Juan A Ballesteros-Cánovas

Clara Rodríguez-Morata

Virginia Garófano- Gómez

Juan M Rubiales
email: jm.rubiales@upm.es

Raúl Sánchez-Salguero

Markus Stoffel
email: markus.stoffel@dendrolab.ch

Dr Juan Antonio

Ballesteros Cánovas
email: juan.ballesteros@dendrolab.ch

Unravelling past flash flood activity in a forested mountain catchment of the Spanish Central System

Keywords: paleohydrology, flash flood, tree ring, dendrogeomorphology, hydrometeorological triggers, Sierra de Guadarrama

Flash floods represent one of the most common natural hazards in mountain catchments, and are frequent in Mediterranean environments. As a result of the widespread lack of reliable data on past events, the understanding of their spatio-temporal occurrence and their climatic triggers remains rather limited. Here, we present a dendrogeomorphic reconstruction of past flash flood activity in the Arroyo de los Puentes stream (Sierra de Guadarrama, Spanish Central System). We analyze a total of 287 increment cores from 178 disturbed Scots pine trees (Pinus sylvestris L.) which yielded indications on 212 growth disturbances related to past flash flood impact. In combination with local archives, meteorological data, annual forest management records and highly-resolved terrestrial data (i.e., LiDAR data and aerial imagery), the dendrogeomorphic time series allowed dating 25 flash floods over the last three centuries, with a major event leaving an intense geomorphic footprint throughout the catchment in 1936. The analysis of meteorological records suggests that the rainfall thresholds of flash floods vary with the seasonality of events. Dated flash floods in the 20 th century were primarily related with synoptic troughs owing to the arrival of air masses from north and west on the Iberian Peninsula during negative indices of the North Atlantic Oscillation. The results of this study contribute considerably to a better understanding of hazards related with hydrogeomorphic processes in central Spain in general and in the Sierra de Guadarrama National Park in particular.

INTRODUCTION

Flash floods represent one of the most common natural hazards in mountain environments [START_REF] Borga | Surveying flash flood response: gauging the ungauged extremes[END_REF][START_REF] Borga | Hydrogeomorphic response to extreme rainfall in headwater systems: flash floods and debris flows[END_REF]. The process is characterized by high flow velocities and important sediment charge causing severe damage and socio-economic losses, especially along the channels and on alluvial fans. In headwater catchments, torrential processes are moreover the main geomorphic agent responsible for landscape evolution [START_REF] Carling | The Noon Hill flash floods, July 17th 1983. Hydrological and geomorphological aspects of a major formative event in an upland landscape[END_REF][START_REF] Foulds | Geomorphological records of extreme floods and their relationship to decadal-scale climate change[END_REF]. From a scientific perspective, the understanding of the temporal dimension of these processes as well as their climatic triggers and subsequent effects on the environment are still a challenge worldwide [START_REF] Schneuwly-Bollschweiler | Hydrometeorological triggers of periglacial debris flows in the Zermatt valley (Switzerland) since 1864[END_REF], probably owing to the frequent occurrence of these processes in sparsely populated areas where archival data and systematic records are usually scarce [START_REF] Mayer | Frequency and spread of debris floods on fans: a dendrogeomorphic case study from a dolomite catchment in the Austrian Alps[END_REF] or unrepresentative (Ayala-Carcedo, 2002). The lack of data on past activity therefore hampers the analysis of flash flood processes and calls for the use of alternative geomorphic approaches [START_REF] Ibsen | The nature, use and problems of historical archives for the temporal occurrence of landslides, with specific reference to the south coast of Britain, Ventnor, Isle of Wight[END_REF][START_REF] Jakob | Debris flow hazard analysis[END_REF]. In this regard, paleohydrologic techniques allow to track the history of past (flash) flood events in ungauged catchments (Baker et al., 2002;Benito et al., 2003;Baker 2008), and consequently, to improve links between process dynamics, climatic conditions and related hazards [START_REF] Kingston | Large-Scale Climatic Controls on New England River Flow[END_REF][START_REF] De | Perspective on climate impacts, snow seasonality and hydrological change in mountain environments[END_REF][START_REF] Merz | Floods and climate: emerging perspectives for flood risk assessment and management[END_REF].

In mountain catchments, trees are frequently present next to torrential channels and on their banks, and can thus be used to reconstruct past flood activity with dendrogeomorphic methods (Stoffel et al., 2010a). The systematic analysis of growth-ring series from trees disturbed by hydrogeomorphic processes [START_REF] Stoffel | Hydrogeomorphic processes and vegetation: disturbance, process histories, dependencies and interactions[END_REF] typically yields valuable records of past events in ungauged mountain catchments with very high spatial and temporal accuracy [START_REF] Shroder | Dendrogeomorphology: review and new techniques of tree-ring dating[END_REF][START_REF] Stoffel | Tree rings and natural hazards: A state-of-the-art[END_REF][START_REF] Stoffel | Dendroecological dating of geomorphic disturbance in trees[END_REF]. Dendrogeomorphical techniques have first been applied in fluvial geomorphology [START_REF] Sigafoos | Botanical evidences of floods and flood-plain deposition[END_REF][START_REF] Sigafoos | Botanical evidence of the modern history of Nisqually Glacier[END_REF]. However, despite that the utility of botanical evidence in paleohydrology has been recognized by many researchers (i.e. Baker 1987;[START_REF] Hupp | Plant ecological aspects of flood geomorphology and paleoflood history[END_REF], tree rings have been used much less frequently as compared to other lines of evidence of paleofloods (e.g., slackwater deposits; Benito and Thorndycraft, 2004). The large potential of dendrogeomorphic tools for the assessment of frequency and magnitude of past events has been demonstrated in recent works (e.g., Ballesteros et al., 2011;[START_REF] Gottesfeld | Flood-plain dynamics of a wandering river: Dendrochronology of the Maurice River, British Columbia, Canada[END_REF][START_REF] Gottesfeld | British Columbia flood scars: maximum flood-stage indicator[END_REF][START_REF] Ruiz-Villanueva | Dendrogeomorphic analysis of flash floods in a small ungauged mountain catchment (Central Spain)[END_REF][START_REF] Schraml | Estimation of debris flood magnitudes based on dendrogeomorphic data and semi-empirical relationships[END_REF][START_REF] St | Paleofloods records for the Red River, Manitoba, Canada derived from anatomical tree-ring signatures[END_REF]. Further research has focused on changes in the spatial geomorphic patterns of processes (Arbellay et al., 2010;[START_REF] Bollschweiler | Dynamics in debris-flow activity on a forested cone -a case study using different dendroecological approaches[END_REF]Stoffel et al., 2008) and on the seasonality and related climatic drivers of hydrogeomorphic processes [START_REF] Schneuwly-Bollschweiler | Hydrometeorological triggers of periglacial debris flows in the Zermatt valley (Switzerland) since 1864[END_REF][START_REF] Stoffel | Rainfall characteristics for periglacial debris flows in the Swiss Alps: past incidences -potential future evolutions[END_REF][START_REF] Stoffel | Dendroecological dating of geomorphic disturbance in trees[END_REF].

Here, we present a case study focusing on the spatio-temporal reconstruction of past flash flood activity in the Arroyo de los Puentes stream (Sierra de Guadarrama National Park, Spanish Central System). We analyze 178 Scots pine trees (Pinus sylvestris L.) disturbed by past flash flood events and couple this data with a large local historical forest management and climatic dataset of the study site to (i) report on the flash flood history of the stream during the last 212 years, and to (ii) identify local meteorological conditions which most likely acted as triggers of flash flood events during the past 83 years for which meteorological records exist locally.

STUDY SITE

The study site is located in the catchment known as Arroyo de los Puentes and its tributaries, located on the northern slope of the Guadarrama Mountains (Sierra de Guadarrama National Park, Spanish Central System, 40° 47´ 37"N, 3° 55´ 14"O; Fig. 1). The catchment covers approximately 2.5Km 2 and extends from the Bola del Mundo at 2258 masl to the alluvial fan area at 1500 masl. The average slope in the main channel is 9° (range: 7° -18°).

(FIGURE 1)

The upper part of the catchment is occupied by extensive accumulations of unconsolidated gneissic materials prone to gelifraction. Water circulation in the source area is mostly sub-superficial, although several well-defined channels exist. The most characteristic geomorphic features in the central part of the catchment (between 1800 and 1600 masl) are linked with torrential activity, such as levees, lobes and well-defined avulsion channels. In this channel segment (at 1660 masl), Majabarca stream joins the main channel and exhibits several trees affected by floods. In the lower part of the catchment (from 1550 to 1480 masl) the valley opens and the gradient decreases, creating an alluvial fan at the confluence of Arroyo de los Puentes with Arroyo de las Pintadas stream. The fan covers an area of approximately 0.03 km 2 ha and its surface is crossed by several channels.

The region is dominated by a Mediterranean climate with continental influence that could be considered as 'humid continental with warm summers' (type 'Dsb': Köppen-Geiger classification, [START_REF] Peel | Updated world map of the Köppen-Geiger climate classification[END_REF]. The study zone also has some Atlantic influence in the regime of rainstorms and is characterized by mild summers and long, cold winters. Average annual precipitation is 1326 mm with maximum rainfall in April, May, October, November and December [START_REF] Aemet | Atlas climático ibérico[END_REF]. Mean annual temperature is 6.5°C at 1890 masl with mean monthly temperatures ranging from 2.9°C in winter to 9.9°C in summer, reaching up to 31.8° C in August and -20.3°C in December, as extreme values.

The study area is located in Montes de Valsaín (hereafter, Valsaín Forest), an extensive, managed P. sylvestris L. forest (10700 ha). The forest is unique in the Mediterranean context owing to the extremely detailed record of forest management interventions started at the end of the 18 th century and its management in general over more than eight centuries (Donés and Garrido, 2001).

MATERIALS AND METHODS

Geomorphic mapping and sampling strategy

A geomorphic characterization of all of the features related to hydrogeomorphic processes in the studied catchment was carried out by combining aerial imagery (cell size: 0.25 m), LiDAR data (cell size: 1 m) and field surveys. All geomorphic features were digitized using ArcGIS TM version 9.3 (ESRI, Redlands, CA, USA, 2009). Disturbed trees -namely wounded, tilted, decapitated or buried trees -located along the channel banks and/or on the fan were sampled following standard procedures in dendrogeomophic studies as described in [START_REF] Stoffel | Dendroecological dating of geomorphic disturbance in trees[END_REF]. Trees with possible disturbances by any process other than hydrogeomorphic (such as rockfall or human activities) were not been included in analysis.

Trees were sampled with increment borers and two increment cores per tree were extracted with sampling positions being chosen according to the nature of the disturbance (Stoffel et al., 2005a). Samples were taken at the contact between the scar edge and the intact wood tissue to make sure that the entire tree-ring record was obtained (Schneuwly et al., 2009a(Schneuwly et al., , 2009b)). In parallel, undisturbed trees were also sampled in the upper and lower parts of the catchment to build a reference chronology and to identify pointer years for a reliable and precise crossdating with disturbed trees [START_REF] Touchan | Dendrochronology course in Valsaín Forest, Segovia, Spain[END_REF]. Tree-ring widths were converted into width indices by standardizing raw data using ARSTAN software [START_REF] Cook | A time series analysis approach to tree-ring standardization[END_REF].

Tree-ring analysis and flash flood chronology reconstruction

Samples were prepared and measured following standard dendrochronological procedures [START_REF] Bollschweiler | Dynamics in debris-flow activity on a forested cone -a case study using different dendroecological approaches[END_REF]2009). Individual growth series were obtained for each tree and cross-dated with the reference chronology, both visually and through statistical procedures (Stoffel et al., 2005a(Stoffel et al., , 2005b)). Signatures of past flash flood activity were then identified on the increment cores and included injuries, callus tissues, compression wood, abrupt growth increase and/or growth suppression. Because Pinus spp. do not form tangential rows of traumatic resin ducts (or TRD, [START_REF] Stoffel | Dating past geomorphic processes with tangential rows of traumatic resin ducts[END_REF], the seasonality of flash flood events had to be based on the position of wound borders within the increment rings.

For the separation of flood signals from noise related with other external processes affecting trees, we applied the weighted index value (W it ) as defined by [START_REF] Kogelnig-Mayer | Possibilities and limitations of dendrogeomorphic time-series reconstructions on sites influenced by debris flows and frequent snow avalanche activity[END_REF]. This index considers the number and the intensity of growth disturbances (GDs) within each tree-ring series and the total number of trees available for the flash flood reconstruction.

The W it is calculated for each year of the reconstruction and considers differences in the intensity of tree reactions to mechanic disturbance related with past events. We screened recent publications to define appropriate thresholds for the identification of past hydrogeomorphic events [START_REF] Corona | Defining optimal sample size, sampling design and thresholds for dendrogeomorphic landslide sampling[END_REF][START_REF] Schneuwly-Bollschweiler | How to improve dating quality and reduce noise in tree-ring based debris-flow reconstructions[END_REF]Stoffel et al., 2011). In addition to the thresholds, we also visually analysed the spatial distribution of the affected trees along the channel and their relationship with geomorphic features [START_REF] Lopez-Saez | Debris-flow activity in abandoned channels of the Manival torrent reconstructed with LiDAR and treering data[END_REF].

Owing to the fact that the anthropogenic influence in the Valsaín Forest is far from being negligible, we also took account of the very detailed and precise official register of forestry works carried out in the catchment since 1940. This additional information allowed discarding doubtful signals in years when intense forestry interventions were realized in the sampled area. Consequently, criteria to define flash flood events was defined as follows: Agency (AEMET). For those events in which a precise dating was not possible, we considered all potential triggering situations with the highest (or majors) rainfall precipitation totals (>30 mm) recorded throughout the year, but we also considered potential rain-on-snow processes.

• W it > 1: EVENT • 1> W it > 0.
Subsequently, 1-, 3-, and 5-day cumulated rainfall totals were analyzed to determine the meteorological characteristics related with flash flood events at the study site. We also evaluated the monthly North Atlantic Oscillation Index (NAO) and synoptic situations prevailing before the occurrence of events.

Hydrometeorological data exists at Navacerrada, a meteorological station located next to the source area of the study site (40° 47´N, 4° 00´W, 1894 masl; Fig. 1) for the period 1933-2011. In addition, we consulted the time series from Segovia's meteorological station (40° 56´N, 4° 10´W, 1005 masl) for which daily data are available since 1894 [START_REF] Aemet | Atlas climático ibérico[END_REF]. Finally, for the definition of the large scale atmospheric flow and weather circulation patterns related with the reconstructed flash flood events, the Hess and Brezowsky Grosswetterlagen (GWL) weather classification system [START_REF] Gestengabe | Katalog der Grosswetterlagen Europas[END_REF][START_REF] Parajka | Seasonal characteristics of flood regimes across the Alpine-Carpathian range[END_REF] was used.

This catalogue provides daily information of patterns over Europe for the period 1881-2004, based on the mean air pressure distribution at both, surface level and 500 hPa level. In this study, we focused on 29 weather types classified into six different groups (, Gestengabe and

RESULTS

Age structure of the forest stand

The reference chronologies developed for the upper and lower parts of the study area covered the periods 1784-2011 and 1748-2011, respectively. The tree-ring growth indices showed for both sites narrow rings in 1829, 1889, 1950, 1972[START_REF] Carling | The Noon Hill flash floods, July 17th 1983. Hydrological and geomorphological aspects of a major formative event in an upland landscape[END_REF]and 1996 (Fig. 2 (Fig. 2). These narrow rings were used for the visual cross-dating of reference trees (undisturbed) with disturbed trees. Data on the pith age at breast height indicated that the 178 trees sampled at Arroyo de los Puentes stream were on average 150 years old (± 77 years). The oldest tree selected for analysis attained sampling height in 1680, whereas the youngest tree reached sampling position in 2000. At least 37 trees were living in 1800, hence the time period considered is well represented by the sampled depth.

(FIGURE 2)

Flash flood reconstruction

A total of 212 GDs were identified in 287 increment cores, most of them in the form of scars (n=91, 42%). Among the other types of GDs, we observed 50 occurrences of reaction wood (23%), 42 trees exhibited an abrupt growth suppression (n=42, 19%), whereas growth releases (n=18, 8%) and callus tissues (n=11, 5%) were much less frequent. From all GDs, 43% were classified as very strong, 34 % as intermediate and 23% as weak reactions.

Based on the number and intensity of reactions as well as the sample depth available for each year of the reconstruction, Fig. 3 shows the computed W it index, which ranks between 0.04 and 36.3. In 28 years, the W it index exceeded the combined threshold of W it >0.5 and GDs ≥3 (Table 1). The historical forest management database indicated that the sampled area, divided into five forest management units, was subjected to periodically forestry interventions, starting in the early 1940s. Table 2 shows the match between forestry activities as well as the reasons for the interventions. At least 3 years (1945, 1950, and 1996), categorized as potential flood events in Table 1, match with years of forestry interventions following "extraordinary events" in the study reach. In the case of 1945 and 1996, these "extraordinary events" were related to large snow falls which provoked widespread forest damage; but the nature of the extraordinary 1950 event is unknown.

As a consequence of these records, the reconstructed flash flood chronology has been reduced to 25 reconstructed events at the study site. The most significant events, in terms of threshold values with a W it > 5, took place in 1936 (W it =36.3), 1954 (W it =10.9) and 1999 (W it =5.18). In addition, three years showed events with indices ranging from 5 ≥ W it ≥ 2, these were in 1813 (W it =3.7), 1966 (W it =2.2) and 1956 (W it =2.0). Nine years showed a W it comprised between 1 and 2 (i.e., 1802, 1869, 1880, 1900, 1926, 1952, 1973, 1984 and 2004), whereas in 10 cases the W it ranked between 0.5 and 1 (i.e., 1830, 1853, 1858, 1871, 1873, 1876, 1884, 1906, 1933, and 1947). Considering the flash flood activity in the Arroyo de los Puentes catchment since 1802, the reconstructed frequency was 0.12 events/year. The partial frequency increased in 1870s with 5 flash floods in 10 years and between 1940s and 1950s with 4 events between 1947 and 1956 (Fig. 4).

(FIGURE 3) 4) Fig. 5A illustrates the three main patterns of GDs distribution at the study site. During more than two-thirds of all events, GDs were limited to trees growing along the Arroyo de los Puentes stream, whereas disturbance signals were missing completely in its tributary (70% of all events, i.e., 1802, 1813, 1830, 1853, 1858, 1869, 1871, 1874, 1876, 1880, 1884, 1900, 1906, 1936, 1945, 1947, 1950, 1956, 1973, 1984, 1996 and 2004). In four cases (17%; 1926, 1933, 1966 and 1944), GDs were observed in the Arroyo de los Puentes and the Majabarca streams.

(TABLE 1) (TABLE 2) (FIGURE
In three cases (13%, i.e. 1952, 1954 and 1999) GDs were only recorded in Majabarca stream.

(FIGURE 5)

Hydrometeorological triggers and synoptic situation

The assessment of the seasonality of flash floods and/or the calendar-dating of reconstructed events was based on the intra-annual position of the damages in the tree-ring records, historical archives, contemporary newspaper reports, as well as on literature existing on flash floods in the wider study region. This approach allowed a realistic calendar-dating of seven 20 th century flash floods (i.e., 1900, 1936, 1952, 1954, 1956, 1984 and 1999). By contrast, local documentation or indirect evidence of events could not be found for the events which occurred in 1906, 1926, 1933, 1947, 1966, 1973 and 2004. In the case of the calendar-dated flash floods, Table 3 shows the main hydrometeorological variables analyzed as well as the distribution of events throughout the year. It becomes obvious that flash floods are mostly limited to the moist seasons with reconstructed flash flood events occurring mainly in autumn-winter (October 1900, February 1936, March 1952, November 1954, January 1956), and a unique flash flood in spring (May 1984) and another one during summer (September 1999). For events without exact monthly date (Table 4), the most likely seasonality were October-November (7 potential cases), May-April (5 potential cases) and March (1947). With exception of only two cases (1954 and 2004), all the events were related to a negative monthly NAO-index.

(TABLE 3) (TABLE 4)
Analysis of rainfall thresholds related with the likely triggers of reconstructed events (period with available data: n = 14) showed differences in precipitation which can likely be related to the seasonality of storms (Fig. 6). Despite the small sample used for analysis, we observed that autumn and winter flash floods were clearly characterized by higher rainfall totals than those occurring in spring or summer. The thresholds for winter events were 54.2 mm, 87.7 mm, and 143.0 mm for the largest 1-, 3-, and 5-day rainfall, respectively. Spring flash floods were characterized by rainfall thresholds of 32.4, 56.6, and 67.6 mm for the same durations, whereas events in autumn had thresholds of 67.9, 139.8, and 146.6 mm for the largest 1-, 3-, and 5-day rainfall totals, respectively. For the case of the unique summer event (1999), the triggering rainfall showed totals of 27, 27.3, and 35.8 mm for the 1-, 3-, and 5-day rainfall, respectively. The very small 1-and 3-day rainfall totals recorded at Navacerrada station located next to the source areas of flash floods may point to a short-lived, but intense storm as the trigger of the flash flood.

The GWL weather classification system and analysis of the isobaric maps at 500hPa and surface level associated with the precipitation records triggering the identified events suggested different weather patterns (Figure 7). Meridional patterns were identified in three cases (i.e. 1952, 1984 and 1999), whereas zonal patterns were observed in two cases (i.e. 1900,and 1954). Two more cases showed a deep low pressure system in Europe (i.e. 1936) and a mixed pattern affected by a north-west air masses (i.e. 1956). Therefore, the identified meteorological causes probably related to the reconstructed events are the following: frontal systems (69% of cases), combination of frontal systems and snowmelt processes (23 % of cases) and convective storms (8% of cases), respectively.

(FIGURE 6) (FIGURE 7)

DISCUSSION

In this study we provide a 212-year chronology of flash floods based on tree-rings reconstructions performed in the ungauged Arroyo de los Puentes catchment located in the Sierra de Guadarrama National Park (Spanish Central System) and characterize the meteorological triggers of those past flash floods. We report results derived from the analysis of 287 increment cores of 178 living P. sylvestris L. trees. The growth-ring records of these trees allowed the detection and dating of 212 GDs attributed to 25 flash flood events which occurred since 1802.

The reconstructed flash flood chronology complements the existing flood record for the contributing highland catchment of the Eresma River and represents the longest and most continuous (annually-resolved), non-systematic record of flash flood activity and related climatic factors in this area. Results are therefore expected to greatly support the definition of flash flood hazard zones in an area characterized by intense tourist activity (e.g., the recently founded Sierra de Guadarrama National Park) and by a concentration of linear transport infrastructures (i.e. railway lines, roads, power lines and hiking paths).

Our reconstruction points to a flash flood frequency of 0.12 events/year in the Arroyo de los Puentes stream for the last 212 years. Individual event years match with the chronology of recorded (i.e., gauged) flood events in the Eresma River. We have successfully identified evidences of the two major lowland floods of 1956 and 1966 in the trees from one of the contributing highland catchments. Moreover, existing archives and local contemporary newspapers reports provide further evidence for extreme hydroclimatic conditions in the Eresma River and other parts of central Spain in 1936Spain in , 1956Spain in , 1966Spain in , 1984Spain in , and 1996 (see Table 5). These years also appear, without any exception, in the dendrogeomorphic flash flood reconstruction in the Arroyo de los Puentes stream, and thereby indirectly confirm the robustness of the time series reconstructed in this study.

(TABLE 5)

Our results showed two periods with intense flash flood activity between 1870 and 1884 and between 1947 and 1956. On the basis of historical literature from the 19 th century, the first period with very frequent flash floods can be related to meteorological anomalies characterized by significant oscillations in atmospheric pressure, high snow cover in winter and an excessive number (103) of intense precipitation events including storms [START_REF] Breñosa | Guía y descripción del Real Sitio de San Ildefonso[END_REF]. The reconstruction from Arroyo de los Puentes is in agreement with observations from the Tagus basin (central Spain) where Benito et al., (2003aBenito et al., ( , 2003b) observed a similar increase in the frequency of flood events. The most intense event was detected in 1936 (W it = 36.38 and GDs = 25). This observation is also in agreement with other studies in the Spanish Central System [START_REF] Potenciano | Las inundaciones históricas en el centro-sur de la Península Ibérica. Condicionantes geomorfológicos y climáticos[END_REF][START_REF] Ruiz-Villanueva | Characterization of flash floods in small ungauged mountain basins of central Spain using an integrated approach[END_REF], and hence confirm the extraordinary nature of this flash flood.

In methodological terms, the preferred use of the W it index value (Kogelnig- [START_REF] Kogelnig-Mayer | Possibilities and limitations of dendrogeomorphic time-series reconstructions on sites influenced by debris flows and frequent snow avalanche activity[END_REF] -instead of the Shroder index [START_REF] Shroder | Dendrogeomorphology: review and new techniques of tree-ring dating[END_REF] -was not only justified by the large amount of flood scars found in the trees along the channel (42%), but also was needed in view of the hypothetical influence of intense forest management on accidental growth anomalies in trees. The preferential weighting of strong signals (i.e., flood scars), as suggested by [START_REF] Stoffel | Dendroecological dating of geomorphic disturbance in trees[END_REF], did not only filter low-frequency noise in the time series more efficiently, but also yielded a reconstruction for which local flash flood events seem to converge with the dates of regional floods. Despite the convincing match between local flash floods and regional flood events, it is necessary to emphasize that dendrogeomorphic reconstructions of past geomorphic activity (of any nature) will always remain as minimum frequencies [START_REF] Bollschweiler | Debris-flood reconstruction in a pre-alpine catchment in Switzerland based on tree-ring records of coniferous and broadleaved trees[END_REF]. Therefore, some events could remain missed even if efforts are undertaken to optimize sample depth to minimize noise and to maximize signals (e.g., [START_REF] Corona | How much of the real avalanche activity can be captured with tree rings? An evaluation of classic dendrogeomorphic approaches and comparison with historical archives[END_REF][START_REF] Corona | Seven centuries of avalanche activity at Echalp (Queyras massif, southern French Alps) as inferred from tree rings[END_REF][START_REF] Corona | Defining optimal sample size, sampling design and thresholds for dendrogeomorphic landslide sampling[END_REF][START_REF] Schneuwly-Bollschweiler | How to improve dating quality and reduce noise in tree-ring based debris-flow reconstructions[END_REF][START_REF] Stoffel | Mass movements and tree rings: A guide to dendrogeomorphic field sampling and dating[END_REF]. This assertion is further underlined by the fact that internal scars (i.e., blurred or hidden injuries which can no longer be seen on the stem surface) cannot be detected easily in P. sylvestris L., as this speciesunlike other conifers -does not form TRD around wounds [START_REF] Stoffel | Dating past geomorphic processes with tangential rows of traumatic resin ducts[END_REF], thus hampering the indirect dating of scars with seasonal precision [START_REF] Stoffel | Snow avalanche and rockfall impacts leave different anatomical signatures in tree rings of Larix decidua[END_REF]. Our study is also limited by the fact that we were unable to show the occurrence of repeated flash floods in the same year [START_REF] Schneuwly-Bollschweiler | Hydrometeorological triggers of periglacial debris flows in the Zermatt valley (Switzerland) since 1864[END_REF].

Analysis of hydrometeorological triggers was based on 14 flash floods covering the span of the available precipitation records. The period of flash flood occurrence could be narrowed down in 7 out of 14 cases using historical archives and/or the intra-seasonal position of scars in trees. In 1956 and 1966, the exact date of events has been documented and could thus been used directly. The intra-annual position of injuries was most helpful for the precise dating of the 1936 and 1999 events. For the remaining events we used newspaper reports and local archives (see above). In case that the seasonality of flash floods was known, we assumed that events were triggered by the greatest rainfall recorded within the time window suggested by the intra-annual position of scars in trees and/or by archival records [START_REF] Schneuwly-Bollschweiler | Hydrometeorological triggers of periglacial debris flows in the Zermatt valley (Switzerland) since 1864[END_REF]. For the 6 remaining cases, neither direct nor indirect references to floods or intense rainfalls could be found, so that the analysis of potential triggers was extended to the entire year. As a result of the proximity of the Navacerrada station (only 900 m away from the source area of flash floods and located at a similar altitude) and the good quality of the data, the climate dataset analyzed in this study is highly reliable. Because the resolution of the precipitation time series was daily (and not hourly), we could not perform an intensityduration-frequency (IDF) analysis for each event, but still believe that the obtained 1-, 3-, and 5-day rainfall thresholds presented in this study are reliable enough to understand average rainfall thresholds involved in the triggering of flash floods at the study site [START_REF] Kundzewicz | Floods at the northern foothills of the Tatra Mountains -A Polish-Swiss research project[END_REF].

Our results suggest that daily rainfall threshold strongly depend on the season of flash floods. In that sense, flash floods occurring in spring were typically related to lower rainfall totals as compared to flash floods in autumn and winter. This observation also clearly points to the predominant role of snowmelt processes and soil moisture structure triggering flash floods in winter and spring [START_REF] De | Perspective on climate impacts, snow seasonality and hydrological change in mountain environments[END_REF], confirming the conclusion provided by [START_REF] Marchi | Characterisation of selected extreme flash floods in Europe and implications for flood risk management[END_REF]. In the case of the largest event on record (1936), we also demonstrate that, depending on the seasonality of flash floods, events will not necessarily be triggered by the most important hydrometeorological event and/or the largest precipitation recorded during a year.

The time series presented in this study also highlight the importance of independent, seasonal flood-frequency analyses to obtain reliable flood hazard assessments (Baratti et al., 2012;[START_REF] Merz | Floods and climate: emerging perspectives for flood risk assessment and management[END_REF].

At the study site, flash flood events typically took place between autumn and spring, which is in agreement with findings from Ruiz-Villanueva et al. ( 2013) who obtained similar conclusions for headwater catchments in the western part of the Spanish Central System. The main synoptic situations related with the triggering of flash floods were cold Atlantic and continental air masses transported to the study site by low pressure systems generating synoptic troughs over the Iberian Peninsula [START_REF] Tomás | Circulation weather types and cloud-to-ground flash density over the Iberian Peninsula[END_REF]. In 84% of the cases, flash floods were related to negative monthly NAO indices, showing in six cases moderate to intense negative values (with NAO < -1). Our observations are in agreement with results from other mountain catchments in central Spain (e.g., Benito et al., 2003aBenito et al., , 2003b[START_REF] Benito | Impactos sobre los riesgos naturales de origen climático[END_REF][START_REF] Benito | Palaeoflood and floodplain records from Spain: Evidence for long-term climate variability and environmental changes[END_REF] and support the idea that intense hydrometeorological activity in central Spain is related to the antecedent NAO mode [START_REF] Cortesi | High resolution reconstruction of monthly precipitation of Iberian Peninsula using circulation weather types[END_REF][START_REF] Salgueiro | Flood magnitudes in the Tagus River (Iberian Peninsula) and its stochastic relationship with daily North Atlantic Oscillation since mid-19th Century[END_REF][START_REF] Trigo | North Atlantic oscillation influence on precipitation, river flow and water resources in the Iberian Peninsula[END_REF]. The weather circulation patterns related with flash floods dated in this contribution confirm the casuistry associated with flash flood events in the study area, which was slightly weighted toward Meridional circulation patterns. The mixed Central European (TM) pattern was in addition responsible in the triggering of the most intense event defined in this study. Our results areconfirmed by observations obtained by [START_REF] Parajka | Seasonal characteristics of flood regimes across the Alpine-Carpathian range[END_REF] who related Meridional circulation patterns with precipitation events and related floods in northern Italy. However, the lower number of event defined in this study is presumably due to the nature of this study (local catchment rather regional study) and thus avoids to draw robust conclusions in this regard

CONCLUSIONS

Flash flood activity in the Arroyo de los Puentes catchment (northern slope of Sierra Guadarrama) was reconstructed during the last 212 years using GDs from a large number of disturbed trees (P. sylvestris L.). The reliability of the presented flash-flood chronology was, therefore, supported by a large climatic and historical forest dataset. Despite the potential limitation related with the use of indirect proxies, this study is supporting the usefulness of the dendrogeomorphic approach for the understanding of past climate-floods in managed mountain areas. Our results describe the complexity of processes in the catchment, as both snowmelt and soil humidity have a major role as triggers of past geomorphic events. In addition, we detected that periods with more intense and frequent flash flood activity were related to negative NAO phases. Therefore, this study provide crucial data to gain a better understanding about the spatio-temporal patterns of flash flood occurrence and their links with climate in the Sierra de Guadarrama National Park.

Arbellay, E., Stoffel, M., Bollschweiler, M., 2010. Dendrogeomorphic reconstruction of past debris-flow activity using injured broad-leaved trees. Earth Surf. Processes. 35, 399-406. Ayala-Carcedo, F.J., 2002. Análisis de riesgos por movimientos de ladera, in: Ayala-Carcedo, F.J., Olcina, J. (Eds.), Riesgos naturales. Ed Ariel, Barcelona, pp. 379-407. Baker, V.R., 1987. Paleoflood hydrology and extraordinary flood events. J. Hydrol. 96, 79-99. Baratti, E., Montanari, A., Castellarin, A., Salinas, J.L., Viglione, A., Bezzi, A., 2012. Estimating the flood frequency distribution at seasonal and annual time scales. Hydrol. Earth Syst. Sci. 16, 4651-4660. Benito, G., Thorndycraft, V.R., 2004. Systematic, palaeo-flood (1829, 1889, 1950, 1972, 1986 and 1996).

LIST OF FIGURES

  5 without forestry work at the study area and more than 3 GDs: EVENT • 1> W it > 0.5 with forestry work at the study area: NO EVENT • W it < 0.5: NO EVENT 3.3. Hydrometeorological analysis of reconstructed flash floods Different historical archives have been screened to refine the (sub-) annually reconstructed dates of flash floods to the month or even the exact date. The sources analyzed included the archives from the Valsaín Forest documentary centre, national and local newspapers, local literature and meteorological ephemerides from the Spanish Meteorological

Werner, 2005 ;

 2005 Parajka et al., 2010): 1. Zonal West: WA, WZ, WS, WW 2. Mixed: SWA, SWZ, NWA, NWZ, HM, BM 3. Mixed Central Europe (CE): TM 4. Meridional North (N): NA, NZ, HNA, HNZ, HB, TRM 5. Meridional Northeast and East (NE, E): NEA, NEZ, HFA, HFZ, HNFA, HNFZ 6. Meridional Southeast and South (SE, E): SEA, SEZ, SA, SZ, TB, TRW

  Baker, V.R., 2008. Paleoflood hydrology: Origin, progress, prospects. Geomorphology 101, 1-13. Baker, V.R., Webb, R.H., House, P.K., 2002. The scientific and societal value of paleoflood hydrology, in: House, P.K., Webb, R.H., Baker, V.R., Levish, D.R. (Eds.), Ancient Floods, Modern Hazards: Principles and Applications of Paleoflood Hydrology, Water Science and Application, vol. 5. Am. Geophys. U., Washington, D.C., pp. 1-19. Ballesteros, J.A., Eguibar, M., Bodoque, J.M., Díez, A., Stoffel, M., Gutiérrez, I., 2011. Estimating flash flood discharge in an ungauged mountain catchment with 2D hydraulic models and dendrogeomorphic paleostage indicators. Hydrol. Process. 25, 970-979.

  and historical data for the improvement of flood risk estimation. CSIC, Madrid.Benito, G., Sopeña, A., Sánchez-Moya, Y., Machado, M.J.,Pérez-González, A., 2003a. Palaeoflood record of the Tagus River (central Spain) during the Late Pleistocene and Holocene. Quaternary Sci.Rev. 22, 1737-1756. Benito, G., Díez-Herrero, A., Fernández De Villalta, M., 2003b. Magnitude and frequency of flooding in the Tagus basin (central Spain) over the last millennium. ClimateChange. 58,

Fig. 1 .

 1 Fig. 1. General scheme of the study area showing the distribution of the anthropogenic

Fig. 2 .

 2 Fig. 2. Reference chronologies of Pinus sylvestris L. for Arroyo de los Puentes stream, obtained

Fig. 3 .

 3 Fig.3. Event-response histograms showing flash-flood induced growth responses from the

Fig. 4 .

 4 Fig. 4. Flash flood chronology and 5-year moving average obtained at Arroyo de los Puentes

Fig. 5 .

 5 Fig. 5. Reconstructed patterns of past flash flood events. (A) Maps showing the location of

Fig. 7 .

 7 Fig. 7. Synoptic situation (at surface level and 500hPa level) related with dated flash flood

Table 2 .

 2 Match between forestry actions, reasons for the works and responses in trees. Grey shading highlight years in which reconstructed disturbances in trees match with forestry work. Letters indicate the causes of the felling of trees during forestry actions in the catchment (Eextraordinary, R-regular, D-dead trees and I-insects,).

	Year	Sample depth	Growth disturbance Medium Weak intensity High intensity intensity	W it	Cause
	1940	112		3	0	0.40	R
	1945	120	3	0	0	0.53	E
	1950	126	3	0	0	0.50	E
	1962	130	1	0	0	0.05	R
	1972	149	0	3	0	0.30	D
	1996	155	0	1	4	0.55	E
	2000	155	2	0	1	0.33	E and I

Table 3 .

 3 Hydrometeorological triggers (in grey analysis based on Segovia station) related to reconstructed events matching with historical archives.

	Year	Date	Temp °C	Rainfall	Rainfall	Rainfall 5-	NAO-	GWL	References/Newspaper/indicator
			(max/min)	1-day	3-day	day	index	system
	1900	13/10	15° / 7°	73	123	123	-0.57	WZ	La Vanguardia-14/10/1900-pp2
	1936	29/02	5° / -°	36	47	91	-2.44	TM	Seasonality based on anatomical interpretation
	1952	31/03	-° / 1.3°	101.8	143.5	276.9	-1.14	HNZ	El Adelantado de Segovia-01/04/1952
	1954	8/11	5.8° / 5°	55	145.7	149.3	1.54	WZ	ABC-07/11/1954-pp 039
	1956	18/01	-° / 1.2°	39	104.2	112.1	-0.76	NWZ	Díez-Herrero et al., (2008)
	1984	14/05	2.2° / -°	27	38.5	38.5	-2.34	HFZ	Morales and Ortega ( 2002)
	1999	18/09	12° / 7°	27	27.3	35.8	-0.54	TRW	Seasonality based on anatomical interpretation

Table 4 .

 4 Potential hydrometeorological triggers (in grey analysis based on Segovia station) related to reconstructed events not recorded on historical archives.

	Year	Date	Rainfall	Temp °C (max/min) Rainfall	Rainfall	NAO-index
			1-day		3-day	5-day
	1906	12/05	27	15°/ 6°	39	79	-1.75
		12/11	30	8°/ 5°	49	64	0.7
	1926	23/10	27	17°/ 2°	52	52	-4.16
	1933	20/04	40	5°/ 3°	40	40	-0.23
		9/10	54 (12h)	-°/ 3°	54	54	-2.03
		21/10	74 (12h)	10°/ 0°	82	82	-2.03
	1947	05/03	40	6.4°/ 1°	57	92	-1.2
	1966	03/10	104.9	9.2°/9°	132.5	136	-2.89
		07/11	94.7	0.8°/0.2°	302.5	321.1	-0.07
	1973	05/11	46.7	6°/ 3.8°	141.5	162.8	-0.26
		03/05	30.2	4° / 0.6°	73.6	73.6	0.37
		19/05	41.6	3° / 0°	106.8	128.7	0.37
	2004	23/05	31	10.2° / 6.6°	49.5	49.5	0.19

Acknowledgments

This study has been funded by the projects IDEA-GesPPNN-163/2010 (OAPN -Organismo Autónomo de Parques Nacionales, Ministerio de Agricultura, Alimentación y Medio Ambiente) and CGL2010-19274 (projects MAS Dendro-Avenidas ) of the Spanish Ministry of Economy and Competitiveness. The authors sincerely acknowledge CENEAM, Montes de Valsaín and OAPN for providing information for this study. Specifically, they want to thank Javier Donés (Director of Montes de Valsaín), Andres Díez-Herrero, Ramzi Touchan and Dave Meko for their support during the study, as well as all participants of the Tree-ring summer school carried out in 2012 in Valsaín (Touchan et al. 2013). R. Sánchez-Salguero thanks the financial support from University of Córdoba-Campus de Excelencia Internacional-ceiA3.

LIST OF TABLES

• Paleohydrology improve the understanding of flash-floods in ungauged catchment

• Dendrogeomorphology is the most suitable approach in mountain forested areas

• We report the flash flood history of the mountain stream during the last 212 years

• We identify the meteorological conditions acting as triggers of flash flood events


Lucas Chesnel
email: lucas.chesnel@cmap.polytechnique.fr

AND Nuutti Hyvönen
email: nuutti.hyvonen@aalto.fi

Stratos Staboulis
email: stratos.staboulis@helsinki.fi

CONSTRUCTION OF INVISIBLE CONDUCTIVITY PERTURBATIONS FOR THE POINT ELECTRODE MODEL IN ELECTRICAL IMPEDANCE TOMOGRAPHY

Keywords: 2010 Mathematics Subject Classification. 35J25, 65N20, 35R30 Electrical impedance tomography, point electrode model, invisibility, localized current sources, elliptic boundary value problems, complete electrode model

come

Construction of indistinguishable conductivity perturbations for the point electrode model in electrical impedance tomography

Lucas Chesnel, Nuutti Hyvönen, Stratos Staboulis

Introduction

Electrical impedance tomography (EIT) is a noninvasive imaging technique with applications, e.g., in medical imaging, process tomography, and nondestructive testing of materials [START_REF] Borcea | Electrical impedance tomography[END_REF][START_REF] Cheney | Electrical impedance tomography[END_REF][START_REF] Uhlmann | Electrical impedance tomography and Calderón's problem[END_REF]. The aim of EIT is to reconstruct the conductivity distribution inside the examined physical body D ⊂ R d , d ≥ 2, from boundary measurements of current and potential. From the purely theoretical standpoint, the inverse problem of EIT, also known as the inverse conductivity problem, corresponds to determining the strictly positive conductivity σ : D → R in the elliptic equation [START_REF] Astala | Calderón's inverse problem for anisotropic conductivity in the plane[END_REF] div(σ∇u) = 0 in D from knowledge of the corresponding Neumann-to-Dirichlet (or Dirichlet-to-Neumann) map at the object boundary ∂D. This formulation corresponds to the idealized continuum model (CM), which is mathematically attractive in its simplicity but lacks a straightforward connection to practical EIT measurements that are always performed with a finite number of contact electrodes. On the other hand, it is widely acknowledged that the most accurate model for real-life EIT is the complete electrode model (CEM), which takes into account electrode shapes and contact resistances at electrode-object interfaces [START_REF] Cheng | Electrode models for electric current computed tomography[END_REF][START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF]. Instead of the CM or the CEM, in this work we employ the point electrode model (PEM) that treats the current-feeding electrodes of EIT as delta-like boundary current sources and models the potential measurements as pointwise evaluations of the corresponding solution to (1) at the electrode locations [START_REF] Hanke | Justification of point electrode models in electrical impedance tomography[END_REF]. (To make such a model well defined, relative potential measurements need to be considered due to the singularities induced by the localized current inputs; see Section 2 for the details.) The PEM is a reasonably good model for practical EIT measurements -in particular, superior to the CM (cf. [START_REF] Hyvönen | Approximating idealized boundary data of electric impedance tomography by electrode measurements[END_REF]) -if the electrodes are small compared to the size of the imaged body: The discrepancy between the PEM and the CEM is of the order O(h 2 ), with h > 0 being the maximal diameter of the electrodes [START_REF] Hanke | Justification of point electrode models in electrical impedance tomography[END_REF]. In this work we develop a method, introduced by [START_REF] Bonnet-Ben Dhia | Obstacles in acoustic waveguides becoming "invisible" at given frequencies[END_REF] in the context of acoustic waveguides (see also [START_REF] Nazarov | Asymptotic expansions of eigenvalues in the continuous spectrum of a regularly perturbed quantum waveguide[END_REF] for the original idea), for building invisible perturbations of a given reference conductivity for an arbitrary but fixed electrode configuration modeled by the PEM.

To put our work into perspective, let us briefly review the main results on the unique solvability of the inverse conductivity problem, which was proposed to the mathematical community by Calderón in [START_REF] Calderón | On an inverse boundary value problem[END_REF]. For d ≥ 3, the global uniqueness for C 2 -conductivities was proven in [START_REF] Sylvester | A global uniqueness theorem for an inverse boundary value problem[END_REF]; an extension of this result to the case Lipschitz conductivities can be found in [START_REF] Haberman | Uniqueness in Calderón's problem with Lipschitz conductivities[END_REF]. The first global uniqueness proof in two spatial dimensions was given by [START_REF] Nachman | Global uniqueness for a two-dimensional inverse boundary value problem[END_REF] for C 2 -conductivities; subsequently, [START_REF] Astala | Calderón's inverse conductivity problem in the plane[END_REF] proved uniqueness for general L ∞ -conductivities. All these articles assume the Cauchy data for [START_REF] Astala | Calderón's inverse problem for anisotropic conductivity in the plane[END_REF] are known on all of ∂D, but the partial data problem of having access only to some subset(s) of ∂D has also been tackled by many mathematicians; see, e.g., [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF][START_REF] Kenig | The Calderón problem with partial data[END_REF] and the references therein. From our view point, the most important uniqueness results are presented in [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF][START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF] where it is shown that any two-dimensional conductivity that equals a known constant in some interior neighborhood of ∂D is uniquely determined by the PEM measurements with countably infinite number of electrodes. This manuscript complements [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF][START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF] by constructively showing that for an arbitrary, but fixed and finite, electrode configuration there exists a perturbation of the unit conductivity that is invisible for the EIT measurements in the framework of the PEM. For completeness, it should also be mentioned that recently the nonuniqueness of the inverse conductivity problem has been studied for finite element method (FEM) discretizations with piecewise linear basis functions [START_REF] Lassas | Inverse problems and invisibility cloaking for FEM models and resistor networks[END_REF]. (Take note that we restrict our attention exclusively to isotropic conductivities because it is well known that the inverse conductivity problem is not uniquely solvable for anisotropic conductivities; see, e.g., [START_REF] Astala | Calderón's inverse problem for anisotropic conductivity in the plane[END_REF][START_REF] Sylvester | An anisotropic inverse boundary value problem[END_REF] and the references therein.)

Our constructive proof for the existence of 'invisible conductivities' is based on introducing a small perturbation of a given reference conductivity and applying a suitable fixed point iteration. The general idea originates from [START_REF] Cardone | Asymptotic behaviour of an eigenvalue in the continuous spectrum of a narrowed waveguide[END_REF][START_REF] Nazarov | Asymptotic expansions of eigenvalues in the continuous spectrum of a regularly perturbed quantum waveguide[END_REF][START_REF] Nazarov | Eigenvalues of the laplace operator with the Neumann conditions at regular perturbed walls of a waveguide[END_REF][START_REF] Nazarov | Trapped waves in a cranked waveguide with hard walls[END_REF][START_REF] Nazarov | Enforced stability of an eigenvalue in the continuous spectrum of a waveguide with an obstacle[END_REF][START_REF] Nazarov | Enforced stability of a simple eigenvalue in the continuous spectrum of a waveguide[END_REF], where the authors develop a method for constructing small regular and singular perturbations of a waveguide that preserve the multiplicity of the point spectrum on a given interval of the continuous spectrum. Subsequently, the same approach has been adapted in [START_REF] Bonnet-Ben Dhia | Obstacles in acoustic waveguides becoming "invisible" at given frequencies[END_REF] (see also [START_REF] Bonnet-Ben Dhia | Underwater topography "invisible" for surface waves at given frequencies[END_REF] for an application to a water wave problem) to demonstrate the existence of regular perturbations of a waveguide allowing waves at given frequencies to pass through without any distortion or with only a phase shift. Recently, in [START_REF] Bonnet-Ben Dhia | Transmission eigenvalues and far field invisibility for a finite number of incident/scattering directions[END_REF] the methodology has also been used in inverse obstacle scattering to construct defects in a reference medium which are invisible to a finite number of far field measurements. Although our technique is, in principle, applicable in any spatial dimension d ≥ 2 and for an arbitrary smooth enough reference conductivity, a necessary intermediate result about linear independence of certain auxiliary functions (cf. Lemma 3.2) is proved here only for d = 2 and the unit reference conductivity, making the main theoretical result of this work two-dimensional. We implement the constructive proof as a numerical algorithm that is capable of producing conductivities that are indistinguishable from the unit conductivity for a given set of point electrodes. We also numerically demonstrate that such conductivities are almost invisible for small electrodes within the CEM as well.

It should be emphasized that the existence of invisible conductivity perturbations for a finite number of electrodes is not very surprising: The space of admissible conductivities is infinite dimensional whereas EIT measurements with N + 1 electrodes include only N (N + 1)/2 degrees of freedom. Be that as it may, we are not aware of any previous works providing examples on this nonuniqueness that is of importance in practical EIT.

This text is organized as follows. Section 2 introduces the setting for our analysis and a representation of the relative PEM measurements with respect to a certain simple basis of electrode current patterns. The general scheme for constructing invisible conductivity perturbations is presented in Section 3 and, subsequently, Section 4 adds the missing theoretical piece in two dimensions. The main result of the article is presented in Theorem 4.1. The algorithmic implementation of the constructive existence proof is considered in Section 5 and, finally, the numerical examples are presented in Section 6.

Setting

Let D ⊂ R d , d ≥ 2, be a simply connected and bounded domain with a C ∞ -boundary. Introduce a real valued reference conductivity σ 0 ∈ C ∞ (D) such that σ 0 ≥ c > 0 in D. (Throughout this text C, c > 0 denote generic constants that may change from one occurrence to the next.) Consider the Neumann boundary value problem [START_REF] Astala | Calderón's inverse conductivity problem in the plane[END_REF] div(σ

0 ∇u 0 ) = 0 in D, ν • σ 0 ∇u 0 = f on ∂D for a current density f in (3) H s ⋄ (∂D) := {g ∈ H s (∂D) | g, 1 ∂D = 0}, with some s ∈ R.
Here and in what follows, ν denotes the unit normal vector of ∂D orientated to the exterior of D. We observe that the dual of H s ⋄ (∂D) is realized by (4) H -s (∂D)/R := H -s (∂D)/span{1}, s ∈ R.

It follows from standard theory of elliptic boundary value problems that (2) has a unique solution u 0 ∈ H s+3/2 (D)/R, which is smooth away from ∂D due to the interior regularity for elliptic equations.

In particular, for any compactly embedded domain Ω ⋐ D it holds that (5)

c u 0 H 1 (Ω)/R ≤ u 0 H s+3/2 (D)/R ≤ C f H s ⋄ (∂D)
, where c = c(s, Ω) and C = C(s) [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]. Let us next consider a perturbed conductivity [START_REF] Borcea | Electrical impedance tomography[END_REF] σ ε := σ 0 + εκ, where κ ∈ L ∞ (D) is compactly supported in D and ε > 0 is such that σ ε ≥ c > 0 in D. The corresponding perturbed Neumann boundary value problem

(7) div(σ ε ∇u ε ) = 0 in D, ν • σ ε ∇u ε = f on ∂D,
with f ∈ H s ⋄ (∂D) for some s ∈ R, has a unique solution in (H s+3/2 (D) ∩ H 1 loc (D))/R 1 (cf., e.g., [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF][START_REF] Hanke | Convex backscattering support in electric impedance tomography[END_REF]).

We note that the relative Neumann-to-Dirichlet map

Λ ε -Λ 0 : f → (u ε -u 0 )| ∂D , D ′ ⋄ (∂D) → D(∂D)/R is well-defined.
Here, the mean-free distributions D ′ ⋄ (∂D) and the quotient space of smooth boundary potentials D(∂D)/R are defined in accordance with (3) and ( 4). This regularity result can be deduced from standard elliptic theory (cf., e.g., [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF][START_REF] Hanke | Convex backscattering support in electric impedance tomography[END_REF]) and follows from the fact that Λ ε and Λ 0 are pseudodifferential operators with the same symbol on ∂D because σ ε -σ 0 vanishes in some interior neighborhood of ∂D. One can also prove the following symmetry property for the relative Neumann-to-Dirichlet map (see, e.g., [START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF]Theorem 2.1]).

Proposition 2.1. For all ϕ, ϕ ′ ∈ D ′ ⋄ (∂D), we have ϕ, (Λ ε -Λ 0 )ϕ ′ ∂D = ϕ ′ , (Λ ε -Λ 0 )ϕ ∂D = D (σ 0 -σ ε )∇u ε ϕ • ∇u 0 ϕ ′ dx,
where u ε ϕ , u 0 ϕ ′ denote respectively the solutions of ( 7), [START_REF] Astala | Calderón's inverse conductivity problem in the plane[END_REF] with f equal to ϕ, ϕ ′ . Now, consider an observer who can impose currents between pairs of (small) electrodes located at x 0 , . . . , x N ∈ ∂D and measure the resulting voltages for both conductivities σ 0 and σ ε . According to the PEM [START_REF] Hanke | Justification of point electrode models in electrical impedance tomography[END_REF], this corresponds to knowing all elements in the matrix of relative measurements

M (σ ε ) ∈ R N ×N defined via (8) M ij (σ ε ) = δ i -δ 0 , (Λ ε -Λ 0 )(δ j -δ 0 ) ∂D , i, j = 1, . . . , N,
where δ n ∈ H -(d-1)/2-η (∂D), η > 0, stands for the Dirac delta distribution supported at x n . It follows from Proposition 2.1 that M (σ ε ) is symmetric and, furthermore,

M ij (σ ε ) = D (σ 0 -σ ε )∇u ε i • ∇u 0 j dx, where u 0 n , u ε n ∈ (H -(d-4)/2-η (D) ∩ H 1 loc (D))/R, η > 0, (9)
are the solutions of the Neumann problems [START_REF] Cheng | Electrode models for electric current computed tomography[END_REF] div(σ

0 ∇u 0 ) = 0 in D, ν • σ 0 ∇u 0 = δ n -δ 0 on ∂D and (11) div(σ ε ∇u ε ) = 0 in D, ν • σ ε ∇u ε = δ n -δ 0 on ∂D,
respectively. The goal of this work is to find σ ε ≡ σ 0 such that M (σ ε ) vanishes.

Remark 2.1. In the framework of the PEM, the relative EIT measurements corresponding to the electrode locations x 0 , . . . , x N ∈ ∂D are, arguably, more intuitively described by the measurement map

(12) M (σ ε ) : I → (u ε I -u 0 I )(x n ) N n=0 , R N +1 ⋄ → R N +1 /R ≃ R N +1 ⋄ , where R N +1
⋄ denotes the mean-free subspace of R N +1 , and u ε I and u ε I are the solutions of (2) and (7), respectively, for the boundary current density

f = N n=0 I n δ n ∈ H -(d-1)/2-η ⋄ (∂D), η > 0.
In other words, M (σ ε )I ∈ R N +1 /R contains the relative potentials at the electrodes (uniquely defined up to the ground level of potential) when the net currents I ∈ R N +1 ⋄ are driven through the boundary 1 Here,

H 1 loc (D) stands for the space of distributions v ∈ D ′ (D) verifying v|ω ∈ H 1 (ω) for all domains ω ⋐ D. points x 0 , . . . , x N ∈ ∂D. It is easy to see that the matrix M (σ ε ) ∈ R N ×N is the representation of M (σ ε ) : R N +1 ⋄ → R N +1
⋄ with respect to the basis e n -e 0 , n = 1, . . . , N, where e n denotes the nth Cartesian basis vector of R N +1 . In particular, M (σ ε ) is the null matrix if and only if M (σ ε ) is the null operator, that is, if and only if σ ε cannot be distinguished from σ 0 by PEM measurements.

General scheme

Let Ω ⋐ D ⊂ R d be a nonempty Lipschitz domain. We recall the decomposition σ ε = σ 0 + εκ and require that κ ∈ L ∞ (D) satisfies supp(κ) ⊂ Ω. Our leading idea, originating from [START_REF] Bonnet-Ben Dhia | Transmission eigenvalues and far field invisibility for a finite number of incident/scattering directions[END_REF][START_REF] Bonnet-Ben Dhia | Obstacles in acoustic waveguides becoming "invisible" at given frequencies[END_REF][START_REF] Nazarov | Asymptotic expansions of eigenvalues in the continuous spectrum of a regularly perturbed quantum waveguide[END_REF], is to consider a small enough ε > 0 so that we can compute an asymptotic expansion for M (σ ε ). To this end, let us first write ( 13)

u ε n = u 0 n + εũ ε n where u ε n , u 0 n ∈ (H -(d-4)/2-η (D) ∩ H 1 loc (D)
)/R are defined by ( 11) and [START_REF] Cheng | Electrode models for electric current computed tomography[END_REF], respectively. Plugging ( 13) in [START_REF] Haberman | Uniqueness in Calderón's problem with Lipschitz conductivities[END_REF], we see that ũε n must satisfy the Neumann problem ( 14)

div(σ ε ∇ũ ε n ) = -div(κ∇u 0 n ) in D, ν • σ 0 ∇ũ ε = 0 on ∂D.
Since div(κ∇u 0 n ) defines a compactly supported source in H -1 (D), by applying the Lax-Milgram lemma to the variational formulation of ( 14), we see that ( 14) has a unique solution, i.e. ũε n , in H 1 (D)/R. 2 Moreover, the Lax-Milgram lemma and ( 5) also imply that [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators. I: Distribution Theory and Fourier Analysis[END_REF] ũε

n H 1 (D)/R ≤ C (σ ε ) -1 L ∞ (D) κ L ∞ (Ω) u 0 n H 1 (Ω)/R ≤ C (σ ε ) -1 L ∞ (D) κ L ∞ (Ω)
where the latter C > 0 depends on Ω.

Inserting ( 13) in ( 9), we deduce that ( 16)

M ij (σ ε ) = -ε D κ ∇u 0 i • ∇u 0 j dx -ε 2 D κ ∇ũ ε i • ∇u 0 j dx.
The representation [START_REF] Hyvönen | Approximating idealized boundary data of electric impedance tomography by electrode measurements[END_REF] demonstrates that M (σ ε ) is of the order ε, which was to be expected as M (σ 0 ) is the null matrix and σ ε -σ 0 = εκ. More interestingly, the first term in the asymptotic expansion of M (σ ε ) has a simple linear dependence on κ, which makes it relatively easy to find a nontrivial κ such that M (σ ε ) is of the order ε 2 . Unfortunately, the higher order terms in ε depend less explicitly on κ. To cope with this difficulty, we will next introduce a suitable fixed point problem. We redecompose κ in the form

(17) κ = κ 0 + N j=1 j i=1 τ ij κ ij
where τ ij ∈ R are free parameters. Moreover, κ 0 , κ ij ∈ L ∞ (D) are fixed functions that are supported in Ω and assumed to satisfy the conditions ( 18)

D κ 0 ∇u 0 i ′ • ∇u 0 j ′ dx = 0 for 1 ≤ i ′ ≤ j ′ ≤ N and ( 19
) D κ ij ∇u 0 i ′ • ∇u 0 j ′ dx = 1 if (i, j) = (i ′ , j ′ ) or (j, i) = (i ′ , j ′ ), 0 else.
The construction of such κ ij will be considered in Lemma 3.2 and Section 4 below, but meanwhile we just assume they exist. We remark that it seems reasonable to assume that the N (N + 1)/2 free parameters τ ij in the perturbation κ are enough to cancel out the symmetric matrix M (σ ε ) ∈ R N ×N (cf. Proposition 2.1). Substituting [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF] in [START_REF] Hyvönen | Approximating idealized boundary data of electric impedance tomography by electrode measurements[END_REF] and using ( 18)-( 19), we obtain the expansion

(20) M (σ ε ) = -ε τ -ε 2 M ε (τ ) where (21) M ε ij (τ ) = M ε ji (τ ) = D κ ∇ũ ε i • ∇u 0 j dx, 1 ≤ i ≤ j ≤ N,
and τ ∈ R N ×N is the symmetric matrix defined by the parameters τ ij .

From [START_REF] Kenig | The Calderón problem with partial data[END_REF] it is obvious that imposing M (σ ε ) = 0 in the chosen setting is equivalent to solving the following fixed point problem: [START_REF] Lassas | Inverse problems and invisibility cloaking for FEM models and resistor networks[END_REF] Find τ ∈ S N such that τ = F ε (τ ), where S N denotes the space of symmetric N × N matrices and F ε : S N → S N is defined by

(23) F ε (τ ) = -ε M ε (τ ).
Lemma 3.1 below ensures that for any fixed parameter γ > 0 and a small enough ε > 0, the map F ε has the invariant set (24)

B γ := {τ ∈ S N |τ | ≤ γ}
on which it is a contraction. In [START_REF] Nachman | Global uniqueness for a two-dimensional inverse boundary value problem[END_REF], | • | denotes an arbitrary norm of S N . In consequence, the Banach fixed point theorem guarantees the existence of ε 0 = ε 0 (γ) > 0 such that for all ε ∈ (0; ε 0 ], the fixed point problem ( 22) has a unique solution in B γ , enabling the construction of σ ε for which M (σ ε ) = 0 via ( 17) and ( 6).

It is important to notice that the constructed κ defined by ( 17) verifies κ ≡ 0 whenever κ 0 ≡ 0 by virtue of the orthogonality conditions ( 18)- [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. In other words, if κ 0 ≡ 0, then also σ ε ≡ σ 0 as required. Moreover, for given σ 0 , κ 0 , κ ij and γ > 0, the upper bound ε 0 > 0 can be tuned to ensure that the perturbed conductivity corresponding to κ of ( 17) satisfies σ ε ≥ c > 0 for all ε ∈ [0; ε 0 ] and τ ∈ B γ , which guarantees that all conductivities involved in the fixed point iteration are admissible.

In particular, we have proved the following result (modulo Lemma 3.1):

Proposition 3.1.
Let Ω ⋐ D be a Lipschitz domain. Assume that there are functions κ 0 , κ ij ∈ L ∞ (D) supported in Ω satisfying ( 18)- [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. Then, there exists a conductivity

σ ε ∈ L ∞ (D), with σ ε ≥ c > 0, such that σ ε -σ 0 ≡ 0, supp(σ ε -σ 0 ) ⊂ Ω and M ij (σ ε ) = δ i -δ 0 , (Λ ε -Λ 0 )(δ j -δ 0 ) ∂D = 0
for all i, j = 1, . . . , N . The following lemma shows that F ε is a contraction as required by the analysis preceding Proposition 3.1. Lemma 3.1. Let γ > 0 be a fixed parameter. Then, there exists ε 0 > 0 such that for all ε ∈ (0; ε 0 ], the map F ε is a contraction on the invariant set B γ defined by [START_REF] Nachman | Global uniqueness for a two-dimensional inverse boundary value problem[END_REF].

Proof. Assume γ > 0 is given. Let τ, τ ′ ∈ B γ be arbitrary, κ, κ ′ ∈ L ∞ (Ω) the corresponding perturbations defined by [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF], and ũε n , ũε n ′ ∈ H 1 (D)/R the associated solutions to [START_REF] Hanke | Convex backscattering support in electric impedance tomography[END_REF], with σ ε = σ 0 + εκ and σ ε = σ 0 + εκ ′ , respectively. By subtracting the equations defining ũε n and ũε n ′ and employing the short hand notation w = ũε n -ũε n ′ , we have

(25) div σ ε ∇w = -div (κ -κ ′ )∇(u 0 n + εũ ε n ′ ) in D, ν • σ 0 ∇w = 0 on ∂D.
From ( 5), ( 15), ( 6) and ( 17) it easily follows that

u 0 n + εũ ε n ′ H 1 (Ω)/R ≤ C uniformly for all ε ∈ [0; ε 0 ] if ε 0 > 0 is chosen small enough.
Subsequently, an application of the Lax-Milgram lemma to [START_REF] Nazarov | Asymptotic expansions of eigenvalues in the continuous spectrum of a regularly perturbed quantum waveguide[END_REF] results in the estimate

ũε n -ũε n ′ H 1 (Ω)/R ≤ C κ -κ ′ L ∞ (Ω) u 0 n + εũ ε n ′ H 1 (Ω)/R ≤ C κ -κ ′ L ∞ (Ω) ≤ C|τ -τ ′ | where all occurrences of C > 0 are independent of ε ∈ [0; ε 0 ]. In particular, M ε (τ ) -M ε (τ ′ ) ≤ C|τ -τ ′ |
by virtue of ( 21), the Cauchy-Schwarz inequality, ( 5) and [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators. I: Distribution Theory and Fourier Analysis[END_REF]. The definition [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] then implies that

(26) |F ε (τ ) -F ε (τ ′ )| ≤ C ε |τ -τ ′ | for all τ, τ ′ ∈ B γ ,
where 21), ( 5) and ( 15), we deduce from (26) that also 26) and ( 27) finally show that the map F ε is a contraction on the invariant set B γ for all ε ∈ [0; ε 0 ].

C > 0 is independent of ε ∈ [0; ε 0 ]. Moreover, noting that |F ε (0)| ≤ C ε due to (
(27) |F ε (τ )| ≤ C ε for all τ ∈ B γ . Reducing ε 0 if necessary, (
Let us denote by τ sol ∈ B γ the unique solution of the problem [START_REF] Lassas | Inverse problems and invisibility cloaking for FEM models and resistor networks[END_REF]. As a side product, the previous proof ensures the existence of a constant C 0 > 0 independent of ε ∈ (0; ε 0 ] such that (28)

|τ sol | = |F ε (τ sol )| ≤ C 0 ε for all ε ∈ (0; ε 0 ].
Combined with [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF], this tells us that κ is equal to κ 0 at the first order. At this stage, the remaining job consists in showing that there are functions κ 0 , κ ij supported in Ω satisfying ( 18)- [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. The following lemma is a classical result on dual basis (see, e.g., [21, Lemma 4.14]); we present its proof here because it will be utilized in the algorithm of Section 5.

Lemma 3.2.

There are functions κ 0 , κ ij ∈ L ∞ (D) supported in Ω and satisfying ( 18)-( 19) if and only if {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is a family of linearly independent functions on Ω. Proof. ⋆ Assume there are functions κ ij verifying [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. If α ij are real coefficients such that

N j=1 j i=1 α ij ∇u 0 i • ∇u 0 j = 0 in Ω,
then multiplying by κ i ′ j ′ and integrating over D, we find that α i ′ j ′ = 0. This allows to show that α ij = 0 for all 1 ≤ i ≤ j ≤ N and proves that the family {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is linearly independent. ⋆ Now, assume that {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is a family of linearly independent functions on Ω. To simplify the notation, we introduce the auxiliary functions

ψ k ∈ L 2 loc (D), k = 1, . . . , K = N (N +1)/2, such that ψ 1 = ∇u 0 1 • ∇u 0 1 , ψ 2 = ∇u 0 2 • ∇u 0 1 , ψ 3 = ∇u 0 2 • ∇u 0 2 , ψ 4 = ∇u 0 3 • ∇u 0 1 , ψ 5 = ∇u 0 3 • ∇u 0 2 , . . . ,
as well as the symmetric matrix A ∈ R K×K given elementwise by

A kk ′ = Ω ψ k ψ k ′ dx, k, k ′ = 1, . . . , K. If a = (α 1 , . . . , α K ) ⊤ ∈ ker A, defining v = K k=1 α k ψ k , we find that Ω v 2 dx = a ⊤ A a = 0.
We deduce that v = 0 in Ω, and thus α 1 = • • • = α K = 0 because, by assumption, the family {ψ k } 1≤k≤K is linearly independent. As a consequence, A is invertible. Consider the functions κ1 , . . . , κK defined via

(29) κk = K k ′ =1 A -1 kk ′ ψ k ′ | Ω in Ω and κk = 0 in D \ Ω.
In particular, κk ∈ L ∞ (D) because ψ k ′ are smooth in a neighborhood of Ω. We have

D κk ψ l dx = K k ′ =1 A -1 kk ′ A k ′ l = 1 if k = l, 0 else.
By renumbering κk , this proves the existence of functions κ ij verifying the conditions [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. Finally, we construct κ 0 ≡ 0 that satisfies (18):

(30)

κ 0 = κ # 0 - N j=1 j i=1 D κ # 0 κ ij dx κ ij ,
where κ # 0 is an arbitrary L ∞ -function such that supp(κ # 0 ) ⊂ Ω and κ # 0 / ∈ span{κ ij } 1≤i≤j≤N .

Construction of invisible conductivities in two dimensions

In this section, we study the two-dimensional case with σ 0 ≡ 1.

In order to complement Proposition 3.1, our goal is to demonstrate that there exist functions κ 0 , κ ij ∈ L ∞ (D) supported in a given Lipschitz domain Ω ⋐ D and satisfying ( 18)- [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. According to Lemma 3.2, this is equivalent to showing that {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is a family of linearly independent functions on Ω. Before setting to work, we remind the reader that the definition of u 0 n can be found in [START_REF] Cheng | Electrode models for electric current computed tomography[END_REF] and that x 0 , . . . , x N are N + 1 distinct points located on ∂D and corresponding to the positions of the electrodes. Proof. In this simple geometry, it is known that (see, e.g., [START_REF] Hakula | Sweep data of electrical impedance tomography[END_REF])

(31) u 0 n (x) = v n (x) -v 0 (x), x ∈ D, n = 1, . . . , N,
where v n is defined by

(32) v n (x) = - 1 π ln |x -x n |, x ∈ D. Let α ij ∈ R be such that (33) N j=1 j i=1 α ij ∇u 0 i • ∇u 0 j = 0 in Ω,
which can be written out explicitly:

(34) N j=1 j i=1 α ij x -x i |x -x i | 2 - x -x 0 |x -x 0 | 2 • x -x j |x -x j | 2 - x -x 0 |x -x 0 | 2 = 0 in Ω.
By analyticity, [START_REF] Sylvester | A global uniqueness theorem for an inverse boundary value problem[END_REF] holds in fact in all of R 2 \ E , where E := ∪ N n=0 {x n }. Multiplying (34) by |x -x j | 2 and letting x tend to x j , we see that α jj = 0 for all j = 1, . . . , N . Moreover, multiplying by |x -x 0 | 2 and letting x go to x 0 , we obtain the relation ( 35)

N j=2 j-1 i=1 α ij = 0.
Let us now introduce [START_REF] Uhlmann | Electrical impedance tomography and Calderón's problem[END_REF] w :=

N j=2 j-1 i=1 α ij u 0 i u 0 j = N j=2 j-1 i=1 α ij (v i v j -v 0 (v i + v j )) in R 2 \ E ,
where the second equality is a consequence of ( 31) and [START_REF] Taylor | Partial differential equations I: Basic theory[END_REF]. Using [START_REF] Sylvester | An anisotropic inverse boundary value problem[END_REF] and the fact that v n defined in [START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF] is harmonic in R 2 \ {x n }, we find that ∆w = 0 in R 2 \ E . Observing that v n is a multiple of the fundamental solution for the Laplacian centered at x n , a standard computation then yields ∆w = -2

N j=2 j-1 i=1 α ij (v j (x i ) -v 0 (x i )) δ i + (v i (x j ) -v 0 (x j )) δ j -(v i (x 0 ) + v j (x 0 )) δ 0 in R 2 .
Motivated by this expression, we introduce another auxiliary function, namely

(37) w = N j=2 j-1 i=1 α ij (v j (x i ) -v 0 (x i )) v i + (v i (x j ) -v 0 (x j )) v j -(v i (x 0 ) + v j (x 0 )) v 0 , which obviously satisfies ∆ w = ∆w in R 2 . Since v n belongs to C ∞ (R 2 \ {x n }), it is clear that w -w ∈ H 1-η loc (R 2
) for all η > 0. Furthermore, at infinity, it holds that w -w = o(|x|). This allows to prove that w -w is a tempered distribution (see, e.g., [15, Chapter VII]) on R 2 .

A classical extension of the Liouville theorem for harmonic functions ([35, Chapter 3, Proposition 4.6]) indicates that the only harmonic tempered distributions in R d are the harmonic polynomials; in particular, any harmonic tempered distribution that behaves as o(|x|) at infinity is a constant. This implies that w -w = C in R 2 for some C ∈ R. On the other hand, reordering the terms in ( 36) and [START_REF] Vauhkonen | Electrical impedance tomography with prior information[END_REF], we obtain the representation

w -w = v N N -1 i=1 α iN (v i -v i (x N )) -(v 0 -v 0 (x N )) + ŵN
where ŵN is a function that is analytic in a neighborhood of x N . We conclude that

v N N -1 i=1 α iN (v i -v i (x N )) -(v 0 -v 0 (x N ))
must also be analytic about x N . Because v N admits a logarithmic singularity at x N , this is possible if and only if

(38) N -1 i=1 α iN ((v i -v i (x N )) -(v 0 -v 0 (x N ))) = 0
in a neighborhood of x N . By analytic continuation, (38) must hold on all of R 2 \ E . Due to the singular behavior of the function v i at x i , this means that α iN = 0 for all i = 1, . . . , N -1.

Studying the behavior of w -w successively at x N -1 , . . . , x 2 , one shows analogously that α ij = 0 for all 1 ≤ i < j ≤ N -2. As the diagonal elements α jj , j = 1, . . . , N , were deduced to vanish already at the beginning of the proof, this shows that {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is a family of linearly independent functions on Ω. As ∂D is smooth, Φ also defines a smooth diffeomorphism of ∂D onto ∂B [START_REF] Pommerenke | Boundary Behaviour of Conformal Maps[END_REF]. Let û0 n ∈ H 1-η (B)/R, η > 0, be the unique solution of ∆û 0 n = 0 in B, ν • ∇û 0 n = δn -δ0 on ∂B, where δi ∈ H -1/2-η (∂B), i = 0, . . . , N , denotes the Dirac delta distribution supported at Φ(x n ) ∈ ∂B. According to [12, Proof of Theorem 3.2], it holds that (39)

u 0 n = û0 n • Φ in D for all n = 1, . . . , N . Let α ij ∈ R be such that N j=1 j i=1 α ij ∇u 0 i • ∇u 0 j = 0 in Ω.
Due to the harmonicity of u 0 n , n = 1, . . . , N , in D, we obtain

∆ N j=1 j i=1 α ij û0 i û0 j • Φ = ∆ N j=1 j i=1 α ij u 0 i u 0 j = 0 in Ω.
Since a composition with the conformal map Φ -1 retains harmonicity, it holds that

2 N j=1 j i=1 α ij ∇û 0 i • ∇û 0 j = ∆ N j=1 j i=1 α ij û0 i û0 j = 0 in Φ(Ω).
Because the family {∇û 0 i •∇û 0 j } 1≤i≤j≤N is linearly independent on Φ(Ω) by Proposition 4.1, it follows that α ij = 0 for all 1 ≤ i ≤ j ≤ N and the proof is complete.

Combining Proposition 3.1, Lemma 3.2 and Proposition 4.2, we obtain the main result of this article.

Theorem 4.1. Assume that σ 0 ≡ 1. Let D ⊂ R 2 be a simply connected and bounded domain with a C ∞ -boundary. For any compactly embedded Lipschitz domain Ω ⋐ D, there exists a conductivity

σ ε ∈ L ∞ (D), with σ ε ≥ c > 0, such that σ ε -σ 0 ≡ 0, supp(σ ε -σ 0 ) ⊂ Ω and M i,j (σ ε ) = δ i -δ 0 , (Λ ε -Λ 0 )(δ j -δ 0 ) ∂D = 0
for all i, j = 1, . . . , N . In other words,

M (σ ε )I = 0 for all I ∈ R N +1 ⋄ ,
where M (σ ε ) : R N +1 ⋄ → R N +1 /R is the relative PEM measurement map from [START_REF] Hakula | Sweep data of electrical impedance tomography[END_REF] corresponding to the electrode locations x 0 , . . . , x N ∈ ∂D. Remark 4.1. In [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF][START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF] it has been shown that M (σ ε ) = 0 if and only if σ ε = σ 0 (≡ 1) when the number of electrodes is countably infinite in two dimensions. Here, we have demonstrated that one can constructively find σ ε ≡ σ 0 such that M (σ ε ) = 0 in case the number of electrodes is finite.

Algorithmic implementation

For a given Lipschitz domain Ω such that Ω ⊂ D, functions κ 0 , κ ij satisfying ( 18)-( 19) can be constructed by following the line of reasoning in the second part of the proof of Lemma 3.2, assuming the potentials u 0 n verifying (10) are available. Notice that the bounded function κ # 0 with supp(κ # 0 ) ⊂ Ω appearing in (30) can be chosen arbitrarily as long as it does not belong to span{κ ij } 1≤i≤j≤N . As a consequence, finding κ # 0 is almost trivial and allows a lot of freedom. We denote by τ k , κ k , u ε,k n , ũε,k n the realizations of τ , κ, u ε n , ũε n at iteration k ≥ 0 of the algorithm; the aforementioned entities are introduced in ( 21), ( 17), ( 11), [START_REF] Hanke | Justification of point electrode models in electrical impedance tomography[END_REF], respectively. Moreover, we set σ ε,k = σ 0 + εκ k . Using formulas ( 21) and ( 23), we recursively define

τ k+1 ij = -ε D κ k ∇ũ ε,k i • ∇u 0 j dx, k ≥ 0, 1 ≤ i ≤ j ≤ N.
Since εũ ε,k n = u ε,k n -u 0 n by [START_REF] Hanke | Justification of point electrode models in electrical impedance tomography[END_REF], we obtain (40)

τ k+1 ij = τ k ij - D κ k ∇u ε,k i • ∇u 0 j dx, k ≥ 0, 1 ≤ i ≤ j ≤ N,
due to [START_REF] Hyvönen | Optimizing electrode positions in electrical impedance tomography[END_REF] and [START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF]. In particular, remark that

|τ k+1 ij -τ k ij | ≤ η implies M ij (σ ε,k
) ≤ ε η for any η > 0 by virtue of (9). To sum up, our algorithm for computing invisible conductivity perturbations is as follows: Algorithm 1. Assume that the potentials u 0 n , n = 1, . . . , N , are available. Construct {κ ij } 1≤i≤j≤N , choose κ # 0 / ∈ span{κ ij } 1≤i≤j≤N and compute κ 0 . Select ε > 0 and τ 0 ∈ S N . Run the fixed point iteration (40) until the desired stopping criterion is met. Construct the invisible conductivity perturbation via [START_REF] Hyvönen | Point measurements for a Neumann-to-Dirichlet map and the Calderón problem in the plane[END_REF]. If divergent behavior is observed, decrease ε.

Numerical experiments

In this section, we implement Algorithm 1 by resorting to finite element (FE) approximations of the PEM. We are interested in validating the convergence of the fixed point iteration (40) and visualizing the output conductivity σ ε (cf. ( 6)) for different subdomains Ω, 'initial guesses' κ # 0 and (point) electrode configurations. It turns out that especially the choice of κ # 0 has a considerable effect on the output, which means that Algorithm 1 can straightforwardly be used to construct several invisible conductivity perturbations for a given measurement configuration and Ω. The degree of indistinguishability of the produced σ ε compared to the background conductivity σ 0 is also tested in the framework of the more realistic CEM: According to our simulations, the relative CEM measurements corresponding to smallish electrodes and the constructed conductivities fall below any reasonable measurement noise level.

We only consider the homogeneous reference conductivity σ 0 ≡ 1 and choose D ⊂ R 2 to be the unit disk, meaning that the functions u 0 n satisfying [START_REF] Haberman | Uniqueness in Calderón's problem with Lipschitz conductivities[END_REF] are explicitly given by [START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF]. Due to the Riemann mapping theorem, this geometric simplification does not severely reduce the generality of The functions u ε,k n needed when running Algorithm 1 are computed at each iteration as follows: First, we (numerically) solve [START_REF] Hanke | Convex backscattering support in electric impedance tomography[END_REF] with σ ε = σ ε,k = 1 + εκ k and κ = κ k to obtain ũε,k n . We remind the reader that ũε,k n belongs to H 1 (D)/R so it can be approximated using standard FE methods. Then, we set u ε,k n = u 0 n + εũ ε,k n . When computing ũε,k n , we employ piecewise quadratic polynomial FE basis (Lagrange P 2 ); the number of triangular elements is around 20 000 in all tests. The conductivity perturbation 'shape functions' κ ij are evaluated (and interpolated) on the FE mesh with the help of ( 29) and [START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF]. Without exception, we use the stopping rule Table 1. The parameters ε and κ # 0 for (a)-(c) together with the respective relative CEM discrepancies E ε CEM (I) defined by (42). The CEM potentials are simulated with the 'trigonometric current basis', the universal electrode width π/32 and contact resistances of magnitude 0.01 [START_REF] Cheney | Electrical impedance tomography[END_REF][START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF].

for Algorithm 1, that is, the difference between consecutive iterates in the scheme (40) is monitored. The choice of the threshold value 10 -8 in (41) is a slight overkill since for practical EIT the noise level in relative measurements is typically over 1% [START_REF] Cheney | Electrical impedance tomography[END_REF], and (41) corresponds to a considerably smaller error for all (relative) PEM data we have simulated. As the starting value for the iteration (40), we choose the null matrix τ 0 = 0 in all experiments. The numerical results are presented in Figures 1234.

Figure 1 illustrates a first convergence test for Algorithm 1, investigating the effect of the choice for the free parameter ε. There are four electrodes and a concentric disk of radius 1/2 serves as the inclusion Ω. The parameter κ # 0 is set to κ # 0 ≡ 1. The convergence rate of the algorithm decreases as a function of ε, but the fixed point scheme remains convergent for all ε in the interval (0, 6] (and beyond).

Figure 2 illustrates that increasing the number of electrodes has a significant impact on the output conductivity σ ε of Algorithm 1. In particular, the spatial frequency of the angular oscillations in σ ε seems to be directly linked to the number of electrodes. Moreover, at least with the chosen simple inclusion shape, the deviations of σ ε from the unit background become less significant as the number of electrodes increases even if ε stays the same. The inclusion shape and the choice of the initial perturbation κ # 0 have considerable effects on the output of the algorithm as shown in Figure 3. A large Ω lying close to ∂D typically yields significant deviations from the unit background in σ ε , even with a relatively high number of electrodes. Moreover, asymmetries in κ # 0 lead to an asymmetric σ ε even in a geometrically symmetric setting. The indistinguishability of the output conductivities in Figure 3 from the unit background is evaluated using FE approximations of the CEM [START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF][START_REF] Vauhkonen | Electrical impedance tomography with prior information[END_REF] with about 10 000 quadratic triangular elements, a fixed electrode width of π/32 and an underlying contact resistance parameter equal to 0.01 on all electrodes. As a measure of the relative discrepancy for simulated (noiseless) CEM voltages, we employ (42) E ε CEM (I) = U ε (I) -U 0 (I) U 0 (I) .

Here, U ε (I), U 0 (I) ∈ R N (N +1) consist of the (stacked) CEM electrode potentials corresponding to the 'trigonometric current basis' I = {I (j) } 1≤j≤N ⊂ R N +1 ⋄ [START_REF] Cheney | Electrical impedance tomography[END_REF] and σ ε , σ 0 , respectively; cf., e.g., [START_REF] Hyvönen | Optimizing electrode positions in electrical impedance tomography[END_REF]. In practice, the relative noise level in EIT measurements is significantly above 0.1% (cf. [START_REF] Cheney | Electrical impedance tomography[END_REF]). Hence, Table 1 indicates that all conductivities shown in Figure 3 are practically indistinguishable from the unit background in the framework of the CEM with smallish electrodes at the depicted locations (cf. [START_REF] Hanke | Justification of point electrode models in electrical impedance tomography[END_REF]).

Finally, Figure 4 illustrates a single simulated relative CEM potential for the conductivity σ ε in Figure 3(a) -to be precise, both the interior electric potential and the electrode potentials are considered (cf. [START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF]). The electrode widths and contact resistances are the same as in Figure 4. It is noteworthy that both the Dirichlet boundary value of the relative interior electric potential and the relative electrode potentials do vanish according to a visual inspection, but the same does not hold for the relative electric potential in the interior of D. Remark 6.1. It is possible that the technique developed in this article could also be adapted to construct invisible conductivities directly for the CEM. However, instead of the functions u 0 n defined by [START_REF] Cheng | Electrode models for electric current computed tomography[END_REF], one would need to work with potentials corresponding to unit net currents between electrodes of finite size modeled by the CEM (cf. [START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF]). This would make proving the convergence of the fixed point iteration as well as implementing the numerical algorithm more technical.

Concluding remarks

We have introduced a fixed point scheme for computing invisible conductivity perturbations (of the two-dimensional unit conductivity) for a given electrode configuration in the framework of the PEM for EIT. Our numerical experiments demonstrate that the constructed conductivities are practically indistinguishable from the background also if the measurements are modeled by the more realistic CEM. In particular, take note that the conductivities shown in Figures 1234are as good solutions as the unit conductivity to the reconstruction problem of EIT with the considered sets of electrodes if no prior information on the behavior of the conductivity is available.

Remark 3 . 1 .

 31 If κ 0 , κ ij supported in Ω and satisfying (18)-[START_REF] Imanuvilov | The Calderón problem with partial data in two dimensions[END_REF] exist, Proposition 3.1 indicates that the conductivities σ ε and σ 0 are indistinguishable by EIT measurements with electrodes at x 0 , . . . , x N modeled by the PEM. In other words, the relative PEM measurement map M (σ ε ) defined by[START_REF] Hakula | Sweep data of electrical impedance tomography[END_REF] satisfies M (σ ε )I = 0 for all electrode current patterns I ∈ R N +1 ⋄ .

4. 1 .Proposition 4 . 1 .

 141 The case of disk. Assume that σ 0 ≡ 1 and let D ⊂ R 2 be the open unit disk. Then, {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is a family of linearly independent functions on any nonempty Lipschitz domain Ω ⋐ D.

4. 2 .Proposition 4 . 2 .

 242 General case in two dimensions. With the help of conformal mappings, Proposition 4.1 can be generalized to the case of arbitrary smooth two-dimensional domains. Assume that σ 0 ≡ 1. Let D ⊂ R 2 be a simply connected and bounded domain with a C ∞ -boundary. Then, {∇u 0 i • ∇u 0 j } 1≤i≤j≤N is a family of linearly independent functions on any nonempty Lipschitz domain Ω ⋐ D. Proof. Denote the open unit disk by B ⊂ R 2 and let Φ : D → B be a conformal map of D onto B.

Figure 1 .

 1 Figure 1. Top row: Outputs σ ε of Algorithm 1 with different values of the perturbation size parameter ε > 0. The point-like electrodes are located at the polar angles θ = 1 • , 91 • , 181 • , 271 • . Bottom left: The discrepancy (41) between consecutive iterates of Algorithm 1. Bottom right: A surface plot of σ ε corresponding to ε = 6.0. our (two-dimensional) numerical experiments: In any smooth simply connected domain D ⊂ R 2 the potentials u 0n employed in the construction of κ ij can be computed using (39) as long as a conformal mapping sending D onto the unit disk is available.The functions u ε,k n needed when running Algorithm 1 are computed at each iteration as follows: First, we (numerically) solve (14) with σ ε = σ ε,k = 1 + εκ k and κ = κ k to obtain ũε,k n . We remind the reader that ũε,k n belongs to H 1 (D)/R so it can be approximated using standard FE methods. Then, we set u ε,k n = u 0 n + εũ ε,k n . When computing ũε,k n , we employ piecewise quadratic polynomial FE basis (Lagrange P 2 ); the number of triangular elements is around 20 000 in all tests. The conductivity perturbation 'shape functions' κ ij are evaluated (and interpolated) on the FE mesh with the help of (29) and[START_REF] Seiskari | Point electrode problems in piecewise smooth plane domains[END_REF]. Without exception, we use the stopping rule

Figure 2 .

 2 Figure 2. Outputs of Algorithm 1 for N = 5, 7, 9, 11 (recall that N + 1 is the number of electrodes). The point-like electrodes are located at the polar angles θ j = 1 • + j N +1 360 • , j = 0, 1, . . . , N . The values of ε used in the computations are 14, 12, 8 and 8, respectively.

4 × 3 (

 43 10 -4 (b) 2.0 exp(-(x + 0.5) 2 -y 2 1.1 × 10 -

Figure 3 .

 3 Figure 3. (a)-(f): Perturbed conductivities produced by Algorithm 1 for the shown electrode configurations. Bottom left: A surface plot of the conductivity (a).Table 1:The parameters ε and κ # 0 for (a)-(c) together with the respective relative CEM discrepancies E ε CEM (I) defined by (42). The CEM potentials are simulated with the 'trigonometric current basis', the universal electrode width π/32 and contact resistances of magnitude 0.01[START_REF] Cheney | Electrical impedance tomography[END_REF][START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF].

10 Figure 4 .

 104 Figure 4. An example of CEM potentials corresponding to σ ε of Figure3(a); see[START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF]. Left: The Dirichlet boundary value of the absolute CEM interior potential (black line), the absolute electrode voltages (black squares), the Dirichlet boundary value of the relative CEM interior potential (red dashed line) and the relative electrode potentials (red squares). The electrode-wise net input currents are also displayed (blue circles). The electrode width and contact resistances are as in Figure3. Right: The relative CEM interior potential. Equipotential lines and electric field arrows are displayed.

Table 1 :

 1

Note that ∇ • L 2 (D) is equivalent to the standard quotient norm of H 1 (D)/R by the Poincaré inequality.