Leakage dynamics of fault zones: experimental and analytical study with application to CO2 storage

Abstract Fault zones have the potential to act as leakage pathways through low permeability structural seals in geological reservoirs. Faults may facilitate migration of groundwater contaminants and stored anthropogenic carbon dioxide (CO$_2$), where the waste fluids would otherwise remain securely trapped. We present an analytical model that describes the dynamics of leakage through a fault zone cutting multiple aquifers and seals. Current analytical models for a buoyant plume in a semi-infinite porous media are combined with models for a leaking gravity current and a new model motivated by experimental observation, to account for increased pressure gradients within the fault due to an increase in Darcy velocity directly above the fault. In contrast to previous analytical fault models, we verify our results using a series of analogous porous medium tank experiments, with good matching of observed leakage rates and fluid distribution. We demonstrate the utility of the model for the assessment of CO$_2$ storage security, by application to a naturally occurring CO$_2$ reservoir, showing the dependence of the leakage rates and fluid distribution on the fault/aquifer permeability contrast. The framework developed within this study can be used for quick assessment of fluid leakage through fault zones, given a set of input parameters relating to properties of the fault, aquifer and fluids, and can be incorporated into basin-scale models to improve computational efficiency. The results show the utility of using analytical methods and reduced-order modelling in complex geological systems, as well as the value of laboratory porous medium experiments to verify results.


Introduction
Buoyancy-driven flows in porous media have many important practical applications, linked to the movement of fluids in the subsurface. Examples include the characterisation of geothermal systems (Cheng 1979;Mahmoudi, Hooman & Vafai 2019), predicting the movement of groundwater contaminating non-aqueous phase liquids such as chlorinated organic solvents (Taylor et al. 2001;Bear & Cheng 2010) and monitoring the motion and storage security of carbon dioxide (CO 2 ) linked with large-scale carbon capture and storage projects (Huppert & Neufeld 2014). The latter have recently received much attention with the growing consensus that large quantities of anthropogenic CO 2 emissions will need to be captured and stored in the subsurface to meet global emissions targets (IPCC 2018;IEA 2020).
The release of a contaminant or the injection of CO 2 results in fluid travelling through the reservoir as a result of gravitational forces until it reaches a confining layer. In typical CO 2 injection scenarios, the buoyant CO 2 rises until it reaches a structural seal, a low permeability rock layer such as a shale, anhydrite or salt, which prevents the fluid from migrating beyond the storage reservoir. For geological carbon storage to be successful, it is vital that the CO 2 remains securely trapped in the subsurface over long time scales (Metz et al. 2005). Defects within the seal may allow stored CO 2 to leak into overlying aquifers and eventually to the surface. Therefore, it is important to understand how these defects may contribute to the migration of CO 2 . Likewise, for contaminants in the subsurface the influence of defects in the reservoir seal may impact the dispersal of the contaminant.
Faults, which are localised zones of brittle deformation, are a common feature in geological reservoirs, and may act as leakage pathways for trapped fluids. There is a large body of work on the structure and fluid flow properties of fault zones (Caine, Evans & Forster 1996;Faulkner et al. 2010;Nicol et al. 2017). A simple model for the structure of a fault zone is a fault core, generally consisting of very fine grained crushed rock (gouge) and larger fragments of broken up host rock, surrounded by a heavily fractured damage zone (Wibberley, Yielding & Toro 2008). The fault core and surrounding damage zone have very different hydraulic properties. Laboratory measurements on samples of fault core show that the permeability can be reduced by two to three orders of magnitude compared with the unfaulted host rock (Zhang & Tullis 1998;Shipton et al. 2002;Shipton, Evans & Thompson 2005). In contrast, the presence of fracture networks in the fault damage zone can increase the permeability of the host rock by two to three orders of magnitude (Simpson, Guéguen & Schneider 2001;Oda, Takemura & Aoki 2002;Mitchell & Faulkner 2008). The resultant model for fluid flow in fault zones is a barrier-conduit system, where the fault core acts as a barrier to across-fault flow and the fracture damage zone channels flow parallel to the fault plane (Caine et al. 1996;Balsamo et al. 2010). Studies have applied this model to explore vertical exchange flows through faults between multiple aquifers (Woods et al. 2015) and upward migration of CO 2 into overlying permeable formations (Chang, Minkoff & Bryant 2008;Kang et al. 2014).
Faults are relatively small features in basin-scale models, so can be computationally challenging and expensive to model accurately using standard multi-scale numerical flow simulations (Class et al. 2009;Nordbotten et al. 2009). One option is to develop analytical models that describe fault behaviour, which are then integrated into the larger-scale models (Nordbotten & Celia 2011;Kang et al. 2014). These models are constrained by the assumptions needed to solve the mathematical system of equations, but solutions can be obtained in seconds as opposed to hours/days, and developing these models also provides insight into the dominant physical processes in the system.
A number of studies have focused on analytical solutions for buoyancy-driven flows or on gravity currents in porous media with leakage from the current. This includes gravity currents in unconfined porous media leaking steadily through a permeable boundary (Acton, Huppert & Worster 2001;Pritchard, Woods & Hogg 2001;Pritchard & Hogg 2002), or leaking through discrete fractures, where the leakage flux is driven by the hydrostatic pressure of the underlying less dense fluid (Pritchard 2007) or with the added effect of the buoyancy of fluid in the fault (Neufeld, Vella & Huppert 2009). Previous studies have also considered leakage in confined aquifers (Avci 1994;Nordbotten, Celia & Bachu 2004 where leakage is driven by an increased pressure due to injection or where a background pressure gradient between multiple aquifers contributes to driving leakage through fractures (Pegler, Huppert & Neufeld 2014b). In all of these studies, the leakage flux is driven by the Darcy velocity in the leaking boundary or fracture, which is a product of the mobility of the injected fluid in the layer or fracture and the pressure gradient across it.
The movement of fluids through porous rocks and small-scale fractures can be modelled by considering miscible flows in porous media, where the flow is governed by Darcy's law. Miscible flows in porous media can be characterised by their Péclet number, Pe = d 0 Uτ/D d , which describes the relative importance of advection to diffusion in fluid transport. Here, d 0 is the mean grain diameter, U is the characteristic velocity of the flow, τ is the tortuosity of the porous medium and D d is the molecular diffusion coefficient. A composite transport coefficient D can model the diffusive and dispersive contributions to movement of fluid, where (Houseworth 1984;Delgado 2007). In diffusion dominated flows, where Pe 1, the transport coefficient D D d /τ is constant whereas in advection dominated flows, where Pe 1, D d 0 U is dependent on the flow speed. As a buoyant fluid is injected into a reservoir or escapes through a fracture into an adjacent aquifer, it initially forms a buoyant plume. The behaviour of a buoyant plume in porous media was first studied by Wooding (1963), who mathematically modelled the dynamics of a rectilinear plume in porous media for Darcy flow. He derived a similarity solution that predicts the amount of entrainment from the ambient and therefore describes the increasing volume flux of the plume with height. Sahu & Flynn (2015) derived a new similarity solution for Darcy plumes with large Péclet numbers to obtain expressions for the plume volume flux and mean reduced gravity as functions of vertical distance from the source.
In this study we develop an analytical model to describe the dynamics of leakage through a fault zone cutting multiple aquifers and seals, which we test against a new set of laboratory experiments. To model flow in faults cross-cutting multiple aquifers, we combine current analytical models for a buoyant plume in a semi-infinite porous medium (Sahu & Flynn 2015) and a leaking gravity current (Pritchard 2007;Neufeld et al. 2009) with a new model for fault leakage which accounts for increased pressure gradients within the fault due to an increase in Darcy velocity directly above the fault. Previous studies providing analytical solutions to fault leakage problems have used numerical solvers to verify their models (Kang et al. 2014). Here, we test the results with a novel set of porous media tank experiments which show a good fit to the model. The results of the modelling are illustrated by application to a naturally occurring CO 2 -charged aquifer system at Green River, Utah, using an extension of the model to calculate the fluid distribution across multiple vertically stacked aquifers, cross-cut by a fault.
In § 2 we formulate a theoretical model for the half-space plume in which the plume volume flux and mean reduced gravity are used as inputs in a model for a leaking gravity current and where leakage through the fault is driven by the hydrostatic pressure within the underlying gravity current and the buoyancy of the fluid in the fault. In § 3 we show numerical results for the gravity current shape and leakage rates. Observations from experiments show that an increase in Darcy velocity within the secondary plume above the fault leads to enhanced pressure gradients within the fault. In § 4 an expression for the pressure above the fault is coupled with the leakage model presented in § 2, resulting in a new model for leakage through the fault. The model is matched against results from a set of analogue porous medium tank experiment in § 5, and in § 6 we demonstrate the application of the model in CO 2 storage by applying an extension of the model to predict CO 2 leakage across multiple aquifers. In § 7 we consider the validity of our model and discuss potential limitations and extensions including the explicit inclusion of multiphase flow relationships in the plume, gravity current and fault leakage models. Finally, we present the conclusions from this study in § 8.

Theoretical model
In this section, we present a model for calculating the leakage flux from an aquifer intersected by a fault. To achieve this, we couple a model for a buoyant plume in a semi-infinite porous medium with a model for a two-dimensional gravity current spreading under an impermeable base containing a fault of finite gap width and thickness, and known permeability and porosity.

Buoyant plume in a semi-infinite porous aquifer
Here, we derive a solution for the steady flow of a two-dimensional buoyant plume in a semi-infinite porous medium with unit thickness in the third dimension. A constant input flux q 0 of fluid with density ρ 0 is injected into a porous medium with permeability k and porosity φ saturated with a denser ambient fluid of density ρ a . Both fluids are assumed to be miscible, with equal viscosities such that there are no capillary effects during the flow. The implications of this assumption for application to CO 2 -water systems is discussed in § 7. The injected fluid mixes with the ambient fluid and forms a plume with a volume flux Q(z) and density ρ(x, z) at a given height z (figure 1a). We follow the analysis of Sahu & Flynn (2015) who considered an unconfined plume in a porous medium for Darcy flow with Pe 1, in contrast to Wooding (1963). Our analysis differs by considering a half-space model, with an impermeable vertical boundary contacting the injection location. Assuming steady Boussinesq flow (shown in previous studies to be a useful approximation in CO 2 -water systems (Soltanian et al. 2016;Amooie, Soltanian & Moortgat 2018)), the governing equations based on mass continuity, momentum continuity, solute transport and a linear equation of state are where u and w are the fluid velocities in the x and z directions respectively, P is the fluid pressure, ν is the kinematic viscosity, g is the acceleration due to gravity, C is the solute concentration, β is the solutal expansion coefficient and D L and D T are the longitudinal and transverse dispersion coefficients, respectively. Since the plume is long and thin we may neglect vertical variations in the horizontal velocity (∂w/∂x ∂u/∂z) and longitudinal dispersion. Hence, the combined momentum equations and the solute transport equation become, (2.7) We now introduce a streamfunction ψ such that w = ∂ψ/∂x and u = −(∂ψ/∂z). For Pe 1, the transverse dispersion coefficient may be approximated as where α is the transverse dispersivity (Delgado 2007 For a steady plume, the buoyancy flux, or equivalently the solute mass flux, remains constant with height, as the horizontally entraining ambient fluid only increases the volume of the plume but does not alter the solute mass. The buoyancy flux of the plume per unit thickness is conserved with height and is where g = g(ρ a − ρ)/ρ a ≡ gβC is the reduced gravity. A scaling analysis of (2.9), (2.10) and (2.11) suggests that ψ (2.12c) Equation (2.12b) motivates us to define a self-similar horizontal length scale of the plume where we note that z = 0 is the point of injection and therefore entrainment of the ambient fluid and the corresponding plume equations are applicable only for z > 0. It is also worth noting that the plume is defined as the region of the flow field where the concentration C > 0, or equivalently g > 0. Now, on considering (2.12a), (2.12c) and (2.13) together, we can define the streamfunction and concentration of the plume in the forms of similarity functions F (η) and G (η) as and respectively. These expressions for ψ and C are related by (2.9) such that F (η) = G (η).
Similarly solute concentration, (2.10) implies that We solve (2.16) subject to the conditions that the solute concentration cannot be negative, and tends to the background concentration in the far field, C(x, z) = 0 and F (η) = 0 as x, η → ∞, whereas inside the plume C(x, z) > 0 and F (η) > 0. The transverse velocity against the fault is zero, such that the value of the streamfunction ψ(0, z) = 0 so F (0) = 0. With these conditions, it can be shown that where c = √ 8/π is a constant of integration which is obtained by applying the solutions for w (= ∂ψ/∂x) and g (= gβC) to (2.11). Therefore, the expression for the volume flux per unit thickness of the plume is This result is a factor of √ 2 smaller than the volume flux of the unconfined plume given by Sahu & Flynn (2015). Using (2.18), the average reduced gravity across the plume can be calculated as a function of height, The equations here are for an ideal plume formed purely by a buoyancy flux, where the volume flux Q → 0 as z → 0. For a non-ideal plume with a finite volume flux these assumptions do not hold. We can correct for this by extrapolating the flow to negative where q 0 is the volume flux into the system, and the buoyancy flux F 0 = q 0 g 0 where g 0 is the reduced gravity of the injected fluid. Given these corrections for source conditions, expressions for the plume volume flux per unit thickness and mean reduced gravity are (2.22) 2.2. Gravity current with leaking fault at the origin Here, we consider the behaviour of a two-dimensional density-driven flow of a fluid in a porous medium of permeability k and porosity φ saturated with an ambient fluid of higher density ρ a and bounded on one side by an impermeable vertical boundary (figure 1b). The injected fluid spreads below an impermeable horizontal baffle containing a fault of gap width d f , thickness h f , permeability k f and porosity φ f which abuts the boundary. At the horizontal baffle, a gravity current forms due to fluid input from a rising plume. The fluid density ρ 1 and input rate q into the gravity current are calculated using the expressions for the plume volume flux and mean reduced gravity derived in (2.21) and (2.22). The values of q andḡ are evaluated at a vertical distance H from the initial injection point, which corresponds with the level of the base of the impermeable baffle. We assume that the depth of the ambient fluid is large compared with the thickness of the current (H h). This means we can neglect the effects of flow in the ambient and assume the values of q andḡ remain constant with time.
Using (2.21) and (2.22), we obtain expressions for the input flux into the current and reduced gravity of the current, (2.24) At the impermeable baffle, the injected fluid forms a gravity current, with some fluid leaking through the fault. We assume that the current is long and thin, and the velocity is predominantly parallel to the baffle so that the pressure in the current is hydrostatic. The fault is modelled using the barrier-conduit system described in § 1, with the fault damage zone within the baffle modelled using an average permeability k f and width d f . Neufeld et al. (2009) formulated a model for drainage through a fissure of given permeability and width, which incorporates both flow driven by hydrostatic pressure as well as the buoyancy of the fluid in the fault itself, . (2.25) Here, k f is the permeability of the fault, h f is the length of the fault, μ is the dynamic viscosity and ν is the kinematic viscosity of the leaking fluid, h 0 is the thickness of the current at x = 0 (h 0 = h(0, t)) and g is the reduced gravity of the current. Note that, despite the width of the fault, d f > 0, we assume that its effect on the gravity current is localised to the point x = 0, and the leakage is driven by the thickness there. A sharp interface between the injected and ambient fluids is described by a thickness h(x, t) below the baffle at z = H. Dispersion will occur predominantly at the edges of the gravity current as it propagates into the reservoir and is expected to alter the shape of the current, the implications of which are discussed during comparison with the experimental results in § 5.3. However, where the plume feeds the current directly below the fault, dispersion will be less pronounced. Hence, we neglect dispersion in the gravity current as it does not significantly affect leakage through the fault. The flow in the gravity current is driven by gradients in the hydrostatic pressure, and the evolution of the height is determined by the divergence of the fluid flux, (2.26) The current forms when the plume impacts the impermeable baffle, and has initial condition Subsequently the gravity current is fed by the plume and has boundary conditions, which describe the input flux at the origin (equal to the plume input flux q, minus the fault leakage flux q F ), a no flux condition through the nose of the current and zero thickness at the nose of the current. The current satisfies global conservation of mass given by, where the final term comes from (2.25) and is equal to the total volume of fluid that has leaked through the fault. Note that time t = 0 in (2.29) is the instance the plume first reaches the impermeable baffle at z = H.

Non-dimensionalisation
Based on (2.23a,b), (2.26) and (2.29) we define the following dimensionless variables, (2.30c) By substituting (2.30a-c), (2.26)-(2.29) become, We introduce the dimensionless parameter that characterises the strength of leakage through the fault due to the buoyancy of fluid in the fault so that (2.33) now has the form, Here, λ = 0 describes the case where the leakage through the fault is only driven by the hydrostatic pressure within the underlying gravity current, and not by the buoyancy of fluid in the fault.

Numerical solutions
We solve for the full time-dependent behaviour of the gravity current numerically using a finite difference scheme. The thickness profile of the current is initiallyh(x, 0) = 0, and the subsequent evolution is described by (2.31)-(2.35). The thickness atx = 0 is obtained by solving (2.35) at each time step, using the thickness profile obtained from solving (2.31) and calculating the new leaked volume of fluid usingh 0 from the previous time step.
Results obtained are shown in figures 2-4. Figure 2(a) illustrates the solution for the gravity current shape as a function of time. The thickness profiles of the current are plotted fromt = 0 tot = 50 at intervals of 5 with λ = 0, which describes the case where leakage through the fault is only driven by the hydrostatic pressure within the current below the fault. Figure 2(b) shows three height profiles att = 0.5,t = 5 andt = 50 (solid lines), where the extent and height of the profiles have been normalised by the maximum extentx N and the thickness of the current atx = 0, h 0 . The self-similar solution for the propagation of a gravity current through a porous medium with a constant input flux and with zero leakage (Huppert & Woods 1995) is also plotted (dashed line). The solutions for the leaking gravity current deviate from the zero leakage case over time, demonstrating a non-self-similar behaviour of the leaking current.
In figure 3(a), the evolution of the horizontal and vertical extent of the leaking gravity current is plotted for λ = 0 (solid lines). At early times, the horizontal extent evolves following at 2/3 power law relationship and the vertical extent of the current evolves following at 1/3 power law relationship which agrees with the solutions for a gravity current with no leakage (Huppert & Woods 1995). At late times, the horizontal and vertical extent evolution deviate from these power law relationships. The vertical extent of the current tends to a constant value, resulting in a constant rate of leakage from the current. Late time asymptotic behaviour of (2.28) and (2.29) indicates that the thickness at the origin h 0 asymptotes to a constant as the gravity current nose tends towards infinity. It is possible to show (by considering a small perturbation to this equilibrium point) that the nose position approaches infinity like x N ∼ (qt) 1/2 in this late time regime (see Appendix A). This late time behaviour can be seen if the gradients of the log-log graph for the vertical and horizontal extent of the current are plotted as a function of time (figure 3b).  The total injected volume of fluid into the system increases linearly with time and partitions between the gravity current and any leaked volume through the fault. A quantitative understanding of this partitioning is a key metric for the storage security of injected fluids in the subsurface. The total injection volume, gravity current volume and total leaked volume, represented by the termst, Initially, the volume of fluid going into the gravity current is greater than the volume of the fluid leaking through the fault. However, as the height of the current below the fault increases, a larger hydrostatic pressure drives more fluid through the fault and the volume of fluid leaking becomes larger than the volume of fluid going into the gravity current. This transition can be seen when the fluxes of fluid going into the gravity current

A simple model for the pressure above the baffle
The theoretical model described in § 2 accounts for leakage due to the hydrostatic pressure in the underlying gravity current and the buoyancy of the fluid in the fault. In several of our laboratory experiments (presented in § 5), we observe a thinning of the secondary plume that forms above the baffle close to the fault (figure 5c). The thinning may be caused by an acceleration of the secondary plume, leading to differential negative pressures near the fault which increase the total pressure gradient across the fault. Transport of concentration in the secondary plume is similar to that described in § 2.1, but the length scale over which this dispersion plays an important role is much larger than the length scale for the thinning of the plume (figure 11b), hence we do not include these effects. In this section we create a model for the pressure above the baffle by considering an asymptotic expansion in terms of a small deviation to the plume inlet velocity.

Thinning plume model
The scenario we consider is illustrated in figure 5(a). Note the redefined position of z = 0. The secondary plume z ≥ 0 is fed by a vertical inlet velocity w f which is slightly smaller than the buoyancy velocity w b = kΔρg/μ, such that w f = w b (1 − ) for some small parameter 1. As a result, the width of the plume, which we denote d(z), must reduce from an initial value d f to a far-field value d b = w f d f /w b . Note that it is possible for to be negative, in which case the width of the secondary plume would increase. We model the flow into the secondary plume using Darcy's law and mass conservation, along with boundary conditions corresponding to constant inflow across the fault, impermeability at the vertical left-hand wall and far-field velocity conditions respectively, where u and w are the fluid velocities in the x and z directions. At the edge of the steady plume, we impose a kinematic boundary condition and we enforce that the pressure equals the ambient hydrostatic The pressure at the edge of the plume is hydrostatic (4.5), since it must match with the pressure in the adjacent static ambient fluid. We therefore consider an asymptotic solution of the form u = û(x, z) + · · · , (4.6) w = w b + ŵ(x, z) + · · · , (4.7) p = p a − ρ a gz + p(x, z) + · · · , (4.8) for 1. Clearly in the limit → 0 the leading-order terms in (4.6)-(4.9) satisfy (4.1)-(4.5), as expected. At first order in , the pressure must satisfy the linear system 14) Conservation of mass dictates that the pressure must satisfy Laplace's equation, which since it is linear results in (4.10). Likewise, the boundary conditions (4.2)-(4.4) and the hydrostatic condition in (4.5b) are all linear equations. Hence, they must apply to leading-and first-order pressure terms alike. The leading-order pressure solution (first part of (4.8)) satisfies all of these trivially. The only nonlinear equation is (4.5a), also known as the kinematic condition. In this case, one must expand out the variables, keeping only first-order terms, which gives (4.15). Conveniently, (4.10)-(4.14) can be solved independently of (4.15), indicating that the plume shape and the pressure are decoupled. The pressure solution is calculated by separation of variables, givinĝ Hence, the pressure (including both leading-and first-order terms) at x = z = 0, which we denote p − , is given by where G = ∞ n=0 (−1) n /(2n + 1) 2 ≈ 0.9160 is Catalan's constant. The plume shape is calculated by evaluating (4.15) and (4.16), which converges to the differential equation We solve (4.18) with boundary conditiond(0) = 0 and plot the solution in figure 6(a). Clearly,d → −d f as z → ∞, such that the the total plume shape d f + d asymptotes to the far-field plume width d b = d f w f /w b , as required. We note that both the pressure p − and the inlet velocity w f are unknown in (4.17). Hence, to close the system we require a second equation relating these two quantities. This is given by considering the flow through the fault −h f < z ≤ 0. In particular, as illustrated in figure 5, the flow through the fault is driven by a difference in pressure p + − p − , where p + is set by the thickness of the gravity current p + = p a + Δρgh 0 + ρ a gh f . (4.19) Hence, by approximating Darcy's law across the fault we arrive at a second relationship, Rearranging (4.17) and (4.20), we arrive at dimensionless equations for w f /w b and p − ,  The pressure p − (written in dimensionless form ( p − − p a )/Δρgd f ) is illustrated in figures 6(b) and 6(c) using contour plots. Largest negative pressures are observed for small values of the gravity current thicknessh 0 , or large values ofh f andk (corresponding to small velocity ratios w f /w b ). We note, however, that we do not expect our asymptotic solution to be valid for w f /w b 1. In such situations a full numerical simulation is required to determine p − . Nevertheless, for the parameter range of the current study (w f /w b ∈ [0.5, 1.2]) we expect (4.22) to remain a good approximation.

Laboratory study
We conducted a series of laboratory experiments to test two aspects of our theoretical predictions. First, the spatial distribution of the gravity current was measured as a function of time. Second, the partitioning of fluxes between the gravity current and the fault was measured by calibrating the image intensity and dye concentration. As it was easier to model experimentally, we considered a system in which the injected fluid is more dense with respect to the ambient. Given the Boussinesq approximation, this system behaves in a symmetrical manner to a less dense fluid rising in an aquifer saturated with a denser fluid. Exp.
Symbol The experimental images presented are rotated by 180 • for ready comparison between the analytical and experimental results.

Experimental set-up
The experiments were performed within a Perspex cell of length 40 cm, height 70 cm and internal thickness 1 cm (figure 7). The cell was filled with glass ballotini which formed a porous layer. The glass ballotini filling the majority of the cell had diameter b 0 = 3.1 ± 0.2 mm except for the region adjacent to a central plastic spacer that lies across the cell where the ballotini diameter b f = 1.0 ± 0.1 mm. The porosity of 3 mm glass ballotini in a cell of width 1 cm was measured by Sahu & Neufeld (2020) and found to be φ = 0.41 ± 0.01, which is slightly higher than the value of φ = 0.37 for randomly close-packed beads as the width of the cell is comparable to the bead size. The permeability of the larger porous medium filling the majority of the cell was estimated using the Kozeny-Carman relation (5.1) The Kozeny-Carman relationship breaks down when applied to a small volume of beads, hence the permeability of the region with smaller ballotini was treated as an unknown parameter and fitted for (see § 5.3). The entire cell submerged in a large water tank with dimensions 90 cm × 80 cm × 80 cm filled with fresh water of density ρ a = 0.998 g cm −3 . The cell was closed on three sides, with one side open but lined with a permeable mesh so that this unconfined side created a hydrostatic pressure boundary condition at the right-hand side of the cell. A total of 15 experiments were conducted, using four different fault geometries (A-D), as summarised in table 1. During the experimental run, a peristaltic pump was used to inject aqueous solutions of sodium chloride (brine) dyed with red food colouring into the cell at a constant rate q 0 which ranged from 0.039-0.245 cm 3 s −1 , and was calculated by measuring the mass  of the brine container over the course of the experiment. The injected brine solutions had an initial dye concentration of c 0 = 3.00 ± 0.05 g l −1 and an initial density ρ 0 = 1.030 or 1.070 g cm −3 , assuming a water temperature of 20 • C (Green & Southard 2019). The kinematic viscosity of the water was 0.01 cm 2 s −1 and the transverse dispersivity of the glass ballotini was assumed to be α = 0.005 cm (Delgado 2007). The experimental parameters were selected so experiments spanned a range of λ values. We used a Nikon D5000 DSLR camera with a resolution of 4288 × 2848 pixels to capture images over the course of the experiment with a time gap which ranged from 4 to 10 s. The cell was backlit using a LED light panel with a Perspex diffuser to ensure uniform illumination.

Post-processing scheme
A set of photographs for experiment B2 is shown in figure 8, rotated by 180 • for ready comparison between the analytical and experimental results. The theoretical predictions for the position of the gravity current interface h (x, t) in the bottom portion of the cell is plotted for the thinning plume model, along with the interface for the zero leakage model (Huppert & Woods 1995). The time t = 0 is defined as the point where the plume first makes contact with the bottom of the central spacer. The horizontal extent of the gravity current was measured as a function of time by picking the furthest front of the current  (Huppert & Woods 1995). Images are rotated 180 • to allow comparison with theoretical results (see § 5).
above a threshold value. The error range of the measurements was set by calculating the extent for a range of threshold values. The results were verified by manually checking the front location from photographs and showed excellent agreement. The height of the current was more difficult to interpret due to dispersion from the plume making the top of the current uncertain.
To measure the flux of fluid leaking through the fault, the concentration of the injected fluid throughout the cell must be determined, and so a series of calibration experiments were performed to determine the functional relationship between the image intensity and dye concentration. In each calibration experiment, the cell was uniformly saturated with a red dye solution with concentration C d . The light intensity for the calibration image was calculated by subtracting the light intensity from a reference image containing no dye. A total of 20 different concentrations between 0 and 3.00 ± 0.05 g l −1 were used to construct calibration curves for the green and blue colour channels (figure 9). The dye calibration is more sensitive to different colour channels at different dye concentrations, so a hybrid calibration curve was used which weighted contributions from the green and blue curves. The full function form of the hybrid calibration curve is given in Appendix B.
The mass of dye across the cell is calculated by summing the average concentration within 0.5 cm × 0.5 cm sized bins, with the leaked mass equal to the total mass below the top of the baffle. The total error in measurement is the sum of two errors. First, the leaked dye mass shows a linear behaviour at late times. A root-mean-square deviation from a regression line was calculated, with the deviations likely caused by small variations in light intensity between each photo, for example due to trapped air bubbles. Second, the total measured mass of dye across the cell was compared with the known input dye mass for each experiment. The measured mass increases linearly as a function of time but the post-processing recovers 83 %-97 % of the input dye mass across experiments. The partial measurement of the input dye is likely due to sensitivity of the calibration curve to smaller concentrations. For each experiment, the measured dye mass is scaled to the known input dye mass and the difference is defined as the error in measurement.  Figure 10(a) shows the total leaked dye mass (black circles) and leakage flux (pink circles) as a function of time for experiment A3. The leaked flux is calculated by taking the gradient of a quadratic polynomial fitted over a seven element moving window. The errors bars represent a 95 % confidence interval for the regression coefficients.

Experimental results and comparison with theory
The experimental results are compared with theoretical results from three different models presented in § § 2 and 4. The first model only considers the contribution from the hydrostatic pressure within the underlying current (λ = 0). The second model considers contributions from the underlying current and the buoyancy of fluid in the fault (λ / = 0), and the third model includes both of these effects but also considers the contribution of flow-enhancing pressure deviations directly above the fault (thinning plume model). There is very little difference between the solutions for the second and third models, suggesting that the effects of pressure deviations above the baffle due to thinning of the secondary plume are minimal. This agrees with experimental observations (figure 10b), which show that the width of the secondary plume is the same as the width of the fault (d b /d f ≈ 1).
Due to the difficulty in estimating the fault permeability k f , this was calculated by minimising the misfit between the experimental data and the thinning plume model using a least-squares regression for one of the experiments. Experiment A3 (figure 10) was selected to fit the fault permeability due to the agreement between the thinning plume model and the λ / = 0 model. A value of k f = k 0 /3.9 ≈ 2.7 ± 0.3 × 10 −9 m 2 was obtained and this value was used to calculate the theoretical results across the other experiments. It is possible that there is some variation in the fault permeability across different experiments, but this should be small as the faults are a similar size and were all packed using the same method. The sensitivity of the thinning plume model to changes in fault permeability is shown by the grey dashed lines in figure 10(a), in which the model solutions are displayed for k f ± 20 %. Figure 11(a) shows the total leaked dye mass (black circles) and leakage flux (pink circles) as a function of time for experiment D4, along with theoretical results for the three models. There is a clear difference between the model which allows for a thinning plume above the baffle and the model only considering contributions from the underlying current and buoyancy in the fault, suggesting that acceleration and thinning of the secondary plume is leading to enhanced pressure gradients within the fault. Experimental observations are in agreement (figure 11b), where significant thinning of the secondary plume can be seen above the fault (d b /d f < 1). Note that the model considering buoyancy in the fault assumes that the fault is initially full of the injected fluid. Furthermore, the thinning plume model assumes that a secondary plume has reached a quasi-steady state profile above the fault. In reality, at early times the fault remains filled with the ambient fluid so the only driving force is the hydrostatic pressure in the underlying current. To account for this, a breakthrough time is introduced, defined as the time it takes for the injected fluid to fill the fault and breakthrough into the upper reservoir, which is obtained from experimental observations. Prior to the breakthrough time, all three models only consider the contribution from the underlying current. After the breakthrough time, the other driving forces are taken into account. Further details on how the breakthrough time is calculated and the sensitivity of the model to the breakthrough time is given in the supplementary material available at https:// doi.org/10.1017/jfm.2021.970. This modification to the theory results in a discontinuous leakage flux profile (pink dashes, figures 10(a), 11(a)), where the leakage flux increases rapidly at the breakthrough time. These predictions broadly agree with the behaviour seen in the experimental data, where initial leakage rates are slow as the fault fills with the injected fluid before increasing rapidly to a higher and constant rate. The thinning plume model shows good agreement with the experimental data in both total leaked dye mass and leakage flux. The horizontal extent of the underlying gravity current along with solutions from the thinning plume model is plotted in figure 13 for experimental set-ups A-D. In general the theoretical model shows good agreement with the experimental data, although it tends to overpredict the horizontal extent of the current. A possible explanation is that the gravity current entrains ambient fluid as it moves into the reservoir which is not captured in the model. This decreases the reduced gravity of the current, and hence reduces the distance it travels. This effect would also cause the thickness of the current to increase, and this can be seen in figure 8, where the experimental current has a greater thickness than the thinning plume model predicts, but appears more dispersed towards its outer edge.

Application to a CO 2 storage reservoir
In reservoirs with multiple faulted aquifers and seals, a fraction of the injected fluid will flow into the aquifers at each level and a diminishing amount will continue to leak through the fault. When applied to geological CO 2 sequestration, this has important implications for providing an initial estimate of the security of a potential field-scale storage site, especially where numerical reservoir simulations are unable to include the details of faults. The critical properties of the faults such as their damage zone widths and permeabilities may be estimated from scaling relationships, analogue surface outcrops or drill cores if available.
To address the case of a fault cutting through multiple layers, we now extend the model developed in § § 2 and 4 and briefly discuss the results. Three aquifers and seals are stacked vertically with a fault cutting through the layers (figure 14a). The same physical parameters are used for each layer, and are modelled on a natural CO 2 -charged aquifer at Green River,  (Kampman et al. 2014). The system is initially saturated with water of density ρ a = 1000 kg m −3 . Less dense CO 2 is injected at a constant rate of 0.1 Mt yr −1 at the bottom of layer 1, and assumed to have constant density ρ = 790 kg m −3 and dynamic viscosity μ = 6.9 × 10 −5 Pa s, calculated at typical storage reservoir conditions of 15 MPa and 40 • C (Dubacq, Bickle & Evans 2013). The mass flux of injected CO 2 is converted into a two-dimensional input flux assuming an along-fault system length of 1 km. These parameters give an initial value of λ = 0.26. The total CO 2 mass in each layer as well as the total leaked CO 2 (defined as the leakage from layer 3) is plotted as a function of time for three different values of the fault permeability (figure 14b). When the fault permeability is comparable to or less than the reservoir permeability, the majority of the CO 2 remains trapped within layer 1 after a 1000 year time scale. However, when the fault permeability is larger than the reservoir permeability, a significant fraction of the injected CO 2 leaks all the way through the system. For example, in the two cases where k f /k = 5 and k f /k = 10, layers 2 and 3 do not accumulate significant amounts of CO 2 . However, even with k f /k = 10 % ∼ 55 % of the CO 2 remains trapped after 1000 years. Note that this is the simple case where all layers have equal permeabilities. In a more realistic system where the permeability varies across layers, it is likely that larger accumulations of CO 2 would be seen in more permeable aquifers, regardless of their position. This system represents the worst case scenario for CO 2 storage as trapping mechanisms such as capillary trapping and dissolution trapping are neglected. In reality, as the CO 2 rises and spreads out into different aquifers these other mechanisms would contribute towards long-term storage of the CO 2 .

Discussion
We have presented a new analytical model that describes the dynamics of leakage through a fault zone given properties relating to fault, aquifer and fluids. This new model has been tested against results from analogue porous medium tank experiments, from which we obtain good fits with the experimental data. This model has then been applied to predict the CO 2 distribution and leakage rates from a naturally occurring CO 2 reservoir. We find that when the permeability of the fault is comparable to the reservoir, the majority (>96 %) of the CO 2 remains trapped after 1000 years. However, it is important to highlight the limitations of the current model and discuss to what extent it can be applied directly in its present form, and in which ways it might be fruitfully enhanced as part of a future study. The gravity current model uses a sharp-interface approximation (Huppert & Woods 1995;Hesse et al. 2007), originally developed for miscible flows in porous media. This is valid for settings such as saline intrusions, but extra multiphase effects need to be considered when applied to immiscible settings such as CO 2 -water systems (Golding et al. 2011). When CO 2 is injected into water-saturated reservoirs, capillary forces play a key role in determining the saturation distribution and flow properties through the relative permeability and capillary pressure relationships (Nordbotten & Celia 2011). When the fault and reservoir have relatively uniform permeability, we expect the behaviour of the system to scale in a similar manner to the miscible case, with the effective permeability instead replaced by the product of intrinsic and relative permeability. However, geologically heterogeneous reservoirs may have significant gradients in capillary pressure, which tend to rearrange the CO 2 saturation into high permeability regions, thereby accelerating plume migration (Jackson et al. 2018;Benham, Bickle & Neufeld 2021). An extension of the present model to fault zones with significant capillary heterogeneity therefore remains an outstanding and important question in predicting the flux through such systems.
Another potential refinement to the model is that the fault zone will have a capillary entry pressure that needs to be overcome before leakage can occur and so the underlying current may have to build up to a large thickness before breakthrough can occur. It is also possible that the CO 2 will be constrained to more localised flow paths with lower entry pressures which, depending on the entry pressure and permeability of the fault and intervening reservoirs, may affect the vertical migration of CO 2 . These may be simply accommodated by the addition of a critical height the gravity current must overcome to achieve breakthrough in the fault (Woods & Farcas 2009;Sayag & Neufeld 2016).
The present gravity current model also assumes that the aquifer is unconfined (h H). This means that there is negligible flow of the ambient fluid and so fluid propagation is independent of the viscosity ratio between the fluids (Huppert & Woods 1995;Vella & Huppert 2006). Pegler, Huppert & Neufeld (2014a) showed that background pressure-driven flow due to confinement starts to play a role when ratio of the current depth and aquifer depth (h/H) is comparable to the viscosity ratio of the injected and ambient fluid. In the case of CO 2 and water, the viscosity ratio μ CO 2 /μ w 0.1, hence the assumption that the size of the reservoir is much larger than the thickness of the current is valid for currents up to ∼10 % of the height of the aquifer, as is the case at many CO 2 storage sites. As the experiments performed in this study were between fluids of near equal viscosity, the thickness of the currents across all experiments are well within the limit for the unconfined approximation. For cases in which h CO 2 0.1H, Pegler et al. (2014b) found that confinement causes greater leakage at earlier times due to the introduction of background pressure, but at later times, as the CO 2 fills the entire depth of the aquifer, the maximum hydrostatic head below the fracture is limited and hence the leakage rate is also limited. Incorporating the effects of confinement on the gravity current therefore also presents an additional extension of the present work.
There are other important fluid dynamical processes which stabilise the storage of CO 2 , including for example the dissolution of CO 2 into the brine. Dissolution depends on the relative flows of CO 2 and brine over the small (e.g. centimetre) length scales related to CO 2 diffusion, presenting significant modelling challenges. However, the mixing of CO 2 and brine within the rising plumes and fault zone is likely to enhance dissolution of the CO 2 (cf. Kampman et al. 2014).
Constraining the hydraulic properties of fault zones in the field remains a much studied area with many different factors affecting the flow of fluids, such as the host rock composition (Wibberley & Shimamoto 2003), fracture density (Mitchell & Faulkner 2012) and fracture connectivity (Saevik & Nixon 2017) and their contributions to the overall permeability. The model presented in this study constrains the sensitivity of CO 2 leakage rates to bulk fault properties such as permeability, width and thickness. By observing the resultant plume behaviour around complex fault zones, this simplified model could be used to provide an estimate of these bulk properties. When applied in a suitable geological context, the model presented here can be used to characterise the flow dynamics of fault zones.

Summary and conclusions
We have presented an analytical model that describes the dynamics of leakage through a fault zone, motivated by the potential risk to storage security of anthropogenic CO 2 . It comprises a two-dimensional gravity current in a porous medium, fed by a buoyant plume and spreading under a horizontal impermeable baffle. The medium is bounded by an impermeable vertical boundary and the horizontal baffle contains a fault through which the current is leaking. This system constitutes a reduced-order model of a faulted caprock, whereby an impermeable fault core is surrounded by a high permeability fracture zone (Wibberley et al. 2008) through which leakage occurs.
In the experiments, we observed thinning of the secondary plume above the fault due to an increase in Darcy velocity, leading to increased pressure gradients across the fault. A new model for leakage through the fault was derived, factoring in this effect. In contrast to previous analytical fault models which use numerical solvers to verify their models (Kang et al. 2014), we matched results using a series of analogous porous medium tank experiments, which supported the dependence of the fault width, fault height, input flux and input density on leakage rates as derived by the model. Crucially, this study showed how these parameters control the ratio of fluid flux into the aquifer compared with fluid leaking through the fault, which is significant for storage efficiency.
We demonstrated the utility of the model for assessment of storage security of anthropogenic CO 2 by application to a naturally occurring CO 2 -charged aquifer at Green River, Utah, using an extension of the model to calculate the fluid distribution and leakage rates across multiple vertically stacked aquifers, cross-cut by a fault. This showed the dependence of leakage rates on the fault/aquifer permeability contrast, with significant leakage occurring when the fault permeability is larger than the reservoir permeability. A detailed discussion of the limitations of the model in application to CO 2 -water systems was provided in § 7.
It is important to note that while other trapping mechanisms such as dissolution trapping, residual trapping and mineral trapping are important in CO 2 storage, they occur on longer time scales and so structural trapping remains the principal mechanism for storage of CO 2 . Understanding to what extent faults within caprocks can act as potential leakage pathways is crucial for ensuring safe storage of CO 2 .
Supplementary material. Supplementary material is available at https://doi.org/10.1017/jfm.2021.970. which is obtained from (A2). We consider the current thickness at x = 0, where δ(t) h 0∞ is a small perturbation. On equating (A1) and (A3) and applying the relationships (A4) and (A6), By considering the leading-order terms, we find that which indicates that the current extent follows a t 1/2 power law relationship at late times.

Appendix B. Functional form of calibration curve
A series of calibration experiments were performed to determine the functional relationship between the image intensity and dye concentration. In each calibration experiment, the cell was uniformly saturated with a red dye solution with concentration C d . A total of 20 different concentrations between 0 and 3.00 ± 0.05 g l −1 were used to construct the calibration curve. The light intensity for each dye concentration was calculated by subtracting the average green or blue colour channel across the image from the average green or blue channel across a reference image where C d = 0 g l −1 . The green and blue channel produce unique calibration curves (figure 9) which are sensitive to changes in image intensity at different concentrations. The functional forms of the calibration curves for the green and blue channels are and C B = aI B , I ≤ 106, (b + cI B ) −1/d , I > 106, where I G and I B are differences in the green or blue image intensity between the reference and calibration images and a, b, c and d are fitting coefficients. The calibration curve for the blue channel is sensitive to small changes in dye concentration for concentrations up to ∼0.9 g l −1 , but insensitive at higher concentrations due to image saturation. In comparison, the green channel is less sensitive at low concentrations but has greater sensitivity at dye concentrations above 0.9 g l −1 . The best calibration results are obtained when using a weighted average of the two curves to convert image intensity to dye concentration. The effective concentration is determined by the function where C B and C G are the concentrations calculated from the blue and green calibration curves, C avg = (C B + C G )/2, τ is the concentration at which C eff is a result of equal contributions from C B and C G (grey solid line, figure 9) and δ sets the width of the region both curves contribute significantly towards C eff (grey dashed lines, figure 9).