Different Couplings Between Precipitation and Deformation at the Same Site: A Case Study at Central Pyrenees (Spain)

Impact of hydrology‐associated deformation spans from the study of the rock hydrological properties and aquifer dynamics to the identification of geological disaster precursors. We analyze a comprehensive set of hydrology‐associated deformation signals, recorded by two underground 70 m long strainmeters (GAL16 and LAB780) operating under the Central Pyrenees (Spain) from the end of 2011 to the end of 2018. The two strainmeters measure extensional strain at the picostrain (10−12) level along two nearly orthogonal horizontal directions. Sign (compression or extension along any of the two directions), size, and time evolution of hydrology‐associated deformation vary considerably from season to season. We identify four classes of signals, related to summer rainfalls (producing compression along both GAL16 and LAB780), late‐spring and autumn rainfalls (producing extension along GAL16 and compression along LAB780), snow melting (producing opposite diurnal cycles along GAL16 and LAB780), and winter snowfalls (producing extension along both GAL16 and LAB780). Similar precipitations of the same class can give rise to differently sized signals; differences may be due to infiltration and/or surface load inhomogeneities as well as to complexities of the underground hydrogeological set up, which can lead to nonlinear, intrinsically time‐dependent, responses. Our modeling sheds light on the diverse effects associated with different environmental conditions. The variety of hydrology‐associated deformation we report is warning for the interpretation of all kinds of deformation data, as for their difficult modeling, the causative physical mechanisms, and the extreme care that is needed when removing hydrologically induced effects. Our insights are particularly useful in mountain areas.

Several reported deformation signals are linked to significant rainfalls and have been modeled through semi-analytical or numerical approaches, considering a flat half-space or a 3-D domain and assuming elastic or poroelastic rheology. Commonly analyzed processes are water load on the topographic surface and/or water accumulation in fractures and karst channels (e.g., Braitenberg et al., 2019;Devoti et al., 2015;Jacob et al., 2010;Lesparre et al., 2017). Others (e.g., Amoruso et al., 2014) investigated the role of different fractures (fracture network and main fractures/faults) in connecting deformation with the seasonal recharge/ discharge cycle of groundwater flow and its long-term changes, using elasto-plastic rheology. As for poroelastic rheology, for example, Barbour and Wyatt (2014) and Jahr et al. (2008) used linear poroelasticity to study the effects of fluid injection and extraction; Silverii et al. (2019) investigated regional-scale transient poroelastic deformation correlated to seasonal and multiyear hydrological variability.
It is very difficult to model hydrology-associated deformation in complex environments and literature does not include evidences of signatures of different processes (snowfalls, rainfalls, and snow melting) at the same site. In this study, we analyze a comprehensive set of hydrology-associated deformation signals, which were recorded by two underground geodetic laser interferometers (strainmeters) installed at about 350 m depth under the Central Pyrenees (Spain). Hydrology-associated deformation varies considerably from season to season, as for sign (compression or extension along one or both directions), size, and time evolution. Our modeling sheds light on the diverse effects associated with different environmental conditions.

Study Area
The study area is located in the Spanish side of the Central Pyrenees (Figures 1a and 1b). The Pyrenees Hercynian orogeny (ca. 300-250 Ma; chronostratigraphic classification is hereinafter from ICS, 2020) generated folds and thrusts (Vergés et al., 2002). The subsequent Alpine orogeny basically affected the Mesozoic-Cenozoic sedimentary cover (ca. 252-34 Ma) and poorly penetrated into the Paleozoic basement (ca. 541-252 Ma), which maintains its Hercynian structures. The underground infrastructure hosting the strainmeters is embedded in the Paleozoic basement of the local axial zone of the Pyrenees chain (black and griotte limestone) and, almost up to the surface, schist of the Pic Lariste series (late Devonian-middle Carboniferous, ca. 383-315 Ma;ESHYG, 1999;IGME, 2011;ITGE, 1989ITGE, , 1994Figures 1c and 1f). The strainmeters are close to two Hercynian faults ( Figure 1f): the Pico de Ladrones fault, which is about vertically dipping and E-W striking, also drawn in Figure 1c, and the Casa de la Cuca fault, which is covered by a dejection fan and is not drawn in Figure 1c (ESHYG, 1999;IGME, 2011;ITGE, 1989ITGE, , 1994. Several faults (sub-vertically dipping, NNE-SSW striking, and perpendicular to the WNW-ESE structural direction) start from the Paleozoic basement and affect the entire stratigraphic series (CHE, 2021b;ESHYG, 1999;IGME, 2011;ITGE, 1989ITGE, , 1994. The study area includes high reliefs which belong to the Axial Pyrenees and other parts characterized by glacial erosion (e.g., the Tortiellas glacial circus zone), river erosion (by the Aragón river and its tributaries), and karst action (originating several sinkholes; ESHYG, 1999;García-Ruiz et al., 2011;IGME, 2011;ITGE, 1994;Jeurissen, 1968). The infrastructure hosting the strainmeters is close to the Aragón river valley, which is about 1,500 m lower than surrounding heights ( Figure 1d).
As for the hydrogeological setting, the study area belongs to the WNW-ESE Ezcuarre-Peña Telera karst aquifer (Figure 1b). Middle Devonian-Lower Carboniferous unit is permeable by fissuring and to a lesser extent karstification, strongly compartmentalized and with local confinements. Its impermeable base level is defined by the Atxerite layers (Lower Devonian;CHE, 2021b). Groundwater flow in the Paleozoic limestone is driven by karst structures, mainly developed along fractures and bedding planes which strike N70°E-N90°E, Figure 1c (ESHYG, 1999).
Close to the strainmeters, permeability of the black and griotte limestone Paleozoic formations is due to cracking and karstification (ESHYG, 1999). Karst fractures are likely to follow the sub-vertical stratification , showing surface geology (Web Map Services https://mapas.igme.es/ and http://infoterre.brgm.fr/) and topography (EU-DEM v.1.1, map generated using GMT v. 5.4.5), respectively; panel (e) zooms the red rectangle inside panel (d) and gives geometry of the two strainmeters (to scale); panel (f) is the geological cross-section along the highway tunnel stretch shown (magenta line) in panels modified from ESHYG, 1999). Strainmeter location is given by red empty circles in panels (c and f). Tunnels (magenta line, highway; green line, railway) are explicitly labeled in panels (c-f). Stations N002, N003, P076, PZ32, A271, and 9198x are indicated by labeled red dots in panel (c). UTM WGS84 30N coordinates.

10.1029/2021WR031081
4 of 20 of the black and griotte limestone immediately South of the Casa de la Cuca fault and North of the Pico de Ladrones fault (Figure 1f), extend down to at least the tunnels depth, and strike around N70°E (ES-HYG, 1999;INECO, 2011), consistently with the local direction of the synclinal fold axis (ITGE, 1989(ITGE, , 1994. Water circulation within these fractures is probably very fast. Schist of the Pic Lariste series is impermeable and acts as an aquitard. Impermeable formations extend up to the surface for about 500 m south of the Pico de Ladrones fault (anticline fold). Where schist does not reach the surface, both karst fractures and groundwater inflows related to rain and snow have been reported as occurring in the black limestone overlying the schist formation, Figure 1f (ESHYG, 1999;INECO, 2011). A hydro-geological investigation of the highway tunnel was performed in 1996-1997; hydrographs of the Aragón river at the Rioseta intake, the Rioseta springs (1,400 m above sea level, asl in what follows, considered a clear example of groundwater behavior), and the highway and railway tunnel drainages exhibit huge peaks with short recession times after precipitation events as well as diurnal cycles during snow melting (ESHYG, 1999).
Glacial erosion in the Tortiellas glacial circus area exposed Cretaceous limestone and calcarenite strata, characterized by medium permeability values, and originated a catchment area feeding a stream which disappears into the Tortiellas sinkhole and lake; the lake forms during snowmelt periods, maybe because the drainage conduit is still, at least partially, blocked with ice. The Tortiellas sinkhole-lake is probably linked to NNE-SSW-striking fractures, Figure 1c (ESHYG, 1999;García-Ruiz et al., 2011;IGME, 2011;Jeurissen, 1968). Main drainage of Cretaceous and Paleogene units (ca. 65-23 Ma) shows karst behavior (very strong flow variations and short recession times) which is restricted to the dynamic zone above the regional discharge levels imposed by the rivers. Vadose zones are several hundred meters thick; below, the saturated zone shows a rapid loss of permeability due to a decrease in karstification and sealing of discontinuities (CHE, 2021b). SW of the strainmeters, groundwater flow is SSW-directed and driven by the Lecherines karst structures, which are mainly located along the NNE-SSW-striking fractures cutting the Paleogene limestones, Figures 1c and 1f (IGME, 2011;ITGE, 1989;Sevil, 2019). Most sinkholes are at about 2,000 m asl. The groundwater flow area extends for about 5.5 km from the Tortiellas sinkhole lake through the Lecherines karst system and beyond (CHE, 2021a;ITGE, 1994), and includes a hydrologically active main fault, which strikes around N21°E and dips sub-vertically down to a depth of approximately 1,000 m (Izquierdo Llavall et al., 2013), and other sub-vertical faults. The main fault, the other nearby faults, and the karst structures act as an outflow route for rainwater; drainage is at around 950 m altitude near the Esjamundo Cave (about 400 m west of piezometer PZ32, Figures 1b and 1c) close to the Villanúa village (Sevil, 2019).
As for meteo-climatic parameters, temperature and precipitation are very heterogeneous because of the Pyrenees morphology and proximity to the Atlantic Ocean to the West and the Mediterranean Sea to the East. The western part (extending eastward as far as the headwaters of the Gállego and Ara rivers, several kilometers east of the study area) is dominated by the Atlantic Ocean weather conditions and most precipitations occur between December and March, whereas in the eastern part most precipitations occur in spring and autumn (April-June, September-November; e.g., López-Moreno et al., 2009). Snow mostly accumulates during winter in the western part (where heavy snowfalls most likely occur) and during late autumn, winter, and spring in the eastern part (Navarro-Serrano & López-Moreno, 2017). Precipitations occur as snow at altitudes higher than 1,500-1,600 m over the entire mountain range from late autumn onwards (López-Moreno & García-Ruiz, 2004). Snow persists above 1,600 m during the whole cold period and between 1,600 and 1,300 m during the coldest months (usually December-February; López-Moreno et al., 2011). Snowfall spatial distribution is not clearly related to elevation and geographic position (Navarro-Serrano & López-Moreno, 2017).

Underground Strain
Two high-sensitivity geodetic laser interferometers (strainmeters), hereinafter called LAB780 and GAL16, have been monitoring extension along two almost-perpendicular azimuthal directions (around N31°W and N77°E respectively, Figure 1e) at a depth of 350 m below the Earth surface from the end of 2011 to the end of 2018. GAL16 is in one of the bypasses that join a highway tunnel and a decommissioned railway tunnel connecting Spain to France; LAB780 is inside two narrow side halls parallel to and built at the same time as the railway tunnel ( Figure 1e). Both interferometers operate in the unequal-arm Michelson configuration, as detailed in Amoruso and Crescentini (2009) and Amoruso et al. (2018). They record extensional strain ΔL/L, where L is the strainmeter length and ΔL its change over time. We give strain in units of 10 −9 ; since each strainmeter is 70 m long, ΔL/L = 10 −9 is equivalent to an increase of its length by 0.07 μm. Nominal sensitivity is a few 10 −13 , frequency band ranges from dc up to hundreds of Hz, recording rate was 600 Hz until November 2017 and 1,000 Hz after then. The instruments became operational in winter 2011, but data acquisition suffered several interruptions (empty parts of the top plots in Figures 2a-2c) mainly caused by laser tube replacements. A long break occurred from December 2015 to December 2016.
We apply a low-pass anti-aliasing filter before decimating strain data at 600 s sampling time. Although strain noise amplitude is inversely proportional to frequency, we do not filter low frequencies out, but for easier visualization, we detrend strain data whenever we plot them. On time-scales longer than tens of seconds and shorter than several days, the largest signal on a strain record is the Earth (body and ocean loading) tides and occasionally a large earthquake or a hydrology-associated deformation event (e.g., Agnew, 2007); loading by nonlinear shallow-water tides in the Bay of Biscay also gives non-negligible deformation in the study area (Amoruso & Crescentini, 2016). Therefore, we subtract linear and nonlinear tides from the strain time series before further analysis (Amoruso & Crescentini, 2020).

Other Data
We also examine meteorological and hydrological data from various stations ( Figure 1c). N002 (around 2,080 m asl) and N003 (around 1,970 m asl) are equipped with snow, temperature, and pressure gauges; P076 and 9198x with pluviometers and temperature gauges; PZ32 with a piezometer (the only one continuously monitoring the Ezcuarre-Peña Telera aquifer) whose steel casing reaches a depth of 122 m, crossing calcarenite (from the surface down to about 40 m) and limestone (CHE, 2021a); A271 with a pluviometer as well as temperature and stream gauges. A small reservoir (Embalse de Canfranc) is located along the Aragón river; the dam is about 2.5 km South of the strainmeters (Figure 1c). We could not obtain data related to the reservoir storage change, however, we will see that dam operation strongly affects A271 stream data, but does not affect any other record visibly. The most indicative time series are shown in Figure 2; unfortunately, not all the records cover the whole investigated period. We use satellite optical images (Landsat 8 and Sentinel 2) to check snow extent. Figure 2 shows the scalogram plot obtained from the continuous wavelet transform (CWT) analysis of strain data. Colors represent the squared magnitude of the wavelet coefficients (i.e., the power spectrum of strain) in arbitrary units, logarithmic scale, as a function of time and frequency. The scalogram is plotted in the frequency range from 0.5 to 4 cycles per day (cpd), where energy related to earthquakes and other high-frequency sources is negligible and the instrumental noise level is exceptionally low (Amoruso & Crescentini, 2016, 2020.

Observations
Figure 2 also includes precipitations (P076), temperature (N003), snow thickness (N002 and N003), and water table level (PZ32). Unfortunately, all the available time series exhibit occasional breaks. Temperature and precipitations at the other stations are not shown for clarity. Precipitations (including rain, sleet, snow, and hail) occur during the whole year; summer is the driest season, with occasional heavy rains. Water table level at PZ32 is usually low (around 952 m asl) in summer, but it increases abruptly during heavy rains (up to around 957 m asl, close to the ground level) and then returns to its previous value in few days. Water table raises at the beginnings of autumn, usually remains high in winter, when snow accumulates, and is often particularly high (around 957 m asl) in spring, during snow melting. Peaks concurrent with precipitations are visible throughout the whole year.
Water level surface elevation (i.e., stage) of the Aragón river (A271) is not shown in Figure 2 because sometimes it is visibly affected by the Canfranc dam operation. Barometric pressure is not shown because on the investigated time scales it does not evidence anomalies during the occurrence of apparent large deformation signals. More generally, we could never find any clear, significant effect of barometric pressure on LAB780 and GAL16 data (Amoruso et al., 2018). Neither GAL16 nor LAB780 temperature is shown because yearly cycles are around 2°C and 0.7°C respectively; diurnal cycles are smaller than 0.05°C in GAL16 and below Comparison between scalograms of strain and time series of precipitation, temperature, snow thickness, and water table level. Panels (a-c) refer to 2012-2013, 2014-2015, and 2017-2018. From top to bottom, each panel shows: occurrence of the clearest hydrology-associated deformation events (dark green vertical arrows, snowfalls; red vertical arrows, late spring and autumn rainfalls; black vertical arrows, summer rainfalls) and the clearest periods of diurnal strain cycles due to snow melting (horizontal cyan segments); scalograms (colors give power spectrum on logarithmic scale in arbitrary units; frequency is plotted on logarithmic scale) of available strain data (LAB780 and GAL16) computed by using Matlab® v. R2020b (Wavelet Toolbox, CWT function, Morlet wavelet); daily precipitation at P076; daily mean temperature at N003; snow thickness at N002 (magenta) and N003 (violet); water table level at PZ32. Dashed blue vertical lines connect the clearest broad-frequency peaks in the scalogram to related precipitation events. Meteorological and hydrological data courteously provided by Confederación Hidrográfica del Ebro (CH Ebro, http://www.chebro.es) and Agencia Estatal de Meteorología, Gobierno de España (AEMET, http:// www.aemet.es). the detection limit in LAB780 (Amoruso et al., 2018). Thus, thermal strain is expected to be negligible on the same time scales as hydrology-associated strain events discussed here (hours to very few days). This conclusion is supported by strain tidal analysis; amplitude of measured S1 strain tide (the harmonics most contaminated by temperature, 1 cpd) is few 10 −10 , thus comparable to although somewhat larger than computed genuine S1 strain tide (Amoruso & Crescentini, 2020).
All the clearest broad-frequency peaks in the scalogram coincide with precipitation events (see dashed blue vertical lines in Figure 2), except the LAB780 peak on February 27-28, 2013 and the GAL16 peak on July 7, 2015; all the apparent peaks around 1 cpd are coeval with periods of snow thickness decrease. Thus, we attribute scalogram peaks to hydrology-associated deformation.
It is practically impossible to show plots of more than five years of strain and precipitation data at a time scale which allows to find out all concurrent events. We carefully examined all data, found concurrent events which are too faint or too slow to produce evident peaks in the scalogram, and finally recognized four classes of hydrology-associated deformation signals with different polarities (compression or extension) and duration. Occurrence of selected hydrology-associated signals is marked by downward vertical arrows and horizontal line segments on the top of Figures 2a-2c. Some of them do not exhibit clear peaks in the scalogram.
The four classes are characterized by: (a) extension along GAL16 and compression along LAB780 in coincidence with rainfalls in late spring and autumn (downward red arrows in Figure 2, example in Figure 3a); (b) opposite diurnal cycles along GAL16 and LAB780 during periods of snow melting in spring and after the first snowfalls in autumn (horizontal cyan segments in Figure 2, example in Figure 3b); (c) compression along both GAL16 and LAB780 in coincidence with rainfalls in summer (downward black arrows in Figure 2, example in Figure 3c); and (d) extension along both GAL16 and LAB780 in coincidence with snowfalls in winter (downward green arrows in Figure 2, example in Figure 3d). Since the two strainmeters are nearly orthogonal to each other, summing their strain gives a reasonable approximation of areal strain. Thus, winter snowfalls cause areal dilatation and summer rainfalls cause areal compression, whereas snow melting and late-spring/autumn rainfalls cause relatively small areal strain.
The main features of typical events belonging to the four classes of strain signals are described in the following subsections. Of course, also hybrid signals caused by coeval events, e.g., rainfall over low-altitude areas and snowfall over high(mid)-altitude ones, are sometimes visible in the strain records.

Late-Spring and Autumn Rainfalls
Red arrows on the top of Figures 2a-2c mark the occurrence of the clearest signals associated with latespring and autumn rainfalls in the strain records. Rainfalls cause extension on GAL16 and contraction on LAB780 when high-altitude areas are snowy, as confirmed by satellite images (see Figure S1 as for the events in Figure 3a). Sometimes it is difficult to estimate the duration of the strain recovery phase robustly; however, when the recovery phase is clear, it lasts about 0.5-1 days (Figure 3a). The water table rises, but in spring the increase is not large because the water table is permanently high. Effects of the Canfranc dam operation on the Aragón river stage are evident in Figure 3a, but there is no visible effect on either strain or water table level at PZ32. Lack of visible effects on strain is consistent with computations (analytical Boussinesq approach, rock rigidity 6 × 10 9 Pa; Terzaghi, 1943), which show that 1 m change of the reservoir level causes as small strain as few 10 −10 at the strainmeter site. The Aragón river stage shows clear transient signals when unaffected by the Canfranc dam activity, as shown in Figure S2; a related satellite image showing the snow line altitude is in Figure S3. Coeval diurnal cycles make it hard to estimate recovery time of the water table level (Figures 3a and S2); nevertheless, it seems longer than the strain one, which in turn is approximately as long as recovery time of the Aragón river stage ( Figure S2).

Snow Melting
Cyan horizontal segments on the top of Figures 2a-2c mark periods when diurnal cycles are clearly visible in strain records. They last several days and are frequently recorded during snow melting in autumn and mainly in spring, when temperature increases almost steadily. For most of those periods, excess energy around 1 cpd is clearly visible on the scalogram of strain. During postmeridian hours (noon to 6 p.m. typically) GAL16 measures extension, whereas LAB780 measures contraction; the opposite occurs during the rest of the day (Figure 3b). Deformation signals related to rainfalls sometimes overlap the diurnal cycles. Satellite images suggest that diurnal strain cycles usually occur when snow line altitude is between 1,500 and 2,000 m (see Figure S4 as for the snow melting period shown in Figure 3b). Water table level at PZ32 shows similar diurnal cycles, few tens of centimeters in amplitude; snow thickness decreases like a staircase function (Figure 3b). Effects of the Canfranc dam activity on the Aragón river stage are evident in Figure 3b, but, as in Figure 3a, there is no visible effect on either strain or water table level at PZ32. Figure S5 shows an example of apparent diurnal cycles in the Aragón river stage during snow melting, with no visible anomaly related to the Canfranc dam activity; a satellite image showing the snow line altitude is in Figure S6. Aragón river stage increases when LAB780 strain decreases (compression) and GAL16 strain increases (extension); water table level is shifted by a couple of hours. Diurnal cycles in the Aragón river stage and in the PZ32 water table level overlap an overall increase in their time series (Figure 3b and Figure S5).

Summer Rainfalls
Occurrence of the clearest strain signals related to summer rainfalls is indicated by black arrows on the top of Figures 2a-2c. CWT time-frequency analysis does not produce clear signatures in the scalogram, because signals are very slow. Summer rainfalls are infrequent and often occur as strong showers. Unfortunately, most of them caused power failures and consequent breaks in the strain acquisition. Available records show that signals related to rainfalls are compressive for both GAL16 and LAB780 when either the whole study area is free from snow (i.e., in summer, Figure 3c) or only the highest reliefs are covered by snow. The Aragón river stage evidences the effects of the Canfranc dam activity in Figure 3c; as in Figures 3a and 3b, such activity does not have any visible effect on either strain or water table level at PZ32. The Aragón river stage shows large transient signals when unaffected by the Canfranc dam activity, as shown in Figure S7; a satellite image showing the snow line altitude is in Figure S8. Duration of the strain recovery phase is often difficult to assess; however, it seems longer than for late-spring rainfalls. This also happens for water table level (Figures 3c and S7) and Aragón stage ( Figure S7).

Winter Snowfalls
Green arrows on the top of Figures 2a-2c mark the occurrence of the clearest strain signals associated with winter snowfalls. Excess energy is clearly visible on the scalogram of strain. In winter, it often snows down to low altitudes; when heavy snowfalls occur, both strainmeters record as large extensive strain as several 10 −7 . Signals manifest as transients, with steep rise and slower recovery, usually approximately as long as for strain associated with late-spring/autumn rainfalls ( Figure 3d); sometimes final strain is compressive. Snow thickness at N002 and N003 may increase by up to tens of cm, sometimes asynchronously. The Aragón river stage rises by tens of cm. Water table level at PZ32 usually reaches its maximum value (around 957 m asl) regardless of its starting value; recovery may be incomplete and is somewhat slower than recovery of stage, which in turn is approximately as long as recovery of strain ( Figure 3d). Unfortunately, we did not find any cloudless satellite image for the signal in Figure 3d; however, Figure S9 shows snow extent the day after the snowfall discussed in Section 6.2.2.2.

Proposed Driving Processes
The class-by-class different responses of GAL16 and LAB780 indicate that observed deformation is generated by different processes occurring under different environmental conditions. In the following subsections, we propose a simplified view of the processes that might drive observed hydrology-associated deformation.

Deformation Associated With Rainfalls
Because of the high permeability of the most local rock, deformation related to rainfalls is probably due to water infiltration. The two strainmeters are located inside a W-E anticline at the bottom of a schist stratum, which is about 400 m thick (Figures 1c and 1f) and has low hydraulic conductivity. Erosion has exposed the schist anticline in the strainmeter area and few kms eastward, but its W-E extension is probably continuous and longer. Limestone has higher hydraulic conductivity and shows karstic features, mainly at the borders of the limestone strata (Section 2). Rainwater moves fast into limestone, causing a sudden increase of the local water table level, as confirmed by the abundant groundwater inflows observed during the railway tunnel excavation concurrently with heavy rain episodes; such a sudden increase does not occur in the schist layer, which acts as an aquitard (ESHYG, 1999;INECO, 2011).
We think that difference between schist and limestone water table levels, as well as limestone fracture dilatation, cause compressional forces (pressure) on the southern and northern faces of the schist unit, which is approximately W-E oriented (Figure 1c), until the water table level returns to its pre-precipitation value everywhere. Compressional forces on the schist unit occur where limestone (or other overlying permeable rock) reaches the surface and is not covered by snow, i.e., where infiltration is easy. Thus, the strainmeters are well inside the potentially SN-compressed area in summer and close to its western end in spring and autumn (Figures 1c and 4a, and satellite images in Figures S1, S3, S4, S6, S8, and S8). Also, the Lecherines karst and the Tortiella sinkhole (Figure 1c) might play an important active role in producing deformation at the strainmeters when karst conduits are not blocked with ice (approximately June-October, Bartolomé, 2019).

Deformation Associated With Snow Melting
Also deformation recorded during snow melting is probably due to water infiltration. We do not ascribe diurnal strain cycles to surface unloading because (a) surface loading (unloading) causes extension (compression) on both LAB780 and GAL16 (Section 5.3), thus the two strainmeters should record in-phase diurnal strain cycles, whereas they are in phase opposition (GAL16 exhibits extension when LAB780 exhibits  Figures S4 and S6) and possible unloading at altitudes higher than about 1,800 m asl seems to be ineffectual in producing either deformation at the two strainmeters or hydrological effects (Section 5.3).
Snow line altitude is approximately between 1,500 m and 2,000 m asl also when late-spring/autumn rainfalls occur (Figures S1 and S3) and recorded strain is negative (compressive) for LAB780 and positive (extensive) for GAL16. Analogously, diurnal strain cycles exhibit LAB780 compression and GAL16 extension in the afternoon, when infiltration because of snowmelt is expected to be large. Thus, we attribute observed strain signals during snow melting to the same infiltration process as late-spring/autumn rainfalls (Section 5.1). Recovery time of deformation associated with late-spring/autumn rainfalls is consistent with the observation of diurnal cycles, because they would be filtered out if the characteristic time were much longer.

Deformation Associated With Snowfalls
When snowfalls occur at low altitude (down to about 1,100-1,200 m asl, see Figure S9), both strainmeters record extension (Section 4.4); when snow line altitude is approximately between 1,500 m and 2,000 m asl, precipitation causes extension along GAL16 and compression along LAB780 (Section 4.1). Thus, precipitation causes different strain responses in spring/autumn and winter because of the different processes occurring at altitudes approximately between 1,200 and 1,800 m. Infiltration is expected to be small because of snow accumulation also at low altitude, and, in any case, infiltration would cause compression along LAB780 (Section 5.1). We think that recorded deformation is due to snow load on the surface; however, rock response to the load is not elastic, otherwise, the strain time history would follow the snow thickness (or, better, the snow water equivalent SWE) time history. Durations of the recovery phase of the strain signal and of the recession limb of the Aragón river stage hydrograph are quite similar (e.g., Figure 3d) thus suggesting that both of them may be driven by the same process; surface load induces uneven rock volumetric deformation, which in turn causes groundwater migration within the rock and to the outside, through springs, until a new equilibrium is reached. Groundwater migration affects strain and river stage.
Detailing this process in a fractured medium is a very difficult task; however, in the simpler case of a poroelastic medium, a sudden surface load induces an initial instantaneous undrained elastic response of saturated rock; afterward deformation changes over time predictably (poroelastic evolution), and at long times it is given by a drained elastic response. The two elastic responses are characterized by the same rock rigidity (shear modulus) but different Poisson ratios ν u and ν d , with ν u > ν d ; the time scale for the completion of the process depends on the solid and fluid material properties, as well as on the size and geometry of a specific case (Wang, 2000). Local snow-load-induced deformation depends on depth to the surface (thus, the interferometer asl altitude), nearby topography, and load extension and location.
When it snows between 1,200 m and 1,800 m asl, it may also snow at higher altitudes. For example, snow thickness increased by about 70 cm at N002 and 40 cm at N003 during the time period shown in Figure 3d. About 25% of the snow thickness increase occurred on December 26, 2013, around noon, and a concurrent abrupt increase of snow thickness (not shown in Figure 3d for clarity) was measured at stations N001 (about 23 km NW of the strainmeters) and N004 (about 24 km ENE of the strainmeters). Thus, a heavy snowfall occurred over a large area on December 26, 2013, around noon, but no concomitant signal is visible as for deformation at the two strainmeters, water table level at PZ32, river stage at A271, precipitation at P076. An even clearer example occurred on February 6, 2013-when a large increase of snow thickness was recorded at N001 (60 cm), N002 (60 cm), N003 (60 cm), and N004 (90 cm) without precipitation at low altitudes in (1,180 m asl) for 10 kPa test pressure. Panels (b and c), strain along LAB780 and GAL16, respectively, during spring and autumn because of groundwater push on the northern and southern faces of the parallelepiped. Panels (d and e), strain along LAB780 and GAL16, respectively, during summer because of groundwater push on the northern and southern faces of the parallelepiped. Panels (f and g), strain along LAB780 and GAL16, respectively, during summer because of groundwater pressure inside the effective N21E fault. Red lines in panels (b-g) give strainmeter (LAB780 or GAL16) geometry to scale; red numbers give computed strain (10 −9 units) at the strainmeter centers. Rock rigidity is 6 × 10 9 Pa, Poisson ratio is 0.33. Note that the color scale is different for strain along LAB780 and GAL16. the study area (stations P076 and A271)-but again no concomitant signal is visible as for deformation at the two strainmeters, water table level at PZ32, Aragón stage at A271 ( Figure S11). Consequently, snow load acting over about 1,800 m asl seems ineffectual in producing either deformation at the two strainmeters or hydrological effects (water table level at PZ32, Aragón river stage at A271).

Modeling
Realistic quantitative modeling of the observed hydrology-associated deformation signals is quite hard because, for example, the study area lacks spatio-temporal groundwater observations and adequate knowledge of underground properties, precipitation is unevenly distributed, karstic features are very complex. In the following subsections, we propose effective simplified models.

Methods
We model deformation using finite element method software COMSOL Multiphysics® v. 5.2. Geometry of the computational domain depends on the physical process under investigation, namely elastic deformation or poroelastic strain evolution. In both cases, we fix normal displacements to zero on the outermost lateral boundaries and on the bottom; ground surface is under stress free condition where it is not loaded by snow.
We compute elastic deformation in order to model strain induced by infiltration and undrained/drained strain (poroelastic evolution endpoints) induced by surface load. The computational domain has a squared surface projection (20 × 20 km 2 ) and extends down to 4,500 m below sea level (see sketch in the Figure S10); thus, it is much larger than the study area shown in Figure 1c. Topography is from the European Union Digital Elevation Model (EU-DEM) v. 1.1.
When modeling infiltration-induced strain (Sections 5.1 and 5.2), we approximate schist subjected to compressional forces from the karst aquifer as a WE-oriented rectangular parallelepiped, 500 m wide and extending from the surface down to 1,000 m asl ( Figure S10); length is different in summer and spring/ autumn because of the different extent of the infiltration area. Surface expression of the northern and southern walls of the parallelepiped is indicated in Figure 4a. Uniform test pressure (10 kPa) is applied on the northern and southern walls of the schist parallelepiped ( Figure 4a) from 1,000 m asl up to 1,500 m asl in spring/autumn and 1,300 m asl in summer, because difference between minimum discharge level and upper overflow at the Rioseta springs reaches almost 100 m in elevation (CHE, 2021b) and groundwater is continuously drained from the tunnels, but flows exhibit large seasonality (ESHYG, 1999). In summer and early autumn, Lecherines karst could also contribute to deformation at the interferometers, because rainfalls increase water pressure inside faults and karst channels. We model Lecherines karst effects through deformation caused by applying uniform test pressure (10 kPa) to the walls of an effective vertical N21E fault (Figure 4a).
The karst aquifer is probably very heterogeneous, but the schist unit may be considered approximately homogeneous; we expect mechanical anisotropy to be of little influence on strainmeter data because compressional forces act horizontally and schist layering is almost horizontal (Figure 1f). Moreover, as for larger-scale phenomena, tidal strain (diurnal band to at least 8 cycles per day) recorded by the two strainmeters does not evidence any anomaly attributable to mechanical heterogeneity or anisotropy (Amoruso & Crescentini, 2016). Thus, for simplicity, we assume homogeneous isotropic conditions in the whole computational domain.
As for undrained and drained snow-load-induced strain, ground surface is subject to test load (1 m SWE) between 1,200 and 1,800 m and stress-free elsewhere (Section 5.3). No pressure is applied on either the southern and northern faces of the parallelepiped or the walls of the N21E fault. Again, we assume homogeneous isotropic conditions in the whole computational domain.
As previously mentioned, undrained/drained elastic deformation represents the endpoints of poroelastic evolution. We model snowfall-associated deformation (Section 5.3) following a three-step procedure. As first step, we test whether and to what extent topography, depth to the surface, and position with respect to the surface load affect strain poroelastic evolution after an impulsive snowfall with no subsequent melting (sudden surface load). We consider a layered poroelastic half-space (a homogeneous aquifer overlying a homogeneous aquitard) with a simple axisymmetric topography (Figure 5a). The outermost lateral boundaries and the bottom of the computational domain also act as confining surfaces. The uppermost area of the computational domain is subjected to a sudden surface load (impulsive snowfall with no subsequent melting) which is 0 before the snowfall (time t < 0) and 9.8 kPa (i.e., we test 1 m of SWE) after the snowfall (t ≥ 0). Initial hydraulic head with respect to the plane surface (z = 0 in Figure 5a) is null. The five computation points mimic strainmeter location. Since in the first step, we obtain a quasi-exponential strain evolution from the initial elastic undrained response to the final elastic drained response, as second step, we check if an exponential strain evolution caused by an impulsive snowfall can reproduce the observed strain signal when considering the snowfall rate history of a real event. If a real finite-duration snowfall is considered, the linear relationship between stress and strain implies that discrete-time strain value y i depends on both the instantaneous values of snowfall rate and snow accumulation (u i−n and v i−n respectively, where n is a possible delay) and also on the previous value of y, through a simple linear input-output polynomial model with three parameters (mathematical details are in Text S1), i.e., y i = ay i−1 + bu i−n + cv i−n where a is linked to the characteristic time of the exponential transition, b to the instantaneous undrained response, and c to the equilibrium drained response. We optimize parameters a, b, c, and n through least squares fitting to various clear snowfall-associated strain signals, using Matlab® v. R2020b (System Identification Toolbox, ARX function) and assuming that snow melting is negligible meanwhile. Lastly, as third step, we compute , axisymmetric computational domain, the dotted-dashed magenta line at r = 0 giving the symmetry axis; superficial load (9.8 kPa, starting at t = 0) is restricted to the uppermost area, as indicated by the violet vertical arrows, whereas the remaining ground surface is under stress free condition; normal displacements are zero on the outermost lateral boundaries and on the bottom; initial hydraulic head with respect to z = 0 is null. Panel (b and c), time history of tangential (radial) strain for the five test points close to the slope shown in panel (a), same colors. Region A (aquifer), undrained Poisson ratio, 0.33; drained Poisson ratio, 0.26; rigidity, 6 × 10 9 Pa; hydraulic diffusivity, 20 m 2 s −1 . Region B (aquitard), undrained and drained Poisson ratio, 0.25; rigidity, 8 × 10 9 Pa; hydraulic diffusivity, 0.13 m 2 s −1 .
snow-load-induced undrained/drained elastic deformation at the interferometers using real topography of the study area.

Deformation Associated With Infiltration (Rainfalls and Snow Melting)
When applying pressure on the northern and southern walls of the schist parallelepiped (Section 6.1), computed strain is negative (compressive) for LAB780 and positive (extensive) for GAL16 in spring and autumn (Figures 4b and 4c); it is negative (compressive) for both strainmeters in summer (Figures 4d and 4e). Several tests, not shown here for brevity, indicate that strain sign (negative/positive) is not affected by different choices of pre-precipitation water table level nearby the strainmeters or the parallelepiped dimensions, but numbers may change.
Although sign of computed deformation agrees with observations, LAB780-to-GAL16 amplitude ratio of recorded signals is usually lower than model values (i.e., computed GAL16 strain is too small in amplitude with respect to computed LAB780 strain) especially in summer (e.g., Figures 3a-3c, S2, S5, and S7) when modeled GAL16 deformation is very small. This discrepancy is probably due to the approximations in our model; however, water pressure increase in the Lecherines karst causes negligible LAB780 deformation and negative (compressive) GAL16 deformation (Figures 4f and 4g); the latter is larger than and adds to compression produced by pressure acting on the southern and northern faces of the schist parallelepiped, thus decreasing the mismatch between computed and observed LAB780-to-GAL16 amplitude ratio.
Modeling does not aim at quantitative comparison between modeled and observed strain, because very simple conditions (uniform force distribution, rectangular schist walls, homogeneous rock properties, etc.) are assumed and the summer relative contribution of Lecherines karst and compressional forces on schist depends (at least) on rainfall distribution. Nevertheless, sign and order-of-magnitude of the observed strain signals are correctly obtained.

Poroelastic Strain Evolution for an Impulsive Snowfall
When sudden uniform surface load is applied on the uppermost part of the computational domain (Section 6.1), computed strain evolution follows a quasi-exponential law at all the five test points. However, initial undrained strain and final drained strain are very different from test point to test point (Figures 5b  and 5c), although test points are close to each other; initial strain is always positive, but final strain may be either positive or negative.
From dimensional analysis and numerical tests, not shown here for brevity, the characteristic time of strain evolution scales as H/d 2 , where H is the aquifer thickness and d is hydraulic diffusivity. Characteristic time is as short as a fraction of a day (as in Figure 3a) when hydraulic diffusivity is of the order of tens of m 2 s −1 (20 m 2 s −1 in Figures 5b and 5c). Such diffusivity is large for competent rock; however, it is comparable with the effective diffusivity obtained elsewhere for fractured rock (e.g., as for the Gran Sasso, Italy, fractured aquifer, Adinolfi Falcone et al., 2008;Amoruso et al., 2014). Hydraulic conductivity, which corresponds to model poroelastic parameters, is 2 × 10 −5 m s −1 , at the lower limit of a "good aquifer" after Bear (1988).

From Impulsive to Real Snowfall
Applying the linear input-output polynomial model (Section 6.1) to a selected event, we obtain the results shown in Figure 6; characteristic time of the undrained-to-drained phase transition is around 0.5 days, lag is around 2 h although the elastic undrained response should be instantaneous. Other clear snowfall-associated strain signals give consistent results. Thus, poroelastic response to snow load satisfies snowfall-associated strain evolution, but for an initial lag. Figure 7 shows computed undrained/drained deformation along the directions of the two strainmeters at the nominal strainmeter altitude (1,180 m asl) in an area 4 × 4 km 2 , when uniform surface load is applied on areas between 1,200 m and 1,800 m asl (Section 6.1). Consistently with observations, at the strainmeter site undrained (initial) deformation is extensive along both strainmeters and about equally large; drained (final) deformation is extensive (but very small) along LAB780 and compressive along GAL16. To check model sensitivity to the loaded altitude interval, we change its endpoints by up to 200 m independently. Although computed values (not shown here for brevity) vary from test to test, undrained (initial) deformation is constantly extensive along both strainmeters and about equally large. Drained (final) deformation may be extensive or compressive but is constantly small along LAB780, whereas it is more compressive (with respect to LAB780) along GAL16.

Endpoints of Poroelastic Evolution for Real Topography
For completeness, we also compute deformation caused by snow load acting above 1,800 m asl only (Figure S12). At the interferometer site, computed undrained and drained elastic deformation is small along LAB780, but large (more than three times as large as deformation caused by snow load acting between 1,200 m and 1,800 m asl) and compressive along GAL16. Figure 6. Linear input-output polynomial model of a snowfall-induced strain signal. Panel (a) shows observed and modeled strain signals (black and red lines respectively), precipitation, water table level, and Aragón stage for a selected snowfall episode occurred during a period of persistent low temperature. Strain has been arbitrarily shifted for clarity. Panel (b) shows modeled strain response (red lines) to an impulsive precipitation (blue bars). A simple linear inputoutput polynomial model with three parameters has been used for modeling (see text and Text S1 for details). Figures S9 shows how low snow line altitude is.

Discussion and Conclusions
A preliminary analysis of four hydrology-associated deformation signals recorded by LAB780 and GAL16 is in Díaz et al. (2014), which focused on seismic monitoring of the Aragón river. Here, we exploit the whole set of available strain records and a more comprehensive set of environmental data, not at our disposal when Díaz et al. (2014) was written; thus, interpretation of strain data is more accurate. The strain signal, shown in Figure 5 of Díaz et al. (2014) and labeled as severe storm, was related to a regional-scale exceptional rainfall and flood event occurred on October 19-20, 2012. Unfortunately, during this event only LAB780 data are available, measured strain being extensive as after winter snowfalls; thus, strain might have been caused by rainwater load, analogously to winter snow load. As for the other strain signals shown in Díaz  2014), we confidently ascribe the moderate rainfall episode in Figure 3 of Díaz et al. (2014) (November 26, 2012) to the late-spring/autumn rainfall class, and the signals shown in Figure 6 of Díaz et al. (2014) (December 14, 2012and January 19, 2013 to winter snowfalls. In this work, we interpret hydrology-associated strain signals basing on simple models. We assume uniform force distributions (e.g., snow surface load, push on the schist unit, overpressure inside the effective N21°E fault), homogeneous or layered elastic/poroelastic medium, simple plane geometries (e.g., for the schist unit and the effective N21°E fault). Nature is much more complex, for example, the aquifer is compartmentalized and rock is heterogeneous and karstic, with severe effects on both mechanical and hydrological underground properties. Environmental time series (e.g., water table level, precipitation, snow thickness) are, at best, indicative or, at worst, suggestive of local behaviors. Thus, realistic quantitative modeling is actually unachievable because available data are scattered and insufficient to characterize the study area; nevertheless, our simple modeling seems able to capture basics of the processes underpinning the hydrology-associated strain signals we observe.
The most prominent signals in the raw strain records and in the scalograms are linked to heavy snowfalls occurring during very cold days. Recorded deformation (Figures 3d and 6a) may be explained through the rock poroelastic response to snow load acting over areas between about 1,200 m and 1,800 m asl (Section 5.3). Surface load induces uneven rock deformation, which in turn causes groundwater flow through the Ezcuarre-Peña Telera aquifer. This flow manifests as changes of the water table level at station PZ32 and the Aragón river stage at station A271 (Figures 3d and 6a).
Actually, data show that snow load acting over about 1,800 m asl is ineffectual in producing either deformation at the two strainmeters or hydrological effects (water table level at PZ32, Aragón river stage at A271; Figure S11). This lack of effects contrasts with what we obtain when we compute deformation caused by snow load acting above 1,800 m asl (Section 6.2.2.3). The discrepancy between model results (small effects on LAB780, large GAL16 compression) and observations (negligible effects on both LAB780 and GAL16) might be due, at least in part, to the presence of vertical faults which can accommodate horizontal deformation, as well as to the complex relations linking surface load, internal stress, and deformation in the vadose and epiphreatic karstic zones which exist at high altitudes (Section 2). Similar unmodeled relations might underpin the lag (around 2 h) between deformation and snowfalls ( Figure 6b).
Late-spring and autumn strain signals are smaller than those caused by heavy snowfalls. We model deformation effects invoking water infiltration close to the strainmeters and consequent N-S compressional forces on the schist unit from the karst aquifer (Section 5.1). Because of the different directions of the two strainmeters, those forces cause extension to GAL16 and compression to LAB780 (Figures 3a, 4b, and 4c). We invoke the same mechanism to explain diurnal strain cycles (Figure 3b) which occur when snow melts at altitudes between approximately 1,500 m-2,000 m asl. Diurnal strain cycles are not clearly visible in early summer, when snow melts above about 2,000 m asl. In the west, such areas are close to the strainmeters, but belong to the sinkhole hydrographic basin and feed the Tortiellas lake; in the east, they are a couple of kms from the strainmeters (Figures 1c, 1d and 4a). Thus, in both cases, snowmelt is not expected to feed limestone close to the strainmeters.
During snowmelt-induced diurnal cycles, LAB780 compression phases are coeval with GAL16 extension ones and typically occur in the afternoon (Figure 3b). Such timing and periodicity are consistent with hydrograph measurements (e.g., Singh et al., 2000) and with dynamics of snowmelt in Alpine regions.
Summer and early autumn rainfalls often affect large areas; related deformation is quite evident in raw data but often unclear in scalograms (Figures 2 and 3c). Similarly to late-spring/autumn rainfalls, we ascribe deformation to water infiltration close to the strainmeters and consequent compressional forces on the schist unit from the karst aquifer (Section 5.1). In summer, pushing probably extends further west with respect to the periods when the western zone is snowy, thus compressing both strainmeters. In addition, faults and karst west of the interferometers (Tortiellas sinkhole and Lecherines area) probably contribute (Figures 4f and 4g), because temperature inside the karst conduits is higher than 0°C (Bartolomé, 2019) allowing water flux.
Our modeling implies the same source of both LAB780 and GAL16 deformation as for spring/autumn rainfalls and diurnal cycles during snow melting, that is, groundwater push on the northern and southern faces of the schist unit; as for summer rainfalls, GAL16 deformation is also due to water pressure inside the faults and karst channels of the Lecherines area. This conclusion is supported by data; LAB780 and GAL16 strain time histories are similar (although opposite in sign) when late-spring/autumn rainfalls occur and during snow melting, but may be quite different (although equal in sign) when summer rainfalls occur (Figures 3a-3c, S2, S5, and S7). Moreover, when summer rainfalls occur, LAB780-to-GAL16 amplitude ratio is particularly variable (Figure 3c), probably because of the inconstant relative contributions to GAL16 deformation from the two deformation sources (push on the schist unit and overpressure in the Lecherines karst).
Each observed hydrology-associated deformation signal can be confidently inserted into one of the four classes-the class-by-class different response of GAL16 and LAB780 being its distinctive signature-and classification is fully consistent with what we expect, based on when it occurred. Our modeling is able to justify the basic characteristics of each class, including the strain order of magnitude; however, as for signal details, we could not find any robust correlation between observed signal amplitude and precipitation size. Also, details of the strain time history may vary. Differences may be due to infiltration and/or surface load inhomogeneities, as well as to the complex underground hydrogeological set up, which can lead to nonlinear, intrinsically time-dependent, responses. Variety of signal details makes it very difficult to foresee system evolution and fitting strain data accurately, and the above-mentioned limitations are not peculiar to our study site.
This work refers to a particular mountain zone, the Central Pyrenees, which is considered one of the most active seismic areas in western Europe. Nevertheless, our insights on hydrology-associated deformation are potentially useful also in other mountain seismic areas (e.g., Italy, Balkans, Turkey, and China). Generally speaking, the variety of hydrology-associated deformation we report warns all kinds of deformation data, as for causative physical mechanisms, their difficult modeling, and the extreme care that is needed when removing hydrology-associated effects, e.g., before attributing apparent deformation anomalies to geological disaster precursors; unfortunately, careful analysis is seldom carried out in the literature.