Leakage From Coexisting Geologic Forcing and Injection‐Induced Pressurization: A Semi‐Analytical Solution for Multilayered Aquifers With Multiple Wells

Abnormal fluid pressures (above or below hydrostatic pressure) can develop and persist in sedimentary basins. The common occurrence of abnormal pressures may cause challenges for project permitting of geological carbon sequestration (GCS), particularly in reservoirs with pre‐injection overpressure. The leaky wells that may exist in some sedimentary basins can provide flow paths between deep brine aquifers and shallower freshwater aquifers. Pre‐injection relative overpressures can cause brine leakage through leaky wells even before any injection occurs. The tendency for flow through leaky wells is coupled with the process of pressure dissipation that occurs through aquitards. Specifically, with non‐zero permeability, aquitards can dissipate pressure over large areal extents and thereby reduce leakage rates through leaky wells. This study presents development of a semi‐analytical solution for hydraulic head and brine leakage in multilayered aquifer–aquitard systems with geologic pressure forcing. The geologic forcing that causes abnormal pressures in the multilayered system can coexist with any number of injection, extraction, and leaky wells that also affect fluid pressure. The semi‐analytical model is applied to explore how leakage through leaky wells varies as functions of pressurization rate, along with aquitard and leaky well properties in an overpressured multilayered system. The results show that although injection‐induced pressures can dissipate rapidly through suitably permeable aquitards, coexisting geologic forcing may create sustained rates of brine leakage into freshwater aquifers through leaky wells. In GCS, a very low‐permeability aquitard with high capillary entry pressure to free‐phase CO2 is desired to serve as the caprock to prevent leakage of CO2 from the storage reservoir. Nevertheless, the results from this study show that the brine leakage impact to shallow freshwater aquifers through leaky wells might substantially decrease with increasing aquitard permeability values, as long as small aquitard permeability and high capillary entry pressure serve to prevent CO2 leakage.

Abandoned wells (e.g., unused/unsuccessful exploration wells, or improperly plugged water or oil/gas wells) that are not sealing and/or whose sealing components have degraded over time can create hydraulic communication through caprocks (aquitards), thereby creating flow paths between deep brine aquifers and shallower freshwater aquifers referred to in the U.S. as underground sources of drinking water (USDWs) (IPCC, 2005;Nicot, 2009;Nordbotten et al., 2004;Young, 1992). Leaky wells are therefore potential conduits for fluid leakage and groundwater contamination related to deep injection projects that must be identified, evaluated, and potentially remediated prior to permitting of any fluid injection operations (Avci, 1994;Gass et al., 1977;Javandel et al., 1988;Lesage et al., 1991). In the case of an aquifer with well penetrations targeted for fluid injection that is overpressured relative to overlying freshwater aquifer(s), there is potential for well leakage before any injection occurs when any existing wells suffer from well-sealing component failure(s). We refer to the potential driving force created by natural hydraulic gradients before injection as pre-injection relative overpressure.
The aim of this work is to determine how leakage through leaky wells varies as a function of aquitard-aquifer-well properties and natural pressurization rates in abnormally pressured multilayered systems. We use a semi-analytical modeling approach for this purpose. Analytical solutions, when available, can be computationally very efficient for problems containing many injection and leaky wells in multilayered systems, compared to numerical simulations requiring local mesh refinement around each well to obtain accurate results (Azzolina et al., 2013;Birkholzer et al., 2012;Celia et al., 2011;Cihan et al., 2013). The literature includes a variety of analytical solutions for leakage through multilayered aquifer systems (e.g., Cihan et al., 2011;Das & Hassanzadeh, 2021;Dejam & Hassanzadeh, 2018;Zeidouni, 2012). However, development of these solutions commonly includes the assumption of initially hydrostatic pressures in the entire subsurface system. Cihan et al. (2011) presented semi-analytical solutions for initially hydrostatic multilayered aquifer systems and coded them into a FORTRAN program called ASLMA (Analytical Solution for Leakage in Multilayered Aquifers). Modifications to ASLMA were made to generate the results presented by Oldenburg et al. (2016) to focus on leakage-only solutions (i.e., well leakage through an impermeable aquitard) for initially overpressured multilayered aquifers. The results were used to assess incremental brine leakage into a USDW for a CO 2 injection scenario compared to a no-injection scenario. Similarly, Burton-Kelly et al. (2021) recently employed this same modified ASLMA code to assess brine leakage risk in a CO 2 storage project in North Dakota. These example applications in overpressured subsurface systems assume impermeable aquitards and no ongoing pressurization source. Aquitards are known to have non-zero permeability that can vary over a wide range from 10 −15 to 10 −23 m 2 , as demonstrated through core, well, and regional-scale studies (Deming, 1993;Neuzil, 1994;Yang & Aplin, 2007, 2010. Regional-scale permeabilities can be orders of magnitude higher than core-scale permeabilities (Hart et al., 2006;Neuzil, 1986). With non-zero permeability, aquitards can dissipate injection-induced pressures over large areal extents and thereby may play a significant role for reducing leakage rates through leaky wells Cihan et al., 2011). As demonstrated by Neuzil (1995), overpressures cannot persist for geologic times without extremely low-permeability aquitards and/or an ongoing pressurizing process. Neuzil (1995) estimated the geologic forcing rates from the different natural pressurization processes listed above to range between 10 −15 s −1 and 10 −13 s −1 . The latter value may represent an upper limit of the over pressurization rate. In the presence of geologic forcing at these estimated rates, overpressures may persist even if the aquitards have relatively high permeability values (e.g., 10 −17 m 2 ). Consequently, although the injection-induced pressures can dissipate rapidly through aquitards, coexisting geologic forcing may add to the manmade forcing to create sustained rates of brine leakage into freshwater aquifers through leaky wells.
This study presents an extension of the semi-analytical solutions for non-hydrostatic multilayered subsurface systems coded in ASLMA to include hydraulic storage and Darcy flow in aquitards, over-or under-pressurization rates in aquifers and aquitards, and multiple injection, extraction and leaky wells. The ASLMA code with these enhancements is now called SALSA (Semi-Analytical Leakage Solutions for Aquifers). Following the problem description, Section 3 presents the development of the solution using a methodology based on Laplace Transform and Matrix Calculus. Section 4 shows applications of the semi-analytical model using a similar problem to that of Oldenburg et al. (2016)

Problem Formulation
The problem is formulated by considering a similar multilayered system and notation as given in Cihan et al. (2011). As depicted in Figure 1 adapted from Cihan et al. (2011), the multilayered system consists of a sequence of aquifers and aquitards of infinite lateral extent. Aquifers are numbered from the bottom i = 1 to the top i = N. The specific example in Figure 1 has N + 1 aquitards (globally numbered from the bottom i = 0 to the top i = N) and includes aquitards at both the bottom and the top. However, other configurations are also possible, that is, the number of the aquitards can range from N − 1 to N + 1, depending on whether the top or the bottom of the multilayered system is an aquitard or an aquifer. We assume horizontal flow in aquifers and vertical flow through aquitards. These assumptions can be justified as long as the ratio of hydraulic conductivity between the aquifers and the aquitards is larger than 100 (Neuman & Witherspoon, 1969). Each aquifer i is homogeneous and isotropic with horizontal hydraulic conductivity (L T −1 ), specific storativity (L −1 ) and uniform thickness (L). Each aquitard i is also homogeneous and isotropic with vertical hydraulic conductivity and uniform specific storativity and thickness. We define a local aquitard numbering and coordinate system relative to each aquifer i with the overlying ( = +) or underlying aquitards ( = −) (Figure 1b). The aquitard properties are denoted by vertical hydraulic conductivity , specific storativity and thickness .
The multilayered system considered can involve any number of injection (or extraction) and leaky wells that may penetrate multiple aquifers. As illustrated in Figure 1c, the increased hydraulic head in aquifer i can drive fluid into the wellbore of the leaky well l (l = 1,…,N L where N L is the number of leaky wells) with a radius of , and the fluid is then transported via the wellbore into the overlying or the underlying aquifer(s). The horizontal flow rate between aquifer i and leaky well l is denoted by (L 3 T −1 ). The leakage from aquifer i through well l is the sum of the rates of leakage through the overlying ( + ) and underlying ( − ) well-aquitard segments, assuming there is no storage effect in the wellbore segment. If the well-aquifer segment is impervious, for example, a cased wellbore, = 0 .

of 18
Although the governing equations and the solutions in this manuscript are presented for a horizontally isotropic system, the solutions can be applied for horizontally anisotropic aquifers by changing coordinate, well radii, and hydraulic conductivity variables as shown in Cihan et al. (2014).

Governing Equations
The governing equations for single-phase radial Darcy flow in aquifers and one-dimensional vertical Darcy flow in aquitards can be expressed in terms of hydraulic heads as (Bear, 1972;Cheng & Morohunfola, 1993;Hemker, 1985;Hunt, 1985;Maas, 1987;Moench, 1985;Zhou et al., 2009): where ( ) and = is the dimensionless local vertical coordinate, where = 0 at the interface between aquifer i and aquitard (i, α) and = 1 at the interface between aquifer + and aquitard (i, α) (see Figure 1b). Note that i + α = i + 1 for α = +, while + = − 1 for α = −. Following Neuzil (1995), we lump the different complex geologic forcing phenomena into pressurization source terms that can exist in individual aquifer and aquitard layers. Γ and Γ (T −1 ) represent the geologic forcing rates in aquifer i and aquitard (i, α). We also assume that the hydraulic properties of aquifers and aquitards remain the same regardless of the pressurizations.
The radial flow equations in aquifers are coupled through the vertical flows in aquitards that are represented by in Equation 1. (L T −1 ), the flow rate per unit area of the aquifer-aquitard interface from aquifer i into overlying ( = +) or underlying ( = −) aquitard, can be expressed as: (3)

Initial and Boundary Conditions
Geologic forcing processes that may exist throughout geologic time in the multilayered system cause initial head differences in the aquifers and the aquitards. Here we assume that t = 0 is the time when well leakage starts to occur through leaky well(s). While the initial heads in each aquifer are uniform over its thickness, the initial heads in the aquitards are assumed to change as functions of the distance from the aquifer-aquitard interface ( ), pressurization rates in the aquitards at t < 0, and the initial heads in the overlying and the underlying aquifers (ℎ 0) . The initial conditions for the heads can be written as follows: Equation 5 is a steady state solution of Equation 2 for head in aquitard ( + ) . Γ 0 represents an average constant pressurization rate that existed in the aquitard at t < 0. Initial head profiles other than the quadratic profile in Equation 5 may be used to solve Equation 2; and in that case the coupling variables f and g and the nonhomogeneous term R defined in Section 3 need to be modified. Irrespective of the chosen initial head profiles in the aquitards, the head must be continuous at aquifer-aquitard interfaces, which is expressed by: CIHAN ET AL. 10.1029/2022WR032343

of 18
We employ the following boundary conditions for the lateral boundary, with the assumption that the perturbed region enclosing wells (injection, extraction or leaky wells) is very far from the actual lateral boundary of the domain: For a multilayered system containing aquitards at the lowermost or uppermost layers, the bottom and the top boundaries of the system may have either prescribed head (fixed or time-dependent) or no-flow conditions: For the no-flow boundary conditions specified at the bottom or the top boundaries of the aquitards, the initial condition for aquitard head given in Equation 5 is replaced with ) .

Boundary Condition at Leaky Wells
The boundary condition at leaky wells is formulated based on Cihan et al. (2011). The rate of fluid exchange between aquifer i and a leaky well (l) depends on the vertical hydraulic head gradient and the flow resistance within the well. For a cylindrical well geometry, the boundary condition at the radial wall of leaky well l in aquifer i is given by where is the hydraulic head due to the unknown well leakage rate ( ) at aquifer i-leaky well l segment. The continuity equation for the net leakage rate into or out of aquifer i through a leaky well l connecting three consecutive aquifers ( − 1, , + 1) (see Figure 1c) is given by where h is the total head (transient actual head based on the superposition) evaluated at leaky well l-aquifer i segment under the effects of all injection, extraction, and leaky wells. Ω (TL −2 ) is the resistance to flow through the well-aquitard segment expressed as Ω = ∕ . Readers are referred to Cihan et al. (2011), (2013) for a more detailed description of the boundary conditions including leaky wells with cased well-aquifer and plugged well-aquitard segments.

Boundary Condition at Injection or Extraction Wells
The boundary condition at the radial wall of a well m screened in aquifer i that injects or extracts fluid with known constant or time-dependent rate is given by (Cihan et al., 2011): where ( ) is the rate (+ for injection and -for extraction) through well m into aquifer i (Figure 1a), r wi,m is the radius of the well at the well m-aquifer i segment, and is the head in aquifer i due to injection at well m. For the wells screened at multiple aquifers, individual injection or extraction rates can differ at each screened well-aquifer segment. If the entire reservoir system is initially at hydrostatic condition, because flow resistances along well-aquifer and well-aquitard segments are very low compared to those of the aquifers, the individual rates may approximately be calculated using , where Q m is the total rate and N s is the total number of screened aquifer layers. However, this approximation cannot be used for initially non-hydrostatic 6 of 18 multilayered systems. For the initially non-hydrostatic systems, the individual rates are unknown and must be calculated following a similar approach as for leaky wells in Section 2.2.1. The boundary condition at the radial wall of a well screened at multiple aquifers can be expressed as where ( ) is the unknown flow rate at well m -aquifer i segment. If the aquifer i is the first screened zone from the top, ( ) can be expressed as follows: where is the total head evaluated at well m-aquifer i segment under the effects of all injection, extraction and leaky wells. The flow rates at the rest of the screened well-aquifer segments are expressed using the continuity equations as in Equation 10 for any injection or extraction well.

Solution Development by Use of Laplace Transform and Matrix Calculus
The Laplace transform of the flow equations (Equations 1 and 2), subject to the initial conditions in Equations 4 and 5, leads to where p is the Laplace transform variable, and the overbar denotes the Laplace transform of a variable. By employing the Laplace transform of the boundary conditions in Equation 6, the solution of Equation 15 for the head in an aquitard (i, α) can be expressed as where ( ) = √ ∕ . The solution given above can also be applied if the lowermost or the uppermost layers are aquitards with the prescribed head boundary conditions at the bottom and top of the multilayered system (see Equation 8), by changing ℎ + in Equation 16 to ℎ (for = 1 and = − ) or ℎ (for = and = + ). Substituting Equation 16 into the Laplace transform of Equation 3, the flow rate per unit area of the aquifer-aquitard interface in Laplace domain, is expressed as where CIHAN ET AL.
10.1029/2022WR032343 7 of 18 For a multilayered system with a no-flow boundary condition specified at the bottom or the top boundary (Equation 8), the head(s) in the corresponding aquitard (s) and at aquitard-aquifer interfaces are expressed as By substituting Equation 17 into Equation 14, we obtain the following flow equations for aquifers: Different from the equations solved in Cihan et al. (2011), Equation 20 forms a coupled system of non-homogeneous ordinary differential equations (ODE). The non-homogeneous term is expressed as Note that the non-homogeneous terms R 1 and R N for Aquifer 1 and Aquifer N must be modified by using Equation 19 for aquifers underlaid or overlaid with an aquitard that has a no-flow boundary condition specified at the bottom or the top boundary.
Equation 20 can be expressed in matrix notation: where matrix A, [a i, j ] N×N , referred to as the diffuse-leakage-coupling matrix (Cihan et al., 2011), is a tri-diagonal matrix, with its components: The system of ODEs, Equation 22, can be decoupled by finding the eigenvalues λ and the eigenvectors ξ for the following eigenvalue system (Cheng & Morohunfola, 1993;Cihan et al., 2011;Hemker & Maas, 1987;Hunt, 1985): where I is a unit diagonal matrix. The matrix A is in general a non-symmetric N × N matrix because the transmissivities (T i ) of each aquifer may be different. Following Maas (1986) and as in Cihan et al. (2011), matrix A is transformed into a symmetrical matrix, Aʹ , using A' = T 1/2 A T −1/2 and ξ' = T 1/2 ξ, where T is the diagonal transmissivity matrix with components T ii = T i . The eigenvalue analysis is done for Aʹ to obtain N real positive eigenvalues λ i with N real independent normalized eigenvectors ′ . Then, the eigenvectors ξ are computed by ξ = −1∕2 ′ , which are orthonormal relative to T (ξ −1 = ξ T T) as shown by Hunt (1985).
Multiplying both sides of Equation 25 by ξ −1 and setting * = −1 , we obtain: The solution of Equation 26 subject to the boundary conditions in Equation 7 can be expressed in matrix notation as * where K is N × N diagonal matrix with the nonzero elements K0 ( √ ) for i = 1,…,N. K 0 is the zeroth-order modified Bessel function of the second kind. c is a column matrix that includes the coefficients obtained from the boundary conditions (e.g., Equations 9 and 11 or 12) at the wall of a single well as a function of the Laplace variable p. By transformation = * , a solution for is given by In the summation form, the hydraulic head in Laplace domain for each aquifer i is expressed as Cihan et al. (2011) showed the calculation methods of c j coefficients with detailed examples both for known rates of injection or extraction at well-aquifer segments and for unknown leakage rates through leaky wells. The readers are referred to Section 3.3 in Cihan et al. (2011) for the calculation details, which are not repeated here.
For the unknown leakage rates for leaky wells or unknown flow rates for multi-segment injection (or extraction) wells, a system of algebraic equations (formed based on Equations 9, 10, 12 and 13) is solved for calculating the unknown N × (N L + N I ) coefficients including for N L leaky wells (l = 1,..,N L ) and for N I injection (or extraction) wells (m = 1,..,N I ).
A general solution for heads in aquifers, based on the superposition principle, for an abnormally pressured multilayered system containing N aquifers, N I injection and N L leaky wells can be expressed as where is the distance between an observation point at a global horizontal coordinate (x, y) and injection (or extraction) well m located at ( , ) , and is the distance between the observation point and leaky well l at ( , ) , which are calculated by The solutions presented here (i.e., for heads in aquifers and aquitards and leakage rates from and into aquifers) in Laplace domain are numerically inverted into the real time domain using the Stehfest numerical Laplace Inversion method (Stehfest, 1970a(Stehfest, , 1970b, as described in the Appendix of Cihan et al. (2011). For the application in the next section, we have incorporated the extended solution into a new version of the ASLMA code we call SALSA (Semi-Analytical Leakage Solutions for Aquifers) that includes computations of the eigenvalues and eigenvectors, and coefficients and the numerical Laplace inversion.

Application for Assessing Brine Leakage Risks in Overpressured CO 2 Storage Reservoirs
In this work, we revisit the problem in Oldenburg et al. (2016) that mimicked a research CO 2 storage project at a mid-continent storage site in the US. The problem here involves a two-aquifer-three-aquitard system (Figure 2), instead of the two-aquifer-one-aquitard system in Oldenburg et al. (2016). The storage aquifer (Aquifer 1) at the bottom is overpressured relative to the top aquifer (Aquifer 2) identified as USDW. We assume that the source of the overpressure is the lower aquitard, from which higher pressures generated migrate to the overlying aquifers and aquitards. Our focus in this study is on brine leakage at distances beyond the CO 2 plume region. Earlier works have showed that both head changes and brine flow rates outside of the CO 2 plume can be described with reasonable accuracy by single-phase flow calculations (e.g., Cihan et al., 2013;Nicot, 2008). We analyze the impacts of overpressure on brine leakage for potential geologic forcing rates of 10 −15 s −1 to 10 −13 s −1 , based on Neuzil (1995). No-flow boundary conditions are used at the bottom of the domain. We also assume that hydrostatic pressure conditions prevail due to fresh water recharge at depths above the upper aquitard, which are represented by using a fixed head boundary condition at the top boundary of the domain. The system includes a high-conductivity leaky well that penetrates the entire domain and starts leaking brine from Aquifer 1 to Aquifer 2 at t = 0. According to the scenario in Oldenburg et al. (2016), the leakage through the well occurs due to pre-injection relative overpressure for an arbitrary pre-injection period of 50 years with a decreasing rate. Then, the leakage rate increases due to pressurization caused by the start of CO 2 injection at t = 50 years through an injection well near the existing leaky well. The injection occurs at a rate of 835.32 m 3 /d (represented as equivalent-volume injection of brine) over 4 years. The radius of both wells is 0.15 m. It is important to explore how additional pressure increases induced by injection within an already overpressurized zone may cause increase in fluid leakage through the hypothetical leaky well. This is relevant to the EPA guidance for CO 2 injection permit applications in overpressured aquifers (Oldenburg et al., 2016;US EPA, 2013). 10 of 18 Each of the aquifers is 50 m thick with a permeability of 3 × 10 −14 m 2 . However, there are slight differences in hydraulic conductivity and specific storativity of the aquifers due to changes in salinity and temperature ( Table  1 in Oldenburg et al., 2016). Aquifer 1 has a hydraulic conductivity of 2.98 × 10 −2 m/d and a specific storativity of 2.11 × 10 −6 m −1 . Aquifer 2 has a hydraulic conductivity of 2.75 × 10 −2 m/d and a specific storativity of 2.04 × 10 −6 m −1 . Each aquitard layer has a thickness of 186 m and a specific storativity of 5.32 × 10 −6 m −1 . In contrast to Oldenburg et al. (2016), here we consider that the aquitards have a finite non-zero permeability (10 −19 m 2 ≤ k aqt ≤ 10 −17 m 2 or 1.03 × 10 −7 m/d ≤ K′ ≤ 1.03 × 10 −5 m/d) and investigate how the geologic forcing and pressure propagation through aquitards can affect the well leakage rates before and after the injection. For this finite range of aquitard permeability values, the capillary entry pressures are expected to be high enough for preventing CO 2 phase flow through aquitards.

Verification of the Semi-Analytical Solution
We use a numerical model to verify the semi-analytical solution presented in Sections 2 and 3 for non-hydrostatic multilayered systems. The numerical model, described in detail in Supporting Information, solves the three-coupled nonlinear partial differential equations for pressure, salt mass fraction, and temperature using the Finite Volume method. Unlike the semi-analytical model, the numerical model takes into account density and viscosity variations as a function of pressure, temperature, and salinity. For the verification, we considered the leakage through a leaky well for the two-aquifer-three-aquitard problem with no injection, as shown in Figure 2.
We solved the numerical model in cylindrical coordinates for the single leaky well (see Supporting Information), assuming Darcy flow through the well and constant rate of pressurization in the lower aquitard. We compared the predictions of the analytical model for brine leakage into the top aquifer (USDW) with those of the numerical model. Figure 3 shows that the leakage rate predictions of the models match very well for different pressurization rates that might exist in natural systems based on Neuzil (1995). Even though the salinity migration and attenuation in the aquifer cannot be modeled using the semi-analytical solution, this comparison shows the applicability of the semi-analytical model to accurately estimate leakage rates. As seen in Figure 3, the leakage rates decrease with time, and the rates change significantly for different pressurization rates. Decreasing leakage rates with time occur because of decreasing potential gradient between the aquifers. Additional model verifications for different aquitard permeability values are presented in Supporting Information S1.

Effects of Geologic Forcing Rates and Pressure Dissipation Through Aquitards on Well Leakage
The modeled scenario in this section includes a 50-year pre-injection leakage period, a 4-year injection period and a long term post-injection period. The injection well is located at 2 km away from the leaky well. The initial head distribution at the onset of well leakage (at t = 0) is used from the steady-state solutions (Γ′ = Γ0′ for Αquitard 1) without the leaky well. Figure 4 shows examples of calculated head profiles that have evolved from initially hydrostatic conditions (h-h top = 0) to steady state over geologic times due to pressurization in the lower aquitard (Γ′ = 10 −14 s −1 ) . The time to reach the steady state decreases with increasing aquitard permeability. For the same pressurization rate of Γ′ = 10 −14 s −1 , the calculated steady state heads as a function of depth take the smallest values for = 10 −17 m 2 (Figure 4a). The steady state head increases (h-h top ) are directly proportional to pressurization rate and inversely proportional to aquitard permeability. As a result, the same steady-state profiles are obtained when the aquitard permeability and pressurization rate are changed at the same rate. For example, the steady state head profile that is computed for = 10 −18 m 2 and Γ′ = 10 −14 s −1 (Figure 4b) is identical to the steady state head profile for a system with = 10 −17 m 2 and Γ′ = 10 −13 s −1 (not shown here). Figure 5 shows that the leakage rates increase proportionally with increasing pressurization rates and decreasing aquitard permeability values. The leakage rates rapidly decrease after the initiation of well leakage because the flow through the high-conductivity leaky well ( = 10 −7 m 2 ) used in this scenario quickly reduces the head gradient between the aquifers. The head gradient increases again during the injection period, resulting in spikes in leakage rates between 50 and 54 years. Due to continuous pressurization, the leakage rates eventually reach steady-state values. Comparing the plots in Figure 5, it is seen that the leakage rates for = 10 −17 m 2 and Γ′ = 10 −13 s −1 (Figure 5a) are similar to the leakage rates for = 10 −18 m 2 and Γ′ = 10 −14 s −1 (Figure 5b) and for = 10 −19 m 2 and Γ′ = 10 −15 s −1 (Figure 5c).
CIHAN ET AL.
10.1029/2022WR032343 11 of 18  12 of 18 Figure 6 presents variation of leakage rates with (Γ′ ≠ 0) and without (Γ′ = 0) ongoing source of pressurization. Γ′ = 0 denotes a scenario in which the geologic forcing process might have ended in very recent geologic times. In Figure 6, the initial head distributions are the same for all the cases when the well leakage starts at t = 0. Without an ongoing source of pressurization, the well leakage rates go to zero more rapidly for = 10 −17 m 2 because the initial overpressures quickly dissipate for such high aquitard permeabilities. As shown by Neuzil (1995), this quick pressure dissipation also implies that overpressures in natural systems may not persist for geological times without either extremely low-permeability aquitards and/or continuous overpressurization forcings, such as, for example, by sediment deposition, tectonic compression, hydrocarbon generation, etc. The effective leaky well permeability of 10 −7 m 2 that is used in the previous figures and Figure 6a may represent a case of Pressurization takes place only in the lower Aquitard. Initial head distribution at the onset of well leakage (at t = 0) are used from the steady-state solutions (Γ′ = Γ0′ for Αquitard 1) without a leaky well. Well leakage starts at t = 0, and injection occurs between t = 50 years and t = 54 years. The distance between the injection well and the leaky well is 2 km. Leaky well permeability = 10 −7 m 2 .

of 18
leakage through an open-well bore with minimal flow resistance. However, the effective permeability can be much smaller for wells with leaky paths through the annulus or the damaged/degraded cement zone outside of the casing (Burton-Kelly et al., 2021). Figure 6b shows that as the leaky well permeability decreases from 10 −7 m 2 to 10 −10 m 2 , the leakage rates decrease more than tenfold. Figure 7 shows head profiles calculated for an observation well located 10 m away from the leaky well and 1,990 m away from the injection well. As seen in Figure 7, the initial head differences observed in the aquifers decrease over time before and after the injection as the flow through the leaky well occurs from the lower aquifer to the upper aquifer. For Γ′ = 0 , the leakage continues with a decreasing rate in the post-injection period, until the system equilibrates relative to the hydrostatic head set at the top boundary of the domain (Figure 7a). Figure 7b shows that for Γ′10 −14 s −1 = 10 −14 s −1 , the system reaches a new equilibrium state for head distribution about one thousand years after the onset of well leakage.

Assessment of Leakage Enhancement Due To Injection
Figure 8 presents change of the incremental cumulative well leakage at 100 years due to the 4-year injection as a function of the distance between the injection well and the leaky well. The incremental leakage is defined as the cumulative leakage minus the baseline cumulative leakage with no injection. Any injection-induced head increases in the storage aquifer decrease with distance from the injection well as shown in Figure 9. Because the well leakage is driven by the head difference between the aquifers, the incremental cumulative leakage also decreases as the distance to the leaky well increases. Figures 9a and 9b show that the head increases around the leaky well located at 2 km distance from the injection well are much greater than the head increases around the leaky well at 6 km distance. The leakage through the well at 2 km causes slightly smaller head increases at the injection well screened at Aquifer 1 and creates greater head changes in Aquifer 2 (USDW).
In GCS, a very low-permeability aquitard with high capillary entry pressure to free-phase CO 2 is desired to serve as the caprock to prevent leakage of CO 2 from the storage reservoir. Nevertheless, Figure 8 shows that the brine leakage impact to USDW through leaky wells might substantially decrease with increasing aquitard permeability values. As demonstrated by , pressure dissipation through aquitards may be significant for aquitard permeabilities as low as 10 −18 m 2 . Consequently, this "pressure bleed-off" in the storage aquifer helps reduce brine well leakage, as long as the aquitard permeability is low enough (and for free-phase CO 2 , the capillary entry pressure is high enough) to prevent CO 2 leakage. Figure 8 also shows the incremental cumulative well leakage drops more rapidly with the distance from the injection well when the aquitard permeability is higher. Figure 8b presents the incremental cumulative leakages for a smaller value of leaky well permeability (10 −10 m 2 ). The incremental cumulative leakage volumes decrease more than about an order of magnitude with a decrease in leaky well permeability from 10 −7 m 2 to 10 −10 m 2 . For the 10 −10 m 2 leaky well permeability and Γ′ = 10 −14 s −1 , the incremental leakage for leaky wells at distances greater than 5 km is less than 1,000 m 3 for > 10 −17 m 2 . In practice, the calculation of incremental cumulative leakage changes as presented here for overpressured storage systems may help in the delineation of Area of Review for site screening of potential leaky paths, before CO 2 injection in the case of pre-injection relative overpressure (Oldenburg et al., 2016). In addition to calculations that include dissipation by the aquitards, other natural attenuation processes in USDWs that can dilute brine concentrations should be considered to determine the degree to which cumulative leakage from an overpressured system may be acceptable, and simulating these processes requires use of numerical models (e.g., Figure 2).

Conclusions
Abnormal fluid pressures (above or below hydrostatic) are commonly observed in sedimentary basins consisting of aquifers alternating with aquitards. In the context of GCS, while caprocks or aquitards are relied upon to contain injected buoyant free-phase CO 2 and limit inter-aquifer brine flow, the large number of wells that could be present in some sedimentary basins creates concerns for potential leakage and groundwater contamination. This is because leaky wells may provide flow paths between deep brine aquifers and shallower freshwater aquifers that can sustain fluid leakage. Simply put, pre-injection relative overpressures in the deep aquifers may cause fluid leakage into the freshwater aquifers through leaky wells even before any injection occurs.
This study presents development of a semi-analytical solution for head changes and leakage in abnormally pressured multilayered aquifer-aquitard systems with geologic forcing sources and any number of injection, extraction and leaky wells. The presented semi-analytical solution, built into a computer code named SALSA (Semi-Analytical Leakage Solutions for Aquifers), was applied to a problem that mimics a research CO 2 storage project at a low-permeability mid-continent US site. The problem involves a two-aquifer-three-aquitard system with a storage aquifer overpressured relative to a shallower freshwater aquifer, due to pressurization originating in the deep subsurface. Leakage through a leaky well, located near an injection well, is assumed to start in a pre-injection time due to the relative overpressure and continue during injection and post-injection periods.
The model results show that without an ongoing source of pressurization, the initial overpressures quickly dissipate for high aquitard permeabilities. Our results on the quick pressure dissipation confirm the earlier findings (e.g., in Neuzil, 1995) that overpressures in natural systems may not persist for geological times without either extremely low-permeability aquitards and/or an active overpressurization source. The pre and post-injection model results show that pressure dissipation over large areal extents provided by non-zero permeability in overand/or underlying aquitards can potentially reduce leakage rates through leaky wells. However, presence of ongoing pressurization due to some geologic forcing process may create sustained rates of brine leakage into freshwater aquifers through leaky wells. For a fixed rate of pressurization, any injection-induced incremental cumulative  (b) with ongoing pressurization at 0 (Γ0′ = Γ′10 −14 s −1 ) . Initial head distributions are the same at the start of well leakage at t = 0. The distance between the injection well and the leaky well is 2 km. Leaky well permeability = 10 −7 m 2 and aquitard permeability = 10 −18 m 2 .
15 of 18 leakage, defined as the cumulative leakage minus the baseline cumulative leakage with no injection, decreases with increasing leaky well distance from the injection well and increasing aquitard permeability. In GCS, a very low-permeability aquitard with high capillary entry pressure to free-phase CO 2 serving as a caprock is desired for preventing leakage of CO 2 through the caprock. Nevertheless, the model results show that the effect of brine leakage caused by injection into shallow freshwater aquifers via leaky wells might substantially decrease with increasing aquitard permeability values, as long as the aquitard permeability is low enough (and for free-phase CO 2 , the capillary entry pressure is high enough) to prevent CO 2 leakage.