A Simple Relation to Constrain Groundwater Models Using Surface Deformation.

A simple relation between pore pressure change and one-dimensional surface deformation is presented. The relation is for pore pressure change in a confined aquifer that causes surface deformation. It can be applied to groundwater models of any discretization and is computationally efficient. The estimated surface deformation from model results can be compared to observed surface deformation through geodetic techniques such as Differential Interferometric Synthetic Aperture Radar. Model parameters then are constrained using the observed surface deformation. The validity of this relation is shown through constraint of model parameters for surface uplift due to pore pressure increase caused by wastewater disposal injection.

Introduction

Constraining hydrologic parameters in groundwater modeling is crucial to model results that are representative of the actual system. Obtaining direct measurements of hydrologic parameters for groundwater models is often difficult. Deep confined aquifers, especially those used for wastewater disposal of produced water from oil and gas extraction, have limited accessibility and few direct measurements of hydrologic parameters are available. Groundwater models in locations of wastewater injection are important to determine the cause of induced seismicity and plan mitigation strategies. In these models, literature values of hydrogeologic parameters are frequently used for the specific aquifer formation or type of formation (e.g., Keranen et al. 2014; Brown and Liu 2016). Some injection‐induced seismicity areas may have core samples available for laboratory testing (e.g., Brown et al. 2017). In selected cases, induced seismicity sites may have well tests for the aquifer of interest that can be used to estimate hydrologic properties (e.g., Hornbach et al. 2015; Ogwari and Horton 2016; Brown et al. 2017).

Surface deformation can be a physical manifestation of pore pressure changes in the subsurface (Teatini et al. 2011). When pore pressure increases within a confined aquifer, the aquifer skeleton expands resulting in surface uplift. Surface deformation is observed using a number of techniques including leveling (e.g., Bell et al. 2002), global positioning system (GPS) (e.g., Ishitsuka et al. 2017), tilt‐meters (e.g. Jahr et al. 2008), and remote sensing like Differential Interferometric Synthetic Aperture Radar (DInSAR) (e.g., Shirzaei et al. 2016; Barba‐Sevilla et al. 2018). Uplift on the order of millimeters to centimeters has been observed using DInSAR in multiple areas (e.g., Shirzaei et al. 2016; Barba‐Sevilla et al. 2018; Kim and Lu 2018; Loesch and Sagan 2018) of wastewater disposal injection related to oil and gas activities. Interferograms are created by co‐registering two SAR images of similar imaging geometries. DInSAR relies on calculating phase differences from two or more passes of the SAR satellite to quantify surface displacement. DInSAR can provide both cumulative deformation and time series of deformation over a broad area with centimeter or better resolution.

In areas of wastewater disposal, oil and gas production in the vicinity must also be considered as a cause of surface deformation. However, if the production is mostly by enhanced oil recovery, where fluids are injected into the subsurface to promote the hydrocarbons flowing into production wells, then the surface deformation from production can be assumed to be negligible. A goal of the enhanced oil recovery is to keep the pressure in the production zone constant (Rubinstein and Mahani 2015). The amount of fluids being injected into the production zone is assumed to be almost equal to the amount being extracted and consequently there is little net fluid injection. Therefore, the surface deformation is caused primarily by the wastewater injection (Figure 1).

Figure 1

Conceptual diagram depicting surface uplift related to pore pressure increase. The enhanced oil recovery (EOR) within the producing zone does not cause surface deformation due to the net zero pore pressure change. The wastewater disposal injection (right well) does cause surface uplift due to the large increase of pore pressure and deformation in the confined disposal zone.

Surface deformation has been used to estimate hydrologic parameters for aquifer systems where in‐situ aquifer test data are scarce (Liu and Helm 2008; Hu et al. 2018). Rinaldi et al. (2017) used inverse modeling of surface uplift to determine mechanical and hydrologic properties at the CO2 injection site at In Salah, Algeria.

Here, we present a simple relation that relates pore pressure changes to one‐dimensional (1D) compaction or expansion, and therefore, surface deformation by hydrologic parameters. This work is motivated by the new applications of DInSAR in areas of wastewater injection (e.g., Shirzaei et al. 2016; Barba‐Sevilla et al. 2018; Kim and Lu 2018; Loesch and Sagan 2018). The relation can be used for constraining hydrologic parameters in groundwater models, numerical or analytical. The simplicity of the relation allows for broad accessibility to modelers of all levels without being computationally expensive.

Relation of Modeled Change in Hydraulic Head to Surface Deformation

Modeled change in hydraulic head in a three‐dimensional (3D) groundwater model can be related to surface deformation by the hydrologic parameters. Vertical deformation due to compaction or expansion of the geologic material in the groundwater system can be calculated, as explained by Hoffman et al. (2003), based on the Terzaghi (1925) principle of coupling between sediment compaction and changes in hydraulic head.

Hydraulic head (h) can be expressed by the sum of elevation head () and pressure head: where is pore pressure (M L−1T−2), (ML−3) is fluid density, and g (L T−2) is gravitational acceleration. Elevation head does not change with time at any particular location; therefore, the change in hydraulic head can be expressed by the change in pressure head or the pore pressure: .

In the most general form, the effective stress tensor () is calculated from: where (M L−1T−2) is total stress tensor. The Kronecker delta is defined as 1 if i = j and 0 if i ≠ j. Equation 2 becomes the following when only the vertical component (zz) is considered: Considering that the change in hydraulic head can be expressed by the change in pressure head (simplified Equation 1), the vertical total stress () remains constant through time, and changes in fluid density caused by the compression or expansion of water are negligible, the change in vertical effective stress can be expressed by substituting Equation 1 solved for pore pressure into Equation 3: where (L) is the change in hydraulic head.

We now relate effective stress change to head change through material compressibility that is defined as: where (L T2M−1) is the compressibility of the material, (L3) is the change in volume of a control volume with initial volume V (L3), and (M L−1T−2) is the change in effective stress. When considering only the vertical direction, Equation 5 becomes: where (L T2M−1) is the vertical compressibility of the material and (L) is the change in thickness of a control volume with initial thickness b (L). If assuming only changes in pore pressure are responsible for changes in effective stress, Equations 4 and 6 can be combined: This formulation is the basis of the USGS Subsidence and Aquifer‐System Compaction (referred to as SUB) package (Hoffman et al. 2003) in MODFLOW‐2005. The surface deformation from wastewater injection is often assumed to be caused by the expansion of the rocks due to increases in pore pressure (e.g., Teatini et al. 2011). We assume that pore pressure increase and deformation are contemporaneous, that is, delays due to fluid draining or pressure dissipation in this process are negligible. We also assume that the deformation is elastic.

For any 3D groundwater model of a confined aquifer, the total surface deformation can be calculated using a variation of Equation 7, integrating deformation over all layers affected by a pore pressure change. The compressibility of the material is taken into account through the specific storage, which is a measure of the ability of the aquifer to store or release water per unit volume per unit change in hydraulic head. Specific storage () (L−1) is a function of the compressibility of the aquifer material () and the compressibility of water (): where is porosity [1] of the aquifer material and β the compressibility of water (L T2M−1). Generally, the product of can be considered negligible compared with . We should note that there are certain lithologies where (Freeze and Cherry 1979); however, these are unlikely to be injection intervals in sedimentary rocks. Therefore, the change in thickness () (L) can be expressed by substituting Equation 8 into Equation 7: Using Equation 9, the total surface deformation can be calculated for a model of discretization index i, j, k (columns, rows, layers), by summing over the model layers between the injection layer and the land surface: where D (L), the total surface deformation as a function of model index location i, j (column, row) and time t, is the sum of the individual change in thickness of N model layers (index k). Each model layer's change in thickness is calculated using the layer's specific storage, initial thickness, and change in hydraulic head. Figure 2 presents a schematic of calculating the surface deformation based on the model discretization. The total surface deformation estimated from the groundwater model can then be compared to the observed surface deformation results calculated by DInSAR or other methods like GPS and leveling.

Figure 2

Schematic of calculating surface deformation from groundwater model results. Cumulative vertical deformation at the red location can be calculated by adding Δb using the thickness b , specific storage , and change in hydraulic head, of N = 3 layers in Equation 10.

Figure 2 shows a uniform groundwater model grid in vertical and lateral directions, but Equation 10 can be applied to models of any discretization. Once the estimates of surface deformation due to modeled hydraulic head (pore pressure) change are completed, a properly georeferenced map can be generated by interpolating the deformation estimates at the center points of each cell. This allows for comparison between the observed surface deformation and the modeled surface deformation.

Example Application of Technique

To demonstrate this method, we use the DInSAR results of Kim and Lu (2018) in Winkler, Texas. Kim and Lu (2018) calculated a vertical uplift between late 2014 and April 2017 of approximately 5 to 6 cm near two wastewater disposal wells (API No. 42‐495‐33675 and 42‐495‐30150). The InSAR data were collected from the European Space Agency's Sentinel‐1 A/B satellites and processed using the multi‐dimensional small baseline subset (MSBAS) method. The majority of the uplift (∼5.0 cm) occurred after the disposal well API No. 42‐495‐33675 started injection in January 2016. The maximum uplift occurred close to the disposal well API No. 42‐495‐33675. A time series of the uplift at approximately 270 m from the well was created by Kim and Lu (2018). The overall uplift was of an irregular elliptic shape that decreased to less than 1 cm approximately 2 km from the disposal well. To illustrate the applicability of the relation in constraining parameters in groundwater models, we constrain the specific storage of the aquifer using Equation 10. This example focuses on the one location of the DInSAR time series of Kim and Lu (2018); however, this technique has the advantage of providing data at numerous locations.

Parameter Estimation Demonstration Using Theis Solution

We estimate the change in hydraulic head caused by wastewater injection as a function of distance from the well, r (m), and time since injection began, t (days), using the analytical Theis (1935) solution: where Q is the injection rate (m3/day), T (m2/day) is transmissivity, is dimensionless time parameter, S [1] is storativity, and x is the variable of integration. Transmissivity of an aquifer is a measure of the ability of the aquifer to transmit fluid and is calculated by multiplying the hydraulic conductivity of an aquifer by the saturated thickness. The storativity (specific storage multiplied by the saturated thickness) of an aquifer is a measure of how much water can be stored or released for a unit volume of the aquifer. In confined aquifers, storativity is a function of the aquifer compressibility and the compressibility of water.

The Theis solution assumes a homogeneous, isotropic confined aquifer of infinite lateral extent. Therefore, this estimate does not take into consideration any pore pressure diffusion into the formations above or below the injection interval. Rather, deformation is calculated as only a result of pressure changes in the injection formation. Any deformation that may occur above the injection interval is assumed negligible because of the minimal pressure changes.

Kim and Lu (2018) indicated that the wastewater is injected into the 70 m thick Bell Canyon Formation sandstones; however, according to the H‐10 forms publicly available online from the Texas Railroad Commission, that is only true of well API No. 42‐495‐30150. The second well (API No. 42‐495‐33675) closest to the observed deformation centroid injects into 730 m of the Bell Canyon Formation and Bone Springs Formation. Kim and Lu (2018) state that the injection zone is confined above by very low permeability limestone, as illustrated schematically on Figure 1. The average permeability of the Bell Canyon Formation sandstones is 40 mD (3.95 × 10−14 m2) (Dutton et al. 2003) and the permeability of the Bone Springs Formation is approximately 2 mD (1.97 × 10−15 m2) (Montgomery 1997). We calculated the average horizontal permeability ( using the arithmetic average: where k 1 and k 2 are the permeabilities (mD) and b and b are the thicknesses (m) of the Bell Canyon Formation (70 m) and Bone Springs Formation (660 m) respectively. The average permeability of the injection interval is 5.64 mD (5.57 × 10−15 m2). This permeability converts to a hydraulic conductivity of 6.14 × 10−8 m/s; transmissivity is calculated by multiplying the hydraulic conductivity by the aquifer thickness for a value of 4.48 × 10−5 m2/s. We estimate a constant injection rate over the 486 days between January 2016 and April 2017, by using the average injection rate, 1265 m3/day, for the disposal well. The injection volumes are publically available (H‐10 forms) online through the Texas Railroad Commission. The change of hydraulic head is calculated for radii between 1 and 2000 m from the well, with a particular interest at 270 m from the well, the approximate location of the DInSAR time series. We start with an estimated initial storativity of 10−4 [1] (Figure 3).

Figure 3

Uplift estimates (cm) from the Theis solution. (a) Results after 486 days of constant injection for multiple storativity values. The distance from the well increases along the x‐axis. The approximate location of the DInSAR time series (Kim and Lu 2018, Figure 2a), 270 m from the injection well, is marked by the dashed black line. (b) Time series of uplift estimates at approximately 270 m from the injection well with a storativity of 3.5 × 10−4. The Theis solution uplift estimates are in the black line/circles, and the DInSAR time series data (Kim and Lu 2018) are indicated by gray triangles.

With the initial storativity estimate, the change in hydraulic head after 486 days of constant injection at radii of 270 and 2000 m from the well is approximately 165 and 63 m respectively. Using Equation 10, the estimated surface deformation is 1.65 and 0.63 cm (Figure 3). We then constrain the storativity of the aquifer by matching the measured surface deformation (∼5 cm since well 42‐495‐33675 began injection) to the estimated surface deformation using Equation 10. For approximately 5 cm of uplift vertical surface deformation at the time series location (approximately 270 m from the well), the storativity must be approximately 3 × 10−4 to 5 × 10−4 [1] (Figure 3). Figure 3b compares the DInSAR time series from Kim and Lu (2018) to the uplift estimated based on the head change from the Theis solution. The two datasets will not match perfectly because we used a constant injection rate in the Theis solution head calculation. However, the general trend in uplift is of the correct order of magnitude and follows the observed DInSAR uplift.

Parameter Estimation Demonstration Using Groundwater Flow Model

To further demonstrate the relation, a groundwater model of the above problem was created. We created a 3D groundwater flow model using MODFLOW‐2005 (Harbaugh 2005) with a domain of 10 km × 10 km × 730 m. The model simulated the injection interval which is the approximately 70 m thick Bell Canyon Formation overlying approximately 660 m of the Bone Springs Formation. The Bell Canyon Formation layer was assigned a hydraulic conductivity of 4.35 × 10−7 m/s (transmissivity of 3.04 × 10−5 m2/s) and the Bone Springs Formation a hydraulic conductivity of 2.18 × 10−8 m/s (transmissivity of 1.44 × 10−5 m2/s). The injection well was placed in the center of the model and injected uniformly over the entire injection interval. The injection rate followed the reported values from January 2016 through April 2017 on the H‐10 forms. General‐head boundaries were set on each side boundary to virtually extend the size of the model domain. General‐head boundaries allow for head dependent flux across the boundary by setting a hydraulic head at some distance beyond the boundary. The General‐head boundaries specify a hydraulic head of zero 10 km beyond the model domain boundaries. The top and bottom of the model were assigned no‐flow boundaries. Note we did not include any heterogeneity or structures, which are likely present in the subsurface. An initial storativity of 3.5 × 10−4 [1] was used based on the average injection rate and Theis solution estimates, which led to a total uplift at 270 m from the injection well of approximately 4.20 cm (Figure 4). Specific storage and hydraulic conductivity of the two formations was adjusted using a systematic ad‐hoc approach within the well‐defined range of K for the formations and an appropriate range of S to match the head change caused surface deformation to the DInSAR surface deformation time series near the injection well (Figure 4b). Future users of this methodology would benefit from an incorporation of an optimization algorithm for parameter estimation.

Figure 4

Uplift estimates from groundwater model. (a) Map view of uplift (cm) estimated from modeled change in hydraulic head after 486 days of injection using the reported injection rates with a storativity (S) of 3.5× 10−4 [1]. The approximate location of the DInSAR time series (Kim and Lu 2018, in Figure 2a), 270 m from the injection well, is shown at the black dot. (b) Time series of uplift estimates at approximately 270 m from the injection well (black dot in part (a)). The modeled uplift (cm) estimates are in the black line/circles, and the DInSAR time series data (Kim and Lu 2018) are indicated by gray triangles. (c) Model parameters for each modeled uplift result.

With the variable injection rate, storativity values between approximately 3 × 10−3 and 3 × 10−4 [1] produced surface deformation that is within the observed values. Figure 4b compares the time series of estimated uplift from the groundwater model to the DInSAR time series (Kim and Lu 2018). The overall pattern matches well and is of the same order of magnitude; however, the observed DInSAR uplift has a more gradual uplift than the uplift estimated by the groundwater model results. The discrepancies between the modeled uplift and the DInSAR time series are likely largely due to the lack of heterogeneity in the model and some of the assumptions for this method as discussed in the following section. The DInSAR data from Kim and Lu (2018) shows an elongation of the uplift to the south of the wells which indicates heterogeneity and/or anisotropy. In particular, there are two distinct elongations approximately 60° apart which may indicate subsurface preferential pathways possibly due to faulting or fractures. This was not modeled for the demonstration. This method assumes instantaneous uplift with the pore pressure change without loss in magnitude between the disposal zone and the surface and no lateral deformation.

For the Theis comparison, we only varied one parameter, specific storage, to constrain the model. Varying multiple parameters, as with the groundwater model demonstration, to optimize the match is best when using surface deformation data. Specific storage of approximately 10−7 to 10−6 m−1 (storativity ∼10−4 to 10−3) translates to an aquifer compressibility of approximately 1 × 10−11 to 1 × 10−10 Pa−1 using Equation 8. This aquifer compressibility is consistent with competent bedrock (Freeze and Cherry 1979).

Discussion and Conclusions

Relating 1D compaction to pore pressure changes is not a new concept (e.g., Terzaghi 1925). However, the concept has not been widely applied to use surface deformation as a constraint on models of pore pressure change from injection. Complicated inverse modeling has been conducted to constrain mechanical and hydrologic parameters (e.g., Rinaldi et al. 2017). While advanced inverse modeling approaches achieve well constrained models with uncertainty estimates, they are computationally expensive and require a high level of modeling experience and time. The simple physics based relation we present here allows for first‐order constraints on model parameters regardless of model discretization and with low computational needs. Deformation caused by pore pressure change can be estimated using simple analytical solutions, or, if numerical models are used to obtain pore pressure in three dimensions, the deformation from each model layer can be integrated vertically to obtain the surface uplift. In areas with little constraint on hydraulic parameters, like areas of wastewater disposal in deep basal aquifers, this relation could be extremely useful. Models of any type of injection into confined aquifers can be constrained using Equation 10. This approach can also be applied to study sites where artificial recharge is being used or proposed to control subsidence (e.g., Zhang et al. 2015) or store water for later use (e.g., Farid et al. 2018). In addition, other natural processes where fluid pressure causes a surface deformation could benefit from use of this relation, for example, to constrain changes in permeability and rock properties where surface uplift are observed following earthquakes or around natural circulation of hydrothermal fluids around large calderas (Hurwitz et al. 2007).

We noted that this method has limitations. Surface deformation is the expression of possibly multiple processes occurring in the subsurface. When using surface deformation to constrain a model the following certain assumptions must be made. If other information about the deformation is not available, one major assumption is that the surface uplift is caused entirely by pore pressure changes. In addition, this relation assumes that the deformation is elastic and instantaneous. There are examples of pore pressure induced deformation that are delayed from days to years (Ireland et al. 1984; Liu and Helm 2008; Chen et al. 2017). This often involves clay lenses within the aquifer which are inconsistent with the deep permeable aquifers (often sandstone) that are used for wastewater disposal. Since the relation is based on 1D‐compaction, deformation is assumed to be predominantly vertical. These assumptions can lead to a less than ideal fit between model estimated uplift and remote sensing inferred uplift. However, that non‐ideal fit can also hold information about the system, including evidence for heterogeneity and indication that other processes are also affecting uplift.

If other processes besides injection are causing surface deformation, the deformation recorded by DInSAR is either over or under estimating the deformation caused by the injection. If the other processes are causing subsidence or uplift, DInSAR underestimates deformation due to injection or DInSAR overestimates deformation due to injection, respectively. If a significant component of deformation is not elastic, then over time the DInSAR inferred deformation will underestimate the uplift due to injection as any compaction that occurs will not be recovered. Since uplift is not instantaneous, there will likely be a slight shift between modeled results and DInSAR deformation. Finally, if there is any lateral deformation, then the vertical DInSAR uplift will be underestimating the deformation caused by injection.

Use of DInSAR data also has its limitations. DInSAR does not always have high temporal resolution and often has problems with decorrelation (when pixels in the SAR images do not match between images) due to vegetation, snow cover, shifting sediments, and more (Massonnet and Feigl 1998). It is important to note that SAR satellites are side‐looking that is they do not look perpendicular to the ground but at some angle off the perpendicular (look angle). Displacement determined by DInSAR is in the line of sight, i.e., in the direction the satellite is looking when the SAR images were acquired. Vertical deformation can only be resolved with knowledge about the lateral motion at the site and/or ascending and descending orbit SAR images. If there is lateral motion captured by DInSAR, then the lateral deformation would need to be removed and only the vertical component included when using Equation 10. In addition, the corrections (e.g., atmospheric, topographic, etc.) that are used in DInSAR processing are based on models or limited data, and therefore may introduce uncertainty into the deformation results (Massonnet and Feigl 1998). There are other methods of measuring that can acquire surface deformation. GPS has a high temporal resolution but would need to be considered a point measurement for this method. Several or more GPS locations would need to be available to have reasonable constraints on a groundwater model. Using GPS to ground‐truth the surface deformation observed via DInSAR, although beyond the scope of this study, would be a further step to take for using this relation.

While there are limitations to the use of this relation, Equation 10 offers a simple formula for estimating surface uplift from groundwater models. This allows a first‐order constraint on the hydrologic parameters of the confined aquifer, especially in cases where parameter data are limited and injection is occurring.

Authors' Note

The authors do not have any conflicts of interest or financial disclosures to report.