Simulation of the Infiltration of Fractured Rock in the Unsaturated Zone

We study infiltration of rainwater into fractured rock and the accompanying capillary exchange processes between fractures and matrix, hereafter referred to as fracture–matrix transfer (FMT). Its influence on the velocity of the wetting front for uniform and variable aperture fractures is of prime interest because it determines the penetration depth of infiltration pulses. FMT is modelled explicitly in a discrete fracture and matrix (DFM) framework realised using a hybrid finite element–finite volume discretisation with internal boundaries. The latter separate the fracture mesh from the rock matrix mesh with the benefit that the flow that occurs within the minute fracture subvolume can be tracked with great accuracy. A local interface solver deals with the transient nonlinear aspects of FMT, including spontaneous imbibition of the rock matrix. Two- and three-dimensional heuristic test cases are used to illustrate how FMT affects infiltration. For the investigated scenario, we find that—beyond a critical fracture aperture around 5–10-mm—infiltration rate is no longer affected by FMT. Fracture aperture variations promote in-fracture-plane fingering, with counter-current flow of water (downward) and air (upward). Fracture flow interacts with FMT in a complex fashion. For systems with a small fracture porosity (≤0.01%), our results suggest that intense, hour-long rainfall events can give rise to tens-of-meter-deep infiltration, depending on fracture/matrix properties and initial saturation of the fractured rock mass.


Introduction
The thicker the vadose zone, the more important are water infiltration dynamics for contaminant ingress into the saturated zone. In brittle sedimentary rocks or volcanic sequences, fractures have a major impact on this flow (e.g., Berkowitz [1], Rosenbom and Jakobsen [2]). It follows that the accurate prediction of unsaturated-zone fracture flow, and infiltration rate and depth during individual rain fall events is key for developing contamination management and treatment strategies (Robins [3], Focazio et al. [4]).
In fractured porous media, the rock matrix normally has a much higher matric potential (capillary pressure) than the fractures. Thus, upon fracture infiltration from the earth's surface, water is sucked into the rock matrix. This process is known as spontaneous imbibition.
Examining borehole data from Yucca Mountain, Nevada, USA, Thordarson [20] asserted that capillary forces immobilise infiltrating water in the rock matrix even if the fractured rock is highly saturated. On the basis of the results of dual continuum fracture and matrix simulations, Wang and Narasimhan [21] inferred that for fracture flow to occur, the rock matrix must be saturated. Their results seemed to indicate that for most parts of the vadose zone at Yucca Mountain, the impact of the fractures on infiltration could be neglected. However, Scott et al. [22], Roseboom [23], Montazer and Wilson [24] challenged this on the basis of field observations, arguing that infiltrating rain water rapidly flows down fractures and open channels with little or no lateral diversion. In their conceptual model, water only saturates the immediate vicinity of fracture walls and therefore keeps flowing downward. Following their reasoning, deep infiltration would not be limited to the scenario of a fully saturated rock matrix, and, due to the low fracture volume, even for a low precipitation rate, rain water may transgress a thick vadose zone, reaching the water table. Experiments conducted by Rangel-German and Kovscek [25], Wood et al. [26], Cheng et al. [27] show that flow rate influences fracture flow in more intricate ways. Water can either imbibe the matrix, such that the wetted halo keeps up with the saturation front in the fracture(s) or water may shoot far ahead through interconnected fractures. Fingering within the fracture plane may also contribute to the fast movement of the water front in the fractures (Pruess [28], Glass [29]). Due to the highly nonlinear process interactions, underpinning this behaviour, predicting the velocity and penetration depth of infiltration fronts in the fractured vadose zone is challenging.
Modelling this rich set of potential behaviours by numeric simulation without oversimplification was the objective giving rise to this paper. We present a new hybrid model that is validated with physical experiments. Beyond this numeric simulation research and software engineering, we apply the new scheme to a few idealised models of fractured rock to obtain first insights on the systematics of fracture infiltration in the vadose zone.
Although capillary transfer can be regarded as a nonlinear diffusion process that is more important on the small scale, our domain of interest, the vadose zone, can be hundreds of meters thick. Considering the frequency of major rainfall events in semi-arid regions, the time scales of interest are hours to tens of years. These constraints imply that physical experiments are very difficult to perform, and numeric simulation is the preferred tool to study this system. However, realistic simulation of fracture infiltration remains challenging (e.g., Makedonska et al. [30]).
In this study, we apply the alternative discrete-fracture and matrix (DFM) simulation approach (Kim and Deo [42], Paluszny et al. [43]). Lower-dimensional and higherdimensional finite elements are used to represent individual fractures and matrix blocks, respectively. Therefore, fractures in 2D models are represented by line elements and by surface elements in 3D models. In this mixed-dimensional approach described in detail by Juanes et al. [44], fracture properties are surface attributes, permitting, for instance, to change fracture aperture without modifying the mesh. This offers flexibility and greatly reduces model complexity. Such DFM models have been applied successfully to investigate single phase flow in fractured rock (Király [13], Kim and Deo [45], and Juanes et al. [44]). Kim and Deo [45] first applied this method to two-phase flow and Matthai et al. [46], Matthai and Bazrafkan [47] extended it to complex field-data based fracture models. Bazrafkan et al. [48] introduced support for saturation discontinuities at material interfaces. However, their model did not yet support the decoupling of phase pressure across material interfaces because, in this hybrid finite element -finite volume method, pressure is discretised on the finite-element nodes and these are shared by elements joined at material interfaces. This missing degree of freedom of this so-called FECFVM method can be overcome through the application of an extended version of the embedded discontinuity method of Nick and Matthäi [49,50] with additional pressure degrees of freedom along interfaces and a non-linear iterative method to compute the interface fluxes, see Tran et al. [51].
Tran's embedded discontinuity method is the point of departure of this study, and here we present a further extension of it, supporting unique phase-pressure values for the lower-dimensional fracture and the rock matrix on either side. Capillary pressure is solved implicitly across the different subdomains of the model separated by interfaces. This embedded discontinuity scheme differs from the many formulations that compute local capillary pressure derivatives only (e.g., Hoteit and Firoozabadi [52]) and often assume that the saturation on either side of the interface reflects capillary equilibrium (e.g., Helmig and Huber [53], Reichenberger et al. [54]). That the latter is not necessarily the case during fracture flow will be illustrated by answering the leading questions of this study.
The article proceeds with the description of the new DFM discretisation via nodematched mesh patches and their nonlinear coupling strategy. Next, we introduce several idealized DFM models and describe their simulation setup. The following result section analyses the nature and propagation velocity of saturation fronts. A discussion section precedes the conclusions of this article.

Mathematical Model
In earlier work, Tran et al. [51], presented a new simulation approach for the dynamic two-phase flow behaviour of material interfaces, including filtration and capillary barrier effects that come into play below the capillary threshold pressure. This paper extends this scheme, presenting a new discretisation of fractured porous media and treatment of spontaneous imbibition. To start with, we briefly recall the mathematical foundation of Tran et al.'s ( [51]) method.

Governing Equations
For simplicity, this description is restricted to immiscible, incompressible two-phase flow. No-slip boundary conditions are applied to the fracture walls, and no Brinkmann terms are used to couple fracture and matrix flow. The mass conservation equations for phase α ∈ {n, w} read Here variable φ is the porosity, S α is the phase saturation, and q α represents volumetric source/sink terms. Phase velocity v α is evaluated from phase pressure via extended where k is the intrinsic permeability tensor, p α is the phase pressure, ρ α the phase density, and g the gravity vector. Phase mobility λ α is defined as the ratio of phase relative permeability k rα and phase dynamic viscosity µ α . Phase saturation must satisfy the constraint The phase pressure and phase saturation are coupled via capillary pressure Capillary pressure is assumed to be a strictly monotonic function, increasing with non-wetting phase saturation.

Pressure Equations and Saturation Equation
An implicit pressure explicit saturation scheme (IMPES) resolves pressure and saturation sequentially. Despite its time-step limitations, this scheme has shown its effectiveness for dealing with highly dynamic systems, see Aziz and Settari [55], Chen et al. [56], and Monteagudo and Firoozabadi [57]. To weaken the coupling between phase pressure and phase saturation, the global pressure formulation is used. This eliminates the capillary pressure term from the implicitly solved pressure equation, Chavent and Jaffre [58].
The total velocity is taken as the sum of the phase velocities v t = v n + v w . In this form of extended Darcy's law, the total velocity is given by v t = −kλ t ∇p + k(λ n ρ n + λ w ρ w )g (5) where p is the global pressure. Volume conservation is expressed as The total volume source term q t is the sum of the phase volume source terms q α . The global pressure formulation is suitable for simulation regardless of whether the primary transport variable is the saturation of the wetting or non-wetting phase. The wetting-phase saturation is updated after the transport step by applying the fractional flow formulation where f α = λ α / λ t is the fractional flow, andλ = λ n λ w / λ t is mobility product. This scheme is wide spread in the simulation of two-phase flow in porous media (Helmig and Huber [53], Geiger et al. [59], Hoteit and Firoozabadi [60]). While this formulation is applied to the rock matrix, fracture flow and fracture-matrix transfer receive a special treatment.

Split-Boundary and Fracture Network Representation
The linear FEFVM collocated mesh method employs a contiguous finite-element mesh to accurately discretise different rock-matrix domains, allowing material properties to vary from element to element in a piecewise constant fashion. However, node-wise linear variables are continuous across element boundaries. As a result, this basic FEFVM scheme cannot capture discontinuities of saturation and pressure at material interfaces. To remove this shortcoming, a mesh splitting operator was applied to the input finite element mesh by Tran et al. [51]. It detects nodes located at material interfaces and multiplicates them. As many collocated nodes created as materials are joined at this mesh coordinate. These interface nodes are reconnected to the corresponding mesh patches (Figure 1). In the new development presented here, the mesh splitting operator now also recognises fractures as stand-alone mesh regions. These originally lower-dimensional mesh domains are now retained as isolated regions within split boundaries, supporting unique properties of such discontinuities in the model. To connect the fractures with the rock matrix, the fracture nodes are added to the node manifolds that identify the multiplicated nodes of the split boundaries and the dissected node-centered finite volumes that they belong to. In this way, split mesh patches representing the fractures can now communicate with each other and with the surrounding rock matrix via the finite-volume discretisation that is collocated with the finite-element mesh.

Manifold Matrix
Fracture Node Figure 1. Triangular finite elements discretising two material patches M 1 and M 2 , and node-matched line elements representing three fractures F 1 , F 2 and F 3 . Manifolds are constructed at multiplicated nodes. This also applies to lower-dimensional fracture manifolds.

Extended Interface Condition
As already mentioned in the introduction, in the extended DFEFVM, both pressure and saturation can be discontinuous at matrix-matrix and matrix-fracture interfaces. Therefore, two additional conditions are required to close the system of equations. The first condition concerns volume conservation. It imposes a constraint on the pressure equation at multiplicated nodes. The second condition is the extended pressure condition, which dictates the capillary value at the interface, Van Duijn et al. [61].
where, S * w stands for the threshold saturation at which capillary pressure in the fracture reaches the entry pressure of the matrix, Figure 2. If the wetting-phase saturation in the fracture is below this threshold value, the capillary pressure in the fracture is higher than the entry pressure of the matrix, and the first condition of (8) is enforced so that capillary pressure is continuous across the interface. When the fracture saturation is higher than the threshold value and capillary pressure in the fracture is smaller than the entry pressure of the matrix, the non-wetting phase cannot leave the fracture. However, this capillary barrier does not stop the flow in the opposite direction. In this situation, capillary pressure is discontinuous across the interface. If matrix and fracture share the same entry pressure, capillary pressure is always continuous across the interface.
For permeable coarse-grained sandstones and sub-mm fracture apertures, the capillary equilibrium condition may be satisfied. However, for larger fractures, with linear relative permeability curves and capillary pressure conditions imposed by local fracture aperture, this does not apply. In particular, when the maximum capillary pressure in the fracture p f max < p m d , the capillary entry pressure, capillary equilibrium between the fracture and the rock matrix is not possible.

Solving the Pressure and Transport Equations
The pressure Equation (6) is solved by applying the standard linear finite element method. The discretised form of this equation leads to definite integrals for each element (9) Here, N is the vector of finite element basis functions. An assembly process iterates over all the finite elements accumulating the algebraic system of equations. Solving this system delivers the global pressure at the finite element nodes. Whether the pressure across multiplicated nodes needs to be coupled or not, depends on whether there is crossflow over the corresponding interfaces. Several methods can be employed to solve this constrained system of equations, for instance, Lagrange multipliers, augmented Lagrange, and-or penalty methods, see Oden and Kikuchi [63]. Like Nick and Matthäi [49], we modify off-diagonal terms of the global conductance matrix and right-hand vector to collapse nodal degrees of freedom as necessary. Once the global pressure has been found, the total velocity is retrieved from Darcy's law in a post-processing operation. Next, the finite volume method is applied to solve Equation (7). The forward Euler time marching scheme advances the wetting phase saturation.
where ∆V is the pore volume of the control volume. The wetting-phase flux over a finite volume facet is computed by This method works everywhere in the model except at material interfaces where fluxes between disjoint mesh patches need to be computed.

Solving Matrix-Fracture Interface Saturation
Consider a manifold and its corresponding truncated control volumes (Figure 3), volume conservation written for phase fluxes exchanged between fracture and matrix without sink or source terms is Tran et al. [51] utilize a saturationbased iteration solver, where the saturation value in matrix finite volume sectors relates to the value in fracture sectors via the condition (8). This method particularly suites the drainage process in which the capillary curves are flat. However, at the saturation front, for spontaneous imbibition, the sharp gradient of capillary pressure significantly compromises solver efficiency. This motivated us to develop a new capillary pressure-based interface solver with more flexibility and a wider trust region. In this approach, capillary pressure is the primary variable. Inserting the inverse of capillary pressure, achieves these objectives. It can be solved by iterative methods such as Newton-Raphson, Regula Falsi or Brent's method [64]. While updating capillary pressure, these methods automatically evaluate the saturation values at the interface, enforcing condition (8).
The flux over the fracture stays between two material patches R 1 and R 2 . Velocity field v are constant within elements; Nodal saturation values of matrix and fracture are generally not the same.
As the water saturation decreases to its residual value, matrix capillary pressure approaches infinity, causing numerical instability. We eliminate this issue by using the linear extension of the curve to a finite value approximating a dry matrix region, see Fink [65], Campbell et al. [66], Webb [67].
Importantly, by contrast to the common Richardson equation-based simulation of flow in the vadose zone, we apply the described two-phase flow model to model water infiltration and air flow in fractured unsaturated rock. Beyond the benefit that the volume and flow of air is resolved explicitly, this approach has the major advantage that the resulting flow and transport equations are better behaved than when the water-content dependent hydraulic conductivity is used.

Conceptual Fracture Modelling
Back of the envelope calculations provide some first insights on how fast water might infiltrate an impervious rock mass intersected by vertical continuous fractures of a given aperture. The fluid is assumed to be Newtonian, and capillary (pinning) and inertia effects are ignored. It is further assumed that the fractures are fully saturated, i.e., that fluid supply does not constrain flow rate. For film flow that might be common in centimetre-aperture fractures Tokunaga and Wan [68] show that flow velocities might be up to 1000 times faster than for saturated flow. This phenomenon will not be considered here. The very much simplified scenario is intended for a preliminary verification of the more sophisticated numeric model.
For isolated opening-mode fractures (penny-shaped cracks), aperture is correlated with fracture length (Pollard and Segall [69]), following linear or sublinear trends established by many field studies (e.g., Olson [70]). Given the power-law length-frequency statistics of many natural fracture patterns (e.g., Bonnet et al., 2001), short fractures often are orders of magnitude more common than long ones, with implications for aperture statistics. An example are aperture measurements in the cross-drifts of the exploratory research facility at Yucca Mountain ( Figure 4). While the average fracture aperture is only 0.1 mm, Mogano [71]) observed sixty five 1 mm aperture fractures for every 1-centimetre aperture fracture, and mapped tens of 2 cm aperture fractures and even larger dilatant jogs. A wide-spread field observation that also applies at Yucca Mountain is that longer fractures tend to be better interconnected than shorter ones (e.g., Bonnet et al. [72]). As a consequence only a small fraction of the overall fracture population plays an active role for fluid flow (Seol et al. [73]). These systematics indicate that the flow properties of fractured rocks in the vadose zone are highly heterogeneous, promoting the development of highly localised flow paths. From the fracture aperture, a, an upper bound on the discharge through the fracture conduit, q(m 3 /s), can be calculated using the parallel plate law (e.g., Snow, 1965), see Zimmermann and Bodvarsson (1996). Here, ρ and µ are the fluid density (kg/m 3 ) and dynamic viscosity (Pa · s), and g is the acceleration of gravity (m/s 2 ) acting in the vertical direction I. We used this relationship to calculate the gravity-driven advance of the wetting front in vertical fractures as a function of fracture aperture, for verification of our DFM simulator for the case of an impermeable solid rock matrix.
Regarding fracture-matrix transfer in the case of a permeable porous rock matrix, the cubic increase of discharge with aperture suggests that there should be a threshold aperture above which, the flow rate in the fracture should no longer be influenced by the transfer. This important threshold cannot be determined analytically because it depends on the nonlinear saturation-capillary pressure relationships of the matrix and its initial water saturation. We will evaluate this further by numeric simulation assuming that co-and counter-current imbibition can be approximated by Van Genuchten saturation functions.
Apart from computing frontal advance in individual fractures, we investigate how multiple fractures with the distribution of apertures given in Figure 4, drain a shared fluid reservoir at their top. This reservoir might be a water saturated coarse sand or gravel overlying the fractured bedrock. The objective of the multi-fracture calculation is to get an idea of the maximum infiltration rate that the fractured bedrock can sustain. This rate will also be the minimum rate required to prevent film flow in the larger fractures.

Single-Fracture Physical Experiment for Model Validation (Model CUBE)
An immiscible displacement experiment conducted by Rangel-German and Kovscek [25] on a homogeneous sandstone block was selected for the validation of the new fracturematrix discretisation and nonlinear transfer scheme because (1) its geometry and setup are simple and fully documented, (2) the properties of the rock matrix and the fluids are given, and (3) there are detailed measurements that we can compare our results with.
The model CUBE is shown in Figure 5. It consists of a 5 × 5 × 5 cm sample of highly permeable Berea sandstone. Except for the bottom boundary that was used to model the fracture, the cube was coated on all sides with epoxy resin. A small gap (0.1 mm) between the block and the impermeable core holder was used to model the constant aperture fracture. The rock sample was dried in a vacuum oven. Then, water was continuously injected through the injection port on the bottom left of the model at a constant flow rate of 0.1 cc/min. Excess fluid was collected from the production port in the model's bottom right corner. To simulate the imbibition experiment, a mixed-dimensional triangular FE mesh was created (Figure 5b). Additional line elements along the bottom side of the block represent the fracture.
Since Rangel-German selected a very small fracture aperture, its value is affected by the grains and pores of the free surface of the sandstone sample. To capture related aperture variations, the Brooks-Corey model was chosen to parameterise the fracture, albeit with a lower lambda-parameter ( Table 1). As in the physical experiment, a constant-rate water source was applied on the fracture's left end, and constant atmospheric pressure on its right end. All other boundaries were treated as natural (no-flow boundary condition).

Fractured Volcanic Tuff (Model TUFF2D)
An idealised two-dimensional DFM model with vertical fractures intersecting a porous volcanic tuff with largely uniform micro-pores and a milli-Darcy permeability (cf., Dinçer and Bostancı [74]) serves as the second test case for the new DFM embedded-discontinuity solver ( Figure 6, Table 2). A 20 cm thick layer of highly permeable and porous sandy soil covers this rock (0.4 −11 m 2 , φ = 0.28). The model dimensions are 4.0 × 3.5 m. The fractured tuff beneath the soil layer accounts for the largest portion of the pore volume of the model. Its permeability is 0.5·10 −15 m 2 and entry pressure p d = 24.4 kPa. Parallel fractures penetrate the rock layer. Their hydraulic aperture alternates between 0.5 and 1.1 mm. The parallel-plate law determines fracture permeability: k f = b 2 / 12, results in permeability of k f 1 = 2.083 · 10 −8 m 2 and k f 2 = 1.008 · 10 −7 m 2 respectively, see Zimmerman and Bodvarsson [75]. A triangular mesh is used to discretise the domain with 3294 elements and 2002 nodes. The split operator introduces 231 new line elements representing the vertical fractures.
The van Genuchten relative permeability model and capillary pressure curves are used for the air phase in both fracture and matrix (Figures 7 and 8).
The model setup consists of a constant pressure bottom boundary condition modelling a water table. The left and right sides of the model are treated as impermeable boundaries.
Two different precipitation rates are explored to study their effect on water movement in the model. In the first scenario (Case A), a surface source is applied at the top to simulate a small rainfall event at an average precipitation rate of 10.5 mm/h. In the second scenario (Case B), the precipitation is much higher, at 73 mm/h, to simulate a vigorous rainfall event. We assume that all the rainwater imbibes into the flat model top. Thus, runoff is ignored as it will be in the scope of another study.
In each rainfall scenario, we investigate the effect of residual saturation on the response of the system. For this purpose, we run two simulations corresponding to different residual water saturation for the matrix and the soil layer. The first run uses the same irreducible water saturation of 0.08 for the soil and the rock layers. The values for the soil and the rock layers in the second run are 0.15 and 0.45, respectively.

Three-Dimensional Model of Fractured Volcanic Tuff (Model TUFF3D)
To study the effect of fracture aperture variations on the saturation front advancing in the fracture plane, a three-dimensional model was constructed (Figure 9a, Table 3). The dimensions of this model are 0.5 × 0.5 × 2 m. Similar to the previous cases, a cover layer of soil stays at the top model, and a fracture goes through a lower permeable rock layer. The finite element mesh divides the matrix into 249,750 tetra elements and the fracture surface into 17,258 triangle elements (Figure 9b).  We employ a Perlin noise generator to create a spatially correlated random variable field used as a possible realisation of fracture aperture. Then other features of the fracture are assessed by the local thickness, stored as the cell's attributes. As a result, fracture properties spread smoothly from cell to cell. The average aperture generated is 0.5 mm, and the ratio between the maximum and minimum aperture is ten folds (Figure 10). Fluid properties are the same as in the previous test case, and the media properties are shown in Table 2. Unlike the previous models, we use mixed relative permeability models for this 3D model. For the matrix and the soil layer, the van Genuchten model describes the relative permeability of the aqueous phase and capillary pressure, while the Brooks -Corey model is used to describe the relative permeability of the aerial phase. And the van Genuchten's model describes both the relative permeability and the capillary pressure for the fracture (Figure 11). To setup the model, a 125 mm/(m 2 · day) water volume source is applied at the top of the model, and the lateral boundaries are treated as impermeable. The bottom boundary is assigned a constant pressure.
A second constant aperture version of this model (a = 0.5 mm) was created to facilitate a comparison with the variable aperture one.   Figure 11. Relative permeability curves assigned to rock matrix and fracture in model TUFF3D.

Results
We start with the verification of the new fracture discretisation and nonlinear FMT algorithm with the analytic calculations introduced earlier. We then proceed with its validation with the physical experiment. Subsequently, the results from the two-dimensional fractured-tuff models are presented, investigating the role of capillary forces on the infiltration of rainwater into unsaturated rocks. The impacts of infiltration rate and initial water saturation of the rock matrix are quantified. Finally the three-dimensional single-fracture model is presented to analyse the impact of fracture-aperture variations on the shape of the wetting front in the fracture plane.

Hand Calculations of Infiltration Front Velocity
The parallel plate law (Equation (13)) gives an average wetting front velocity of 0.01 m/s in a fracture with a uniform aperture of 0.1 mm, 0.2 m/s for an aperture of 0.5 mm, and 90 m/s for a 1-cm aperture fracture. All these values are for gravity-driven laminar flow of water displacing air, and ignoring capillary effects. To create an approximate multi-fracture model of the Tiva Canyon Tuff, the fracture aperture-frequency distribution (Figure 4), fitting the powerlaw, f = 37.2824 a −1.5715 (14) is In the aperture-statistics based multi-fracture model, wider, ≥3-mm-aperture fractures determine the transport capacity in spite of their relatively rare occurrence. The maximum rate at which water can infiltrate is 0.0016 m 3 /s. This means that the amount of water needed to saturate the fracture network was available at the historic peak precipitation rate irrespective of fracture intensity. An excess volume of fluid would have been available where surface runoff collects into depressions, canyons or washes. These backof-the-envelope calculations further suggest that due to the very small fracture porosity of the modelled system (≤10 −6 ) and amplified by preferential flow in the larger fractures, drainage of the fully water saturated soil layer downward into the fractures would support the ≥100-m-deep infiltration of ≥1-cm aperture fractures. For these calculations to be applicable, flowing fractures would need to become fully saturated, flow would have to remain laminar, and imbibition of the rock matrix should not occur. Inducing peak fracture flow velocities of metres per second, a multi-hour rainfall event might even provoke penetration down to the base of the vadose zone. Yet this is just a simplistic back-of-the-envelope calculation for idealised vertical fractures. In the following we investigate this with the more realistic albeit still highly idealised discrete fracture and matrix model TUFF2D, which contains only fractures with narrow apertures of 0.5 and 1.1-mm so that laminar creeping flow is a reasonable assumption.

Validation of Discrete Fracture and Matrix Model
X-ray tomography images of saturation in the physical experiment are shown alongside with our simulation results in Figure 12. An excellent match is observed between the position of the wetting front as a function of time. Since the matrix is completely dry, spontaneous imbibition is strong and most of the injected water imbibes the matrix instead of filling the fracture ( Figure 13). Subsequently, the wetting front in the matrix moves at a similar speed as that in the fracture. Rangel-German and Kovscek [25] termed this behaviour 'filling' regime, to contrast it from that which typified faster injection rates where the wetting front in the fracture overtakes that in the rock matrix. Figure 14 quantifies matrix imbibition over time. The average matrix saturation curves have a distinct 'S' shape. Their first part represents an early time response when water enters the fracture. The middle part describes the spontaneous imbibition process during which FMT scales linearly with the square root of time. The plateau of the final part illustrates 'late-time' behaviour that begins when the wetting front reaches the sealed sample boundaries, and the rest of the air is trapped in the rock sample. The late-time response is characterised by a drop in FMT and nearly constant saturation gradients in the rock matrix. The outlet port produces mostly water.
Our simulation captures these three evolutionary stages albeit with a small but constant offset in time. In the simulation, FMT ramps up earlier than in the physical experiment. There may be several reasons for this mismatch: Candidates are the constitutive relationships used for the fracture (the experimental study did not establish these and we just used Brooks-Corey functions), the geometry of the inlet fluid port that was not captured exactly or perhaps a minor delay between when the pump was switched on and when the fluid pulse arrived in the sample. A further possibility relates to the regression of the experimental results (cf., Schmid and Geiger [39]). However, the similarities in saturation patterns, average water saturation over time, and evolutionary stages, strongly suggest that the proposed fracture discretisation and nonlinear FMT scheme can realistically capture the prominent physics at play in the experiment. We therefore suggest that these results constitute a model validation.

Infiltration Simulations
For a fully water-saturated rock matrix, the velocity in stand-alone vertical fractures only depends on aperture and fluid availability. The maximum velocity observed in the simulations conducted with model TUFF2D in its 1.1 and 0.5-mm aperture fractures at the low rainfall intensity are 0.006 and 0.0029 m/s, respectively. At the high intensity, these velocities rise to 0.041 m/s and 0.020 m/s, but are still almost an order of magnitude lower than those predicted by the back-of-the-envelope calculations. Apart from the consideration of FMT, this difference arises because the numeric model also considers the multiphase fluid flow in the soil layer at its top while the back-of-the-envelope model assumes an unlimited fluid supply rate to each fracture.
The velocity field of the fully water-saturated model is shown in Figure 15. The 1.1-mm aperture fractures become the preferential flow paths regardless of rainfall rate. Snapshots of an initially dry version of model TUFF2D, after three hours of lowintensity rainfall are shown in Figure 16a. This saturation map shows saturation jump discontinuities at the interface between the soil and the bedrock as well as between the fractures and the rock matrix. Water absorbs into 76% of the domain area. While the singlephase velocity analysis suggests that the wetting front might advance very quickly, at least in the cm-aperture fractures, the results for the initially dry system and the moderate rainfall intensity show a much reduced wetting front velocity. This is the effect of spontaneous imbibition of water into the rock matrix that dominates the behaviour of the system. The infiltrating water imbibes creating a large and smooth halo surrounding each fracture and there are hardly any aperture-related differences. While fracture water imbibes the matrix, matrix air drains out in a counter-current flow, rises upward, and discharges into the soil layer where the fracture terminates. Due to the limited rainwater supply, water saturation is higher in the 0.5-mm-wide fractures than in the mm-wide ones, which remain essentially air saturated. This difference is due also to the lower capillary pressure within the mm-wide fractures. Within each fracture, water and air move by counter-current flow, while the dry rock matrix rapidly absorbs most of the water expelling its air into the fracture. This composite behaviour explains the relative slowdown of the wetting front in the mm-wide fractures. For the high rainfall intensity (Case B), the front moves much faster through the initially dry model TUFF2D, reaching the same infiltration depth as in the low-intensity Case A in just 10 min Figure 17a. Now, the water penetrates the mm-wide fractures multiple times faster than the narrower ones. The wetted fracture matrix halos are correspondingly narrower spanning only about 35% of the total model area. At this early time in the evolution, less air has a chance to leak out into the soil layer through the fracture system.
(a) (b) Figure 17. Air distribution in model TUFF2D after ten minutes corresponding to a high precipitation rate scenario (case B). (a) low initial water saturation, (b) high initial water saturation. Figure 18 displays co-and counter-current flows in detail, revealing that the flow pattern is essentially the same for the 1.1 and 0.5-mm wide fractures: air and water move in the same direction ahead of the wetting front. Behind the wetting front counter-current flow of water and air occurs.

Effect of the Initial Saturation of the Rock Matrix
Initial water saturation has a pronounced effect: high initial matrix saturation implies that fracture-matrix transfer is reduced, accelerating fracture infiltration. This can be seen in Figures 16b and 17b. For both rainfall intensities, the partial saturation of the rock matrix results in a faster downward penetration of the wetting front in the fractures as compared with the dry scenario.

Simulation of Wetting Front Retardation by Fracture-Matrix Transfer
All results obtained from model TUFF2D show wetting-front velocities lower than those estimated by the hand calculations. The velocity differences mark the significant retardation of the front by fracture-matrix transfer. The front velocity diagrams Figure 19 reveals two distinct behaviours. Early time behaviour is marked by steep slopes, indicating rapid penetration of water into the fracture system. Subsequent front advance in the fracture system is up to three times slower. A similar behaviour is seen in the experiments conducted by Rangel-German and Kovscek [25], and referred to by these authors as earlyand late-time response. For low fluid supply rates, capillary diffusion dominates the system, leading to extensive fracture-matrix transfer and front retardation. For this scenario, the aperture-related difference in wetting front velocity in model TUFF2D is small (Figure 19a, Table 4). Thus, during late-time behaviour, the delayed velocity of the wetting front is only one-tenth of the value obtained by the hand calculation.
By contrast, in the high rainfall intensity scenario, gravity-driven fracture flow dominates the response of the system. Capillary-driven fracture-matrix transfer is relatively small. Consequently, front retardation is small and front velocity reaches about one third of the hand-calculation value (Figure 19b, Table 5). Also, the curves for the different fractures diverge more.
As the last case that was considered for model TUFF2D, an elevated initial saturation of the rock matrix eliminates front retardation.  The reference depth 0.0m-is that of the soil-bedrock interface. Note that the front moves much faster in the model at a matrix water saturation of 0.45 than at a hypothetical residual water saturation (moisture retention capacity) taken as 0.08. Clearly, a two-dimensional model cannot resolve processes that operate in the plane of fractures, like the counter-current flow of water and air which was nevertheless indicated by the results obtained from model TUFF2D. Moreover, spatial aperture variations in the fracture cannot be represented.
The three-dimensional constant aperture model TUFF3D displays a behaviour largely similar to TUFF2D ( Figure 20): water spontaneously imbibes into the soil layer before infiltrating further downward through the fractures. The wetting front in the constant aperture fracture is sharp and stable, resulting in a kind of sheet flow.
By contrast, in the variable aperture model TUFF3D-VA, an unstable saturation front evolves in the fracture (Figure 20). Water preferentially flows along the constrictions, trapping air pockets between these flow zones. Further downward, the water eventually also occupies fracture segments with a wider aperture, developing lobes that infiltrate faster than the average velocity of the wetting front. In constant-aperture TUFF3D, the average wetting front velocity is 1.75 mm/s, but in the variable aperture model TUFF3D-VA it reaches 1.93 mm/s. Co-and counter-current flow patterns in the rock matrix are similar to the ones seen in the low rainfall intensity scenario of the two-dimensional model, (Figure 21).

Discussion
The simulation results confirm that aperture is of central importance for the transport velocity in fractures embedded in a permeable rock matrix (Duguid and Lee [79], Guzman and Aziz [80], Bogdanov et al. [81]). Our work highlights that the velocity of the wetting front during gravity-driven flow in small-aperture fractures strongly depends on the magnitude of capillary-driven fracture-matrix transfer. Our study corroborates Cey et al. [82] finding that the matrix hydraulic conductivity and initial moisture content are most influential for this exchange. However, this control becomes less and less important with increasing fracture aperture. Above a certain threshold aperture flow velocity matches the analytic solution for gravity-dominated flow in an impermeable matrix.
While only narrow ≤1-mm apertures were investigated herein, we observed that for the heavy to violent rain that typifies thunderstorms, gravitational forces began to dominate over viscous in the one-mm aperture fractures underneath the saturated sandy soil layer. This seems to indicate that fracture-matrix transfer will have little effect on the wetting front velocity in cm-aperture fractures. This is the scope of ongoing work on the larger scale which will be the subject of a further publication.
Wang and Narasimhan [21] argued that the rock matrix needs to be fully water saturated before there is fluid flow in the fractures. This does not apply for larger fractures. Earlier on, Montazer and Wilson [24] had already proposed that it only takes a wet fracture wall for water to infiltrate a fracture. The experimental work on film flow by Tokunaga and Wan [68] and the numeric simulations presented herein confirm this suggestion: While water spontaneously imbibes the rock surrounding the fracture, this infiltration is restricted to a relatively small fracture selvage and once this wet zone is established, the wetting front in the fracture advances freely. Consistent with the capillary interface model, spontaneous imbibition occurs rapidly directly adjacent to the fracture because of the sharp initial capillary pressure gradient, but the exchange rate quickly drops asympotically as the selvage widens. Transfer is marginal when the rock matrix already is at a saturation that reflects the plateau of the capillary pressure curve at the time when fracture infiltration occurs, and it ceases completely as the matrix approaches saturation.
Aperture variations induce unique flow patterns within the fracture(s). If they are slight (<10% of the overall aperture), sheet flow remains prominent, in agreement with the conclusion of Cey et al. [82]. When they are larger, fingering destabilises the wetting front and water flow in the individual fingers is accelerated (Glass [29], Kawamoto and Miyazaki [83], Nicholl and Glass [84], Cueto-Felgueroso et al. [85]).

Conclusions
Rainwater infiltration into open fractures in porous rock is investigated by means of discrete fracture and matrix (DFM) simulation, using a hybrid finite element-finite volume method with lower dimensional fracture representations. Fracture nodes with pressure and saturation degrees of freedom are separated from the rock matrix by "split boundaries." Complex fracture-matrix interface conditions and fluid transfer are modelled using a specific nonlinear solver embedded into an operator-splitting IMPES approach. This permits simulation of spontaneous imbibition of rainwater into the rock matrix as well as capillary sealing. The new scheme also resolves the movement of water and air in the fractures, which is demonstrated with two-and three-dimensional models with variable fracture aperture. As expected, the results indicate that spontaneous imbibition retards the wetting front significantly. However, this effect depends on rainfall intensity, fracture aperture, and rock matrix water saturation. Matrix imbibition retards the wetting front in fractures only up to a critical fracture aperture, above which the wetting front velocity matches analytic solutions for saturated fracture flow. In all cases, front velocity is constant unless the fluid supply runs out. Multiple interacting fractures have yet to be investigated. Regarding rainfall events in semi arid regions, vertically continuous fractures are shown to rapidly conduct water far below the earth's surface, suggesting that the saturated-zone aquifers below can get recharged by these events. Fracture aperture variation within individual fractures destabilises wetting fronts. Thus, even without consideration of complex wetting front physics like pinning, sheet flow can give way to fingering as narrow, water-saturated down-flow zones alternate with air-filled up-flow regions within the fracture.