An Analytical Solution for Pressure-Induced Deformation of Anisotropic Multilayered Subsurface

Abstract

We present a generalized Geertsma solution that can consider any number of finite-thickness layers in the subsurface whose mechanical properties are different from layer to layer. In addition, each layer can be assumed either isotropic or anisotropic. The accuracy of the generalized solution is validated against a numerical reference solution. The generalized Geertsma solution is further extended by a linear superposition framework that enables a response simulation due to an arbitrarily-distributed non-uniform pressure anomaly. The linear superposition approach is tested and validated by solving a realistic synthetic model based on the In Salah CO₂ storage model and compared with a full 3D finite element solution. Finally, by means of a simple inversion exercise (based on the linear superposition approach), we learn that the stiffnesses of cap rock and reservoir are the most influencing parameter on the inversion result for a given layering geometry, suggesting that it is very important to estimate high-confidence mechanical properties of both cap rock and reservoir.

1. Introduction

Disturbing underground pore pressure by fluid injection or production may induce significant changes in stress and strain fields in a subsurface reservoir, which can then propagate to neighboring layers both down- and upwards and even to ground surface or seabed [1]. Such a mechanical behavior can be well described by the poro-elasticity theory (e.g., [2,3]) and have been used for characterization and monitoring of subsurface geomechanical behavior in many different areas of engineering and geoscience (e.g., [4,5,6]). The same mechanical behavior is relevant for understanding geomechanics related to geological CO₂ storage sites. In particular, the interferometric synthetic aperture radar (InSAR) data acquired at the In Salah CO₂ storage project in Algeria have been successfully applied to understand the pressure front evolution and hydraulic communication network in the subsurface during the whole period (around 7 years) of CO₂ injection (e.g., [7,8,9,10]). In addition, Hu et al. [11] showed how the change in temperature and pressure may influence CO storage capacity, which is also one of the main parameters together with the geomechanical integrity in the context of CO₂ storage characterization and monitoring.

When relating surface deformation to the change in pore pressure, stress, and strain in the subsurface, we need a well-developed solution-frame that can calculate the quantitative dependency between the surface deformation and disturbed pore pressure under realistic geological environment e.g., heterogeneous stratigraphy, arbitrarily-distributed pressure anomaly. The Geertsma’s solution [2] is well known to provide good insights (as a first-order approximation) on how ground deformation behaves, when a fluid is injected into or produced yet only from an isotropic homogeneous subsurface. In addition, the Geertsma’s solution is relevant only for a disk-shaped reservoir with a uniform pressure perturbation. Du and Olson [12] extended the Geertsma’s solution by numerically integrating the nucleus of poro-elastic strain over a rectangular prism so that the spatial change of pore pressure can be taken into account, followed by superposing results from each rectangular prism. Fokker and Orlic [13] presented an advanced semi-analytical solution that is based on a boundary element method and can consider a multi-layered visco-elastic subsurface. The semi-analytical solution can be faster than e.g., finite element approach, but is still rather slow compared to the analytical Geertsma’s solution. Tempone et al. [14] extended the Geertsma’s solution a bit further by adding a rigid layer beneath a compacting reservoir, and also reported a correction to an error found in [15]. Lately, Mehrabian and Abousleiman [16] not only provided a good overview of the analytical solutions developed further to address the limitations of the Geertsma’s solution, but also presented a poro-elasticity framework that can simulate an isotropic layered-stratigraphy.

The current study makes further improvements in such a Geertsma type problem by deriving a more-generalized analytical solution that can consider any number of finite-thickness layers, each of which can have not only different but also anisotropic stiffnesses. These novelties are not found in any of the previous studies mentioned above. The analytical solution requires a numerical Hankel transform, whose accuracy and efficiency are also discussed. The accuracy of the analytical solution is validated in comparison with a finite element solution. Then, we discuss the importance of the heterogeneity in the layers through a set of synthetic models. In addition, we discuss the possibility of applying the analytical solution for a spatially-varying pressure anomaly, not only a single cylinder-shaped constant pressure. The approach is based on a linear superposition and should be efficient, since it avoids the numerical integration presented by Du and Olson [12]. At the end, we apply the analytical solution to analyze a realistic synthetic surface heave data that is made based on the In Salah CO₂ storage site in Algeria.

2. Analytical Solution for Anisotropic Layered Subsurface

Here, we present an analytical solution for displacement field for an anisotropic layered subsurface subjected to fluid-induced pore pressure disturbance in a reservoir layer. The anisotropy herein refers to a particular case of vertically-transverse isotropic (VTI) stiffness, as shown in Figure 1. Note that the anisotropy (horizontal-to-vertical velocity ratio) can be applied individually to S- and P-wave velocities and each layer may have different velocity ratios. The analytical solution is derived for a constant pressure change within a reservoir layer, and the pressure change is of cylinder shape i.e., axis-symmetric problem. The reservoir thickness can be any value, not only very small as in [2]. In addition, any number of VTI or isotropic layers can be taken into account. Therefore, the analytical solution in the present study can be considered as a generalized Geertsma solution (GGS).

To derive the analytical solution for the VTI medium shown in Figure 1, we apply the following axis-symmetric governing equation in cylindrical coordinate (r, z) and Hankel transforms with k being the transform parameter (or wavenumber).

Governing equation:

$$\begin{gathered} \left[\left\{\begin{array}{cc}\lambda +2G& \cdot \\ \cdot & G_t\end{array}\right\}\frac{{\partial }^2}{\partial r^2}+\right.\left\{\begin{array}{cc}{G}_t& \cdot \\ \cdot & {\lambda }_t+2G_t\end{array}\right\}\frac{{\partial }^2}{\partial z^2}+\left\{\begin{array}{cc}\cdot & {\lambda }_t+G_t\\ {\lambda }_t+G_t& \cdot \end{array}\right\}\frac{{\partial }^2}{\partial r\partial z} \\ +\left\{\begin{array}{cc}\lambda +2G& \cdot \\ \cdot & G_t\end{array}\right\}\frac{1}{r}\frac{\partial }{\partial r}+\left\{\begin{array}{cc}\cdot & \cdot \\ {\lambda }_t+G_t& \cdot \end{array}\right\}\frac{1}{r}\frac{\partial }{\partial z}\left.-\left\{\begin{array}{cc}\lambda +2G& \cdot \\ \cdot & \cdot \end{array}\right\}\frac{1}{r^2}\right]\left\{\begin{array}{c}{u}_r\\ u_z\end{array}\right\}=-\nabla \left(\alpha p\right) \end{gathered}$$

(1)

Hankel transforms:

$$\begin{array}{cc}{u}_r={\int }_0^{\infty }{U}_1\left(z\right)k{J}_1\left(kr\right)dk \to \hfill & :\mathrm{radial}\mathrm{displacement}\hfill \\ u_z={\int }_0^{\infty }{U}_3\left(z\right)k{J}_0\left(kr\right)dk \to \hfill & :\mathrm{vertical}\mathrm{displacement}\hfill \\ p={\int }_0^{\infty }P\left(z\right)k{J}_0\left(kr\right)dk \to \hfill & :\mathrm{pore}\mathrm{pressure}\hfill \end{array}$$

(2)

In the equations above, u_r and u_z are the radial and vertical displacements, respectively; p is the pressure anomaly of radius R; ∇ is the gradient operator; (λ, G) and (λ_t, G_t) are the two pairs of the Lamé first parameter and shear modulus for the horizontal (or radial) and vertical directions, respectively; α is the Biot coefficient. It is also noted that the Lamé and shear moduli are related to the elastic P- and S-wave velocities as below by Park and Kaynia [17]:

$$V_p^2=\frac{\lambda +2G}{\rho },V_s^2=\frac{G}{\rho },V_{pt}^2=\frac{{\lambda }_t+2G_t}{\rho },V_{st}^2=\frac{G_t}{\rho }$$

(3)

where ρ, V_s, V_st, V_p, and V_pt are, respectively, mass density, radial/horizontal and vertical S-wave velocities, radial/horizontal and vertical P-wave velocities of each layer. Formally applying the Hankel transforms defined above to the governing equation, we obtain the governing equation in the k-z domain as below:

$$\begin{gathered} -\left(\lambda +2G\right)k^2{U}_1+G_t\frac{d^2{U}_1}{d{z}^2}-k\left({\lambda }_t+G_t\right)\frac{d{U}_3}{dz}=-k\alpha P \\ -G_t{k}^2{U}_3+\left({\lambda }_t+2G_t\right)\frac{d^2{U}_3}{d{z}^2}+\left({\lambda }_t+G_t\right)k\frac{d{U}_1}{dz}=\frac{\partial \alpha P}{\partial z} \end{gathered}$$

(4)

The solution to this linear system of ordinary differential equations can be obtained through a standard algebra and the final results are given below:

$$\begin{array}{cc}\hfill U_1& =\frac{1}{k\left(\rho V_p^2\right)}\alpha P+A{e}^{k_{1}z}+B{e}^{-k_{1}z}+C{e}^{k_{2}z}+D{e}^{-k_{2}z}\hfill \\ \hfill U_3& =A{\varphi }_1{e}^{k_{1}z}-B{\varphi }_1{e}^{-k_{1}z}+C{\varphi }_2{e}^{k_{2}z}-D{\varphi }_2{e}^{-k_{2}z}\hfill \end{array}$$

(5)

where:

$$\begin{array}{cc}\hfill {\varphi }_{\pm 1,\pm 2}& ={\displaystyle \frac{\left(V_{pt}^2-V_{st}^2\right)\left(\frac{k_{1,2}}{k}\right)}{V_{st}^2-V_{pt}^2{\left(\frac{k_{1,2}}{k}\right)}^2}}\hfill \\ \hfill \left({\displaystyle \frac{k_{1,2}}{k}}\right)& =\pm {\displaystyle \sqrt{\frac{1}{2}\frac{\left[V_p^2-\left(V_{pt}^2-2V_{st}^2\right)\right]}{V_{st}^2}\pm \sqrt{{\left(\frac{\left[V_p^2-\left(V_{pt}^2-2V_{st}^2\right)\right]}{V_{st}^2}\right)}^2-4\frac{V_p^2}{V_{pt}^2}}}}\hfill \end{array}$$

(6)

In addition, there are four unknown coefficients of A, B, C, and D for each layer, which can be determined by satisfying the interface conditions between adjacent layers (i.e., displacement continuity and traction equilibrium). Furthermore, by imposing the zero-displacement boundary conditions at the bottom of the last finite-thickness layer, we can also simulate the rigid layer beneath a compacting reservoir studied by Tempone et al. [14]. The related two traction components (S_rz and S_zz) at the layer interfaces can also be written in the k-z domain as below:

$$\begin{gathered} S_{rz}=\rho V_{st}^2\left[\left(+k_1-k{\varphi }_1\right)A{e}^{k_{1}z}-\left(k_1-k{\varphi }_1\right)B{e}^{-k_{1}z}+\left(k_2-k{\varphi }_2\right)C{e}^{k_{2}z}-\left(k_2-k{\varphi }_2\right)D{e}^{-k_{2}z}\right] \\ {S}_{zz}=+\rho \left[\left(V_{pt}^2-2V_{st}^2\right)k+k_1{\varphi }_1{V}_{pt}^2\right]\left(A{e}^{k_{1}z}+B{e}^{-k_{1}z}\right)+\rho \left[\left(V_{pt}^2-2V_{st}^2\right)k+k_2{\varphi }_2{V}_{pt}^2\right]\left(C{e}^{k_{2}z}+D{e}^{-k_{2}z}\right)-\frac{\alpha P\left[V_p^2-\left(V_{pt}^2-2V_{st}^2\right)\right]}{V_p^2} \end{gathered}$$

(7)

Once we determine the four coefficients of A, B, C, and D for each layer, we can then obtain the final displacement field in the full-space (r, z) domain by performing numerically the Hankel transforms defined above. In theory, the integrations above should be done from 0 to infinite along parameter k. In practice, however, we should define the upper limit of parameter k (denoted as k_max) so that the numerical calculation can be terminated. Throughout an auxiliary numerical experiment, it is found that k_max is inversely proportional to the total thickness of the overburden (denoted as H_ob), and the criterion of $k_{\max }\ge 10/H_{ob}$ is suggested.

It is useful to notice that the horizontal S-wave velocity (V_s) is not showing up in the final expression of the analytical solution. Therefore, if we know the vertical S-wave velocity (V_st), we need to measure the anisotropic factor only for the P-wave velocity (V_p, V_pt). It is also important to remember that the solution derived above is valid only for VTI layers, not for an exact isotropic layer (i.e., V_p = V_pt), because the linear independence among the four unknown coefficients are not valid for the exact isotropic case. For the latter, we need to use the corresponding solution for the isotropic medium, which is presented by Mehrabian and Abousleiman [16]. Therefore, if both anisotropy and isotropy exist in a subsurface media, the two solutions are used jointly. Alternatively, to simulate an isotropic layer, we may specify the ratio of horizontal-to-vertical velocities as ~1 but not exactly 1.0, e.g., 1.001, which also produces almost the same result as that by the exact isotropic layer solution.

Correction to Mehrabian and Abousleiman

In addition to the analytical solution for the VTI medium derived above, we also write down the solution for the isotropic medium derived by Mehrabian and Abousleiman [16], for the completeness of the current manuscript. It should also be noted that a few critical typos in Mehrabian and Abousleiman [16] are found through [18]. The corrected solutions for U₁, U₃, S_rz, and S_zz for the isotropic layer are given in the following form, using the two same parameters of a and c_m as in [16]:

$$\begin{array}{cc}\hfill U_1& =\frac{c_{m}P}{k}+az\left(A{e}^{kz}-B{e}^{-kz}\right)+C{e}^{kz}+D{e}^{-kz}\hfill \\ \hfill U_3& =\left(\frac{a+1}{k}-az\right)A{e}^{kz}-\left(\frac{a+1}{k}+az\right)B{e}^{-kz}-C{e}^{kz}+D{e}^{-kz}\hfill \\ \hfill S_{rz}& =2G\left[\left(akz-\frac{1}{2}\right)A{e}^{kz}+\left(akz+\frac{1}{2}\right)B{e}^{-kz}+kC{e}^{kz}-kD{e}^{-kz}\right]\hfill \\ \hfill S_{zz}& =2G\left[\left(1-akz+\frac{\nu }{1-2\nu }\right)A{e}^{kz}+\left(1+akz+\frac{\nu }{1-2\nu }\right)B{e}^{-kz}-kC{e}^{kz}-kD{e}^{-kz}-c_{m}P\right]\hfill \end{array}$$

(8)

where $a=\frac{1}{2\left(1-2\nu \right)}$, $c_m=\frac{\alpha \left(1-2\nu \right)}{2G\nu }$ (P modulus), and ν is the Poisson’s ratio. Then, following the same Hankel transform procedure, we can calculate the displacement field for an isotropic layered subsurface subjected to the same type of axis-symmetric pore pressure disturbance of finite radius and thickness.

3. Validation

In order to validate the analytical solutions for the VTI anisotropic layered medium developed in this study, we solve a set of numerical examples and compare the results with a finite element solution via COMSOL Multiphysics™ (Solids Mechanics Module, version 5.3.0.316).

Table 1. Layering and material properties for the reference isotropic model for validation. Note that we apply the velocity ratio of horizontal to vertical for the VTI anisotropic models as specified in Figure 2.

Layer.	Thickness [m]	Mass Density [kg/m3]	z-dir. S-Wave Velocity [m/s]	z-dir. P-Wave Velocity [m/s]
1	1500	1000	1500	2598
2	100	1500	1000	1732
3	∞	2000	2000	3464

shows the reference isotropic model consisting of 3 layers. Layer 2 is the reservoir where we apply a 500 m-radius pore pressure anomaly of 10 MPa (fluid injection scenario). Note that a constant Poisson’s ratio is applied to all the three layers of the reference isotropic model. For simplicity, we apply a constant horizontal-to-vertical velocity ratio not only to all the layers, but also both of P- and S-waves. However, we consider two different velocity anisotropy ratios of 0.8 and 1.2 and compare the results with those obtained for the reference isotropic model. Figure 2 shows the results by comparing the radial and vertical displacements at surface obtained from the analytical (denoted as GG) and FE solutions for the three different horizontal-to-vertical velocity ratios. First, we notice that the agreement between the analytical and FE solutions is almost perfect, which confirms that the analytical solution derived in this study and the implementation are correct. In addition, we see that the anisotropic models with ±20% horizontal-to-vertical velocity contrast (blue and yellow curves/dots) behaves quite differently in comparison to the reference isotropic model (red curves/dots). Therefore, it is suggested to apply the VTI medium model, when relevant.

4. Representation of Arbitrarily-Distributed Pressure Anomaly

The analytical solutions presented above provide the surface deformations for a simple case where the pressure anomaly is uniform and cylindrical. However, the pressure in a real reservoir is typically non-uniformly and varies in space (Figure 3a). In practice, such distribution can be described via a series of square-cuboid grids. Such a grid can be represented approximately by an equivalent cylinder (or drum) whose volume is the same as the square-cuboid grid (See Figure 3b). Then, we can describe arbitrarily-distributed pressure anomaly of any shape by applying the linear superposition of the contribution from as many square-cuboid grids as needed, for each of which we apply the analytical solution presented in the current study and using the equivalent radius (R_eq). The linear superposition can be expressed as below:

$$u_{z,i}={\displaystyle \sum }_{j=1}^J{g}_{z,ij}{p}_j$$

(9)

where u_z,i is the total vertical displacement at point i, g_z,ij is a vertical displacement at point i due to a unit-magnitude pressure anomaly at grid j (i.e., so-called fundamental or Green’s function), and p_j is the pressure magnitude at grid j. When we have the displacement at multiple points (Figure 3a), then we can also write the linear superposition expression in the following matrix form:

$$U_z=G_{z}P$$

(10)

with:

$$\begin{gathered} U_z={\left[u_{z,1},u_{z,2},\dots u_{z,i},\dots ,u_{z,I}\right]}^T, \\ P={\left[p_1,p_2,\dots ,p_j,\dots ,p_J\right]}^T, \\ {G}_z=\left[\begin{array}{ccc}{g}_{z,11}& \cdots & g_{z,1J}\\ ⋮& \ddots & ⋮\\ g_{z,I1}& \cdots & g_{z,IJ}\end{array}\right] \end{gathered}$$

(11)

Note that the linear superposition expression in matrix form can also be used for an inversion problem where we invert known (or measured) surface displacement (U_z) to estimate pressure anomaly distribution (P). In this case, it is assumed that the subsurface model (layering and stiffness or the Green’s function matrix G_z) is known.

Figure 3. (a) Schematic description of arbitrarily-distributed pressure anomaly. We assume that the pressure anomaly is discretized into J number of square-cuboid grids, to each of which a value of pressure p_j is specified. The vertical displacement at a point (say u_z,i) on the surface can be obtained by linearly summing up the contribution from all the pressure grids. (b) Illustration of an equivalent cylinder (radius R_eq = L_re/√π) to a square-cuboid of length L_re. Note both of the square-cuboid and the equivalent cylinder have the same height of H_re.

Figure 3. (a) Schematic description of arbitrarily-distributed pressure anomaly. We assume that the pressure anomaly is discretized into J number of square-cuboid grids, to each of which a value of pressure p_j is specified. The vertical displacement at a point (say u_z,i) on the surface can be obtained by linearly summing up the contribution from all the pressure grids. (b) Illustration of an equivalent cylinder (radius R_eq = L_re/√π) to a square-cuboid of length L_re. Note both of the square-cuboid and the equivalent cylinder have the same height of H_re.

Since the equivalent-radius representation is an approximation, we need to know its performance and limitation so that we can decide R_eq to be used without introducing any significant error in the calculation. To investigate this issue, we perform a parameter study by using FE simulation (COMSOL Multiphysics™, version 5.3.0.316). Namely, we compare two FE solutions obtained (1) by the exact square-cuboid shape pressure anomaly and (2) by the equivalent-radius cylinder shape pressure anomaly, respectively, with varying the lateral length of reservoir pressure anomaly (length, L_re). For the simplicity in calculation, we simulate an isotropic homogeneous half-space of 1 GPa shear modulus and 0.25 Poisson’s ratio. Figure 4 shows the results of the parameter study by plotting the relative difference between the two FE solutions versus the ratio of the radial coordinate to the overburden thickness (H_ob). Note that in the parameter study we fix the pressurized reservoir thickness (H_re) to be 100 m, while we consider four different overburden thicknesses of 500 m, 1000 m, 1500 m, and 2000 m. In addition, for each overburden thickness, we vary the ratio of L_re/H_ob from 0.5 to 2. As shown in Figure 4, the relative difference between the two FE solutions depends on all the parameters used in the current sensitivity. It is also noticed that the dependency is most influenced by the ratio of L_re/H_ob, and the cylinder-shape approximation for a square-cuboid-shape constant pressure anomaly works very well (<1% relative error) as long as L_re/H_ob ≤ 1.4–1.5 (depending on reservoir depth). For the practical purpose, we choose the condition of L_re/H_ob ≤ 1.0 (i.e., R_eq/H_ob ≤ 1/√π) in order to guarantee the accurate displacement calculation at the top surface for a square-cuboid shape pressure via the equivalent cylinder shape pressure. The purple dashed lines in Figure 4 correspond to the upper limit of the criterion (i.e., L_re/H_ob = 1.0). However, note that the real size for L_re (i.e., size of square-cuboid grids) to be used is also depending on the spatial variation of pore pressure disturbance in the subsurface e.g., sharp variation requiring smaller L_re.

5. Synthetic Model Study Inspired by in Salah CO₂ Storage Dataset

Here, we apply the linear superposition approach presented in the previous section to a realistic synthetic model and evaluate the performance by comparing with a full 3D FE solution. Furthermore, we also utilize the same superposition approach in the context of the simple inversion problem described earlier to explore the influence of the subsurface properties on the inverted pressure. For these purposes, we use an available FE model that was made for the history-matching study of the injection data and surface heave of the In Salah storage site [3]. To make things simpler, we remove the vertical fault (F12) at Well KB502, and the pressure build-up within the fault from the FE simulation. Namely, we take only the pressure build-up in the reservoir layer of 20 m thickness from the FE simulation results, and produce a synthetic surface deformation data by means of the simplified FE model. Then, using the superposition approach combined with the analytical solution derived in the present study, we invert for the reservoir pressure.

Table 2. Layering and material properties for In Salah-inspired isotropic synthetic model (after Bjørnarå et al. [3]).

Layer	Thickness [m]	Young’s Modulus [GPa]	Poisson’s Ratio [-]	Remark
1	900	3	0.25	Shallow aquifer (Cretaceous)
2	750	5	0.30	Cap rock (Visean mudstone)
3	130	2	0.30	Lower cap rock
4	20	20	0.25	Tight sandstone
5	20	9	0.15	Reservoir
6	∞	15	0.30	Devonian (underburden)

shows the layering and material properties for the In Salah model simplified from [3]. In addition, Figure 5 shows the reservoir pressure change calculated by Bjørnarå et al. [3] in a grid setting of 1000 m × 1000 m × 20 m (i.e., L_re = 1 km and H_re = 20 m for each square-cuboid grid), covering a large area of 25 km × 25 km.

First, we re-calculate the surface deformation resulting from the reservoir pressure build-up shown in Figure 5 yet now by using the linear superposition approach with the analytical solution derived earlier in order to check if the linear superposition approach produces the same results as the simplified FE model. The comparison of the results is shown in Figure 6. It is clearly seen that the two results look almost identical to each other. Note that the surface deformation in Figure 6 is plotted in a grid setting of 125 m × 125 m.

The next task is to invert the surface deformation calculated by the FE model (i.e., shown in Figure 6a) by means of the inversion idea explained earlier. In our inversion study here, we introduce some changes in the layer stiffness to see the impact of such stiffness change on the inverted reservoir pressure build-up. Namely, we increase the stiffness of each layer by 10% in the three different following scenarios. In Scenario 1, we increase the stiffnesses of Layers 1–3, while keeping the other layers’ stiffnesses unchanged. In Scenario 2, we increase the stiffnesses of only Layers 4 and 5 and in Scenario 3, for Layer 6 only. In addition, we also run the inversion without changing any stiffness as Scenario 0. The inversion results for all the four scenarios are plotted in Figure 7 in terms of inverted pressure (left panels) and difference between inverted and true (right panels). It is clearly shown that the stiffness change has noticeable impact on the inverted pressure, and the impact is the most significant (more than 1 MPa or 10%) when the stiffnesses of cap rock and reservoir are changed (Scenario 2 and shown in Figure 7e,f). This observation illustrates that it is very important to estimate the high-confidence stiffnesses of cap rock and reservoir in comparison to the other layers in the subsurface.

6. Conclusions

We have presented a generalized Geertsma solution that can consider any number of layers in the subsurface whose properties and thicknesses can be different from layer to layer. In addition, each layer can be VTI anisotropic, with which more realistic stiffness can be simulated depending on condition of in situ stress or sediment bedding plane. The accuracy of the newly-derived analytical solution is validated with respect to an FE-based reference solution. Furthermore, it is found that the Hankel transform to be numerically performed can be done via a definite integral whose upper limit can be determined optimally as inversely proportional to the overburden thickness (i.e., $k_{\max }\ge 10/H_{ob}$). The generalized Geertsma solution is then applied to a linear superposition framework via square-cuboid-shaped grids so that we can calculate arbitrarily-distributed pressure anomaly cases, not only the simple case of a constant pressure anomaly of cylinder shape. Furthermore, we have also defined the condition that the lateral length of a square-cuboid should be shorter than the overburden thickness (i.e., L_re/H_ob ≤ 1.0) in order to guarantee the accurate displacement calculation at the top surface for a square-cuboid shape pressure via the equivalent cylinder shape pressure. The performance of this linear superposition approach is tested by solving a realistic synthetic model based on the In Salah CO₂ storage site, simplified by removing the vertical fault in Well KB502. Finally, we look at the influence of subsurface stiffness on inverted pressure by means of a simple inversion based on the linear superposition approach. It is learned that the stiffnesses of cap rock and reservoir are the most influencing parameters on the inversion result, suggesting that it is very important to estimate the mechanical properties of cap rock and reservoir with high accuracy (i.e., low uncertainty).

The next step to follow up among many possibilities would be that we apply the solution approach presented in the current study to real data inversion of the In Salah storage with taking into account known fault structure such as F12. For this purpose, we may have to improve a limitation of the current approach that the subsurface should be a horizontally-layered medium and cannot consider any tilted structure (e.g., faults with different properties from layers). Tackling this limitation is our immediate task to pursue. Finally, since the analytical solution approach is fast, its application to machine learning-based inversion is also a desired task to follow up.