Multiscale simulation of inelastic creep deformation for geological rocks

Subsurface geological formations provide giant capacities for large-scale (TWh) storage of renewable energy, once this energy (e.g. from solar and wind power plants) is converted to green gases, e.g. hydrogen. The critical aspects of developing this technology to full- scale will involve estimation of storage capacity, safety, and eﬃciency of a subsurface formation. Geological formations are often highly heterogeneous and, when utilized for cyclic energy storage, entail complex nonlinear rock deformation physics. In this work, we present a novel computational framework to study rock deformation under cyclic loading, in presence of nonlinear time-dependent creep physics. Both classical and relaxation creep methodologies are employed to analyze the variation of the total strain in the specimen over time. Implicit time-integration scheme is employed to preserve numerical stability, due to the nonlinear process. Once the computational framework is consistently deﬁned using ﬁnite element method on the ﬁne scale, a multiscale strategy is developed to represent the nonlinear deformation not only at ﬁne but also coarser scales. This is achieved by developing locally computed ﬁnite element basis functions at coarse scale. The developed multiscale method also allows for iterative error reduction to any desired level, after being paired with a ﬁne-scale smoother. Numerical test cases are studied to investigate various aspects of the developed computational workﬂow, from benchmarking with experiments to analysing the impact of nonlinear deformation for a ﬁeld-scale relevant environment. Results indicate the applicability of the developed multiscale method in order to employ nonlinear physics in their laboratory-based scale of relevance (i.e., ﬁne scale), yet perform ﬁeld-relevant simulations. The developed simulator is made publicly available at https://gitlab .tudelft .nl /ADMIRE _Public /mechanics. Numerical results are presented to study various aspects of the developed computational workﬂow, from benchmarking with experiments to analysing the impact of nonlinear deformation for a ﬁeld-scale relevant environment. Results indicate the applicability of the developed multiscale method in order to employ nonlinear physics in their laboratory-based scale of relevance (i.e. ﬁne scale), yet perform ﬁeld-relevant simulation.


Introduction
A successful transition towards cleaner energy relies on the availability of large scale renewable energy storage. To achieve this, renewable energy can be converted to green fuels such as hydrogen [1] or green methane to be stored in large-scale quantities in geological subsurface [2]. Compressed air energy storage [3], energy stored in the form of heat [4,5] are other alternatives to store energy in the subsurface, which are cyclic in nature too. Depending on the supply and demand, the energy stored can vary in the range of different time scales and different capacities of storage volumes [6].
Clayed rock Yang et al. [26] ε cr = s−D b s cr Subsurface storage technology has been implemented worldwide in many countries including the Netherlands [7]. To safely utilize the underground storage capacities (e.g. in porous rocks and salt caverns [8][9][10]), it is critical to understand the underlying physics involved when a cyclic load is imposed. Cyclic loading has been studied in mining sciences to study hydraulic fracturing, rock cutting, forecasting volcanic hazard, designing tunnels against earthquakes (see e.g. [11][12][13][14]). Very few researchers have studied the applications of cyclic loading on oil and gas storage both numerically and experimentally [15][16][17]. Relevant hysteric rock and fluid constitutive governing equations can be implemented to capture the important physics in the cyclic subsurface storage technology [18][19][20]. However, there is a significant need for a better understanding of the complexities involved in the subsurface storage physics to implement it safely, efficiently, and reliably on a long term basis. Especially this becomes a crucially important field of study when storage cities become closer to the urban areas.
Of particular interest is to develop a scalable (i.e., multiscale) mechanical deformation modelling concept to include not only the elastic but also time-dependent nonlinear creep physics. As geological domains span large length scales, having inelastic (nonlinear) and heterogeneous properties, development of a multiscale strategy which allows to capture these finescale physics and parameters in the coarse-scale systems is crucial. Such a multiscale strategy would also allow for accurate quantifications of the consequences of the pore-pressure dynamics in the subsurface geological domains. This motivates the development of the current study.
Complex thermal, geo-mechanical, and chemical effects might occur simultaneously making it difficult to predict the feasibility of the energy storage [27]. This work is focused on studying the influence of in-elastic geo-mechanical effects. These effects include creep, plasticity and viscoplasticity. Few researchers have studied plasticity and visco-plasticity extensively in porous media [28][29][30]. In this study we focus on nonlinear time-dependent creep deformation.
Creep deformation depends on the type of material and the type of loading involved. Examples can be found in deep resource mining [31] and radioactive nuclear waste disposal [25], in geological reservoirs such as shale rocks [32,33] to name a few. Several constitutive laws have been introduced in the past to study creep and this still is an active field of research in the rock physics labs and geophysical field studies. Table 1 show different creep formulations employed in the past for different materials. In this work, the constitutive equation for modelling inelastic creep deformation is chosen based on the existing literature which is elaborated in section 3.
The present work introduces a multiscale 3D non-linear finite element method (FEM) to study the creep mechanics that is involved in the compression of rocks under normal and cyclic loading. Classical creep and relaxation creep are the two methodologies that are used to analyze and quantify the creep. Both are considered in this work. The computational framework allows using both explicit and implicit time integration scheme. However, implicit time integration is employed to guarantee convergence irrespective of the size of the time steps. For field-scale applications, one has to employ the physical constitutive laws which are validated in the laboratory settings (cm-scale) only.
Applications of these experimental laws for large (several tens of meters) grid blocks, having heterogeneous elastic and nonelastic properties, are questionable. To avoid excessive upscaling of these nonlinear terms and highly heterogeneous elasticity coefficients, a multiscale procedure is also introduced. The multiscale method is built on finite element (hence, MSFE) to obtain coarse-scale systems using locally-supported finite-element coarse-scale basis functions. The multiscale formulation of this type, i.e. local basis function-based, has been applied and developed by various researchers in the past. ( [34][35][36][37][38][39][40][41][42][43]). Researchers have employed multiscale strategy to study flow in porous media [44][45][46], geo-mechanics in porous media [47][48][49][50][51] and other non-linear problems [52][53][54]. Note that the rich domain of multiscale modelling for inelastic deformations also includes non-basis-function-based approaches, such as FE 2 for composite materials [55], FEM-DEM (Discrete Element Method) for granular materials [56] and high-porosity rocks [57]. The multiscale strategy developed in this work aims to (1) develop local coarse-scale basis functions, (2) resolve fine-scale heterogeneities in coarse system, (3) allow for iterative error reduction to the machine accuracy, (4) allow for algebraic formulation and implementation. In these aspects, it stays novel compared to the existing literature of multiscale FEM for inelastic deformation modelling. Note that in order to guarantee convergence of the iterative multiscale error reduction strategy, the developed multiscale procedure is paired with a fine-scale smoother [58,59].
Numerical results are presented to study various aspects of the developed computational workflow, from benchmarking with experiments to analysing the impact of nonlinear deformation for a field-scale relevant environment. Results indicate the applicability of the developed multiscale method in order to employ nonlinear physics in their laboratory-based scale of relevance (i.e. fine scale), yet perform field-relevant simulation.
The paper is structured as follows. Firstly, FEM method for simulation of linear elasticity is briefly revisited in section 2. Different creep mechanisms, as inelastic deformation, are then introduced in section 3. Constitutive equations employed to model creep are presented. The algorithms employed to model creep are described in detail. Next, the algebraic MSFE procedure for nonlinear deformation mechanism is developed to allow for field-scale analyses without relying on excessive upscaling of the nonlinear creep physics and heterogeneous material properties in section 4. MSFE formulation is further extended by employing iterative MSFE. Numerical results are then presented in section 5 starting with (i) Uniaxial compression test case for rock salt for both fine-scale and algebraic MSFE. Next, (ii) triaxial loading of coal is studied using both fine-scale and iterative multiscale strategy. Next, (iii) tri-axial cyclic loading test case is studied for sandstone using finescale and iterative multiscale strategy. For the above test cases, the numerical results are compared with experimental data by using fitting parameters to understand the variation of deformation over time. Next, the influence of creep on subsidence is studied in the northern Italian field [60] by assuming plain strain. Finally, in section 6, concluding remarks are presented.

Linear elasticity
The momentum balance for drained solid porous medium is given by where pore pressure is constant. The stress is given by σ = C : ε, in which the elasticity coefficient matrix and strain tensor are denoted by C and ε, respectively. Here f is the volumetric force expressed as [N/m 3 ]. The governing equation when expressed in terms of deformation reads Here, ∇ s u = 1 2 ∇u + ∇ T u is the gradient symmetric of displacement vector. The elasticity tensor is expressed using Lamé parameters λ and μ, i.e., , where E is Youngs modulus and ν is Poisson's ratio. Using Voigt's notation and transpose of divergence operator, the final set of equations in 3D are

FEM formulation of the system
To solve equations (4a)-(4c), finite element formulation is followed. Using minimization of weighted residual Galerkin approach, the displacement of each element is interpolated within an element using fine-scale shape functions. The displacement inside an element (ũ e ) can be approximated from nodal displacements (û e,i ) as Here, û e,i can be written as displacements (u x , u y , u z ) for each node i. N i,e are the shape functions used to interpolate solution between the nodes [59]. The shape functions are expressed in terms of local coordinates (ξ 1 , ξ 2 , ξ 3 ) by  [61,22].
Using the Gauss divergence theorem, the final discrete equation for an element reads The above volume integrals are approximated using Gauss quadrature rule by transforming it from global to local coordinate system. Here, D e = di v T N e , the applied stress on the boundary along the normal is t (Neumann bc) which is integrated over area σ , the volumetric force is f , integrated over an elemental volume e . The operator di v T is given by By assembling the stiffness matrix for each element, the global linear system is formed as A mechũ =f . (10) After computing elastic deformation, stresses and strains are computed using constitutive law from Equation (1). Next, the nonlinear elastic deformation (creep) is revisited.

Creep formulation
Creep is the tendency of solid material to deform permanently under certain load that depend on time and temperature. Typical creep process has three phases which are primary creep (creep rate decreasing over time), secondary creep (constant creep rate) and tertiary creep (increasing creep rate until failure) as shown in Fig. 1. In this work, tertiary creep is not modelled. When the strains are infinitesimal, the total strain ε t can be split as ε t = ε el + ε cr + ε th . (11) Here, ε el is elastic strain, ε th is thermal strain and ε cr is the creep strain. In this work, time independent plastic strains like viscoplasticity and thermal strains are not included. The modified creep equation which involves power law and Bingham body are [26,24] ε cr = X t P + Y B.  Depending on the type of material and formulation either power law ( X = 1, Y = 0) or Burgers model ( X = 0, Y = 1) can be chosen. Here, Material constants are A, n, m in the power model and in the Burgers model the parameters G 1 , G 2 are the shear modulus, η 1 , η 2 is the viscosity coefficient and t is the time. Here s is deviatoric stress and σ vM is von Mises stress. The Burgers model for 3D stress state is obtained as series of elements which represent Hookean body (elastic, spring), Kelvin body (damper and spring in parallel) and Newton body (dashpot/damper). So the stresses of the individual elements are given by The non-linear strains are incorporated in the governing equation by plugging Equation (11) in Equation (2), The deformation is computed by Using the principal of virtual displacements, which states that sum of external and internal work is 0 [21], the final discrete equation can be written as modified form of Equation (7), i.e., Equation (17) can be rewritten as Where, F cr is the fictitious creep forces and F b is the body forces.

Numerical methodology
The creep strain ε cr can be modelled using classical creep and relaxation creep methodology [62]. Classical creep is a modelling approach, when the stress is assumed to be constant and the creep strain accumulates with respect to time. If the stress is constant, the elastic strain remains constant which would mean total strain is accumulated over time. Whereas, the relaxation creep is an approach when the total strain is assumed to be constant, and as creep strain accumulates, elastic strain ε el reduces over time. This would imply that the stress reduces with time. Fig. 2 shows the schematic diagrams of classical and relaxation creep. These numerical methodologies are similar in nature to the stress and strain controlled compression experiments on the rocks, where either stress or strain is controlled on the rock sample when the experiments are performed.
The algorithms for classical and relaxation creep are depicted in Algorithm 1 and Algorithm 2 respectively. The developed computational framework allows to use both implicit and explicit formulation. The results shown here are obtained from implicit formulation which ensures computational stability even for higher time step sizes. The non-linear creep formulation is discretised implicitly as, Calculate creep strain at time step i : ε i Compute total strain at time step i : ε i t = ε e + ε i cr 8: Calculate creep strain at time step i : ε i Calculate elastic strain ε i el = ε t − ε i cr 8: Calculate stress σ i = C : ε i el 9: In the classical creep formulation, the stress remains constant with respect to time. However, in relaxation the stress is dependent on time. In general, the stress term in the creep model is obtained based on the new time step (i + 1). In this case, Newton-Raphson iterative approach is employed to solve for the displacement iteratively. The nonlinear residual reads as To find u i+1 the above equation is further linearized, i.e., The iterations ν are repeated until the residual norm falls below the prescribed threshold, i.e., || ν || 2 < r to achieve convergence. Algorithm 3 describes the implicit formulation to model creep. Set iteration counter ν = 0, creep forces F i+1 : while || ν || 2 < r and ν < max(ν) do

Multiscale finite element method for nonlinear deformation under cyclic loading
The approximate fine scale displacement ū h is computed from coarse scale system ū H using Prolongation operator P h H , i.e., The local conditions are solved in each coarse element u i for basis functions N H i of i coarse block. Here δ ij is the Kronecker delta and x j is the position vector of coarse mesh node associated to j-th unknown. ∇ s and ∇ s || denote approximate divergence and symmetric gradient operators acting in tangential direction of coarse elements boundary (local coordinate system) to solve local boundary conditions at ∂ u for faces and edges, respectively.
Based on the wire basket decomposition of fine scale system [63], a permutation matrix W is employed to rearrange all the nodes in the order of interior, faces, edges and vertex. Based on that the fine scale system can be written aŝ A hû =f h . (25) Denoting I, F, E and V as interior, face, edge and vertex node, relative to the coarse mesh configuration, allows for re-stating the permuted system aŝ The above matrix is reduced by performing Gaussian elimination twice which gives ⎡ ⎢ ⎢ ⎣Â Here,Ŝ 7 By applying the reduced order boundary conditions it can be stated that The matrices are computed by modifying the divergence and gradient operators for faces and edges respectively and solving Equation (24). The information from the faces is interpolated to faces, edges and vertices by reduced boundary conditions, and information from edges is interpolated to edges and vertices. This meanŝ By implementing Equation (24), the modified system is ⎡ ⎢ ⎢ ⎣Â The restriction operator is the transpose of prolongation matrix, i.e., Using the coarse scale deformation and employing Equation (23), fine scale displacements are computed.  Fig. 4 shows the basis functions at the vertex, edge, face and interior coarse element. The coarsening ratio along each axis is given by The MSFE procedure for inelastic creep deformation is shown in Algorithm 4. The normalized average error between the fine-scale and multiscale solution is given by The solutions obtained from MS strategy employ fewer degrees of freedom which leads to the reduction of computational effort compared to finescale strategy.

Iterative multiscale solver
When creep strains are included, the solutions obtained from multiscale formulation are approximate. To reduce the error to a desired tolerance, an iterative procedure is used, employing a fine scale smoother like ILU or GMRES [64,65]. In this work, ILU(0) and GMRES is employed for few test cases to compare the higher computational efficiency and also improve MS solutions to the desired tolerance. The stopping criteria used here e r is given by Where e r is the average error between fine-scale and multiscale solutions. In realistic field cases, fine scale solution will not be available. In those cases, the stopping criterion can be the relative residual, i.e., The procedure for iterative formulation is described in Algorithm 5. The smoothing stage (stage -2) can be applied n s number of times, depending on the complexity of the problem and the computational efficiency.

Convergence test
To check the consistency of the fine-scale 3D FEM formulation, a synthetic test case for mechanical equilibrium is developed which is shown in Equation (40)   Table 2.

Synthetic shear test case
In this subsection, a synthetic test case is modelled with creep to show the capabilities of the computational framework and understand the influence of creep on the body deformation. Rock salt material properties are chosen for modelling this test case. Fine-scale and multiscale strategies are employed to this model. The parameters used in this test case are shown in Table 2, see also Fig. 6 for the schematic diagram. The parameters are chosen artificially, to demonstrate the deformation caused due to creep.
A horizontal shearing force is applied on the top face, the bottom face is fixed and rest of the faces are traction free. Fig. 7b shows the snapshot of the deformed body at time = 0, obtained from both fine-scale and multiscale strategies. Fig. 7b shows the finescale and multiscale solutions obtained after time = 1.5 years. This deformation occurs due to creep.
The error between the FS and MS results as calculated from Equation (37) is less than 10%. Iterative MS is not further employed in this test case. The solution obtained from MS strategy is not accurate because the reduced boundary conditions have not provided an accurate localisation condition. fixed. Fine-scale grid has 25 × 25 × 25 cells, and the coarsening ratio is 5 × 5 × 5. Fig. 7a shows the results at time = 0 for fine-scale (Surface) and MS (Wire-frame). Fig. 7b shows the results at time = dt for fine-scale (Surface) and MS (Wire-frame) and white outline is the fine-scale result at time = 0.
The parameters used in this test case are presented in Table 2.

Uni-axial compression of rock salt
Uni-axial compression of rock salt is investigated in this subsection. A distributed load is imposed on the top face and all the four faces (north, east, west and south) of the domain are subjected to roller constraints. The bottom face is fixed. The schematic is shown in Fig. 9a. A distributed load P of 0.1-0.3 MPa is imposed for a period of 8 months. The results for P = 0.2 MPa are compared with experimental results [68]. Least square fit is used on the strain rate to obtain the parameters of the creep formulation. Equation (13a) is used here as a constitutive equation. The parameters used to model creep are presented in Table 3. The experiments are conducted on cylindrical specimens, however the effect of curvature is ignored in this case and only the dimensions are taken.
Firstly, classical and relaxation creep methodology is employed to understand uni-axial compression. Fig. 8 shows the results of both classical and relaxation creep. Fig. 8a and Fig. 8b show the variation of absolute stress (σ zz ), absolute strains (ε t,zz , ε el,zz and ε cr,zz ) with time for classical creep. Here, we can see that the stress is constant with time which means that elastic strain is constant with time. The absolute of the total strain is increasing due to increasing creep strain. Fig. 8c and Fig. 8d show the variation of absolute stress (σ zz ), absolute strains (ε t,zz , ε el,zz and ε cr,zz ) with time for relaxation creep.
In this case, the total strain is constant with time, since the absolute of the creep strain is increasing the elastic strain is reducing with time. Since the elastic strain is reducing with time, the absolute stress is reducing with time. These plots depict the same methodology as shown in Fig. 2.
The results are shown in Fig. 9 are obtained from employing classical creep. All the results are plotted for a node that is located in the centre of the top face. The variation of strain rate with time is shown in Fig. 9b. It can be seen that the strain rate increases when the imposed load increases. Similar observation is seen in Fig. 9c. Least square fit is conducted on strain rate and accordingly the strain is computed. There is some difference between the experimental and numerical results for variation of strain with time Fig. 9c. This is because it is assumed that the properties of the rock salt are assumed to be homogeneous. However, Young's modulus and Poisson ratio vary depending on the salt composition. The presence of salt of varying compositions would affect the fitting parameters of the creep constitutive equation. Investigating the effects of the composition of rock salt on the overall behaviour of the numerical model is beyond the scope of this work. Fig. 9d shows the variation of deformation with time for different loading conditions shown for both fine-scale and multiscale. It can be seen that the results of fine-scale and multiscale overlap. As the imposed load increases, the elastic deformation and deformation accumulated over time due to creep increases. This is a situation in which the boundary condition imposed on local basis functions is very accurate, leading to accurate multi-scale solutions. Note that, the multiscale solution accuracy depends highly on the accuracy of the local basis functions. When the reduced boundary condition has provided an accurate localisation condition, the MS solutions are accurate.   Table 3.  (Fig. 10b), deformation (Fig. 10c) with time for coal specimen when it undergoes creep. The experimental results [26] are compared with numerical results using parameters shown in Table 4. The deformation results are compared for FS, MS and iterative MS. The convergence plot obtained by using iterative MS is shown in Fig. 10d for deformation (u z ). The convergence plot is shown for stage 2 of the iterative multiscale solver for different number of smoothing stages using only ILU(0) n s = 1, 2, 3 or using only GMRES without any preconditioner.

Tri-axial compression
In this subsection, influence of creep on tri-axial loading is studied. Fig. 10a shows the schematic of tri-axial loading where a confining stress is imposed on all the four faces (north, east, south and west), the bottom is fixed and a vertical load is imposed on the top face. Bingham body formulation is employed here, as in Equation (13b), to study the influence of creep on coal specimen. The parameters chosen to study this formulation are presented in Table 4. The experimental results [26] are compared with the numerical results. The experiment is conducted for a period of 45 minutes. The numerical results are compared for fine-scale, multiscale and iterative multiscale simulations. Fig. 10b shows the variation of creep strain and creep strain rate with time. Here it can be seen that the experimental data fits well with the numerical data for the parameters chosen. However, the experimental data also shows the tertiary creep phase which is not captured by the numerical formulation. To model tertiary creep, a different damage constitutive formulation has to be employed that can explain the coupling of highly non-linear effects which leads to material failure. Though it is of critical importance, it is beyond the scope of this work. The variation of deformation (u z ) with time is shown in Fig. 10c for a node that is located in the centre of the top face. It can be seen that the deformation increases with time for the imposed load. The primary and secondary stage of creep can also be seen in this formulation. The normalized error computed from Equation (37) is less than 2% for the chosen coarsening ratio.
The results of fine-scale and multiscale are slightly different. Multiscale solution is further improved by employing iterative multiscale strategy. ILU(0) and GMRES are compared for their efficiency of smoothing in i-MS strategy. Fig. 10d shows the variation of error for deformation (u z ) with number of iterations to achieve convergence. The tolerance was set to 1e-10.  11. Tri-axial cyclic loading: Fig. 11a shows the variation of total strain with time compared with the experimental results [69] for sandstone rock, and Fig. 11b shows the variation of deformation (u z ) with time for both fine scale and multiscale results. The parameters used for this formulation are presented in Table 5.
The convergence plot is compared for ILU (0) with varying smoothing stages n s = 1, 2, 3 and GMRES without any preconditioner. It can be seen that GMRES without any preconditioner converges to the solution in less than 5 iterations, at a faster rate compared to ILU (0). The deformation obtained from MS is not very different from FS. However, when the properties of rock are heterogeneous and different loading condition the convergence rate of the iterative solvers will vary.

Tri-axial cyclic loading of sandstone
In this subsection, a cyclic load is imposed when the rock is confined by a stress in all the 4 sides. The experimental results [69] are compared with numerical results for sandstone. The parameters chosen for this formulation are presented in Table 5. Investigation of the effect of varying material properties on the fitting parameters of creep is beyond the scope of this work. The simulation is conducted for a period of 210 minutes with 4 cycles for a vertical load of 265 MPa. Equation (13a) is used by employing classical creep in this study. Fig. 11 shows the variation of total longitudinal strain and longitudinal deformation for fine-scale and multiscale strategies. Fig. 11a compares the experimental and numerical total strain data. Least square fit is employed on the experimental data only on the first cycle. The fitting parameter m = 1.005 is very close to 1. This would mean that the strain rate is almost not dependent on time. In addition, it can be seen that the experimental and numerical data of strain match for more than 2 cycles. The observed error between the numerical and experimental data is slightly higher in the last cycle after 300 hours. We think that this might be because due to cyclic loading there might be opening and closing of micro-cracks or expansion of voids. To incorporate this, the fitting parameters of the creep constitutive formulation would depend on the type of loading like frequency and amplitude of loads. It is of critical importance to further investigate this effect. Due to lack of literature, this is beyond the scope of this work. Fig. 11b compares the deformation obtained from fine-scale and multiscale strategy. The multiscale solutions are obtained using much fewer degrees of freedom. As such it clearly outperforms the fine-scale simulation strategy. The variation of the deformation of the node that is located in the centre of the top face is shown here. It can be seen that, the approximate solutions obtained from multiscale solutions are very close to fine-scale solutions. Hence for this test case, iterative multiscale strategy is not further employed.

Influence of creep on plain strain subsidence
Subsurface energy storage technology will involve cyclic injection and production of the reservoir. Due to this, subsidence and uplift of the earths top layer might occur. To comprehend this, the influence of creep on subsidence and uplift has to be studied. In this work, the influence of non-linear deformation on land subsidence which is caused due to injection and production of reservoir is modelled in heterogeneous media by assuming plain strain in the domain (x − z plane). [70,60] gives you the properties of the geological formation inside the domain. The domain mainly consists of clayed rock. The  Fig. 12a shows the schematic of the plain strain subsidence test case where the pressurized H 2 gas is stored in the reservoir (small box) inside the geological domain. The domain is constrained by roller boundary condition at north, east, south and bottom. Fig. 12b depicts the distribution of Young's modulus in the x-z plane within the geological domain as described by Equation (44). The green box shows the location of the reservoir.
reservoir is modelled in the span of 10 km and for a height of 3 km. The elasticity modulus is obtained from vertical uni-axial compressibility (c M ) which is given by [71], i.e., Here c M and σ Y are expressed in [bar −1 ] and [bar] respectively. Here σ Y is the vertical effective stress. The vertical effective stress is obtained as superposition of total vertical stress σ Y and hydrostatic pressure p, i.e., The Young's modulus can be expressed as from [71], i.e., The schematic of this test case is shown in Fig. 12a. Here, rollers are subjected to the bottom, west and east boundary and the north boundary is traction free. The hydrogen gas has to be stored at right temperature and pressure for feasibility. In the natural state of the earth, lithostatic pressure would be present inside the geometrical domain. However, it is assumed that the lithostatic pressure would not cause nonlinear time dependent deformation. A deviation from the equilibrium state, which involves injection and production of the reservoir would undergo nonlinear time dependent deformation. Assume that the natural state of the reservoir has a pressure of P 0 = 150 bar without any subsidence. The effective pressure P is given by In this work, the pressure P is varied between 50 and 110 bar in the form of cyclic load as shown in Fig. 12b. If P = 0 bar, then P = P 0 which means that the earth is in equilibrium and there is no subsidence. For P = 50 bar ( P = −100 bar) higher subsidence is expected compared to P = 110 bar ( P = −40 bar) because of the higher pressure drop. The temperature is chosen to be 373 K [72], with minimum microbial and bio-chemical reactions [73]. Elasticity modulus degradation due to possible reservoir rock in-elasticity is neglected. The rock material is assumed to be isotropic. Classical creep methodology is employed here.
The fine-scale grid contains 110 × 110 cells, while the multiscale solver employs 22 × 22 cells, having the coarsening ratio of 5 × 5. The parameters for creep formulation are obtained from previous test case as shown in Table 5. The time-step (dt) size is chosen to be 4 months and the entire simulation is conducted for a period of 10 years (30 ×dt). The time-step size is assumed to be of similar order as Maxwell time which is given by ratio of shear viscosity to its shear modulus. The uncertainty in viscosity due to different grain sizes is not considered here. Poisson ratio is ν = 0.3 and Lame's parameters are obtained from Youngs modulus in Equation (44). Fig. 14 shows the results obtained for plain strain subsidence test case when nonlinear time-dependent deformation is incorporated. Fig. 13a shows the variation of subsidence inside the domain for different effective pressures P = −40, −100 bar when time = 0. The multiscale solutions are obtained using much fewer degrees of freedom. This reduces the computational effort and outperforms the fine-scale simulation strategy. It can be seen that the multiscale solution and fine-scale solution have similar magnitude except for the smooth kink which is not obtained in the multiscale solution. This is because the number of cells occupied by the reservoir is close to the number of fine-scale elements in a coarse element. If required, an iterative MS strategy can be employed to improve the solution. Overall, the quality of the multiscale solution is good. It can be seen that the lower is the effective pressure, the higher is the magnitude of the subsidence due to higher pressure drop compared to the natural state. The influence of creep is studied by injection and production inside the reservoir which varies the pressures as shown in Fig. 13b in the form of constant load and cyclic load for 10 years.   Table 5. The fine scale grid is 110 × 110. The multiscale solutions are shown for a coarsening ratio of 5 × 5. Constant pressure ( P = −100 bar) and cyclic pressure ( P = −100 to P = −40 bar) is imposed inside the reservoir as shown in (Fig. 13b) and the influence of creep on subsidence is shown in Fig. 14a and Fig. 14b respectively. The variation of the average error (Equation (37)) between fine-scale and multi-scale solutions for u x , u z at each time step is shown in Fig. 14c and Fig. 14d. Fig. 14a shows the subsidence inside the domain at t = 0 (without creep) and t = 10 years (with creep) at a constant load for effective pressure of P = −100 bar. It can be seen that the magnitude of subsidence has increased at t = 10 years compared to t = 0. Due to creep, the deformation increases which also increases the magnitude of the subsidence caused. The deviation between fine-scale and multiscale solutions has increased which is because the basis functions are constructed based on linear elasticity and the magnitude of the fictitious creep forces F cr increases with time. Fig. 14b shows the variation of subsidence inside the domain when the reservoir is injected and produced cyclically. The magnitude of subsidence caused only due to creep in cyclic load will be lesser than constant load (of higher magnitude). This is because a constant load of higher magnitude will result in higher magnitudes of stress which further results in the accumulation of higher magnitude of creep strain and deformation. It can be seen that the magnitude of error between FS and MS solutions for P = −40 bar at t = 0 (as seen in Fig. 13a) and at t = 10 years in cyclic loading (as seen in Fig. 14a) is lower. The error between the FS and MS solutions is shown in Fig. 14c and Fig. 14d for constant and cyclic load respectively. The normalized error is less than 13% for constant load and 11% for cyclic load. As time increases creep strain increases which leads to increase in the error. The error for u x increases at a faster rate than u z . The MS solutions can be further improved if the iterative MS strategy is used only near and around the reservoir as a block iterative MS solver. Overall the MS solutions are satisfactory, hence iterative MS strategy is not further employed in this case.

Conclusion
In this work, an iterative multiscale nonlinear modelling strategy was proposed which can model time-dependent nonlinear geomechanical deformation using finite element method (FEM). The developed computational framework would be of critical importance to model subsurface energy storage which involves nonlinear deformation. Constitutive formulae were chosen to model the creep strains of different rock materials from the available literature. Employing the least square fit on the experimental data, the fitting parameters of creep strain were computed. Classical creep and relaxation creep were both modelled incorporated. Due to the computational stability, an implicit time discretization strategy was followed. Multiscale finite element (MSFE) method, formulated and implemented algebraically, was then developed based on locally computed elastic basis functions. To further improve the approximate multiscale solution, 2 stage iterative strategy is developed to combine the multiscale system with a smoother at fine scale.
Several numerical experiments were conducted with synthetic and realistic material parameters, to understand the nonlinear time-dependent deformation. First, creep deformation of rock salt under uni-axial deformation is studied under both classical and relaxation creep. It was found that the fine-scale and multiscale solutions match well. Employing Bingham creep formulation, tri-axial compression of coal is studied. The difference between fine-scale multiscale solution was found to be less than 2%. Different iterative multiscale smoothing strategies were compared to obtain fine-scale solution of desired accuracy. Since subsurface storage technology in reservoirs involves cyclic loading, as such tri-axial cyclic loading on sandstone was also studied. In agreement with the previous cases, the difference between fine-scale and multiscale solutions was less than 1%. The influence of non-linear deformation was further studied on realistic reservoirs which involved land subsidence in heterogeneous fields. Without creep, the multiscale solutions show a similar magnitude as fine-scale solutions. For the studied realistic case, the quality of multiscale solution (with no iterations) was satisfactory even when creep is incorporated.
This study allows for quantification of the inelastic time-dependent creep deformation for geoscience applications, specially those relevant to subsurface storage. It also makes it possible to represent the complex nonlinear and heterogeneous systems with much fewer degrees of freedom than the fine-scale classical approaches. As such, it sheds new lights in advancing the accuracy and speedup of the quantitative analyses of complex natural porous media behaviour.
Future work includes poro-inelastic systems to study the influence of fluctuating pore pressure on the cyclic behaviour of the rock materials. It was also found in subsection 5.5 that as the number of cycles increases, the experimental data starts to deviate from the fitting data. This needs to be investigated further when the time-scale is in years. Propagation of faults and fractures will be incorporated inside the domain to study their influence on inelastic deformation when a cyclic load is imposed.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.