Attenuation of Seismic Waves in Partially Saturated Berea Sandstone as a Function of Frequency and Confining Pressure

Frequency-dependent attenuation (1/Q) should be used as a seismic attribute to improve the accuracy of seismic methods and imaging of the subsurface. In rocks, 1/Q is highly sensitive to the presence of saturating fluids. Thus, 1/Q could be crucial to monitor volcanic and hydrothermal domains and to explore hydrocarbon and water reservoirs. The experimental determination of seismic and teleseismic attenuation (i.e., for frequencies < 100 Hz) is challenging, and as a consequence, 1/Q is still uncertain for a broad range of lithologies and experimental conditions. Moreover, the physics of elastic energy absorption (i.e., 1/Q) is often poorly constrained and understood. Here, we provide a series of measurements of seismic wave attenuation and dynamic Young’s modulus for dry and partially saturated Berea sandstone in the 1–100 Hz bandwidth and for confining pressure ranging between 0 and 20 MPa. We present systematic relationships between the frequency-dependent 1/Q and the liquid saturation, and the confining pressure. Data in the seismic bandwidth are compared to phenomenological models, ultrasonic elastic properties and theoretical models for wave-induced-fluid-flow (i.e., squirt-flow and patchy-saturation). The analysis suggests that the observed frequency-dependent attenuation is caused by wave-induced-fluid-flow but also that the physics behind this attenuation mechanism is not yet fully determined. We also show, that as predicted by wave-induced-fluid-flow theories, attenuation is strongly dependent on confining pressure. Our results can help to interpret data for near-surface geophysics to improve the imaging of the subsurface.

Laboratory studies offer the opportunity to understand the physics of attenuation mechanisms and to validate theories (e.g., Tisato and Quintal, 2013). However, the few lithologies and physical conditions investigated to date are insufficient to achieve such a goal, and some studies report measurements which are not yet fully understood (e.g., Tisato and Madonna, 2012).
The present contribution reports a series of seismic wave attenuation measurements in extensional mode (Green et al., 1967(Green et al., , 1974a as a function of frequency and confining pressure for two partially saturated Berea sandstone specimens. Samples were saturated with water, glycerin, and air. Saturation levels were chosen to mimic different saturation conditions as proxies of subsurface reservoir conditions. Moreover, an additional sample is used to estimate pressure-dependent ultrasonic wave velocities for dry Berea sandstone. We employ phenomenological models to fit the dataset confirming the relation of causality between modulus dispersion and attenuation in our measurements (i.e., Kramers-Kronig relation). Further, squirt-flow and patchy-saturation models and ultrasonic wave speed measurements indicate that WIFF can cause the observed frequency-dependent seismic wave attenuation.

Samples and Measurements
We present measurements of physical properties for three Berea sandstone samples named: BS2, BS3, and BS5s. The three cylindrical samples were acquired from Cleveland Quarries 1 and had dimensions (length x diameter) of 250 x 75.1 mm, 250 x 75.1 mm, and 59 x 25.6 mm, respectively. The samples' length and diameter were measured with an uncertainty of ± 0.025 mm using a 0.01 mm division caliper. Sample BS3 has hydraulic permeability of ∼10 −12 m 2 , while BS2 and BS5s, which originated from the same block, have hydraulic permeability of ∼0.75x10 −12 m 2 . Permeability values were provided by Cleveland Quarries. Berea sandstone contains bedding planes controlling elastic properties anisotropy (Lo et al., 1986;Sayers et al., 1990). The axes of our cylindrical samples were perpendicular to such bedding planes.
For samples BS2 and BS3, we measured the frequencydependent Young's modulus (E) and seismic wave attenuation, i.e., 1/Q E that here we call 1/Q for the sake of simplicity. Measurements were performed at room temperature (∼20 • C) in the bandwidth 1-100 Hz and at confining pressures (P c ) ranging 0 -20 MPa. Samples were oven-dried at 105 • C for more than 24 h, then let cool at room temperature and humidity for several hours, jacketed for testing, and finally loaded inside the testing apparatus. Data are presented for dry and partially saturated samples. Under the aforementioned conditions, we also measured the ultrasonic longitudinal wave speeds of sample BS3. Such measurements were performed at a frequency (f ) of ∼100 kHz.
Sample BS5s was smaller than BS2 and BS3, and used to estimate the longitudinal and transverse ultrasonic (f ∼1 MHz) velocities of dry samples. The density (ρ) of the dry samples is 2110 ± 31 kg/m 3 , and it was obtained from the sample volume and weight, which was measured with a 0.1 g resolution scale.
Porosity ( ) of sample BS5s is 0.205, and it was estimated using a helium pycnometer (AccuPyc 1330, Micromeritics). We use such a porosity value as a proxy for the porosity of BS2 and BS3.
We saturated the samples with air, water, and a waterglycerin mixture. The latter is called GWmix and contains 55 and 45 volume % of glycerin and distilled water, respectively. In general, the subscript W and GW refer to water and GWmix, respectively. GWmix viscosity (η GW = 11 ± 0.25 mPa s) and density (ρ GW = 1.12 ± 0.02 g/cm 3 ) were measured with a falling ball viscometer, a 0.1 mg precision scale and a 0.1 ml precision vial, respectively. The bulk modulus of GWmix (K GW = 2.4 ± 0.1 GPa) was estimated by analyzing the pressurevolume curve obtained with a high-pressure syringe pump (ISCO 260D). Sample BS2 was saturated with air and water, while BS3 was saturated with air and GWmix.
While inside the pressure vessel, water or GWmix (i.e., liquid) was pumped into the sample through the sample bottomend-face, letting air flow through the sample top-end-face. Saturation was calculated according to the quantity of injected liquid and the porosity of the sample. The injected liquid was drained from a graduated burette, and any excess liquid flowing through the top end-faces was collected on a second graduated burette. Thus, the volume of injected liquid was estimated as the difference between the liquid levels in the two graduated burettes. Some of the injected fluid was used to fill the hydraulic circuit and was not considered in the saturating fluid volume. In particular, before recording the initial volume of fluid in the burettes, we prefilled the injection circuit except for the bottom sample holder to avoid sample imbibition. We also estimated the internal volumes of the sample holders and draining circuit with an accuracy of +/−1 ml representing less than 0.5% of BS2 or BS3 pore space volume. Once the wanted saturation was reached, the sample was isolated from the hydraulic circuit by closing the inlet and outlet valve. All the experiments were performed with fluid pressure equal to room pressure (∼0.1 MPa), and P c represents a proxy for effective stress (σ').
Each saturated sample effective density (ρr) was estimated considering the rock porosity, saturation, and fluid density.

Low-Frequency Testing
Attenuation and complex Young's modulus in the bandwidth 1 -100 Hz are obtained employing the sub-resonance method (McKavanagh and Stacey, 1974). Attenuation is: where ϕ is the phase shift between the stress (σ) imposed by the seismic wave and the strain (ε) measured across the whole sample.
Q is the quality factor. The σ to ε ratio equals the complex Young's  modulus (E z ), so that 1/Q is also: The real part of E z [i.e., Re(E z )] or simply E SF is given by: where SF stands for "seismic-frequencies, " σ pp and ε pp are the peak-to-peak values of the σ and ε, respectively (Lakes, 2009). ε pp was kept around 1.2x10 −6 , which is similar to the strain caused by the propagation of a seismic wave (Karato and Spetzler, 1990). The extensional seismic wave attenuation (1/Q) and complex Young's modulus (E Z ) of sample BS2 and BS3 were measured using the Burlini attenuation vessel (BBAV) (Figure 1; Tisato and Madonna, 2012). We calculated i) the stress from the load-cell measurement (i.e., force) and the sample cross-section area, and ii) the strain from the displacement sensor measurement (i.e., sample shortening) and sample length. Force and shortening time series are fit with sinusoidal functions. Amplitude and phase of the fitting functions yield the parameters to calculate 1/Q and E SF . The fitting functions are assumed to represent the true values. We calculate two standard deviations between the true values and the measurements to estimate the 95% confidence bounds on both the amplitude and the phase of the signals. Such uncertainties and those related to the sample cross area and length provide a proxy for accuracy on Young's modulus and attenuation. Precision (i.e., repeatability) is estimated by calculating the standard deviation of five measurements performed at each frequency. Accuracy and precision are summed to provide uncertainty on the measurements that are provided as error bars in Figures 2-8.
The curved face of the samples was covered with an aluminum foil, which was glued onto the rock, and a jacket made of Fluorinated Ethylene Propylene (FEP) shrink-tubing. These two layers do not affect the stiffness of the sample, isolate the sample from the confining medium (i.e., silicon oil) and avoid that the curved face of the sample acts as a free-flow boundary by limiting the saturating fluid to escape radially while the sample is dynamically stressed (Gardner, 1962;Dunn, 1986) (see Supplementary Material). In addition, the specimen endfaces were machined on a grinder to obtain a parallelism tolerance of ∼40 µm to meet the required experimental accuracy (Paffenholz and Burkhardt, 1989).
A uniaxial press, which is driven by a servo-electric uniaxial actuator, hosts the BBAV and applies the vertical stress on the sample (σ 1 ). During testing, the stress field applied to the specimen was: σ 1 = P c + 2.8 MPa and σ 2 = σ 3 = P c . The vertical stress (σ 1 ) was intentionally kept 2.8 MPa higher than the confining pressure to ensure good coupling between the horizontal interfaces of the sample assembly.  (Liu et al., 1976). The low-frequency seismic wave was generated as a sinusoidal variation (σ) of the vertical stress. A piezoelectric motor, which is controlled by a high voltage amplifier, generates such a sinusoidal variation of σ 1 that is necessary to estimate the complex Young's modulus and attenuation. σ and ε are estimated for the measurements performed using a load cell and a cantilever, respectively. Calibration and additional details about the BBAV apparatus are reported in Tisato and Madonna (2012).

Ultrasonic Testing
Longitudinal (V P ) and transverse (V S ) ultrasonic velocities were obtained using the pulse transmission method (Birch, 1960). V P of BS3 was measured during low-frequency testing inside the BBAV, and the stress field applied to the specimen was: σ 1 = P c + 2.8 MPa and σ 2 = σ 3 = P c . Uncertainties were estimated combining: (i) the uncertainty on picking the first arrival, (ii) the uncertainty on the delay introduced by the apparatus (26.1 ± 2 µs), (iii) the uncertainty on the sample length (1 mm), and (iv) the variability between two consecutive sets of measurements. The two sets of measurements were respectively collected for consecutive increasing and then decreasing P c. For each confining pressure, the velocity was computed as the average between such two measurements. For each seismogram, the first arrival uncertainty was assessed by picking a minimum (t min ) and a maximum (t max ) time value to the first arrival. Such values were always around the wavelet onset, and they typically differ of ∼1 µs. The first arrival and its uncertainty were then calculated as t min + t max 2 and t min −t max 2 , respectively.
On the other hand, V P and V S of BS5s were measured employing another pressure vessel named Geneva-Rig in which the stress field applied to the specimen was hydrostatic, i.e., σ 1 = σ 2 = σ 3 = P c (Zappone et al., 2000;Kästner et al., 2020). We calculated uncertainties combining: (i) the uncertainty on picking the P-and S-wave arrivals, (ii) the uncertainties on the delay introduced by the apparatus (∼1.06 ± 0.4 µs and ∼10.5 ± 0.8 µs for P and S-waves, respectively) and (iii) the uncertainty on the sample length. Two tests were conducted at increasing and decreasing confining pressures between 1 and 40 MPa. The results, i.e., velocities vs confining pressure, were fit with empirical relations that are described by the following equation: where A, B, C, and D are the best-fit parameters estimated through least square minimization (Eberhart-Phillips et al., 1989). Ultrasonic shear (µ HF ) and bulk (K HF ) moduli were estimated as: where HF stands for "high-frequency" and ρ r is the effective sample density. For dry conditions, we refer to the ultrasonic bulk and shear modulus as: K HF(Sa = 1) and µ HF(Sa = 1) , respectively; where S a is the gas (air) saturation. Then, ultrasonic Young's modulus (E HF ) was estimated as: (Mavko et al., 2009).

Phenomenological Models
Modulus dispersion and attenuation are linked through the Kramers-Kronig relation or causality (Lakes, 2009). Experimental data have been previously verified through the Kramers-Kronig relation (e.g., Mikhaltsevitch et al., 2016). Here, to verify that our results honor causality, we use two phenomenological models that obey the Kramers-Kronig relation. For frequencyindependent measurements, we utilize the nearly constant Q model (NCM) (Liu et al., 1976). For frequency-dependent measurements, we use the sum of the NCM and the standard linear solid model (SLS) (Carcione, 2007) (see Supplementary Material). We assume that the sum of attenuation is an appropriate operation as Johnston et al. (1979) and Tisato and Quintal (2013) show that the frame-related attenuation and the fluid-related attenuation are independent and can be summed to obtain the total attenuation of the saturated sample.
In particular, we assume that the frame related attenuation is constant, regardless of the saturation level, and the fluid-related attenuation is a function of saturation and can be represented by the SLS rheology. Such fitting is possible as wave-induced fluid flow mechanisms and SLS produce a low-frequency asymptote that scale as Q −1 ∝ f (Mueller et al., 2010).  (Liu et al., 1976;Carcione, 2007).  (Liu et al., 1976;Carcione, 2007).
where κ is the rock permeability, K s is the bulk modulus of the saturated rock, L is the characteristic size of the patches (White, 1975). For the present experiments, we can assume κ = 10 −12 m 2 , L ranging between 1 and 10 cm, η f ranging between η W and η GW and K s ranging between 10 and 20 GPa. The selected range of K s agrees with the literature, and the ultrasonic measurements will later confirm it. According to this parameter space, f patchy ranges between ∼28 Hz and ∼63 kHz. Jones (1986) provides a formula for the squirt-flow characteristic frequency, which is: where α is the characteristic crack aspect ratio for the rock, and K g is the bulk modulus of the rock frame. For the present experiments, we assume α ranging between 10 −4 and 10 −3 (Zimmerman, 1991), η f ranging between η W and η GW and K g = 37 GPa, as Berea sandstone is mainly composed of quartz. According to this parameter space, f squirt ranges between ∼3 Hz and 37 kHz. The two wave-induced-fluid-flow mechanisms operate in the bandwidth of investigation (i.e., 1-100 Hz); therefore, both are considered as potential attenuation mechanisms in our experiments. We also rule out the possibility that the observed frequency-dependent attenuation is caused by global flow as the characteristic frequency of such attenuation mechanism ranges above 30 kHz, which is well above the upper limit of the investigated bandwidth (Biot, 1956).
To further investigate the occurrence of wave-inducedfluid-flow we compared ultrasonic and low-frequency elastic properties to effective medium theories: namely Biot-Gassmann-Hill (Johnson, 2001) and Gassmann fluid substitution (Gassmann, 1951), respectively; And, we also compared the phenomenological model fittings with analytical solutions for squirt-flow and patchy-saturation (White, 1975;Dutta and Odé, 1979;Gurevich et al., 2010).

Analytical Solution for Squirt-Flow
We use the model proposed by Gurevich et al. (2010) to estimate complex frequency-dependent Young's modulus due to squirtflow. Such a model deals with fully saturated rocks while our samples are partially saturated. Thus, we assume that squirt-flow operates only in the liquid saturated portion of the sample and that the dry portion of the sample presents negligible attenuation. Under this assumption, we approximate the imaginary part of the sample Young's modulus by multiplying the imaginary part of the calculated Young's modulus and the saturation. Attenuation is then computed as the ratio between the imaginary and the real part of Young's modulus. The procedure and the equations that are used to model squirt-flow are reported in the Supplementary Material.
The model requires the knowledge of: (i) Bulk and shear modulus of the dry rock; (ii) Bulk and shear modulus of the mineral making up the rock (K g and µ g ); (iii) Crack porosity ( c ); and (iv) Mean crack aspect ratio (α). For the saturating fluid, the required properties are: (i) Density (ρ w or ρ GW ); (ii) Viscosity (η w or η GW ); and (iii) Bulk modulus (K w or K GW ). In addition, the model needs the effective pressure at which the cracks can be considered closed (P h ) and the related rock bulk modulus (K h ).

Analytical Solution for Patchy-Saturation
Complex Young's modulus and attenuation related to patchysaturation are calculated according to the model proposed by White (1975) and Dutta and Odé (1979). To model patchy saturation, we use the Matlab code "patchw" provided by Mavko et al. (2009), and that can be downloaded at https://pangea.stanford.edu/departments/geophysics/dropbox/ SRB/public/data/RPHtools.htm.
The model assumes that the rock is saturated with two fluid phases, i.e., liquid and gas. In addition to the petrophysical parameters that are listed in the previous paragraph, and with the exception of c , α, P h , and K h , the model requires (i) Rock porosity ( ); (ii) Hydraulic permeability (κ); (iii) Gas density (ρ a ); (iv) Gas viscosity (η a ); (v) Gas bulk modulus (K a ); and (vi) The radius representing the average size of the liquid patches that are surrounded by air (r p ).
The portion of pore volume that was not saturated with water or GWmix was saturated with air. Thus, we express air saturation as: Where ∨ means or. We present Young's modulus and attenuation values for case i through v. For each case we report curves for eight different confining pressures ranging between 0 and 20 MPa. Each curve consists of 30 to 36 measurements taken at frequencies between 1 and 100 Hz. Therefore, we present ∼1440 measurements of Young's modulus and attenuation. Each attenuation measurement is calculated as the average of five measurements acquired at the same frequency. The standard deviation of the five repetitions is used to estimate precision that is shown as half-height of the error bar.
For dry Berea sandstone (case i), attenuation is frequencyindependent, and around 0.01, and Young's modulus increases  (Liu et al., 1976). (B) Second step of the fitting procedure. The residuals between the NCM fitting and the measurements are interpolated with the standard linear solid model (SLS) (Carcione, 2007). Finally, NCM and SLS fittings are summed.
logarithmically with frequency (Figure 2). Data are in agreement with the NCM model (Figure 3).
For each test of case i, we averaged Young's modulus and the attenuation over the frequency bandwidth. As a general trend, we observe that as confining pressure increases, the averaged E increases from ∼17 to ∼26 GPa, and 1/Q decreases from ∼0.01 to ∼0.005.
Besides, measurements performed on partially saturated samples (case ii, iii, iv, and v) show two distinct regimes: a. At confining pressures above 14 MPa, attenuation is frequency-independent with values similar to those measured in case I (i.e., at dry conditions). Moreover, Young's modulus increases logarithmically with frequency, which is in agreement with the NCM model; b. At confining pressures below 14 MPa, attenuation and Young's modulus are frequency-dependent (Figures 4-7) and can fit the sum of the NCM and SLS model (Figure 8).
SLS fitting provides an estimate for corner frequency (f c ), i.e., the frequency at which attenuation is maximum. Such corner frequencies range between 114 and 579 Hz ( Table 2). For example, for case ii (S W = 60%) and confining pressure 8 MPa, 1/Q increases from ∼0.01 to ∼0.08 as frequency increases from 1 to 100 Hz (Figure 4 curve 5). Similarly, measurements conducted on the samples saturated with GWmix (case iii, iv and v) are frequency-dependent, and 1/Q varies with confining pressure. For case iii (S GW = 53.5%), 1/Q values are statistically lower than those measured for case iv and v (S GW = 80.2% and 86.4%, see Such values are reported in Figure 9, highlighting hoŵ E increases with confining pressure and how dispersion and attenuation decrease for confining pressures > 12 MPa. Longitudinal (V p ) and transverse (V s ) ultrasonic velocities of dry BS5s increase with confining pressure (0 -40 MPa) from 2426 to 3752 m/s and from 1880 to 2429 m/s, respectively. Curves that show velocities as a function of confining pressure are fit according to Eq. 4, yielding the following two equations: V S m s = 2413 + 0.86 P C (MPa) − 533 e −0.08 P C (MPa) , (12) V p was also measured for sample BS3. Figure 10 shows velocities curves as a function of confining pressure for different cases. Saturated BS3 exhibits higher V p than those measured for dry BS3. For example, for case v (S GW = 86.4%) and vi (S GW = 35.6%), and confining pressure comprised between 2 and 25 MPa, V p ranges between 3328 and 3822 m/s and between 3073 and 3652 m/s, respectively. These results agree well with the literature (e.g, Christensen and Wang, 1985;Winkler, 1985;Prasad and Manghnani, 1997).

MODELING AND DISCUSSION
Frequency-independent attenuation is observed for dry samples and samples confined at P c ≥ 14 MPa, and it could be caused by grain-to-grain friction (Walsh, 1966). Measurements show that as confining pressure increases, elastic moduli increase, and 1/Q decrease ( Table 2). Similar results for ultrasonic waves were presented, for instance, by Molyneux and Schmitt (1999), Prasad and Manghnani (1997), and Toksöz et al. (1979). Elastic moduli increase as confining pressure closes cracks increasing grain-to-grain traction (Zimmerman, 1991). A decreasing 1/Q might be explained by pressure-dependent crack closure and grainto-grain traction increment (Molyneux and Schmitt, 1999). The contact between two grains occurs only on a limited area throughout the contact of several asperities populating the grains' surfaces (Archard, 1953). Hulikal et al. (2015) suggest that each single asperity contact rheology could be modeled by a standard linear solid exhibiting a specific relaxation time (i.e., corner frequency). Relaxation time is a function of the asperity state and the petrophysical properties of the material. A large number of asperity contacts, making up the whole population of grain-to-grain contacts of a rock sample, will eventually exhibit a broad spectrum of corner frequencies. The superposition of all the corner frequencies would eventually produce a macroscopic rheology that might be described by a nearly constant attenuation model (Liu et al., 1976). We hypothesize that a decrease in attenuation, as P c increases, could be the result of contact stiffening and increase of the grain contacts areas (Archard, 1957). In fact, contact stiffening would change the strain partitioning between contacts and solid grains, hindering the overall viscous behavior due to contact deformation . Attenuation measurements are frequency-dependent for partially saturated Berea sandstone and confining pressure < 14 MPa. Measurements agree with the sum of the nearly constant Q model (NCM) and the standard solid model (SLS), confirming that attenuation and Young's modulus obey the Kramers-Kronig relation (Lakes, 2009 ; Figures 4-7). As already suggested by Kuteynikova et al. (2014) and Quintal (2013, 2014) in partially saturated Berea sandstone attenuation can result from the sum of frequency-independent and frequency-dependent mechanisms. Here, we assume that also in saturated samples, friction causes the frequency-independent attenuation, while patchy-saturation and/or squirt-flow cause the frequency-dependent attenuation (White, 1975;Dutta and Odé, 1979;Mavko and Jizba, 1991) (White, 1975;Dutta and Odé, 1979;Mavko and Jizba, 1991). In the next paragraphs, we will test such a hypothesis by employing some theoretical models to fit the observations. The petrophysical parameters employed in the modeling are listed in Table 3. FIGURE 10 | Longitudinal (V p ) and transverse (V s ) wave speeds measured for samples BS5s and BS3 as a function of vertical stress (σ 1 ). Ultrasonic wave speeds were measured for all the cases except case ii. However, for the sake of visualization, we show only a few cases. Error bars are calculated based on manual-picking uncertainty, sample length uncertainty and apparatus calibration.

Evidences for Patchy-Saturation
Patchy-saturation characteristic frequency f patchy sets the limit between the flow or "relaxed" limit and no-flow or "unrelaxed" limit. The two limits are valid for f << f patchy and f >>f patchy , respectively. As already mentioned in par. 2.5, in our case patchy-saturation might operate at frequencies between 28 Hz and 63 kHz. Thus, ultrasonic measurements (E HF ), which have been collected at f > 100 kHz, represent a good proxy for the "unrelaxed" limit.
The unrelaxed Young's modulus (E U ) can be calculated as: where K U is the unrelaxed bulk modulus calculated according to the Biot-Gasmmann-Hill limit (Johnson, 2001): where K 1 or K 2 (i.e., K 1 ∨ 2 ) are: where K f is the liquid bulk modulus, K a is the air bulk modulus, K dry is the dry bulk modulus measured at ultrasonic frequencies, i.e., K HF (S a = 1) . At low frequencies, we observe the left flanks of attenuation peaks whose corner frequencies are typically above 100 Hz. Thus, E SF at f < 100 Hz represents a good proxy for the "relaxed" limit. We calculate Young's modulus for the "relaxed" limit (E LF ) by averaging E SF measurements at frequencies < 10 Hz.
The relaxed Young's modulus (E R ) is then calculated as: (Mavko et al., 2009), where K R is the relaxed bulk modulus for a specific fluid saturation (i.e., K R (S a < 1) ) and can be estimated by means of the Gassmann fluid substitution: (Gassmann, 1951). Because we have no direct measurements of bulk modulus at low frequencies and dispersion is negligible at seismic frequency, in Eq. 17 we use K dry . Dry Young's modulus calculated from K HF and µ HF agree very well with low frequency dry Young's modulus indicating negligible dispersion. In Eq. 17 K f is the effective fluid bulk modulus that can be estimated as: where K a and K l are the bulk moduli of the gas and the liquid, respectively (Reuss, 1929). The liquid can be either water or GWmix. According to Gassmann (1951) shear modulus is not affected by saturation, therefore we assume that the relaxed shear modulus µ R is: E HF and E LF at different saturations and confining pressures (i.e., 2, 6, 10, and 14 MPa) are then compared to E U and E R (Figure 11). E LF and E HF follow the two theoretical limits and the dispersion -i.e., the vertical distance between E LF and E HF decreases with increasing confining pressure, similarly, to the dispersion between E R and E U , suggesting again that the wave-induced-fluid-flow mechanism is more active at low confining pressures. Both theoretical limits do not predict the decrease in Young's modulus between dry and saturated conditions. Such a behavior is common in sandstones and, according to the literature, in rocks like Berea sandstone, it could be related to the softening of the cement or the swelling of the clay minerals at grain boundaries (Cadoret, 1993).

Evidences for Squirt Flow
Squirt-flow could also explain the observed frequency-dependent attenuation curves (Mavko and Jizba, 1991;Gurevich et al., 2010). Here we show that the observed corner frequencies are compatible with squirt-flow theory and typical values of crack aspect ratios for Berea sandstone. We estimate the range for crack aspect ratio by means of Eq. 9 and by substituting f squirt with f c . This, provides a crack aspect ratio α ranging between ∼1.4x10 −4 and ∼4.5x10 −4 , which is in agreement with the literature (Zimmerman, 1991).
On the other hand, ultrasonic measurements suggest that crack closure occurs at confining pressure ∼40 MPa, where ultrasonic velocities vs. confining pressure curves tend to flatten (Figure 10). Walsh (1965) indicates that elliptical cross-section cracks under plane-stress conditions close at confining pressure P cc : by substituting P cc with 40 MPa and E g with 95 GPa (i.e., quartz Young's modulus), we obtain a value of crack aspect ratio α of ∼8x10 −4 . The fact that α and α , which are independently estimated, are similar suggests that our samples bear cracks whose geometry could have generated squirt-flow in a bandwidth close to 100 Hz.

Squirt-Flow vs. Patchy-Saturation Modeling
Squirt flow and Patchy saturation have been calculated for the conditions of case ii, iii, and v and confining pressures of 2, 8 and 12 MPa, i.e., we generated nine attenuation and dispersion curves (Figure 12). Precisely, Figures 12A,B show the comparison between the modeling results and the best fitting curves from the phenomenological models.
Patchy saturation was calculated for an interval of patchy radii comprised between 3 and 30 mm. Squirt flow was calculated for an interval of P h comprised between 20 and 100 MPa. Such values were chosen to span an order of magnitude and according to the sample dimensions and the increase of ultrasonic wave speed as a function of confining pressure. As a consequence, in Figure 12, the modeling results appear as shaded areas where the two attenuation mechanisms operate.
From the comparison, we observe that patchy saturation accurately describes attenuation in 3 cases - Figure 12, curves 4, 7 and 9 -and slightly underestimates the measurements in one case - Figure 12, curve 8. However, in 5 cases, patchy saturation strongly underestimates the measurements presenting corner frequencies much higher than those observed in the laboratory. Similarly, squirt flow accurately describes attenuation in 3 cases - Figure 12, curves 7, 8, and 9. However, in 6 cases, squirt flow strongly underestimates the measurements presenting corner frequencies much higher than those observed in the laboratory. Although the modeled Young's moduli overestimated the measurements of 1 to 4 GPa, the modeled increase of Young's modulus with confining pressure is similar to the experimental results FIGURE 11 | Young's modulus measured at low frequencies (1-10 Hz, E LF ) and at ultrasonic frequencies (0.1 MHz, E HF ) compared to the Biot-Gassmann-Hill limit (E U ) and the Gassmann fluid substitution (E R ) (Gassmann, 1951;Johnson, 2001). The former represents a proxy for patchy-saturation model (Mavko et al., 2009). (Figure 12), suggesting a nearly constant offset between the measurements and the models.
From the present analysis, it is difficult to assess whether patchy saturation or squirt flow is more appropriate to describe the observed frequency-dependent attenuation. In particular, the largest discrepancies between modeling and experimental data occur at 8 MPa confining pressure suggesting that both the models do not accurately capture the evolution of the petrophysical parameters with confining pressure.
The models are less accurate in describing the watersaturated case than the GWmix saturated cases. In particular, they largely overestimate the corner frequencies of case ii but less those of case iii and v. For instance, although case ii and iii have similar saturation and probably patchy sizes, the models describe better case iii. This might suggest that the characteristic frequency of patchy-saturation and squirtflow might not be a simple inverse relation with viscosity. The fluid flow of low-viscosity fluids in porous media might be turbulent, increasing the dragging forces and limiting the fluid flow (Reynolds, 1883). Such a process would result in reduced mobility (Batzle et al., 2006). Thus, we speculate that under certain conditions, wave-induced-fluid-flow might not be laminar as assumed in squirt-flow and patchy-saturation models manifesting as a nonlinear relation between corner frequencies and petrophysical parameters.

Final Remarks and Future Work
Models for both patchy-saturation and squirt-flow predict that attenuation fades as confining pressure increases. Such behavior is also present in our experimental data. Nevertheless, during experiments, the sample was deformed of a fixed strain amount. Thus, as confining pressure increases and cracks close, the strain must redistribute differently within the sample. At low confining pressures, a change in vertical stress mainly causes the splitting of vertical cracks and the closure of horizontal cracks. In contrast, at higher confining pressure, cracks tend to be closed, and a change of vertical stress causes shear strain along inclined faults (Jaeger et al., 2007). In the former case, crack deformation creates the local overpressures needed to generate the wave-induced-fluid-flow at the microscopic scale (i.e., squirt-flow). Nevertheless, we assume that overpressures generated in cracks also control the mesoscopic overpressure leading to patchy-saturation. The disappearance of frequencydependent attenuation for P c ≥ 14 MPa could be related to the closure of the cracks and the strain distribution variation affecting both the attenuation mechanisms analyzed.
The modeling that we present here should be furthered in future work. Petrophysical parameters could be estimated more precisely following the literature, and different models could be employed. For example, we could follow de Paula et al. (2012) to better estimate squirt-flow parameters, Buckingham FIGURE 12 | Comparison between (i) measured attenuation (A -symbols) and Young's modulus (B -symbols), and (ii) calculated attenuation and Young's modulus for squirt-flow (red) and patchy-saturation (green) models. Squirt-flow is modeled for a minimum and a maximum value of P h , 20 and 100 MPa, respectively. Patchy saturation is modeled for a minimum and maximum value of r p of 0.003 and 0.03 m, respectively.
(1999) to verify the hypothesis that frequency-independent attenuation in dry samples is related to grain friction, and use Vogelaar et al. (2010) to model the patchy saturation. Since the observed attenuation might result from the combination of squirt flow and patchy saturation, we could model our data using the double double-porosity model proposed by Ba et al. (2017). Besides, we could test analog materials whose physical properties could be thoroughly controlled (e.g., Xie et al., 2019). Nevertheless, results and modeling presented here suggest that in near-surface sedimentary terrains, where the effective stress is low -e.g., < 20 MPa -and partial saturation is likely to be present, attenuation and dispersion might strongly affect the propagation of seismic waves like suggested by field observations (e.g., Abercrombie, 1997).

CONCLUSION
We have measured low-frequency seismic wave attenuation and Young's modulus, and ultrasonic wave speeds for three dry and partially saturated Berea sandstone samples. Measurements have been conducted at confining pressure between 0 and 40 MPa, and frequency ranges 1-100 Hz, and 0.1-1 MHz. Seismic wave attenuation in samples partially saturated and confined with pressures below 14 MPa are frequency-dependent and as high as 0.08. On the other hand, the dry sample and samples confined above 14 MPa exhibit negligible and frequency independent attenuation.
Low-frequency attenuation and dispersion curves obey to the Kramers-Kronig relation. Low frequency and ultrasonic Young's moduli are compatible with effective-medium theories that are calculated for relaxed and unrelaxed conditions. This suggests that the observed frequency-dependent attenuation is caused by a wave-induced-fluid-flow type of mechanism and that attenuation could be significant in the near-surface (<∼1 km) where effective stress is low and partial saturation might be present. Nevertheless, it is not yet clear whether the observed frequency-dependent attenuation is the result of squirt-flow or patchy-saturation as analytical modeling provides similar results for these two attenuation mechanisms. Such and other aspects, like the role of non-laminar fluid-flow in controlling attenuation, are not yet clear and demand further investigation.

DATA AVAILABILITY STATEMENT
The full datasets including data not shown in figures can be found here: https://zenodo.org/record/4517286#.YCCtDXlMGUk.

AUTHOR CONTRIBUTIONS
NT conceived and performed the experiments, analyzed and modeled the results, wrote the manuscript, and prepared the figures and tables. CM assisted in preparing and performing the experiments. ES provided funding and supervised the study. All authors contributed to the article, edited the manuscript, and approved the submitted version.

FUNDING
This project has been partially funded by the Swiss Innovation Promotion Agency (KTI) project nr. 22701.