The effect of ﬂuid compressibility and elastic rock properties on deformation of geothermal reservoirs

A geothermal reservoir deforms when the extraction of pore ﬂuid exceeds reservoir recharge, causing a decrease in pore pressure. The magnitude of this deformation is related to the amount of pore ﬂuid that is extracted. Assuming compressible material properties in a homogeneous reservoir, we derive an expression for the ratio of reservoir volume change per extracted ﬂuid mass. We show that this ratio depends on a number of parameters, notably the compressibilities of reservoir rock and pore ﬂuid. We apply the obtained relationship to three different geothermal areas (Hellisheidi, Reykjanes and The Geysers) to illustrate under which circumstances the relation between reservoir deformation and the amount of extracted ﬂuid is able to help us learn more about reservoir conditions. We ﬁnd that the ﬂuid compressibility, depending on whether the system is single-phase or two-phase, may explain large differences in estimates of reservoir volume changes per mass of extracted ﬂuid.


I N T RO D U C T I O N
Extraction of geothermal fluids can cause significant surface deformation, sometimes with magnitudes reaching the metre scale as observed at the Wairakei field in New Zealand (e.g. Hatton 1970;Allis 2000). Satellite-based geodesy has made it possible to detect millimetre-scale surface deformation, resulting in new opportunities to study deformation of reservoir rocks. Many studies have been published in recent years highlighting deformation caused by geothermal fluid extraction and injection observed using In-SAR and/or GNSS technology (e.g. Mossop & Segall 1997Eysteinsson 2000;Fialko & Simons 2000;Vasco et al. 2002;Foxall & Vasco 2003;Keiding et al. 2010;Sarychikhina et al. 2011;Ali et al. 2016; Barbour et al. 2016;Drouin et al. 2017;Juncu et al. 2017). These geodetic observations help us to quantify the deformation around fluid reservoirs. By using simple analytical deformation models, we can estimate a reservoir's volume change due to contraction or dilatation, as well as the depth and size of the reservoir. The volume change can be compared to measured quantities such as pressure or temperature changes, and to the amount of fluid extraction/injection. In this study, we focus on the relationship between fluid extraction and the volume change of the reservoir for cases where deformation is driven by changes in pore pressure.
When comparing known volumes of fluid extraction or injection to geodetically estimated changes in reservoir volumes, many studies find that the former volume change exceeds the latter, sometimes by orders of magnitude (e.g. Eysteinsson 2000;Keiding et al. 2010;Sarychikhina et al. 2011;Juncu et al. 2017;Liu et al. 2018). To examine what causes this discrepancy, we employ the framework of poroelasticity to investigate how fluid and reservoir volume changes are related and which parameters control this relation. We describe a simple methodology to help examine plausible parameter combinations and explain the magnitude of deformation for a given amount of extracted fluid. To this end, we derive an equation that shows that the ratio of fluid-to-reservoir volume change depends on various reservoir properties, like the bulk compressibility of the reservoir rock and the compressibility of the pore fluid, as well as the fluid exchange with the surrounding formations (i.e. recharge in case of fluid extraction, leakage in case of fluid injection).
The reservoir properties on which the magnitude of deformation depends, however, are often poorly constrained. The fluid compressibility in a geothermal reservoir, for instance, can have a wide range of values depending on whether the pore fluid is single-phase steam or water, or consists of both phases (see Section 2.1.1). In a two-phase system, the fluid compressibility depends on the relative amounts of steam and water, which is often not well known. The bulk rock compressibility at field scale can be difficult to constrain, as well. Tabulated values found in well-known literature (e.g. Turcotte & Schubert 2002) are often high, in the tens of GPa. But such values are usually based on seismological data and may not always be suitable for large-scale deformation problems. This may be partly because seismologically derived elastic moduli are typically larger Reservoir deformation and material properties 123 than static moduli that prevail in deformation problems (e.g. Heuze 1980;van Heerden 1987;Eissa & Kazi 1988), but also because of local factors that affect rock strength and rigidity. Elastic moduli depend on effective stress, (Zimmerman et al. 1986), implying a depth dependence that needs to be considered for deformation stemming from shallow reservoirs. Hydrothermal alterations of minerals can lower the rigidity of rock formations Frolova et al. 2014), which may affect deformation due to geothermal fluid production. The effect of mineral alterations and effective stress on elastic moduli could be an indicator that elastic moduli of shallow formations in geothermal areas may be relatively low but, unfortunately, estimates from laboratory tests are only rarely available to confirm this. For the Wairakei-Tauhara geothermal field in New Zealand where very large amounts of subsidence have been observed (up to 15 m over 5 decades ;Bromley et al. 2013), Pender et al. (2013) report values of the constrained modulus K v (see eq. B3) obtained by core sample compression tests as low as 0.1-1 GPa for samples from depths of around 400 m. When values of elastic moduli are required as input for deformation models-or used to interpret modelling results-of various other geothermal fields, they are often assumed to exceed 10 GPa (e.g. Fialko & Simons 2000;Keiding et al. 2010;Juncu et al. 2017). This way, it seems that possible values for elastic rock properties can range over orders of magnitude, posing a challenge for deformation modelling in geothermal areas.
Using examples from different geothermal fields, we investigate the link between the magnitude of deformation of geothermal reservoirs and the reservoirs' properties in order to identify which combinations of property values are able to explain a given magnitude of deformation. We show that it is possible to narrow down parameter ranges using our methodology. Considering the aforementioned problem in constraining the reservoir properties using other methods, for example, laboratory measurements, we can view man-made reservoir deformation as an opportunity, a large-scale experiment, that we can use to learn more about the mechanical behaviour of porous and fractured rocks at field scale. The advantage of the reservoir setting stems from a priori knowledge about changes in pore fluid volume and pore pressure. This is information that is generally not available in cases of natural deformation, for example, driven by magmatic or tectonic processes. Measurements of changes in pore pressure and the amount of extracted fluid can be used to better constrain reservoir conditions and estimate material properties, thus improving our understanding of reservoir processes.

P O R E F L U I D A N D R E S E RV O I R V O L U M E C H A N G E S
Fluids in rock formations are stored in pores and fractures, where the fluid's pressure, or pore pressure, supports the rock matrix, counteracting compressional normal stresses. A change in pore pressure can be achieved by a change in volume of pore fluid, which, according to the theory of poroelasticity (Biot 1941), produces stress and deformation in the rock formation. The theory also implies that the change in stress in the rock in turn produces a change in pore pressure. This two-way coupling makes it difficult to obtain analytical solutions to the equations of poroelasticity. To obtain an analytical solution the two-way coupling can be turned into one-way coupling (changes in pore pressure produce changes in stress, but not vice versa) under the assumption of uniaxial strain and constant vertical stress. Using this approach, Geertsma (1966Geertsma ( , 1973 derived an expression that relates changes in pore pressure within a fluid reservoir to surface deformation. The reservoir is assumed to be a flat axisymmetric finite structure separated from its surroundings by an impermeable barrier, see Geertsma (1973), and it is assumed that material properties are homogeneous (i.e. the properties of reservoir rock and surrounding formations are the same). Fig. 1 depicts a sketch of the concept. For more detailed explanations on poroelastic theory we refer the reader to Wang (2000) whose notation we follow in this study.

Model equation
We use the analytical solution of Geertsma (1966Geertsma ( , 1973 for the vertical compaction of a fluid reservoir, and retain the same set of assumptions to obtain a relation for fluid and reservoir volume change. According to Geertsma's work, reservoir volume change V r of a reservoir with volume V r , under the condition of uniaxial (vertical), constrained strain (Eshelby 1957), depends on a change in pore pressure p, such that (following Geertsma 1966): where c m = α/K v is the uniaxial poroelastic expansion coefficient, K v is the uniaxial drained bulk modulus (also known as the constrained modulus) of the reservoir and α is Biot's poroelastic expansion coefficient (Geertsma 1966(Geertsma , 1973Wang 2000). Material properties are assumed to be homogeneous and isotropic. Next we consider the relation between change of pore fluid volume V f and pore pressure change p. The amount of fluid that is moved from/into a control volume given a change in pore pressure depends on both solid and fluid compressibilities. As noted by Geertsma (1957), for a porous solid, three different compressibilities need to be distinguished: rock matrix compressibility 1/K S , rock bulk compressibility 1/K and pore compressibility 1/K φ . They are linked through a number of poroelastic coefficients, see for example, eq. (B2; Appendix B). The combined effect of the different compressibilites can be expressed through the specific storage coefficient, S. Specific storage is defined as the volume of released water per bulk volume and change in hydraulic head (Wang 2000). S can be defined in terms of pore pressure change for a volume V r as (see Green & Wang 1990;Wang 2000): Recasting eq. (2) in terms of the change in fluid volume gives By taking the ratio of eqs (1) and (3) we are able to eliminate the V r p term and obtain an expression for the ratio V f / V r : Retaining the assumption of uniaxial strain and constant vertical stress, the uniaxial specific storage coefficient can be written in terms of material compressibilites as (Wang 2000, chapter 3.3.3): where 1/K φ is the pore compressibility, 1/K f is the compressibility of the pore fluid, α is Biot's expansion coefficient and φ is the porosity. We can simplify this equation by assuming that the pore compressibility 1/K φ equals that of the rock matrix 1/K S (eq. B2). This assumption holds for the case of a solid phase consisting of a single constituent (Berge & Berryman 1995;Wang 2000), displacements. V f, t and M f, t are the total (gross) volume and total mass of extracted fluid. V r is the volume change of the reservoir due to deformation. V f is the net volume of extracted fluid, that is, V f, t minus recharge. K is the reservoir bulk modulus, K f is the bulk modulus of the pore fluid and K S is the bulk modulus of the solid matrix. α, ν and φ are Biot's coefficient, Poisson's ratio and porosity, respectively. and is commonly applied to idealize numerical models (see e.g. Zienkiewicz & Shiomi 1984). Using this assumption and inserting eqs (5) and (B3) into eq. (4), we obtain where ν is Poisson's ratio. A problem with applying eq. (6) is that the net change in pore fluid volume ( V f ) is generally not known.
The value of V f depends on the fluid balance in the reservoir, that is, how much fluid is extracted and how much of the extracted fluid is replenished from surrounding formations or by injection. What we know instead, is the total volume of extracted fluid V f, t , and that V f equals V f, t minus recharge into the reservoir. Hence, the value of V f represents only an effective fluid volume change, which we can relate to the known total extracted fluid volume V f, t by introducing a factor χ , so that where 0 < χ ≤ 1 denotes the fraction of V f, t that is causing the observed deformation. Within the assumptions of our model, homogeneous material parameters and small vertical extent of the reservoir, this χ essentially depends only on the recharge into the reservoir, which means that in this case, 1 − χ is the fraction of V f, t that is recharged. Combining eqs (6) and (7) can be written as This shows that the volume ratio of fluid flow and reservoir deformation depends on six parameters. Under a set of simplifying assumptions we can illustrate a limiting behaviour of eq. (8). If pore fluid and solid matrix are incompressible, or their compressibilities are much smaller than the bulk compressibility of the reservoir (K K f , K << K S , α = 1), the term in the square brackets in eq. (8) becomes unity so that the ratio of volume changes depends only on χ : V f, t / V r = 1/χ . In most cases, however, this will not be very realistic, in particular because K K f is not a common occurrence.
In high-temperature geothermal reservoirs, we have to consider the additional problem of two-phase conditions: water and steam may both be present within the reservoir while their relative fractions of the total fluid volume are unknown. Therefore, it may not be possible to estimate the volume of extracted fluid V f, t . We can circumvent this problem if we use mass instead of volume to quantify the fluid loss from the reservoir. The fluid volume in eq. (8) can be replaced by the total mass of extracted fluid M f, t by where ρ f is the bulk fluid density within the reservoir, which depends on the relative fractions of water and steam. We can take eq. (9) and rearrange eq. (8) in terms of the reservoir volume change per total mass extraction V r / M f, t : (10) To find an upper limit for this ratio, we can consider the case of an incompressible fluid (K f = ∞), incompressible solid grains (K S >> K , α = 1, see Section 2.1.2), no recharge (χ = 1) and a fluid density of ρ f = 1000 kg m −3 , for which we obtain V r / M f, t = 10 −3 m −3 kg −1 . Any values lower than this are caused by the combined effect of the parameters on the right-hand side of eq. (10).

Two-phase fluid compressibility
While for a single-phase fluid the compressibility (the inverse of the fluid's bulk modulus K f ) is readily obtained, for a two-phase fluid this is, not surprisingly, more difficult. Grant & Sorey (1979) showed a convenient way of simplification, however: in a two-phase water-steam system, where both phases are in thermal equilibrium with each other and the surrounding matrix, changes in pressure lead to phase changes that result in significant volume change due to the large density difference between liquid and gas phases. This way, the effective fluid compressibility is dominated by the mass exchange between the phases and exceeds the effect of the respective singlephase compressibilities, which we can thus neglect (Grant & Sorey 1979). The value of the effective fluid compressibility depends on the thermal and mechanical properties of rock and pore fluid, the effect of latent heat and the relative fractions of water and steam in the system. Grant & Sorey (1979) give the following expression for the two-phase compressibility (see also Grant 2011;Watson 2013):

125
where ρC is the reservoir heat capacity, defined as (Moench & Atkinson 1978): Here, ρ denotes density, c p specific heat, L the latent heat of vaporization, f l the volume fraction of liquid and T the temperature (in Kelvin). Subscripts s, w and r denote steam phase, water phase and rock, respectively. The density of the water-steam system is a function of the relative fractions of liquid and steam, given as (Grant & Sorey 1979): To illustrate the difference in compressibility values we can consider an example case of a fluid at 250 The fluid is liquid water at a pressure of 5 MPa and has a compressibility of 0.001 MPa −1 , it is steam at 3 MPa and has a compressibility of 0.5 MPa −1 , and using eqs (11)-(13), a water-steam system on the saturation curve at 4 MPa has the highest compressibility of 6 MPa −1 . We can plug eqs (11)-(13) into eq. (10) and obtain an expression for the reservoir volume change per total mass extraction V r / M f, t in terms of the liquid fraction f l as well as material properties of rock mass and pore fluid:

Model parameters
Eqs (8) and (10) contain a number of variables that depend on the physical properties of the reservoir. Here, we give a short overview of these parameters as well as the range of values that they commonly assume. For properties of water and steam that depend on temperature and pressure (density, compressibility, heat capacity) we use the tables given in Wagner & Kretzschmar (2007).
Reservoir bulk modulus K. The reservoir bulk modulus (sometimes referred to as wet rock bulk modulus or drained bulk modulus) can vary approximately from the order of ∼10 −1 to the order of ∼10 1 GPa. Its value depends on rock type, temperature, degree of fracturing, confining pressure and degree and type of hydrothermal alterations (e.g. Zimmerman et al. 1986;Lynne et al. 2013;Hassanzadegan et al. 2014). It can be estimated by laboratory tests on core samples (e.g. Blöcher et al. 2014) or at field scale, using geodetic data of surface deformation if reservoir pressure changes can be constrained (see Section 3.1.1).
Pore fluid's bulk modulus K f . The fluid's bulk modulus is known in cases where the fluid composition is known. For the case of a two-phase geothermal field where both water and steam occur, K f can be calculated in terms of the fraction of liquid it contains (see Section 2.1.1). At a pressure of 5 MPa, K f of liquid water is around 1000 MPa between temperatures of 200 and 250 • C. For steam at 5 MPa, K f is much lower, around 4 MPa at temperatures between 275 and 325 • C. For a two-phase system at 5 MPa and 264 • C (which is on the boiling point curve), it is the lowest, approximately 0.3 MPa, calculated using eq. (11) for 0 < f l < 1, φ = 0.1, ρ r = 2500 kg m −3 , c p, r = 1.0 kJ kg −1 • C −1 , and the properties of water and steam according to the given temperature and pressure (Wagner & Kretzschmar 2007).
Biot's coefficient α. Biot's expansion coefficient α controls how much of the pore pressure is being transferred onto the solid matrix. It can be expressed in terms of the reservoir bulk modulus K and the solid matrix bulk modulus K S , the relation is given in eq. (B2). Theoretically, it can take any value between 0 and 1, where the limiting cases are an incompressible (relative to the bulk reservoir strength) matrix (K << K S ⇒ α = 1) and the reservoir bulk modulus being equal to the matrix bulk modulus (K = K S ⇒ α = 0). It is commonly assumed, however, that the value of α falls between 0.5 and 1 (e.g. Zienkiewicz & Shiomi 1984). We also know that the bulk moduli of the mineral components, which control the value of K S , are with few exceptions above 30 GPa (e.g. Mavko et al. 2009). This can be used to limit the range of possible values of α. A value of K < 3 GPa, for instance, in many cases requires that α > 0.9.
Porosity φ. The porosity of rocks varies with rock type and confining pressure (i.e. depth) and typically ranges from 0.01 to 0.5 (see Manger 1963). For hydrocarbon reservoir sandstones and carbonates at 1000 m depth, these values are generally below 0.35 (Ehrenberg & Nadeau 2005). We assume that geothermal reservoirs made up of, for example, volcanoclastic or fractured crystalline rocks do not have higher porosities. For geothermal reservoirs in Iceland, φ = 0.1 is a common assumption for modelling reservoir processes (see e.g. Axelsson et al. 2015). These values can be better constrained by geological information from drill cores or other means.
Poisson's ratio ν. Poisson's ratio is the ratio of transverse strain to axial strain in an elastic medium. For rocks, the values of ν fall most commonly between 0.1 and 0.4 (Gercek 2007).
Fraction of effective fluid volume change χ . The parameter χ quantifies the ratio of the net fluid loss V f in the part of the reservoir exhibiting detectable deformation, to the total volume of extracted fluid V f, t (eq. 7). In this study, we assume that we are able to detect deformation from the whole reservoir, which means that χ only relates to reservoir recharge, which is 1 − χ . The value of χ ranges from 0 (complete recharge i.e. no deformation) to 1 (no recharge). It is possible to constrain reservoir recharge using gravity measurements to infer mass changes in the reservoir (e.g. Guðnason et al. 2015).
Liquid volume fraction f l . The liquid volume fraction is limited by f l = 0 (pure steam) and f l = 1 (pure water). In a geothermal twophase reservoir, it is often not well constrained, because of natural and induced variations in pressure and temperature in space and time. In our model this parameter is, as are all other parameters, averaged over the whole reservoir.

E X A M P L E S : D E F O R M AT I O N A N D F L U I D E X T R A C T I O N AT G E O T H E R M A L F I E L D S
We can apply the equations we derived in Section 2 to a set of examples where both fluid extraction rates and reservoir volume change are known. Due to parameter trade-off, eq. (10) requires additional constraints to allow any kind of quantitative interpretation. While for some parameters, constraints may be derived from prior information (e.g. examination of drill cores), for others, researchers or reservoir engineers have to rely on experience-based assumptions or adopt values from comparable cases.
We will use three examples to demonstrate how the ratio of deformation per fluid extraction can be used to constrain reservoir parameters and to learn more about reservoir conditions. To this end we employ the random walk Metropolis algorithm (see Metropolis et al. 1953;Gelman et al. 2014) to sample from the Bayesian posterior probability distribution of parameters. To this end, we assume that uncertainties on V r / M f, t are normally distributed with a standard deviation of 30 % of our estimate (based on an approximation for the Hellisheidi results, see Supporting Information Fig. S3; neglecting uncertainties in the mass extraction value). We take samples from uniform probability distributions over parameter ranges discussed in Section 2.1.2, but use parameter values (or ranges) from the specific study area if these are available. Temperatureand pressure-dependent water and steam properties are taken from Wagner & Kretzschmar (2007). Given the simplifying assumptions that imply homogeneous material parameters, it has to be noted that the resulting range of values depicts an average range in both space and time for the part of the reservoir that is deforming. Here, we study three examples where production of geothermal fluids caused observable surface deformation. In all presented cases, it has been assumed that fluid pressure changes are driving the deformation of the reservoir and thermoelastic deformation has been neglected.

Hellisheidi, Iceland
The study of Juncu et al. (2017) reports surface deformation due to production at the Hellisheidi high-temperature geothermal field in southwest Iceland between 2012 and 2015. For this case, we know the extraction rates from the reservoir in terms of mass rate, which was on an average M f, t = 38 × 10 9 kg yr −1 during the same time interval (Gunnlaugsson 2016). Of the extracted mass, ∼60 % was reinjected back into the crust, although it is not certain how much of this fraction flowed back into the reservoir (Juncu et al. 2018).
In the study of Juncu et al. (2017), a shear modulus of 10 GPa and a Poisson's ratio of 0.25 was assumed, which implies a bulk modulus of 17 GPa. However, pressure changes in the Hellisheidi reservoir have been inferred from borehole measurements (Haraldsdóttir 2014; Juncu et al. 2017) and can be used as an upper bound for an average reservoir pressure decrease. Therefore, we are able to estimate the formation's bulk modulus using the geodetic data reported by Juncu et al. (2017). We also invert for the reservoir volume change using a different reservoir model than in the original study, which we argue has a more suitable geometry. For all parameters that we invert for we estimate the posterior probability density functions, shown in the Supporting Information.

Inversion of displacement data using Geertsma's poroelastic model
Using GNSS and InSAR surface displacement data (2012)(2013)(2014)(2015) from the Hellisheidi geothermal area and observed reservoir pore pressure changes (Juncu et al. 2017) we can perform a joint optimization of the two displacement data sets in order to constrain the bulk modulus of the formation. We use the disk-shaped reservoir model given by Geertsma (1973;see also Segall 1992). Surface displacements for such a reservoir of radius R, vertical thickness H, depth D, experiencing a pressure change p at a point at radial distance r from the reservoir, are given as (Geertsma 1973): where J 0 and J 1 are Bessel functions of the first kind, of zero and first order, respectively. A rigorous derivation for eq. (15) can be found in Segall (1992). We solve the integrals in eq. (15) numerically using Matlab's integral (Shampine 2008) and besselj (Amos 1986) functions. As an upper bound for the average reservoir pressure change we use the average observed pressure decrease range at Hellisheidi of 0.3 MPa yr −1 (Juncu et al. 2017). We use a Bayesian Markov chain Monte Carlo type inversion approach (the Catmip algorithm, Minson et al. 2013) to find the optimal reservoir source parameters as well as probability distributions of the parameter values ( Fig. 2 and Supporting Information Figs S1 and S2). As in Juncu et al. (2017) we correct the data for the plate spreading signal and invert for three local deformation sources, two for the geothermal fields at Hellisheidi and Nesjavellir and a deep source in eastern Hengill. Results for all deformation sources are given in the Supporting Information. For the purpose of this study, however, we are only concerned with the Hellisheidi source. For Hellisheidi, we obtain a 90 % confidence interval (CI) for the bulk modulus K of 1.3-6.4 GPa for a Biot-coefficient range of 0.9 ≤ α ≤ 1. We find that the deformation originates from a depth of around 1.3 km (90 % CI ranging from 0.7 to 2.0 km) with a source thickness of 0.7 km (90 % CI: 0.5-1.0 km). The reservoir volume change is V r ≈ 2.1 × 10 5 m 3 yr −1 (90 % CI: 1-4 × 10 5 m 3 yr −1 ), which results in a ratio of reservoir volume change per mass extraction of V r / M f, t ≈ 6 × 10 −6 m 3 kg −1 , or, using the CI, gives a 90 % confidence range of 3-10 × 10 −6 m 3 kg −1 (neglecting uncertainties in the mass extraction rate). The optimal result for Hellisheidi in the original study (Juncu et al. 2017), using a spheroid model, is a vertical thickness (at the source centre) of 2.5 km (90 % CI: 2.0-3.0 km) and a reservoir volume decrease of 4 × 10 5 m 3 yr −1 , which agrees with the 90 % confidence limits of the reservoir volume change estimated here. A key difference between source geometries used in this study and in that of Juncu et al. (2017), is that in the latter study the source model was a spheroid with the long axis being orientated near horizontally, meaning that the vertical and the short horizontal axis need to be of same length. The geodetic data, however, constrain the horizontal extent of the source better than the vertical, so that this can lead to a bias of the source thickness for near horizontal spheroidal models. The results for source thickness and source volume presented here are about one third of the result from Juncu et al. (2017) and we believe the new result is a more realistic representation of the vertical extent of the part of the geothermal reservoir that causes measurable surface deformation.

Relating V r / M f, t to fluid properties
Having obtained constraints on the reservoir's bulk modulus we can explore the possible parameter combinations that yield V r / M f, t ≈ 6 × 10 −6 m 3 kg −1 . The pressure in the central region of the Hellisheidi reservoir is approximately 7 MPa at depths between 1 and 1.5 km (see Haraldsdóttir 2014) and the reservoir temperature is at the boiling point curve, which is 286 • C at this pressure. We expect the system to be two-phase, because a steam cap has formed in the top region of the reservoir due to the pressure drop (Gunnarsson et al. 2011). However, we test the cases of pure steam and pure water, as well. We use eq. (10) and a random walk Metropolis algorithm to find the parameter combinations that yield values of V r / M f, t that fall in the 90 % confidence range of 3 × 10 −6 -10 × 10 −6 m 3 kg −1 , where parameter combinations that result in values close to the optimal value of V r / M f, t ≈ 6 × 10 −6 m 3 kg −1 are sampled more frequently.
We first test the single-phase case for pure liquid and pure steam just below and above the boiling point curve. For the pure liquid case, we assume a density of ρ f = 750 kg m −3 and a fluid bulk modulus of K f = 500 MPa; for pure steam we use ρ f = 35 kg m −3 and K f = 5 MPa. To constrain the value of Biot's expansion coefficient α, we make use of the results from the inversion of geodetic data (Section 3.1.1), where we obtained a 90 % CI of K of 1.3-6.4 GPa. For such low values of K, we assume that 0.8 < α < 1 (see Section 2.1.2). We also use a range of 0.01 ≤ φ ≤ 0.2, 0 ≤ χ ≤ 1 and 0.1 ≤ ν ≤ 0.4. We find that these conditions require a very low value of χ , that is, recharge fractions of greater than ∼80 % for steam or ∼90 % for water (for K ≤ 7 GPa), see Fig. 3.
Next, we test the assumption of a two-phase pore fluid. For this case, the parameter values for these conditions are shown in Table 1. Parameter value ranges of α, K, φ, χ and ν are the same as above. Using eq. (11), we can calculate a possible range of K f of approximately 0.1-3 MPa (depending on porosity φ and the volume fraction of liquid f l ). We show in Fig. 4 how the reservoir bulk modulus K and f l correlate, and find that f l is very sensitive to K and practically unconstrained for K < 1 GPa. However, for K = 3 GPa, the optimal value determined through the inversion of geodetic data (Section 3.1.1), we find that most of the values fall around f l ≈ 0.15, which translates to a liquid mass fraction of 80 %.
Given the very high required recharge fraction for the singlephase cases and that we know that both water and steam occur in the reservoir, it seems more likely that the reservoir behaves like a two-phase system. It also shows that with a more accurate estimate of K it may be possible to constrain the ratio of water-to-steam in the deforming part of the reservoir. Keiding et al. (2010) describe deformation between 2005 and 2008 at the Reykjanes geothermal field, located around 80 km southwest of the Hellisheidi field and commissioned in 2006. For this time interval they report a total mass extraction of around 58 Mton and find an estimated reservoir volume change of about 2.1 × 10 6 m 3 . Hence, the extracted mass per reservoir volume change is V r / M f, t ≈ 4 × 10 −5 m 3 kg −1 .

Reykjanes, Iceland
It has been suggested that due to production-induced pressure drop a steam cap developed in the reservoir . Since the deformation study of Keiding et al. (2010) covers only the beginning of the geothermal fluid production, it is not clear whether or not the amount of steam in the reservoir is significant at that time. We test both the effect of a purely liquid and a two-phase pore fluid. Given the low estimate of K, for both cases we assume that α > 0.5, see eq. (B2). For simplicity, we assume the pore fluid's liquid phase is pure water, although in reality it is hydrothermally altered sea water . First, we test the case of a purely liquid pore fluid, given material parameters as discussed above and a porosity of 0.01 ≤ φ ≤ 0.35. We find that V r / M f, t = 4 × 10 −5 m 3 kg −1 can be explained by a wide range of combinations of χ and K, as shown in Fig. 5(a). For geothermal reservoirs in Iceland a common assumption is φ = 0.1, which we can apply to see if it is possible to better constrain the parameters. The trade-off between K and χ after adding the constraints is shown in Fig. 5(b). We can compare the range of χ to a recharge estimate by Guðnason et al. (2015) for the time interval 2008-2010, which is 30-50 % (i.e. 0.5 < χ < 0.7, the area between the dashed black lines in Fig. 5). For this range of χ , we can see that for the case of unconstrained φ (Fig. 5a) the results have a high likelihood for K > 5 GPa, while for φ = 0.1 the regions of higher likelihood are where recharge rates are higher (or χ is lower) than the estimate by Guðnason et al. (2015).
We test whether it is possible to achieve V r / M f, t of ∼4 × 10 −5 m 3 kg −1 by two-phase reservoir conditions. We assume that average reservoir pressure is 10 MPa and that reservoir conditions are on the boiling point curve (temperature of 310 • C). The parameters for these conditions are shown in Table 1. The range of plausible combinations of liquid volume fraction f l and reservoir bulk modulus K are any 0 < f l < 1 and K < 4 GPa. Since the studied time interval of 2005-2008 only covers the beginning of production at Reykjanes, it may be likely that the amount of steam is low in the reservoir. For a liquid mass fraction of at least 90 %, a volume fraction f l > 0.4 for reservoir conditions at Reykjanes, our results indicate that K needs to be below 1 GPa (Fig. 6).
These results indicate that both two-phase and pure liquid compressibilities can possibly explain the observed deformation, albeit in both cases with restrictions. A possible reason for the ambiguity in results is that the reservoir conditions may have changed from liquid to two-phase during the studied time interval, see Section 4.

The Geysers, California, USA
Large-scale geothermal power production at The Geysers began in 1960, with production being strongly increased from 1972 to the late 1980s (Barker et al. 1992). Mossop & Segall (1997 published two studies presenting deformation obtained by GPS and levelling  1 (b). Colour scale depicts the relative probability of parameter combinations. The vertical dashed black lines limit 0.5 < χ < 0.7 (recharge between 30 and 50 %, see Guðnason et al. 2015). Figure 6. Bivariate histogram of liquid saturation f l and reservoir bulk modulus K for a two phase pore fluid. Calculated using eq. (10) for the deformation at the Reykjanes geothermal field (Keiding et al. 2010). Colour scale depicts the relative probability of parameter combinations. observations at The Geysers between 1977 and 1996. From their estimates of volumetric strain and observations of reservoir pressure changes, they estimated K ≤ 4.6 GPa (Mossop & Segall 1999). We can derive an estimate range for values of α using eq. (B2) and their estimate for K and the value they assumed for the solid grain bulk modulus, K S ≈37 GPa. This gives 0.9 < α < 1, which we will apply here. Gunderson (1992) and Barker et al. (1992) report reservoir matrix and fracture porosities of 1-6 and 1-2 %, respectively. According to this, we assume the gross reservoir porosity to be 0.01 ≤ φ ≤ 0.1. We assume that the average reservoir pressure is approximately 2 MPa (see Barker et al. 1992). For Poisson's ratio of the reservoir we use a range of 0.1 ≤ ν ≤ 0.4.
The total extracted mass for the investigated time interval is 1.4 × 10 12 kg and the reservoir volume change is V r ≈ 9 × 10 7 m 3 (Mossop & Segall 1997), giving a ratio for volume change per extracted mass of V r / M f, t of 6 × 10 −5 m 3 kg −1 .
First we test the case of two-phase compressibility. As for the previous examples, we can use eq. (11) to calculate K f , which for a temperature of 212 • C and a pressure of 2 MPa (saturation conditions, see Table 1) yields values of K f < 0.1 MPa. Such low values of K f produce low values of V r / M f, t , considerably lower than 6 × 10 −5 m 3 kg −1 . Only in combination with very low values of K < 0.1 GPa such high values of 6 × 10 −5 m 3 kg −1 can be reached, which may be unlikely.
Next, we investigate single-phase conditions. For the pure water case, we assume a temperature of 210 • C and pressure of 2 MPa. According to these conditions we assume a density of water of 850 kg m −3 and a fluid bulk modulus of K f = 1000 MPa. We can apply eq. (10) to determine the most likely combinations of χ and K (see Fig. 7a). The results imply that under the assumption that we observe deformation from the whole reservoir, the case of pure water as pore fluid requires a minimum recharge (1 − χ ) of 90 %. For the pure vapour case, we assume a temperature of 240 • C and a pressure of 2 MPa. For these conditions, we use a bulk modulus of vapour of K f = 2 MPa and a density of 10 kg m −3 . The results allow a wider range of χ than the pure water case, with a minimum recharge of 50 % for K < 5 GPa (see Fig. 7b).
It seems to be generally assumed that reservoir conditions at The Geysers are vapour-dominated two-phase (Williamson 1992;Allis & Shook 1999;Mossop & Segall 1999). Our results indicate that this can only serve as an explanation for the magnitude of deformation per mass of fluid extracted if the bulk modulus is very low, with K < 0.5 GPa. Pure steam compressibility on the other hand can plausibly explain the observed deformation for a wider range of K, see Fig. 7(b). Considering the long-time interval being studied it is possible that the reservoir evolved through different stages during the time of observation, which would make the estimated amount of reservoir volume change an average of two-phase and steam conditions.

D I S C U S S I O N
In a fluid reservoir, the relation between the volume or mass of extracted fluid to the volume change of the reservoir due to contraction depends on elastic and poroelastic parameters of fluid and solid (K, K f , ν, α, φ) as well as the amount of recharge. The equations derived in this study show that the reservoir volume change  (Mossop & Segall 1997, both for pure liquid (a) and pure vapour phase pore fluid (b). The dashed black line shows K = 4.6 GPa, the estimated maximum by (Mossop & Segall 1999). Colour scale depicts the relative probability of parameter combinations. and the fluid volume change cannot be expected to have the same value, except for the idealized case of both the pore fluid and solid matrix being incompressible. Expected values should be V r / V f, t < 1, or V r / M f, t < 10 −3 m 3 kg −1 in case of fluid mass, where the magnitude primarily depends on the ratio of the reservoir bulk modulus and fluid bulk modulus and the amount of recharge into the reservoir.

Difference in V r / M f, t between different geothermal reservoirs
Three examples of geothermal systems in Iceland and the USA, all of which exhibit fluid pressure decrease, show different magnitudes of deformation despite all of them being fractured rock two-phase systems. The Hellisheidi system exhibits considerably lower values of deformation per extraction ( V r / M f, t of ∼6 × 10 −6 m 3 kg −1 ) than the Reykjanes system and The Geysers ( V r / M f, t of ∼4-6 × 10 −5 m 3 kg −1 ). Investigating the different examples indicate that the difference of this value is strongly affected by the compressibility of the pore fluid, the value of which being controlled by whether the fluid is in a single-phase or of two-phase state. The two-phase compressibility of a water-steam system can be very high, higher than either of the single-phase counterparts (see Section 2.1.1). This means that in a two-phase system the pressure change, as well as the deformation, will be significantly smaller even if extraction rate, reservoir volume and reservoir properties are the same. Considering how geothermal reservoirs evolve over time of exploitation and what the circumstances at the respective geothermal fields are, fluid conditions may play an important role in controlling the observed values of V r / M f, t . Due to production-induced pressure decrease, high-temperature geothermal reservoirs can go through different stages during their lifetime, from liquid dominated to two-phase to vapour-dominated. This way, order of magnitude differences in V r / M f, t between two reservoirs could be explained by different fluid compressibilities, which depends on reservoir pressure and temperature, and production-induced changes in reservoir pressure and/or temperature could cause large changes in V r / M f, t over the time for a single reservoir.
At Hellisheidi we find that the relatively low V r / M f, t can very well be explained by two-phase conditions in the reservoir. The results agree well with the estimate for the bulk modulus K (1-6 GPa) that we obtained by the inversion of geodetic data (Section 3.1.1) and yield realistic values of the liquid volume fraction f l . At Reykjanes and The Geysers, single-phase conditions seem to allow a wider range of parameter values, while two-phase conditions can only serve as explanations with very low bulk modulus values, see Figs 5-7. For two-phase conditions at Reykjanes, values of K of above 1 GPa can only be obtained for very low values of f l , at The Geysers K needs to be below 0.1 GPa for two-phase fluid compressibility.
The Reykjanes power plant had just been commissioned at the time when Keiding et al. (2010) studied the deformation (2005)(2006)(2007)(2008). Hence, it is likely that the reservoir was still dominated by liquid during this time period, explaining the high ratio of deformation-to-mass extraction (∼4 × 10 −5 m 3 kg −1 ) we estimate from their study. A recently published study on Reykjanes by Parks et al. (2018) found that for the time interval 2009-2016 the total deformation for their preferred model is 1 × 10 6 m 3 . This means that for a total extraction mass of approximately 1.1 × 10 11 kg (Parks et al. 2018) during the same time period, the deformation per mass of extracted fluid had decreased to ∼9 × 10 −6 m 3 kg −1 (close to that at Hellisheidi of ∼6 × 10 −6 m 3 kg −1 ). The decrease in this ratio might indicate the transition from an initial liquid state of the fluid to two-phase. This effect would also be observable as a stabilization of pressure in the reservoir (see Grant 2011, chapter 3.4.2), which seems to occur at Reykjanes from approximately 2007 onwards . Another possibility is that the natural recharge increased due to the decrease in reservoir pressure, which we have not yet tested for the more recent ratio of V r / M f, t at Reykjanes. At The Geysers, the production had been going on for several years before the time interval of 1977-1996 that was studied by Mossop & Segall (1997. In such a long-time interval it is possible that a change in prevailing fluid phase has occurred. Mossop & Segall (1999) break down the long-time interval into shorter ones and find that the yearly volume changes vary between the shorter time intervals. In addition, they report that volume strain scales differently with pressure in different regions within the reservoir. This way both Reykjanes and The Geysers may experience temporal and spatial variations of reservoir (and, in particular, fluid) properties that are being averaged in the volume change estimate that can lead to ambiguity in the parameter estimates. Another factor that could affect the results is the model for two-phase compressibility, which assumes thermal equilibrium between the phases. In a reservoir where water and steam are not in thermal equilibrium, this model may have to be adjusted.
For the examples discussed the parameters are only very little, if at all, constrained by a priori information, which makes it difficult to reach conclusive results. There may be cases, however, where the elastic properties or the state of reservoir fluid is better known in which case the methodology of this paper can lead to better constrained results. Nevertheless, our study indicates that even if reservoir properties cannot be conclusively determined, the ratio V r / M f, t can be an indicator for fluid conditions in the reservoir. The large differences in compressibility between singleand two-phase fluids may explain order-of-magnitude differences in V r / M f, t . High values of this ratio can be explained by singlephase compressibility, low values by two-phase compressibility. If the state of the fluid is known a priori, on the other hand, eq. (10) or eq. (14) can help to constrain the bulk modulus of the reservoir rock mass, reservoir recharge or other model parameters.

Deformation caused by reservoir temperature changes
While this study focuses on deformation caused by reservoir pressure changes, it should be noted that changes in temperatures can contribute to the deformation field, as well. The three geothermal areas we investigate appear to be dominated by pressure-driven deformation (as suggested by Mossop & Segall 1997;Keiding et al. 2010;Juncu et al. 2017), which is why we did not consider the effect of thermoelastic deformation in this study. If there was significant thermoelastic deformation, however, it should be included in the modelling because it would have an effect on parameter estimates (e.g. for the bulk modulus).
An example of deforming geothermal reservoirs due to temperature decrease may be the Krafla geothermal field in North Iceland. Drouin et al. (2017) report that pressure levels have been constant but temperatures have been decreasing, and based on their modelling they reason that the temperature decrease caused by boiling of incoming recharge water is the driving cause of deformation. Interestingly, Drouin et al. (2017) find a value of V r / M f, t very similar to the Hellisheidi example, which in their case is ∼8 × 10 −6 m 3 kg −1 . For the thermoelastic case, the deformation behaviour could be studied similarly to the poroelastic analysis presented in this study.

Low values of K in geothermal areas
The values of the reservoir bulk modulus K that we estimate for Hellisheidi using geodetic data, pressure changes and extraction rates, are quite low when compared to values typically obtained from laboratory experiments or seismic data (see e.g. Turcotte & Schubert 2002;Pimienta et al. 2017). There are a number of possible explanations for this discrepancy that can be found in the literature: seismic wavelengths greatly exceed the pore and fracture length scales, scales of laboratory experiments being much below field scale (e.g. Heuze 1980;Ciccotti & Mulargia 2004) and rock characteristics that are specific to the site's location. Those characteristics can be rock type, degree of fracturing, confining pressures, degree and type of hydrothermal alterations, and temperature (e.g. Zimmerman et al. 1986;Lynne et al. 2013;Frolova et al. 2014;Hassanzadegan et al. 2014). Of those, in particular hydrothermal alterations can have a significant effect on rock strength and rigidity, as shown, for example, at hydrothermal fields in New Zealand  and Russia (Frolova et al. 2014). Lynne et al. (2013) reported that for samples from the Tauhara geothermal field in New Zealand that values of the constrained modulus (uniaxial bulk modulus) could vary from below 100 MPa to over 1 GPa depending on whether clay minerals are present as a result of hydrothermal alteration. The Hellisheidi field exhibits a large range of hydrothermal alterations (Helgadóttir et al. 2010) and their possibly degrading effect on the reservoir rock could explain low bulk moduli.

C O N C L U S I O N S
We derive an equation (eq. 10) that relates rock volume changes estimated from deformation studies to the mass of extracted fluid from a flat reservoir, for a set of common simplifying assumptions. The equation shows that the change in reservoir volume, given the mass of extracted fluid, depends on fluid compressibility and elastic rock properties, as well as the rate of fluid recharge. We use examples from three different geothermal fields, where the ratio of deformation volume per mass of extracted fluid varies by one order of magnitude, to investigate how the obtained relation between deformation, extraction and material properties can be applied. Using a stochastic approach we identify how varying parameter combinations affect the difference in the volume change per fluid mass extraction ratio, and find that the compressibilities of pore fluid and reservoir rock have a significant effect. For the Hellisheidi field in Iceland, for example, a combination of two-phase compressibility and low reservoir bulk modulus, K 10 GPa, agree well with the observed deformation per fluid mass extracted, and with the deformation-based estimate of K of 1-6 GPa. For Reykjanes and The Geysers, results are less conclusive but indicate that deformation per extracted fluid mass may be best explained by pure liquid compressibility at Reykjanes, and by pure steam compressibility at The Geysers. As a general result, our calculations indicate that large, order-of-magnitude differences in the ratio of reservoir volume change per mass of extracted fluid (either between different geothermal systems, or for changes over the time in a single reservoir) can be well explained by differences in fluid compressibility, which is determined by the reservoir fluid conditions (two-phase, liquid or steam). While the fluid compressibility may not be the only possible cause of such differences, we argue that it should be considered a plausible candidate that can be tested employing the equations and methods we used in this study.

A C K N O W L E D G E M E N T S
We would like to thank Andrew Hooper (University of Leeds), Björn Lund (Uppsala University), Sigurjón Jónsson (King Abdullah University of Science and Technology),Ólafur Flóvenz (Iceland Geo-Survey) and Freysteinn Sigmundsson (University of Iceland) for their helpful comments on the manuscript. Comments from the editor (Michael Ritzwoller) and two anonymous reviewers were highly appreciated and helped to improve this paper. COMET is the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics. This work was supported in part by grants from the Iceland Research Fund (grant numbers 130371-051/052/053 and

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Bivariate histograms of the source parameters for the Hellisheidi geothermal field in the model optimization using geodetic data. Bulk modulus K with Poisson's ratio ν (a), Biot's coefficient α (b), depth (c), radius (d), thickness (e) and pressure decrease (f) of the reservoir. The colour scale gives the frequency of results within a bin, from higher frequencies depicted by yellow to lower frequencies depicted by blue. Figure S2. Bivariate histograms of the source parameters for the Hellisheidi geothermal field in the model optimization using geodetic data. Reservoir volume change V r with Poisson's ratio ν (a), Biot's coefficient α (b), depth (c), radius (d), thickness (e), pressure decrease (f) and bulk modulus K (g) of the reservoir. The colour scale gives the frequency of results within a bin, from higher frequencies depicted by yellow to lower frequencies depicted by blue. Figure S3. Histogram of the estimated posterior probability of V r at Hellisheidi as black bars, and normal distribution with standard deviation of 30 % of the mean of 2.1 × 10 5 m 3 as blue line. Table S1. Estimated parameters of the deformation sources and their 90 % CI from joint inversion of GNSS and InSAR data. The depths and the coordinates of the sources refer to the centre of the bodies. For the geothermal sources, the semi-major axis is the half the long axis of the elliptic cross-section, orientated along its strike direction that is measured from north.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.