A Three-Dimensional Analytical Solution for Reservoir Expansion, Surface Uplift and Caprock Stress Due to Pressurized Reservoirs

An analytical solution is presented for the displacement, strain and stress of a three-dimensional poro-elastic model with three layers, where the three layers are an underburden, a reservoir with a given fluid pressure, and an overburden. The fluid pressure in the reservoir is assumed symmetrical around the z-axis and represented by a Fourier cosine series. The poro-elastic solution is expressed as a superposition of the solutions for each term in the Fourier series. It is shown that the bulk strain in the reservoir layer is proportional to the fluid pressure and that the bulk strain in the underburden and overburden is zero. Using these properties of the bulk strain, a solution is derived for the three-layer model where the fluid flow and mechanics are fully coupled. A particular aim of the model is to study the surface uplift from a given reservoir pressure. The expansion of the reservoir and the uplift of the surface are studied in terms of the wavelengths in the Fourier representation of the pressure. It is shown that the surface uplift can be written in a similar form to the 1D vertical expansion of the reservoir layer, but where the fluid pressure is based on the Fourier series. It is shown that the amplitudes with average wavelengths longer than 2π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\pi $$\end{document} times the thickness of the reservoir give expansion of the reservoir, but average wavelengths much shorter than this limit do not. Similarly, amplitudes with average wavelengths much longer than 2π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\pi $$\end{document} times the thickness of the overburden produce surface uplift, but wavelengths much shorter do not. The stress in the overburden, which is generated by the reservoir fluid pressure, is also analysed in terms of the wavelengths. A case is given where the analytical uplift is compared with the results of a numerical simulation and the agreement is excellent.


Introduction
The global annual emissions of CO 2 were 32 Gt in 2016 (International Energy Agency 2016), and the concentration of CO 2 in the atmosphere has increased from around 280 ppm in 1900 to more than 400 ppm in 2017. CO 2 is a greenhouse gas and its increasing concentration in the atmosphere has been found to be the reason for global warming (Bryant 1997). The injection of large volumes of CO 2 into aquifers or depleted oil and gas reservoirs is considered a promising way to reduce CO 2 emissions and thereby reduce global warming (Bachu 2008;Benson and Cole 2008;Bickle 2009). CO 2 injection leads to a pressure build-up and an expansion of the reservoir, and in turn to surface uplift. An example of subsurface CO 2 injection and surface uplift is the In Salah storage project, where 3.8 Mt has been injected since its start in 2004, and where a total uplift of roughly 15 mm has been measured around three injection wells (Rutqvist 2012). The uplift measured at In Salah has attracted a considerable amount of scientific interest, and a number of modelling studies have been published (Vasco et al. 2008;Rutqvist et al. 2010;Rutqvist 2012;Bissell et al. 2011;Verdon et al. 2011;Zhou et al. 2010;Shi et al. 2012;Rinaldi and Rutqvist 2013;Rucci et al. 2013;White et al. 2014;Vilarrasa et al. 2015Vilarrasa et al. , 2017Rinaldi et al. 2017). Similarly, seabed uplift is expected over offshore CO 2 sites, although it is more difficult to measure than land uplift. Land-based water pumping from aquifers and land-based oil and gas production are known to produce surface subsidence as the fluid pressure is lowered. Geertsma (1973) developed an early model to compute land subsidence, which was applied to the Groningen gas field. The model gives the poro-elastic subsidence above a horizontal disk-shaped reservoir with a constant thickness and a constant pressure. The reservoir is placed at a finite depth in the infinite half-space and has a stress-free surface. The model is based on the nucleus-of-strain concept in the half-space, which was introduced by Mindlin and Cheng (1950) and Sen (1950) as a method to solve thermoelastic problems. A particular feature of the model of Geertsma (1973) is that a constant pressure reduction gives a maximum subsidence that is approximately 1.5 times larger than the one-dimensional (vertical) compaction of the reservoir layer. The model gives a maximum horizontal displacement (contraction) that is also approximately 1.5 times the one-dimensional compaction.
The classical model of Geertsma (1973) has recently been modified by Tempone et al. (2010): the infinite half-space was replaced by a rigid basement underneath the reservoir. Tempone et al. (2010) compare the displacement, stress and strain of their model with the classical Geertsma model (Geertsma 1973). The comparisons show that a rigid basement close to the reservoir gives more uplift than the classical model, because less of the reservoir compaction goes vertically downward into the layer underneath the reservoir.
Another approach is taken by Fokker and Orlic (2006). They model the poro-elastic subsidence by a semi-analytical approach, which allows for an arbitrary number of layers. They use a numerical method to fit a set of analytical solutions to approximate the boundary condition in an optimal manner. Selvadurai (2009), Kim and Selvadurai (2015), Selvadurai and Kim (2016) and Niu et al. (2017) have developed analytical models for different configurations of a reservoir and a caprock. Selvadurai and Kim (2016) present analytical poro-elastic solutions for a storage aquifer with a caprock; the solutions are given by an integral representation. The model was used to study how the surface displacement is controlled by the radius and the depth of the injection region. Finally, Chang and Segall (2016) studied how fluid pressure variations from injection or production in reservoirs could induce poroelastic stress changes and fault slip in the basement.
In the present article, a different approach is taken, by first solving the poro-elastic equations for a three-dimensional model with three layers when the reservoir pressure is just one term in a Fourier series. The same approach was taken by Wangen et al. (2018) with a two-dimensional model of a reservoir with an overburden, but no underburden. The three layers in the present model are an overburden, a reservoir and an underburden, where the underburden is placed on a rigid basement. Each layer in the model is of infinite lateral extent and of constant thickness. A periodic pressure distribution can be represented by a Fourier series, and the displacement, strain and stress becomes the superposition of the effect of each term in the Fourier series. The model has a stress-free surface, except for the shear stress xz and yz components, which turn out to be negligible for normal configurations of the reservoir and the overburden. The displacement field and the stress, except for the x x and yy stress components, are continuous across the internal layer boundaries, which are the top and the base of the reservoir. A particular feature of this model is that it allows a decoupling of the equation for fluid flow from the equations of the displacement field. This decoupling simplifies considerably the derivation of analytical solutions of the fully coupled problem of fluid flow and mechanics. This model is well suited for a study of how the displacement, stress and uplift depend on the wavelengths of the reservoir pressure and the layer thicknesses. 'Short' and 'long' wavelengths in the pressure distribution produce different displacement fields and stresses, where 'short' and 'long' are with respect to layer thicknesses.
The present paper is organized as follows: The poro-elastic assumption and notation are introduced, and the expansion of the one-dimensional vertical poro-elastic layer is established. The geometry of the three-layer model is explained, and the boundary conditions are discussed. The solution for the displacement field is presented, and then the limit as the wavelength goes to infinity is investigated. How the fluid pressure can be decomposed in a double Fourier series is demonstrated before a Fourier solution for the fully coupled problem of fluid flow and mechanics is derived. The uplift is studied in terms of the wavelengths of the Fourier decomposition, and a special solution for the displacement for a reservoir on top of a rigid basement is presented. The analytical solution for the uplift is compared with a numerical finite element solution. An example of the analytical solution is given, which demonstrates the displacement field and the stress field. The stress in the caprock is investigated with respect to the wavelength of the Fourier components in the reservoir pressure.

Poro-elasticity
The stress state σ i j in the model is given by the sum of the initial stress σ (0) i j and the change σ (1) i j in the poro-elastic stress caused by a change in the fluid pressure in the reservoir where the indices i and j indicate the three spatial dimensions x = 1, y = 2 and z = 3. The stress state σ i j fulfills the equilibrium equations where ρ is the rock bulk density, g is the acceleration due to gravity, and δ iz is the Kronecker delta, that is, δ iz = 1 for i = z and δ iz = 0 for i = z. Einstein's summation convention is used in the equilibrium Eq.
(2), where summation is understood for each pair of the same indices.
The initial stress state also fulfills the equilibrium Eq.
(2), which implies that the stress produced by poro-elastic changes fulfills an equilibrium equation without body forces The effective stress can also be decomposed into the initial effective stress and the effective stress produced by changes in the fluid pressure in the reservoir The initial effective stress is given by where p (0) is the initial fluid pressure and α is the Biot coefficient. Similarly, the difference in the poro-elastic effective stress is given by where p (1) is the change in the reservoir pressure. The change τ (1) i j in the effective stress produces poro-elastic deformations where Λ is the Lamé parameter, G is the shear modulus, ε i j = 1 2 (u i, j + u j,i ) is the strain tensor, and u i is the displacement vector (i = x, y, z). An upper case Λ is used for the Lamé parameter, since the lower case λ will be used for the wavelength. Notice that normal effective stress is positive when the rock is tensile and negative when it is compressive. The elastic moduli are for drained conditions. The equations for the displacement field u i are obtained by inserting the stress from Eq. (7) into the equilibrium Eq. (3), which yields where γ = Λ/G. For simplicity, the layers in the system are taken to have the same geomechanical parameters Λ, G and α. The gradient of the pressure change appears as an internal load on the left-hand side of the equilibrium equations. Therefore, Eq. (8) gives the displacement field when the pressure is known. The displacement field has a feedback on the pressure. The equation for the fluid pressure that includes poro-elastic deformations is where K f is the modulus of fluid compressibility, K s is the modulus of the solid constituting the rock matrix, φ is the porosity, κ is the permeability, μ is the pore fluid viscosity, ε = ∂u i /∂ x i is the bulk strain and t is the time. The notation is simplified by denoting the pressure change by p = p (0) . The notation for the displacement field is also simplified by writing u = u x , v = u y and w = u z . For instance, the bulk strain can now be written ε = ∂u i /∂ x i = ∂u/∂ x + ∂v/∂y + ∂w/∂z, where it is noted that summation is understood over pairs of equal indices. In the following, the fluid overpressure is assumed known and the displacement equations are solved with the given overpressure.

One-dimensional Vertical Uplift
The one-dimensional vertical poro-elastic expansion of a reservoir layer due to a constant increase p 0 in the overpressure is an important reference model. The equations for the displacement field (8) become for one-dimensional vertical displacements, when u = 0 and v = 0. The vertical displacement w(z) is obtained by integrating Eq. (10) twice where it is assumed that the base of the reservoir layer at z = z 0 is fixed (w = 0). The surface of the reservoir layer has a maximum displacement given by Fig. 1 The three-layer model and the coefficients of the homogeneous solution for each layer. The base layer is at z = 0, the base of the reservoir is at z = z 1 , the top of the reservoir is at z = z 2 and the top surface is at z = z 3 where h is the thickness of the reservoir layer. In case the reservoir layer has an overburden, the overburden is uplifted by w when the base is fixed. The one-dimensional reservoir expansion (12) turns out to be an important reference for the uplift.

The Three-Layer Model
The model consists of three layers of infinite extent, as shown in Fig. 1. The middle layer is the reservoir and it has an overburden and an underburden. To simplify the analytical solution for the displacement field, the three layers are assigned the same mechanical properties. Only the reservoir layer is considered permeable and has a change in the fluid pressure. The z-axis is positive upwards, and z = 0 is the base of the model (the underburden). The base of the reservoir layer is at z = z 1 , the top of the reservoir layer is at z = z 2 and the top of the model (the surface) is at z = z 3 . The thicknesses of the underburden, the reservoir and the overburden are h 1 , h 2 and h 3 , respectively, and the layer boundaries become z 1 = h 1 , z 2 = h 1 + h 2 and z 3 = h 1 + h 2 + h 3 .

Boundary Conditions
The layers are of infinite extent, and the boundary conditions are at the base and at the top surface of the model. The base surface has zero vertical displacement, and the top surface has zero lateral displacement. These boundary conditions can be written as The boundary condition w(z = 0) = 0 implies that the underburden is resting on an absolutely rigid base. The top surface may move vertically, but it is prevented from lateral displacements, since it is assumed that u(z = z 3 ) = 0 and v(z = z 3 ) = 0. It will be seen that this boundary condition leads to an approximately stress-free top surface when the overburden is much thicker than the reservoir. This approximation works fine for 'long' wavelengths in the reservoir pressure, because they give locally almost one-dimensional vertical surface deformations. Furthermore, it is seen that 'short' wavelengths in the reservoir pressure produce deformations that decrease to almost zero towards the surface.

The Solution for the Displacement Field
The displacement field is first found for a periodic reservoir fluid overpressure of the form which is one term in a Fourier representation of a fluid pressure that is symmetrical around the z-axis. The overpressure (14) is independent of the depth inside the reservoir layer. The amplitude of the fluid pressure is p 0 , and the wavenumbers in the xand y-directions are k 1 and k 2 , respectively. These wavenumbers correspond to the wavelengths respectively. The full solution for the displacement field of the three-layer model consists of three parts: one for each of the three layers comprising the underburden, the reservoir and the overburden. The solution gives a displacement field that is continuous across the layer boundaries. The derivation of the solution is presented in the "Appendix". The underburden (0 ≤ z ≤ z 1 ) has the following solution for the displacement field where the average wavenumber is and where and are two coefficients. The reservoir (z 1 ≤ z ≤ z 2 ) has the displacement field where and are two functions of z. Finally, the overburden (z 2 ≤ z ≤ z 3 ) has the solution where is a coefficient. It is seen that all parts of the displacement field are proportional to the coefficient D 0 and the pressure amplitude p 0 . Another observation is that the functions F 2 (z) and F 3 (z) have the following properties dF 2 dz = cF 3 (z) and dF 3 dz = cF 2 (z).
These relations are useful in the expressions for the strain and for verifying that the stress fulfills the equilibrium equations. The strain ε i j = 1 2 (u i, j +u j,i ) is straightforward to compute from the displacement field, and the strain tensor for each layer is given in the "Appendix". The effective stress follows from the strain by Eq. (7) and the effective stress tensor τ (1) i j for each layer is also given in the "Appendix". The stress σ (1) i j is equal to the effective stress in the overburden and the underburden where the fluid overpressure is zero. In the reservoir layer, the normal stress is obtained from the effective stress by subtracting the reservoir pressure multiplied by Biot's coefficient; see Eq. (14). Finally, it is straightforward to verify that the stress σ (1) i j fulfills the equilibrium Eq. (3). It is also straightforward to verify that the displacement field (u, v, w) and the stress components σ xy are continuous across the internal boundaries at the base and at the top of the reservoir. In the same way, it can be seen that the stress components σ (1) yy are discontinuous across the horizontal layer interfaces. The stress is σ (1) x x = σ (1) yy = ΛD 0 p 0 cos(k 1 x) cos(k 2 y) larger in the reservoir than right above in the caprock or right below in the underburden. The following two relations are useful for showing the continuity of the stress field across the internal layers where the first relation applies to the base of the reservoir and the second to the top of the reservoir.
7 The Limit as k 1 → 0 and k 2 → 0 An important aspect of the solution (16)-(29) is what happens in the limit as both k 1 → 0 and k 2 → 0. This is the limit of infinite wavelengths in both the x-and ydirections, where the pressure (14) becomes laterally constant with an amplitude p 0 . Notice that the four functions F 1 to F 4 depend on the arguments cz 1 , cz 2 and cz 3 , and that c also goes to zero when both k 1 and k 2 do so. The limit k 1 → 0 and k 2 → 0 can be studied by using k 1 = k 2 = k and c = √ 2k, and letting the wavenumber k → 0. A fixed value of x gives that and similarly for fixed values of y and z. This implies that the factor k/c 2 disappears in all expressions for u (x, y, z) and v(x, y, z), which are Eqs.
and (28). Next, it is seen that which implies that the lateral displacement field goes to zero in the limit k 1 , k 2 → 0. This is as expected for a constant pressure distribution, because it is a one-dimensional vertical model. The same argument as above implies that the vertical displacement field of the underburden goes to zero in the limit k 1 , k 2 → 0. The vertical displacement field of the reservoir depends on the function F 3 in the limit c → 0, and it can be seen that F 3 → (cz − cz 2 ). Therefore, the vertical displacement field of the reservoir (24) has the limit w → D 0 p 0 · (z − z 1 ) when k 1 , k 2 → 0, which is the solution for the one-dimensional vertical reservoir expansion.
In the overburden, the vertical displacement field depends on the function F 4 , which has the limit (cz 2 − cz 1 ) when c → 0. Therefore, the vertical displacement of the overburden has the limit w → D 0 p 0 · (z 2 − z 1 ), which is, as expected, the maximum vertical expansion of the reservoir layer.
Furthermore, it follows that the strain and the effective stress tensors are zero in the limit k 1 , k 2 → 0 for all layers, with exceptions for the strain ε zz = αp 0 /(Λ + 2G) and effective stress τ (1) zz = αp 0 for the reservoir layer.

Fourier Representation of the Reservoir Pressure
The reservoir pressure can be represented as a double Fourier cosine series, assuming that the pressure is symmetrical around the z-axis A 0,n cos(k n y) where the Fourier coefficients are given by and where the wavenumbers for indices m and n are respectively; see Kreyszig (2011). The domain sizes in the x-and y-directions are L 1 and L 2 , respectively. A Fourier series gives an exact representation of any periodic and well-behaved function, with periods L 1 and L 2 in the x-and y-directions, respectively, in the limit N 1 → ∞ and N 2 → ∞. The domain size is assumed sufficiently large for the pressure to go to zero before the edges of the domain. A finite representation, where N 1 and N 2 are of order 100, normally gives an excellent approximation of the reservoir overpressure. It is convenient to write the series (40) compactly as where the coefficients A mn are given by The displacement field is linear in the reservoir overpressure, which also implies that the strain and the stress are linear in the overpressure. Since the model is linear, the displacement, stress and strain from a pressure distribution written as a Fourier series become the superposition of the displacement, stress and strain from each term in the Fourier series. For instance, if S(x, y, z) is a property such as the displacement, the stress or the strain, it implies that where S (k m x, k n y, cz) is the contribution to the property from one Fourier component, and where p 0 is replaced by the pressure A m,n . The function S (k m x, k n y, cz) is zero for m = 0 and n = 0, except for a displacement w in the vertical direction, the strain ε zz and the effective stress τ (1) zz , as already shown.

Solution of the Coupled Biot Equations
The strain tensor in the underburden, reservoir and the overburden is straightforward to compute from the displacement field; see the Appendix. It may be seen from the Appendix that the bulk strain is zero in the underburden and the overburden. On the other hand, the bulk strain in the reservoir layer is proportional to the pressure, ε = D 0 p, when the pressure is represented by a Fourier series. This allows for a decoupling of the pressure Eq. (9) from the displacement Eq. (8). The pressure equation for the reservoir layer can now be written as where is the characteristic time associated with the average wave length λ = 2π/(k 2 m +k 2 n ) 1/2 . The Fourier coefficients in Eq. (46) give the initial condition, which decays to zero with time. The solution (46) provides a Fourier series that solves for the case of fully coupled fluid pressure and deformations.

Uplift as a Function of Wavelengths and Layer Thicknesses
A particular aim with this model is to study which wavelengths (wavenumbers) in a Fourier series of the reservoir pressure give surface uplift. Therefore, the vertical deformations are studied for one term in the Fourier series, as represented by the fluid pressure (14). The surface uplift is given by the vertical displacement (29) for z = z 3 , and it can be written as and where The first of the three factors in the uplift (48) is the one-dimensional reservoir expansion produced by a constant reservoir pressure p 0 . The second factor, f surf , is a dimensionless surface uplift: it turns out to be a number between 0 and 1. The dimensionless function depends on the three parameters ch 1 , ch 2 and ch 3 , where c = (k 2 1 + k 2 2 ) 1/2 and λ = 2π c .
are the average wavenumber and the average wavelength, respectively. The dimensionless surface uplift f surf describes how the uplift depends on the layer thicknesses and the average wavenumber, or, alternatively, the average wavelength. Figure 2 shows the dimensionless uplift as a function of ch 1 and ch 2 for the three different underburden thicknesses given as ch 1 = 0.1, 1 and 10. For all three cases, the uplift is close to one for ch 2 1 and ch 3 1. Another way to state this condition is that the average wavelength must be larger than both the reservoir thickness and the overburden thickness multiplied by 2π , i.e. λ 2π h 2 and λ 2π h 3 . This is a generalization of the condition for uplift found with the two-dimensional plain-strain model of Wangen et al. (2018), which has a reservoir layer and an overburden, but no underburden. The dimensionless uplift (49) gives exactly the same dimensionless uplift as in Wangen et al. (2018) under plain strain conditions (k 2 = 0) and with an underburden of zero thickness (h 1 = 0).  (29) at z = z 2 , because the displacement field is continuous at the layer interfaces, and it can be written as and where w 0 = cosh(ch 3 ) − cosh(ch 2 + ch 3 ) cosh(ch 1 + ch 2 + ch 3 ) .
(53) Figure 3 shows the dimensionless uplift of the top reservoir, Eq. (52). It lies between 1/2 and 1 for ch 2 1, or alternatively λ 2π h 2 . When ch 2 1 and ch 3 1, the dimensionless uplift is ≈ 1, and when ch 2 1 and ch 3 1, it is ≈ 1/2, unless the overburden is 'thin', ch 1 1, in which case it is ≈ 1. In the opposite regime, ch 2 1 (λ 2π h 2 ), it is almost zero. Wavelengths much shorter than 2π h 2 do not produce an uplift of the top of the reservoir. The vertical movement of the base of the reservoir is given by Eqs. (18) and (24) . 4 The dimensionless base reservoir uplift as a function of ch 1 , ch 2 and ch 3 where the dimensionless vertical displacement for the base of the reservoir is The dimensionless displacement f base is plotted in Fig. 4, which shows that the base of the reservoir has a negative displacement. This means that it is pressed down into the underburden. Figure 4a shows that the conditions ch 2 1 and ch 3 1 give small negative displacements of the base reservoir. Furthermore, it can be seen from Fig. 4 that the maximum downward displacement of the base of the reservoir takes place for ch 2 1 and ch 1 , ch 3 1. This condition can also be written as λ 2π h 2 and λ 2π h 1 , 2π h 3 . The condition says that the average wavelength must be 'long' compared to the thickness of the reservoir and at the same time 'short' compared with the thickness of the underburden and the overburden.
Finally, the dimensionless expansion of the reservoir is given by the difference between the vertical displacements of its top and base The dimensionless expansion is plotted in Fig. 5. It can be seen that the condition for expansion of the reservoir is that ch 2 1 or, alternatively, that λ 2π h 2 , regardless  Fig. 2 shows that the expansion of the reservoir reaches the surface only when ch 3 1 or λ 2π h 3 . Since the overburden is normally much thicker than the reservoir, h 3 h 2 , the condition for surface uplift becomes λ 2π h 3 , in terms of the average wavelength.

Reservoir and Overburden Above a Rigid Basement
The model is simplified considerably if it is assumed that the reservoir rests on an absolutely rigid basement. Then, the underburden thickness can be set to zero (h 1 = 0). In this case, the displacement field for the reservoir layer (0 ≤ z ≤ z 2 ) becomes and for the overburden (z 2 ≤ z ≤ z 3 ) The vertical displacement (58) gives the dimensionless uplift for this case: which becomes the same as the dimensionless uplift reported for the two-layer plainstrain model by Wangen et al. (2018) by setting either k 1 = 0 or k 2 = 0. It should be noted that the dimensionless uplift is given by a different expression by Wangen et al. (2018), but it is easily rewritten to be as in Eq. (59).
where Q is the injection rate, h 2 is, as before, the thickness of the aquifer, r = (x 2 + y 2 ) 1/2 is the radius from the injection well and t is the time. The example has a reservoir thickness h 2 = 50 m, porosity φ = 0.1, fluid compressibility 1/K f = 5 × 10 −8 Pa −1 , solid matrix compressibility 1/K s = 0 Pa −1 , aquifer permeability κ = 1 × 10 −12 m 2 and injection rate Q = 0.278 m 3 s −1 . The overburden thickness is h 3 = 1000 m, and the reservoir layer is placed on a rigid basement, which gives h 1 = 0 m. Young's modulus is E = 15 GPa, and the Poisson ratio is ν = 0.2.
Since E K f , the fluid compressibility dominates the compressibility of the rock, and the right-hand side of Eq. (9) can be ignored.
How the pressure plume spreads out from the origin is plotted in Fig. 6a by showing the pressure at 10, 100 and 1000 days. In Fig. 6a, the solid curve shows the Theis solution, and the markers show the pressure computed from the Fourier series. The uplift for one Fourier component is expressed by Eq. (48). Since the model is linear, the uplift becomes the superposition of uplift from each term in the Fourier series. Therefore, the uplift from the pressure plume can be expressed as where p up (x, y), the uplift pressure, is defined by x-coordinate [m] The surface uplift (61) has the same form as the one-dimensional reservoir expansion (12) when it is expressed using the uplift pressure (62). The uplift pressure is the Fourier representation of the reservoir pressure filtered with the factor f surf , which is a number between 0 and 1. Figure 6b shows the uplift pressure corresponding to the pressure plume in Fig. 6a. The uplift pressure is seen have roughly the same width as the corresponding pressure plume, but it has much less height. These observations are explained by the factor f surf , which reduces 'short' wavelengths in the Fourier series. The Fourier coefficients as a function of the average Fourier number, n av = (m 2 + n 2 ) 1/2 , are plotted in Fig. 7a. Figure 7a shows that the first Fourier coefficients become larger as the plume spreads out. The first three Fourier coefficients for the pressure plume at 1000 days are considerably larger than the others. Figure 7a shows that the largest Fourier coefficients are found for an average Fourier number less than 2. These Fourier numbers correspond to an average wavelength that is longer than 10 km, as seen from Fig. 7b. The pressure plumes at times 10 and 100 days are not wide enough to have noticeable Fourier coefficients for the lowest Fourier numbers.
The Fourier coefficients are multiplied with the f surf -function, which lies in the range 0.8-1 for Fourier numbers up to 2. For Fourier numbers larger than 4, the f surffunction is less than 0.5, as it converges towards 0. The large Fourier coefficients of large wavelengths, which dominate the pressure plume at time 1000 days, are not reduced very much by the f surf -function, and therefore they contribute to the surface uplift.
The actual uplift at the three times is plotted in Fig. 8a-c. These curves show the uplift computed by expression (62) (the red circles) and the results of a finite element computation (the blue curve). The agreement between the analytical expression and the numerical FE computation is excellent. Figure 8a-c also show the one-dimensional reservoir expansion (12), when applied locally, as the green curve. The uplift from the pressure plume approaches the one-dimensional reservoir expansion as it gets wider  Fig. 8 a Surface uplift at time 10 days. b Surface uplift at time 100 days. c Surface uplift at time 1000 days and therefore has Fourier components with longer wavelengths. The one-dimensional reservoir expansion becomes a reasonable approximation for the surface uplift when the width of the plume exceeds 2π h 3 ≈ 6000 m.
Expressions similar to those for the surface uplift (61), where the uplift has the same form as the one-dimensional reservoir expansion, can be made for the displacement at the base of the reservoir and the top of the reservoir.

An Example of Displacement, Strain and Stress from Reservoir Pressure
Another version of the preceding example is shown in Fig. 9. The reservoir layer with the reservoir pressure and the overburden is the same, but a 2000 m underburden is added. Figure 9 shows a two-dimensional xz-cross section of the model through the z-axis. The fluid overpressure is limited to the reservoir layer, as seen from Fig. 9a, and Fig. 9b shows that the displacement field u is symmetric around the origin. The displacement field u is also restricted to the surroundings of the reservoir layer.  Figure 10 shows the stress corresponding to the overpressure and the displacement field in Fig. 9. Recall that the stress is the effective stress minus the pore pressure, Eq. (6), and that it is the effective stress that is responsible for deformations, Eq. (7). A negative value of σ (1) x x and σ (1) zz indicates compression. A positive pore pressure acts to expand the rock, but there is a compressive stress that holds the expansion back. Otherwise, the rock would have been maximally expanded and stress-free. The deformations are much larger in the vertical direction than in the x-direction, which implies that there is less compressive stress in the z-direction than in the x-direction. The shear stress σ (1) xz is located near the centre of the pressure plume in the underburden, the reservoir and the overburden. The absolute value of the shear stress is also two orders of magnitude less than the overpressure.

Stress in the Overburden
The surface of the model (z = z 3 ) is stress-free, except for the shear stress components σ (1) xz and σ (1) yz , as seen from Eqs. (178)-(182). It turns out that the coefficient F 4 can be used to obtain a condition for when these two shear stress components are close to zero. The shear stress σ (1) xz at the surface (z = z 3 ) can be constrained in the following way, since k 1 /c ≤ 1, α ≤ 1 and γ > 0. How F 4 = F 4 (ch 1 , ch 2 , ch 3 ) depends on ch 1 , ch 2 and ch 3 is shown in Fig. 11, where it can be seen that the shear stress σ (1) xz at the surface is negligible if one of following two conditions is fulfilled. The first is that ch 3 1, and the second is that ch 2 1 and ch 3 ch 2 . If it is assumed that the overburden is much thicker than the reservoir layer, the second condition becomes just ch 2 1. These two inequalities can be rewritten in terms of the average wavelength The coefficient F 4 is plotted as a function of ch 2 and ch 3 for four values of ch 1 as λ 2π h 3 and 2π h 2 λ. Since the overburden is assumed thicker than the reservoir layer, these intervals overlap, and the condition for negligible shear stress at the surface is fulfilled for all wavelengths. Exactly the same argument shows that σ (1) yz is close to zero at the surface.
The stress in the caprock, Eqs. (178)- (183), shows that all the components of the stress tensor are proportional to the coefficient F 4 . The components of the stress tensor depend on the z-position in the caprock by either the factor cosh(cz 3 − cz) or sinh(cz 3 − cz). This implies that the stress increases nearly exponentially, from the surface towards the base of the overburden, when ch 3 1, because 0 ≤ cz 3 − cz ≤ ch 3 . In the other regime, where ch 3 1, these two factors can be approximated by cosh(cz 3 −cz) ≈ 1 and sinh(cz 3 −cz) ≈ 0. It can therefore be concluded that the short wavelengths in the fluid pressure, ch 3 1 or λ 2π h 3 , produce stress in the caprock that is maximal towards its base. Long wavelengths, ch 3 1 or λ 2π h 3 , produce a weak stress in the caprock for normal basin configurations where the overburden is much thicker than the reservoir (h 2 h 3 ).

Conclusions
An analytical solution is presented for the displacement, strain and stress of the poroelastic equations for a three-dimensional model of three horizontal layers with the same rock properties. The three layers are an underburden on a rigid basement, a reservoir with a pore pressure change, and an overburden. The overburden and the underburden are assumed impermeable with no fluid pressure change. The boundary conditions are zero vertical displacement at the base of the underburden and zero horizontal displacement at the top surface. It has been shown that the top surface is stress-free, except for the shear stress components xz and yz, which are almost zero for common basin configurations where the overburden is much thicker than the reservoir. The pressure distribution in the reservoir is assumed symmetrical around the zaxis, in which case the pressure can be represented by a double Fourier cosine series. The model provides a solution for the displacement, the strain and the stress as a superposition of contributions from each term in the Fourier series. It is shown that the bulk strain in the overburden and the underburden is zero, and that the bulk strain in the reservoir layer is proportional to the fluid overpressure. By means of this property of the bulk strain, a Fourier series solution is derived where the fluid flow and mechanics are fully coupled.
The contribution from each term in the Fourier series is proportional to the amplitude and a factor of proportionality that is a function of the average wavelength and the layer thicknesses. These factors provide insight into how the fluid pressure deforms the layers with respect to wavelengths and layer thicknesses. It is seen that only terms in the Fourier series with average wavelengths larger than 2π times the thickness of the reservoir (λ 2π h 2 ) produce an expansion of the reservoir. A second condition must be fulfilled for the expansion of the reservoir to produce a surface uplift, which is that the average wavelength must be larger than 2π times the thickness of the overburden (λ 2π h 3 ). When the second condition is not fulfilled, the expansion of the reservoir does not reach the surface. The expansion of the reservoir can also press down the underburden in this case.
The surface uplift can be written in the same form as the one-dimensional vertical reservoir expansion, where the fluid pressure is based on the Fourier series for the reservoir pressure. This expression shows that the surface uplift in this model can never be larger than the expansion of the reservoir. The analytical expressions for the surface uplift are compared with the numerically computed surface uplift for a case where the reservoir is resting on a rigid basement, and the match is excellent.
The model predicts a discontinuity in the x x and yy components of the stress tensor at the top and the base of the reservoir. Furthermore, 'short' wavelengths (λ 2π h 3 ) in the Fourier series for the fluid pressure produce a stress at the base of the overburden with nearly the same strength as the amplitude in the Fourier series, but the stress does not extend far upwards into the overburden. Terms with 'long' average wavelengths (λ 2π h 3 ) produce negligible stress in the overburden. It remains to study how the model eventually can be extended to layers with different mechanical properties.
The system of linear Eq. (73) has a nonzero solution for the vector [u 1 , v 1 , w 1 ] when det(A) = 0. If the determinant of A is zero, then and the linear Eq. (73) implies that the coefficients u 1 , v 1 and w 1 are related by The same procedure can be applied for a homogeneous solution of the displacement field of the alternative form u = u 2 exp(−cz) sin(k 1 x) cos(k 2 y), (76) v = v 2 exp(−cz) cos(k 1 x) sin(k 2 y), (77) w = w 2 exp(−cz) cos(k 1 x) cos(k 2 y), where det(A) = 0 once more implies c 2 = k 2 1 + k 2 2 , and where now implies that k 1 u 2 + k 2 v 2 − cw 2 = 0.
These two pairs of expressions for w 3 and w 4 are useful in deriving these coefficients from either w 0 or w 5 . Finally, these two alternative expressions give that The coefficients w 0 and w 5 give w 3 and w 4 from the pair of equations (124)-(125) or (126)-(127). Then, the coefficients u i and v i for i = 1, 2, 3, 4 are obtained from w 3 and w 4 using relations (118)-(121). The coefficients u 5 and v 5 are found from interface relations (107) and (107), respectively, and finally, u 6 , v 6 and w 6 are obtained from Eqs. (104)-(106), respectively.