Pressure-varying CO 2 distribution a ﬀ ects the ultrasonic velocities of synthetic sandstones

We performed a novel experiment in which three synthetic sandstones – manufactured using a common method but having di ﬀ erent porosities – were saturated with brine and progressively ﬂ ooded with CO 2 under constant con ﬁ ning pressure. The ﬂ uid pressure was varied around the critical pressure of CO 2 and repeated measurements were made of resistivity, in order to assess the saturation, and elastic wave velocity during the ﬂ ood. The measured saturated bulk moduli were higher than those predicted by the Gassmann – Wood theory, but were consistent with behaviour described by a recently derived poroelastic model which combines “ patch ” and “ squirt ” e ﬀ ects. Measurements on two of the samples followed a patch-based model while those on the highest porosity sample showed evidence of squirt- ﬂ ow behaviour. Our analysis suggests that the appropriate ﬂ uid mixing law is pressure dependent, which is consistent with the notion that the e ﬀ ective patch size decreases as ﬂ uid pressure is increased. We derive simple empirical models for the patch dependence from ﬂ uid pressure which may be used in seismic modelling and interpretation exercises relevant to monitoring of CO 2 injection.


Introduction
Estimation of CO 2 saturation through seismic methods is a central component of seismic monitoring geological carbon storage projects (Chadwick et al., 2006). Our ability to perform this task reliably and effectively depends on an understanding of the fundamental physics of wave propagation through rocks saturated with multiple fluids. Unfortunately, this area is not yet completely understood.
In the literature of CO 2 monitoring there is considerable discussion of the merits of "patchy" versus "uniform" saturation (e.g. Eid et al., 2015). Even though difficult to measure in the field, patchy saturation can be assessed in the controlled environment of the laboratory. In this regard, Lebedev et al. (2009) showed direct experimental evidence of fluid patches and assessed their influence in measured seismic velocities and Caspari et al. (2011) inverted patch-size from sonic logs.
Laboratory rock physics measurements have the clear potential to shed light on these issues, but unfortunately interpretation of the measurements is not always straightforward. Most patchy saturation models associate the existence of patches to a characteristic patch size (for example Dutta and Odé, 1979;Toms-Stewart et al., 2009). This patch size influences the wave propagation causing the wave to disperse. Many authors (e.g. Endres and Knight, 1997;Mavko and Jizba, 1991;King et al., 2000) have noted that ultrasonic velocity measurements are influenced by a different dispersive mechanism related most probably, to squirt flow.
Recently, Papageorgiou and Chapman (2017) proposed a model in which the two mechanisms of squirt flow and patchy saturation, are intertwined (for a different approach, see also Mayr and Burkhardt, 2006). The effect of patches persists at the low frequency limit of the model but also affects the squirt flow characteristic frequency in the ultrasonic range. This is manifested through an "effective patch" parameter which constrains both the low and high frequency range of the theory.
A sequence of ultrasonic experiments performed on synthetic rocks by Falcon-Suarez et al. (2016 in the rock physics CO 2 lab of the National Oceanography Centre in Southampton, have demonstrated that there may be a complicated relation between the pressure dependence of the CO 2 phase, the effect of pore pressure in the samples and the variation in ultrasonic velocities when the samples are partially saturated. In these experiments, synthetic samples are partially saturated with CO 2 of varying pressure at conditions close to the critical point, providing an ideal setting to test the exact interdependence of the dispersive mechanisms dominant in patchy saturation and ultrasonic squirt flow. In this work we use ultrasonic velocity data from two such experiments (Falcon-Suarez et al., 2016, one performed with a synthetic sandstone of 26% and one of 38% porosity. In addition we present new velocity data of a synthetic sandstone of high porosity (45%), obtained in a similar CO 2 -brine flow-through experiment.
We show that all three datasets can be explained by a single theoretical model. Two of the samples are adequately modelled under a low frequency assumption but the highest porosity sample shows behaviour consistent with an intermediate frequency. These results lead to a new fluid mixing law which, in contrast to previous work, is dependent on pressure rather than patch-size.
We begin by introducing the experimental setup used in these measurements and review the theoretical background for the model we use. We then analyse the experimental results and discuss the implications of our work.

Experimental setup and results
The experimental results in this work consist of CO 2 injection tests on three different synthetic sandstone samples: one of ≃26% porosity (sample 1), one of ≃38% porosity (sample 2) and a ≃45% porosity sandstone (sample 3). The data obtained from sample 3 are original and shown here for the first time (Table 1) whereas the data for samples 1, 2 have been presented in Falcon-Suarez et al. (2016 where the experimental setup is also described in detail.
Samples were manufactured by mixing quartz-sand, kaolinite and sodium-silica cement (> 90% silica) using the methodology proposed by Tillotson et al. (2012). After removal from the manufacturing mould, the resulting specimens were flushed with deionized water, and then cut and ground flat and parallel to within ± 0.01 mm. The physical properties of the samples are discussed in Falcon-Suarez et al. (2016.
The experiments were conducted on the CCS experimental rig at the National Oceanography Centre, Southampton (Falcon-Suarez et al., 2016). The rig has been designed to simultaneously measure geophysical and hydromechanical properties of rock samples exposed to the co-injection of up to two fluid phases, at controlled pressure and temperature conditions matching those of realistic shallow North Sea-like CO 2 storage reservoirs. The rig allows independent application of axial (σ 1 ) and lateral (σ 2 = σ 3 ) confining stresses up to 65 MPa, using a high pressure, high accuracy EX-100D (ISCO) dual-controller. Similarly, pore pressure (Pp) is also controlled by a second EX-100D. The sample is inserted into a triaxial cell core holder, with an inner sleeve equipped with 16 electrodes for bulk electrical resistivity tomography (ERT) measurements (North et al., 2013). Under our test conditions, the resistivity measurement error is 5% (at frequencies 1-500 Hz). The ultrasonic data are obtained using the pulse-echo technique of McCann and Sothcott (1992), providing P-and S-wave velocities (V p and V s ) with accuracies of ± 0.3% (Best, 1992). Further information about the features and multi-flow configurations of the experimental rig are presented in Falcon-Suarez et al. (2016.
The experimental procedure followed in the three tests consisted of a steady state flooding configuration in which brine-CO 2 fractional flows were increased by 20% stepwise from 100% (35 g/L NaCl) synthetic brine to 100% CO 2 . For each fractional flow, the pore pressure (Pp) was increased by 1 MPa stepwise from 7 to 12 MPa and back to 7 MPa, while keeping constant the confining pressure P c at 16.4 MPa under hydrostatic conditions (i.e. σ 1 = σ 2 = σ 3 = P c ). The maximum differential pressure (P d = P c − P p ) applied was P d ∼ 9.4 MPa. These values were selected in order to replicate the conditions of Utsira sand formation in the Sleipner field, North Sea (Chadwick et al., 2010;Eiken et al., 2011). The fluid temperature was held constant at 32°C for sample 2 and 35°C for samples 1, 3. Further details about fluid mixture and flow conditions can be found in Falcon-Suarez et al. (2016. At the end of each differential pressure/fractional flow step, ultrasonic velocity and electrical resistivity were measured in short Table 1 V p and V s measurements for sample 3. The dry measurements have been calibrated to the same state of stress as the corresponding wet measurements. At the time of publication, the dataset for sample 3 was in the process of becoming publicly available via the UKCCS data and information archive. succession. Resistivity was transformed into degree of saturation using Archie's law (Archie et al., 1942), particularly modified for CO 2 -brine porous media (Carrigan et al., 2013;Falcon-Suarez et al., 2016Nakatsuka et al., 2010).

Theory
In rock physics modelling the most established tool for performing fluid substitution in rocks is traditionally the model of Gassmann (1951). In Gassmann's model the saturated matrix bulk modulus K sat. is given in terms of the dry matrix, grain and fluid moduli K d , K m , K f and the porosity ϕ respectively, by the formula: Often, when the rock is partially saturated, this model is used in conjunction with Wood's law (see Domenico, 1974) in what is known as Gassmann-Domenico or Gassmann-Wood model. Strictly speaking, this model is valid only when the fluids are in pressure equilibrium and connected across the pore space. Under these assumptions, and taking the matrix to be saturated by CO 2 and brine, one simply replaces the fluid stiffness 1/K f by a volume averaged compliance = + in the above formula, where we have denoted the volume fraction of brine (saturation) in the pore space to be S w . The Gassmann-Wood model is a valuable theoretical tool and provides the low-frequency lowest bound for the composite elasticity of a partially saturated rock but the assumptions under which it is valid are often too restrictive. Alternatively a low-frequency upper bound for a partially saturated matrix is achieved by assuming an iso-strain condition for the fluids. The iso strain condition suggests volume averaging the fluid moduli of the fluids and replacing in Eq. (1) in what is often referred to as "patchy-saturation" model although a more precise description of patchy saturation is the so-called Hill-average discussed in Mavko and Mukerji (1998). In a recent theoretical work, Papageorgiou et al. (2016) established that these two boundsthe Wood and iso-strain limitscan be achieved by varying a constant of proportionality q between the pressure differentials of brine and CO 2 : One reason for departing from the pressure equilibrium given by q = 1 is indeed having a patchy distribution of fluids but other justifications, including capillary pressure or membrane stress effects, are also possible. In Papageorgiou and Chapman (2015) the relation of squirt flow to hysteresis phenomena is made more explicit but in any case, the effective fluid modulus arising from Eq. (2) has the form: where the patchy and Wood's bounds are reproduced at the endpoints = q K Kw CO2 and q = 1 respectively. The parameter q affects the partial saturation behaviour of a rock in a way shown in Fig. 1.
There are several reasons for introducing the parameter q. First of all, using Gassmann's model with the fluid given by Eq. (3) can consistently model situations where the partially saturated rock behaviour is between the patchy and Gassmann-Wood bound. In Papageorgiou et al. (2016) it was shown that the modulus of Eq. (3) reproduces the empirical model proposed by Brie et al. (1995) which has been observed in both lab and field data (i.e. Perozzi et al., 2016;Eid et al., 2015;Carcione et al., 2006;Kim et al., 2011). Papageorgiou and Chapman (2017) used Eq.
(2) to derive a fully frequency dependent model including the effects of squirt flow (Jakobsen and Chapman, 2009) with multiple fluids. They found the characteristic frequency of squirt flow is coupled to the parameter q of Eq.
(2) as well as the relative permeability of the matrix.
The characteristic frequency ω c scales relative to a reference frequency ω 0 assumed to be the squirt flow characteristic frequency of fully brine saturated rock: CO2 the viscosities and k k , w CO2 the relative permeabilities of CO 2 and brine respectively and q͠ is given in Eq. (3).
With the definition above, the frequency-dependent bulk modulus introduced in Papageorgiou and Chapman (2017) incorporates the effect of patchy saturation via the definition of the fluid (3) and characteristic frequency (4). The frequency dependent bulk modulus has the form: The parameters used in this expression are as follows: the symbols K m , ν refer to the mineral bulk modulus and the effective medium Poisson's ratio respectively. The symbols ϕ p , ϕ c are the spherical and microcrack porosity where the latter is given in terms of the microcrack aspect ratio r and crack density ϵ . The schematic behaviour of the bulk modulus when the experiment frequency is near the squirt flow frequency for brine is shown in Fig. 2 where the low frequency limits of Fig. 1 are also shown for comparison.
An important feature of the squirt flow model presented of Eq. (5) is the limit ω → 0, where the model reduces to Gassmann's model with the effective fluid introduced in Eq. (3) meaning the two equations can be used interchangeably when the frequency driving the fluid to squirt is assumed low.

Data analysis
To isolate the effect of the fluid solely on the rock physics parameters we go through a series of data processing and analysis steps aiming to make this effect more transparent. First of all we preprocess the dry V p , V s velocities of all three samples according to: . 1. The effective patch parameter affects the partial saturation behaviour of the bulk modulus but stays within well defined low-frequency bounds (ω = 0).
G. Papageorgiou et al. International Journal of Greenhouse Gas Control 74 (2018) 1-8 to obtain the dry bulk K d and shear G moduli. These dry measurements were performed with 100% air saturated samples under hydrostatic conditions. The confining pressure was chosen so that it matches the range of differential pressure states of the saturated measurements. We have plotted the variation of the dry moduli with confining pressure in Fig. 3 where we have also annotated the x-axis with a saturated equivalent "pore pressure". The dry moduli of all three samples show little sensitivity to confining pressure indicating that the confining pressure may be too small to introduce closure of microcracks. This is apparent from the small standard deviations around the mean dry moduli shown in Table 2 where we also document the solid and brine input parameters used.
Note that sample 3 has a higher shear than bulk dry modulus which may be due to its abnormally high porosity and in fact, the sample lies above the critical porosity value for sandstones quoted in Nur et al. (1998) which suggests that it may have the structure of pumice or similar volcanic rock. Sample 3 also has an unusually low, at least for sandstones, Poisson's ratio of ∼0.03.
Similarly we process the saturated measurements using Eq. (7) to obtain the saturated bulk modulus from the measurements of V p , V s . This time, the CO 2 density variation is taken into account in Eq. (7) and for each water saturation S w and pore pressure we use the composite density: The solid and brine densities ρ ρ , s w are assumed constant with pore pressure but the CO 2 density ρ CO2 changes according to the curves of Fig. 4 which are based on data from the National Institute of Standards and Technology (NIST) database. We have taken into account the temperature variability of these experiments as sample 2 was measured in slightly lower temperature to match the conditions of the Utsira reservoir.
The pressure dependence of CO 2 properties is also considered when calculating the composite fluid modulus K f . A constant modulus for brine K w (see Table 2) is mixed with CO 2 modulus K CO2 whose pressure dependence is taken from the NIST database and shown in Fig. 5. We use the mixing law of Eq. (3) with the appropriate values for each pressure step as the final input parameter. As mentioned, the mixing law of Eq. (3) incorporates the iso-stress (Wood's law used by Domenico, 1974) and iso-strain states as extremes of the parameter q which, in the current context will be a constrained fitting parameter. As argued in Papageorgiou and Chapman (2017), this parameter may itself be dependent on saturation or the wetting properties of the matrix. However, in this work we assume it to be constant in what can be thought as a zeroth order approach.
The porosity, dry and mineral modulus and a composite fluid modulus are the necessary input parameters we use to model the dependence of the partially saturated bulk modulus for each pressure step. Although care was taken so that the dry measurements are performed under the same state of differential pressure as the saturated measurements, in many cases it was not possible to adjust the pore pressure exactly. To account for this variability we apply a boxcar filter of width 1.3 MPa to the data centred around the reference pore pressure. In what follows, when we refer to fits corresponding to a particular pressure P 0 Fig. 2. The combination of varying the effective patch parameter q and the crack density ϵ affects the partial saturation behaviour of the bulk modulus but can lead to an increase beyond the low-frequency bounds (dashed lines) if the characteristic frequency takes intermediate values.

Fig. 3. Dry bulk (filled symbols) and shear (open symbols) moduli as a function
of confining pressure (the pore pressure equivalent axis is shown to assist comparison with the saturated measurements). There is only slight pressure dependence of these moduli justifying the use of a single crack density to model the squirt flow behaviour of each rock. Note that in the 45% porosity rock, the calculated shear modulus is higher than the bulk modulus.

Table 2
Measured rock physics parameters and water parameters used in the modelling of the sandstones (we have shown only the mean and standard deviation of the dry measurements of K d , G as their dependence from confining stress is small). The CO 2 properties are variable with pore pressure and are not shown here.  it should be understood that we refer to data measured within P 0 ± 0.65 MPa. Bearing in mind that the velocity measurements were performed with waves of frequency 600 kHz, we distinguish between two possibilities for the modelling of the data. One possibility is there is no squirt flow taking place between pores of different compliance in the samples. In this case, the characteristic squirt flow frequency of water ω 0 is assumed much larger than the experimental frequency and the behaviour of the bulk modulus is adequately captured by Gassmann's model with the composite fluid modulus of Eq. (3).
We can also assume the squirt flow effect is dominant and set the experiment frequency equal to the characteristic squirt flow frequency of water ω = ω 0 . In this case the fully frequency dependent squirt flow model of Section 3 will be used.
In both cases, the parameter q is seen as a fitting parameter. Since this parameter captures the effect of fluid patches on the saturated bulk modulus, we can investigate the correlation of q with the varying fluid pressure bearing in mind that values of q closer to unity correspond to uniform saturation and lower values correspond to more patchy fluid distribution. Note that the so-called patchy saturation limit depends on the pore pressure because of the variation in the elastic properties of CO 2 .
The characteristic squirt flow frequency scales inversely with viscosity so when the matrix is fully saturated with CO 2 we expect the squirt flow frequency to be ten to forty times higher. Therefore, irrespectively of what happens at intermediate saturations, we can hypothesise that the experiment frequency of 600 kHz is in the low frequency regime when the matrix is fully CO 2 saturated.
With this in mind, we use Gassmann's fluid substitution to calculate the expected CO 2 saturated moduli from the dry (air saturated) measurements since in either case the low-frequency assumption is expected to apply to the CO 2 saturated matrix. Noting that the CO 2 modulus is in the order of tens to hundreds of MPa, whereas the grain and dry moduli are of the order of GPa, the calculated values for the CO 2saturated moduli are close to the dry measurements. Their pore pressure dependence is dominated by the change in CO 2 properties rather than the stress dependence of the dry rock.
In any case, it should be noted that in calculating the moduli of the fully CO 2 saturated matrix, the stress dependence of the dry measurementshowever smallis taken into account. These calculated CO 2 saturated moduli are used in conjunction with the measured data points of the partially CO 2 saturated bulk modulus and constitute the dataset that is being fitted in the two cases below.

Low frequency case: ω = 0
In the low frequency case, Eq. (5) reduces to Gassmann's model of Eq. (1) with the composite fluid modulus given in Eq. (3). This model is fitted to the partially saturated data and inverted for the parameter q for each pressure step using a constrained Nelder-Mead nonlinear minimisation scheme.
The results, shown in Fig. 6a-c against the data, differ significantly between samples. The inverted values of the parameter q are shown against the pore pressure steps in Fig. 9 and suggest an exponential correlation.
More explicitly, in sample 1, the model performs well against the data for all pressure steps (Fig. 6a) and the inverted q correlates exponentially with pore pressure as is demonstrated in Fig. 9. This can be understood as a lowering of patch size with increasing fluid pressure.
For sample 2, the cluster of data points corresponding to 100% water saturation lie well below the predicted water-saturated Gassmann modulus for all pressures (Fig. 6b). The parameter q again correlates exponentially with pore pressure. A possible explanation for the 100% data/model mismatch may be the error-prone method of measuring  (3) with a q parameter fitted for different pore pressures. The variability at 0% water saturation (CO 2 saturated modulus) is dominated by the change of CO 2 properties with pressure. Hollow points indicate modelled rather than directly measured data. The light shaded area defines the range that can be covered under the uniform saturation assumption as the CO 2 properties change with fluid pressure. The darker shading defines the range that can be covered only when q < 1 which indicates patchy saturation. The black cross in the lower right hand corner shows the errors associated with the measurement of resistivity and velocity translated to saturation and bulk modulus. saturation using resistivity where the error in the resistivity measurement translates to 5% error in the measurement of the saturation whereas the measurement for the velocity is accurate to 0.3% of its value.
With sample 3, the situation is quite different (Fig. 6c). The inversion returns the minimum value for q corresponding to a Voigt fluid mixturea theoretical upper boundand still all data points lie well above that value. As the error in the measurement of the bulk modulus is very small, the measured increase indicates the existence of a different mechanism.
4.2. Squirt flow: ω = ω 0 When the experiment frequency is set to the characteristic squirt flow frequency of water, squirt flow effects become important and the fully frequency dependent Eq. (5) is used. Remembering that the definition of the dry modulus in (5) is different to that of (1), we need to make sure the two models are consistent with each other.
The model (5) is based on the Eshelby inclusion problem (Eshelby, 1957) and one of its assumptions is that each pore can be thought of as an isolated inclusion. This, however, is not the case when one considers high porosity sandstones. To remedy this issue, we solve for the Poisson's ratio of the effective medium defined in Eq. (5) so that Gassmann's Eq. (1) is satisfied identically. In essence, this implies that the effective medium of the matrix has a different Poisson's ratio for every pressure and porositya necessary concession when modelling sandstones with high porosity using inclusion models.
The values for the effective medium Poisson's ratio derived to ensure consistency of the squirt flow model with Gassmann's model together with additional input parameters for the model(crack density and aspect ratio) are given in Table 3.
The reasoning behind the choices of crack density is that no squirt flow is observed in samples 1, 2 hence it is expected they have relatively low crack density (∼1%), but sample 3 has a higher crack density to justify the abnormally high velocity measurements. Unfortunately it is not possible to invert the values for crack density without knowledge of the characteristic squirt flow frequency when the dry moduli are pressure insensitive so direct evidence (via, say, micro CT imaging) would be required to make more accurate estimates.
As was noted in Papageorgiou and Chapman (2017), a necessary ingredient for calculating the characteristic squirt flow frequency ω c of the partially saturated matrix is the relative permeabilities of brine/CO 2 as per Eq. (4). Here, we choose a simple symmetric relative permeability model This is a simplistic choice but consistent with standard relative permeability curves for CO 2 -brine systems published in reference literature (i.e. Benson et al., 2013). With this choice the characteristic frequency variation ω c at arbitrary saturation, relative to the reference frequency ω 0 (at full water saturation) can be written: In what follows, we fix ω 0 = 1 assuming the frequency f 0 of brine is the frequency of the experiment i.e. 600 kHz with a characteristic time constant τ 0 = 1.7 × 10 −6 s. Even though at the brine saturated end the characteristic frequency tends to one irrespective of pore pressure conditions, under partial saturation the changes in CO 2 viscosity affect the model. We therefore use appropriate viscosity curves of CO 2 shown in Fig. 7 with data taken from the NIST database.
With these considerations we fit Eq. (5) to the partially saturated data and invert for the parameter q for each pressure step. We use a constrained, nonlinear Nelder-Mead minimisation scheme.
The fitted bulk moduli corresponding to these results are shown in Fig. 8a-c. The low frequency limits discussed in Section 4.1 are shown at the same plots in grey for comparison. The bulk modulus variation in the dispersive model is almost indistinguishable from its low frequency counterpart for samples 1, 2.
To the contrary, when sample 3 is modelled on the frequency dependent approach (Fig. 8c) the fitted model performs significantly better when the frequency-dependent squirt flow mechanism is taken into account. Although this does not exclude other justifications for the abnormally high partially saturated modulus observed in this sample, it does provide a strong indication that the squirt flow effect may be responsible for this stiffening.
The inversion results for the parameter q as a function of pore pressure are compared in Fig. 9 for each of the three samples. The inversion for sample 3 lies almost exclusively at the patchy saturation limit but there appears to be a pressure dependence of the patch size for samples 1, 2. Irrespective of whether squirt flow is dominant or not, this dependence follows the same asymptotic behaviour. We have therefore postulated a dependence of the form and fit the variation of the patch parameter q from pressure using this model. The resulting fits are shown in Fig. 9 and the parameters a, b in Table 4 for each case.
Note there is a quantitative difference in the inverted patch size of Fig. 9. The parameter q differs little between the low and high frequency approach for sample 2 but for sample 1 the stiffening is due to the squirt flow when ω = 1, rather than fluid patches. But the residual difference between fit and data is so similar that it does not allow to distinguish which mechanismsquirt flow or patchy saturationis better suited to describe the measurement corresponding to that sample.

Discussion
Laboratory observations of the bulk modulus of rocks saturated by multiple fluids are often higher than the corresponding Gassmann-Wood predictions. Explanations of this stiffening are typically related to unevenness in the induced fluid pressures, either Table 3 Demanding consistency of Eq. (5) with Gassmann's model requires a Poisson's ratio for the effective medium that is dependent on porosity and the input parameters crack density and aspect ratio shown above.  G. Papageorgiou et al. International Journal of Greenhouse Gas Control 74 (2018) 1-8 between different parts of the pore-space or between the different fluids. Differences in induced pressures between the fluids are most commonly explained by assuming fluid patches, but following recent theoretical developments we capture the effect via an "effective patch parameter" q which provides a scaling constant between induced fluid pressures (Eq. (2)). We measured elastic wave velocity and resistivity during CO 2 flooding of brine saturated synthetic sandstone plugs and resistivity was converted to degree of saturation. Errors were small for velocity measurements but in the conversion to saturation the resistivity measurement introduced an error ranging from ∼3% for the fully saturated, to ∼5% for the partially saturated measurements. Great care was taken to fully saturate the samples by continuously flushing them with brine while simultaneously subjecting them to pressure cycles to dissolve any remaining air (a similar approach has been taken in Blake and Faulkner, 2016, for example). Nonetheless, the samples were constructed under vacuum conditions and the presence of empty, disconnected porosity cannot be ruled out. We also measured attenuation of the three samples, and found high (0.05-0.2) and generally fluid independent values consistent with an interpretation in terms of scattering. As the attenuation due to squirt flow peaks to around 0.07, the attenuation due to scattering dominates and any attempt to model this attenuation would be overly ambitious.
Our work emphasises the importance of treating squirt and patch phenomena in a unified way, both for understanding measurements and extrapolation to seismic frequencies. The influence of fluid mobility on squirt effects is controlled by relative permeability. This leads, at intermediate frequencies and high water saturations, to a potentially counterintuitive decrease in bulk modulus with increasing water saturation. We believe there is some weak evidence for this in the data corresponding to sample 3 where a small reduction in the bulk modulus is observed when the saturation is increased to 100% water. This reduction is captured by the fitted squirt flow model.
Recent discussion (Thomsen, 2017;Duranti, 2017) has focused on the potential importance of the deviation between the approaches to fluid substitution taken by Gassmann and Brown and Korringa (1975), based on the potential of an independent "pore space compressibility" parameter. We considered whether the discrepancy could provide an explanation for the behaviour of the highest porosity sample. The Brown-Korringa pore space compressibility can lead to a higher predicted bulk modulus but this is only pronounced when the dry and mineral moduli are of similar magnitude. We concluded that this concept was not able to explain the measurements of our highest porosity sample.
The outcome of our work is an empirical fluid mixing law for brine-CO 2 saturation in clean sandstones which, in principle, is valid in the low frequency limit. In contrast to previous work, behaviour is controlled by fluid pressure rather than a patch size. We find it intuitive that increasing fluid pressure leads to more uniform mixing of fluids.
As discussed by Papageorgiou and Chapman (2017), the effective patch parameter, q, is likely to be influenced by wetting phenomena (3) with a q parameter fitted for different pore pressures. The variability at 0% water saturation (CO 2 saturated modulus) is dominated by the change of CO 2 properties with pressure. Hollow points indicate modelled rather than directly measured data. The light shaded area defines the range that can be covered under the uniform saturation assumption as the CO 2 properties change with fluid pressure. The darker shading defines the range that can be covered only when q < 1 which indicates patchy saturation. The black cross in the lower right hand corner shows the errors associated with the measurement of resistivity and velocity translated to saturation and bulk modulus. Fig. 9. Inversion of the parameter q of Eq. (3) and the fit of samples 1, 2 (denoted s − 1, s − 2) to an exponential pore pressure model. The fit for ω = 0 is with a dashed, whereas for ω = ω 0 with a solid line. The fitted empirical model describing variation of q with pore pressure is constructed by observing the asymptotic behaviour of the data. Sample 3 (denoted s − 3) was not fitted as its inversion for q coincides with the theoretical minimum curve K Kw CO2 as CO 2 modulus changes with fluid pressure. Table 4 Fitting parameters a, b of the empirical model = + ( ) q P ( ) exp a P b P 3 shown in Fig. 9.

Conclusions
Ultrasonic measurements of the bulk modulus of three synthetic sandstones of different porosities saturated by brine-CO 2 mixtures can be understood in terms of coupled patch and squirt mechanisms. The two lowest porosity sandstones are well modelled by a patch model in the low frequency limit where squirt flow does not operate. Moduli of the higher porosity sample are higher than predicted by low frequency theories and consistent with squirt effects. In all cases, account must be taken of the pressure dependence of the modulus and viscosity of CO 2 .
The results indicate that the effective patch parameter follows a simple exponential law with respect to the fluid pressure. This suggest that in many practical situations the deviation from the Gassmann-Wood theory will be more pronounced for low fluid pressures. We suggest that interrogation of such laboratory measurements with a coupled patch and squirt model will point the way to improved practical methodologies for brine-CO 2 fluid substitution.