Estimate of a spatially variable reservoir compressibility by assimilation of ground surface displacement data

Fluid extraction from producing hydrocarbon reservoirs can cause anthropogenic land subsidence. In this work, a 3-D finite-element (FE) geomechanical model is used to predict the land surface displacements above a gas field where displacement observations are available. An ensemble-based data assimilation (DA) algorithm is implemented that incorporates these observations into the response of the FE geomechanical model, thus reducing the uncertainty on the geomechanical parameters of the sedimentary basin embedding the reservoir. The calibration focuses on the uniaxial vertical compressibility cM, which is often the geomechanical parameter to which the model response is most sensitive. The partition of the reservoir into blocks delimited by faults motivates the assumption of a heterogeneous spatial distribution of cM within the reservoir. A preliminary synthetic test case is here used to evaluate the effectiveness of the DA algorithm in reducing the parameter uncertainty associated with a heterogeneous cM distribution. A significant improvement in matching the observed data is obtained with respect to the case in which a homogeneous cM is hypothesized. These preliminary results are quite encouraging and call for the application of the procedure to real gas fields.


Introduction
Fluid extraction from aquifer systems and hydrocarbon reservoirs are among the most frequent causes of anthropogenic land subsidence.In the framework of a sustainable development of energy resources, the availability of numerical models able to reproduce the monitoring data and to predict the future development of the land settlement is nowadays of paramount importance.In this study, an ensemble-based DA method is used to infer the geomechanical parameters characterizing the rock formation of a deep gas reservoir, thus reducing the prior uncertainties of the geomechanical model response.The DA framework essentially requires three main ingredients: (i) a model to simulate the physical process of interest, (ii) a set of observation data and (iii) a suitable algorithm to incorporate these data into the model response.In this work, a 3-D finite-element (FE) geomechanical model is used to predict the land surface displacements above a gas field where displacement observations are available.The cal-ibration focuses on the uniaxial vertical compressibility c M , which mostly influences the occurrence of land subsidence.Partitioning the reservoir into blocks by faults provides a basis for assuming a homogeneous c M within each block, and a heterogeneous c M from one block to another.The effectiveness of the DA algorithm to reduce the parameter uncertainty associated with the block-to-block c M heterogeneity is evaluated for a synthetic test case.Significant improvements are obtained with respect to the assumption of a homogeneous c M field throughout the model domain.In this work, an Ensemble Smoother (ES) is the DA algorithm used to estimate the compressibility by inversion of land surface displacement data.The ES is deemed adequate for this purpose because its implementation (i) avoids the simulation restart necessary to apply other more common data assimilation techniques, such as the Ensemble Kalman Filter (EnKF); (ii) significantly reduces the overall computational cost required by geomechanical model and the inversion procedure.In this paper, the description of the geomechanical model and its implementation Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.are presented is Sect.2, whereas Sect. 3 reviews the basic concepts of the ES approach.In Sects. 4 and 5, the generation of the prior uncertain parameters and the synthetic case results are presented and discussed with both a homogeneous and a heterogeneous c M .

Model description
The subsurface deformation results from the pore pressure change in space and time due to fluid injection into, or extraction from, deep reservoirs.In order to simulate the deformation up to the ground surface, we need to solve both the governing flow and the structural partial differential equations (PDEs).Geomechanical simulations are performed using a finite-element (FE) poro-elasto-plastic numerical model (e.g.Gambolati et al., 2001).In this study, an isotropic stress-strain constitutive law is used with the vertical uniaxial compressibility c M that depends on the stress state according to the hypo-plastic hysteretic model developed by Baù et al. (2002) andFerronato et al. (2013).The uncertain compressibility c M is calibrated by introducing a spatially variable multiplicative factor f CM , which allows scaling c M values in the regions where fluid pressure changes occur.Poisson's coefficient ν is assumed to be known and equal to 0.3.

Model setup
A one-way coupling approach is followed with the flow model first run and the outcome subsequently used as input data for the geomechanical model.The geomechanical FE grid comprises 320 901 nodes and 1 824 768 tetrahedral elements.The model domain covers an area of 52 km × 49 km and extends to about 5 km depth.Zero-displacements conditions are prescribed on the lateral and bottom boundaries.The top of the domain is a traction-free boundary.The simulations span one year, during which the reservoir experiences a fluid extraction.The pressure data are synthetic.Figure 1 shows the reservoir domain and the grid used for the geomechanical simulations.

Ensemble smoother
The ES algorithm consists of Monte Carlo stochastic simulations based extension to nonlinear problems of the classic Kalman Filter (KF) (Kalman, 1960).The ES algorithm follows a two-step forecast-update process.The forecast step involves simulating an ensemble of model states X based upon the solution of the geomechanical FE model , which depends upon uncertain parameters P and forcing terms q (e.g.pore-pressure): X = (P, q).
(1) In these simulations, each model state X is represented by land surface displacements from a subsidence distribution map.Each realization of the ensemble is run forward in time using random sets of the uncertain geomechanical parameters P, thus creating an ensemble X f of model states X.The model results at any given location in the simulated domain are spread over a range of values, representing the uncertainty in the surface displacement prediction.In the update step, the set of measurements z collected to-date, i.e. point measurements of land displacement at a number of locations, is perturbed to account for measurement errors and assimilated into the forecast system state X f to produce the updated state ensemble: Matrix d includes the perturbed measurements, and the matrix H contains binary constants (0 or 1) that map model results at measurement locations.The matrix d − H • X f incorporates the residual at these locations between the measured and the predicted data.The matrix K is called the "Kalman Gain" matrix (Kalman, 1960), and has the following structure: where C f is the forecast error covariance matrix associated with the model forecast X f , and R is the measurement error covariance matrix associated with the perturbed measurements d.The matrix K plays the dual role of: (a) spreading information from measurement locations to adjacent areas, allowing for the measurements to correct the predicted values throughout the model domain; and (b) acting as a weight that scales the correction terms according to model and measurement errors.As R approaches zero, which means lowerror measurements, the influence of K increases and the residual is weighted more heavily, so that the model forecast approaches the measurements.In contrast, as C f approaches zero, which indicates a relative agreement among model realizations, the influence of K decreases, and the residual is weighted less heavily, so that the model forecast receives little or no correction from measurements.Within the ES algorithm, any variable incorporated into the system state matrix X f can be corrected by assimilating measurement data if a spatial correlation exists between these variables and the data (van Leeuwen, 2001).In Eq. ( 1), since the geomechanical parameters P dictate the behavior of the model response , all uncertain parameters used in the forecast step can be included into the system state matrix X f , and conditioned by land surface movements in Eq. ( 3).This conditioning may provide updates for the geomechanical parameters P that should approach those of the rock formation.Compared to other techniques for characterizing subsurface systems, the ES algorithm is quite attractive because of its low computational burden and the ability to run entirely independent of the simulation model (Bailey and Baù, 2010a, b).

Prior distribution of uncertain parameters
In the present study, c M is assumed to be the only uncertain geomechanical parameter.Because the ES relies on a Monte Carlo approach, a prior probability distribution function (PDF) is needed to sample the prior ensemble of the multiplicative calibration factor f CM .In this section, the generation of the prior PDFs in two test cases is described: (1) f CM is uniform within the reservoir, and (2) f CM is spatially distributed.In the latter, the heterogeneity on f CM occurs only in the area shown in the enlargement of Fig. 2 where the pressure has changed due to fluid extraction.

Homogeneous compressibility (test case 1)
The calibration factor f CM is assumed to be spatially uniform within the area of Fig. 2. The values of the prior f CM ensemble are randomly sampled between 1 and 10 from a uniform PDF: The selected range is based on the outcome of a sensitivity analysis (not shown here) carried out to investigate the possible interval of the f CM variation.Figure 3 shows the spatial distribution of the mean and the standard deviation (σ z ) of the vertical displacements (u z ) from the forecast ensemble obtained by performing 100 Monte Carlo geomechanical simulations using the prior f CM ensemble.

Heterogeneous compressibility (test case 2)
In test case 2, the f CM is spatially distributed within the same area used in test case 1.This area is (14 × 10) km 2 wide and subdivided into 140 square cells.Each cell is characterized by a different f CM .A categorical indicator algorithm, that creates random realizations of a heterogeneous f CM field according to a given covariance model, is used.The f CM values are drawn from a discrete uniform distribution with the prescribed ten categories ranging from 1 to 10  on the distance between these cells and a prescribed correlation length, λ.The λ value has a direct influence on the degree of heterogeneity assigned to the spatial distribution of f CM .Preliminary simulations have suggested the choice of λ = 4000 m. Figure 4 shows one of the 100 realizations of the generated prior ensemble of the f CM field.Obviously, the mean over the ensemble in each grid block is equal to 5. As for test 1, the mean and the standard deviation of the u z forecast ensemble are shown in Fig. 5.

Synthetic land subsidence data
Ground-surface displacements data are used to infer the model state and the geomechanical parameters.The observations are collected from the land subsidence map shown in Fig. 6, obtained from a geomechanical reference simulation with a prescribed and "known" compressibility distribution.The f CM field is assigned on the basis of a plausible reservoir partition derived from the presence of faults and thrusts (Fig. 6).The synthetic data locations are uniformly distributed over the reservoir (Fig. 6) in the area with σ z greater than 0.001 (see Figs. 3b and 5b).This is to account for the reservoir behavior where pressure variation occurs causing land surface deformation.

Update of f CM
In test case 1, a significant reduction of the uncertainty associated with the prior f CM distribution is evident by comparing the prior and the posterior cumulative distribution functions (CDFs) plotted in Fig. 7. Defining the ensemble spread as the average absolute difference between the ensemble members and the ensemble mean, the Average Ensemble Spread (AES) for n MC Monte Carlo realizations is calculated as:

AES
where x i and x are the f CM value of the ith ensemble member and the ensemble mean, respectively.The prior AES, i.e. before the assimilation of data, equals 2.41 and reduces by about 93 % after assimilation.Therefore, the posterior CDF is highly constrained compared to the prior CDF with the most probable estimate value for f CM equal to 1.89, corresponding to the mean, or expected value, of the updated ensemble.
In test case 2, the ES performance is evaluated by calculating the AES index over each grid block of Fig. 2. The prior AES index ranges between 2.0 and 2.9 over the domain.After assimilation, the spread of the updated ensemble significantly reduces over the area where data points are collected, while higher AES values are found in the surrounding area (see Fig. 8b).Although data are collected only over a limited portion of the domain, a sensitivity analysis (not shown here) reveals that collecting data over the entire domain does not improve the assimilation outcome.Indeed, these observa-Proc.IAHS, 372, 351-356, 2015 proc-iahs.net/372/351/2015/tions cannot yield enough information to infer the model parameters because the deep reservoir experiences small pressure variation with a negligible influence on the land surface deformation of outer regions.The average AES over all grid blocks reduces by about 33 % after assimilation.Despite a lower relative reduction of AES compared to test 1, the inverse problem outcome is much better verified in terms of ground surface displacements (Sect.6.2). Figure 8a depicts the spatial distribution of the mean from the updated f CM field.

Forecast of ground surface displacement
The quality of the parameter estimation is validated by executing the posterior geomechanical simulations for both uniform and spatially variable c M .The f CM model input is the mean of the updated f CM ensemble from the outcome of test cases 1 and 2. The results of this simulation are compared with the synthetic land surface data of Fig. 6 using the Nor- where u zi , sim and u zi , obs are the simulated and observed land ground displacement at the ith assimilation location, respectively; u zobs,max and u zobs , min are the maximum and minimum observation values, respectively, and N is the total number of assimilation locations, i.e. 60 data points.Test case 1 and 2 gives NRMSE values equal to 15 and 3 %, respectively.Hence, the predicted vertical displacements provides lower data mismatch in test case 2 and the assumption of a heterogeneous compressibility field allows for a better fit of the synthetic observations (Fig. 9).The fitting of the model predictions to the reference data is very accurate with the heterogeneous c M , while a relevant underestimation is found in the homogeneous test case for the largest subsidence values, i.e. u z greater than 1 mm.

Conclusions
In this study, an ensemble-based DA approach, i.e. the ES, is used to infer the compressibility of the geomechanical model of a producing hydrocarbon reservoir by assimilating ground surface displacements.The methodology is applied herein to investigate two different conceptual models assuming (i) a homogeneous c M within the reservoir and (ii) a heterogeneous reservoir compressibility.The latter assumption is made to account for the reservoir compartmentalization due to the presence of a complex fault system.Thus, c M spatially varies within the reservoir and calls for the calibration of a heterogeneous c M field.A test case is used to assess the validity of the proposed methodology on a reservoir with a synthetic pressure variation.The assimilation of synthetic ground-surface displacements, i.e. simulated data obtained with a priori known c M distribution, provides satisfactory results.In particular, the heterogeneous c M field gives a much better fitting than the homogeneous c M .Hence, a more satisfactory application to a real case is expected by a DA approach with a spatially variable c M .

Figure 1 .
Figure 1.Axonometric view of the 3-D FE grid of the geomechanical model.The reservoir is embedded within the grid with different colors distinguishing the producing layers.The vertical exaggeration is 5.

Figure 2 .
Figure 2. 2-D view of the geomechanical FE grid of Fig. 1.The enlargement view within the red rectangle refers to the area affected by the pore pressure change.

Figure 3 .
Figure 3. Test case 1: (a) mean and (b) standard deviation (σ z ) from the forecast ensemble of the vertical displacement (u z ).

Figure 4 .
Figure 4. Test case 2: one realization out of the 2-D f CM field ensemble.The compressibility varies in the area corresponding to the enlargement of Fig. 2.

Figure 5 .
Figure 5. Test case 2: (a) mean and (b) standard deviation (σ z ) from the forecast ensemble of the vertical displacement (u z ).

Figure 6 .
Figure 6.Spatial distribution of the synthetic land subsidence data (u z ).Assimilation data are collected at the 60 measurement locations displayed in the map.The red lines represent the trace of the faults subdividing the reservoir into blocks with the prescribed reference f CM values.

Figure 7 .
Figure 7. Test case 1: prior and Posterior CDFs from the updated f CM ensemble.

Figure 8 .
Figure 8. Test case 2: spatial distribution of (a) the mean and (b) the performance index AES from the updated f CM ensemble after data assimilation.

Figure 9 .
Figure 9. Simulated vs. observed values of u z in test case 1 (homogeneous c M ) and test case 2 (heterogeneous c M ).