Hydrogeologic Property Estimation in Plate Boundary Observatory Boreholes Using Tidal Response Analysis

Because of the in ﬂ uence pore pressures have on e ﬀ ective stress, understanding hydrogeologic properties that control ﬂ uid ﬂ ow and pressure distribution is important in characterizing earthquake and deformation processes. Here, we utilize borehole pressure changes in response to earth tides to determine hydrogeologic properties and their time variations for 17 boreholes within the NSF Earthscope ’ s Plate Boundary Observatory (PBO) network along the San Andreas fault and Cascadia subduction zone. Our analysis considers solutions for both con ﬁ ned and uncon ﬁ ned aquiares. Resulting permeability and hydraulic di ﬀ usivity values range from 6 : 4 × 10 − 16 – 8 : 4 × 10 − 14 m 2 and 1 × 10 − 4 – 9 × 10 − 1 m 2 s − 1 , respectively, whereas speci ﬁ c storage values are generally ~ 1 × 10 − 6 m − 1 . The values are fairly consistent through time, reasonable given lithology, and are comparable to other regional studies. For one borehole, values are also comparable to those determined with traditional aquifer test data. In contrast with previous determinations of the high-frequency poroelastic response to seismic waves, no obvious spatial trends in hydrogeologic properties determined from long-wavelength tidal perturbations are observed. Within the recurring time-series estimates, only one borehole exhibits clear permeability enhancement by earthquakes, whereas nearby boreholes with similar lithology and hydrogeologic property values do not. This highlights the variable susceptibility of rocks to permeability enhancement. Together, these results provide quantitative constraints useful for models of large-scale groundwater ﬂ ow around large fault systems and the potential hydrologic in ﬂ uence on deformation and fault slip behavior. local geology.


Introduction
Pore fluid pressures within rocks and sediments have a direct impact on the slip behavior of faults and deformation processes through their influence on effective stress [1,2]. The subsurface distribution of pore fluid pressures and how they may change over time is largely controlled by the hydrogeologic properties of the rocks and sediments [3,4]. Characterizations of hydrogeologic properties within fault systems are therefore particularly of interest in efforts to understand processes controlling earthquakes and deformation [4].
Hydrogeologic properties control the distribution of pore fluid pressures through the pore pressure diffusion equation: where P in units of pascals is excess pressure defined as the excess formation pore pressure above hydrostatic equilibrium conditions [5]. Permeability, k, with units of length squared (m 2 ), describes how easily water flows through rock in the presence of a gradient in excess pressure. Specific storage, S s , with units of inverse length (m -1 ), describes the volume of water released by a volume of rock per unit change in hydraulic head, which is defined as the pore fluid pressure divided by gravity and fluid density. Table 1 lists and defines the variables used throughout this study. Specific storage is a function of the elastic properties of the formation matrix and fluid by where α and β are the matrix and fluid compressibilities, respectively, in units of Pa -1 , and n is porosity. Together, 2 Geofluids specific storage and permeability control the diffusion of excess pressure through hydraulic diffusivity with units of length squared per time, which is defined by the terms in front of the right-hand side of Equation (1) as A challenge in characterizing hydrogeologic properties, especially in tectonically active areas, is that the presence of faults and fracture networks at multiple scales limit what can be determined from laboratory-derived estimates on rock [6][7][8]. In situ measurements that interrogate large volumes are therefore necessary to characterize the controls on large-scale fluid flow and pressure distribution. Thus, in situ measurements of hydrogeologic properties are typically acquired through aquifer tests that require active perturbations to the subsurface hydrogeology. These tests generally involve adding or removing water within a borehole and monitoring the hydraulic recovery [9][10][11][12].
Alternatively, hydrogeologic properties can also be determined passively by monitoring the borehole pressure response to subsurface perturbations caused by solid earth tides [13][14][15]. This method requires longer monitoring times on the order of months to years, but with modern datalogging equipment, allows for many boreholes to be assessed over a large network simultaneously, and does not require active disturbance on the part of the investigator. In addition, borehole tidal response analysis allows for time series of hydrogeologic property estimates to be made. Hydrogeologic properties are not necessarily constants through time, and permeability enhancement and slow recovery from damage and healing processes have been shown to occur in response to earthquakes [13,14,16].
Here, we use borehole tidal response analysis to compute hydrogeologic properties of boreholes along the San Andreas fault system and Cascadia subduction margin from Vancouver Island in the north to the Salton Sea in the south ( Figure 1). These locations allow us to constrain hydrogeologic properties within the shallow crust along tectonically active margins where there is particular interest for understanding the poroelastic response to earthquakes and interpreting deformation signals associated with fault [17][18][19]. Our analysis focuses on determining the hydrogeologic properties of each borehole, assessing spatial variations within regions, and investigating potential temporal changes in response to regional earthquakes. We then compare and contrast our results to similar studies in these regions, standard values based on lithologies, and a traditional aquifer test performed on one borehole in Washington state, B001.

Materials and Methods
2.1. Tidal Response Analysis. Tidal response analysis allows for the passive determination of hydrogeologic properties by monitoring borehole pore pressure or water level within a well or borehole. Just as the gravitational pull of the sun and the moon cause ocean tides, they also impose strains on the solid earth [20]. Earth tides are observable within gravity measurements and borehole strainmeters at various tidal frequencies. The four largest amplitude tidal forcings occur roughly diurnally, K 1 and O 1 , or semidiurnally, S 2 and M 2 ( Table 2).
The strain from these forcings can cause undrained pore pressure changes within the formation dependent upon the compressibilities of the fluid and rock matrix. Because it takes time for the water to flow from the formation into a borehole, the pressure response observed within a well or borehole generally has a time delay from the tidal strain dependent upon the borehole dimensions and surrounding hydrogeologic properties. Similarly, the amplitude of the borehole response depends upon the amplitude of the tidal forcing at a particular frequency and the hydrogeologic and poroelastic properties of the borehole and formation [21][22][23].
Borehole tidal responses are also dependent upon groundwater flow characteristics, such as whether the screened or open interval is within a confined or unconfined aquifer [21,22]. However, with known borehole design characteristics, theoretical or observed tidal strain estimates, and reasonable assumptions regarding groundwater flow scenarios, determinations of the amplitude response and phase lag of tidal signals within wells and boreholes can be used to estimate the hydrogeologic properties of the surrounding formation.

Flow Models.
We consider two hydrogeologic flow scenarios for our analysis which reflect different hydrogeologic conditions. The first assumes that the borehole screen interval is situated within a confined aquifer and thus flow to and from the borehole is dominated by horizontal fluid flow [21]. The second assumes an unconfined aquifer which incorporates a free surface boundary and vertical fluid flow [22].
Both solutions describe the borehole pressure response to a periodic volumetric strain in the formation in terms of amplitude ratio and phase lag. Similar to Allegre et al., 2016, we define the amplitude ratio A ′ between the borehole pore pressure and strain as where A Pp is the amplitude of the borehole pore pressure at a given frequency, and A ε is the corresponding amplitude of the strain. Phase lag η is defined as where φ Pp and φ ε are the observed phase of the borehole pore pressure signal and strain phase at a given frequency, respectively. By this definition, a negative phase lag occurs when the strain leads the pore pressure signal.

Horizontal Flow within a Confined Aquifer from Periodic
Forcing. The first flow solution we consider predicts the borehole tidal response to periodic forcings within the formation 3 Geofluids assuming a single, laterally extensive, homogenous, isotropic, confined aquifer (Equations (15) and (16)) [21].
Amplitude ratio, A ′ , and phase lag, η, are related to hydrogeologic properties and borehole dimensions by where E and F are defined as Here, ω is the frequency of periodic forcing, r c is the radius of the cased section, μ is the viscosity, ρ is the density, g is the gravitational acceleration on the surface of the earth, b w is the length of the screened interval, Kei and Ker are the imaginary and real parts of the zeroth-order Kelvin function, k is the permeability, and α ω is defined as This model has been successfully utilized to determine values of hydrogeologic properties from borehole tidal response analysis in many settings [13-15, 21, 24-26]. Figure 2(a) shows the range of phase lag solutions from Equation (7) as a function of specific storage values and permeability considering the M 2 tidal frequency and casing radius, r c , the screened interval radius, r s , and screen depth b w specific to borehole B001. Phase lags are represented by color with large values of phase lag in yellow. The figure panel illustrates how phase lags for this flow scenario are particularly dependent upon permeability, with larger phase lags corresponding to lower permeabilities. In this example, and for most of the boreholes in this study, observable phase lags at the M 2 frequency are able to resolve permeabilities between~5 × 10 −15 m 2 and~1 × 10 −13 m 2 . Similar to Figure 2(a), Figure 2(b) shows the corresponding solutions for amplitude ratio (Equation (6) Figure 1: Map of study area with boreholes indicated by red dot boreholes clustered into northern, central, and southern groupings. 4 Geofluids of specific storage and permeability. In this panel, amplitude ratios are shown by the color, again with higher values in yellow. Figure 2(b) illustrates that amplitude ratios are particularly sensitive to specific storage, which is controlled by the elastic properties of the formation (Equation (2)).

Flow within an Unconfined Aquifer from Periodic
Forcing. The second flow solution we consider describes the borehole pressure response expected from a periodic forcing within an unconfined aquifer (Equations 6.62 and 6.64) [22]. Although originally formulated to describe the response to a periodic surface load, this model is also applicable for periodic volumetric forcings and has been used in a number of previous borehole tidal response analysis studies to account for situations in which the fluid flow to the well or borehole has strong vertical flow components [15,24,27]. For this scenario, the solution for amplitude ratio is where S s is specific storage, z is the depth from the drained pore pressure boundary, and δ is defined as   where D is the hydraulic diffusivity and ω is the angular frequency of oscillation.
The corresponding phase lag solution is which predicts phase lags ranging from -1 to 45 degrees [22]. The positive phase lags, or phase leads, result from the poroelastic response with the free surface and retains causality in terms of the increment of fluid mass, which is the dependent variable rather than borehole pressure. Similar to panels a and b of Figure 2, panels c and d show the range of solutions for phase lag and amplitude ratio as a function of permeability and specific storage for an unconfined aquifer considering borehole dimensions specific to B001. For the unconfined aquifer model, since phase lags and amplitude ratios are dependent on hydraulic diffusivity, which is a function of both permeability and specific storage (Equations (11) and (13)), the contours of equal phase lag and amplitude ratio run diagonally in Figures 2(c) and 2(d) [22]. However, unique values of permeability and specific storage can be resolved from inversions using a particular amplitude ratio/phase lag pair.

PBO Network/Data
Used. For our analysis, we use data collected from the borehole array component of the NSF Earthscope's Plate Boundary Observatory (PBO) network. We analyze data from 17 of the 23 boreholes, which have long time series data over many years without large data gaps or considerable data quality issues. These boreholes are situated on-land from Vancouver Island in the North to the Salton Sea in the South ( Figure 1).
The boreholes are constructed with 2 ″ PVC casing and communicate with formation fluid pressures through a screened interval typically 3-9 m long at depths ranging between 95 and 226 m below the surface. Each borehole is sealed at or near the top to isolate pressure measurements from and minimize the influence of atmospheric pressure and barometric fluctuations ( Figure 3; Table 3). The boreholes are each instrumented with Digiquartz depth sensors which record the borehole fluid pressure at 1 Hz frequency ( Figure 4). Barometric pressure sensors and rainfall meters are also located at each borehole site or within close proximity.
In contrast with most wells and boreholes in which borehole tidal response analysis has been performed, the PBO boreholes are also each instrumented with borehole Gladwin Tensor strainmeters that measure geologic strain every 300 s. Our analysis focuses on the time period from January 2009, when high-frequency borehole pore pressure sampling began, through December 2017.
Strain data from this network have been used to study earthquakes and fault slip behavior [28]. Pore pressure data have also been used in coordination with strain data to solve for poroelastic response variations in the region [29,30].
2.6. Constraints on Earth Tide Strain. Determination of hydrogeologic properties using the borehole response to earth tides requires knowledge of the amplitude and phase of the earth tide strain within the formation. Although strainmeter data are available for each borehole site in this study, direct observations of geologic strain are not commonly available. Instead, predicted values for a given location are commonly used [14,15,24,25,31] based on synthetic tidal loading models [32,33].
For comparison to the directly observed strain data, we also compute results using modeled strain generated using the software package Some Programs On Tidal Loading (SPOTL) [33]. The modeled strain time series in this study incorporate the effects of solid earth tides and the effects of ocean loading determined by global and local tidal models available for the US west, San Francisco Bay, and Gulf of California. Estimates of hydrogeologic properties calculated from both strainmeter data and modeled strain allow us to assess the implications of using modeled strain on hydrogeologic property estimates when direct observations are not available.

Isolating Tidal Signals from Barometric and Other
Biasing Effects. When analyzing signals for borehole tidal response analysis, it is important to ensure that the signals being interpreted are not greatly affected by processes beyond the tidal forcing. In order to avoid potential aliasing from short and long period signals outside of the diurnal and semidiurnal tidal range, we apply an acausal second order Butterworth filter with a bandpass from 0.8 to 2.2 cycles per day to both borehole pressure measurements and strain data [14].  6 Geofluids In addition, we consider the potential biasing influence of surface loading due to barometric pressure variations at similar diurnal and semidiurnal tidal frequencies [27]. Previous borehole tidal response analysis studies within wells and boreholes have largely been able to avoid the influence of barometric effects by analyzing the M 2 tidal frequency, which has a large earth tide strain amplitude, is often resolvable in borehole pore pressure or water level data, and typically does not appear as a strong component within barometric pressure [13-15, 21, 24, 27].
For our analysis, we assess the resolvability of tidal signals and the potential biasing effects of barometric pressure variations at various tidal periods by analyzing the Fourier amplitude spectra of both the borehole pore pressure and barometric pressure for each borehole ( Figure 5). From this analysis, we find that the barometric pressure data consistently have strong S 2 and K 1 signatures, which could bias the analysis, whereas M 2 and O 1 do not. In our subsequent analysis, we focus on the M 2 frequency and use time windows of around one month which allows us to adequately separate the M 2 and S 2 signals.
2.8. Tidal Amplitude and Phase Determination. In order to determine the phase and amplitude of the M 2 tidal component within the borehole pore pressure and strain data and calculate the amplitude ratio and phase lags, we simultaneously fit the four dominate frequencies within the filtered data at once using a least-squares fit. We do this analysis in a 29.6-day moving window with 80% overlap following previous studies [14,24]. The choice of a 29.6-day moving window allows for the semidiurnal M 2 and S 2 to be distinguished for one another as well as the diurnal O 1 and K 1 tidal constituents.
Next, using the relationship that where A = ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi a 2 + b 2 p and φ = atan2ð−b, a Þ, we can obtain the phase and amplitude of our tidal constituent signals using least-squares regression.
The kernel constructed by Equation (15) allows us to solve for a set of equations for multiple frequencies with the functional dependence of Equation (14): where ω 1 -ω n represent tidal frequencies of interest and t 1 t m represent the time vectors of interest. The model parameters m to be solved correspond with the constants in where f is a constant mean offset to the data, and a 1a n represent the a value and b 1b n represent the b value for each frequency of consideration. We then take the kernel G, model fit parameters m, and the data vector, d, and calculate the least-squares fit inversion by where m represents the fit properties and d represents the data vector of interest, either strain or water level data, and T reflects the transpose of a matrix. Error in fit, e, is defined as And the posterior total error, E, and standard deviation for each fit are calculated following the example of Menke and Menke, 2016 by where N is the length of d and M is the number of model parameters m [34].
, and (13) developed for each borehole specific properties (e.g., Figure 2). Although both amplitude ratio and phase lag are both dependent on permeability and specific storage (Equations (6), (7), (11), and (13), Figure 2), each pair defines a solution. This can be thought of as two contours of the range of solutions for amplitude ratio and phase lag that overlap and match to one solution of permeability and specific storage.
To assess the resolvability and potential uncertainty in our determinations of hydrogeologic properties, we propagate errors in phase lag and amplitude ratio by using the following relationship: where σ ktot is the total uncertainty in permeability due to the uncertainty in permeability from amplitude ratio, σ kA , and uncertainty in permeability from uncertainty in phase lag, σ kz . This error propagation relationship is used for both the low and high estimates of specific storage and permeability.
For B028, B073, B079, and B086, the observed pore pressures do not have to contain clear M 2 signals such that hydrogeologic property estimates are not able to be determined. There seems to be no clear distinction in lithology, borehole design, or geologic setting for these boreholes that could explain the lack of this tidal component. If S s is comparable to values at other boreholes, this lack of resolvable M 2 signal could be explained by low values of permeability for a confined aquifer setting (Figures 2(a) and 2(b)) or high values of permeability for an unconfined aquifer setting   9 Geofluids (Figures 2(c) and 2(d)), such that the poroelastic response to earth tides is not able to transmit formation pressure to the borehole within the tidal period or the response diffuses away too quickly due to the sensitivity of free surface, respectively. Alternatively, the lack of a strong M 2 signal could also reflect large specific storage values which result in lower formation pressure changes in response to changes in strain (Figures 2(b) and 2(d)).

Comparing Results from Modeled Strain versus
Observed Strain. Our analyses utilized both collocated strainmeter data which provides a direct measurement of borehole strain and modeled strain that is more commonly used in borehole tidal response analysis.
For a number of boreholes, particularly those in the Cascadia and Parkfield regions, the phase lags when using strainmeter data and modeled strain differ by as much as 38  Figure 6: Median values of hydrologic properties. Wells are arranged descending to the south. Inversion 1 uses filtered borehole pore pressure and filtered modeled strain for analysis. Inversion 2 uses filtered borehole pore pressure and filtered strainmeter data for analysis. Background colors represent regions determined by geographic locality and similarity of local geology. 10 Geofluids degrees ( Figure 7) which likely reflects limitations in the tidal models, perhaps because of local ocean loading within Puget Sound not fully considered within the model. In general, when using modeled strain near a body of water such as a large bay, follow-up with either observed strain or gravity measurements would provide a useful check to see if the model is reasonable in this location. Whereas constraints based on the directly observed strain data are ideal, the modeled strain can be used when these data are missing. For some boreholes, hydrogeologic properties were not able to be determined using the observed or modeled strain, or both, because the phase lags were outside of the predicted range of either flow model. This suggests that for these boreholes, in particular boreholes B005, B011, and B078, flow conditions are likely more complex than the end-member confined and unconfined scenarios considered. Where both observed and modeled strain can be used, the difference in solutions is visible as separation in the blue and red dots in Figure 6. For some boreholes such as B081, these differences in estimates of permeability can be nearly an order of magnitude. When accuracy is of particular importance, these results highlight the importance of independent constraints on the phase tidal strains.

Comparison to Traditional Aquifer
Testing. For one borehole, B001, data from traditional slug tests are available which allow us to determine estimates of permeability and specific storage independent of the borehole tidal response analysis determinations (Figure 8). A similar comparison between passive and active methods for in situ hydrogeologic property determination was performed by Allegre et al. (2016), which found comparable estimates using the two methods.

Geofluids
For a slug test, a volume or "slug" of water is rapidly introduced into a borehole with a starting hydraulic head of H 0 and at a later time withdrawn, and the hydraulic head, H, is monitored for both tests where hydraulic head H is defined by Equation (22).
We interpret the slug test data using the Cooper et al. [11] solution: where γ is defined as and ϵ is defined as where t is the value of time starting from the initial perturbation. It should be noted that for both slug tests, the normalized head does not recover fully and curve matching is more ambiguous and limiting the confidence in the resulting property estimates (Figure 8).
For the slug insertion test, the resulting estimates are a permeability of 5:0 × 10 −15 m 2 and specific storage of 1:9 × 10 −6 m −1 . For the slug withdraw test, the estimated value of permeability is 1:5 × 10 −15 m 2 and 1:9 × 10 −6 m −1 for specific storage once again. When compared with the median specific storage value of 6:1 × 10 −7 m −1 from the borehole tidal analysis using strainmeter data, the estimates closely agree. When compared with the median permeability of 5:0 × 10 −14 m 2 obtained from borehole tidal response analysis for this borehole (Table 4), these results are comparable; however, lower permeabilities are obtained using the slug test. This may reflect differences in the volume of formation that is being interrogated by the different methods. Alternatively, as Allegre et al., 2016, notes, the variety of models used for traditional aquifer testing can yield variations of more than an order of magnitude in property determinations, further suggesting that the results between methods are comparable [15].

Temporal Changes in Hydrogeologic Properties.
Because of the continuous cycle of the earth tide forcing over time, borehole tidal response analysis can allow for repeated determinations of hydrogeologic properties and observation of their changes over time (Figures 9 and 10). Although traditional aquifer tests have been used to assess permeability evolve over time in an active fault zone from sparse repeat measurements [35][36][37], borehole tidal response analysis has allowed for long time-series estimates of hydrogeologic properties which have revealed temporal changes in permeability in response to local and distant earthquakes both inside fault zones and in regular country rock [14,25,26,38].
Earthquake-related changes in hydrogeologic systems from earthquakes have long been observed and studied, including changes in water levels in wells [13,[39][40][41], increased discharge of streams and springs [42][43][44], generation of seiches [45,46], and elevated fluid flow rates through faults [4,47,48]. With our time series data for permeability, we are able to look for enhancements in permeability that may be related to earthquake-induced damage. Changes in water level or discharge after earthquakes can result from multiple processes, including the poroelastic response to static or dynamic stresses and strains, and do not require changes in permeability or hydrogeologic properties [17,18]. Interestingly, Lai et al. [25] find that many wells around the epicenter of the 2008 Wenchuan earthquake exhibited increased permeability afterwards, irrespective of whether water level changes went up or down from poroelastic strain change. Permeability enhancement is of interest, however, because it can greatly affect fluid flow and pore pressure distribution. Because permeability of rocks can range over several orders of magnitude and is highly sensitive to pore and fracture connectivity, changes due to damage and disturbance are likely resolvable within time series data [49,50]. Additionally, permeability at different spatial locations will likely vary, as the lithologies and fracture densities are varied [5]. Conversely, specific storage is not expected to vary over time or across sites, as even though some rocks may be more frac-tured; the compressibility of rock matrix, α, (10 -8 -~10 -10 Pa -1 ) has a narrow range of values among consolidated materials and is relatively less sensitive to damage [5].
Of the 17 boreholes analyzed, only B084 shows clear changes in permeability in response to earthquakes. Six earthquakes produced changes in permeability, including the August 3 rd , 2009, M w 6.9 Gulf of California Earthquake, the April 4 th , 2010, M w 7.2 El Mayor-Cucapah, and smaller, more local earthquakes. While a change in permeability would likely be instantaneous in nature, our moving window of 29.6 days used for analysis causes a convolution of the signal and a gradual rise over a month, as seen in Figure 11. The 13 Geofluids permeability then appears to decrease back to the long period trend after a period of weeks to months.
To evaluate the sensitivity of permeability from dynamic shaking, we analyze earthquakes with seismic energy densities, e, estimated to be greater than 10 -3 J m -3 [44]. Seismic energy density is calculated using the empirical relationship: where r is the distance between the earthquake and borehole and M w is the magnitude of the earthquake. The earthquakes that meet these criteria are plotted overlaying the permeabil-ity and specific storage time series data (Figures 9 and 10). To aid in this, we apply a 10 s seventh order Butterworth high pass filter to the borehole pore pressure data, as earthquakes will appear as spikes ( Figure 11) [29].

Analysis and Discussion
3.2.1. Hydrogeologic Properties. We compare our values of hydrogeologic properties obtained from tidal response to standard values for the lithology surrounding to check for consistency [5,51]. 14 Geofluids to 6:3 × 10 −14 m 2 . Values lower than 1 × 10 −14 m 2 for a permeable basalt are unusual [5], but for B022 in Mendocino, the permeable basalt overlays a sandstone unit within 10 m of the screened interval, complicating interpretation. For boreholes that are surrounded by sandstones, B022 and B078, permeabilities of 2 × 10 −14 m 2 and 6:7 × 10 −14 m 2 are consistent with standard values for these lithologies. Finally, all other boreholes analyzed are set in granites and granitic rocks along the San Andreas fault system, ranging from 1:7 × 10 −15 m 2 to 8:6 × 10 −14 m 2 , within standard values for fractured igneous and metamorphic rocks. We see diffusivity values vary over more than three orders of magnitude across systems, likely due to differences in lithology [5]. Since most values of specific storage are similar, the differences in diffusivity largely result from differences in permeability. Prior analysis by Barbour [29] has investigated the poroelastic response of these boreholes to dynamic strains associated with seismic waves. Barbour describes the relationship between areal strain, E a , and pore pressure, P, by the relationship [29] where the variable d that allows for responses to be scaled by drainage where 1 is undrained and 0 is drained, and c, which is comparable to Skempton's coefficient times bulk shear modulus when d equals 1, describes the undrained pore pressure response to strain. Because d relates to drainage and pressure diffusion, we investigate whether there is a correlation with our permeability and hydraulic diffusivity determinations (Figure 12). Similarly, since c is a measure of the pore pressure response to strain, we investigate any correlation with ρg/S s which has the same units (Pa), but whose inverse describes the volumetric strain due to a pressure change ( Figure 12). The results, however, do not show any strong correlations between variables. We see from panel a in Figure 12 that ρg/S s values seem to slightly increase with increasing c value but this correlation is weak, meaning that c cannot be thought of as equivalent to S s or not explained by elastic properties alone. Additionally in panels b and c, we see that permeability and diffusivity appear to decrease as d increases, a counter-intuitive result, as one might expect that the higher the permeability and diffusivity, the more easily fluids would be allowed to drain through the system, indicating a higher d value. Inversion 1 Inversion 2 Figure 11: Raw borehole pore pressure, high-pass filtered borehole pore pressure, and permeability plotted to look for changes in permeability caused by earthquakes. Earthquakes with a seismic energy density above 10 -3 J m -3 are plotted as magenta dashed lines. Inversion 1 uses filtered borehole pore pressure and filtered modeled strain for analysis. Inversion 2 uses filtered borehole pore pressure and filtered strainmeter data for analysis.

Geofluids
Barbour [29] also found spatial variations in the c and d for these boreholes as a function of distance from large faults in the Parkfield and Anza clusters. Similarly, for a series of boreholes in central California north of Parkfield immediately adjacent to the San Andreas Fault, Xue et al. [24] found increased permeability and decreased specific storage nearest the fault such that diffusivity was relatively constant. One would expect that increased fracture density near faults may lead to higher permeability and perhaps lower specific storage.
In panels d though f of Figure 12, we plot our determinations of permeability, specific storage, and hydraulic diffusivity for the California boreholes as a function of distance from the nearest large fault within the San Andreas Fault system. No clear trends are observed, except perhaps a slight decrease in specific storage with increasing distance, which is counter to the hypothesis. The lack of clear trends may suggest that lithologic differences may have a greater influence on differences between boreholes than fracture density and fault proximity. When compared to the trends observed for c and d, the difference likely reflects the large frequency difference between seismic waves and the M 2 earth tide, the corresponding differences in the scale of volume being measured and flow processes involved.

Temporal Changes in Hydrogeologic Properties.
Our results show only one borehole, B084 in southern Califor-nia, has clear permeability enhancements from earthquakes. When earthquakes that generated a response in permeability are plotted with those that did not in Figure 13, it appears as if not only did larger, farther away events affect permeability, but as did some nearby, smaller magnitude events. Additionally, the aftershocks of the larger events of sizable seismic energy density may also increase permeability, but it is impossible to distinguish these from the effects of the main shock. B084 is surrounded by a lithology of granite, but nearby boreholes B081, B086, and B087 are as well, so the difference in responsivity cannot be explained by a difference in lithology.
It is noteworthy that B084 is the only borehole that reliably exhibits permeability enhancements from dynamic shaking, suggesting that some boreholes may be more responsive to dynamic shaking than others. Barbour [29] notes that boreholes farther from the fault tend to exhibit more of a response in dynamic pore pressure to dynamic shaking caused by earthquakes, but B088 is farther away from the fault then B084 by about 6 km and does not exhibit permeability enhancements [29]. Additionally, it is observed that some earthquakes over the minimum seismic energy density caused no increase in permeability. Further investigations are needed into the factors influencing why some formations and boreholes are more susceptible to permeability enhancement due to dynamic shaking and the controls on variations in healing time.

Conclusions
Our analysis of the borehole pore pressure response to earth tides provides time-series estimates of permeability, specific storage, and hydraulic diffusivity for 17 boreholes in the NSF's Plate Boundary Observatory (PBO) network along the San Andreas fault and Cascadia subduction zone margins between 2009 and 2018. Median permeability values over time range from 6:4 × 10 −16 m 2 to 6:2 × 10 −14 m 2 ; specific storage is relatively similar among boreholes, with a value of~1 × 10 −6 m −1 , and shows little variation over time.
Values of hydrogeologic properties are generally consistent through time and are reasonable given lithologies, comparable to other studies, and for one borehole, B001, in Washington state, comparable to a traditional aquifer test. When using modeled strain for analysis, changes in properties are visible, but observed strain or other verifications of tidal phases are necessary to best constrain the overall magnitude of hydrogeologic properties. B084 in southern California contains changes in permeability that are comparable with damage and healing processes following some earthquakes with sufficient dynamic shaking. Although B084 shows evidence of permeability enhancement to earthquakes, other boreholes in the vicinity do not. It is not clear why this borehole responds differently, as B084 contains similar lithologies and hydrogeologic conditions to surrounding boreholes, highlighting the uncertainties in understanding the variable susceptibility to permeability enhancement from earthquakes. Results show no obvious trends with distance from large faults, in contrast with values of poroelastic properties previously constrained by analyzing the pressure response to seismic waves, but are consistent among regions and sides of a fault. Overall, our results provide constraints on hydrogeologic properties that can be useful for efforts to understand hydrologic influence on deformation and fault slip behavior.

Data Availability
Pore pressure and areal strainmeter data are available at https://www.unavco.org/data/data.html and environmental data are available at https://service.iris.edu/irisws/ timeseries/docs/1/builder/ with network PB, station number, and channel codes RRO for rainfall and LDO for atmospheric data, respectively.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.  Figure 13: (a) Catalogued earthquakes for B084. Seismic energy density plotted as contours, with all earthquakes larger than 10 -3 J m -3 plotted as larger markers. Red events represent main shocks that caused changes in permeability; black events are aftershocks from these events, while blue events caused no change in permeability. (b) Magnitude of events over time. Red events represent main shocks that caused changes in permeability; black events are aftershocks from these events, while blue events caused no change in permeability.