A Simultaneous Model for Simulating the Propagation of Hydraulic Head in Aquifers with Leakage Pathways

A three-dimensional simultaneous solution model was developed for analysis of variable transient flows, specifically leakage associated with arbitrary groundwater aquifers. The model simultaneously calculates leakage rates using the hydraulic gradients between coupled leakage points at two leaky aquifers and evaluates the propagation of hydraulic heads in multiple aquifers resulting from that leakage. Important considerations for leakage simulation are variant leakage rates with time and leakage starting time. Two types of leakage pathways are specified in terms of the generated time, including (1) pre-existing leakage pathways and (2) abruptly-induced leakage pathways at specific times. The governing equation with a leakage term is composed of three finite difference equations, and its form depends on whether leakage pathways pre-exist or are induced at specified times according to known or assumed aquifer histories. The developed numerical code was validated by comparison with the results of the TOUGH2/EOS1 program for a three-dimensional conceptual domain.


Introduction
We developed a simultaneous solution model for the three-dimensional analysis of transient flow in arbitrary multilayered groundwater aquifers with leakage.The developed model uses the Finite Difference Method (FDM) and is designed to understand the propagation of groundwater through several possible leakage zones in porous media, such as fractures and abandoned wells.We view this research as an applied single-phase flow analysis in a research program focused on storage of CO 2 in geological formations and associated potential leakage.In the geologic carbon storage (GCS) project, CO 2 injected into deep storage reservoirs saturated with brine may induce pressure buildup and migration of brine.The existence of fractures and faults or abandoned pre-existing wells in the storage system can cause the leakage of pressurized brine into adjacent formations.Before CO 2 leakage occurs, the adjacent aquifers saturated with brine act as single-phase reservoirs; therefore, the single-phase code can be an effective tool in evaluating the brine leakage features in multiphase flow systems (Lee et al., in press).In addition, Cihan et al. (2013) demonstrated that the analytical solution model for single-phase fluids is efficient in evaluating the regional pressure buildup in a two-phase flow system.This study was thus conducted as part of fundamental research to realize pressure anomalies induced by CO 2 leakage in the GCS.
Researchers have recently developed the methodologies to quantify leakage behaviors and pressure anomalies induced by leakage in multilayered aquifer systems: Nordbotten et al. (2004) developed an analytical methodology to simulate the leakage in multiple aquifers.Zhou et al. (2009) applied both one-dimensional radial flow equations and vertical flow equations for the semi-analytical leakage modeling.Cihan et al. (2011) studied semi-analytical solutions for pressure anomalies induced by leakage in a multilayered aquifer system.This method was later applied to simulating pressure perturbations resulting from leakage pathways of various sizes and permeability values by Cihan et al. (2013).Veling and Maas (2009) also studied the method of solution for both horizontal and vertical flow induced by leakage in multilayered aquifers.The commercial MLU (Multi-Layer Unsteady State) code was developed for leakage diffuse due to pumping and injection in multiple aquifers (Hemker & Post, 2013).Those existing studies focused on the methodology to model leakage resulting from pre-existing leakage pathways but not abruptly-induced leakage pathways which can be induced by an external effect such as a micro-seismic event caused by overpressure from water (or CO 2 ).Therefore, in this study, we describe abruptly induced leakage pathways at specific times as well as pre-existing leakage pathways for more realistic predictions of leakage behavior and pressure perturbation in nature leaky multiple aquifers.This study provides an improvement in methods discussed above to solve the flows in porous media.
The developed code simulates the transient release of leakage into adjacent aquifers along abandoned wells and fractures penetrating confining layers from the confined aquifers with higher hydraulic heads (or pressures) by water injection and anomalies of hydraulic heads in the aquifers caused by the transient release of leakage.The leakage simulations can be performed using a three-dimensional flow equation with a leakage term intended to represent leakage.Leakage rates expressed as Darcy's flow equation depend on the hydraulic conductivity of the leakage pathway, the cross-sectional area of the leakage pathway, and the hydraulic gradient between the leakage aquifers overlying and underlying the confining bed.The hydraulic gradient between the leakage aquifers is calculated from two leakage points on the leakage pathway.Thus, the inflow and outflow rates along the leakage pathway through the confining bed are calculated based on the hydraulic conductivity and the cross-sectional area of the leakage pathway, and the hydraulic heads of the confined leakage aquifers caused by leakages can be modeled simultaneously.Section 2 introduces the basic theory and the methodology of the leakage modeling.In Section 3, a synthetic model domain and parameters are described.In addition, Section 3 illustrates the simulation results from different leakage applications to both the pre-existing and the abruptly-induced leakage pathways in the model domain, and the developed model is verified by comparison with the results of a commercial numerical code for pre-existing leakage scenarios.

Methodology
As suggested previously, this study is associated with potential CO 2 leakage in the storage of CO 2 in geological formations.Cihan et al. (2013) suggested that the methodology of solving for the groundwater (or single-phase fluid) governing equation can be applied to simulations of mobile brine/CO 2 entering or exiting the confined aquifers through leakage pathways if equivalent volumes of water are used to mimic the CO 2 injection in the CO 2 storage reservoirs.Lee et al. (in press) also determined that the inverse simulator for a single-phase fluid is applicable for leakage pathway estimation in the multilayered system of two-phase fluids (brine and CO 2 ).However, those studies were applied to the analysis of pressure induced by a pre-existing leakage pathway.Increased pressure from water (or CO 2 ) injection can cause fractures or cracks in confining beds with vulnerable areas (for example, incompletely plugged abandoned wells).In addition, leakage pathways can be abruptly induced.In this study, the important consideration is the abruptly-induced leakage at specific times.However, this study does not consider specific mechanisms for how leakage pathways are generated in the confining beds by increased pressure.Rather, these analyses focus only on the migration of leakage along leakage pathways under the assumption that the leakage pathways are generated in the confining beds.Two types of leakage pathways are specified in terms of the generated time: (1) pre-existing leakage pathways and (2) induced leakage pathways at specific times.First, the pre-existing leakage pathway corresponds to a leakage pathway that is completely generated and saturated before water injection.The instantaneous inflow rate into the leakage pathway is assumed to be equivalent to the outflow rate from the leakage pathway, corresponding to the fundamental continuity equation.Second, the induced leakage pathway corresponds to a leakage pathway that is generated at an arbitrary time through injection-induced increased pressure or other mechanisms.The model releases leakage from the pathway into an adjacent confined aquifer at the time when the pathway is generated.By assumption (see the next section), the inflow rate into the leakage pathway and the outflow rate from the leakage pathway are simultaneous and equivalent.

Governing Equation and Methodology of Leakage Simulation
The governing equation is the mass conservation equation for the three-dimensional movement of groundwater of constant density through porous media.The governing equation may be described by a partial differential equation as follows: where , , and are values of hydraulic conductivity (L/T) along each direction; ℎ is the hydraulic head (L); ); is sources and/or sinks of water (a volumetric flux per unit volume) (T -1 ); is the leakage term, the leakage rate per unit volume (T -1 ) ; Ss is the specific storage coefficient of the porous media (L -1 ); and t is time (T).
Equation ( 1) describes groundwater flow with a leakage term for a heterogeneous and anisotropic medium.The leakage induces changes of hydraulic heads in the aquifers, which is referred to as "leakage effects" in this study.
The leakage term has been added in the governing equation to consider the leakage effects, and the Finite Difference Method (FDM) is used to solve for the governing equation.In this study, hydraulic heads are calculated at discrete points in space, which are nodes in the centers of meshes; therefore, assigning hydraulic conductivities and specific storage coefficients for the cells is straightforward.In addition, multiple aquifers should be considered to simulate a confined system with possible leakage pathways (Cihan et al., 2011).Figure 1 shows a schematic of the leakage confined system, with a confining bed (or caprock) and two confined aquifers overlying and underlying the confining bed.In this study, the underlying aquifer in which water is injected is referred to as the "storage formation" and the aquifer above the caprock is referred to as the "overlying formation." where the leakage rate ( ) in Equation ( 1) is expressed to , , per unit volume.I is the leakage pathway number, , , is the cross-sectional leakage pathway area, ( ) is the length of I-th leakage pathway, , , is the z-directional hydraulic conductivity of the leakage pathway, ( ) is the z-coordinate at a storage formation of I-th leakage, and ( ) is the z-coordinate at an overlying formation of I-th leakage (Figure 1 depicts how the leakage pathways and parameters are related).It is assumed that there are no lateral flows at the interface between the leakage pathway and the confining bed, i.e., leakage migrates vertically in confining beds along specific leakage pathways.In addition, the leakage pathway is not deformable and the continuity equation is satisfied in the leakage pathway; therefore, the leakage rates in all cross-sections of a leakage pathway are instantaneously identical.
Consistent with Darcy's Law, the leakage rates can change with time because the leakage affects the hydraulic heads in the overlying and storage formations.Thus, the leakage rates should be simultaneously modeled using the hydraulic heads between the overlying and storage aquifers; the transient flow in these two aquifers caused by the leakage must also be calculated.The leakage simulations are performed by solving the hydraulic heads at the "coupled leakage points" within each leakage pathway (see Figure 1).The methodology uses three finite difference equations for leakage simulations: two finite difference equations for coupled leakage points and one finite difference equation for no leakage points (these finite difference equations are discussed in the subsequent section).The hydraulic head at one leakage point is determined using hydraulic heads at seven discrete nodes, which consist of six discrete adjacent FDM nodes and one corresponding leakage point.The hydraulic head at the corresponding leakage point is also determined in the same manner.We can simultaneously calculate hydraulic heads at coupled leakage points and leakage rates, which are based on the difference of hydraulic heads between the coupled leakage points.In general, the leakage simulations using the conventional method are conducted by calculating hydraulic heads (or pressures) at the inner nodes of an explicitly meshed leakage pathway in the confining layers.However, this methodology does not require the geometry of the leakage pathways to be specified by meshes because the leakage term implies the geometry of the leakage pathways, thereby providing an advantage by reducing the leakage simulation computations.

Discretization of the Governing Equation
Figure 2 illustrates x-directional cross-sections through three cells and the numerical approximation of derivatives of hydraulic head under anisotropic and heterogeneous conditions (Bennett, 1976).In terms of node (i, j, k), the notation f indicates the region in which water flows into node (i, j, k) from the upstream node, and the notation b indicates the region in which water flows out from node (i, j, k) to the downstream node.
The effective hydraulic conductivity is calculated as a weighted harmonic mean, as described by Collins (1961).For example, The time derivative of the heads can be approximated with the backward difference method (Burden & Faires, 2005): where n is the time step index.This finite difference formula has a local truncation of order O(∆t).
The x-directional 1 st order partial derivative at node (i, j, k) can be approximated by the arithmetic mean of the fluxes of sections (i) and (ii): , , ( 5 ) The 2 nd order approximation for the x-directional second order partial derivative of head at node (i, j, k) can be expressed as: ( 6 ) These finite difference formulas for 2 nd order partial derivatives have the local truncation error of order O(∆ ).
The finite difference equations of , , and , , can be determined in the same manner.
The three finite difference equations have been specifically derived from the governing equation with the leakage term.One of three finite difference equations is applied to the nodes in a model depending on whether the nodes exhibit the leakage pathways at the specific time.The first finite difference equation is In the modeling scenario, the leakage occurs at three pathways with different properties (see Table 1).The first pathway (LP1), located at node (x, y) = (5300 m, 4100 m), has leakage at time zero (this is a pre-existing pathway).
The second pathway (LP2) at node (5300 m, 5700 m) and third pathway (LP3) at node (5300 m, 5100 m) are induced at time 3.171 years and 4.76 years after the simulation begins, respectively.The three pathways are arranged in a straight line at coordinate (x) = (5300 m).It is assumed that LP2 and LP3 (abruptly-induced leakage pathways) will arise at the confining beds by an external forces such as increased pressure from water injection.Therefore, the simulation using the new model was performed using three leakage conditions with time intervals.

Validation of the Developed Code and Simulation Results
Figure 4 illustrates the simulated pressure distributions at each leak point of the three leakage pathways in the 3 rd and 9 th layers from both the newly developed code (dashed lines) and TOUGH2 code (solid lines).Pressure is a dependent variable of TOUGH2 to describe variable density groundwater flow, and hydraulic head is a dependent variable of the new code to describe constant density groundwater flow; therefore, the results of the new code were converted from the hydraulic head (m) to pressure (Pa) using = ( − ), where P is the pressure, H is the hydraulic head and z is the elevation.
From results of the new code in Figure 4 (a), the leak point of LP1 in the overlying formation reflects gradual changes of hydraulic heads as soon as water is injected into the storage formation caused by the increased gradient of hydraulic heads between coupled leak points of LP1.After the LP2 pathway is generated at 3.171 years, the leak point of LP2 in the 3 rd layer immediately shows a rapid rise of hydraulic heads, and it causes the propagation of hydraulic head anomalies throughout the overlying formation; therefore, the leak points of LP1 and LP3 also exhibit the sudden rise of hydraulic heads.The effect of leakage induced by LP3 at 4.76 years is similar to that of LP2, but the largest increase of hydraulic head occurs from LP3 because of the larger size (0.5 m × 0.5 m) of the LP3 and the substantial hydraulic head gradient between coupled leak points of LP3, which is the closest leakage pathway to the injection well.In Figures 4 (a  Figure 5 illustrates the relative errors between both models at (x) = (5300 m) in the 3 rd and 9 th layers of the YZ-plane, in which LP1 is located at 3.170 years before the abruptly-induced leakage occurs.The maximum relative errors in the 3 rd layer and in the 9 th layer are approximately 4.65E-7 (0.0000465%) at (y) = (4100 m) (i.e., the leak point of LP1) (Figure 5 (a)) and 4.46E-4 (0.0446%) at (y) = (5100 m) (closest to the injection well) (Figure 5 (b)), respectively.The increased errors in the 9 th layer might be caused by the effects of variable water density on changes in pressure.As discussed previously, TOUGH2 calculates pressures based on variable density, while the new code employs hydraulic heads in terms of the constant density; therefore, the change of density caused by an increased pressure buildup in the storage formation can lead to a slight increase in relative errors between both models.
Figure 6 illustrates the hydraulic head propagation in the XY-plane of the third layer over a period of 10 years.Figure 6 (a) represents the hydraulic head distribution at time 3.06 years before LP2 is induced, which slowly evolves at LP1, located at (x, y) = (5300 m, 4100 m).The hydraulic head begins to rapidly increase at LP2 at the coordinate (x, y) = (5300 m, 5700 m) after time 3.171 years when LP2 is generated (Figure 6 (b)).The hydraulic head in the overlying formation is consistently and substantially increased by LP2 (Figure 6 (c)).The occurrence of LP3 causes the sudden increase of hydraulic head in the vicinity of node (5300 m, 5100 m) (Figs. 6 (d) and (e)), and the hydraulic head anomalies are continuously propagated in the 3 rd layer (Figure 6 (f)).As discussed previously, the increasing amount of hydraulic head from LP3 is more substantial because of the larger cross-sectional size of LP3 and the greater hydraulic gradient between the coupled leak points of LP3 than those of other pathways.
The results for the time period prior to 3.170 years demonstrate that the methodology of the newly developed model using the coupled leak points originating from a leakage term in the governing equation can be applied for the transient flow simulations caused by the pre-existing leakage pathway.In addition, the methodology can efficiently recognize hydraulic head anomalies caused by the abruptly induced leakage pathways because it does not require the use of meshes to specify leakage pathway geometries and hydrogeological properties, which is different from the conventional models such as TOUGH2 (i.e., the methodology can readily parameterize/generate the leakage pathway features at a specific time without the meshes).This methodology also provides an advantage in terms of the number of computations because the number of cells decreases in the leakage simulations.Finally, the new model can be used as a tool to evaluate the flow patterns caused by potential leakage in porous media.

Conclusions
This study focused on the analysis of transient flows caused by leakage in porous media.A developed model provides the three-dimensional analysis of transient flows caused by leakage in the arbitrary groundwater aquifers.The leakage effects are considered when the leakage term is added in the governing equation using Darcy's Law.In the model, the leakage rates are simultaneously simulated by using parameterized hydrogeological properties of leakage pathways and hydraulic gradients between multiple aquifers, and the transient flow in two aquifers caused by the leakage migration is calculated.The hydraulic gradient is determined using coupled leakage points and parameterizing the geometries and the hydrogeological properties of leakage pathways.The hydraulic head at one of the coupled leakage points is calculated from aggregate hydraulic heads at six discrete adjacent nodes of the leakage point and one corresponding leakage point.The hydraulic head at the corresponding leakage point is also calculated in the same manner.Thus, the transient/dynamic flow caused by leakage can be effectively simulated without the use of a meshed leakage pathway in the confining layers, so the methodology provides an advantage by reducing the number of computations for the leakage simulations.In particular, unlike the existing studies, the method enables the simulation of the effect of the leakage migration through induced leakage pathways at specific times as well as pre-existing leakage pathways.The new model was validated through application of the TOUGH2/EOS1 model for a three-dimensional conceptual domain.The maximum relative errors in the 3 rd and 9 th layers in the domain between both models were 4.65E-7 and 4.46E-4, respectively.
By the given assumptions, our method does not allow the lateral flow from the leakage pathways into the caprock, so cannot simulate the leakage diffuse through the caprock from the leakage pathway.In the future, our study will solve this problem for a more realistic modeling of the transient flow caused by the leakage in multiple aquifer systems.In addition, this study will be applied to develop a leakage detection system using the inverse analysis.The future study will proceed to be helpful in evaluating/reducing the risk of CO 2 leaks so that it can be applied to the CO 2 storage project in deep subsurface formations with leakage issues.

Figure 1 .
Figure 1.Schematic showing leakage pathways and relation parameters ) and (b), TOUGH2 is simulated from time 0 to 3.170 years (0.99 × 10 8 sec) before LP2 is induced because the TOUGH2 model could not directly simulate pressure anomalies through the abruptly induced leakage pathways.The results from the newly developed code and TOUGH2 for the pre-existing leakage condition (by LP1) through 3.170 years are thus compared to validate the new code.The simulated pressure distributions through 3.170 years from the new code are generally analogous to those from TOUGH2, as shown in Figures4 (a) and (b).

Figure 4 .Figure 6 .
Figure 4. Pressure distributions at each leak point in the (a) 3 rd layer and (b) 9 th layer of three leakage pathways from both models

Table 1 .
Specification of the conceptual model