# A multimethod Global Sensitivity Analysis to aid the calibration of geomechanical models via time-lapse seismic data

A multimethod Global Sensitivity Analysis to aid the calibration of geomechanical models via... Summary Time-lapse seismic attributes are used extensively in the history matching of production simulator models. However, although proven to contain information regarding production induced stress change, it is typically only loosely (i.e. qualitatively) used to calibrate geomechanical models. In this study we conduct a multimethod Global Sensitivity Analysis (GSA) to assess the feasibility and aid the quantitative calibration of geomechanical models via near-offset time-lapse seismic data. Specifically, the calibration of mechanical properties of the overburden. Via the GSA, we analyse the near-offset overburden seismic traveltimes from over 4000 perturbations of a Finite Element (FE) geomechanical model of a typical High Pressure High Temperature (HPHT) reservoir in the North Sea. We find that, out of an initially large set of material properties, the near-offset overburden traveltimes are primarily affected by Young's modulus and the effective stress (i.e. Biot) coefficient. The unexpected significance of the Biot coefficient highlights the importance of modelling fluid flow and pore pressure outside of the reservoir. The FE model is complex and highly nonlinear. Multiple combinations of model parameters can yield equally possible model realizations. Consequently, numerical calibration via a large number of random model perturbations is unfeasible. However, the significant differences in traveltime results suggest that more sophisticated calibration methods could potentially be feasible for finding numerous suitable solutions. The results of the time-varying GSA demonstrate how acquiring multiple vintages of time-lapse seismic data can be advantageous. However, they also suggest that significant overburden near-offset seismic time-shifts, useful for model calibration, may take up to 3 yrs after the start of production to manifest. Due to the nonlinearity of the model behaviour, similar uncertainty in the reservoir mechanical properties appears to influence overburden traveltime to a much greater extent. Therefore, reservoir properties must be known to a suitable degree of accuracy before the calibration of the overburden can be considered. Geomechanics, Numerical modeling, Statistical Methods 1 INTRODUCTION A reduction in reservoir pore pressure, as a result of production, can lead to significant compaction. A decrease in reservoir thickness results in the straining of the surrounding rock mass and thus a change to the in situ stress state (e.g. Herwanger & Horne 2009; Fjær & Kristiansen 2009). Failure to anticipate and manage geomechanical issues associated with production can have a detrimental effect on the economic performance of a field. For example, wellbore failure due to fault reactivation in the overburden (e.g. Vudovich et al. 1989; Bruno 2002) and platform instability due to seabed subsidence (e.g. Kvendseth 1988) are all but some of the consequences reported in literature over the past few decades. Therefore, geomechanical modelling can be an essential reservoir monitoring tool. Understanding and anticipating the spatial changes to the in situ stress field can aid production and help maintain well and platform integrity. Estimating production induced stress change typically requires the coupling of a geomechanical model to a reservoir simulator (e.g. Minkoff et al. 2003; Samier et al. 2003; De Gennaro et al. 2010). The coupled model allows the modelling of pore pressure and saturation change within the reservoir (in the fluid domain) while additionally simulating changes in stress and strain (in the mechanical domain) to the reservoir and surrounding rock mass. Often the coupled model is constrained with the use of production data. However, time-lapse seismic data have been found to contain information about reservoir saturation (e.g. Landrø 2001; Calvert 2005; Gonzalez-Carballo et al. 2006; Alerini et al. 2014) and production induced stress change (e.g. Guilbot & Smith 2002). As a result, time-lapse seismic data have been used extensively to aid the calibration of both the reservoir (e.g. Staples et al. 2005) and geomechanical model (e.g. Hawkins et al. 2007; Staples et al. 2007; De Gennaro et al. 2008; Angus et al. 2015). However, unlike advanced calibration procedures applied to production simulation models, such as, advanced history matching (e.g. Gosselin et al. 2003; Emerick & Reynolds 2012), the geomechanical model is typically loosely calibrated, that is, both model and data show similar patterns of anomalies (e.g. Hatchel & Bourne 2005). By not utilizing numerical matching techniques, are we making the most out of the geomechanical information stored in time-lapse seismic data? To investigate, we conduct a multimethod Global Sensitivity Analysis (GSA) on over 4000 model perturbations of a Finite Element (FE) geomechanical model of a typical High Pressure High Temperature (HPHT) reservoir in the North Sea. The GSA results are used to investigate the behaviour of geomechanical models (i.e. evaluate the potential geomechanical information content of time-lapse seismic data), assess the feasibility of seismic calibration and aid future numerical matching studies. We focus our attention primarily to the overburden where the mechanical properties are often as unknown as those of the reservoir rock and considered to be equally as important (e.g. Vudovich et al. 1989; Bruno 2002). A suitable example being the uncertain mechanical properties of overburden chalks seen in North Sea reservoirs. However, overburden time-lapse seismic anomalies are not complicated by fluid effects (i.e. saturation change as seen within the reservoir). Therefore, they are assumed to be a purely mechanical consequence and thus a suitable focus for this study. Specifically, via the GSA, we examine the sensitivity of overburden seismic time-shifts to the various properties of a single unknown overburden layer in a geomechanical model. We attempt to screen model parameters with negligible influence (i.e. those with little influence on time-lapse seismic data), rank those that are most influential (i.e. those with the most influence over time-lapse seismic data and hence the main focus of calibration) and map the input parameter space (i.e. can the uncertainty in influential parameters be reduced by seismic data). 2 THE GEOMECHANICAL MODEL We use the FE modelling software ELFEN (Rockfield Software Ltd) to construct a geomechanical model typical of an HPHT reservoir in the North Sea. We also utilize ELFEN's single-phase production simulator to model hydromechanical coupling (described in greater detail in Appendix A). We simplify our model to 2-D and isotropic due to the large number of model perturbations that need to be computed for the GSA. Although ELFEN has the potential to model anisotropy, we assume an isotropic scenario to reduce the number of model parameters for this initial study. The model domain consists of a 20 × 9 km2 subsurface region, with a reservoir interval located at a depth of approximately 5 km. This depth is similar to that of HPHT fields in the North Sea. Multiple overburden layers and faults are included to increase model complexity. However, for this particular study the faults are not initialized and are left as simple lithological discontinuities. This helps reduce the number of model parameters analysed in the GSA. Fig. 1 shows the complete model geometry. The reservoir interval is overpressured to 110 MPa prior to production and depleted by roughly 50 per cent over a period of 20 yrs. The rate of depletion is shown in Fig. 2 and is based on the production profile of the HPHT reservoir given by Hawkins et al. (2007). Figure 1. View largeDownload slide Model geometry, with reservoir region highlighted in blue and faults in red. The chalk layer whose physical properties are deemed uncertain shaded in grey. Also shown are the three locations in which we estimate vertical traveltimes used in the Global Sensitivity Analysis (GSA). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 1. View largeDownload slide Model geometry, with reservoir region highlighted in blue and faults in red. The chalk layer whose physical properties are deemed uncertain shaded in grey. Also shown are the three locations in which we estimate vertical traveltimes used in the Global Sensitivity Analysis (GSA). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 2. View largeDownload slide Reservoir pore pressure reduction given in normalized units, 1 = 110 MPa. Figure 2. View largeDownload slide Reservoir pore pressure reduction given in normalized units, 1 = 110 MPa. The in situ stress and pore pressure state, prior to production, is determined through a period of settling time steps in which the model is loaded under gravity. Any disequilibrium caused by the loading on the geometry is given time to relax, avoiding numerical oscillations (i.e. unwanted dynamic effects) in the modelling results. The phreatic surface (pore pressure origin) is set 50 m above the surface of the model to mimic a shallow North Sea environment. The model mesh is unstructured with triangular elements approximately 70 m in size. Element size was determined such that the model produced stable results (free from unwanted dynamic effects), had acceptable computation time, and yielded suitable resolution (considering typical vertical and lateral resolution of seismic data between 1 and 5 km depth). Fig. 3 shows the CPU time for models with varying sized meshes. Note the exponential increase in CPU time per relatively small increments in mesh size. Figure 3. View largeDownload slide Model runtime versus mesh element size. Figure 3. View largeDownload slide Model runtime versus mesh element size. We design nine different homogeneous isotropic elastoplastic materials, one for each of the nine lithological layers of our model (Fig. 1). In keeping with North Sea geology, we assume a sandstone reservoir with an impermeable shale caprock and a thick layer of stiff chalk in the overburden. To simplify model behaviour, we substitute a complex salt underburden (typically seen in the North Sea) with a simple mechanically strong rock, typical of unfractured limestone/dolostone. To define an elastoplastic material in ELFEN requires the parametrization of a state boundary surface, nonlinear elastic and permeability relationships and a suite of consolidation properties. A detailed description of all these material property relationships and associated parameters can be found in Appendix B, whilst a definitive summary in Table 1. Material properties are determined using both the extensive literature available on North Sea geology and generic materials designed by Rockfield in their generic material database (e.g. Rockfield Software Limited 2012). Fig. 4 shows logs of the final geomechanical model properties, post-geostatic initialization and prior reservoir production. Using Young's Modulus, Poisson's ratio and bulk density, we estimate seismic P-wave velocities. Note that no static to dynamic conversion was used when calculating the P-wave velocity as the static mechanical properties give rise to credible dynamic velocities. Figure 4. View largeDownload slide Logs of Young's modulus, Poisson's ratio, bulk density and pore pressure through the final geomechanical model, post-geostatic initialization and prior to reservoir production (i.e. 0 yr). These logs are then used to calculate the P-wave Vp velocity profile. Layer boundaries are marked via dotted horizontal lines. The lithostatic and hydrostatic gradients are plotted in red on the pore pressure log. Reservoir layer is shaded. Note that no static to dynamic conversion was used when calculating Vp as the static mechanical properties give rise to credible velocities. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 4. View largeDownload slide Logs of Young's modulus, Poisson's ratio, bulk density and pore pressure through the final geomechanical model, post-geostatic initialization and prior to reservoir production (i.e. 0 yr). These logs are then used to calculate the P-wave Vp velocity profile. Layer boundaries are marked via dotted horizontal lines. The lithostatic and hydrostatic gradients are plotted in red on the pore pressure log. Reservoir layer is shaded. Note that no static to dynamic conversion was used when calculating Vp as the static mechanical properties give rise to credible velocities. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Table 1. Material properties for each layer in Fig. 1. The layer number represents increasing depth from the surface. Layer  Lithology  Grain &  SR3 model  Yield Surface Evolution  Elastic Properties  Permeability Properties    No.    fluid densities  (Equations B1-B3)  (Equations B4-B7)  (Equations B8 & B9)  (Equation B10)  Consolidation Properties      ρg  ρf  pc  pt  β  n  ψ  κ  λ  ϕref  ϕinit  Eref  A  B  n  υmin  υmax  m  c  K0  x  y  Kf  Ks  α  η  Kxy  Pp      (g cc−1)  (g cc−1)  (Mpa)  (Mpa)                (MPa)                (m2)      (MPa)  (MPa)    (MPa s−1)    (MPa)  1  Shale  2.69  1.02  0.17  −1  60  1.3  51  0.012  0.106  0.55  0.38  910  −0.2758  −0.2758  0.085  0.31  0.31  1  −0.9  4.20 × 10−19  3  2  2400  36000  1  1  0.6  0  2  Sandstone  2.65  1.02  0.28  −1  66  1.7  56  0.012  0.174  0.55  0.33  1850  −0.2758  −0.2758  0.15  0.31  0.31  1  −0.5  6.61 × 10−11  5  2  2400  26000  1  1  0.6  0  3  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.35  0.2  410  −0.2758  −0.2758  0.38  0.31  0.31  1  −2  1 × 10−22  3  2  2400  130000  1  1  0.6  0  4  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.2  0.105  30000  −0.2758  −0.2758  0.038  0.31  0.31  1  −0.015  1 × 10−22  3  2  2400  130000  1  1  0.6  5  5  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.06  30500  −0.2758  −0.2758  0.02  0.33  0.33  1  −0.1  1 × 10−22  3  2  2400  130000  1  1  0.6  20  6  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.12  6000  −0.2758  −0.2758  0.27  0.31  0.31  1  −0.13  1 × 10−22  3  2  2400  130000  1  1  0.6  0  7  Shale  2.71  1.02  0.6  −1  56  1.3  47  0.012  0.077  0.29  0.03  3800  −0.2758  −0.2758  0.059  0.33  0.33  1  −0.21  1 × 10−22  3  2  2400  36000  1  1  0.6  0  8  Sandstone  2.65  1.02  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  395  −0.2758  −0.2758  0.15  0.13  0.13  1  −1.2  6.61x10−11  5  2  2100  36000  1  1  0.6  set to 110  8  Sandstone (Pay)  2.65  0.81  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  225  −0.2758  −0.2758  0.23  0.13  0.13  1  −1.28  6.61x10−11  5  2  2400  36000  1  0.638  0.6  set to 110  9  Dolostone  2.811  1.02  2  −100  67  1.76  55  0.008  0.06  0.005  0.01  59500  −0.2758  −0.2758  0  0.25  0.25  1  −0.03  4.20x10−19  5  2  2400  200000  1  1  0.6  0  Layer  Lithology  Grain &  SR3 model  Yield Surface Evolution  Elastic Properties  Permeability Properties    No.    fluid densities  (Equations B1-B3)  (Equations B4-B7)  (Equations B8 & B9)  (Equation B10)  Consolidation Properties      ρg  ρf  pc  pt  β  n  ψ  κ  λ  ϕref  ϕinit  Eref  A  B  n  υmin  υmax  m  c  K0  x  y  Kf  Ks  α  η  Kxy  Pp      (g cc−1)  (g cc−1)  (Mpa)  (Mpa)                (MPa)                (m2)      (MPa)  (MPa)    (MPa s−1)    (MPa)  1  Shale  2.69  1.02  0.17  −1  60  1.3  51  0.012  0.106  0.55  0.38  910  −0.2758  −0.2758  0.085  0.31  0.31  1  −0.9  4.20 × 10−19  3  2  2400  36000  1  1  0.6  0  2  Sandstone  2.65  1.02  0.28  −1  66  1.7  56  0.012  0.174  0.55  0.33  1850  −0.2758  −0.2758  0.15  0.31  0.31  1  −0.5  6.61 × 10−11  5  2  2400  26000  1  1  0.6  0  3  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.35  0.2  410  −0.2758  −0.2758  0.38  0.31  0.31  1  −2  1 × 10−22  3  2  2400  130000  1  1  0.6  0  4  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.2  0.105  30000  −0.2758  −0.2758  0.038  0.31  0.31  1  −0.015  1 × 10−22  3  2  2400  130000  1  1  0.6  5  5  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.06  30500  −0.2758  −0.2758  0.02  0.33  0.33  1  −0.1  1 × 10−22  3  2  2400  130000  1  1  0.6  20  6  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.12  6000  −0.2758  −0.2758  0.27  0.31  0.31  1  −0.13  1 × 10−22  3  2  2400  130000  1  1  0.6  0  7  Shale  2.71  1.02  0.6  −1  56  1.3  47  0.012  0.077  0.29  0.03  3800  −0.2758  −0.2758  0.059  0.33  0.33  1  −0.21  1 × 10−22  3  2  2400  36000  1  1  0.6  0  8  Sandstone  2.65  1.02  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  395  −0.2758  −0.2758  0.15  0.13  0.13  1  −1.2  6.61x10−11  5  2  2100  36000  1  1  0.6  set to 110  8  Sandstone (Pay)  2.65  0.81  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  225  −0.2758  −0.2758  0.23  0.13  0.13  1  −1.28  6.61x10−11  5  2  2400  36000  1  0.638  0.6  set to 110  9  Dolostone  2.811  1.02  2  −100  67  1.76  55  0.008  0.06  0.005  0.01  59500  −0.2758  −0.2758  0  0.25  0.25  1  −0.03  4.20x10−19  5  2  2400  200000  1  1  0.6  0  View Large The time-lapse model results for six different production years is shown Fig. 5. Shown are plots of the change in vertical displacement Δz and vertical effective stress $$\Delta \sigma ^{\prime }_{v}$$ between the initial model state (at time 0) and the corresponding production year (i.e. monitor–baseline). Our geomechanical model estimates total reservoir compaction of roughly 0.8 m, surface subsidence of up to 0.2 m and a reduction in overburden vertical effective stress within the range 0–7 MPa. These values are typical of HPHT scenarios and thus render our synthetic model a good representation of a North Sea production scenario. Figure 5. View largeDownload slide Predicted change in vertical displacement Δz (left) in metres and vertical effective stress $$\Delta \sigma ^{\prime }_{v}$$ (right) in MPa from the initial pre-production state (i.e. 0 yr) for six separate production years: 1, 2, 3, 5, 10 and 20. For interpretation of colour in this figure, the reader is referred to the web version of this paper. Figure 5. View largeDownload slide Predicted change in vertical displacement Δz (left) in metres and vertical effective stress $$\Delta \sigma ^{\prime }_{v}$$ (right) in MPa from the initial pre-production state (i.e. 0 yr) for six separate production years: 1, 2, 3, 5, 10 and 20. For interpretation of colour in this figure, the reader is referred to the web version of this paper. We estimate overburden preproduction seismic traveltimes using the in situ P-wave velocity (Fig. 4). Seismic traveltimes are affected by changes in stress and path length. Thus, we calculate production induced vertical time-shifts using the modelled Δz and $$\Delta \sigma ^{\prime }_{v}$$ (Fig. 5). To model the influence of production related stress change on seismic velocity requires a rock physics model (e.g. Angus et al. 2015). These are often calibrated using stress–velocity relationships derived from core data and differ depending on lithology. This is an important step in coupling observed time-lapse anomalies to mechanical subsurface changes and has proven to be a challenging task (e.g. Price et al. 2016). However, for the purposes of this study we simplify this step by assuming each lithological layer in our model has a similar stress–velocity relationship which is known. We set dv/dσ΄ = 0.004 km s−1 MPa−1 which is a fair representation of the (effective) stress–velocity relationship considering typical core-measurements (e.g. Price et al. 2016). Specifically, we calculate changes to individual overburden layer vertical traveltimes Δtv as a result of production induced changes in layer thickness H and seismic (i.e. P-wave) velocity V. The results can be seen in Fig. 6. Also shown in Fig. 6 is each layer's R-factor which describes the proportion of fractional changes in layer thickness H to fractional changes in velocity V (Hatchel & Bourne 2005). Over the entire 20 yrs of production we estimate a total vertical time-shift to top reservoir of approximately 10 ms. Note that we focus on vertical displacements and vertical stress and thus near-offset (i.e. vertical) vertical, seismic traveltimes. This is because accurately modelling offset-dependant time-lapse time-shifts is difficult and currently in debate (e.g. Kudarova & Macbeth 2016). Figure 6. View largeDownload slide Fractional change to overburden layers thickness H and P-wave velocity V over the whole 20 yrs production period (i.e. Monitor-Base). Also shown is each layer's R-factor (i.e. R = [ΔV/V]/[ΔH/H]) and resultant change in each layer's traveltime Δtv. Note that these values are taken at location 2 in Fig. 1, directly above the reservoir region. The reservoir region is shaded. Figure 6. View largeDownload slide Fractional change to overburden layers thickness H and P-wave velocity V over the whole 20 yrs production period (i.e. Monitor-Base). Also shown is each layer's R-factor (i.e. R = [ΔV/V]/[ΔH/H]) and resultant change in each layer's traveltime Δtv. Note that these values are taken at location 2 in Fig. 1, directly above the reservoir region. The reservoir region is shaded. 3 GLOBAL SENSITIVITY ANALYSIS METHODS In this study we use a total of four different GSA methods. Each approach defines and measures sensitivity differently, capturing different aspects of the models response. This results in different, yet complimentary, sensitivity measures for the same input factor. Using a number of different GSA methods allows for a robust analysis of the models response and avoids the potential bias conclusions drawn from using just one specific method. In this section we describe the four GSA methods applied to our model: the Elementary Effects Test (EET), Regional Sensitivity Analysis (RSA), Variance-Based Sensitivity Analysis (VBSA) and a density based approach called PAWN. We also discuss how we modify the VBSA and PAWN methods such that they can be applied to the same data set as that of the RSA. This avoids additional model runs through tailored sampling strategies. Each of these four GSA methods, along with other GSA techniques, is described in greater detail in Petropoulos & Srivastava (2016). 3.1 Elementary Effects Test The Elementary Effects Test (EET; Morris 1991), calculates an effect per input from a one-factor-at-a-time (OAT) sample matrix, x:   \begin{eqnarray} x_{j,k} = {\left(\begin{array} {cccc}x_{1,1} &\quad x_{1,2} &\quad \cdots &\quad x_{1,k} \\ x_{2,1} &\quad x_{2,2} &\quad \cdots &\quad x_{2,k} \\ \vdots &\quad \vdots &\quad \ddots &\quad \vdots \\ x_{j,1} &\quad x_{j,2} &\quad \cdots &\quad x_{j,k} \end{array}\right)} , \end{eqnarray} (1)where, k is equal to the total number of parameters and j = k + 1 representing an independent sample or model run. The sample matrix x is ordered such that its first row (i.e. j = 1) is a randomly sampled set of model parameters whilst its jth row differs in only the (j − 1)th element. An Elementary Effect (EE) is calculated for each parameter k via   \begin{eqnarray} {\rm EE}_{i} = \frac{|Y_{i+1} - Y_{1}|}{(|x_{i+1,i} - x_{1,i}|)} (a_{i}\,-\,b_{i}) \quad {\rm {for}} \quad i=1,....,k , \end{eqnarray} (2)where Y is a 1× j matrix of each independent model result and a and b, 1×k matrices that define the maximum and minimum sample ranges for each parameter k respectively. Repeating this procedure, generating an ensemble n of x matrices, builds a finite distribution of n EE's per parameter k that is, EEn, i. To build a distribution of n elementary effects per input k would require n different x matrices and hence n(k + 1) model runs. A large (absolute) measure of central tendency in these EE distributions indicates an input with an important ‘overall’ influence on the output whilst, a large spread indicates an input with an important nonlinear effect (i.e. it is heavily affected by the values of other inputs and their interactions). This is shown schematically in Fig. 7. Note that the EET requires a tailored sampling strategy for the generation of the sample matrix x. Figure 7. View largeDownload slide Elementary Effect (EE) distributions of three different parameters x1, x2 and x3. A large (absolute) measure of the central tendency (i.e. mean value μ) indicates an input with an important direct influence on the model output, while a large spread (i.e. standard deviation SD) indicates an input with a strong nonlinear effect. Therefore, parameters that fall within the top right-hand section of an EE μ-SD plot are most influential to the model output. Figure 7. View largeDownload slide Elementary Effect (EE) distributions of three different parameters x1, x2 and x3. A large (absolute) measure of the central tendency (i.e. mean value μ) indicates an input with an important direct influence on the model output, while a large spread (i.e. standard deviation SD) indicates an input with a strong nonlinear effect. Therefore, parameters that fall within the top right-hand section of an EE μ-SD plot are most influential to the model output. 3.2 Regional Sensitivity Analysis RSA (Spear & Hornberger 1980) requires the separation of the input space into ‘behavioural’ and ‘non-behavioural’ regions. Formally, the set Xb of behavioural inputs is defined as   $$X_{b} = \lbrace x|y_{i} = f_{i}(x)\le \bar{y_{i}} {\rm \, for \, all}\,\, i \rbrace ,$$ (3)where x = [x1, ⋅⋅⋅, xk] is the vector of all k input parameters, yi either model output or an objective function (i.e. models fit to observed data) and $$\bar{y_{i}}$$ a predefined threshold value. Note that this particular criterion lends itself to assigning behavioural inputs as those that minimize a pre-defined objective function (i.e. distance between measured and observed data). However, different, less harsh criterions can be defined. For example, behavioural inputs can be defined such that they meet only one of many pre-defined threshold values. For this study we define behavioural samples to be those which show absolute differences from the data of less than the average absolute difference seen across the whole ensemble. Once the input sample is decomposed, sensitivity is measured by comparing the marginal Cumulative Density Functions (CDFs) of the two groups. Specifically, sensitivity is defined by means of the Kolmogorov–Smirnov (KS) statistic. The sensitivity index for the kth input factor xk is expressed as:   $$S_k = \max \limits _{x_k} |F_{k}^{B}(x_k)- F_{k}^{\bar{B}}(x_k)|,$$ (4)where $$F_{k}^{B}(x_k)$$ and $$F_{k}^{\bar{B}}(x_k)$$ are the behavioural and non-behavioural CDFs, respectively. A schematic demonstrating the KS statistic of two different CDFs is shown in Fig. 8. The larger the distance between the two CDFs (i.e. the larger the KS statistic) the greater the sensitivity. Note that, unlike the EET, the RSA does not require a tailored sampling strategy but only a generic input–output data set. Figure 8. View largeDownload slide Cumulative Density Functions (CDF's) of behavioural and non-behavioural samples. Different criterion can be used to define behavioural regions of the parameter space. Typically behavioural samples are those which minimize a pre-defined objective function such as the difference between measured and observed data. The Kolmogorov–Smirnov (KS) statistic describes the difference between the two CDF's, which in this study we take to be the maximum difference. The larger the KS statistic the larger the Regional Sensitivity Analysis (RSA) index. Figure 8. View largeDownload slide Cumulative Density Functions (CDF's) of behavioural and non-behavioural samples. Different criterion can be used to define behavioural regions of the parameter space. Typically behavioural samples are those which minimize a pre-defined objective function such as the difference between measured and observed data. The Kolmogorov–Smirnov (KS) statistic describes the difference between the two CDF's, which in this study we take to be the maximum difference. The larger the KS statistic the larger the Regional Sensitivity Analysis (RSA) index. 3.3 Variance–Based Sensitivity Analysis VBSA assigns a sensitivity index to each input parameter based upon its contribution to the variance of the model output (Sobol 1990). The direct contribution of the kth input factor to the variance of the output is defined as   $$S_{k} = \frac{V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]}{V(y)},$$ (5)where E is the expected value, V the variance and x∼ k a vector of all input factors but the kth. Sk can be described as the reduction of the total model output variance V(y) that would be observed on average when the uncertainty about xk would be removed by setting xk to a fixed value (Tarantola 2002). Since an analytic solution to eq. (5) is typically impossible, numerical approximations are often used (e.g. Saltelli et al. 2010) which require tailored sampling techniques. However, Petropoulos & Srivastava (2016) approximate eq. (5) such that it can be used on a generic input–output data set. They approximate $$E_{x_{\sim k}}(y|x_{k})$$ as a linear combination of Radial Basis Functions (RBFs),   $$\hat{E_{k}} = \sum _{j=1}^{n}[a_{j} {\rm {exp}}(-(x_{k}-w_{j})^{2})],$$ (6)where aj and wj are parameters that define the shape of the RBF. The variance $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ in eq. (5) can now be approximated by the sample variance of $$\hat{E_{k}}$$ while V(y) approximated from the variance of the sample output y. To obtain $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ operationally for each input factor the steps are (1) calibrate the regression function of eq. (6) by calculating the best fit parameters aj and wj (in our case we use a linear combination of five RBFs, thus j = 1, …, 5), (2) evaluate the optimized regression function for all values of xk and finally, (3) calculate the sample variance of $$\hat{E_{k}}$$. A schematic example of this methodology is shown in Fig. 9. Figure 9. View largeDownload slide A linear combination of Radial Basis Functions (RBFs) $$\hat{E_{k}}$$ is used as a regression function for the input-output (i.e. xk − Y) data set (left). The optimized regression function is then evaluated at all values of xk and the variance of $$\hat{E_{k}}$$ (right) used to approximate the term $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ in eq. (6). Figure 9. View largeDownload slide A linear combination of Radial Basis Functions (RBFs) $$\hat{E_{k}}$$ is used as a regression function for the input-output (i.e. xk − Y) data set (left). The optimized regression function is then evaluated at all values of xk and the variance of $$\hat{E_{k}}$$ (right) used to approximate the term $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ in eq. (6). 3.4 PAWN sensitivity analysis PAWN (Pianosi & Wagener 2015b) is a density-based method where sensitivity is measured by estimating the variation to the output y distribution when removing the uncertainty in one or more parameters xk. This variation is calculated from the measure of distance between the unconditional (when all inputs vary simultaneously) and conditional (when all inputs vary but xk, which is set to a nominal value) CDFs. The PAWN sensitivity index for the kth input is defined as:   $$S_{k} = \max _{x_{k}}\max _{y}|F_{y}(y)-F_{y|x_k}(y|x_{k})|,$$ (7)where Fy(y) and $$F_{y|x_{k}}(y|x_{k})$$ are the unconditional and conditional CDFs of the output respectively. The inner maximum of eq. (7) defines the maximum absolute difference between the two CDFs approximated via the KS statistic using empirical distribution functions. As the KS statistic will depend on the nominal (i.e. fixed) value of xk, the outer maximum of eq. (7) extracts the maximum KS statistic over all values of xk. If the data set does not contain multiple samples with the same value of xk, that is, a generic input-output data set, conditional distributions can be conditioned on ‘similar’ values of xk. Therefore eq. (7) can be approximated as   $$S_{k} = \max _{j=1,....,n}\max _{y}|F_{y}(y)-F_{y|x_{k}}(y|x_{k} \in \alpha _{j})|,$$ (8)where αj are n (e.g. 10) equally spaced intervals over the range of variation of xk. A schematic example of this method is shown in Fig. 10. Figure 10. View largeDownload slide Red line (left image) indicates the unconditional (when all inputs vary simultaneously) CDF while shaded lines the conditional CDFs (all inputs vary but xk) when xk is fixed at incremental nominal values. The KS statistic (see caption for Fig. 8) is computed for each unconditional–condition CDF pair and the PAWN sensitivity index taken as the maximum KS value for the input xk (right). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 10. View largeDownload slide Red line (left image) indicates the unconditional (when all inputs vary simultaneously) CDF while shaded lines the conditional CDFs (all inputs vary but xk) when xk is fixed at incremental nominal values. The KS statistic (see caption for Fig. 8) is computed for each unconditional–condition CDF pair and the PAWN sensitivity index taken as the maximum KS value for the input xk (right). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. 4 EXPERIMENTAL SETUP: DEFINING MODEL INPUT FACTORS AND OUTPUT We consider a single unknown overburden chalk layer (layer 5 in Table 1) whose material properties are largely uncertain. We assume little a priori knowledge of its properties and thus large, uniform, independent uncertainty distributions for all its physical parameters. The production profile and mechanical properties of the reservoir and underburden are assumed known. In other words, known to greater degree of accuracy (i.e. much smaller uncertainty) that of the overburden chalk. Thus, their properties are kept constant and not considered in the GSA. In total, 21 different input factors are subjected to the GSA, each a material property of the overburden chalk. A summary of these parameters is presented in Table 2 along with their range of uncertainties. These uncertainty ranges chosen are based on typical properties presented in the generic material database of Rockfield Software Limited and are made as wide as possible. Most of these parameters are described in greater detail in the definition of an elastoplastic material found in Appendix B. However, to ease GSA parameter space sampling we consider the poroelastic parameters A and B (eq. B8) to be equal (for simplicity) and a minimum Poisson's ratio (υmin in eq. B9) as a ratio υratio of υmax. Also, we assume that overburden rocks behave elastically during production as production related stress changes in the overburden are generally small compared to the yield strength of the rock. Therefore, we assume those parameters that define the shape of state boundary surface to have no effect on the seismic traveltimes. However, we do vary the yield surface evolution parameters as they will affect the stress dependant porosity parameter in eq. (B8). It should also be noted that for the purpose of this study we assume all unphysical, or algorithm specific (e.g. coupling rate, mesh type and size, etc.) parameters are optimized and are not considered in the GSA. Table 2. Chalk layer physical properties and their parameter sensitivity ranges. No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Kxy  Vertical-Horizontal stress coef    0.6  0.4  1.1  2  Pp (Mpa)  Over/underpressure    20  0  +40  3  ρf (g cc−1)  Fluid density    1.02  1  1.2  4  ρg (g cc−1)  Grain density    2.71  2.6  2.8  5  λ  Cam-clay constant  (B4) and (B5)  0.06  0.02  0.1  6  κ  Cam-clay constant  (B4) and (B5)  0.008  0.002  0.012  7  Eref (Mpa)  Reference Young's Modulus  (B8)  30500  1 × 104  3.3 × 104  8  n  Elastic constant  (B8)  0.02  0.001  0.1  9  A/B  Elastic constant  (B8)  −0.2758  0  −0.5  10  c  Elastic constant  (B8)  −0.1  −0.5  −0.001  11  υmax  Max Poisson's ratio  (B9)  0.33  0.2  0.4  12  υratio  Min Poisson's ratio  (B9)  1  1  1.5  13  m  Elastic constant  (B9)  1  0.01  1  14  α  Biot constant    1  0.5  1  15  K0 (m2)  Permeability constant  (B10)  1 × 10−22  1 × 10−23  1 × 10−18  16  x  Permeability constant  (B10)  3  1  4  17  y  Permeability constant  (B10)  2  1  6  18  Ks (Mpa)  Grain stiffness    2400  2400  5000  19  Kf (Mpa)  Fluid stiffness    13 × 104  8 × 104  15 × 104  20  Φinit  Initial porosity  (B6)  0.06  0.01  0.23  21  Φref  Reference porosity  (B6)  0.3  0.3  0.5  No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Kxy  Vertical-Horizontal stress coef    0.6  0.4  1.1  2  Pp (Mpa)  Over/underpressure    20  0  +40  3  ρf (g cc−1)  Fluid density    1.02  1  1.2  4  ρg (g cc−1)  Grain density    2.71  2.6  2.8  5  λ  Cam-clay constant  (B4) and (B5)  0.06  0.02  0.1  6  κ  Cam-clay constant  (B4) and (B5)  0.008  0.002  0.012  7  Eref (Mpa)  Reference Young's Modulus  (B8)  30500  1 × 104  3.3 × 104  8  n  Elastic constant  (B8)  0.02  0.001  0.1  9  A/B  Elastic constant  (B8)  −0.2758  0  −0.5  10  c  Elastic constant  (B8)  −0.1  −0.5  −0.001  11  υmax  Max Poisson's ratio  (B9)  0.33  0.2  0.4  12  υratio  Min Poisson's ratio  (B9)  1  1  1.5  13  m  Elastic constant  (B9)  1  0.01  1  14  α  Biot constant    1  0.5  1  15  K0 (m2)  Permeability constant  (B10)  1 × 10−22  1 × 10−23  1 × 10−18  16  x  Permeability constant  (B10)  3  1  4  17  y  Permeability constant  (B10)  2  1  6  18  Ks (Mpa)  Grain stiffness    2400  2400  5000  19  Kf (Mpa)  Fluid stiffness    13 × 104  8 × 104  15 × 104  20  Φinit  Initial porosity  (B6)  0.06  0.01  0.23  21  Φref  Reference porosity  (B6)  0.3  0.3  0.5  View Large The outputs analysed by the GSA are the individual overburden layer vertical time-shifts Δtv; specifically, Δtv over each production year (i.e. 1 through to 20 yrs) at three separate locations (highlighted in Fig. 1). Sensitivity indices are generated for each parameter by analysing each Δtv output. In this study, we present global sensitivity indices by averaging combined sets of individual results. For example, the sensitivity of chalk layer 5 to a certain parameter will be the average of the three individual indices calculated at the three locations shown in Fig. 1. Therefore, the sensitivity indices take into account variations in Δtv across the model domain. 5 RESULTS Utilizing over 4000 model perturbations we conduct an in depth multimethod GSA using the SAFE Toolbox of Pianosi et al. (2015a). In our GSA, we initially use the EET on an ensemble of 1540 (i.e. n = 70 EE's per input) model runs to screen those parameters (see Table 2) that have little effect over overburden Δtv. We then create an additional random ensemble of 3000 model runs with a reduced number of variable parameters. The RSA, VBSA and PAWN techniques are then applied to this data set simultaneously to rank the most influential parameters in order of importance. We then extract the most influential parameters (hereby referred to as ‘active’ parameters) and assess their potential seismic calibration by model performance (i.e. mapping their location within the model space). 5.1 Screening model parameters Fig. 11 shows the results of the EET analysis. Taking the three vertical locations shown in Fig. 1 we calculate average EE measures for both the uncertain chalk layer and remaining overburden layers at 1 yr intervals over the total 20 yrs of production. As discussed in Section 3.1, a large (absolute) measure of central tendency (i.e. mean) indicates an input with an important direct influence on the model output. A large spread (i.e. standard deviation) indicates an input with a strong nonlinear effect on the model output. Thus, parameters which show a significant shading have a greater influence over the modelled Δtv. Figure 11. View largeDownload slide Time varying Elementary Effects (EEs) considering the resultant change in layer traveltime Δtv at yearly intervals over the total 20 yrs of production. Results are computed considering the Δtv results of the uncertain chalk layer only and the Δtv result of the other (unchanged) overburden layers at the locations specified in Fig. 1. Figure 11. View largeDownload slide Time varying Elementary Effects (EEs) considering the resultant change in layer traveltime Δtv at yearly intervals over the total 20 yrs of production. Results are computed considering the Δtv results of the uncertain chalk layer only and the Δtv result of the other (unchanged) overburden layers at the locations specified in Fig. 1. It is apparent from Fig. 11 that altering the material properties of a single layer affects Δtv across the whole overburden. This demonstrates a complex nonlinear model behaviour. However, similar sensitivity patterns emerging across all overburden layers suggests the total modelled overburden Δtv could potentially be explained by a handful of model variables. Although Fig. 11 suggest certain model parameters to be more influential than others, it is difficult to conclusively screen a number of parameters as being non-influential. However, it is apparent from these results that the Cam-Clay parameters λ and κ do not affect the modelled Δtv. Both parameters measure zero EE mean and standard deviation. This is not unexpected as they primarily affect the yield surface evolution (eqs B4 and B5) and, as their influence is negligible, it confirms our assumption that the overburden in our model remains elastic during production. However, their zero measure also suggests negligible influence over the stress-dependant porosity parameter in eq. (B8). As a result we screen out the Cam-Clay parameters λ and κ as non-influential but consider all remaining 19 parameters as potentially influential to Δtv. 5.2 Ranking model parameters To further investigate the influence of the remaining uncertain parameters on the modelled Δtv, we run an additional 3000 different model input combinations. We use the same uncertainty ranges as expressed in Table 2, but fix the Cam-Clay constants to their true value (Table 1). We create 3000 ensemble by randomly sampling the parameter space using Latin-hypercube sampling (Forrester et al. 2008). We apply the RSA, VBSA and PAWN sensitivity techniques simultaneously to this input–output data set to give complimentary parameter sensitivity indices. Fig. 12 shows the results of all three GSA methods focused on the final modelled Δtv (i.e. after 20 yrs of production). Note that here we use a similar procedure as that of the EET by summarizing sensitivity indices as average values for the uncertain chalk layer and remaining overburden layers over the locations shown in Fig. 1. Figure 12. View largeDownload slide GSA sensitivity indices of the reduced set of model parameters. Blue circles represent the RSA results, black squares PAWN and grey hollow boxes the VBSA results. Sensitivity indices are computed considering the Δtv results of the uncertain chalk layer only and the results of the remaining (unchanged) overburden layers. These results focus on the final model Δtv that is, after 20 yrs of production, at the locations highlighted in Fig. 1. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 12. View largeDownload slide GSA sensitivity indices of the reduced set of model parameters. Blue circles represent the RSA results, black squares PAWN and grey hollow boxes the VBSA results. Sensitivity indices are computed considering the Δtv results of the uncertain chalk layer only and the results of the remaining (unchanged) overburden layers. These results focus on the final model Δtv that is, after 20 yrs of production, at the locations highlighted in Fig. 1. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. A parameter with a greater sensitivity index indicates one which has a greater direct influence over Δtv. It is apparent from Fig. 12 that, although not giving exactly the same absolute measures of sensitivity, all three GSA techniques provide suitably similar global trends. We observe similar results as that of the EET analysis (Fig. 11) where the parameters that control the nonlinear elastic response (Eref, c, ϕinit or eq. B8) and the Biot coefficient α are notably the most influential across all overburden layers. It is therefore fair to assume these to be the active parameters of our model. If we plot the time-varying sensitivity indices of these parameters when considering Δtv of the chalk, we can compare the influence of the elastic properties to α over production. The results are presented in Fig. 13. Here, we see the elastic parameters to be most influential during earlier production times but they appear to be outweighed by α later in production. Figure 13. View largeDownload slide The time-varying GSA sensitivity indices of the four most influential parameters within the uncertain chalk layer. The elastic parameters are shown in blue whilst the Biot coefficient α in red. For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this paper. Figure 13. View largeDownload slide The time-varying GSA sensitivity indices of the four most influential parameters within the uncertain chalk layer. The elastic parameters are shown in blue whilst the Biot coefficient α in red. For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this paper. 5.3 Mapping data to model space We take the 3000 model ensemble and compare these simulations to the result of our original model run. Assuming the overburden Δtv of the original model (Fig. 6) to be observed data, we compare model residuals in the form of a Δtv Root Mean Square Error (RMSE). As time-lapse seismic data is time consuming and costly to acquire, shooting data at yearly intervals (or less) is generally unfeasible. Therefore, we compute the RMSE considering the Δtv results at just three production time steps. The results of the GSA (Figs 11 and 13) show that the active variables do not start to significantly affect overburden Δtv till 3 yrs of production. They also show that the Biot coefficient α takes approximately 10 yrs to become as influential as the elastic coefficients (e.g. Fig. 13). Taking this into consideration we choose to compute the RMSE using the Δtv result at 3, 5 and 10 yrs. Note that the RMSE is calculated assuming the results from all three vertical locations shown in Fig. 1. Fig. 14 show parallel coordinate plots of the best 5 per cent (i.e. 15) of models whose Δtv measurements most closely resemble that of the data (i.e. our original, truth model). These possess the lowest RMSE and are referred to as behavioural models. Also shown in Fig. 14 are their corresponding Δtv logs (plotted just at location 2 of Fig. 1) after 10 and 20 yrs of production. The results show that each behavioural model Δtv result closely resembles that of the data. Even their forward predicted Δtv values (i.e. after 20 yrs of production) are similar to those observed in the data. However, these models appear randomly scattered along the uncertain parameter range. Thus, models with significantly different active parameters each produce similar, possibly acceptable, solutions. This suggests a complex, possibly ill-posed, model space. If we simplify the objective function to the RMSE of just the uncertain chalk layers Δtv (i.e. ignoring the Δtv of other overburden layers) we see a slightly different result (Fig. 14). This optimization produces a different set of behavioural models which, as expected, do a better job of fitting the data of the uncertain chalk layer. These behavioural models also appear less scattered throughout the model space. Almost all of their active parameters more closely resemble those of the original (i.e. true) model. However, these models contain a significantly lower value of the elastic coefficient Eref (and as a result a lower pre-production Young's Modulus). Thus, significantly different model parameters still produce similar, possibly acceptable, solutions even with a simplified, condensed data space. Figure 14. View largeDownload slide Parallel coordinate plots showing the active parameters of the best 5 per cent (i.e. 15) of models whose Δtv results most closely resemble that of original, that is truth, model (Fig. 6). Also shown are the corresponding models Δtv overburden logs after both 10 and 20 yrs of production. The original model results are highlighted in red, while the closest models in black. The overburden Δtv logs of the whole model ensemble shown in grey. The model residuals were computed by taking the Root Mean Square Error (RMSE) of the whole overburden Δtv results after 3, 5 and 10 yrs of production. Also shown are the results when only the uncertain chalk layers Δtv results are considered in the residual calculation. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 14. View largeDownload slide Parallel coordinate plots showing the active parameters of the best 5 per cent (i.e. 15) of models whose Δtv results most closely resemble that of original, that is truth, model (Fig. 6). Also shown are the corresponding models Δtv overburden logs after both 10 and 20 yrs of production. The original model results are highlighted in red, while the closest models in black. The overburden Δtv logs of the whole model ensemble shown in grey. The model residuals were computed by taking the Root Mean Square Error (RMSE) of the whole overburden Δtv results after 3, 5 and 10 yrs of production. Also shown are the results when only the uncertain chalk layers Δtv results are considered in the residual calculation. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Fig. 14 along with the sensitivity results of the GSA emphasizes the complexity and nonlinearity of the model behaviour. The uncertainty in a single overburden layer's material properties appears to affect the predicted Δtv across the whole overburden domain. Due to significant pore pressure changes, the reservoir undergoes far more extreme mechanical changes than the overburden during production. Thus, it is reasonable to presume that the uncertainty in the mechanical properties of the reservoir may, more heavily, influence overburden Δtv than similar uncertainty in overburden properties. To test this hypothesis, we assume the mechanical properties of our reservoir to be uncertain, whilst holding the properties of the overburden constant (Table 1). We take the four active parameters highlighted by the GSA and assume uniform, independent uncertainty distributions for their values in the reservoir. A summary of their uncertainty ranges is outlined in Table 3. We assume large, uniform, independant uncertainty distributions for each parameter, similar to that used for the overburden chalk in the GSA (Table 2). Using the same Latin Hypercube sampling as that used in the GSA (Forrester et al. 2008), we run 200 different model input combinations and compute their overburden Δtv results. The Δtv values are again computed at the same locations as that used in the GSA (Fig. 1). Table 3. Active parameters of the reservoir layer and their uncertain sensitivity ranges. No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Eref (Mpa)  Reference Young's Modulus  (B8)  225  100  1000  2  c  Elastic constant  (B8)  −1.28  −1.5  −0.001  3  α  Biot constant    1  0.5  1  4  Φinit  Initial porosity  (B6)  0.12  0.01  0.23  No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Eref (Mpa)  Reference Young's Modulus  (B8)  225  100  1000  2  c  Elastic constant  (B8)  −1.28  −1.5  −0.001  3  α  Biot constant    1  0.5  1  4  Φinit  Initial porosity  (B6)  0.12  0.01  0.23  View Large Shown in Fig. 15 are the corresponding overburden logs for each model run after 10 yrs of production. The results confirm the significant affect reservoir uncertainty has on overburden Δtv. Its uncertainty appears to influence overburden Δtv to a much greater extent than suitable similar uncertainty in overburden layers (Fig. 14). It is important to note that the extreme mechanical changes seen in the reservoir could result in plastic deformation. Thus, parameters which govern the rocks state boundary surface could potentially have significant influence. These parameters should be included if an in depth sensitivity study is to be undertaken for the reservoir region. Figure 15. View largeDownload slide Overburden Δtv logs of the uncertain reservoir ensemble (Table 3) after 10 yrs of production shown via grey lines. The, truth model results are shown in red (Fig. 6) while the dotted lines represent the extreme values seen within the original GSA results (Fig. 14). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 15. View largeDownload slide Overburden Δtv logs of the uncertain reservoir ensemble (Table 3) after 10 yrs of production shown via grey lines. The, truth model results are shown in red (Fig. 6) while the dotted lines represent the extreme values seen within the original GSA results (Fig. 14). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. 6 DISCUSSION The results of the GSA highlighted that out of an initial 21 model parameters, the modelled overburden Δtv is most heavily influenced by just 4 (e.g. Figs 11 and 12). The Biot coefficient (α) and the parameters that govern a materials elastic behaviour (Eref, c and ϕinit). The GSA also shows that these variables take 3 yrs (from the start of production) to become significantly influential and that α takes approximately 10 yrs to become as influential as the elastic parameters (e.g. Fig. 13). This time varying sensitivity demonstrates how acquiring multiple vintages of time-lapse seismic data could be advantageous. However, these results also suggest that significant overburden near-offset seismic time-shifts may take time to manifest. Thus, time-lapse seismic data taken shortly after the start of production may not be conclusive enough to aid the advanced numerical calibration of geomechanical models. However, acquiring early seismic data can be beneficial as an early warning system. It can highlight any large discrepancies between model and reality, and suggest the potential need for additional data (e.g. well log or core data). The elastic parameters being influential is not totally unexpected, as seismic traveltimes are affected by changes in stress and path length and thus governed by rock stiffness (i.e. Young's Modulus, eq. B8). However, slightly unexpected is the significant influence of the Biot coefficient α. Typically, the overburden is modelled as an undrained scenario (i.e. no fluid flow) in which you assume that there is no production related pore pressure change. Therefore, you would expect α (i.e. σ΄ = σ − αPp) to have little influence over changes in effective stress and hence Δtv. However, ELFEN implicitly evaluates the pore pressure Pp of the whole model domain (i.e. whole domain coupling) as a function of α and the volumetric strain εv. Therefore, although there is no (or little) fluid flow outside of the reservoir, the pore pressure is affected by mechanical changes in volumetric strain. The GSA demonstrates this and emphasizes the importance of modelling fluid flow and pore pressure outside of the reservoir. Their slight instabilities can cause non-negligible effects to the model output. The consequence of the material properties of a single layer affecting Δtv across the whole overburden demonstrates a complex nonlinear model behaviour. Thus, analysing model activity globally (i.e. across whole modelled domain), as opposed to locally, could be crucial for potential calibration. Although we highlight just four active parameters for the Δtv results, it is also influenced, albeit to a lesser extent, by the remaining uncertainty in other parameters (e.g. Fig. 12). Therefore it must be stressed that changes in overburden Δtv are not entirely determined by changes in these four active variables. Thus, calibration procedures focussing on a condensed model space should also account for variations to model output caused by changes in less sensitive variables. Models with significantly different input parameters produce similar Δtv results (e.g. Fig. 14). This is the case when considering both global Δtv results (i.e. across the whole overburden) and when focused to local model output (i.e. results of uncertain chalk layer only). This highlights the complexity of the model space where a single global solution will most likely not exist. Instead, numerous models, of different input combinations will produce equally acceptable solutions (e.g. Fig. 14). In this study, we assume only the uncertainty of a single overburden layer. The complexity of the model space will undoubtedly increase when we consider the uncertainty in the mechanical properties of other layers. The GSA results do suggest that time-lapse seismic data could potentially be able to distinguish between certain models within our ensemble. The difference in Δtv between certain models being on the order of 2–3 ms (e.g. Fig. 14). However, the size and complexity of the model space makes calibration via Monte Carlo based approaches unfeasible. For example, it is clear from the results of the GSA that 3000 randomly sampled model perturbations are insufficient to determine definitive solutions (e.g. Fig. 14). However, with a more sophisticated method of sampling of the model space, we may be able to determine the most plausible solutions from a more feasible number of model runs. Although still unlikely to determine a single unique solution, you may potentially be able to determine a set of potential scenarios within sensible time frames. This could provide vital additional information for use in conjunction with qualitative analysis. We find that the uncertainty in the mechanical properties of the reservoir heavily influence overburden Δtv (e.g. Fig. 15). Its uncertainty appearing to influence Δtv to a much greater extent than suitable similar uncertainty in the overburden (e.g. Fig. 14). This demonstrates the nonlinearity of the model behaviour and the importance of a suitably accurate reservoir model. Time-lapse seismic calibration of other properties of a geomechanical model will thus only be possible once the reservoir behaviour is known to a suitable degree of accuracy. It is important to note that the sensitivity measurements of the GSA are heavily affected by the uncertainty in model parameters. For example, considering a much smaller uncertainty range in the elastic coefficients would result in their sensitivity being significantly less than seen in this study. Therefore it is always important to cross analyse the results of the GSA with the uncertainty range used. It is also important to highlight that we concentrate our analysis on near-offset that is, vertical traveltimes. If we consider time-shift versus offset estimations we would expect an increased sensitivity in parameters such as the Poisson's ratio and the horizontal-to-vertical stress ratio Kxy. It is also important to mention that in this study we do not include the uncertainty in the stress-dependant rock physics model. We assume the in situ seismic velocity of each rock, and its stress-dependence, is known (i.e. no uncertainty). In reality, this process is highly uncertain, especially in overburden rocks where stress–velocity core data is often not available (e.g. Price et al. 2016). This complex modelling step is beyond the scope of this study. However, any future quantitative calibration of geomechanical models with near offset time-lapse seismic data must include this uncertainty in its calculations. Finally, in this study we have not accounted for random modelling errors (e.g. variations to the implicit and explicit solutions caused by the parametrization of their solvers), which can act as noise or bias to the resulting output distributions. However, since large distributions are seen in Δtv (e.g. Figs 14 and 15) it is safe to assume these random modelling errors to be insignificant. 7 CONCLUSION In this study, we have analysed the potential of near-offset time-lapse seismic data to aid the geomechanical calibration of the overburden. We build a geomechanical model of a typical North Sea HPHT reservoir and utilize over 4000 model perturbations to conduct an in-depth multimethod GSA. We outline that out an initially large set of uncertain material properties, the modelled overburden Δtv are mainly affected by just four ‘active’ parameters. These being the effective stress (i.e. Biot) coefficient (α) and the parameters that govern the material's elastic behaviour (i.e. stiffness). The influence of the Biot coefficient highlights the importance of modelling fluid flow and pore pressure outside of the reservoir. These results show that the model space can be condensed to these active variables for calibration. However, the variations caused by less sensitive variables are non-negligible and so should also be taken into account. We demonstrate how the FE model is complex and highly nonlinear. Altering the material properties of a single layer affects the Δtv results of the whole overburden domain. As a result, constraining the uncertainty in model parameters appears challenging. Multiple combinations of model parameters can yield equally possible model realizations and hence Δtv results. Consequently, calibration via a large number of random model perturbations is unfeasible. However, the significant difference in the Δtv estimates of certain models within our ensemble suggests more sophisticated calibration methods could potentially be feasible. Determining a set of potential solutions using a small number of model runs could be possible given intelligent sampling of the model space. The results of the time-varying GSA suggest that significant overburden time-shifts may take time to manifest. Thus, time-lapse seismic data taken shortly after the start of production may not be conclusive enough to aid the advanced numerical calibration of geomechanical models. We find that similar uncertainty in reservoir mechanical properties (as to overburden properties) appear to influence overburden Δtv to a much greater extent. This demonstrates the complex nonlinear model behaviour and stresses the importance of reservoir characterization. Thus, to calibrate parts of the geomechanical model other than the reservoir, such as the overburden, the reservoir behaviour must be known to a suitable degree of accuracy. This lends credit to small scale sensitivity studies of reservoir uncertainty before overburden calibration is considered. Although this particular study focuses on the modelling software ELFEN. The findings and conclusions can be related to all FE geomechanical modelling software. We stress that, although not straightforward, there is potential to quantitatively calibrate geomechanical models via time-lapse seismic data. Further studies like this could extend models to include anisotropic behaviour and fault properties. Acknowledgements DP is supported by a scholarship funded by Total GRC and the School of Earth and Environment, University of Leeds. The author would like to thank Francesca Pianosa of Department of Civil Engineering, University of Bristol for providing a copy of the SAFE toolbox and for her advice on the use of different sensitivity techniques. The author would also like to thank Rockfield Ltd for providing the University of Leeds with an academic license for their modelling software ELFEN. DA acknowledges the Research Council UK (EP/K035878/1; EP/K021869/1; NE/L000423/1) for financial support. REFERENCES Alerini M., Ayzenberg M., Ek T., Feng T., Hustoft L., Lie E., Liu S., Skjei N., Skjervheim J.A., 2014. Utilization of Time-lapse seismic for reservoir model conditioning, in76th EAGE Conference and Exhibition 2014  EAGE. Angus D.A. et al.  , 2015. Integrated hydro-mechanical and seismic modelling of the Valhall reservoir: a case study of predicting subsidence, AVOA and microseismicity, Geomech. Energy Environ.  2 32– 44. Google Scholar CrossRef Search ADS   Bloch S., Lander R.H., Bonnell L., 2002. Anomalously high porosity and permeability in deeply buried sandstone reservoirs: origin and predictability, AAPG Bull.  86( 2), 301– 328. Bruno M.S., 2002. Geomechanics and decision analysis for migration compaction related casing damage, SPE drilling and completion, SPE  79519 179– 188. Calvert R., 2005. Insights and Methods for 4D Reservoir Monitoring and Characterization  Distinguished Instructor Series No. 8, EAGE. Google Scholar CrossRef Search ADS   Carman P.C., 1937. Fluid flow through a granular bed, Trans. Inst. Chem. Eng.  15 150– 156. Crook A.J.L., Willson S.M., Yu J.G., Owen D.R.J., 2006. Predictive modelling of structure evolution in sandbox experiments, J. Struct. Geol.  28( 5), 729– 744. Google Scholar CrossRef Search ADS   De Gennaro S., Onaisi A., Grandi A., Ben-Brahim L., Neillo V., 2008. 4D reservoir geomechanics: a case study from the HP/HT reservoirs of the Elgin and Franklin fields, First Break  26 53– 59. Google Scholar CrossRef Search ADS   De Gennaro S., Schutjens P., Frumau M., Fuery M., Ita J., Fokker P., 2010. The role of geomechanics in the development of an HPHT field, in 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium  Salt Lake City, UT, American Rock Mechanics Association, pp. 10– 450. Emerick A.A., Reynolds A.C., 2012. History matching time-lapse seismic data using the ensemble Kalman filter with multiple data assimilations, Comput. Geosci.  16( 3), 639– 659. Google Scholar CrossRef Search ADS   Fjær E., Kristiansen T.G., 2009. An Integrated Geomechanics, Rock Physics and Seismic model, in 71st EAGE Conference & Exhibition  Amsterdam. Forrester A., Sobester A., Keane A., 2008. Engineering Design via Surrogate modelling: A Practical Guide  John Wiley & Sons. Garcia A.G., 2011. Dynamic reservoir characterization from overburden time-lapse strains, PhD thesis  Institute of Petroleum Engineering, Heriot-Watt University. Gonzalez-Carballo A., Guyonnet P.Y., Levallois B., Veillerette A., Deboiasne R., 2006. 4D monitoring in Angola and its impact on reservoir understanding and economics, Leading Edge  25 1150– 1159. Google Scholar CrossRef Search ADS   Gosselin O. et al.  , 2003. History matching using time-lapse seismic, in SPE Annual Technical Conference and Exhibition  Society of Petroleum Engineers. Guilbot J., Smith B., 2002. 4-D constrained depth conversion for reservoir compaction estimation: Application to Ekofisk Field, Leading Edge  21 302– 308. Google Scholar CrossRef Search ADS   Hatchell P.J., Bourne S.J., 2005. Rocks under strain: Strain-induced time-lapse time shifts are observed for depleting reservoirs, leading Edge  24( 12), 1222– 1225. Google Scholar CrossRef Search ADS   Hawkins K. et al.  , 2007. Production-induced stresses from time-lapse time shifts: A geomechanical case study from Franklin and Elgin fields, Leading Edge  26 655– 662. Google Scholar CrossRef Search ADS   Herwanger J., Horne S.A., 2009. Linking reservoir geomechanics and time-lapse seismics: predicting anisotropic velocity changes and seismic attributes, Geophysics  74 W13– W33. Google Scholar CrossRef Search ADS   Jackson D., Richardson M., 2007. High-frequency Seafloor Acoustics  Springer Science & Business Media. Google Scholar CrossRef Search ADS   Jaeger J.C., Cook N.G., Zimmerman R., 2009. Fundamentals of Rock Mechanics  John Wiley & Sons. Japsen P., 1998. Regional velocity-depth anomalie, North Sea chalk: a record of overporessure and neogene uplift and erosion, AAPG Bull  82( 11), 2031– 2074. Japsen P., 2000. Investigation of multi-phase erosion using reconstructed shale trends based on sonic data, Sole Pit axis, North Sea, Glob. Planet. Change  24( 3), 189– 210. Google Scholar CrossRef Search ADS   Jones J.C., 2010. Hydrocarbons-physical Properties and Their Relevance to Utilisation  BookBoon, pp. 1– 100. Kozeny J., 1927. Über kapillare leitung des wassers im boden, Sitzungsber. Akad. Wiss. Wien  136 271– 306. Kvendseth S.S., 1988. Giant Discovery, A history of Ekofisk through the First 20 Years  Philips Petroleum Company. Kudarova A.M., Macbeth C., 2016. Offset-dependence of production-related 4D time shifts: Real data examples and modelling, in 87th SEG Annual Meeting  Dallas, TX. Landrø M., 2001. Discrimination between pressure and fluid saturation changes from time-lapse seismic data, Geophysics  66 836– 844. Google Scholar CrossRef Search ADS   Mallon A.J., Swarbrick R.E., 2008. Diagenetic characteristics of low permeability, non-reservoir chalks from the Central North Sea, Mar. Pet. Geol.  25( 10), 1097– 1108. Google Scholar CrossRef Search ADS   Marcussen ø., Maast T.E., Mondol N.H., Jahren J., Bjørlykke K., 2010. Changes in physical properties of a reservoir sandstone as a function of burial depth—the Etive formation, northern North Sea, Mar. Pet. Geol.  27( 8), 1725– 1735. Google Scholar CrossRef Search ADS   Minkoff S.E., Stone C.M., Bryant S., Peszynska M., Wheeler M.F., 2003. Coupled fluid flow and geomechanical deformation modeling, J. Pet. Sci. Eng.  38 37– 56. Google Scholar CrossRef Search ADS   Morris M.D., 1991. Factorial sampling plans for preliminary computational experiments, Technometrics  33( 2), 161– 174. Google Scholar CrossRef Search ADS   Okiongbo K.S., 2011. Petrophysical properties of North Sea shales, Res. J. Appl. Sci. Eng. Technol.  3( 1), 46– 52. Petropoulos G., Srivastava P.K., 2016. Sensitivity Analysis in Earth Observation Modelling  Elsevier Science & Technology Books. Pianosi F., Wagener T., 2015. A simple and efficient method for global sensitivity analysis based on cumulative distribution functions, Environ. Modelling Softw.  67 1– 11. Google Scholar CrossRef Search ADS   Pianosi F., Sarrazin F., Wagener T., 2015. A Matlab toolbox for Global Sensitivity Analysis, Environ. Modelling Softw.  70 80– 85. Google Scholar CrossRef Search ADS   Price D.C., Angus D.A., Garcia A., Fisher Q.J., 2016. Probabilistic analysis and comparison of stress-dependent rock physics models, Geophys. J. Int.  210( 1), 196– 209. Google Scholar CrossRef Search ADS   Rao S.S., 2010. The Finite Element Method in Engineering  Elsevier. Rockfield Software Limited, 2012. ELFEN GeoDB Generic Materials  Swansea. Saltelli A., Annoni P., Azzini I., Campolongo F., Ratto M., 2010. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, Comput. Phys. Commun.  181( 2), 259– 270. Google Scholar CrossRef Search ADS   Samier P., Onaisi A., Fontaine G., 2003. Coupled analysis of geomechanics and fluid flow in reservoir simulation, SPE Reservoir simulation symposium  Houston, USA, 3–5 Febuary, SPE 79698. Schneider F., Potdevin J.L., Wolf S., Faille I., 1996. Mechanical and chemical compaction model for sedimentary basin simulators, Tectonophysics  263( 1), 307– 317. Google Scholar CrossRef Search ADS   Slagstad T., Barróre C., Davidsen B., Ramstad R.K., 2008. Petrophysical and thermal properties of pre-Devonian basement rocks on the Norwegian continental margin, Geol. Surv. Norway Bull.  448 1– 6. Sobol I., 1990. Sensitivity analysis for non-linear mathematical models. Math. Modelling Comput. Exp.  1 407– 424. Spear R.C., Hornberger G.M., 1980. Eutrophication in peel inlet II. Identification of critical uncertainties via generalized sensitivity analysis, Water Res.  14( 1), 43– 49. Google Scholar CrossRef Search ADS   Staples R., Stevens T., Leoutre E., Jolley S., Marshall J., 2005. 4D seismic history matching-the reality, in 67th EAGE Conference and Exhibition  Madrid. Staples R., Ita J., Nash R., Hague P., Burrell R., 2007. Using 4D seismic data and geomechanical modelling to understand pressure depletion in HPHT fields of the Central North Sea, in 69th EAGE Conference and Exhibition  EAGE. Tarantola S., Giglioli N., Jesinghaus J., Saltelli A., 2002. Can global sensitivity analysis steer the implementation of models for environmental assessments and decision-making?, Stoch. Environ. Res. Risk Assess.  16( 1), 63– 76. Google Scholar CrossRef Search ADS   Tarantola A., 2005. Inverse Problem Theorey and Methods for Model Parameter Estimation  SIAM. Google Scholar CrossRef Search ADS   Vudovich A., Chin L.V., Morgan D.R., 1989. Casing deformation in Ekofisk, J. Pet. Technol.  41 729– 734. Google Scholar CrossRef Search ADS   Wright W.A., 1967. Prediction of bulk moduli and pressure-volume-temperature data for petroleum oils, ASLE Trans.  10( 4), 349– 356. Google Scholar CrossRef Search ADS   Zang J., Reederm R.J., 1999. Comparative compressibilities of calcite-structure carbonates: deviations from empirical relations, Am. Mineral.  84( 5–6), 861– 870. Google Scholar CrossRef Search ADS   APPENDIX A: COUPLING Hydro-mechanical modelling in ELFEN involves iterative coupling. Solutions for the mechanical domain are calculated explicitly whereas the fluid domain is solved implicitly via a nonlinear Newton–Raphson approach (e.g. Tarantola 2005). Solutions for both fields are solved simultaneously via transferring the data between the two domains at specified time intervals. The coupling rate for the hydro-mechanical process is similar to that of the implicit time step. This is typically 400–500 times larger than the explicit step. A robust explicit time step Δt can be estimated via:   $$\Delta t = f_{{\rm crit}} \times {\rm min} \Bigg | l^{e} \sqrt{\frac{\rho ^{e}}{E^{e}}} \Bigg |.$$ A1$$l^{e} \sqrt{\rho ^{e}/E^{e}}$$ is defined as the critical time step where Ee, ρe and le are the Young's Modulus, density and characteristic length of the minimum element, respectively. The critical time step is a stability condition which prevents the magnification of round-off errors caused by the explicit scheme (e.g. Rao 2010). Its calculation is approximate and therefore we introduce the factor of critical step fcrit. Rockfield suggests that for 2-D hydro-mechanical models fcrit = 0.9. Using the approximated mechanical time step, we choose a coupling rate (i.e. the rate at which the explicit mechanical and implicit fluid field are coupled) of 0.01 time steps. This is well within the stability threshold for our model while not being too small such that it compromises computational runtime or numerical stability. APPENDIX B: ELASTOPLASTIC MATERIAL PROPERTIES B1 Basic properties We assume all non-reservoir rocks to have a fluid density ρf of 1.02 g cc−1 (Japsen 1998), while the reservoir rock has a fluid density of 0.81 g cc−1 to represent the density of light crude oil in the North Sea (Jones 2010). The grain densities ρg of each shale is assumed to be 2.69 g cc−1 (Okiongbo 2011), chalk assumed to be purely calcite with a grain density of 2.71 g cc−1 (Japsen 1998), each sandstone to be composed mainly of quartz, 2.65 g cc−1, and the stiff underburden a grain density of 2.81 g cc−1 typical of dolostone. B2 State boundary surface The state boundary surface used for each material is the Soft Rock (SR3) model of Crook et al. (2006). The SR3 yield function is a smooth, three-invariant surface that intersects the hydrostatic axis in both tension pt0 and compression pc0. It is defined as   $$\Phi (p)=g(\theta , p)q+(p-p_{t0})\tan \beta \bigg (\frac{p-p_{c0}}{p_{t}-p_{c0}}\bigg )^{1/n},$$ B1where p and q are the effective mean stress and deviatoric stress respectively, β and n material constants and θ the Lode angle. Finally, g(θ, p) is the deviatoric plane correction term that controls the shape of the yield surface in the deviatoric plane. The evolution of the plastic flow is defined by a non-associated flow rule:   $$\dot{\varepsilon ^{p}} = \dot{\lambda } \frac{\partial \Psi }{\partial \sigma },$$ B2where $$\dot{\lambda }$$ is the plastic multiplier and Ψ is the plastic potential defined as:   $$\Psi (p)=g(\theta , p)q+(p-p_{t0}){\rm tan}\psi \bigg (\frac{p-p_{c0}}{p_{t0}-p_{c0}}\bigg )^{1/n}.$$ B3Note that eq. (B3) is of identical form to that of the state boundary surface defined in eq. (B1). However, the plastic potential is defined in terms of the angle tanψ, where ψ is the dilation parameter controlling the shape of the plastic potential surface. The deviatoric plane correction term g is scaled to be 1 such that the strength in triaxial compression directly corresponds to the strength calibrated using compressive triaxial (CTC) tests (Crook et al. 2006). The initial state boundary surface is defined at a reference porosity ϕref (i.e. surface conditions). To define the state boundary surface for each material we use a compilation of test data presented in the generic material database of Rockfield Software Limited. A summary of the final chosen parameter values presented in Table 1. B3 Yield surface evolution The evolution (or hardening) of the primary yield surface is determined through relationships that define pc and pt as a function of volumetric plastic strain $$\varepsilon _{v}^{p}$$:   $$p_{c} = p_{c0} + (p_{c0}-p_{c({\rm resid})})\bigg [{\rm exp}\bigg (-\frac{v\Delta \varepsilon _{v}^{p}}{(\lambda -\kappa )}\bigg )-1\bigg ],$$ B4  \begin{eqnarray} p_{t}^* &=& p_{t0} + (p_{t0}-p_{t({\rm resid})})\bigg [{\rm exp}\bigg (-\frac{v\Delta \varepsilon _{v}^{p}}{(\lambda -\kappa )}\bigg )-1\bigg ],\nonumber \\ p_{t} &=& {\rm max}\left[p_{t0}, p_{t}^* \right]. \end{eqnarray} B5Here, κ and λ are Cam-Clay hardening constants and v is the specific volume. v can be related to porosity via 1/(1 − ϕ). Note that pc(resid) = pc0/100 and pt(resid) = pc0/100 to ensure the yield surface is always of finite size. The volumetric plastic strain $$\varepsilon _{v}^{p}$$ can also be defined in terms of porosity via:   $$\varepsilon _{v}^{p} = {\rm {log}}\bigg [\frac{1-\phi _{{\rm ref}}}{(1-\phi _{{\rm init}})}\bigg ].$$ B6These hardening relationships allow a material characterization defined at surface conditions, with a specific reference porosity ϕref, to be used to generate data suitable for a similar material at greater depth (subjected to compaction) with a different, initialization porosity ϕinit. The shape of the state boundary surface remains unchanged but its size is governed by the scaling of pc0 to pc(init) and pt0 to pt(init). We use the compilation of test data presented in the generic material database of Rockfield Software Limited to define κ and λ. The initialization porosity ϕinit corresponds to the porosity of the material at the start of the simulation. Its value for each material is determined through either porosity-depth relationships for North Sea rocks (Bloch et al. 2002; Mallon & Swarbrick 2008) or estimated using typical North Sea bulk density ρ values (Japsen 1998, 2000; Slagstad et al. 2008; Marcussen et al. 2010; Okiongbo 2011) via:   $$\rho =\rho _g(1-\phi _{{\rm init}})+\rho _f\phi _{{\rm init}}.$$ B7A complete summary of yield surface evolution parameter values are presented in Table 1. B4 Elastic properties Young's modulus E and Poisson's ratio υ are defined in ELFEN as empirical functions of effective mean stress $$\hat{\sigma}$$:   $$E = E_{{\rm ref}}\bigg [\frac{p+A}{B}\bigg ]^n \phi (\hat{\sigma})^{c},$$ B8  $$\upsilon = \upsilon _{{\rm min}} + (\upsilon _{{\rm max}} + \upsilon _{{\rm min}})(1-e^{m\hat{\sigma}}).$$ B9Eref is the effective Young's Modulus while υmax and υmin are the Poisson's ratio values at low and high effective mean stress $$\hat{\sigma}$$. n, c and m are material constants, while A and B are also material constants used to prevent problems near zero values of $$\hat{\sigma}$$. Finally ϕ($$\hat{\sigma}$$) is the porosity, which itself is a function of effective mean stress $$\hat{\sigma}$$. These elastic relationships are calibrated such that the Young's Modulus and Poisson's ratio of each material, prior to simulation, resemble those stated in Garcia 2011 and their nonlinear behaviour calibrated using the generic material database of Rockfield Software Limited, 2012. Again, all elastic parameters can be found in Table 1. Note that we have initially specified a constant Poisson's ratio for this study as we assume negligible changes to Poisson's ratio with effective mean stress. B5 Porous flow properties The porosity dependent intrinsic permeability Kin is based upon the Kozeny–Carman (Kozeny 1927; Carman 1937) relationship   $$K_{{\rm in}}(\phi ) = K_{0} \frac{\phi ^x}{(1-\phi )^y},$$ B10where K0, x and y are material constants. The chosen parameters are provided in Table 1 and based on the values given by Schneider et al. (1996). B6 Consolidation properties Additional material consolidation properties are also required by ELFEN. An effective stress coefficient (i.e. Biot) α of 1 is chosen for all rock types (a parameter usually set to 1 in most modelling scenarios). A fluid viscosity η is also required and set to 1×10−9 MPa   s−1 (typical of pure water) in all non-reservoir layers. Within the reservoir we use η = 0.638 × 10−9 MPa s−1, similar to the viscosity of hydrocarbons from the North Sea Forties field (Jones 2010). Fluid Kf and grain Ks stiffnesses are also required. For these we use typical values found in North Sea literature (Wright 1967; Zang & Reeder 1999; Jackson & Richardsons 2007; Jaeger et al. 2009). Finally a horizontal to vertical stress ratio Kxy is prescribed at 0.6 for all materials. In addition to the required consolidation properties, we include an over/under-pressure parameter Pp. This alters the settled in situ (post-geostatic initialization) pore pressure state of each material. In addition to the overpressured reservoir region (e.g. Fig. 2), we overpressure two chalk layers in the overburden, one by +20 MPa and another by +5 MPa. A summary of all the consolidation properties can be found in Table 1. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# A multimethod Global Sensitivity Analysis to aid the calibration of geomechanical models via time-lapse seismic data

, Volume 212 (3) – Mar 1, 2018
16 pages

/lp/ou_press/a-multimethod-global-sensitivity-analysis-to-aid-the-calibration-of-IUaBexz8Tx
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx516
Publisher site
See Article on Publisher Site

### Abstract

Summary Time-lapse seismic attributes are used extensively in the history matching of production simulator models. However, although proven to contain information regarding production induced stress change, it is typically only loosely (i.e. qualitatively) used to calibrate geomechanical models. In this study we conduct a multimethod Global Sensitivity Analysis (GSA) to assess the feasibility and aid the quantitative calibration of geomechanical models via near-offset time-lapse seismic data. Specifically, the calibration of mechanical properties of the overburden. Via the GSA, we analyse the near-offset overburden seismic traveltimes from over 4000 perturbations of a Finite Element (FE) geomechanical model of a typical High Pressure High Temperature (HPHT) reservoir in the North Sea. We find that, out of an initially large set of material properties, the near-offset overburden traveltimes are primarily affected by Young's modulus and the effective stress (i.e. Biot) coefficient. The unexpected significance of the Biot coefficient highlights the importance of modelling fluid flow and pore pressure outside of the reservoir. The FE model is complex and highly nonlinear. Multiple combinations of model parameters can yield equally possible model realizations. Consequently, numerical calibration via a large number of random model perturbations is unfeasible. However, the significant differences in traveltime results suggest that more sophisticated calibration methods could potentially be feasible for finding numerous suitable solutions. The results of the time-varying GSA demonstrate how acquiring multiple vintages of time-lapse seismic data can be advantageous. However, they also suggest that significant overburden near-offset seismic time-shifts, useful for model calibration, may take up to 3 yrs after the start of production to manifest. Due to the nonlinearity of the model behaviour, similar uncertainty in the reservoir mechanical properties appears to influence overburden traveltime to a much greater extent. Therefore, reservoir properties must be known to a suitable degree of accuracy before the calibration of the overburden can be considered. Geomechanics, Numerical modeling, Statistical Methods 1 INTRODUCTION A reduction in reservoir pore pressure, as a result of production, can lead to significant compaction. A decrease in reservoir thickness results in the straining of the surrounding rock mass and thus a change to the in situ stress state (e.g. Herwanger & Horne 2009; Fjær & Kristiansen 2009). Failure to anticipate and manage geomechanical issues associated with production can have a detrimental effect on the economic performance of a field. For example, wellbore failure due to fault reactivation in the overburden (e.g. Vudovich et al. 1989; Bruno 2002) and platform instability due to seabed subsidence (e.g. Kvendseth 1988) are all but some of the consequences reported in literature over the past few decades. Therefore, geomechanical modelling can be an essential reservoir monitoring tool. Understanding and anticipating the spatial changes to the in situ stress field can aid production and help maintain well and platform integrity. Estimating production induced stress change typically requires the coupling of a geomechanical model to a reservoir simulator (e.g. Minkoff et al. 2003; Samier et al. 2003; De Gennaro et al. 2010). The coupled model allows the modelling of pore pressure and saturation change within the reservoir (in the fluid domain) while additionally simulating changes in stress and strain (in the mechanical domain) to the reservoir and surrounding rock mass. Often the coupled model is constrained with the use of production data. However, time-lapse seismic data have been found to contain information about reservoir saturation (e.g. Landrø 2001; Calvert 2005; Gonzalez-Carballo et al. 2006; Alerini et al. 2014) and production induced stress change (e.g. Guilbot & Smith 2002). As a result, time-lapse seismic data have been used extensively to aid the calibration of both the reservoir (e.g. Staples et al. 2005) and geomechanical model (e.g. Hawkins et al. 2007; Staples et al. 2007; De Gennaro et al. 2008; Angus et al. 2015). However, unlike advanced calibration procedures applied to production simulation models, such as, advanced history matching (e.g. Gosselin et al. 2003; Emerick & Reynolds 2012), the geomechanical model is typically loosely calibrated, that is, both model and data show similar patterns of anomalies (e.g. Hatchel & Bourne 2005). By not utilizing numerical matching techniques, are we making the most out of the geomechanical information stored in time-lapse seismic data? To investigate, we conduct a multimethod Global Sensitivity Analysis (GSA) on over 4000 model perturbations of a Finite Element (FE) geomechanical model of a typical High Pressure High Temperature (HPHT) reservoir in the North Sea. The GSA results are used to investigate the behaviour of geomechanical models (i.e. evaluate the potential geomechanical information content of time-lapse seismic data), assess the feasibility of seismic calibration and aid future numerical matching studies. We focus our attention primarily to the overburden where the mechanical properties are often as unknown as those of the reservoir rock and considered to be equally as important (e.g. Vudovich et al. 1989; Bruno 2002). A suitable example being the uncertain mechanical properties of overburden chalks seen in North Sea reservoirs. However, overburden time-lapse seismic anomalies are not complicated by fluid effects (i.e. saturation change as seen within the reservoir). Therefore, they are assumed to be a purely mechanical consequence and thus a suitable focus for this study. Specifically, via the GSA, we examine the sensitivity of overburden seismic time-shifts to the various properties of a single unknown overburden layer in a geomechanical model. We attempt to screen model parameters with negligible influence (i.e. those with little influence on time-lapse seismic data), rank those that are most influential (i.e. those with the most influence over time-lapse seismic data and hence the main focus of calibration) and map the input parameter space (i.e. can the uncertainty in influential parameters be reduced by seismic data). 2 THE GEOMECHANICAL MODEL We use the FE modelling software ELFEN (Rockfield Software Ltd) to construct a geomechanical model typical of an HPHT reservoir in the North Sea. We also utilize ELFEN's single-phase production simulator to model hydromechanical coupling (described in greater detail in Appendix A). We simplify our model to 2-D and isotropic due to the large number of model perturbations that need to be computed for the GSA. Although ELFEN has the potential to model anisotropy, we assume an isotropic scenario to reduce the number of model parameters for this initial study. The model domain consists of a 20 × 9 km2 subsurface region, with a reservoir interval located at a depth of approximately 5 km. This depth is similar to that of HPHT fields in the North Sea. Multiple overburden layers and faults are included to increase model complexity. However, for this particular study the faults are not initialized and are left as simple lithological discontinuities. This helps reduce the number of model parameters analysed in the GSA. Fig. 1 shows the complete model geometry. The reservoir interval is overpressured to 110 MPa prior to production and depleted by roughly 50 per cent over a period of 20 yrs. The rate of depletion is shown in Fig. 2 and is based on the production profile of the HPHT reservoir given by Hawkins et al. (2007). Figure 1. View largeDownload slide Model geometry, with reservoir region highlighted in blue and faults in red. The chalk layer whose physical properties are deemed uncertain shaded in grey. Also shown are the three locations in which we estimate vertical traveltimes used in the Global Sensitivity Analysis (GSA). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 1. View largeDownload slide Model geometry, with reservoir region highlighted in blue and faults in red. The chalk layer whose physical properties are deemed uncertain shaded in grey. Also shown are the three locations in which we estimate vertical traveltimes used in the Global Sensitivity Analysis (GSA). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 2. View largeDownload slide Reservoir pore pressure reduction given in normalized units, 1 = 110 MPa. Figure 2. View largeDownload slide Reservoir pore pressure reduction given in normalized units, 1 = 110 MPa. The in situ stress and pore pressure state, prior to production, is determined through a period of settling time steps in which the model is loaded under gravity. Any disequilibrium caused by the loading on the geometry is given time to relax, avoiding numerical oscillations (i.e. unwanted dynamic effects) in the modelling results. The phreatic surface (pore pressure origin) is set 50 m above the surface of the model to mimic a shallow North Sea environment. The model mesh is unstructured with triangular elements approximately 70 m in size. Element size was determined such that the model produced stable results (free from unwanted dynamic effects), had acceptable computation time, and yielded suitable resolution (considering typical vertical and lateral resolution of seismic data between 1 and 5 km depth). Fig. 3 shows the CPU time for models with varying sized meshes. Note the exponential increase in CPU time per relatively small increments in mesh size. Figure 3. View largeDownload slide Model runtime versus mesh element size. Figure 3. View largeDownload slide Model runtime versus mesh element size. We design nine different homogeneous isotropic elastoplastic materials, one for each of the nine lithological layers of our model (Fig. 1). In keeping with North Sea geology, we assume a sandstone reservoir with an impermeable shale caprock and a thick layer of stiff chalk in the overburden. To simplify model behaviour, we substitute a complex salt underburden (typically seen in the North Sea) with a simple mechanically strong rock, typical of unfractured limestone/dolostone. To define an elastoplastic material in ELFEN requires the parametrization of a state boundary surface, nonlinear elastic and permeability relationships and a suite of consolidation properties. A detailed description of all these material property relationships and associated parameters can be found in Appendix B, whilst a definitive summary in Table 1. Material properties are determined using both the extensive literature available on North Sea geology and generic materials designed by Rockfield in their generic material database (e.g. Rockfield Software Limited 2012). Fig. 4 shows logs of the final geomechanical model properties, post-geostatic initialization and prior reservoir production. Using Young's Modulus, Poisson's ratio and bulk density, we estimate seismic P-wave velocities. Note that no static to dynamic conversion was used when calculating the P-wave velocity as the static mechanical properties give rise to credible dynamic velocities. Figure 4. View largeDownload slide Logs of Young's modulus, Poisson's ratio, bulk density and pore pressure through the final geomechanical model, post-geostatic initialization and prior to reservoir production (i.e. 0 yr). These logs are then used to calculate the P-wave Vp velocity profile. Layer boundaries are marked via dotted horizontal lines. The lithostatic and hydrostatic gradients are plotted in red on the pore pressure log. Reservoir layer is shaded. Note that no static to dynamic conversion was used when calculating Vp as the static mechanical properties give rise to credible velocities. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 4. View largeDownload slide Logs of Young's modulus, Poisson's ratio, bulk density and pore pressure through the final geomechanical model, post-geostatic initialization and prior to reservoir production (i.e. 0 yr). These logs are then used to calculate the P-wave Vp velocity profile. Layer boundaries are marked via dotted horizontal lines. The lithostatic and hydrostatic gradients are plotted in red on the pore pressure log. Reservoir layer is shaded. Note that no static to dynamic conversion was used when calculating Vp as the static mechanical properties give rise to credible velocities. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Table 1. Material properties for each layer in Fig. 1. The layer number represents increasing depth from the surface. Layer  Lithology  Grain &  SR3 model  Yield Surface Evolution  Elastic Properties  Permeability Properties    No.    fluid densities  (Equations B1-B3)  (Equations B4-B7)  (Equations B8 & B9)  (Equation B10)  Consolidation Properties      ρg  ρf  pc  pt  β  n  ψ  κ  λ  ϕref  ϕinit  Eref  A  B  n  υmin  υmax  m  c  K0  x  y  Kf  Ks  α  η  Kxy  Pp      (g cc−1)  (g cc−1)  (Mpa)  (Mpa)                (MPa)                (m2)      (MPa)  (MPa)    (MPa s−1)    (MPa)  1  Shale  2.69  1.02  0.17  −1  60  1.3  51  0.012  0.106  0.55  0.38  910  −0.2758  −0.2758  0.085  0.31  0.31  1  −0.9  4.20 × 10−19  3  2  2400  36000  1  1  0.6  0  2  Sandstone  2.65  1.02  0.28  −1  66  1.7  56  0.012  0.174  0.55  0.33  1850  −0.2758  −0.2758  0.15  0.31  0.31  1  −0.5  6.61 × 10−11  5  2  2400  26000  1  1  0.6  0  3  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.35  0.2  410  −0.2758  −0.2758  0.38  0.31  0.31  1  −2  1 × 10−22  3  2  2400  130000  1  1  0.6  0  4  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.2  0.105  30000  −0.2758  −0.2758  0.038  0.31  0.31  1  −0.015  1 × 10−22  3  2  2400  130000  1  1  0.6  5  5  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.06  30500  −0.2758  −0.2758  0.02  0.33  0.33  1  −0.1  1 × 10−22  3  2  2400  130000  1  1  0.6  20  6  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.12  6000  −0.2758  −0.2758  0.27  0.31  0.31  1  −0.13  1 × 10−22  3  2  2400  130000  1  1  0.6  0  7  Shale  2.71  1.02  0.6  −1  56  1.3  47  0.012  0.077  0.29  0.03  3800  −0.2758  −0.2758  0.059  0.33  0.33  1  −0.21  1 × 10−22  3  2  2400  36000  1  1  0.6  0  8  Sandstone  2.65  1.02  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  395  −0.2758  −0.2758  0.15  0.13  0.13  1  −1.2  6.61x10−11  5  2  2100  36000  1  1  0.6  set to 110  8  Sandstone (Pay)  2.65  0.81  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  225  −0.2758  −0.2758  0.23  0.13  0.13  1  −1.28  6.61x10−11  5  2  2400  36000  1  0.638  0.6  set to 110  9  Dolostone  2.811  1.02  2  −100  67  1.76  55  0.008  0.06  0.005  0.01  59500  −0.2758  −0.2758  0  0.25  0.25  1  −0.03  4.20x10−19  5  2  2400  200000  1  1  0.6  0  Layer  Lithology  Grain &  SR3 model  Yield Surface Evolution  Elastic Properties  Permeability Properties    No.    fluid densities  (Equations B1-B3)  (Equations B4-B7)  (Equations B8 & B9)  (Equation B10)  Consolidation Properties      ρg  ρf  pc  pt  β  n  ψ  κ  λ  ϕref  ϕinit  Eref  A  B  n  υmin  υmax  m  c  K0  x  y  Kf  Ks  α  η  Kxy  Pp      (g cc−1)  (g cc−1)  (Mpa)  (Mpa)                (MPa)                (m2)      (MPa)  (MPa)    (MPa s−1)    (MPa)  1  Shale  2.69  1.02  0.17  −1  60  1.3  51  0.012  0.106  0.55  0.38  910  −0.2758  −0.2758  0.085  0.31  0.31  1  −0.9  4.20 × 10−19  3  2  2400  36000  1  1  0.6  0  2  Sandstone  2.65  1.02  0.28  −1  66  1.7  56  0.012  0.174  0.55  0.33  1850  −0.2758  −0.2758  0.15  0.31  0.31  1  −0.5  6.61 × 10−11  5  2  2400  26000  1  1  0.6  0  3  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.35  0.2  410  −0.2758  −0.2758  0.38  0.31  0.31  1  −2  1 × 10−22  3  2  2400  130000  1  1  0.6  0  4  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.2  0.105  30000  −0.2758  −0.2758  0.038  0.31  0.31  1  −0.015  1 × 10−22  3  2  2400  130000  1  1  0.6  5  5  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.06  30500  −0.2758  −0.2758  0.02  0.33  0.33  1  −0.1  1 × 10−22  3  2  2400  130000  1  1  0.6  20  6  Chalk  2.71  1.02  2  −100  67  1.76  55  0.008  0.06  0.3  0.12  6000  −0.2758  −0.2758  0.27  0.31  0.31  1  −0.13  1 × 10−22  3  2  2400  130000  1  1  0.6  0  7  Shale  2.71  1.02  0.6  −1  56  1.3  47  0.012  0.077  0.29  0.03  3800  −0.2758  −0.2758  0.059  0.33  0.33  1  −0.21  1 × 10−22  3  2  2400  36000  1  1  0.6  0  8  Sandstone  2.65  1.02  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  395  −0.2758  −0.2758  0.15  0.13  0.13  1  −1.2  6.61x10−11  5  2  2100  36000  1  1  0.6  set to 110  8  Sandstone (Pay)  2.65  0.81  0.65  −1  66  1.7  56  0.012  0.129  0.5  0.12  225  −0.2758  −0.2758  0.23  0.13  0.13  1  −1.28  6.61x10−11  5  2  2400  36000  1  0.638  0.6  set to 110  9  Dolostone  2.811  1.02  2  −100  67  1.76  55  0.008  0.06  0.005  0.01  59500  −0.2758  −0.2758  0  0.25  0.25  1  −0.03  4.20x10−19  5  2  2400  200000  1  1  0.6  0  View Large The time-lapse model results for six different production years is shown Fig. 5. Shown are plots of the change in vertical displacement Δz and vertical effective stress $$\Delta \sigma ^{\prime }_{v}$$ between the initial model state (at time 0) and the corresponding production year (i.e. monitor–baseline). Our geomechanical model estimates total reservoir compaction of roughly 0.8 m, surface subsidence of up to 0.2 m and a reduction in overburden vertical effective stress within the range 0–7 MPa. These values are typical of HPHT scenarios and thus render our synthetic model a good representation of a North Sea production scenario. Figure 5. View largeDownload slide Predicted change in vertical displacement Δz (left) in metres and vertical effective stress $$\Delta \sigma ^{\prime }_{v}$$ (right) in MPa from the initial pre-production state (i.e. 0 yr) for six separate production years: 1, 2, 3, 5, 10 and 20. For interpretation of colour in this figure, the reader is referred to the web version of this paper. Figure 5. View largeDownload slide Predicted change in vertical displacement Δz (left) in metres and vertical effective stress $$\Delta \sigma ^{\prime }_{v}$$ (right) in MPa from the initial pre-production state (i.e. 0 yr) for six separate production years: 1, 2, 3, 5, 10 and 20. For interpretation of colour in this figure, the reader is referred to the web version of this paper. We estimate overburden preproduction seismic traveltimes using the in situ P-wave velocity (Fig. 4). Seismic traveltimes are affected by changes in stress and path length. Thus, we calculate production induced vertical time-shifts using the modelled Δz and $$\Delta \sigma ^{\prime }_{v}$$ (Fig. 5). To model the influence of production related stress change on seismic velocity requires a rock physics model (e.g. Angus et al. 2015). These are often calibrated using stress–velocity relationships derived from core data and differ depending on lithology. This is an important step in coupling observed time-lapse anomalies to mechanical subsurface changes and has proven to be a challenging task (e.g. Price et al. 2016). However, for the purposes of this study we simplify this step by assuming each lithological layer in our model has a similar stress–velocity relationship which is known. We set dv/dσ΄ = 0.004 km s−1 MPa−1 which is a fair representation of the (effective) stress–velocity relationship considering typical core-measurements (e.g. Price et al. 2016). Specifically, we calculate changes to individual overburden layer vertical traveltimes Δtv as a result of production induced changes in layer thickness H and seismic (i.e. P-wave) velocity V. The results can be seen in Fig. 6. Also shown in Fig. 6 is each layer's R-factor which describes the proportion of fractional changes in layer thickness H to fractional changes in velocity V (Hatchel & Bourne 2005). Over the entire 20 yrs of production we estimate a total vertical time-shift to top reservoir of approximately 10 ms. Note that we focus on vertical displacements and vertical stress and thus near-offset (i.e. vertical) vertical, seismic traveltimes. This is because accurately modelling offset-dependant time-lapse time-shifts is difficult and currently in debate (e.g. Kudarova & Macbeth 2016). Figure 6. View largeDownload slide Fractional change to overburden layers thickness H and P-wave velocity V over the whole 20 yrs production period (i.e. Monitor-Base). Also shown is each layer's R-factor (i.e. R = [ΔV/V]/[ΔH/H]) and resultant change in each layer's traveltime Δtv. Note that these values are taken at location 2 in Fig. 1, directly above the reservoir region. The reservoir region is shaded. Figure 6. View largeDownload slide Fractional change to overburden layers thickness H and P-wave velocity V over the whole 20 yrs production period (i.e. Monitor-Base). Also shown is each layer's R-factor (i.e. R = [ΔV/V]/[ΔH/H]) and resultant change in each layer's traveltime Δtv. Note that these values are taken at location 2 in Fig. 1, directly above the reservoir region. The reservoir region is shaded. 3 GLOBAL SENSITIVITY ANALYSIS METHODS In this study we use a total of four different GSA methods. Each approach defines and measures sensitivity differently, capturing different aspects of the models response. This results in different, yet complimentary, sensitivity measures for the same input factor. Using a number of different GSA methods allows for a robust analysis of the models response and avoids the potential bias conclusions drawn from using just one specific method. In this section we describe the four GSA methods applied to our model: the Elementary Effects Test (EET), Regional Sensitivity Analysis (RSA), Variance-Based Sensitivity Analysis (VBSA) and a density based approach called PAWN. We also discuss how we modify the VBSA and PAWN methods such that they can be applied to the same data set as that of the RSA. This avoids additional model runs through tailored sampling strategies. Each of these four GSA methods, along with other GSA techniques, is described in greater detail in Petropoulos & Srivastava (2016). 3.1 Elementary Effects Test The Elementary Effects Test (EET; Morris 1991), calculates an effect per input from a one-factor-at-a-time (OAT) sample matrix, x:   \begin{eqnarray} x_{j,k} = {\left(\begin{array} {cccc}x_{1,1} &\quad x_{1,2} &\quad \cdots &\quad x_{1,k} \\ x_{2,1} &\quad x_{2,2} &\quad \cdots &\quad x_{2,k} \\ \vdots &\quad \vdots &\quad \ddots &\quad \vdots \\ x_{j,1} &\quad x_{j,2} &\quad \cdots &\quad x_{j,k} \end{array}\right)} , \end{eqnarray} (1)where, k is equal to the total number of parameters and j = k + 1 representing an independent sample or model run. The sample matrix x is ordered such that its first row (i.e. j = 1) is a randomly sampled set of model parameters whilst its jth row differs in only the (j − 1)th element. An Elementary Effect (EE) is calculated for each parameter k via   \begin{eqnarray} {\rm EE}_{i} = \frac{|Y_{i+1} - Y_{1}|}{(|x_{i+1,i} - x_{1,i}|)} (a_{i}\,-\,b_{i}) \quad {\rm {for}} \quad i=1,....,k , \end{eqnarray} (2)where Y is a 1× j matrix of each independent model result and a and b, 1×k matrices that define the maximum and minimum sample ranges for each parameter k respectively. Repeating this procedure, generating an ensemble n of x matrices, builds a finite distribution of n EE's per parameter k that is, EEn, i. To build a distribution of n elementary effects per input k would require n different x matrices and hence n(k + 1) model runs. A large (absolute) measure of central tendency in these EE distributions indicates an input with an important ‘overall’ influence on the output whilst, a large spread indicates an input with an important nonlinear effect (i.e. it is heavily affected by the values of other inputs and their interactions). This is shown schematically in Fig. 7. Note that the EET requires a tailored sampling strategy for the generation of the sample matrix x. Figure 7. View largeDownload slide Elementary Effect (EE) distributions of three different parameters x1, x2 and x3. A large (absolute) measure of the central tendency (i.e. mean value μ) indicates an input with an important direct influence on the model output, while a large spread (i.e. standard deviation SD) indicates an input with a strong nonlinear effect. Therefore, parameters that fall within the top right-hand section of an EE μ-SD plot are most influential to the model output. Figure 7. View largeDownload slide Elementary Effect (EE) distributions of three different parameters x1, x2 and x3. A large (absolute) measure of the central tendency (i.e. mean value μ) indicates an input with an important direct influence on the model output, while a large spread (i.e. standard deviation SD) indicates an input with a strong nonlinear effect. Therefore, parameters that fall within the top right-hand section of an EE μ-SD plot are most influential to the model output. 3.2 Regional Sensitivity Analysis RSA (Spear & Hornberger 1980) requires the separation of the input space into ‘behavioural’ and ‘non-behavioural’ regions. Formally, the set Xb of behavioural inputs is defined as   $$X_{b} = \lbrace x|y_{i} = f_{i}(x)\le \bar{y_{i}} {\rm \, for \, all}\,\, i \rbrace ,$$ (3)where x = [x1, ⋅⋅⋅, xk] is the vector of all k input parameters, yi either model output or an objective function (i.e. models fit to observed data) and $$\bar{y_{i}}$$ a predefined threshold value. Note that this particular criterion lends itself to assigning behavioural inputs as those that minimize a pre-defined objective function (i.e. distance between measured and observed data). However, different, less harsh criterions can be defined. For example, behavioural inputs can be defined such that they meet only one of many pre-defined threshold values. For this study we define behavioural samples to be those which show absolute differences from the data of less than the average absolute difference seen across the whole ensemble. Once the input sample is decomposed, sensitivity is measured by comparing the marginal Cumulative Density Functions (CDFs) of the two groups. Specifically, sensitivity is defined by means of the Kolmogorov–Smirnov (KS) statistic. The sensitivity index for the kth input factor xk is expressed as:   $$S_k = \max \limits _{x_k} |F_{k}^{B}(x_k)- F_{k}^{\bar{B}}(x_k)|,$$ (4)where $$F_{k}^{B}(x_k)$$ and $$F_{k}^{\bar{B}}(x_k)$$ are the behavioural and non-behavioural CDFs, respectively. A schematic demonstrating the KS statistic of two different CDFs is shown in Fig. 8. The larger the distance between the two CDFs (i.e. the larger the KS statistic) the greater the sensitivity. Note that, unlike the EET, the RSA does not require a tailored sampling strategy but only a generic input–output data set. Figure 8. View largeDownload slide Cumulative Density Functions (CDF's) of behavioural and non-behavioural samples. Different criterion can be used to define behavioural regions of the parameter space. Typically behavioural samples are those which minimize a pre-defined objective function such as the difference between measured and observed data. The Kolmogorov–Smirnov (KS) statistic describes the difference between the two CDF's, which in this study we take to be the maximum difference. The larger the KS statistic the larger the Regional Sensitivity Analysis (RSA) index. Figure 8. View largeDownload slide Cumulative Density Functions (CDF's) of behavioural and non-behavioural samples. Different criterion can be used to define behavioural regions of the parameter space. Typically behavioural samples are those which minimize a pre-defined objective function such as the difference between measured and observed data. The Kolmogorov–Smirnov (KS) statistic describes the difference between the two CDF's, which in this study we take to be the maximum difference. The larger the KS statistic the larger the Regional Sensitivity Analysis (RSA) index. 3.3 Variance–Based Sensitivity Analysis VBSA assigns a sensitivity index to each input parameter based upon its contribution to the variance of the model output (Sobol 1990). The direct contribution of the kth input factor to the variance of the output is defined as   $$S_{k} = \frac{V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]}{V(y)},$$ (5)where E is the expected value, V the variance and x∼ k a vector of all input factors but the kth. Sk can be described as the reduction of the total model output variance V(y) that would be observed on average when the uncertainty about xk would be removed by setting xk to a fixed value (Tarantola 2002). Since an analytic solution to eq. (5) is typically impossible, numerical approximations are often used (e.g. Saltelli et al. 2010) which require tailored sampling techniques. However, Petropoulos & Srivastava (2016) approximate eq. (5) such that it can be used on a generic input–output data set. They approximate $$E_{x_{\sim k}}(y|x_{k})$$ as a linear combination of Radial Basis Functions (RBFs),   $$\hat{E_{k}} = \sum _{j=1}^{n}[a_{j} {\rm {exp}}(-(x_{k}-w_{j})^{2})],$$ (6)where aj and wj are parameters that define the shape of the RBF. The variance $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ in eq. (5) can now be approximated by the sample variance of $$\hat{E_{k}}$$ while V(y) approximated from the variance of the sample output y. To obtain $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ operationally for each input factor the steps are (1) calibrate the regression function of eq. (6) by calculating the best fit parameters aj and wj (in our case we use a linear combination of five RBFs, thus j = 1, …, 5), (2) evaluate the optimized regression function for all values of xk and finally, (3) calculate the sample variance of $$\hat{E_{k}}$$. A schematic example of this methodology is shown in Fig. 9. Figure 9. View largeDownload slide A linear combination of Radial Basis Functions (RBFs) $$\hat{E_{k}}$$ is used as a regression function for the input-output (i.e. xk − Y) data set (left). The optimized regression function is then evaluated at all values of xk and the variance of $$\hat{E_{k}}$$ (right) used to approximate the term $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ in eq. (6). Figure 9. View largeDownload slide A linear combination of Radial Basis Functions (RBFs) $$\hat{E_{k}}$$ is used as a regression function for the input-output (i.e. xk − Y) data set (left). The optimized regression function is then evaluated at all values of xk and the variance of $$\hat{E_{k}}$$ (right) used to approximate the term $$V_{x_{k}}[E_{x_{\sim k}}(y|x_{k})]$$ in eq. (6). 3.4 PAWN sensitivity analysis PAWN (Pianosi & Wagener 2015b) is a density-based method where sensitivity is measured by estimating the variation to the output y distribution when removing the uncertainty in one or more parameters xk. This variation is calculated from the measure of distance between the unconditional (when all inputs vary simultaneously) and conditional (when all inputs vary but xk, which is set to a nominal value) CDFs. The PAWN sensitivity index for the kth input is defined as:   $$S_{k} = \max _{x_{k}}\max _{y}|F_{y}(y)-F_{y|x_k}(y|x_{k})|,$$ (7)where Fy(y) and $$F_{y|x_{k}}(y|x_{k})$$ are the unconditional and conditional CDFs of the output respectively. The inner maximum of eq. (7) defines the maximum absolute difference between the two CDFs approximated via the KS statistic using empirical distribution functions. As the KS statistic will depend on the nominal (i.e. fixed) value of xk, the outer maximum of eq. (7) extracts the maximum KS statistic over all values of xk. If the data set does not contain multiple samples with the same value of xk, that is, a generic input-output data set, conditional distributions can be conditioned on ‘similar’ values of xk. Therefore eq. (7) can be approximated as   $$S_{k} = \max _{j=1,....,n}\max _{y}|F_{y}(y)-F_{y|x_{k}}(y|x_{k} \in \alpha _{j})|,$$ (8)where αj are n (e.g. 10) equally spaced intervals over the range of variation of xk. A schematic example of this method is shown in Fig. 10. Figure 10. View largeDownload slide Red line (left image) indicates the unconditional (when all inputs vary simultaneously) CDF while shaded lines the conditional CDFs (all inputs vary but xk) when xk is fixed at incremental nominal values. The KS statistic (see caption for Fig. 8) is computed for each unconditional–condition CDF pair and the PAWN sensitivity index taken as the maximum KS value for the input xk (right). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 10. View largeDownload slide Red line (left image) indicates the unconditional (when all inputs vary simultaneously) CDF while shaded lines the conditional CDFs (all inputs vary but xk) when xk is fixed at incremental nominal values. The KS statistic (see caption for Fig. 8) is computed for each unconditional–condition CDF pair and the PAWN sensitivity index taken as the maximum KS value for the input xk (right). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. 4 EXPERIMENTAL SETUP: DEFINING MODEL INPUT FACTORS AND OUTPUT We consider a single unknown overburden chalk layer (layer 5 in Table 1) whose material properties are largely uncertain. We assume little a priori knowledge of its properties and thus large, uniform, independent uncertainty distributions for all its physical parameters. The production profile and mechanical properties of the reservoir and underburden are assumed known. In other words, known to greater degree of accuracy (i.e. much smaller uncertainty) that of the overburden chalk. Thus, their properties are kept constant and not considered in the GSA. In total, 21 different input factors are subjected to the GSA, each a material property of the overburden chalk. A summary of these parameters is presented in Table 2 along with their range of uncertainties. These uncertainty ranges chosen are based on typical properties presented in the generic material database of Rockfield Software Limited and are made as wide as possible. Most of these parameters are described in greater detail in the definition of an elastoplastic material found in Appendix B. However, to ease GSA parameter space sampling we consider the poroelastic parameters A and B (eq. B8) to be equal (for simplicity) and a minimum Poisson's ratio (υmin in eq. B9) as a ratio υratio of υmax. Also, we assume that overburden rocks behave elastically during production as production related stress changes in the overburden are generally small compared to the yield strength of the rock. Therefore, we assume those parameters that define the shape of state boundary surface to have no effect on the seismic traveltimes. However, we do vary the yield surface evolution parameters as they will affect the stress dependant porosity parameter in eq. (B8). It should also be noted that for the purpose of this study we assume all unphysical, or algorithm specific (e.g. coupling rate, mesh type and size, etc.) parameters are optimized and are not considered in the GSA. Table 2. Chalk layer physical properties and their parameter sensitivity ranges. No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Kxy  Vertical-Horizontal stress coef    0.6  0.4  1.1  2  Pp (Mpa)  Over/underpressure    20  0  +40  3  ρf (g cc−1)  Fluid density    1.02  1  1.2  4  ρg (g cc−1)  Grain density    2.71  2.6  2.8  5  λ  Cam-clay constant  (B4) and (B5)  0.06  0.02  0.1  6  κ  Cam-clay constant  (B4) and (B5)  0.008  0.002  0.012  7  Eref (Mpa)  Reference Young's Modulus  (B8)  30500  1 × 104  3.3 × 104  8  n  Elastic constant  (B8)  0.02  0.001  0.1  9  A/B  Elastic constant  (B8)  −0.2758  0  −0.5  10  c  Elastic constant  (B8)  −0.1  −0.5  −0.001  11  υmax  Max Poisson's ratio  (B9)  0.33  0.2  0.4  12  υratio  Min Poisson's ratio  (B9)  1  1  1.5  13  m  Elastic constant  (B9)  1  0.01  1  14  α  Biot constant    1  0.5  1  15  K0 (m2)  Permeability constant  (B10)  1 × 10−22  1 × 10−23  1 × 10−18  16  x  Permeability constant  (B10)  3  1  4  17  y  Permeability constant  (B10)  2  1  6  18  Ks (Mpa)  Grain stiffness    2400  2400  5000  19  Kf (Mpa)  Fluid stiffness    13 × 104  8 × 104  15 × 104  20  Φinit  Initial porosity  (B6)  0.06  0.01  0.23  21  Φref  Reference porosity  (B6)  0.3  0.3  0.5  No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Kxy  Vertical-Horizontal stress coef    0.6  0.4  1.1  2  Pp (Mpa)  Over/underpressure    20  0  +40  3  ρf (g cc−1)  Fluid density    1.02  1  1.2  4  ρg (g cc−1)  Grain density    2.71  2.6  2.8  5  λ  Cam-clay constant  (B4) and (B5)  0.06  0.02  0.1  6  κ  Cam-clay constant  (B4) and (B5)  0.008  0.002  0.012  7  Eref (Mpa)  Reference Young's Modulus  (B8)  30500  1 × 104  3.3 × 104  8  n  Elastic constant  (B8)  0.02  0.001  0.1  9  A/B  Elastic constant  (B8)  −0.2758  0  −0.5  10  c  Elastic constant  (B8)  −0.1  −0.5  −0.001  11  υmax  Max Poisson's ratio  (B9)  0.33  0.2  0.4  12  υratio  Min Poisson's ratio  (B9)  1  1  1.5  13  m  Elastic constant  (B9)  1  0.01  1  14  α  Biot constant    1  0.5  1  15  K0 (m2)  Permeability constant  (B10)  1 × 10−22  1 × 10−23  1 × 10−18  16  x  Permeability constant  (B10)  3  1  4  17  y  Permeability constant  (B10)  2  1  6  18  Ks (Mpa)  Grain stiffness    2400  2400  5000  19  Kf (Mpa)  Fluid stiffness    13 × 104  8 × 104  15 × 104  20  Φinit  Initial porosity  (B6)  0.06  0.01  0.23  21  Φref  Reference porosity  (B6)  0.3  0.3  0.5  View Large The outputs analysed by the GSA are the individual overburden layer vertical time-shifts Δtv; specifically, Δtv over each production year (i.e. 1 through to 20 yrs) at three separate locations (highlighted in Fig. 1). Sensitivity indices are generated for each parameter by analysing each Δtv output. In this study, we present global sensitivity indices by averaging combined sets of individual results. For example, the sensitivity of chalk layer 5 to a certain parameter will be the average of the three individual indices calculated at the three locations shown in Fig. 1. Therefore, the sensitivity indices take into account variations in Δtv across the model domain. 5 RESULTS Utilizing over 4000 model perturbations we conduct an in depth multimethod GSA using the SAFE Toolbox of Pianosi et al. (2015a). In our GSA, we initially use the EET on an ensemble of 1540 (i.e. n = 70 EE's per input) model runs to screen those parameters (see Table 2) that have little effect over overburden Δtv. We then create an additional random ensemble of 3000 model runs with a reduced number of variable parameters. The RSA, VBSA and PAWN techniques are then applied to this data set simultaneously to rank the most influential parameters in order of importance. We then extract the most influential parameters (hereby referred to as ‘active’ parameters) and assess their potential seismic calibration by model performance (i.e. mapping their location within the model space). 5.1 Screening model parameters Fig. 11 shows the results of the EET analysis. Taking the three vertical locations shown in Fig. 1 we calculate average EE measures for both the uncertain chalk layer and remaining overburden layers at 1 yr intervals over the total 20 yrs of production. As discussed in Section 3.1, a large (absolute) measure of central tendency (i.e. mean) indicates an input with an important direct influence on the model output. A large spread (i.e. standard deviation) indicates an input with a strong nonlinear effect on the model output. Thus, parameters which show a significant shading have a greater influence over the modelled Δtv. Figure 11. View largeDownload slide Time varying Elementary Effects (EEs) considering the resultant change in layer traveltime Δtv at yearly intervals over the total 20 yrs of production. Results are computed considering the Δtv results of the uncertain chalk layer only and the Δtv result of the other (unchanged) overburden layers at the locations specified in Fig. 1. Figure 11. View largeDownload slide Time varying Elementary Effects (EEs) considering the resultant change in layer traveltime Δtv at yearly intervals over the total 20 yrs of production. Results are computed considering the Δtv results of the uncertain chalk layer only and the Δtv result of the other (unchanged) overburden layers at the locations specified in Fig. 1. It is apparent from Fig. 11 that altering the material properties of a single layer affects Δtv across the whole overburden. This demonstrates a complex nonlinear model behaviour. However, similar sensitivity patterns emerging across all overburden layers suggests the total modelled overburden Δtv could potentially be explained by a handful of model variables. Although Fig. 11 suggest certain model parameters to be more influential than others, it is difficult to conclusively screen a number of parameters as being non-influential. However, it is apparent from these results that the Cam-Clay parameters λ and κ do not affect the modelled Δtv. Both parameters measure zero EE mean and standard deviation. This is not unexpected as they primarily affect the yield surface evolution (eqs B4 and B5) and, as their influence is negligible, it confirms our assumption that the overburden in our model remains elastic during production. However, their zero measure also suggests negligible influence over the stress-dependant porosity parameter in eq. (B8). As a result we screen out the Cam-Clay parameters λ and κ as non-influential but consider all remaining 19 parameters as potentially influential to Δtv. 5.2 Ranking model parameters To further investigate the influence of the remaining uncertain parameters on the modelled Δtv, we run an additional 3000 different model input combinations. We use the same uncertainty ranges as expressed in Table 2, but fix the Cam-Clay constants to their true value (Table 1). We create 3000 ensemble by randomly sampling the parameter space using Latin-hypercube sampling (Forrester et al. 2008). We apply the RSA, VBSA and PAWN sensitivity techniques simultaneously to this input–output data set to give complimentary parameter sensitivity indices. Fig. 12 shows the results of all three GSA methods focused on the final modelled Δtv (i.e. after 20 yrs of production). Note that here we use a similar procedure as that of the EET by summarizing sensitivity indices as average values for the uncertain chalk layer and remaining overburden layers over the locations shown in Fig. 1. Figure 12. View largeDownload slide GSA sensitivity indices of the reduced set of model parameters. Blue circles represent the RSA results, black squares PAWN and grey hollow boxes the VBSA results. Sensitivity indices are computed considering the Δtv results of the uncertain chalk layer only and the results of the remaining (unchanged) overburden layers. These results focus on the final model Δtv that is, after 20 yrs of production, at the locations highlighted in Fig. 1. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 12. View largeDownload slide GSA sensitivity indices of the reduced set of model parameters. Blue circles represent the RSA results, black squares PAWN and grey hollow boxes the VBSA results. Sensitivity indices are computed considering the Δtv results of the uncertain chalk layer only and the results of the remaining (unchanged) overburden layers. These results focus on the final model Δtv that is, after 20 yrs of production, at the locations highlighted in Fig. 1. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. A parameter with a greater sensitivity index indicates one which has a greater direct influence over Δtv. It is apparent from Fig. 12 that, although not giving exactly the same absolute measures of sensitivity, all three GSA techniques provide suitably similar global trends. We observe similar results as that of the EET analysis (Fig. 11) where the parameters that control the nonlinear elastic response (Eref, c, ϕinit or eq. B8) and the Biot coefficient α are notably the most influential across all overburden layers. It is therefore fair to assume these to be the active parameters of our model. If we plot the time-varying sensitivity indices of these parameters when considering Δtv of the chalk, we can compare the influence of the elastic properties to α over production. The results are presented in Fig. 13. Here, we see the elastic parameters to be most influential during earlier production times but they appear to be outweighed by α later in production. Figure 13. View largeDownload slide The time-varying GSA sensitivity indices of the four most influential parameters within the uncertain chalk layer. The elastic parameters are shown in blue whilst the Biot coefficient α in red. For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this paper. Figure 13. View largeDownload slide The time-varying GSA sensitivity indices of the four most influential parameters within the uncertain chalk layer. The elastic parameters are shown in blue whilst the Biot coefficient α in red. For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this paper. 5.3 Mapping data to model space We take the 3000 model ensemble and compare these simulations to the result of our original model run. Assuming the overburden Δtv of the original model (Fig. 6) to be observed data, we compare model residuals in the form of a Δtv Root Mean Square Error (RMSE). As time-lapse seismic data is time consuming and costly to acquire, shooting data at yearly intervals (or less) is generally unfeasible. Therefore, we compute the RMSE considering the Δtv results at just three production time steps. The results of the GSA (Figs 11 and 13) show that the active variables do not start to significantly affect overburden Δtv till 3 yrs of production. They also show that the Biot coefficient α takes approximately 10 yrs to become as influential as the elastic coefficients (e.g. Fig. 13). Taking this into consideration we choose to compute the RMSE using the Δtv result at 3, 5 and 10 yrs. Note that the RMSE is calculated assuming the results from all three vertical locations shown in Fig. 1. Fig. 14 show parallel coordinate plots of the best 5 per cent (i.e. 15) of models whose Δtv measurements most closely resemble that of the data (i.e. our original, truth model). These possess the lowest RMSE and are referred to as behavioural models. Also shown in Fig. 14 are their corresponding Δtv logs (plotted just at location 2 of Fig. 1) after 10 and 20 yrs of production. The results show that each behavioural model Δtv result closely resembles that of the data. Even their forward predicted Δtv values (i.e. after 20 yrs of production) are similar to those observed in the data. However, these models appear randomly scattered along the uncertain parameter range. Thus, models with significantly different active parameters each produce similar, possibly acceptable, solutions. This suggests a complex, possibly ill-posed, model space. If we simplify the objective function to the RMSE of just the uncertain chalk layers Δtv (i.e. ignoring the Δtv of other overburden layers) we see a slightly different result (Fig. 14). This optimization produces a different set of behavioural models which, as expected, do a better job of fitting the data of the uncertain chalk layer. These behavioural models also appear less scattered throughout the model space. Almost all of their active parameters more closely resemble those of the original (i.e. true) model. However, these models contain a significantly lower value of the elastic coefficient Eref (and as a result a lower pre-production Young's Modulus). Thus, significantly different model parameters still produce similar, possibly acceptable, solutions even with a simplified, condensed data space. Figure 14. View largeDownload slide Parallel coordinate plots showing the active parameters of the best 5 per cent (i.e. 15) of models whose Δtv results most closely resemble that of original, that is truth, model (Fig. 6). Also shown are the corresponding models Δtv overburden logs after both 10 and 20 yrs of production. The original model results are highlighted in red, while the closest models in black. The overburden Δtv logs of the whole model ensemble shown in grey. The model residuals were computed by taking the Root Mean Square Error (RMSE) of the whole overburden Δtv results after 3, 5 and 10 yrs of production. Also shown are the results when only the uncertain chalk layers Δtv results are considered in the residual calculation. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 14. View largeDownload slide Parallel coordinate plots showing the active parameters of the best 5 per cent (i.e. 15) of models whose Δtv results most closely resemble that of original, that is truth, model (Fig. 6). Also shown are the corresponding models Δtv overburden logs after both 10 and 20 yrs of production. The original model results are highlighted in red, while the closest models in black. The overburden Δtv logs of the whole model ensemble shown in grey. The model residuals were computed by taking the Root Mean Square Error (RMSE) of the whole overburden Δtv results after 3, 5 and 10 yrs of production. Also shown are the results when only the uncertain chalk layers Δtv results are considered in the residual calculation. For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Fig. 14 along with the sensitivity results of the GSA emphasizes the complexity and nonlinearity of the model behaviour. The uncertainty in a single overburden layer's material properties appears to affect the predicted Δtv across the whole overburden domain. Due to significant pore pressure changes, the reservoir undergoes far more extreme mechanical changes than the overburden during production. Thus, it is reasonable to presume that the uncertainty in the mechanical properties of the reservoir may, more heavily, influence overburden Δtv than similar uncertainty in overburden properties. To test this hypothesis, we assume the mechanical properties of our reservoir to be uncertain, whilst holding the properties of the overburden constant (Table 1). We take the four active parameters highlighted by the GSA and assume uniform, independent uncertainty distributions for their values in the reservoir. A summary of their uncertainty ranges is outlined in Table 3. We assume large, uniform, independant uncertainty distributions for each parameter, similar to that used for the overburden chalk in the GSA (Table 2). Using the same Latin Hypercube sampling as that used in the GSA (Forrester et al. 2008), we run 200 different model input combinations and compute their overburden Δtv results. The Δtv values are again computed at the same locations as that used in the GSA (Fig. 1). Table 3. Active parameters of the reservoir layer and their uncertain sensitivity ranges. No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Eref (Mpa)  Reference Young's Modulus  (B8)  225  100  1000  2  c  Elastic constant  (B8)  −1.28  −1.5  −0.001  3  α  Biot constant    1  0.5  1  4  Φinit  Initial porosity  (B6)  0.12  0.01  0.23  No.  Parameter  Nomenclature  Equation No  Truth value  Range            Min  Max  1  Eref (Mpa)  Reference Young's Modulus  (B8)  225  100  1000  2  c  Elastic constant  (B8)  −1.28  −1.5  −0.001  3  α  Biot constant    1  0.5  1  4  Φinit  Initial porosity  (B6)  0.12  0.01  0.23  View Large Shown in Fig. 15 are the corresponding overburden logs for each model run after 10 yrs of production. The results confirm the significant affect reservoir uncertainty has on overburden Δtv. Its uncertainty appears to influence overburden Δtv to a much greater extent than suitable similar uncertainty in overburden layers (Fig. 14). It is important to note that the extreme mechanical changes seen in the reservoir could result in plastic deformation. Thus, parameters which govern the rocks state boundary surface could potentially have significant influence. These parameters should be included if an in depth sensitivity study is to be undertaken for the reservoir region. Figure 15. View largeDownload slide Overburden Δtv logs of the uncertain reservoir ensemble (Table 3) after 10 yrs of production shown via grey lines. The, truth model results are shown in red (Fig. 6) while the dotted lines represent the extreme values seen within the original GSA results (Fig. 14). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. Figure 15. View largeDownload slide Overburden Δtv logs of the uncertain reservoir ensemble (Table 3) after 10 yrs of production shown via grey lines. The, truth model results are shown in red (Fig. 6) while the dotted lines represent the extreme values seen within the original GSA results (Fig. 14). For interpretation of the references to colour in this figure, the reader is referred to the web version of this paper. 6 DISCUSSION The results of the GSA highlighted that out of an initial 21 model parameters, the modelled overburden Δtv is most heavily influenced by just 4 (e.g. Figs 11 and 12). The Biot coefficient (α) and the parameters that govern a materials elastic behaviour (Eref, c and ϕinit). The GSA also shows that these variables take 3 yrs (from the start of production) to become significantly influential and that α takes approximately 10 yrs to become as influential as the elastic parameters (e.g. Fig. 13). This time varying sensitivity demonstrates how acquiring multiple vintages of time-lapse seismic data could be advantageous. However, these results also suggest that significant overburden near-offset seismic time-shifts may take time to manifest. Thus, time-lapse seismic data taken shortly after the start of production may not be conclusive enough to aid the advanced numerical calibration of geomechanical models. However, acquiring early seismic data can be beneficial as an early warning system. It can highlight any large discrepancies between model and reality, and suggest the potential need for additional data (e.g. well log or core data). The elastic parameters being influential is not totally unexpected, as seismic traveltimes are affected by changes in stress and path length and thus governed by rock stiffness (i.e. Young's Modulus, eq. B8). However, slightly unexpected is the significant influence of the Biot coefficient α. Typically, the overburden is modelled as an undrained scenario (i.e. no fluid flow) in which you assume that there is no production related pore pressure change. Therefore, you would expect α (i.e. σ΄ = σ − αPp) to have little influence over changes in effective stress and hence Δtv. However, ELFEN implicitly evaluates the pore pressure Pp of the whole model domain (i.e. whole domain coupling) as a function of α and the volumetric strain εv. Therefore, although there is no (or little) fluid flow outside of the reservoir, the pore pressure is affected by mechanical changes in volumetric strain. The GSA demonstrates this and emphasizes the importance of modelling fluid flow and pore pressure outside of the reservoir. Their slight instabilities can cause non-negligible effects to the model output. The consequence of the material properties of a single layer affecting Δtv across the whole overburden demonstrates a complex nonlinear model behaviour. Thus, analysing model activity globally (i.e. across whole modelled domain), as opposed to locally, could be crucial for potential calibration. Although we highlight just four active parameters for the Δtv results, it is also influenced, albeit to a lesser extent, by the remaining uncertainty in other parameters (e.g. Fig. 12). Therefore it must be stressed that changes in overburden Δtv are not entirely determined by changes in these four active variables. Thus, calibration procedures focussing on a condensed model space should also account for variations to model output caused by changes in less sensitive variables. Models with significantly different input parameters produce similar Δtv results (e.g. Fig. 14). This is the case when considering both global Δtv results (i.e. across the whole overburden) and when focused to local model output (i.e. results of uncertain chalk layer only). This highlights the complexity of the model space where a single global solution will most likely not exist. Instead, numerous models, of different input combinations will produce equally acceptable solutions (e.g. Fig. 14). In this study, we assume only the uncertainty of a single overburden layer. The complexity of the model space will undoubtedly increase when we consider the uncertainty in the mechanical properties of other layers. The GSA results do suggest that time-lapse seismic data could potentially be able to distinguish between certain models within our ensemble. The difference in Δtv between certain models being on the order of 2–3 ms (e.g. Fig. 14). However, the size and complexity of the model space makes calibration via Monte Carlo based approaches unfeasible. For example, it is clear from the results of the GSA that 3000 randomly sampled model perturbations are insufficient to determine definitive solutions (e.g. Fig. 14). However, with a more sophisticated method of sampling of the model space, we may be able to determine the most plausible solutions from a more feasible number of model runs. Although still unlikely to determine a single unique solution, you may potentially be able to determine a set of potential scenarios within sensible time frames. This could provide vital additional information for use in conjunction with qualitative analysis. We find that the uncertainty in the mechanical properties of the reservoir heavily influence overburden Δtv (e.g. Fig. 15). Its uncertainty appearing to influence Δtv to a much greater extent than suitable similar uncertainty in the overburden (e.g. Fig. 14). This demonstrates the nonlinearity of the model behaviour and the importance of a suitably accurate reservoir model. Time-lapse seismic calibration of other properties of a geomechanical model will thus only be possible once the reservoir behaviour is known to a suitable degree of accuracy. It is important to note that the sensitivity measurements of the GSA are heavily affected by the uncertainty in model parameters. For example, considering a much smaller uncertainty range in the elastic coefficients would result in their sensitivity being significantly less than seen in this study. Therefore it is always important to cross analyse the results of the GSA with the uncertainty range used. It is also important to highlight that we concentrate our analysis on near-offset that is, vertical traveltimes. If we consider time-shift versus offset estimations we would expect an increased sensitivity in parameters such as the Poisson's ratio and the horizontal-to-vertical stress ratio Kxy. It is also important to mention that in this study we do not include the uncertainty in the stress-dependant rock physics model. We assume the in situ seismic velocity of each rock, and its stress-dependence, is known (i.e. no uncertainty). In reality, this process is highly uncertain, especially in overburden rocks where stress–velocity core data is often not available (e.g. Price et al. 2016). This complex modelling step is beyond the scope of this study. However, any future quantitative calibration of geomechanical models with near offset time-lapse seismic data must include this uncertainty in its calculations. Finally, in this study we have not accounted for random modelling errors (e.g. variations to the implicit and explicit solutions caused by the parametrization of their solvers), which can act as noise or bias to the resulting output distributions. However, since large distributions are seen in Δtv (e.g. Figs 14 and 15) it is safe to assume these random modelling errors to be insignificant. 7 CONCLUSION In this study, we have analysed the potential of near-offset time-lapse seismic data to aid the geomechanical calibration of the overburden. We build a geomechanical model of a typical North Sea HPHT reservoir and utilize over 4000 model perturbations to conduct an in-depth multimethod GSA. We outline that out an initially large set of uncertain material properties, the modelled overburden Δtv are mainly affected by just four ‘active’ parameters. These being the effective stress (i.e. Biot) coefficient (α) and the parameters that govern the material's elastic behaviour (i.e. stiffness). The influence of the Biot coefficient highlights the importance of modelling fluid flow and pore pressure outside of the reservoir. These results show that the model space can be condensed to these active variables for calibration. However, the variations caused by less sensitive variables are non-negligible and so should also be taken into account. We demonstrate how the FE model is complex and highly nonlinear. Altering the material properties of a single layer affects the Δtv results of the whole overburden domain. As a result, constraining the uncertainty in model parameters appears challenging. Multiple combinations of model parameters can yield equally possible model realizations and hence Δtv results. Consequently, calibration via a large number of random model perturbations is unfeasible. However, the significant difference in the Δtv estimates of certain models within our ensemble suggests more sophisticated calibration methods could potentially be feasible. Determining a set of potential solutions using a small number of model runs could be possible given intelligent sampling of the model space. The results of the time-varying GSA suggest that significant overburden time-shifts may take time to manifest. Thus, time-lapse seismic data taken shortly after the start of production may not be conclusive enough to aid the advanced numerical calibration of geomechanical models. We find that similar uncertainty in reservoir mechanical properties (as to overburden properties) appear to influence overburden Δtv to a much greater extent. This demonstrates the complex nonlinear model behaviour and stresses the importance of reservoir characterization. Thus, to calibrate parts of the geomechanical model other than the reservoir, such as the overburden, the reservoir behaviour must be known to a suitable degree of accuracy. This lends credit to small scale sensitivity studies of reservoir uncertainty before overburden calibration is considered. Although this particular study focuses on the modelling software ELFEN. The findings and conclusions can be related to all FE geomechanical modelling software. We stress that, although not straightforward, there is potential to quantitatively calibrate geomechanical models via time-lapse seismic data. Further studies like this could extend models to include anisotropic behaviour and fault properties. Acknowledgements DP is supported by a scholarship funded by Total GRC and the School of Earth and Environment, University of Leeds. The author would like to thank Francesca Pianosa of Department of Civil Engineering, University of Bristol for providing a copy of the SAFE toolbox and for her advice on the use of different sensitivity techniques. The author would also like to thank Rockfield Ltd for providing the University of Leeds with an academic license for their modelling software ELFEN. DA acknowledges the Research Council UK (EP/K035878/1; EP/K021869/1; NE/L000423/1) for financial support. REFERENCES Alerini M., Ayzenberg M., Ek T., Feng T., Hustoft L., Lie E., Liu S., Skjei N., Skjervheim J.A., 2014. Utilization of Time-lapse seismic for reservoir model conditioning, in76th EAGE Conference and Exhibition 2014  EAGE. Angus D.A. et al.  , 2015. Integrated hydro-mechanical and seismic modelling of the Valhall reservoir: a case study of predicting subsidence, AVOA and microseismicity, Geomech. Energy Environ.  2 32– 44. Google Scholar CrossRef Search ADS   Bloch S., Lander R.H., Bonnell L., 2002. Anomalously high porosity and permeability in deeply buried sandstone reservoirs: origin and predictability, AAPG Bull.  86( 2), 301– 328. Bruno M.S., 2002. Geomechanics and decision analysis for migration compaction related casing damage, SPE drilling and completion, SPE  79519 179– 188. Calvert R., 2005. Insights and Methods for 4D Reservoir Monitoring and Characterization  Distinguished Instructor Series No. 8, EAGE. Google Scholar CrossRef Search ADS   Carman P.C., 1937. Fluid flow through a granular bed, Trans. Inst. Chem. Eng.  15 150– 156. Crook A.J.L., Willson S.M., Yu J.G., Owen D.R.J., 2006. Predictive modelling of structure evolution in sandbox experiments, J. Struct. Geol.  28( 5), 729– 744. Google Scholar CrossRef Search ADS   De Gennaro S., Onaisi A., Grandi A., Ben-Brahim L., Neillo V., 2008. 4D reservoir geomechanics: a case study from the HP/HT reservoirs of the Elgin and Franklin fields, First Break  26 53– 59. Google Scholar CrossRef Search ADS   De Gennaro S., Schutjens P., Frumau M., Fuery M., Ita J., Fokker P., 2010. The role of geomechanics in the development of an HPHT field, in 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium  Salt Lake City, UT, American Rock Mechanics Association, pp. 10– 450. Emerick A.A., Reynolds A.C., 2012. History matching time-lapse seismic data using the ensemble Kalman filter with multiple data assimilations, Comput. Geosci.  16( 3), 639– 659. Google Scholar CrossRef Search ADS   Fjær E., Kristiansen T.G., 2009. An Integrated Geomechanics, Rock Physics and Seismic model, in 71st EAGE Conference & Exhibition  Amsterdam. Forrester A., Sobester A., Keane A., 2008. Engineering Design via Surrogate modelling: A Practical Guide  John Wiley & Sons. Garcia A.G., 2011. Dynamic reservoir characterization from overburden time-lapse strains, PhD thesis  Institute of Petroleum Engineering, Heriot-Watt University. Gonzalez-Carballo A., Guyonnet P.Y., Levallois B., Veillerette A., Deboiasne R., 2006. 4D monitoring in Angola and its impact on reservoir understanding and economics, Leading Edge  25 1150– 1159. Google Scholar CrossRef Search ADS   Gosselin O. et al.  , 2003. History matching using time-lapse seismic, in SPE Annual Technical Conference and Exhibition  Society of Petroleum Engineers. Guilbot J., Smith B., 2002. 4-D constrained depth conversion for reservoir compaction estimation: Application to Ekofisk Field, Leading Edge  21 302– 308. Google Scholar CrossRef Search ADS   Hatchell P.J., Bourne S.J., 2005. Rocks under strain: Strain-induced time-lapse time shifts are observed for depleting reservoirs, leading Edge  24( 12), 1222– 1225. Google Scholar CrossRef Search ADS   Hawkins K. et al.  , 2007. Production-induced stresses from time-lapse time shifts: A geomechanical case study from Franklin and Elgin fields, Leading Edge  26 655– 662. Google Scholar CrossRef Search ADS   Herwanger J., Horne S.A., 2009. Linking reservoir geomechanics and time-lapse seismics: predicting anisotropic velocity changes and seismic attributes, Geophysics  74 W13– W33. Google Scholar CrossRef Search ADS   Jackson D., Richardson M., 2007. High-frequency Seafloor Acoustics  Springer Science & Business Media. Google Scholar CrossRef Search ADS   Jaeger J.C., Cook N.G., Zimmerman R., 2009. Fundamentals of Rock Mechanics  John Wiley & Sons. Japsen P., 1998. Regional velocity-depth anomalie, North Sea chalk: a record of overporessure and neogene uplift and erosion, AAPG Bull  82( 11), 2031– 2074. Japsen P., 2000. Investigation of multi-phase erosion using reconstructed shale trends based on sonic data, Sole Pit axis, North Sea, Glob. Planet. Change  24( 3), 189– 210. Google Scholar CrossRef Search ADS   Jones J.C., 2010. Hydrocarbons-physical Properties and Their Relevance to Utilisation  BookBoon, pp. 1– 100. Kozeny J., 1927. Über kapillare leitung des wassers im boden, Sitzungsber. Akad. Wiss. Wien  136 271– 306. Kvendseth S.S., 1988. Giant Discovery, A history of Ekofisk through the First 20 Years  Philips Petroleum Company. Kudarova A.M., Macbeth C., 2016. Offset-dependence of production-related 4D time shifts: Real data examples and modelling, in 87th SEG Annual Meeting  Dallas, TX. Landrø M., 2001. Discrimination between pressure and fluid saturation changes from time-lapse seismic data, Geophysics  66 836– 844. Google Scholar CrossRef Search ADS   Mallon A.J., Swarbrick R.E., 2008. Diagenetic characteristics of low permeability, non-reservoir chalks from the Central North Sea, Mar. Pet. Geol.  25( 10), 1097– 1108. Google Scholar CrossRef Search ADS   Marcussen ø., Maast T.E., Mondol N.H., Jahren J., Bjørlykke K., 2010. Changes in physical properties of a reservoir sandstone as a function of burial depth—the Etive formation, northern North Sea, Mar. Pet. Geol.  27( 8), 1725– 1735. Google Scholar CrossRef Search ADS   Minkoff S.E., Stone C.M., Bryant S., Peszynska M., Wheeler M.F., 2003. Coupled fluid flow and geomechanical deformation modeling, J. Pet. Sci. Eng.  38 37– 56. Google Scholar CrossRef Search ADS   Morris M.D., 1991. Factorial sampling plans for preliminary computational experiments, Technometrics  33( 2), 161– 174. Google Scholar CrossRef Search ADS   Okiongbo K.S., 2011. Petrophysical properties of North Sea shales, Res. J. Appl. Sci. Eng. Technol.  3( 1), 46– 52. Petropoulos G., Srivastava P.K., 2016. Sensitivity Analysis in Earth Observation Modelling  Elsevier Science & Technology Books. Pianosi F., Wagener T., 2015. A simple and efficient method for global sensitivity analysis based on cumulative distribution functions, Environ. Modelling Softw.  67 1– 11. Google Scholar CrossRef Search ADS   Pianosi F., Sarrazin F., Wagener T., 2015. A Matlab toolbox for Global Sensitivity Analysis, Environ. Modelling Softw.  70 80– 85. Google Scholar CrossRef Search ADS   Price D.C., Angus D.A., Garcia A., Fisher Q.J., 2016. Probabilistic analysis and comparison of stress-dependent rock physics models, Geophys. J. Int.  210( 1), 196– 209. Google Scholar CrossRef Search ADS   Rao S.S., 2010. The Finite Element Method in Engineering  Elsevier. Rockfield Software Limited, 2012. ELFEN GeoDB Generic Materials  Swansea. Saltelli A., Annoni P., Azzini I., Campolongo F., Ratto M., 2010. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, Comput. Phys. Commun.  181( 2), 259– 270. Google Scholar CrossRef Search ADS   Samier P., Onaisi A., Fontaine G., 2003. Coupled analysis of geomechanics and fluid flow in reservoir simulation, SPE Reservoir simulation symposium  Houston, USA, 3–5 Febuary, SPE 79698. Schneider F., Potdevin J.L., Wolf S., Faille I., 1996. Mechanical and chemical compaction model for sedimentary basin simulators, Tectonophysics  263( 1), 307– 317. Google Scholar CrossRef Search ADS   Slagstad T., Barróre C., Davidsen B., Ramstad R.K., 2008. Petrophysical and thermal properties of pre-Devonian basement rocks on the Norwegian continental margin, Geol. Surv. Norway Bull.  448 1– 6. Sobol I., 1990. Sensitivity analysis for non-linear mathematical models. Math. Modelling Comput. Exp.  1 407– 424. Spear R.C., Hornberger G.M., 1980. Eutrophication in peel inlet II. Identification of critical uncertainties via generalized sensitivity analysis, Water Res.  14( 1), 43– 49. Google Scholar CrossRef Search ADS   Staples R., Stevens T., Leoutre E., Jolley S., Marshall J., 2005. 4D seismic history matching-the reality, in 67th EAGE Conference and Exhibition  Madrid. Staples R., Ita J., Nash R., Hague P., Burrell R., 2007. Using 4D seismic data and geomechanical modelling to understand pressure depletion in HPHT fields of the Central North Sea, in 69th EAGE Conference and Exhibition  EAGE. Tarantola S., Giglioli N., Jesinghaus J., Saltelli A., 2002. Can global sensitivity analysis steer the implementation of models for environmental assessments and decision-making?, Stoch. Environ. Res. Risk Assess.  16( 1), 63– 76. Google Scholar CrossRef Search ADS   Tarantola A., 2005. Inverse Problem Theorey and Methods for Model Parameter Estimation  SIAM. Google Scholar CrossRef Search ADS   Vudovich A., Chin L.V., Morgan D.R., 1989. Casing deformation in Ekofisk, J. Pet. Technol.  41 729– 734. Google Scholar CrossRef Search ADS   Wright W.A., 1967. Prediction of bulk moduli and pressure-volume-temperature data for petroleum oils, ASLE Trans.  10( 4), 349– 356. Google Scholar CrossRef Search ADS   Zang J., Reederm R.J., 1999. Comparative compressibilities of calcite-structure carbonates: deviations from empirical relations, Am. Mineral.  84( 5–6), 861– 870. Google Scholar CrossRef Search ADS   APPENDIX A: COUPLING Hydro-mechanical modelling in ELFEN involves iterative coupling. Solutions for the mechanical domain are calculated explicitly whereas the fluid domain is solved implicitly via a nonlinear Newton–Raphson approach (e.g. Tarantola 2005). Solutions for both fields are solved simultaneously via transferring the data between the two domains at specified time intervals. The coupling rate for the hydro-mechanical process is similar to that of the implicit time step. This is typically 400–500 times larger than the explicit step. A robust explicit time step Δt can be estimated via:   $$\Delta t = f_{{\rm crit}} \times {\rm min} \Bigg | l^{e} \sqrt{\frac{\rho ^{e}}{E^{e}}} \Bigg |.$$ A1$$l^{e} \sqrt{\rho ^{e}/E^{e}}$$ is defined as the critical time step where Ee, ρe and le are the Young's Modulus, density and characteristic length of the minimum element, respectively. The critical time step is a stability condition which prevents the magnification of round-off errors caused by the explicit scheme (e.g. Rao 2010). Its calculation is approximate and therefore we introduce the factor of critical step fcrit. Rockfield suggests that for 2-D hydro-mechanical models fcrit = 0.9. Using the approximated mechanical time step, we choose a coupling rate (i.e. the rate at which the explicit mechanical and implicit fluid field are coupled) of 0.01 time steps. This is well within the stability threshold for our model while not being too small such that it compromises computational runtime or numerical stability. APPENDIX B: ELASTOPLASTIC MATERIAL PROPERTIES B1 Basic properties We assume all non-reservoir rocks to have a fluid density ρf of 1.02 g cc−1 (Japsen 1998), while the reservoir rock has a fluid density of 0.81 g cc−1 to represent the density of light crude oil in the North Sea (Jones 2010). The grain densities ρg of each shale is assumed to be 2.69 g cc−1 (Okiongbo 2011), chalk assumed to be purely calcite with a grain density of 2.71 g cc−1 (Japsen 1998), each sandstone to be composed mainly of quartz, 2.65 g cc−1, and the stiff underburden a grain density of 2.81 g cc−1 typical of dolostone. B2 State boundary surface The state boundary surface used for each material is the Soft Rock (SR3) model of Crook et al. (2006). The SR3 yield function is a smooth, three-invariant surface that intersects the hydrostatic axis in both tension pt0 and compression pc0. It is defined as   $$\Phi (p)=g(\theta , p)q+(p-p_{t0})\tan \beta \bigg (\frac{p-p_{c0}}{p_{t}-p_{c0}}\bigg )^{1/n},$$ B1where p and q are the effective mean stress and deviatoric stress respectively, β and n material constants and θ the Lode angle. Finally, g(θ, p) is the deviatoric plane correction term that controls the shape of the yield surface in the deviatoric plane. The evolution of the plastic flow is defined by a non-associated flow rule:   $$\dot{\varepsilon ^{p}} = \dot{\lambda } \frac{\partial \Psi }{\partial \sigma },$$ B2where $$\dot{\lambda }$$ is the plastic multiplier and Ψ is the plastic potential defined as:   $$\Psi (p)=g(\theta , p)q+(p-p_{t0}){\rm tan}\psi \bigg (\frac{p-p_{c0}}{p_{t0}-p_{c0}}\bigg )^{1/n}.$$ B3Note that eq. (B3) is of identical form to that of the state boundary surface defined in eq. (B1). However, the plastic potential is defined in terms of the angle tanψ, where ψ is the dilation parameter controlling the shape of the plastic potential surface. The deviatoric plane correction term g is scaled to be 1 such that the strength in triaxial compression directly corresponds to the strength calibrated using compressive triaxial (CTC) tests (Crook et al. 2006). The initial state boundary surface is defined at a reference porosity ϕref (i.e. surface conditions). To define the state boundary surface for each material we use a compilation of test data presented in the generic material database of Rockfield Software Limited. A summary of the final chosen parameter values presented in Table 1. B3 Yield surface evolution The evolution (or hardening) of the primary yield surface is determined through relationships that define pc and pt as a function of volumetric plastic strain $$\varepsilon _{v}^{p}$$:   $$p_{c} = p_{c0} + (p_{c0}-p_{c({\rm resid})})\bigg [{\rm exp}\bigg (-\frac{v\Delta \varepsilon _{v}^{p}}{(\lambda -\kappa )}\bigg )-1\bigg ],$$ B4  \begin{eqnarray} p_{t}^* &=& p_{t0} + (p_{t0}-p_{t({\rm resid})})\bigg [{\rm exp}\bigg (-\frac{v\Delta \varepsilon _{v}^{p}}{(\lambda -\kappa )}\bigg )-1\bigg ],\nonumber \\ p_{t} &=& {\rm max}\left[p_{t0}, p_{t}^* \right]. \end{eqnarray} B5Here, κ and λ are Cam-Clay hardening constants and v is the specific volume. v can be related to porosity via 1/(1 − ϕ). Note that pc(resid) = pc0/100 and pt(resid) = pc0/100 to ensure the yield surface is always of finite size. The volumetric plastic strain $$\varepsilon _{v}^{p}$$ can also be defined in terms of porosity via:   $$\varepsilon _{v}^{p} = {\rm {log}}\bigg [\frac{1-\phi _{{\rm ref}}}{(1-\phi _{{\rm init}})}\bigg ].$$ B6These hardening relationships allow a material characterization defined at surface conditions, with a specific reference porosity ϕref, to be used to generate data suitable for a similar material at greater depth (subjected to compaction) with a different, initialization porosity ϕinit. The shape of the state boundary surface remains unchanged but its size is governed by the scaling of pc0 to pc(init) and pt0 to pt(init). We use the compilation of test data presented in the generic material database of Rockfield Software Limited to define κ and λ. The initialization porosity ϕinit corresponds to the porosity of the material at the start of the simulation. Its value for each material is determined through either porosity-depth relationships for North Sea rocks (Bloch et al. 2002; Mallon & Swarbrick 2008) or estimated using typical North Sea bulk density ρ values (Japsen 1998, 2000; Slagstad et al. 2008; Marcussen et al. 2010; Okiongbo 2011) via:   $$\rho =\rho _g(1-\phi _{{\rm init}})+\rho _f\phi _{{\rm init}}.$$ B7A complete summary of yield surface evolution parameter values are presented in Table 1. B4 Elastic properties Young's modulus E and Poisson's ratio υ are defined in ELFEN as empirical functions of effective mean stress $$\hat{\sigma}$$:   $$E = E_{{\rm ref}}\bigg [\frac{p+A}{B}\bigg ]^n \phi (\hat{\sigma})^{c},$$ B8  $$\upsilon = \upsilon _{{\rm min}} + (\upsilon _{{\rm max}} + \upsilon _{{\rm min}})(1-e^{m\hat{\sigma}}).$$ B9Eref is the effective Young's Modulus while υmax and υmin are the Poisson's ratio values at low and high effective mean stress $$\hat{\sigma}$$. n, c and m are material constants, while A and B are also material constants used to prevent problems near zero values of $$\hat{\sigma}$$. Finally ϕ($$\hat{\sigma}$$) is the porosity, which itself is a function of effective mean stress $$\hat{\sigma}$$. These elastic relationships are calibrated such that the Young's Modulus and Poisson's ratio of each material, prior to simulation, resemble those stated in Garcia 2011 and their nonlinear behaviour calibrated using the generic material database of Rockfield Software Limited, 2012. Again, all elastic parameters can be found in Table 1. Note that we have initially specified a constant Poisson's ratio for this study as we assume negligible changes to Poisson's ratio with effective mean stress. B5 Porous flow properties The porosity dependent intrinsic permeability Kin is based upon the Kozeny–Carman (Kozeny 1927; Carman 1937) relationship   $$K_{{\rm in}}(\phi ) = K_{0} \frac{\phi ^x}{(1-\phi )^y},$$ B10where K0, x and y are material constants. The chosen parameters are provided in Table 1 and based on the values given by Schneider et al. (1996). B6 Consolidation properties Additional material consolidation properties are also required by ELFEN. An effective stress coefficient (i.e. Biot) α of 1 is chosen for all rock types (a parameter usually set to 1 in most modelling scenarios). A fluid viscosity η is also required and set to 1×10−9 MPa   s−1 (typical of pure water) in all non-reservoir layers. Within the reservoir we use η = 0.638 × 10−9 MPa s−1, similar to the viscosity of hydrocarbons from the North Sea Forties field (Jones 2010). Fluid Kf and grain Ks stiffnesses are also required. For these we use typical values found in North Sea literature (Wright 1967; Zang & Reeder 1999; Jackson & Richardsons 2007; Jaeger et al. 2009). Finally a horizontal to vertical stress ratio Kxy is prescribed at 0.6 for all materials. In addition to the required consolidation properties, we include an over/under-pressure parameter Pp. This alters the settled in situ (post-geostatic initialization) pore pressure state of each material. In addition to the overpressured reservoir region (e.g. Fig. 2), we overpressure two chalk layers in the overburden, one by +20 MPa and another by +5 MPa. A summary of all the consolidation properties can be found in Table 1. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations