## Abstract

An efficient, numerical local upscaling technique for estimating elastic geomechanical properties in heterogeneous continua is proposed. The upscaled anisotropic elastic properties are solved locally with various boundary conditions and reproduce the anisotropic geomechanical response of fine-scale simulations of sand–shale sequence models with horizontal and inclined shale bedding planes. The algorithm is automated in a parallel program and can be used to determine optimum upscaling ratios in different regions of the reservoir. The successful application of the proposed upscaling method in a field-scale coupled reservoir–geomechanics simulation demonstrates an improvement in overall computational efficiency while maintaining accuracy in the geomechanical response and reservoir performance.

Heterogeneity in a geological medium is observed at all length scales from nanometres for pores in shales to centimetres for facies in core samples to kilometres for lithologies imaged in seismic data. Geological models developed for hydrocarbon reservoirs that capture this range of heterogeneity are typically created at the centimetre scale to incorporate features that will influence fluid-flow predictions. This results in fine-scale geological models that are computationally expensive for reservoir simulations, but established workflows have been established for upscaling geological models for flow simulation studies (Deutsch 1989; Durlofsky 1991; Pickup *et al.* 1994; Zhang & Okuno 2015). With growing recognition over the last decade that geomechanical processes in reservoirs can impact fluid flow, cap-rock integrity and well integrity, geomechanical behaviour is increasingly integrated into reservoir simulation studies. These reservoir–geomechanics simulations impose significant computational challenges and require similar attention to upscaling fine-scale geological models for geomechanical modelling. However, the development of efficient and robust numerical geomechanical upscaling techniques for reservoir-scale geomechanical simulations has received less attention than flow-based upscaling. Anisotropy in rock properties and the resulting geomechanical response is a significant complicating factor in the context of heterogeneous continua where existing techniques poorly predict upscaled anisotropic properties. The ultimate geomechanical upscaling workflow should be able to consider anisotropic elastic responses, elastoplastic behaviours and failure mechanisms in heterogeneous reservoirs. The proposed methodology focuses on the upscaling of elastic geomechanical properties (Young's modulus, Poisson's ratio and shear modulus), which control the stress–strain response before failure, and represents the initial step in the development of a fully integrated reservoir–geomechanics upscaling workflow.

Elastic upscaling techniques can be classified into three categories targeting different scales and facies distribution: randomly distributed, layering and spatially correlated systems. Randomly distributed models (Fig. 1a) are widely used to represent the distribution of heterogeneities in a representative element volume (REV). REV scales are chosen sufficiently large to ensure a statistically relevant distribution for the microscopic heterogeneities, but small enough to be considered as a material point (local equilibrium) at the reservoir scale. Voigt (1889) and Reuss (1929) bounds are widely accepted as the upper and lower boundaries of the effective elastic responses of an aggregate sample of heterogeneous materials. However, Voigt and Reuss bounds assume a uniform strain or stress field, respectively, which are rare in lithologically heterogeneous material. Hashin & Shtrikman (1962, 1963) developed improved bounds based on variational principles of elasticity because it is not trivial to find an exact analytical solution for the elastic modulus of heterogeneous solids that match with experimental results. In the derivations of Hashin–Shtrikman bounds based on the principle of minimum potential energy, the deformed body is assumed to be infinite and isotropic.

Evolved from the elegant Eshelby solution (Eshelby 1957) of a single isolated ellipsoidal inclusion in an infinite matrix under uniform exterior loading, dilute schemes (Withers *et al.* 1989) have been developed for homogenization of various simply shaped particles in a matrix (e.g. spheres). However, the non-interacting assumption in dilute schemes restricts its application to randomly distributed models with low-volume fractions of inclusions. The self-consistent method (Budiansky 1965; Hill 1965) has been proposed to broaden the application of dilute schemes by assuming each inclusion is embedded in an unknown homogenized medium rather than a matrix. Unfortunately, this leads to an erroneous prediction of the effective modulus in high-porosity porous media. To deal with this issue, the generalized self-consistent method (Christensen & Lo 1979) was developed which assumes that the inclusions are surrounded by a matrix in the homogenized medium. Due to the difficulty of determining an appropriate mixing scale and complicated closed-form solution, the generalized self-consistent method is not as popular as the Mori–Tanaka approach (Mori & Tanaka 1973), which considers the interactions between inclusions via a tensor sensitive to the increased volume fraction of inclusions. For spherical inclusions, the Mori–Tanaka approach matches well with the lower bound of the Hashin–Shtrikman bounds. However, neglecting the interactions between inclusions and the matrix leads to imprecise estimations for a homogenized medium with more than 30% inclusion. For more detailed comparisons between these methods developed for REV determinations of microscopic heterogeneity, readers can refer to Zohdi & Wriggers (2008), Goodarzi *et al.* (2016) and Wetzel *et al.* (2017). In summary, the previous analytical solutions can be used to estimate the effective elastic properties of randomly distributed models at the REV scale with limited volume fractions of inclusions, but the anisotropy within the homogenized continua are usually omitted.

Geological features start to dominant geomechanical behaviour as the scale increases. Depositional or sedimentation processes generally create layered geological systems and, as such, upscaling of the elastic response in layered systems (Fig. 1b) has received significant attention. For perfectly layered porous media with isotropic elastic properties and zero intrinsic energy dissipation between layers, Backus (1962) developed an analytical upscaling method to interpret seismic wave propagation. The upscaled compliance matrix with five linearly independent parameters is calculated for a transversely isotropic layered structure. Salamon (1968) extended the Backus technique to parallel layers with transversely isotropic properties, while Gerrard (1982) extended Salamon's method to layered systems with orthotropic elastic properties. Rijpsma & Zijl (1998) established an upscaling technique for imperfectly layered rocks by introducing a perturbation approximation to adjust the solution. The non-symmetrical upscaled compliance matrix matches the fine-scale results with compression and shear-boundary conditions.

Numerical upscaling methods have also been developed for reservoir-scale geomechanical simulation. Chalon *et al.* (2004) computed the lower and upper bounds of the upscaled compliance matrix for each upscaled geomechanical cell. For a simplified synthetic layered reservoir, Salamon's method matches horizontal and vertical displacements under normal and shear loading conditions. For the spatially correlated systems illustrated in Figure 1c, the exact analytical solutions of the homogenized equations cannot be easily obtained (Zijl *et al.* 2002). Zijl *et al.* (2002) extended Salamon's homogenization method to complex spatially correlated reservoir models with node-based finite-element code with displacement–stress and displacement–energy averaging techniques providing the upper and lower boundaries of the upscaled compliance matrix.

The recovery of oil and gas through thermal methods for oil sands or hydraulic fracturing in shale gas are usually associated with interactions between flow and deformation in the reservoir. Thus, coupled geomechanical and flow simulations have been applied to improve the accuracy when large deformations lead to significant changes in flow response (Tran *et al.* 2005; Kim *et al.* 2011; Cui *et al.* 2017). Correct modelling of upscaled elastic geomechanical properties increases the accuracy of flow prediction (e.g. porosity, absolute and relative permeability) because these properties are impacted by void ratio and pore structure (Touhidi-Baghini 1998). Willis (2013) combined the Backus method with cluster analysis for upscaling a 3D mechanical earth model of the Vaca Muerta Formation to assist in hydraulic fracturing design. Berryman (2011) developed an upscaling method for poroelastic and orthotropic layered media containing fluids in drained and undrained conditions analogous to Backus' work. Ita & Malekzadeh (2015) extended Berryman's method to multiporosity layered poroelastic medium combined with both drained and undrained layers by introducing an indicator for permeable and impermeable layers.

The importance of boundary conditions assigned to a geomechanical cell to compute the effective elastic modulus attracted growing attention when numerical upscaling techniques were applied to complex reservoirs in the coupled reservoir–geomechanics simulations. Couples *et al.* (2003) found that upscaled fluid-flow and geomechanical properties in the matrix + fractures + fluids system are highly non-linear with strong dependencies on each parameter. Noting that simple boundary conditions may result in erroneous upscaled properties, they suggested building libraries with realistic boundary conditions to consider changes to adjacent coarse cells. Adopted from permeability upscaling, an oversampling technique was applied by Zhang & Fu (2010) to calculate the effective elastic modulus of a broader region than the upscaled geomechanical cell for various external loading conditions. The upscaled model reasonably reproduces displacement and pressure in a 2D numerical heterogeneous ground-subsidence problem. Settari *et al.* (2013) proposed a methodology to determine a dynamic equivalent stiffness tensor for subgrid interbedded shales in the coupled simulation of a larger offshore reservoir. The equivalent stiffness tensors are calculated by a function of volume fraction of sand, Young's modulus and pressure depletion from the uniaxial deformation history, which may lead to unexpected error in complex *in situ* stress states. Yang *et al.* (2013) introduced a numerical upscaling technique for permeability and elastic stiffness tensors, and successfully applied it to a heterogeneous 2D reservoir model. Khajeh *et al.* (2012) proposed an elastic upscaling method for transversely isotropic materials by solving local compression tests in *X*, *Y* for all 2D coarse cells, improving accuracy in a 2D synthetic case when compared to traditional averages.

Here, the proposed upscaling technique, termed ‘E-upscale’, is an extension of Khajeh's upscaling technique to 3D orthotropic materials with improved local numerical tests adopting the concept of six linearly independent loading schemes from Zohdi & Wriggers (2008), which were originally designed to estimate the anisotropic elastic response caused by microscopic heterogeneity at the REV scale. An improved estimation of the effective shear modulus is obtained through local shear tests rather than estimation from Young's modulus and Poisson's ratio, as in Khajeh's method. The local numerical testing regimes demonstrate the capability in dealing with complex spatial variability and anisotropic deformations in heterogeneous reservoirs. The non-symmetry of the upscaled compliance matrix is found to be more pronounced for inclined heterolithic stratification (IHS) than layered systems. E-upscale successfully minimizes the potential error in the prediction of global compression and shear response compared to the arithmetic (Voigt), geometric and harmonic (Reuss) averages. The application of the proposed upscaling methodology in a coupled reservoir–geomechanics simulation of the steam-assisted gravity drainage (SAGD) process demonstrates its potential to significantly reduce the computational efforts in coupled simulation while honouring the deformation response of the reservoir.

## Methodology for E-upscale

### Theoretical background

The upscaled elastic compliance matrix for each coarse (upscaled) geomechanical cell is calculated based on Hooke's law. The stress–strain relationship of an anisotropic linear elastic material can be expressed as:*i* direction; *j* when stress is in *i*; and *ij*.

The normal and shear strain of the heterogeneous upscaled geomechanical cell can be written as functions of elastic properties and applied stress:

For each numerical testing regime, equations (4–9) are simplified with a user-defined stress boundary

### Local numerical testing regimes

The local numerical testing regimes were set up for each upscaled geomechanical cell, which consisted of fine-scale geological cells with heterogeneous elastic properties. Figure 2 illustrates the overall workflow for the local numerical test procedure. During production and/or injection operations, the reservoir stress changes due to pressure and temperature fluctuations, which result in variations in the stress boundary conditions for the upscaled geomechanical cell (Fig. 2a) and lead to complex deformations of each fine-scale geological cell (Fig. 2b). The effective strain tensors for the upscaled cell (Fig. 2c) can be computed from the displacements of fine-scale models at the face boundary (Fig. 2b). The upscaled compliance matrix (Fig. 2d) is solved in a series of local numerical tests given the applied stress tensors and corresponding effective strain tensors.

The E-upscale algorithm was demonstrated on a 5 × 5 × 5 upscaled region (Fig. 3c) extracted from a section of the geological model developed for the MacKay River oil sands area of NE Alberta (Fig. 3a). Each cell has unique geomechanical properties modelled by sequential Gaussian simulations. The sand zones have larger Young's modulus, Poisson's ratio and shear modulus than the shale zone (Tables 1 and 2). The orthotropic model in FLAC^{3D} (2012: Itasca Consulting Group, Inc.) was used as the constitutive model for the geomechanical simulation. The boundary conditions of the six independent testing regimes are illustrated in Figure 4. For the first test, a constant stress of 0.1 MPa was applied normal to the *YZ* plane of the upscaled cell. The other faces were free of confinement. The stress boundary can be expressed as:

Substituting equation (10) into equations (3–6) can be rearranged to solve for

The effective strains *YZ* plane after the simulation reaches the criteria to stop, which is set as the equilibrium state (unbalanced stress/applied stress <10^{−5}).

For the second and third tests, constant stress boundaries *σ _{yy}* = 0.1 MPa as the stress tensor for the second compression test, the corresponding Young's modulus and Poisson's ratio can be computed via equations (14–16):

The stress tensor for the third compression test is *σ _{zz}* = 0.1 MPa. Equations (17–19) express the determination of

The upscaled shear modulus is calculated with simple shear tests in each axis (Fig. 4d–f) for each upscaled region. For the first shear test (Fig. 4d), the stress tensor is:

The shear strains are calculated from the average shear displacements of the fine cells at the boundaries where shear stresses are applied. The shear moduli

### Reservoir-scale algorithm

An automated algorithm is required for the E-upscale workflow to upscale all necessary regions in a large-scale reservoir geomechanical simulation model. The algorithm is coded in C++ with data transformation and automation of the local numerical tests applied sequentially for each selected upscaled region. The detailed procedure for E-upscale is explained below.

Step 1: E-upscale generates the local testing module for the selected upscaled region of the fine-scale geological model. The number of fine cells in each upscaled coarse block can be customized in

*X*,*Y*and*Z*, making it possible to optimize the upscaling ratio for different areas of the reservoir (e.g. finer mesh near wells and coarser mesh for the under-/side-burden of the reservoir).Step 2: six local tests are performed sequentially with the necessary stress boundaries (Fig. 4) at user-defined stress levels or an incremental strain boundary. The displacements of each fine cell in the non-uniformly deformed body (Fig. 2b) is computed and recorded at the boundary of the upscaled block (Fig. 2c).

Step 3: the effective strain tensors of the upscaled coarse block are calculated as the average of the fine cells at the stress boundary.

Step 4: the upscaled anisotropic elastic properties for each coarse block are calculated based on the given stress boundaries and corresponding strain tensors obtained in Step 3.

Step 5: the outputs of upscaled elastic properties are organized in an ASCII format file compatible with most geological modelling software.

The upscaling ratio is defined as the number of geological cells that are merged into the coarse geomechanical cell. For the 5 × 5 × 5 upscaled region (Fig. 3c), the upscaling ratio is 125:1. To further improve the efficiency of E-upscale in full-field geological models, the upscaled regions can be split into parallel groups and the numerical tests run simultaneously on different processors.

## Numerical validation of E-upscale

### Model description

A field-scale geological model was built for the Alberta MacKay River oil sands area in SKUA GOCAD (Paradigm). Log and core data were extracted from geoSCOUT software (geoLOGIC Systems). The geological model has five formations, designated Till, Clearwater, Wabiskaw, McMurray and Devonian, which were arranged, respectively, from shallow to deep. A series of numerical geomechanical tests with different boundary conditions were conducted for the 20 × 20 × 20 cell model (Model 1) extracted from the reservoir (Fig. 3b). The geomechanical elastic properties for sand and shale within this model are provided in Table 1, where the variations in density and elastic properties are caused by different effective stress at different depths. The input elastic properties are assumed isotropic, but it will be shown that anisotropic effective elastic properties are caused by subgrid scale heterogeneity of the upscaled region.

The shale layers in Model 1 are almost horizontal. Through intensive geological studies of the Alberta oil sands, Thomas *et al.* (1987) demonstrated the importance of inclined heterolithic stratification (IHS) lithosomes. Consequently, a second model (Model 2) was generated where the shale layers are inclined at 30° (Fig. 5) to demonstrate the difference between a transversely isotropic system (Model 1) and a fully orthotropic system (Model 2). The geomechanical properties for Model 2 are provided in Table 2.

### Traditional analytical bounds of the elastic modulus

To model the effective elastic properties of a heterogeneous coarse region consisting of two lithologies (e.g. sand and shale), the elastic properties, the volume fraction and the spatial correlation of each facies are required. At the REV scale, microscopic heterogeneities are difficult to measure. Thus, Voigt and Reuss bounds (equation 21) are powerful tools to estimate the upper and lower limits of the effective elastic modulus following a power-law average of sand and shale:*M*_{s} and *M*_{sh} are the modulus of sand and shale, respectively; *V*_{s} and *V*_{sh} are the volume fraction of sand and shale, respectively; and *ω* is a constant determined by the geometry of facies and stress–strain distribution, and ranges from −1 to 1.

Voigt assumed a uniform strain field parallel to the displacement field leading to

The geometric average is an intermediate power-law average between the arithmetic and harmonic averages often used in the upscaling of reservoir properties without apparent geological features. Subsurface reservoirs usually have complex depositional environments for different formations. However, the Voigt, Reuss, VRH bound and geometric average assume an isotropic elastic response within the upscaled block. It will be shown subsequently that this is not entirely appropriate for conditions where the upscaled block is highly heterogeneous.

### Anisotropic E-upscaled properties

The proposed E-upscale methodology is designed to reproduce the anisotropic elastic response of an upscaled block under various boundary conditions. Following the algorithm of E-upscale, each 5 × 5 × 5 region in models 1 and 2 is upscaled into a single cell, resulting in 8000 fine-scale cells in each model being upscaled to 64 coarse blocks in the upscaled model. The sand–shale heterogeneity within each coarse block, measured by *V*_{sh}, is used as a metric to characterize the development of anisotropic elastic properties in the upscaled model. For both Model 1 and Model 2, Figure 6 illustrates how the development of anisotropy evolves as *V*_{sh} varies within the coarse blocks. For each model, notable anisotropic behaviour initiates at *V*_{sh} values of <10% and becomes increasingly anisotropic with increasing *V*_{sh}. The upscaled elastic properties for Model 1 comprised of horizontal shale layers becomes transversely isotropic (Fig. 6a–c) and consistent with the work of Backus (1962); while the inclined shale layers of Model 2 produce fully orthotropic elastic properties (Fig. 6d–f).

For Model 1, the arithmetic average (Voigt bound) and harmonic average (Reuss bound) are observed as rigorous upper and lower bounds on E-upscaled Young's and shear modulus. *E _{x}* and

*E*are close to the Voigt bound as the numerical testing regime for E-upscale (loading is parallel with shale layers) is similar to Voigt's iso-strain model.

_{y}*E*is closer to the Reuss bound as the compression direction is perpendicular to the shale layers, as in the Reuss iso-stress model (Fig. 6a). The Voigt bound is a good estimate of

_{z}*G*, while the Reuss bound is close to

_{xy}*G*and

_{yz}*G*(Fig. 6b). The shear modulus is larger when shear stress is in the same orientation as the shale bedding planes. Poisson's ratio

_{xz}*Y*when the loading is in

*X*, is higher than the arithmetic average; while

The arithmetic and harmonic averages also work as rigorous upper and lower bounds on Young's and shear modulus in Model 2. However, they cannot provide an appropriate estimate of the upscaled elastic properties as the shale configurations in each upscaled block become more complex. Because of the 30° shale layer inclination along *X*, upscaled *E _{x}* becomes smaller than

*E*(Fig. 6d). The upscaled Poisson's ratio for Model 2 is not constrained by the arithmetic and harmonic averages, and Figure 6f illustrates the development of three clear clusters of anisotropy.

_{y}### Non-symmetry of the upscaled compliance matrix

For each fine grid, symmetrical strain tensors are assumed if under symmetrical stress tensors; therefore, the compliance matrix should be symmetrical according to Hooke's law. In upscaling problems, the strain tensors are not guaranteed symmetrical due to heterogeneities at the subgrid scale, even under symmetrical stress tensors. As a consequence, symmetry may not hold for the upscaled compliance matrix when there are lithological trends that deviate from perfect layering. Rijpsma & Zijl (1998) provided derivations showing that the thermodynamic foundations of elasticity theory are not violated by a non-symmetrical upscaled compliance matrix.

The symmetry ratio is defined to further investigate the relationship between symmetry of compliance matrix and the geometry of shale configurations:*XY* plane. For *XZ* and *YZ* planes results in the ratios varying from unity. With a substantial departure from planar symmetry in Model 2 as the shale layers start to incline, the compliance matrix becomes more non-symmetrical, as shown by the scatter in symmetry ratios (Fig. 7b). The skewness of the auxiliary compliance tensor is strongly related to the non-symmetry of the depositional body. The previous upscaling method, including Backus' averaging method, cannot correctly capture the non-symmetry of the upscaled compliance matrix which may lead to inaccurate predictions of the upscaled elastic properties for highly heterogeneous systems similar to Model 2.

### Validation of global model response

The successful application of a local upscaling technique should also demonstrate its ability to match or reproduce the global geomechanical response of the fine-scale model. In this study, this implies that under the same loading conditions, the 64-cell upscaled model should, within an acceptable error, provide the same geomechanical response as the 8000-cell fine-scale model. The following sections compare the axial stress, shear stress and volumetric strain for models 1 and 2 for the fine-scale and corresponding upscaled geometries.

#### Global compression response

For Model 1, three uniaxial compression tests in the *X*, *Y* and *Z* directions are conducted on the fine-scale (8000 cells) and the upscaled (64 cells) model by applying a constant displacement rate of 10^{−4} m/step (Fig. 4a–c). Figure 8 shows the results from the numerical *Z*-direction compression test on Model 1. The *Z* direction is chosen to demonstrate the difference in stress–strain behaviour between the different upscaling methods because it is perpendicular to the plane of symmetry in Model 1. E-upscaled model predictions are able to match the axial stress and volumetric strain results from the fine-scale model (Fig. 8a and b). The harmonic average coarse model also performs well for axial stress prediction because the numerical experimental set-up is analogous to the Reuss iso-stress model. Predictions of volumetric strain (Fig. 8b) are less consistent because the Reuss bound is not strictly valid for Poisson's ratio.

As shown in Figure 9a, the Voigt (arithmetic) upscaling approach nearly matched the fine-scale response for compression loading in Model 1 for the *X* and *Y* directions. These results by the Voigt bound approach are expected since the configuration of the numerical experimental set-up is identical to Voigt's iso-strain model; however, the mismatch error exceeds 10% for loading in the *Z* direction. For all directions of loading, the E-upscale model matches the axial stress response of the fine-scale model with <1% error. It is expected that the ability of the E-upscale model to match the global model behaviour depends on the finite-difference discretization scheme. A sensitivity analysis of the upscaling ratio should be conducted to validate the global model responses of small-scale models prior to applying E-upscale to the full-field simulation models. For other loading configurations, the performance of the other homogenization techniques is poor and highly dependent on the direction of loading (Fig. 9a and b), providing additional evidence of the applicability of E-upscale to provide anisotropic elastic properties.

For Model 2, uniaxial compression tests in the *X* and *Z* directions have been chosen to assess upscaling methods when inclined shale layers produce non-uniform deformations. For this case, arithmetic and harmonic averaging are unable to replicate the axial stress behaviour (Fig. 10a and c) because the geometry of the model violates the assumed boundary conditions for the Voigt or Reuss oversimplified model. For these geological settings, E-upscale provides the best estimation of upscaled elastic properties and predictions of axial stress and volumetric strain. Fine-scale volumetric strain can only be approximated by the E-upscaled model (Fig. 10d) as the Voigt and Reuss bounds do not apply for Poisson's ratio. Again, for these compression loading assumptions, E-upscale is able to provide anisotropic elastic properties in highly heterogeneous models with inclined lithologies.

#### Global shear response

The E-upscale methodology is validated for the global model response in axis-parallel shear tests with the stress boundaries shown in Figure 4d–f but with 10 MPa shear stress to show larger deformations for comparison. E-upscale is able to replicate the results of the fine-scale model with an average error of <0.2% (Fig. 9c). Significantly larger errors, in some cases as high as 10%, occur when using the three other upscaling techniques. Subgrid-scale heterogeneity results in an anisotropic deformation response of the model under shear-stress loading boundary conditions that cannot be captured using arithmetic, geometric or harmonic averaging. The inclined shale layers present in Model 2 result in even greater non-uniform deformations (Fig. 11). For an upscaling ratio of 125:1, the E-upscaled model can capture the shape of the deformed model for three different shear-stress boundaries. The mismatch error in shear strain for E-upscale is *c*. 1%, whereas for the arithmetic (Voigt) and harmonic (Reuss) averages, mismatch errors are 5–10 times larger due to oversimplified uniform strain/stress assumptions.

## E-upscale applied to a coupled reservoir–geomechanics SAGD simulation

As shown in the earlier subsection on ‘Global compression response’, upscaled models using E-upscale can capture the geomechanical behaviour of a fine-scale model under both compression and shear-stress loading conditions. To assess the computational efficiency gained in adopting E-upscale, it is applied to a field-scale coupled reservoir–geomechanics SAGD simulation. The interactions of geomechanical and flow response are widely accepted to have a significant impact on reservoir performance and cap-rock integrity in SAGD projects. The increase in pore pressure and thermal expansion associated with steam-chamber development in heterogeneous oil sands reservoirs can lead to complex deformations. Laboratory experiments show that volumetric strain changes due to these deformations alter porosity, absolute permeability and relative permeability (Oldakowski 1994; Touhidi-Baghini 1998). Field observations of SAGD projects also provide evidence that high-porosity and high-permeability zones ahead of steam chambers allow hot fluids to flow into colder regions of the reservoir (Ito *et al.* 2000; Birrell 2001; Collins 2007). Numerical simulations for these scenarios have typically adopted uncoupled or coupled reservoir geomechanical modelling approaches. Geomechanical modelling requirements impose different constraints to reservoir modelling and typically lead to larger upscaling ratios to improve computational efficiency. This process is only effective if the upscaled model can honour the behaviour of the fine-scale geocellular model.

### Fine-scale reservoir–geomechanics SAGD simulation

A fine-scale coupled reservoir–geomechanics simulation was generated based on the geological model (Fig. 3a) of the MacKay River oil sands area. The Wabiskaw and McMurray formations are the zones of interest for the reservoir and geomechanical models. The Till, Clearwater and Devonian formations are only included in the geomechanical model as the over- and under-burden for the reservoir (Fig. 12a). The geomechanical and reservoir properties for the fine-scale reservoir–geomechanics simulation are summarized in Table 3. A steam-injection pressure of 1283 kPa with a corresponding steam-injection temperature of 191°C was chosen for the simulation (Suncor MacKay River Project 2017). Figure 12b shows the location of *XZ* cross-sections within the model that are used to illustrate model heterogeneity and modelling results. The initial horizontal permeability distribution (Fig. 13) for the fine-scale model in *XZ* cross-sections 1, 20, 40 and 60 highlights the degree of heterogeneity within the model.

For the simulation, a preheating period was modelled where steam is circulated within the injection and production wells until Day 150. Steam is then injected in the upper well and fluids are produced from the lower well, under steam trap control (Suncor MacKay River Project 2017). For the preheating period, absolute permeability is updated every 30 days until Day 150. For the production period, sequential coupling between the reservoir and geomechanical processes occurs every 100 days until Day 4760, which is the end of the simulation. The absolute permeability is updated based on a correlation generated for the McMurray Formation oil sands at low confining stresses (Touhidi-Baghini 1998), which is a simplified form of Tortike and Farouq's correlation (Tortike & Farouq 1991) incorporating volumetric strain into the changes of absolute permeability:*i* represents the direction of permeability (*X*, *Y* and *Z*) and *n* is the number of coupling steps. Constants *a* and *b* are tuned to fit laboratory and field data. Constants *c* and *d* are the upper and lower bounds of the updated permeability to avoid spurious calculated permeability. Table 4 provides the sand and shale tuning parameters, and the upper and lower boundaries for the permeability updating function.

The simulations predict the evolution of deformations within the reservoir and surrounding formations as the fluid pressures and temperatures change during steam-chamber development. For *XZ* cross-section 20, Figure 14 illustrates the growth of the steam chamber and the temperature field at days 1060, 1760, 3060 and 4760. For the same time periods, reservoir–geomechanical coupled processes occurring during the development of the steam chambers impart significant changes to the mobility of bitumen (e.g. horizontal permeability: Fig. 15), and these changes exert a significant influence on the geometry of steam-chamber development. These results demonstrate that flow and geomechanical response are closely integrated into SAGD, and it is critical that upscaling workflows adopted for these types of coupled simulations honour the geomechanical response of fine-scale models.

### Flow-based permeability upscaling

To facilitate the application of E-upscale to field-scale modelling problems, the workflow is implemented in the sequentially coupled reservoir-geomechanics simulation platform (Reservoir Geomechanics Research Group, University of Alberta, Canada). Weighted volume averages are used for porosity and saturation upscaling as they are not directionally dependent. The upscaled permeability tensor is calculated by the pressure-field solution method through the conservation of mass from the fine scale to the coarse grid. The upscaling of permeability has been extensively investigated as it is one of the key factors in flow response, according to Darcy's law:*u*, *k*, *μ* and

The effective permeability of the coarse block can be obtained by the integration of flow over the volume of the upscaled region:*V* is the volume of the upscaled cell and *k*_{eff} is the effective permeability tensor. The total flux (*Q*) in the upscaled region (left-hand side of equation 26) is calculated by the summation of flux at all permeable boundaries:*n* is the outward unit normal to the pressure boundary and *S* is the area of the face boundary.

A steady-state flow boundary is assumed in the local flow tests in the *X*, *Y* and *Z* direction. Thus, the effective permeability *k*_{eff} in the direction *i* is calculated by:

### E-upscale in a coupled reservoir–geomechanics simulation

Due to depositional processes, the vertical range of spatial correlation is usually shorter than horizontal (e.g. thin shale layers in oil sands reservoir). The longer horizontal range of continuity in shale bedding planes allows for larger upscaling ratios horizontally without over-mixing of the sand and shale facies. For the upscaled 6 × 6 × 1 model, the horizontal upscaling ratio is fixed to 6 while keeping the vertical grid size. E-upscale, arithmetic, geometric and harmonic averages are used for the upscaling of elastic properties. Flow-based upscaling is performed for all models, leading to the same initial permeability fields.

The cumulative oil production (COP) for fine-scale and upscaled cases are illustrated in Figure 16. The E-upscaled model has the most accurate oil recovery prediction compared to the fine-scale response. The arithmetic (Voigt) and harmonic (Reuss) bounds of the effective elastic properties lead to deviations in predicted oil recovery from the fine-scale model due to the complex stress state that occurs during the SAGD process. As the intermediate average between the Voigt and Reuss bounds, the geometric average has a slightly better prediction in terms of oil recovery as it provides a more accurate deformation response for the upscaled model.

Surface heave is commonly monitored in SAGD operations to avoid potential cap-rock integrity issues. While all upscaled models do not replicate the predicted surface heave from the fine-scale model, the E-upscaled model is the most accurate prediction of magnitude and shape of surface heave (Fig. 17). The displacement field shows that vertical displacements dominate the reservoir deformations (Fig. 18) in this model with confined side and bottom boundary conditions. The side view of the vertical displacement contour for *XZ* Section 60 also demonstrates the improved prediction of E-upscale. The deformation and permeability change of the reservoir is closely related to the steam-chamber propagation of the reservoir. Therefore, E-upscale does improve the overall prediction of the upscaled coupled simulation.

The reduction in computational time is significant for fine-scale coupled simulation. Figure 19 illustrates the strong dependence of simulation time on the number of model cells. The simulation time using the E-upscale workflow is 0.75 h compared to 15.4 h for the fine-scale model.

## Summary and conclusions

A local numerical upscaling technique, E-upscale, is developed to estimate the effective elastic geomechanical properties for heterogeneous continua. For each upscaled cell, the anisotropic Young's modulus, Poisson's ratio and shear modulus are solved to replicate the local strain tensors of the fine-scale model under three compression and three shear-stress boundaries in a parallel program. The global compression and shear response are validated for two different sand–shale sequence systems: Model 1 with horizontal shale bedding planes (transversely isotropic); and Model 2 with 30° inclined shale layers (fully orthotropic). Successful application of E-upscale in reservoir-scale coupled SAGD simulation demonstrates that E-upscale can significantly improve the prediction of surface heave, subsurface displacements and oil production when integrated with flow-based permeability upscaling.

The anisotropic E-upscaled elastic properties are solved with different numerical testing regimes to reproduce the anisotropic geomechanical response caused by subgrid-scale heterogeneities. The orthotropic assumption of the compliance matrix may not be applicable for reservoirs with massive discontinuities (faulting and fractures), where further investigation is required to confirm whether the anisotropic compliance matrix with 21 or 36 constants can deal with upscaling problems in heterogeneous discontinua. E-upscale can help to better accommodate the complex interactions of flow and deformations in reservoir processes. When significant inelastic behaviour is involved, additional simulation time steps are required for the local tests to determine the effective geomechanical properties which can capture the effective/dynamic modulus of elastoplastic behaviour and failure criteria.

As flow parameters are strongly related to deformation response, it is recommended that E-upscale be applied in recovery processes involved with significant pressure and temperature change (e.g. SAGD). The improved predictions of the E-upscaled model can enable a larger upscaling ratio for full-field geological models while maintaining the accuracy of the coupled simulation. The reliable surface and subsurface deformations associated with production performance provided by the E-upscaled model for full-field pilots can provide valuable predictions to make an improved resource recovery plan.

## Acknowledgements

The authors would like to acknowledge Abel Juncal for his valuable suggestions regarding the E-upscale approach. The valuable feedback received from anonymous reviewers of the original manuscript is gratefully acknowledged.

## Funding

This research was funded by the Energi Simulation Industrial Research Consortia in Reservoir Geomechanics.

## Author contributions

BZ: Conceptualization (Lead), Data Curation (Lead), Formal Analysis (Equal), Investigation (Lead), Methodology (Lead), Software (Lead), Validation (Lead), Visualization (Lead), Writing – Original Draft (Lead), Writing – Review & Editing (Lead); ND: Methodology (Supporting), Software (Supporting), Supervision (Supporting), Validation (Equal), Writing – Review & Editing (Supporting); MK: Formal Analysis (Supporting); RC: Conceptualization (Equal), Funding Acquisition (Lead), Supervision (Lead), Writing – Review & Editing (Equal); JB: Conceptualization (Supporting), Investigation (Supporting), Supervision (Lead), Writing – Review & Editing (Equal).

- © 2019 The Author(s). Published by The Geological Society of London for GSL and EAGE. All rights reserved