CO2‐Brine Substitution Effects on Ultrasonic Wave Propagation Through Sandstone With Oblique Fractures

Seismic monitoring of injected CO2 plumes in fractured storage reservoirs relies on accurate knowledge of the physical mechanisms governing elastic wave propagation, as described by appropriate, validated rock physics models. We measured laboratory ultrasonic velocity and attenuation of P and S waves, and electrical resistivity, of a synthetic fractured sandstone with obliquely aligned (penny‐shaped) fractures, undergoing a brine‐CO2 flow‐through test at simulated reservoir pressure and temperature. Our results show systematic differences in the dependence of velocity and attenuation on fluid saturation between imbibition and drainage episodes, which we attribute to uniform and patchy fluid distributions, respectively, and the relative permeability of CO2 and brine in the rock. This behavior is consistent with predictions from a multifluid rock physics model, facilitating the identification of the dispersive mechanisms associated with wave‐induced fluid flow in fractured systems at seismic scales.


Introduction
From the microscopic scale to large faulting systems, cracks are present in almost every rock formation found in the Earth's crust. Anisotropic and potentially frequency-dependent wave propagation are key geophysical signatures describing cracked rocks. Such properties are crucial to assess the suitability of hydrocarbon reservoirs and saline aquifers for geological carbon storage (GCS) (Chiaramonte et al., 2015;Iding & Ringrose, 2010).
Previous studies have analyzed a number of factors affecting the brine-CO 2 saturation dependence in saline sandstone reservoirs-the most suitable geological context for GCS (Michael et al., 2010). These factors include the distinction between pore pressure and pore fluid distribution effects (Falcon-Suarez et al., 2016, the frequency dependence of elastic wave properties and the methodology to upscale information collected in laboratory (Lei & Xue, 2009;Mikhaltsevitch et al., 2014;Nakagawa et al., 2013), or the effect of mineralogical changes in the elastic and transport properties of the rock (Canal et al., 2013;Hangx et al., 2010Hangx et al., , 2015Vialle et al., 2014;Vialle & Vanorio, 2011). However, most of these studies are restricted to nonfractured rocks (Nooraiepour et al., 2018), although the nucleation and reactivation of fractures endanger safe GCS (Rutqvist, 2012;Velcin et al., 2020).
The distribution of different pore fluids is largely conditioned by reservoir heterogeneities, wettability (Al-Khdheeawi et al., 2017), and saturation history (Knight & Nolen-Hoeksema, 1990). To address the uncertainties associated with CO 2 -brine partially saturated fractured rocks, we need to conduct controlled experiments on samples with well-defined physical and structural properties (e.g., Amalokwu et al., 2016;Rathore et al., 1995;Tillotson et al., 2011), to generate data sets to help constrain robust rock physics models. However, fluid-dependent phenomena dominating in laboratory observations such as frequency-dependent behavior or fluid distribution in the pore space are rarely examined in conjunction with their intrinsic anisotropy (Amalokwu et al., 2014;Murphy, 1984).
Few rock physics theories incorporate the subtleties associated with such partially saturated cracked rocks holistically. Recently, Papageorgiou and Chapman (2017) have suggested a theory that combines uniform/patchy fluid distribution (Dutta & Odé, 1979;White, 1975) and squirt flow (Chapman, 2003;Dvorkin et al., 1995) phenomena for partially saturated isotropic rocks. Later, Jin et al. (2018) extended their work to incorporate mesoscopic flow that occurs between fractures and the pore space for anisotropic rocks. The model of Jin et al. (2018) is therefore appropriate to describe frequency-dependent wave velocities of partially saturated, cracked samples. If properly calibrated, this model also could be used to describe seismic velocities obtained during and after CO 2 injection in fractured GCS reservoirs.
Here, we present a brine-CO 2 partial saturation experiment in cracked sandstone with a well-defined fracture network (cracks aligned at oblique angle), under realistic geological conditions of confining (40 MPa) and pore (10 MPa) pressure. During the test, we measured ultrasonic P and two orthogonal S wave velocities and attenuations, together with electrical resistivity that we used to determine the evolution of the degree of CO 2 saturation (up to~60%, including drainage and imbibition saturation paths). We simultaneously fit ultrasonic velocities to constrain the rock physics model of Jin et al. (2018), which explains the observed anisotropic and dispersive properties, and variation with partial saturation, as well as a counterintuitive increase of the velocity despite the decrease of stiff fluid. The transport and mechanical results are assessed in Muñoz-Ibáñez et al. (2019).

Experimental Procedure
We used a~2-cm length,~5-cm diameter synthetic sandstone core plug, containing fractures aligned at 45°f rom its axis. The sample was manufactured using a mixture of sand, kaolinite, and sodium silica gel , but with a predetermined number of 2-mm diameter (δ), 0.2-mm thickness (τ) aluminum discs (Amalokwu et al., 2015;Tillotson et al., 2012). Once baked, the sample was immersed in an acidic bath to remove the aluminum discs, resulting in a silica-cemented sandstone with 45°aligned penny-shaped cracks. The fracture density (ε f ¼ 0.0298 ± 0.0077) and an average fracture aspect ratio (τ/δ ¼ 0.088 ± 0.001) were obtained from X-ray CT scan image processing (Amalokwu et al., 2015), which leads to a porosity fraction attributed to fractures of~10% (Muñoz-Ibáñez et al., 2019).
The sample was subjected to steady-state brine-CO 2 flow-through (BCFT), using the experimental rig for CO 2 storage multiflow tests described in Falcon-Suarez et al. (2017). The rig implements sensors for measuring, simultaneously, electrical resistivity and ultrasonic waves (velocity and attenuation) using the pulse echo technique (Amalokwu, 2016;Best, 1992;McCann & Sothcott, 1992). Here, we used a sensor that incorporates P and S wave transducers in one platen, and a single S wave transducer in the other (all transmitting 400-to 1,000-kHz broadband acoustic pulses). For these transducers, the velocity precision is ±0.1% with an accuracy of ±0.3% (95% confidence), while attenuation accuracy is ±0.2 dB/cm (Best, 1992). The S wave transducer polarizations were aligned parallel and perpendicular to the fracture strike, measuring oblique S 1 (fast) and S 2 (slow) wave velocity and attenuation, respectively, for estimations of the shear wave anisotropy (expressed as (n S1 − n S2 )/n S1 , for a given S wave attribute n).
The BCFT test was conducted under constant hydrostatic confining (P c ¼ 40 MPa) and pore (P p ¼ 10 MPa) pressure, at room temperature (20°C), using 3.5% NaCl synthetic brine and liquid CO 2 , and covered drainage and imbibition episodes. An overpressure limit of 2 MPa was imposed (Muñoz-Ibáñez et al., 2019), which prevents rock mechanical damage (fracture propagation) at the simulated conditions (Velcin et al., 2020). The total flow (Q ¼ Q w + Q CO2 ) was initially regulated in each episode (Q < 0.5 cm 3 min −1 ), and the CO 2 fractional flow (X CO2 ¼ Q CO2 /Q) increased in 0.2 units episodes during the drainage part of the test. After the last drainage episode, the BCFT test concluded with a 100% brine flow-through to force the natural imbibition process that occurs in saline aquifers after ceasing the CO 2 injection (Gaus, 2010).
During the test, axial and radial strains were measured continuously (every second), while ultrasonic P and S 1,2 wave velocities (V P , V S1,2 ) and attenuations (inverse quality factors Q P −1 , Q S1,2 −1 ), and resistivity were acquired, at least, every 0.5 times the sample pore volume (PV~10.46 cm 3 ) of flow-through. Then, using Archie's (1942) first and second laws, the bulk resistivity of porous media partially saturated with brine and CO 2 was transformed into degree of brine saturation (Carrigan et al., 2013;Falcon-Suarez et al., 2016Muñoz-Ibáñez et al., 2019;Nakatsuka et al., 2010). For resistivity, the error increases from <1% up to 5% with resistivity and sample anisotropy (North et al., 2013). The upper value leads to a saturation uncertainty of ±5%, when converting resistivity into degree of saturation using Archie's approach with a saturation exponent n ¼ 2, usually adopted for sandstones (Mavko et al., 2009). However, the dual porosity of our sample influences the rock wettability and the saturation history (Muñoz-Ibáñez et al., 2019). In turn, n might become variable for a given resistivity-saturation relationship, invalidating the use of Archie's law (Suman & Knight, 1997;Zhou & Stenby, 1997). To account for this, we propagate saturation error derived from ±10% variations in n, which covers the deviation from drainage to imbibition in nonuniform pore size distribution in water-wet samples (Suman & Knight, 1997).
After the BCFT test, the sample was cleaned with deionized water and then dried, before being saturated with CO 2 at the test conditions to collect the corresponding ultrasonic properties. Further information about the testing methodology in Falcon-Suarez et al. (2017)

Ultrasonic Waves
The pulse-echo method used in this study (Winkler & Plona, 1982) provides accurate measurements of ultrasonic velocities and attenuation (Amalokwu, 2016;Best, 1992;Best et al., 2007;Tillotson et al., 2011Tillotson et al., , 2012. Figure 1a shows energy dissipation related to multiple arrivals for both the P and S wave signals between top and base reflection peaks (framed within the data-processing windows W1 and W2, respectively). The energy dissipation is independent of the partial saturation for S 1 and S 2 , which trends are almost identical under different partial CO 2 saturations (S CO2 ), while multiple internal reflections distort the signal after the first pulse. P wave also evidences scattering arrivals in the time between the top and base reflections, but varies with the saturation state, which magnifies the overall energy dissipation, as calculated from the reduction in amplitude of the base reflection relative to the top reflection.
The horizontal and vertical cross-sections of the fractures (Figure 1b) show local fracture propagation that in some cases can exceed original fracture length (e.g., Y3). This fracture elongation is preferentially inducing energy dissipation for S waves, which the wavelength (3 mm) is closer to the fracture length (2 mm) than that of P waves (5 mm). Therefore, the attenuation data shown in Figure 2 have to be understood in qualitative terms, with Q P −1 showing significant changes between drainage and imbibition paths, while the Q S1 −1 and Q S2 −1 are largely affected by scattering with no significant changes (i.e., variations lie within the experimental error).

BCFT Test
The brine-CO 2 flow-through (BCFT) test lasted~75 hr, resulting in a total fluid throughput of~63 PV. Figure 2 shows the ultrasonic properties (V P , V S1,2 and Q P −1 , Q S1,2 −1 ; all the results at 600 kHz obtained from Fourier analysis of the broadband signals), plotted together with the electrical resistivity and volumetric strains measured during the BCFT test.
The transition from the first brine: CO 2 fractional flow episode (X CO2 ¼ 0) to the second (X CO2 ¼ 0.2) marks the arrival of the free CO 2 phase in the rock (PV~3.6), and the resistivity increases progressively with the volume of CO 2 passing through the rock. The CO 2 arrival leads to a sharp increase of V P (~1.5%) and Q p −1 (~50%), followed by a progressive decrease of V P with the increasing CO 2 content and the opposite for Q P −1 . The two orthogonal components of the shear wave show similar trends in both the velocities (V S1 , V S2 ) and the attenuations (Q S1 −1 , Q S2 −1 ), with relative differences of~6% between V S1 and V S2 and throughout the BCFT test. V S1 and V S2 drop during X CO2 ¼ 0.2, followed by a soft increasing trend thereafter. Similarly, Q S1 −1 and Q S2 −1 progressively increase during drainage and drop to approximately the original values during the forced imbibition episode (X CO2 ¼ R-0). However, this trend is inconclusive as it occurs within the experimental uncertainty.
The imbibition (R-0) started with a single brine pulse (~0.5 PV at 42.8 PV; BP in Figure 2), with continuous monitoring of the geophysical properties of the sample. Then, continuous brine flow through the sample was progressively replacing the CO 2 from the pore space, as inferred from the resistivity drop (down to~24% 10.1029/2020GL088439

Geophysical Research Letters
FALCON-SUAREZ ET AL.
above the original brine-saturated conditions). The imbibition affected P and S wave properties differently: V P increased and Q P −1 decreased, both up to~2% below and above the original brine-saturated values, respectively, while V S1,2 and Q S1,2 −1 fully recovered by the end of the imbibition. Yin et al. (1992) also found that the attenuation depends not only on the degree of saturation and the frequency of the measurements but also on the saturation history. They reported~7% shift on the degree of saturation for the attenuation peak from drainage to imbibition paths at sonic frequencies (~1,500 Hz). In our case, this shift reaches up to 30% for Q P −1 (Figure 2), because of the interplay of the scattering on heterogeneities (crack-induced internal reflections) and energy dispersive mechanisms related to partial fluid flow in porous media. This render the attenuation measurements unreliable as quantitative indicators of the changing physical properties of the sample due to pore fluid substitution.
The sample evidences mechanical deformation (inflation) after X CO2 > 0.5, which stabilizes during imbibition. This minor inflation (~0.12%) is within the porosity fraction commonly attributed to microcracks (0.7%; Fortin et al., 2007). Thus, it could be related to microcracks reopening due to differential pressure (P diff ¼ P c − P p ) drop during the experiment associated with the increase of inlet-to-outlet P p gradients with the increasing X CO2 (Muñoz-Ibáñez et al., 2019).

Insights From Data Analysis
Theories for wave propagation through partially saturated rocks emphasize two dominant dispersive mechanisms taking place when wave propagation induces local and mesoscale fluid flow: squirt flow (Chapman, 2003;Dvorkin et al., 1995) and patchy saturation (Dutta & Odé, 1979;White, 1975) effects. While the diffusive character of Biot flow (Biot, 1956) is significant in unconsolidated sediments, in consolidated rocks, the predicted dispersion is much smaller and occurs at much higher frequencies than the experimental frequency. We therefore focus on the first two: squirt flow and patchy saturation.
Papageorgiou and Chapman (2017) combined these two effects into a single theory. Later, Jin et al. (2018) extended their theory to include anisotropic effects due to fractures, showing that the moduli of a fractured, partially saturated rock can be interpreted as two coupled standard linear solid (SLS) models (Mavko et al., 2009). These two SLS models have relaxation frequencies linked by the scale ratio of the fracture size to the pore/microcrack size. For a rock with aligned fractures, the SLS with lower relaxation frequency Figure 1. P and S 1,2 wave signals (a) at the saturation states corresponding to fractional flow episodes of pure brine (X CO2 ¼ 0) and CO 2 (X CO2 ¼ 1). W1 and W2 indicate the time windows used for the Fourier analysis (Best, 1992). (b) Thin section analysis from two (mutually orthogonal) slides prepared after the BCFT test. From original images, color channel filtering was applied to emphasize porous (bright) with respect to grains (dark) areas.
corresponds to an anisotropic mechanism that induces dispersion in the fracture normal direction only, whereas the SLS with higher relaxation frequency is isotropic.
In partial saturation, both the combined squirt-patch effect and the effective fluid moduli impact the two SLS relaxation frequencies (Brie et al., 1995;Papageorgiou & Chapman, 2017). The relaxation frequencies depend on relative permeability and a patch parameter q, which lies between 1 at uniform saturation and the ratio of the fluid moduli at patchy saturation (Papageorgiou & Chapman, 2017). The magnitude of each SLS model induced dispersion is controlled by the fracture density (ε f ) and by microcrack density (ε m ), for the SLS with lower relaxation frequency and higher relaxation frequency, respectively.
To calculate the relaxation frequencies at partial saturation, we use the Brooks-Corey relative permeability model, suitable for explaining the two fluid phases' interaction in our sample (Muñoz-Ibáñez et al., 2019), and the fracture and microcrack densities ε f ¼ 0.042 and ε m ¼ 0.02. ε f is derived from the design of synthetic fractures in the sample (Tillotson et al., 2012), while ε m is arbitrarily adopted within the range of values observed in similar synthetic rocks . In the Supporting Information, we show how we derive the V P , V S1 , V S2 of an effective medium with these properties. The increase in V P at low saturation, together with the absence of fluid-dependent shear wave anisotropy, is consistent with the behavior of our model in the frequency regime, where the large fractures are effectively sealed but squirt flow dispersion is still observable isotropically. Current data prevent us from further constraining these frequencies, so we present modeling for a range of relaxation frequencies of the microcrack SLS relative to the experiment frequency, as well as saturation distribution scenarios for patchy and uniform saturated pores.  (Figures 3b-3d). For the fit, we perform a nonlinear Nelder-Mead minimization in two steps for each of the three relaxation frequencies ω i . First, we deduce effective medium matrix parameters λ, μ at full water saturation; that is, we determine λ, μ given constant ε, ε f , ω i such that, when S w ¼ 1, (V P(th) (λ, μ) − V P(ex) ) 2 +(V S1(th) (λ, μ) − V S1(ex) ) 2 +(V S2(th) (λ, μ) − V S2(ex) ) 2 is minimum (subscripts "th" and "ex" for theoretic and experimental, respectively). As a second step, we invert for the patch parameter at partial saturation for each of the imbibition/drainage stages; that is, given λ, μ, ε, ε f , ω i , we determine where S j are the saturation steps for which we have data.
Despite the difference in fluids and incomplete saturation path, our V P results are similar in form to the observations in Knight and Nolen-Hoeksema (1990), where velocities showed strong hysteresis with values consistently higher during drainage. Knight and Nolen-Hoeksema (1990) attribute this discrepancy to a difference in pore-fluid distribution under drainage (patchy saturation) and imbibition (uniform saturation), which is consistent with our observations. However, we note the lack of significant V S hysteresis and the strong impact of the saturation distribution to the velocity dispersion in our work.
The absence of imbibition data with S w ≥ 0.85 prevents saturation history analysis at the high saturation regime and indicates a residual CO 2 trapping of~15% at the experimental conditions. However, the Figure 2. Brine-CO 2 flow-through test in fracture sandstone. P and S 1,2 wave velocities (V P , V S1,2 ), attenuations (Q P −1 , Q S1,2 −1 ) and electrical resistivity, together with volumetric strains for six brine: CO 2 flow rates drainage episodes (X CO2 ¼ 0 to 1) and the forced imbibition (X CO2 ¼ R-0; BP denotes the initial imbibition brine pulse; see text). The ultrasonic properties were measured at a single frequency of 600 kHz (pulse-echo technique), obtained from the Fourier analysis of broad band signals. Error bars displayed every 10 measurements.

Geophysical Research Letters
FALCON-SUAREZ ET AL.
increase observed in porosity (~5%) and permeability (~45%) after the BCFT test might have occurred during drainage. Then, abnormally high values of V S1,2 measured at full CO 2 saturation (marked by stars in Figure 3), together with the elevated V S1,2 and reduced V P values with respect to modeling along the imbibition path, indicate the reduction of the sample's density (i.e., porosity increase) affected its shear modulus minimally, in agreement with Delle Piane and Sarout (2016).
We also observe a striking increase in V P at low CO 2 saturation (S w~0 .9), despite the addition of softer fluid in the system. Given the similar density between the CO 2 and brine, this increase is attributed to an increase in the elastic modulus. This phenomenon is unusual but still identified in previously reported data for sandstones (Batzle et al., 2006), limestones (Cadoret et al., 1995), and shales (Szewczyk et al., 2018), with air and water as pore fluids under atmospheric pressure conditions. In our experiment, the viscosity and moduli of the two fluids differ by only one or two orders of magnitude within the experimental bandwidth, which magnifies the observed effect with respect to the previous observations. This increase is also present in the modeling predictions, and we have excluded CO 2 dissolution into water as a possible cause for this phenomenon as its impact would be of the order of <1%, much smaller than the observation. The reason for this increase lies in the complex interplay between partial fluid mobility and effective fluid modulus: a S w range exists between 0.8 and 0.9 in which the dispersion from frequency-dependent effects overcompensates for the addition of the soft fluid, making the system stiffer at partial than at full saturation. , V S1 and V S2 (c), and S wave velocity anisotropy, SWVA ¼ (V S1 − V S2 )/V S1 × 100 (d), versus degree of brine saturation (S w ). Two modeling scenarios for pore fluid distribution: uniformly (U, with patch parameter q ¼ 1) and patchy (P, with q ¼ 0.2), at frequencies ω 1 ¼ 0.027, ω 2 ¼ 0.17, and ω 3 ¼ 1.0725. Subscripts d and i refer to drainage and imbibition, respectively. The star symbol indicates the CO 2 saturation measurement (see text). Two single cross bars show horizontal and vertical errors for low and high S w .

Discussion
The ultimate goal of our work is to understand the geophysical behavior of fractured reservoirs undergoing CO 2 injection and identify the mechanisms that can be used to effectively monitor GCS projects. We extend the work of Nooraiepour et al. (2018) in considering properties of distributions of fractures, which is closer to the equivalent medium concept required to scale to reservoir conditions. To date, data regarding partial saturation in fractured media are limited to controlled experiments at atmospheric conditions using water and air as nonwetting phase (e.g., Amalokwu, 2016). The air has higher compressibility, lower viscosity, and density 3 orders of magnitude lower than ours at the experimental conditions (i.e., liquid CO 2 ). Therefore, we expect to find discrepancies between our results and those found in the literature related to the nonwetting fluid properties.
The properties of the pore fluids and their mobility through the porous media complicate the distinction between pore fluid distribution and rock fabric in fractured media. Kong et al. (2017) remarked that the fracture-filling fluid can disguise P wave dispersion and attenuation effects due to fracture orientation, when its modulus is much smaller than that of the pore-filling fluid. Amalokwu et al. (2015) found that fractures obliquely oriented with respect to wave propagation can lead to S wave velocity anisotropy (SWVA) drop with increasing water saturation, because V S2 (sensitive to fluid compressibility reduction) decreases, while V S1 , considered independent of the pore fluid (Tillotson et al., 2012), remains constant. Conversely, our results show minor SWVA variation with the CO 2 content, within the bounds imposed by experimental error.
The distribution of the two fluids in our porous medium is conditioned by the wettability of the CO 2brine-rock system, the pore geometry, and the saturation history (Cinar & Riaz, 2014;Krevor et al., 2015). Our experiment simulates different stages of CO 2 injection in fractured saline siliciclastic reservoirs, at realistic geological stress conditions. The test covers fluid substitution during CO 2 plume advance (drainage) and the natural aquifer recharge after ceasing the injection activities (imbibition) with a maximum CO 2 saturation~0.6, in agreement with previous CO 2 -water fluid substitution tests in sandstones (Burnside & Naylor, 2014). The low capillarity forces of the fractures with respect to the pores lead to preferential CO 2 fracture filling (Muñoz-Ibáñez et al., 2019), which magnifies the effect of drainage versus imbibition paths.
Our results indicate that P wave properties (V P , Q P −1 ) are highly sensitive to the saturation history, while S waves show very little differences between drainage and imbibition saturation paths. This finding has important implications for CO 2 storage, since the information inferred from seismic surveys during the injection and post injection activities are significantly different; but this also applies for the differential monitoring of the front (drainage) and the trailing edge (imbibition) of the CO 2 plume. The degree of partial saturation in sandstones with fractures aligned at oblique angles might have some effect in the SWVA, in agreement with Amalokwu et al. (2016), but certainly insufficient to be used as a robust indicator of pore fluid distribution. This finding has direct application on the interpretation of S wave splitting data sets from gas and GCS reservoirs containing or developing fractures during production or injection activities.
The use of ultrasonic waves to explore fractured media is influenced by the characteristics of the fracture, such as fracture length (c), aperture, crack surface roughness, and the saturating fluid, as well as the wavelength (λ w ). If c/λ w < 1, the porous medium appears homogeneous, while for c/λ w > 1, heterogeneities begin to act as energy scattering fronts. In our rock sample, c is slightly lower than λ w for both the P and S waves, which leads to some wave scattering in the recorded signals. In addition, regardless of how homogeneous the medium is, absolute scale also plays a part.
Here, we identify and analyze the mechanisms that are necessary to describe an effective medium that reproduces the dispersive behavior of the rock. The results from our controlled experiment suggest that fluid flow induces dispersive behavior at the pore and fracture scale (with a transition frequency that scales with effective mobility), which is sufficient to describe this fractured medium. Likewise, we observe that despite the single experimental frequency, dispersion varies across saturation consistently with our modeling strategy indicating that the dominant mechanisms-at least with respect to wave-induced fluid flow-have been taken into account in the effective medium. This interpretation is supported by Figure 3 saturation depending on drainage/imbibition paths, in conjunction with the relatively small dispersion in shear wave velocities.
Likewise, we anticipate that a fractured, partially saturated reservoir rock will exhibit dispersive behavior according to these mechanisms at the field scale. To make this fully predictive, we would need to calibrate the model with specific field velocity measurements and have a hydrological characterization of the reservoir formation. Our identification of these mechanisms offers a step toward overcoming the difficulty in distinguishing saturation effects from fluid distribution effects in fractured systems. This problem is dominant in the characterization and monitoring of CO 2 storage reservoirs where faults have been recognized (Chiaramonte et al., 2015;Iding & Ringrose, 2010).

Conclusions
In this work, we record experimental observation suggesting anisotropic rocks with aligned fractures have bulk and shear elastic behavior that depends on the saturation path (imbibition/drainage). During our high-pressure brine-CO 2 flow-through experiment on a synthetic sandstone, stiffening was observed for low (~10%) CO 2 saturation in the drainage phase, despite the apparent fluid softening. We found shear wave velocity anisotropy to be independent of saturation path and were unable to draw conclusions about its dependence on fluid content due to high uncertainties associated with the S wave measurements. Using recently published work, we were able to explain the dependence of the elastic constants on imbibition and drainage. We observe manifestations of different patchiness in the fluid distribution, with more uniform saturation during imbibition and patchy saturation during drainage, in accordance to previously published experiments performed on isotropic rocks. Using the same rock physics model, we were able to capture the increase in modulus observed during drainage as a result from a decrease in effective mobility due to relative permeability effects.