Recurrence of Fault Valve Behavior in a Continental Collision Area: Evidence From Tilt/Strain Measurements in Northern Adria

We analyzed the data recorded by the NE-Italy subsurface tilt and strainmeter network evidencing a coherent transient signal in the recordings of four tiltmeter sites in the 1984–1990 period that produced a tilt along the main fractures. Borrowing from classical seismology techniques, we used the uprise times to locate the transient signal source. The propagation velocity is compatible with a fluid diffusion process that starts from a source located close to the hypocenter of the February 10, 1983 Uccea earthquake, MD = 4.2 at the Italy-Slovenia border, at an estimated depth of 10.8 km. Our results add to the previous interpretation of a transient signal recorded by several global navigation satellite system (GNSS) stations in the 2006–2009 period in terms of fluid diffusion below the Bovec basin (Slovenia). That source was located upon continuation to the northwest of the Ravne fault, few kilometers to the northeast from the present one, and about 6 km from the July 12, 2004 Bovec–Krn earthquake, Mw = 5.1, depth ~6.1 km. These observations suggest that the area is subject to fault valve behavior episodes that released fluids trapped at depth to the surrounding region as pore-pressure bulges. The convergence between Alpine and External Dinarides structures in this area puts highly permeable dolomitic limestones in contact with low-permeable fine-grained limestones and flysch formations. Therefore, the conditions for overpressure generation can be created, whereas fault movements, from time to time, in close relation with seismic events, can enable fluid diffusion in the surroundings. We also estimated the possible fluid influx needed to maintain overpressure and possible discharge across both the faults. The study provides insights on pore–fluid pressure variations related to slow slip events from a context different from subduction or transform margins, i.e., in a continental collision area.


INTRODUCTION
Nowadays, space geodetic observations cover large portions of the Earth, with time series even longer than 20 years, i.e., long enough to allow distinguishing and characterization of periodical or transient variations. The latter consist of variable amplitude signals (from few millimeters to some centimeters, Riel et al., 2014) due to slow slip events often associated with tremors [episodic tremor and slip (ETS)] induced by tectonic or volcanic processes. They have been observed in various parts of the world (e.g., Dragert, 2001;Lowry et al., 2001;a review in Bürgmann, 2018, and references therein).
Independently upon the fluid-origin, pore-pressure variations are often thought to be at the origin of slow slip events, ETS, and seismicity (Rice, 1992;Kodaira, 2004;Skarbek and Rempel, 2016;Cruz-Atienza et al., 2018;Nakajima and Uchida, 2018), but most of the voluminous literature regard oceanic-continental plate subduction zones and, in lesser part, volcanic, hotspot, or transform settings (Khoshmanesh and Shirzaei, 2018). Recently, however, fluid migration related to seismic activity has been reported in the Mediterranean, in different parts of the Adria microplate margins (Rossi et al., , 2017(Rossi et al., , 2018Nespoli et al., 2018;Console et al., 2020).
In particular, in the northern tip of the Adria microplate, an area of continental collision , corrected in 20172018-hereafter: RossiII), analyzing the global navigation satellite system (GNSS) data from the Friuli Regional Deformation Network (FreDNet) (OGS), Marussi (Friuli Venezia Giulia regional council), and the EUREF networks, reported the occurrence of a transient signal propagating through the region that, between 2006 and 2009, induced a tilt parallel to the principal faulting lineaments. They interpreted the signal as a propagating surge in pore fluid pressure (Revil and Cathles, 2002;Lupi et al., 2011), i.e., a porosity wave. Porosity waves are packets of fluid-filled, hydraulically connected cracks that propagate in response to a pressure gradient set by ambient tectonic stresses (Connolly and Podladchikov, 2013). They have been hypothesized by Skarbek and Rempel (2016) as the transport mechanism of deep fluids triggering ETS in subduction zones. Also, porosity waves can originate in the presence of suprahydrostatic gradients in fluid pressure across an active fault (Lupi et al., 2011) that may behave as a fluid pressure-activated valve, the so-called fault valve behavior of Sibson (Sibson, 1992(Sibson, , 2014Fischer et al., 2017;Matsumoto et al., 2021). According to this model, during the interseismic and pre-seismic phases, the fault acts as an impermeable barrier to fluid flow from a region characterized by a supra hydrostatic regime to the overlying layers where the fluid pressure of interconnected pores is hydrostatic ( Figure 1A). When failure condition is met, through further increase in fluid pressure or accumulation of shear stress, an earthquake occurs, and the fault and the barrier break, allowing fluid discharge along the fault ( Figure 1B). The scheme in Figure 1C shows the time behavior of a fault through time: in the pre-seismic phase (to the left of the earthquakes, indicated as EQ), we observe a slow increase in fluid pressure (P f ) and shear stress (τ ). In correspondence with the earthquake, the fault permeability (k) suddenly increases, and the following fluid discharge progressively reduces P f and τ , with opposite behavior of the frictional shear strength of the fault (τ f ). Progressive healing due to hydrothermal processes decreases permeability, restoring the impermeable barrier, and the cycle can start again.
Although much lower than the fractured fault, the permeability across the fault is not null. Fluids can propagate at a lower rate outside the fault zone in the porous and fractured rocks that constitute a deformable matrix. Hence, a fluid flow would lead to a porosity transient that propagates in a waveform that could be detected by GNSS sensors and tiltmeters, as Braitenberg et al. (2019) demonstrated for pressure changes in karst conduits (Figure 2).
In order to support that the transient signal in the different permanent GNSS stations was due to a porosity wave, and, hence, fluid diffusion, RossiI and RossiII inverted the arrival times of the transient, obtaining both propagation velocity and hydraulic diffusivity values consistent with the fluid diffusion hypothesis, and, from these, they inferred permeability. The permeability values turned out to be consistent with the rocks present in the area: highly permeable dolomitic limestones, low permeable finegrained limestones, and flysch formations. Such rock formations are involved in the overlapped thrusts and imbricate structures created by several tectonic phases that affected the region of study. The contact between permeable and poorly permeable rocks in a compressive and transpressive tectonic regime can give rise to overpressure conditions favorable to the fault valve behavior (Sibson, 1992(Sibson, , 2014. RossiI and RossiII, based on the location they inferred for the porosity-wave source, close to Bovec (Slovenia), hypothesized a valve behavior of the Ravne fault, a transpressive fault of the External Dinarides, recognized as being responsible for the Bovec-Krn seismic events (1998 M w = 5.6 at 6.9 km depth; 2004 M w = 5.1 at 6.1 km depth) (Kastelic et al., 2008). In particular, the transient appeared related to the 2004 earthquake. Unfortunately, there are no GNSS data that could confirm such behavior for the 1998 earthquake.
The hypothesis that the phenomenon reported by RossiI and RossiII is not an isolated episode but could have occurred in the past in the region is at the basis of the present research, FIGURE 1 | Schematic representation of the fault valve behavior. (A) Pre-failure: on the left, the pore pressure gradient in the fault zone, above and below the impermeable barrier, i.e., the transition from the hydrostatic to the supra-hydrostatic flux regime (right); (B) post-failure: (right) an earthquake (red star) breaks the barrier, and the fluids rise upwards; on the left, the fluid pressure gradient is nearer to the hydrostatic curve; (C) scheme of fault valve action in relation to earthquake occurrence (EQ): when the impermeable barrier is broken, the local permeability k increases, allowing fluid discharge, decreasing fluid pressure P f , and increasing the frictional strength (τ f ) (from Sibson, 2020, redrawn). P l , lithostatic pressure; P h , hydrostatic pressure; τ , tectonic shear stress (From Sibson, 1992, redrawn).
FIGURE 2 | Conceptual model of a porosity wave originating at a fault (red line) and propagating in the hosting rocks. Geodetic instruments, namely, tiltmeters-horizontal pendulums (installed underground, when possible, in natural caves) and GNSS sensors (on the surface) record the induced deformation at different times, depending upon the path of rays, normal to the propagating front (blue dashed lines).
focused on re-analysis of old tiltmeter data in light of more recent GNSS observations. The region, in fact, has also been continuously monitored for over 50 years by the NE-Italy subsurface tilt and strainmeter network (Braitenberg, 1999). The Grotta Gigante horizontal pendulum station has been active since 1964 (Marussi, 1959;Braitenberg et al., 2006), wherein tiltmeters and strainmeters have been installed starting from 1977. In the 3 years before the 1976 Friuli M w > 6 earthquakes, the Grotta Gigante pendulums recorded "tremors" and an anomalous southward tilt (Chiaruttini and Zadro, 1976;Dragoni et al., 1985). Bonafede et al. (1983) interpreted the tremors in terms of silent earthquakes that originated from the deep non-brittle part of the fault. The frequency is proper for the very low frequency earthquakes (VLFEs) reported by Ito et al. (2007). The suggestion that a transient similar to the 2006-2009 episode could have occurred in the same region in the past came from the observations by Rossi and Zadro (1996) and Zadro and Braitenberg (1999) on long-term variations affecting the signals of the regional tiltstrainmeter network and elastic parameter changes between 1983 and 1986, evidenced by Mao et al. (1989).
Hence, in particular, we focused the analysis on a coherent transient tilt signal, recognizable in the interval 1984-1990, by applying a data-driven approach to the time series. We used the same tomographic approach of RossiI to locate the source of the transient while calculating propagation velocity and hydraulic diffusivity. Finally, we verified whether the associated permeability values are compatible with the rock formation in the region and whether there are conditions for overpressure. To verify the plausibility of the fluid diffusion hypothesis, we also estimated the

THE TECTONIC FRAME AND THE DATA
Several tectonic phases affected the region, originating different tectonic systems interfering with each other. The primary tectonic systems present in the region are those of the Alpine (with S-directed E-W-striking thrusts) and the Dinaric, represented by the External Dinarides (with NW-SE-oriented dextral transpressive faults) (Slejko et al., 1999;Poljak et al., 2000;Bressan et al., 2007). Originated by the collision between the Adria and Eurasia plates and the counter-clockwise rotation of the former, the two systems merge and intersect in the central part of the region. Here, we record the highest seismic activity (Figure 3), characterized by an intricate pattern of thrusts and strike-slip fault mechanisms, with a compressional axis oriented from NW-SE to NNW-SSE and N-S (Bressan et al., 2018b).
The strain field has been continuously monitored for over 50 years by the NE-Italy subsurface tilt and strainmeter network (Braitenberg, 1999). Presently, it consists of two long baseline horizontal pendulums installed in Grotta Gigante (GG), near Trieste (installed in 1959;continuous data:1964-present), and two tiltmeter stations located in the Bus de la Genziana (BG) and Villanova (VI) caves, active since 2005 and 1977, respectively (Figure 3, Table 1; Chiaruttini and Zadro, 1976;Braitenberg, 1999;Zadro and Braitenberg, 1999;Grillo et al., 2018). The VI and BG tiltmeter stations are equipped with two Marussi-type, horizontal pendulums with Zöllner suspension oriented to NS and EW, respectively, with a cast iron conic housing. The nominal angular resolution is 2.5 nrad, and signals are recorded with a sampling rate of 1 h. GG instruments are Zöllnertype horizontal pendulums of exceptional dimensions with an oscillation period of 360 s and an amplification factor of equal to 24,000. VI also hosts three Cambridge wire strainmeters (King and Bilham, 1976) in the N128E, N27E, and N67E directions. The network included at its maximum expansion (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) four further tiltmeter stations in the Piedmont area of Friuli (CE, BA, IN, and GE) with tiltmeters of the same kind of those in VI and BG (in italic in Table 1; Figure 3).
This study considers the continuous recordings of GG, VI, CE, and IN stations that worked continuously between 1979 and 1996 ( Figure 4). Although the raw data have different amplitudes, depending on the instrument and site characteristics (Rossi and Zadro, 1996), some features are recognizable for all stations in the two components (EW and NS). Pluri-annual fluctuations modulate seasonal (annual and semiannual) oscillations. In some cases, a clear linear long-term trend is recognizable (e.g., CE and VI NS-component) (Figure 4). Higher frequency signals (some days) (e.g., in the EW component of VI and GG between the end of 1991 and the beginning of 1992) are correlated to hydrological effects (Braitenberg et al., 2019).
The application of a running average over 2 years evidences long-term variations (Figure 4). In particular, one fluctuation with a transient character and that affected all the recordings between 1984 and 1990 is evident (blue ellipse in Figure 4). For GG, we note another higher amplitude fluctuation in the 1989-1994 period (dashed blue ellipse), which will also be analyzed in the following chapter. The search for  (Figure 5; RossiI). The 15 • spacing enables a first reasonably accurate scan of the direction with the maximum amplitude to be refined later.

DATA ANALYSIS Time Series Analysis
As shown in Figures 4, 5, the oscillations have transients and likely dispersive characters, so that uprise and duration could change in the different sites. Hence, intending to exploit the time series in its whole length and avoiding the risk of the filtering parameters conditioning the results, in this study we applied a different method from the filtering of RossiI and RossiII: the data-driven empirical mode decomposition (EMD) technique, initially proposed by Huang et al. (1998). Since its introduction, EMD has been successfully applied in various fields (Ditommaso et al., 2012;Zhang et al., 2019). EMD decomposes a noisy, complex signal into a variable number (N) of intrinsic mode functions (IMFs) through an iterative sifting process. Based on the detection of local maxima and minima in the time series, sifting produces a (generally limited) number of IMFs, each containing higher frequency oscillations than the following IMF component. Orthogonality is not a necessary condition of the EMD technique that provides a complete and adaptive basis for expanding linear, non-linear, and non-stationary data.
The original signal x(t) is expressed as where r n (t) is the residual. We performed EMD analysis on the projections calculated every 15 • (Figure 5). The first 5 IMFs capture the highest frequency signals, including those induced by rainfall, groundwater variations, and co-seismic displacements (red ellipses in original data, IMF1 and IMF2), whereas IMFs 6 and 7 represent the seasonal, periodic, and semi-periodic components (Figure 6). IMF 8 contains a longer wavelength term, dominated by a transient surge, recognizable for all the stations (blue ellipse). For each site, once we found the direction in which IMF 8 shows the maximum amplitude, we repeat the analysis by moving from directions 5 and 10 • to North and South, to be sure to find the direction of maximum amplitude of our signal, with an uncertainty of only 5 • . The maximum absolute amplitude is recorded in directions N85 • E for VI, N135 • E for CE, N10 • E for IN, and N135 • E for GG (blue arrows in Figure 7). The positive tilt is N85 • E for VI, N45 • W for CE and GG, and N170 • W for IN (darker blue). Due to the time delay between the different sites, the image refers to the end of 1985: GG shows a null tilt (light-blue arrow). Transient oscillation is recognizable in the VI signal in 1983-1987, in CE and IN signals in 1984-1988, with slight differences between the two sites, and in the GG signal in 1988-1992. There is a substantial agreement with both the tectonic lineaments existing in the neighbors of the four sites and the directions in which the lower crust background seismicity tends to align, as evidenced by principal component  The four tiltmeter time series transient uprise times progressively increase, starting with VI, then with CE and IN northwesterly, and with GG southeasterly (Figure 8). Notably, the transient uprise time in GG is consistent with the 1989-1994 oscillation observed in the original time series (dashed ellipse in Figure 4).

Transient Source Parameters
To locate the source of the unknown transient (in space and in time), we applied a standard, simplified earthquake location procedure based on circles, whose radii are estimated by the arrival times of the transient and its velocity of propagation (RossiI). We estimated the average propagation velocity based on various distance/differential uprise time ratios between couples of tilt recordings. Then, we refined localization through a tomographic approach, i.e., ray tracing in a stratified medium characterized by laterally varying velocities (Böhm et al., 1999), minimizing the observed and calculated travel time differences. The technique is depicted in Figure 9. For both the frequencies involved and the limited number of observations, we are aware that we are working at the limits of the method. However, the staggered grid procedure limits the drawbacks of inadequate coverage (Vesnaver and Böhm, 2000). Therefore, we initially chose a very coarse grid, consistent with the broad wavelength of the pulses, to have a rough but reliable velocity distribution.
We started from the RossiI model ( Figure 10A) and performed ray tracing, following the principle of Fermat. The 3D model origin is at 11.8 • E, 45.6 • N, and it extends to 14.6 • E, 46.7 • N, from the Moho surface to the topographic surface, with three surfaces in between at average depths of 2, 4, and 11 km, respectively (Rossi I). We refined the velocity model iteratively, together with the source location. When the difference between the observed and calculated travel times was minimal, we applied to the grid small shifts in x and y directions, and repeated the inversion. At the end of the process, we averaged the results, obtaining a higher resolved image without losing the reliability of the original grid (Figure 9; Vesnaver and Böhm, 2000).
The velocities of the final model range between 6.4 and 27 km/year (Figure 10). Notwithstanding the more limited area covered by the rays, the values are close to the ones obtained by RossiI. In the second and third layers, the propagation velocity does not exceed 16 km/year (±0.5 km/year; Figures 10B,C) in the volumes crossed by the rays. The shallowest layer (layer 1) shows the highest velocity values, up to 27 km/year (±0.9 km/year) ( Figure 10D). The source of the transient (green circle     in Figure 7) is located at 46 • 16'50.50 "N, 13 • 19'0.40" E, hence about 5 km to the East of the VI station, at 6 km depth, with an uncertainty of about 7 km in the plane, and 1 km in depth. Low velocities (∼8 ± 0.2 km/year) characterize the "hypocentral" area ( Figure 10B). Origin time is the end of November 1982 (red star in Figure 8), with an uncertainty of 1.5 months.

Fluid Diffusion
We verified whether the uprise times are compatible with an episode of fluid diffusion in the present case, as in the case reported by RossiI. For this aim, we used the uprise times of the transient as input for hydraulic tomographic inversion to obtain hydraulic diffusivity (Brauchler et al., 2013). We followed the scheme shown in Figure 9, excluding the source location, which is now fixed, and started from the model of hydraulic diffusivity obtained by RossiI. Figure 11 shows the hydraulic diffusivity (D h ) values obtained for the same layers, shown in Figure 10.
The values are in the range 0.1 · 10 −3 -3.4 · 10 −2 m 2 /s (±.02 · 10 −3 m 2 /s). The highest values are in the first layer (Figure 11C), and the lowest is in the second one (Figure 11B), as observed by RossiI. Since permeability is a physical property for which more measurements are available than hydraulic diffusivity, it is more suitable to compare the obtained values with those that can be expected from the lithology present in the region. The equation binding the two quantities is (Talwani et al., 1999).
where β f , β r , η f , and φ are the fluid and rock compressibility, fluid dynamic viscosity, and rock porosity, respectively, the latter available from a borehole and laboratory measurements (Faccenda et al., 2007). The rays mainly cross the Jurassic limestones and the Cenozoic flysch formations. The finely layered limestones of Jurassic are characterized by permeability values in the order of 10 −17 m 2 , whereas the Cenozoic flysch shows higher permeability (∼4 · 10 −16 m 2 ) ( Table 2). The effective stress σ 0 associated with the pulse is (RossiI) where ρ f and ρ g are the density of fluid and rock, g is the gravity acceleration constant, v is the propagation velocity, and D h is the hydraulic diffusivity. We chose for the Jurassic limestones and the Cenozoic flysch the propagation velocity and hydraulic diffusivity values in locations where these values can be assumed to be dominant at the depth crossed by the rays, based on the available stratigraphic information (Faccenda et al., 2007;ViDEPI, 2009).
In the absence of other information, we assumed water density for fluid density. The density of the grains in the rock ρ g values, known from laboratory measurements, are obtained from Table  2 of RossiII. The effective stress obtained using Equation (2) is quantified as about 19 MPa, with a standard deviation of 0.11 MPa. To see whether this value implies an overpressure condition, we  calculated the vertical stress due to lithostatic loading σ v : where P p is the pore pressure, λ is the pore fluid ratio (Terzaghi, 1923;Hubbert and Rubey, 1959), and ρ b is the bulk density of sediment, defined as We considered the same vertical profiles reported by RossiII, since they are also very close to the source of the transient in this study (A,B,C,D,E,F and S2,S4 in Figure 3). In the present case, given that the origin of the transient is shallower, at a depth of about 6 ± 1.0 km, the lithostatic load is about 167 ± 25 MPa with a standard deviation of ∼2.5 MPa. The pore fluid ratio λ is about 0.89 ± 0.01, indicating condition of a supra-hydrostatic state, with pore pressure very close to the lithostatic load.

DISCUSSION
The study of transient variations in deformation requires several years of observations of long-term inter-seismic strain variations, necessary to reveal a particular region or a repetitive behavior of a structure. Since GNSS and interferometric synthetic aperture radar (InSAR) data are limited to the last 20 years, observations on past episodes of slow slip can come from strainmeters, tiltmeters, and dilatometers that had been continuously recorded in the past. This study focused on long-term oscillation of the strain field in NE-Italy observed in the second half of the last century, reported by Rossi and Zadro (1996) and Zadro and Braitenberg (1999), not correlated with meteorological factors, and that has similarities with the recent GNSS observations in the same region (RossiI and RossiII). We chose the data-driven EMD technique to delineate, at best, transient non-stationary signals in the recordings. This choice has the advantage of avoiding any a priori hypothesis and is independent of any parameter choice while exploiting the time series in its whole length, which is essential when considering transients with some years of duration. The IMF 8 resulting from the tiltmeter data analyses revealed a transient pulse, slowly propagating throughout the region between 1983 and 1987 (Figure 8): the uprise times in the various stations slightly increased from VI to IN-to the west of VI-and to GGto the east of VI. There is a substantial agreement between the directions of the tilt induced by the transient (blue arrows in Figure 7), the directions found through the analysis of the GNSS data (red-black arrows; RossiII), and the directions in which the background seismicity tends to align (gray dashed arrows; Bressan et al., 2016Bressan et al., , 2019. We can argue that the propagation follows the primary fractures, presenting the same orientation (Serpelloni et al., 2016;Peterie et al., 2018).
The velocities range between 6.4 and 27 km/year (Figure 8), in good agreement with the results of RossiI. Such velocities strongly suggest pore pressure diffusion, and a hydraulic tomographic inversion of the surge times in various tiltmeter stations provided the hydraulic diffusivity values: the values range between 0.1 · 10 −3 and 3.4 · 10 −2 m 2 /s. The excellent agreement of the derived values for permeability of the Jurassic limestones and Cenozoic sandstones with the values reported in the literature (Hawle et al., 1967;Sibson and Rowland, 2003) strongly supports the hypothesis that the transient signal is the expression of pore fluid pressure diffusion throughout the region.
The origin time (end of November 1982) and the location of the transient source (Figure 7, green circle) indicate a close relationship between fluid surge and the 1983 February 10 M D = 4.2 Uccea earthquake, 6.5 km apart (Figure 7, red star). We can estimate the uncertainty in origin time to about 1.5 months, about 7 km for the location in the plane, and about 1 km in depth. The 2006-2009 origin of the transient (RossiI) (orange pentagon in Figure 7) was located in the northwestward continuation of the Ravne fault, responsible for the 2004 earthquake (M w = 5.1) (orange star in Figure 7), about 6 km from the epicenter of that earthquake (orange star in Figure 7). The earthquake occurred 3.5 months after the calculated origin time of the transient. In that case, uncertainty in origin time was estimated at about 15 days, whereas the location uncertainty was about 2 km in the horizontal dimension and 0.8 km in depth.
The two sources are <20 km apart, both in the area where the Dinaric fault system crosses and merges with the Alpine one, near the border between NE-Italy and W-Slovenia. In both cases, the effective stress that can be calculated as an origin for the pore fluid pressure surge results low (19 MPa for this case study, 22 MPa for the 2006-2009 one), implying pore pressures about 0.9 the lithostatic one.
These observations suggest that the area mentioned above is subject to fault valve behavior episodes (Sibson, 1992), releasing fluids to the surrounding region as pore pressure bulges. The tectonic style of the area, a complex superposition of thrusting (Alpine) and transpressive structures (the Dinaric ones), is most favorable to fault valve behaviors (Sibson, 2020). The model of Sibson implies accumulation of a volume of high-pressure fluid beneath a low-permeability seal at the base of a fault zone (interseismic period); an earthquake would break the seal, and the fluid would surge upward immediately following the earthquake (Figure 1).
Although the two transients relate to two different faults, some characteristics of the 2004 Bovec-Krn M w = 5.1 earthquake sequence support the interpretation that, at least in that case, the Ravne fault acted as a valve. Gentili (2010) recognized seismic quiescence preceding the 2004 M w = 5.1 Bovec-Krn earthquake, which is compatible with an overpressure condition. On the other hand, the depth distribution of the aftershocks that followed the main shock is compatible with upward spreading, as it would be in the case of a fluid discharge (Bressan et al., 2018a). The seismic monitoring network was still limited at the time of the 1983 Uccea earthquake; however, in the preceding year, we observe a reduction in seismic activity close to the forthcoming epicenter. The aftershocks were few and had tendency to become shallower.
Several models have been proposed to relate overpressurized fluids with faulting (Hickman et al., 1995;Faulkner et al., 2010). The model of Rice (1992), formulated to explain the weakening of a fault, implies that a fault zone behaves as an impermeable barrier to the country rocks. A continuous fluid flux from the lower crust maintains high pore pressure. On the contrary, the model of Byerlee (1990) assumes that water from the country rocks saturates a highly permeable fault zone first, and, when the shearing causes fault zone compaction, is released back to the country rocks, until the flow is arrested because of hydrothermal deposition. It results in a complex system of seal-bound compartments and heterogeneous distribution of overpressurized fluids within a fault zone (Khoshmanesh and Shirzaei, 2018).
Some independent observations can help understand which model can be more suitable in describing the phenomenon of study.
As it is known, the state of overpressure can be related to seismic velocities ratio (V p /V s ) anomalies (Dvorkin et al., 1999;Chiarabba et al., 2018). The source of the transient in this study coincides with a higher value (1.91) (Bressan et al., 2012). The same authors report a value of 1.86 in correspondence with the source of the 2006-2009 transient. Although high V p /V s values are proper for limestones and dolomitic limestones, such values can be interpreted as indicative of a state of overpressure. Unfortunately, these values result from the seismic activity analysis over a relatively long period of time, encompassing different phases of the seismic cycle. A different analysis is reported by Mao et al. (1989), who analyzed the seismic velocities in the same region using a sliding time window. They also compared the amplitude of Earth strain tides observed by the three strainmeters located in the VI station with theoretical ones, in 1978-1986, and recovered the elastic parameters. The results from Mao et al. (1989) can be seen in Figure 12, together with the V p /V s that we calculated, starting from the velocity values reported in Mao et al. (1989)'s Table 3. It is noteworthy, that they observed changes in velocities and elastic parameters in the same years with a wavelength similar to our transient. In particular, Mao et al. (1989) focused on changes in seismic velocities (3%) and bulk modulus K (50%), started in 1982 and recovered after the 1983 Uccea earthquake, about a year later. The increase in observed bulk modulus (K O ), compared with the theoretical one (K T ) in correspondence with the earthquake could agree with a state of overpressure at the time of the event. Notably, the V p /V s ratio markedly increased from the end of 1982 to the end of 1986, in agreement with the surge in the tilt anomaly in VI. However, the information is relative for a broader part of the region and not limited to a single fault, as, e.g., Chiarabba et al. (2018), and, hence, may reflect pore fluid diffusion through the rock formations hosting the fault. The increase in ratio between  Table 3 of Mao et al. (1989), with the error bars; (C) the ratio between the observed rigidity (µ O ) and bulk modulus (K O ) and the theoretical rigidity (µ T ) and bulk modulus (K T ), from Mao et al. (1989) redrawn.
the observed shear modulus (µ O ) and the theoretical one (µ T ) that we observe from 1980 to 1986 suggests hardening of the rock matrix with time, interrupted by temporary decrease in the course of 1982, in agreement with the pattern shown in Figure 1. We have, however, to remind that the calculation of the shear modulus is strongly affected by the loading tides, and it is unstable, owing to the closeness in values of the two normal strain components of the observed tidal strains (Mao et al., 1989).
Our observations, together with the variations in elastic parameters reported by Mao et al. (1989), are in agreement with the time variation in the fluid content, possibly related to compaction-driven porosity and lateral and temporal variations in permeability (Byerlee, 1993;Khoshmanesh and Shirzaei, 2018).
The contribution of meteoric water cannot be excluded in the region of study, characterized by karstic phenomena in outcropping and subsurface formations (Cimolino et al., 2010), implying that complex groundwater flow in the rock matrix, fractures, and conduits also inferred from deformations (Devoti et al., 2015;Braitenberg et al., 2019), as in other regions characterized by karst aquifers (Silverii et al., 2016(Silverii et al., , 2019. In the present case, based on all the observations presented, the models of Rice (1992) and Byerlee (1993) appear the most suitable.
In the absence of geochemical information, the possibility of a fluid source from crustal deeper parts, as in the model of Rice (1992), could not be excluded. We can estimate the flux rate necessary to maintain the overpressure at depth and verify the plausibility of the values for both cases: the one treated in the present study, and the one reported by RossiI and RossiII. The fluid influx q is defined by David et al. (1994) as where K 0 is the permeability, and ξ is the normalized flux given by with σ * being the characteristic stress. We used the permeability values shown in Table 2 for the Jurassic limestones, whereas for permeability and density of the Dolomia Principale and Paleozoic sandstones, not crossed by the rays in this case, we referred to RossiII's Table 2.
In this case, we considered the formations that should be present at the depth of the two transient sources. Table 3 of David et al. (1994) gives the values of σ * : 21.3 MPa for Dolomia Principale and the Jurassic limestones and 26.3 MPa for the Paleozoic sandstones. For the two transients, we obtain the values of ξ in the range 0.21-0.33 for the permeable Dolomia Principale, and in the range 0.53-0.63 when Jurassic limestones or Paleozoic sandstones are considered, in agreement with the values reported by David et al. (1994).
It follows that for both transients the value of the fluid influx q is in the range 0.47-0.65 m/year for the permeable Dolomia Principale and about 0.013-0.015 m/year for the Jurassic limestones, whereas it is 0.023-0.027 m/year for the Paleozoic sandstones. If we use the values of the effective stress we calculated for both transients, instead of σ * , we have a variation of about 8% for ξ and q. The fluid influx values are reasonable compared with the measurements of the mantle fluid flow in different tectonic environments from geochemical measurements: the highest range of values, 0.42-2.6 m/year, is obtained for the Alpine Fault in New Zealand (Menzies et al., 2016), whereas 0.004-0.14 m/year characterizes the San Andreas Fault (Kennedy et al., 1997;Kulongoski et al., 2013), 0.13-0.19 m/year the North Anatolian Fault (de Leeuw et al., 2010), and −0.012-0.019 m/year the Karakoram fault (Klemperer et al., 2013).
To calculate the approximate areas involved in the flux discharge according to the fault valve model, we use the wellknown relationships between magnitude and rupture area (Wells and Coppersmith, 1994), considering the February 1983 M D = 4.2 and July 2004 M w = 5.1 earthquakes. These earthquakes occurred two to three months after the origin of the hypothesized transient. The highest flux discharge corresponds to the more permeable Dolomia Principale (7.2 · 10 5 -4.0 · 10 6 m 3 /year), whereas, for the much less permeable Jurassic limestones, nearlithostatic pressure would be maintained even with a flux range of 1.7 . 10 4 -1.1 . 10 5 m 3 /year. In all the cases, we obtain the highest values for the 2006-2009 transient and the greater magnitude of the correlated earthquake.

CONCLUSIONS
Transient variations in crustal deformation data can potentially detect slow pore fluid pressure propagation and associated effective stress fluctuations. The empirical mode decomposition (EMD) method employed here has been efficient in detecting transient tilt signals. Through the present analysis, we demonstrated that the area at the boundary between NE-Italy and W-Slovenia, where the Dinaric fault system crosses and merges with the Alpine one, strongly supports the recurrence of fault valve behavior, releasing fluids from depth to the surrounding region as pore pressure bulges. The pore fluid pressure appears to have migrated from a depth of 6 km upward and through the surrounding as a porosity wave, as inferred from a strong, coherent tilt transient of about 2.5 years in duration in the direction of the main fractures. At the origin of the fluid release, there are overpressure conditions at depth, with pore pressure ∼0.89 the lithostatic one, in agreement with previous and independent observations on elastic parameters and seismic velocity variations. The fluid influx values necessary to maintain such pore pressure values are comparable with other observations in other contexts. Geochemical tracer analysis of the spring water and concretions could help evaluate the budget of meteoric, metamorphic, and mantle fluids. Such information, together with full waveform inversion, more adequate in a near-field case (as for VI and the source of transient), and a pore fluid diffusion model, could better define this complex interaction between fluids, deformation, and seismicity.
Our research provides novel pieces of evidence on the relationship between pore fluid pressure variations and slow slip moderate magnitude events in an area of continental collision. Although based on an almost unique strain data set, our results and observations can help consider phenomena like the ones reported in this study as a possible cause of deformation and seismicity migration patterns in other regions.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available upon request from the authors. Requests to access these datasets should be directed to berg@units.it, apastorutti@units.it.

AUTHOR CONTRIBUTIONS
GR ideated the study, performed data analysis and interpretation of the results, and wrote the manuscript. IN and AP performed the preprocessing of data. CB is responsible for the maintenance of the NE-Italy tilt-strainmeter network of the University of Trieste, and discussed the results and interpretation. SP suggested the EMD technique, designed the EMD software, and participated in EMD data analysis. All the authors read and corrected the manuscript.

FUNDING
We acknowledge funding by the Italian Ministry of University and Research (MUR) Department of Excellence grant, given to the Department of Mathematics and Geosciences of the University of Trieste (Italy) through which AP is supported.

ACKNOWLEDGMENTS
We dedicate this study to the memory of Prof. Maria Zadro, who initiated the NE-Italy tilt-strainmeter network of the University of Trieste, and her scientific curiosity that inspired this research. We thank all the technicians and Dr. Barbara Grillo who, through the years, have helped in the management of the network. We thank Dr. Gianni Bressan for constructive discussions and the relocalization of the Uccea and Bovec-Krn earthquakes. We are grateful to Alexandre Canitano, Weijian Mao, Pierre-Michel Rouleau, and the editor, György Hetényi, for the accurate and thorough revision, stimulating comments, and helpful suggestions.