Numerical and Analytical Modeling of Flow Partitioning in Partially Saturated Fracture Networks

Infiltration processes in fractured‐porous media remain a crucial, yet not very well understood component of recharge and vulnerability assessment. Under partially saturated condition flows in fractures, percolating fracture networks and fault zones contribute to the fastest spectrum of infiltration velocities via preferential pathways. Specifically, the partitioning dynamics at fracture intersections determine the magnitude of flow fragmentation into vertical and horizontal components, hence the bulk flow velocity and dispersion of fracture networks. Here we derive an approximate analytical solution for the partitioning process and validate it using smoothed particle hydrodynamics simulations. The transfer function is conceptually based on simulation results and laboratory experiments carried out in previous works. It allows efficient flow simulation through fracture networks with simple cubic structures and an arbitrary number of fractures and aperture sizes via linear response theory and convolution of a given input signal. We derive a nondimensional bulk flow velocity ( v∼ ) and dispersion coefficient ( D~ ) to characterize fracture networks in terms of dimensionless horizontal and vertical time scales τm and τ0. The dispersion coefficient strongly depends on the horizontal time scale and converges toward a constant value of 0.08 within reasonable fluid and geometrical parameter ranges, while the nondimensional velocity exhibits a characteristic v~∼τm−1/2 scaling. Given that hydraulic information is often only available at limited places within (fractured‐porous) aquifer systems (boreholes or springs), our study intends to provide an analytical concept to potentially reconstruct internal fracture network geometries from external boundary information, such as the dispersive properties of discharge (groundwater level fluctuations).

therefore dispersion of an input signal. The complexity arises from geological heterogeneities that provide continuous pathways on various scales for rapid percolation and transport within fractures. In karst systems, precipitation is commonly partitioned into diffuse and preferential components, where the latter is commonly linked to direct infiltration in the surrounding area of surface depressions, dry valleys, and dolines (Gunn, 1981;Kordilla et al., 2012;Sauter, 1992;Williams, 2008). Fault zones may cut across several geological units and provide catchment-scale preferential flow paths in the form of strongly connected clusters of fractures (Bodvarsson et al., 1997;Flint et al., 2001;H. H. Liu et al., 2004). Tectonically induced stress fields and stress field changes generally promote the formation of local discontinuities, such as fractures, joints, and fault zones in consolidated porous rocks (Ford & Williams, 2013;Neslon, 2001). What sets such features apart from typical pore space geometries is their strong anisotropic character, that is, their length or spatial extent is orders of magnitude larger than their aperture. When fractures are connected, they can form percolating clusters (Adler et al., 2013;Berkowitz & Scher, 1995) that can reach length scales far beyond the thickness of individual geological layers/units and potentially extend across the entire vadose zone. Similar features can be observed in soil systems, a type of material where the heterogeneity is commonly associated with macropores (wormholes), which can also form percolating clusters (Hussain et al., 2019;Jarvis, 1998;Nimmo, 2010).
Assessing recharge dynamics in fractured-porous systems on the field scale is difficult (Scanlon & Cook, 2002). Phreatic zone techniques assess recharge at the water table or at springs (e.g., tracers, water table fluctuations; Cook & Solomon, 1997;Nimmo & Perkins, 2018). Therefore, the estimates can potentially reflect catchment-scale dynamics or at least subcatchment recharge processes within the hydraulic influence area of the measurement point. In contrast, vadose zone techniques rely on measurements above the groundwater table (e.g., lysimeters, Darcy's law, tracers, Chambers et al., 2019;Heppner et al., 2007;Rossman et al., 2014). They allow a rather localized quantification of recharge or water content and only integrate a limited volume above the point of measurement as infiltration commonly occurs nearly vertical. Because most of these methods rely on simple assumptions about internal systems geometry and percolation processes, the predictive power and temporal resolution is often limited (Scanlon & Cook, 2002).
To shed light on complex infiltration processes, laboratory-scale experiments have been a promising addition to the former investigation methods because they allow important processes to be isolated under well-controlled conditions that are often impossible to observe in situ. Small-scale laboratory experiments for gravity-driven partially saturated flow often exhibit erratic or chaotic flow dynamics (Dragila & Weisbrod, 2004;Su et al., 2001;T. Wood & Huang, 2015). In general, flow modes on the wall of wide fractures evolve with increasing flow rates from thin adsorbed films to droplets and rivulets to wavy surface films (Dippenaar & Van Rooy, 2016;Dragila & Wheatcraft, 2001;Ghezzehei, 2004;Jones et al., 2017). Different flow modes may also coexist. Consequently, experimental results are difficult to cast into meaningful frameworks. This especially concerns the complex flow dynamics at fracture intersections, which act as critical relay points controlling: (1) the overall connectivity of fracture networks (Adler et al., 2013); (2) the flow partitioning dynamics between connected fracture elements (Dragila & Weisbrod, 2004;Xue et al., 2020;Yang et al., 2019); and, ultimately, (3) the distribution of flow modes on fracture surfaces (Dippenaar & Van Rooy, 2016;Jones et al., 2017;Shigorina et al., 2019), which, in turn, can affect the interaction between the porous matrix and fracture (Tokunaga, 2009;Tokunaga & Wan, 1997). Here, the term "partitioning" refers to the process of fluid redistribution at a fracture intersection, which depends on the relation between capillary, inertial, and viscous forces  and complexities such as velocity-dependent contact angles (Xue et al., 2020;Yang et al., 2019).
In terms of fracture aperture, numerical and laboratory studies of unsaturated flow in fractures have covered various length scales, from submillimeter scales Ji et al, 2004Ji et al, , 2006 over ranges close to the capillary-inertial transition around 0.7 mm (T. R. Wood et al., 2002Wood et al., , 2005 to apertures well within the inertial-dominated regime (Dragila & Weisbrod, 2004;Tartakovsky & Meakin, 2005a, 2005bTokunaga & Wan, 1997Huang et al., 2005;M. Liu et al., 2007). Studies of free surface flow on a fracture plane without an intersection have been conducted by Shigorina et al. ( Depending on the experimental setup, studies have focused either on the partitioning process at fracture intersections (Dragila & Weisbrod, 2004) or the (long-term) bulk system response (Ebel & Nimmo, 2013;Nimmo, 2010), both of which are integral parts of understanding preferential flow dynamics through fractured systems. For the former case, a single fracture intersection of a horizontal and a vertical fracture is often "constructed," for example, by breaking glass plates, which results in a quasi-two-dimensional setup (Ji et al., 2006). Modifications to this setup include experiments with a slight offset at the fracture intersection or T-shaped intersections at various degrees of rotation (T. R. Wood et al., 2005;Xue et al., 2020;Yang et al., 2019). Intersections resembling an inverted Y-structure have been studied by Dragila and Weisbrod (2004), M. Liu et al. (2007) and Tartakovsky and Meakin (2005a). The combination of several commonly crossshaped fracture intersections allows flow convergence to be studied, that is, the deviation from classical volume-effective diffusive flow dynamics that are typical for nonfracture porous media. Studies of this kind have been conducted by T. R. Wood et al. (2002Wood et al. ( , 2005; T. Wood and Huang (2015); Glass et al. (2003); LaViolette et al. (2003), often reaching timescales of several minutes or days. In the study by Glass et al. (2003), fractures are embedded into an impermeable matrix, while the other authors constructed their fracture networks from geological materials. Incorporation of the porous matrix can be considered another important classification parameter of fracture-scale studies. For large-aperture fractures, that is, inertial-dominated flow systems that limit the contact time between fracture flow, matrix flow, and/or a low-permeable matrix, the effect of matrix storage may be neglected. This can be observed in fractured karst systems, where fractures are often enlarged by dissolution (Benson, 2001;Dijk et al., 2002) or fractured crystalline rocks with extremely low matrix porosity and a severely limited advective potential.
Despite these research efforts, the gap between small-scale process understanding and larger-scale application is still significant. In our recent work (Noffz et al., 2019), we demonstrated how to model breakthrough behavior in terms of discharge at the bottom of arbitrary long stacks of sugar-cube fracture arrays (Barenblatt et al., 1960) via linear response theory and convolution of input signals, whereas the transfer function has been obtained empirically for a given setup of a wide-aperture vertical surface intersected by a horizontal one. However, it is desirable to obtain the form of the transfer function a priori using information about the internal geometry as well as fluid properties and fluid-solid interaction characteristics. Therefore, in this work, we provide an analytical solution for the transfer function and validate it using numerical simulations. The analytical solution describes the horizontal fracture infiltration until critical pressure thresholds trigger the breakthrough. The following dynamics are governed by Washburn-type flows and are conceptually based on the numerical studies and former laboratory studies (Noffz et al., 2019). Vertical flows are approximated by a film flow model. Finally, we employ linear response theory to model flow through arbitrary numbers of fracture intersections with explicit geometry and derive nondimensional dispersion and velocity parameters (   , D v) that depend on the dimensionless horizontal and vertical fracture time scales (τ m , τ 0 ). Flows are shown to converge to a near-constant dispersion coefficient with increasing τ m , while nondimensional velocities scale as within feasible critical Reynolds number ranges.

SPH Model
We use a two-dimensional SPH model to analyze complex flow partitioning at fracture intersections. SPH is a Lagrangian meshless method able to simulate complex flows with highly dynamic interfaces and is especially suited for the simulation of free-surface (pseudo-multiphase) liquid flows with a continuous gas phase, effects of surface tension, and static/dynamic contact angles. We use a two-dimensional version of the massively parallel three-dimensional code of Kordilla et al. (2017) that has been extended with an alternative formulation of the no-slip boundary condition. For a detailed description of the SPH free flow model and its implementation in a parallel code, the reader is referred to Kordilla et al. (2017) and references therein. The SPH equations are summarized in Appendix A.
Here, we validate the SPH code for two classical static and dynamic flow cases that are related to the processes encountered in our application of flow in fractures, Poiseuille flow in a parallel plate system, and capillary rise in a vertical tube. The Poiseuille flow example serves to verify the proper implementation of viscous forces at the solid-fluid interface. These forces are important for both the film flow, where flow is bounded on one side only, and the flow that fully fills horizontal fractures.

Poiseuille Flow
In this section, we demonstrate that for sufficiently large values of the friction coefficient β in the SPH momentum conservation Equations A3 and A6, the SPH method recovers the solution of the NS equations subject to no-slip boundary conditions at the fluid-solid boundary. Specifically, we use the SPH code with β ranging from 1 × 10 −1 to 1 × 10 2 kgm 2 s −1 to simulate a two-dimensional Poiseuille flow problem and validate the SPH solutions for velocity against the analytical solution for the no-slip boundary condition (Sigalotti et al., 2003) where the center is located at y = 0, d = L/2 such that the solid boundaries are located at y = ±d, ν = μ/ρ is the kinematic viscosity, and a f is the component of the body force per unit mass acting in the x direction (the direction of flow).
Flow is simulated using the following parameter set: the interparticle spacing is Δx = 5 × 10 −5 m, L = 200Δx = 1 × 10 −2 m, ρ = 1,000 kgm −3 , and μ = 1.25 × 10 −3 kgm −1 s −1 , and an acceleration of a f = 1.25 × 10 −5 ms −2 is applied parallel to the x direction. Five layers of boundary particles are placed at y = 0 and y = L to assure kernel consistency. For the given parameter set, this yields a Reynolds number of Results indicate that the SPH solution converges to the exact no-slip solution for β > 10 (which corresponds to the artificial slip length λ < L/100) with an error on the order of 1.5% or lower. This holds for all time steps during the initial acceleration of the fluid within the capillary. It should be noted that for much higher values of β, the time steps have to become increasingly smaller according to Equation A7, which was not enforced in these simulations. Therefore, a slight increase in the error can be observed for values above β = 20, a strategy that is employed in this study.

Capillary Rise in a Tube
Here we simulate capillary rise in tubes of varying radii and compare the equilibrium fluid column height to the classical theory of Jurin (1718) and the extended theories of Legait and de Gennes (1984) and Barozzi and Angeli (2014).
The classical theory of capillary rise is based on the parallel plate concept: where r is the radius of the fracture, and h is the height of the triple contact line from the water surface. The total pressure in a two-dimensional system consists of the capillary pressure and the pressure due to the weight of the water column Plugging the total pressure   2 Δ Δ D c h P P P into Equation 2 and for dh/dt = 0, the maximum rise becomes: Because the curvature of the meniscus slightly depends on h, a common extension of Equation 4 is given by Legait & de Gennes (1984) (5) Barozzi and Angeli (2014) extend the solution by adding a correction term that accounts for the additional fluid volume over the apex of the meniscus The SPH simulations are run with an interparticle spacing of Δx = 5 × 10 −5 m, and a density of ρ = 1,000 kgm −3 and a body force of g = 9.81 ms −2 are applied in normal direction to the bottom boundary. The viscosity is μ = 0.001 kgm −1 s −1 and the no-slip condition is enforced with β = 25. The speed of sound is set to c 0 = 3 m.s −1 . Interaction forces are set to s sf = 0.015 and s ff = 0.02, which yields a surface tension of σ = 0.0742 kgs −2 and a static contact angle of θ 0 = 69°. The domain has a width of L x = 800Δx = 4.0 cm. The height of the capillary is L y = 340Δx = 1.7 cm and is placed ΔL y = 60Δx = 3 mm above the bottom boundary. Mirror boundaries are applied in the x-direction. All solid boundaries are five particles thick to assure kernel consistency. Simulations are initiated with a flat fluid surface covering the domain with an initial height of 145Δx = 7.25 mm. The aperture of the capillary is varied in a range of 1.5-3.5 mm.
Simulations are run until an equilibrium is established and the maximum height is reached within the capillary. In order to measure Δh, we determine the minimum height h min of the fluid as the average of the water height 20Δx away from the left and right mirror boundaries (see Figure 1). The maximum height h max of the fluid column is measured at the outer part of the capillary meniscus, and we therefore obtain Δh = h max − h min . The contact angles at equilibrium are obtained from a circle fit using the Pratt method (Pratt, 1987) (Figure 2, left).
Results of the SPH simulations and theoretical results are shown in Figure 2. Numerical results are in good agreement with the theoretical predictions and lie in between the predictions of Jurin (1718), Legait and de Gennes (1984) and Barozzi and Angeli (2014). The latter study takes into account the effect of additional fluid volume above the meniscus apex which becomes increasingly important when the capillary rise h c is on the order of capillary radius r. This is true for the shown example (h c /2r ∼ 0.5-2.2).

Results and Discussion
In the following subsections we (1) conceptualize the flow partitioning at a T-type fracture intersection, (2) derive an analytical transfer function for the partitioning including fluid movement on the vertical surfaces, (3) provide an upscaling solution via convolution and linear response theory, and finally (4) derive expressions for arrival times and bulk dispersion that are then (5) analyzed in nondimensional form to provide a comprehensive picture of the larger-scale infiltration dynamics and its relation to the internal geometry.

Rivulet Flow Partitioning at a Fracture Intersection
In this section, we derive a solution for the partitioning dynamics of rivulet flow down a vertical plane intersected by a horizontal smooth fracture (see Figure 3) and compare it to our SPH model results. We consider the threshold at which critical capillary pressures within the horizontal fracture are high enough to route flow further down onto the vertical surface. At this point, flow in the horizontal fracture transitions from a linear plug-flow type into a Washburn-type flow regime. It should be noted that the conceptual model assumes that flow within the horizontal fracture is capillary-driven; therefore, wide-aperture horizontal fractures (e.g., cave-like structures) where the fluid may establish free-surface flows (low capillary flows) are not taken into account. Furthermore, we assume that flow on the vertical surfaces takes the shape of rivulets. This assumption is often made in soil systems for preferential flows (Bogner & Germann, 2019;Germann & Hensel, 2006;Nimmo, 2010) and has been shown to also hold under laboratory conditions (Noffz et al., 2019), yet, more complex flow-rate-dependent modes (droplet, slugs) may occur (Ghezzehei, 2004;Kordilla et al., 2017) that trigger early partitioning at intersections (T. R. Wood et al., 2005;Xue et al., 2020;Yang et al., 2019) and are different from the sequential approach considered in this work.
The flow rate Q h (t) (m 2 s −1 ) in the horizontal fracture is approximately given by the Darcy law: where t is time from the moment water entered the horizontal fracture, k = a 2 /12 (m 2 ), a is the aperture (m), μ is the viscosity, P f and P in (t) are the pressures at the invading front (point 1 in Figure 3a) and the horizontal KORDILLA ET AL.
10.1029/2020WR028775 6 of 27 fracture entrance (point 2 in Figure 3a), respectively, and l(t) is the distance from the front to the fracture entrance.
From the Young-Laplace law, the pressure at the invading front is where σ is the surface water-air surface tension, P air is the air pressure,   2 cos a R is the front curvature, and θ is the contact angle.
Initially, all flow in the vertical fracture is diverted to (imbibed into) the horizontal fracture, that is, Q h (t) = Q 0 , as shown in Figure 3a. Later, flow partitions, where flow both penetrates the horizontal fracture and flows down the wall of the vertical fracture segment, as depicted in Figure 3b. When flow is partitioned, P in (t) = P air (point 3 in Figure 3b). In the following analysis, we assume that partitioning occurs instantaneously at time t = t c . Then, Equation 7 can be rewritten as The front position in the horizontal fracture at the time of partitioning is obtained by setting   Figure 4 shows l c and t c , computed from Equations 11 and 12 and direct SPH simulations, as a function of the flow rate Q in for three horizontal apertures (2.5, 3, and 3.5 mm). SPH simulations for the fracture aperture 2.5 mm and three different flow rates are shown in Figure 5. In SPH simulations, we use an interparticle spacing of Δx = 5 × 10 −5 m, a density of ρ = 1,000 kgm −3 , and a body force of g = 9.81 ms −2 applied normal to the horizontal fracture plane. The surface tension is σ = 0.0742 kgs −2 with the interaction parameters s ff = 0.015 and s sf = 0.0125 and speed of sound c 0 = 3 ms −1 . The viscosity is slightly increased to μ = 0.005 kgm −1 s −1 to limit the required length of the horizontal fracture (and hence computation time), which was set to L = 0.25 m. The no-slip boundary condition is enforced with β = 25 (see Figure 6). For the flow rates between Q 0 = 3 × 10 −5 m 2 s −1 and 8 × 10 −5 m 2 s −1 and fracture apertures between a = 2.5 and 3.5 mm, the critical penetration length can be observed within the chosen fracture length. To avoid erratic partitioning behavior at the fracture intersections (i.e., bypassing droplets), we initiate the simulations with a rivulet on the upper vertical surface, which is already in contact with the horizontal fracture aperture at the start of KORDILLA ET AL.    inated by rivulets. Figure 4 demonstrates that our SPH simulations are in good agreement with the analytical predictions of Equations 11 and 12 for small fracture apertures (2.5 and 3 mm) and larger Q 0 but slightly deviate for larger apertures (3.5 mm) and smaller Q 0 , resulting in the maximum error of ∼12%. We partially attribute this to the fact that the contact angle θ changes during the water penetration into the horizontal fracture (e.g., Popescu et al., 2008). In our analytical model, we disregard dynamic variations in the contact angle and compute θ as an average of the contact angles right after the onset of fracture penetration and close to t c . Yet, Figure 4 demonstrates that our analytical solutions provide an overall good approximation of the partitioning dynamics.

Analytical Solution for the Transfer Function
We now derive analytical solutions for l(t), the front position (or the depth of penetration) in the horizontal fracture, and Q(t), which is the outflow rate below the horizontal fracture junction where fluid is discharged.
Given that the inflow consists of (1) a linear penetration phase and (2) a Washburn-type penetration period, we obtain an analytical solution as follows is the Washburn-type flow velocity after the critical time t c , and t max is the time at which the horizontal fracture is fully saturated and all fluid is channeled further down into the vertical fracture segment. Figure 5 shows the different flow stages observed in the SPH simulation, and Figure 7 shows the time-dependent penetration l(t) and inflow velocity dl(t)/dt obtained from Equation 14 for an aperture a = 2.5 mm and an inflow rate of Q 0 = 4 × 10 −5 m 2 s −1 to Q 0 = 8 × 10 −5 m 2 s −1 . The early time behavior is characterized by a plug-flow regime and hence l(t) ∼ t, whereas after the critical time t c , the inflow scales as  ( ) l t t (gray lines show both scaling regimes). At the time t c , the analytical solution assumes a sudden jump in velocity due to the flow transitions from pure plug-flow within the horizontal fracture to a Washburn-type flow due to the breakthrough at the fracture intersection. The flow velocity in the horizontal fracture drops because of less fluid volume entering the fracture (which is instead channeled downwards onto the vertical fracture). The analytical solution for the penetration velocity dl/dt can describe both regimes (before and after the critical time t c ) and is in very good agreement with the numerical result. A slight deviation can be observed right after the onset of the Washburn behavior at time t c , where channeling into the lower vertical fracture is initiated. Here, a very brief buildup of fluid at the fracture intersection occurs until a critical contact angle is reached and fluid flows downwards, that is, the breakthrough process is not perfectly instantaneous as assumed in the analytical solution but occurs over a small time window close to t c . Figure 5 shows three partitioning types that have been determined based on a qualitative observation. For smaller t c , the breakthrough process is a rather rapid process ("fast breakthrough" right column, Figure 5), while for larger t c (and for example lower flow rates Q 0 ), the buildup of fluid on the vertical surface slightly disperses the breakthrough ( Figure 5, left and middle column, "slow" and "moderately fast" breakthrough), yet a clear sequential progression of plug flow followed by Washburn flow in the horizontal fracture can be observed. The process of fluid buildup is not explicitly considered in our solution and is likely to induce the small temporary drop of the inflow velocity right after t c (beyond the correct drop in velocity due to fluid entering the lower vertical surface instead of the horizontal fracture). However, at later times the velocity correctly converges toward the l(t) ∼ t 0.5 scaling. The cutoff at t max is not shown here because the simulations are stopped when flow reaches the end of the horizontal fracture such that dl/dt = 0. It should be noted that for very small t c or large Q 0 , flow may not exhibit the clear dynamics of sequential partition-KORDILLA ET AL.   14) are in very good agreement. Note that the fluid front did not fully penetrate the horizontal fracture in the simulations; therefore, the cutoff at t max (l = const., dl/dt = 0) is not visible here. Here, a = 2.5 mm and the inflow rate ranges from Q 0 = 4 × 10 −5 m 2 s −1 (bright red) to Q 0 = 8 × 10 −5 (dark red). Gray lines indicate the two characteristic scaling regimes. ing and a breakthrough can occur right away even before the theoretical time t c due to effects of inertia, which we do not consider. Yet, for the covered range of flow rates, our model is in very good agreement with the theoretical solution.
To model the response of the system to a constant input signal Q 0 , we obtain the outflow rate Q leaving the system: Thus, the dimensionless flow rate is Next, we define the normalized transfer function as The transfer function for a plug-flow type regime followed by a Washburn-type behavior has the form where δ is the Dirac delta function.
To numerically integrate the transfer function, we replace the Dirac delta function in Equation 18 by where Δt = 0.1. Figure 8 shows the normalized outflow rate Q/Q 0 (exact and approximate solution) and its derivative, the normalized transfer function    Q dQ dt 0 1 / . The outflow Q/Q 0 is zero at first (all fluid is filling the horizontal fracture) until the critical time t c , where partitioning sets in and inflow is characterized by a Washburn behavior. Finally, when the horizontal fracture is fully saturated at t max , the outflow Q/Q 0 reaches its maximum value, that is, Q = Q 0 and Q/Q 0 = 1.

Extension of the Transfer Function
In the previous section, we focused on the process of horizontal fracture inflow and partitioning; however, we did not consider the effect of additional vertical surfaces above or below the fracture intersection, which affect the system response and therefore the transfer function. In the following section, we extend the transfer function based on classical Nusselt film flow approximations (Nusselt, 1916), which assume a constant film thickness.  (1/Q 0 )dQ/dt (s -1 ) Q/Q 0 (exact) Q/Q 0 (approx.) where y is the direction normal to the surface and α is the inclination angle from the horizontal. The boundary conditions are established via a no-slip condition at y = 0, that is, v x (0) = 0, and the normal viscous stress being zero at the free surface y = h, The solution of this problem is The volumetric flux down the plane is then calculated as and, hence, the maximum film height for a given Q is The depth-averaged velocity can be obtained as The effect of the upper vertical surfaces is simply a delay in the first arrival, i.e., a positive shift in the transfer function by Δ v up t , given that Q = Q 0 . Using Equation 25, we can then simply compute Δ v up t as where v up L is the total length of the upper vertical surface. On the lower vertical surface, a similar shift in the transfer function is induced; however, here the outflow rate is initially Q c = Q(t c ). It should be noted that here t c is the critical time since the beginning of the fracture penetration. For the sake of simplicity, we neglect the increase in Q after the breakthrough at t c and assume that the flow velocity on the lower vertical surface depends on the breakthrough flow rate Q c . The flow rate at the critical breakthrough is obtained via Equations 14 and 15 as where v c = lim ϵ→0 v f (t = t c + ϵ). We then obtain the time Δ v low t as and define the total shift induced by the upper and lower vertical surfaces as The cutoff at time t max , when the horizontal fracture is fully saturated can be computed by setting l(t max ) = l max , which gives The flow rate at the cutoff time t max can be evaluated using Equation 14 with We are now able to compute the full transfer function, including the residence times on the upper and lower vertical surfaces as well as the cutoff at full fracture saturation. In the next chapter, we extend this analysis to model discharge through arbitrary large stacks of fracture intersections via linear response theory.

Analytical Percolation Model for Fracture Cascades
Following Noffz et al. (2019), we employ the transfer function in the context of linear response theory (Jury et al., 1986) to model the outflow rate Q n (t) at the bottom of the vertical surface intersected by n horizontal fractures. For the sake of comparison, two additional simplified cases are considered: (1) Flow within the horizontal fracture is governed entirely by plug-flow even when t max > t c , and (2) horizontal fractures are inactive due to missing capillary action (e.g., macropores, cave-like structures). The considered geometry and its properties with respect to the transformation of an input signal Q 0 serves as a proxy for consecutive routing through further fracture intersections of similar geometry. The outflow rate can be found as a convolution of the input signal Note that for n = 1, the outflow rate Q 1 = Q is given by Equation 15. . Here, the maximum horizontal fracture length is L max = 0.3 m, the aperture a = 2.5 mm, the static contact angle θ 0 = 69.0°, the density ρ = 1,000 kgm −3 , the surface tension σ = 0.0742 kgs −2 , viscosity μ = 0.005 kgm −1 s −1 , and the inflow rate Q 0 = 5 × 10 −5 m 2 s −1 . The upper and lower vertical surfaces have a length of   0.2 v v up low L L m. The two additional cases exhibit no dispersion because of missing Washburn-type dynamics. If the impact of horizontal fracture is completely ignored, the arrival times are nearly one order of magnitude lower than in the complex case. When the penetration is entirely governed by plug-flow, the first moments of the breakthrough correspond to the first moments of the complex case including Washburn-flow. However, the first arrival in the latter case occurs earlier due to the pressure-induced breakthrough. Figure 10 shows the outflow rate for a system with L max = 0.05 m, where flow is dominated by a plug-flow behavior and t max > t c . As expected, the mean breakthrough velocity is higher and the maximum outflow rate Q 0 KORDILLA ET AL.    is reached faster because of the stronger dispersive effect of deeper horizontal fractures. As in the previous example, a faster breakthrough is observed when horizontal fractures are ignored. However, due to the comparably small length of the horizontal fracture, the effect is less pronounced.

Arrival Times and Dispersion
The distribution of residence times after n horizontal fractures is defined by and its Laplace transform is given by The first and second moments of the travel time are given by for j = 1, 2, or Thus, for the mean and the variance of residence time we obtain This means that the first moment and the variance are given by where h 1 and h 2 are the first and second moments of the residence time for a single fracture. They are given by (see Appendix C) To determine the fluid arrival times after n fractures, we add a constant time shift ΔT to the residence time in a single horizontal fracture. Thus, the quantities h 1 and h 2 are modified as We nondimensionalize time with respect to the critical time t c such that where τ m = t max /t c is found from Equation 30 as The arrival time moments are nondimensionalized accordingly as where the dimensionless τ 0 = ΔT/t c (Equations 12 and 29) is given by The equivalent flow velocity and dispersion coefficients are given in terms of the mean m 1 and variance σ 2 of the arrival times at a plane at z = nΔz, where Δz is the spacing between horizontal fractures, We nondimensionalize lengths by Δz and obtain In the following, we study the behavior of the nondimensional dispersion coefficient  D and the flow velocity  v as functions of the nondimensional times τ m and τ 0 .

Dimensionless Analysis of Flow Through a Fracture Network
To investigate the effect of nondimensional times τ m and τ 0 on the dimensionless flow velocity  v and dispersion coefficient  D (Equations 53 and 54), we conduct a multiparameter study. The minimum and maximum values of τ m and τ 0 are computed for all parameter combinations of Q 0 , l max , θ 0 , a, and Δz. As we carry out our study in 2D, the proposed ranges for the flow rate should be interpreted in terms of the maximum film thick-KORDILLA ET AL.

10.1029/2020WR028775
15 of 27 ness on the vertical surface. We vary Q 0 in the range of 1 × 10 −6 to 1 × 10 −2 m 2 s −1 , which yields a film thickness on the vertical surfaces between 63.7 µm and 1.5 mm according to Equation 24. This range has been observed and modeled in various studies (Dragila & Weisbrod, 2003;Dragila & Wheatcraft, 2001;Patnaik & Perez-Blanco, 1996;Tokunaga & Wan, 1997). However, it should be noted that flows in the upper range on the order of millimeters are known to develop a wavy-laminar structure with low and high amplitude waves traveling as a rolling film (Ghezzehei, 2004;Patnaik & Perez-Blanco, 1996) and therefore are likely to deviate from the simple film flow approximation. The horizontal fracture depth l max ranges from 0.01 to 4 m with an aperture a of 0.5-10 mm. The upper limit of the fracture length range has been chosen following Olson et al. (2007), who demonstrated that cumulative fracture lengths in fractured systems under various stress boundary conditions converge toward a maximum value of ∼4 m. While fracture apertures cover all length scales from micron to centimeter scales (Bonnet et al., 2001), apertures in the vadose zone are more likely to cover the upper portion of this range due to less overburden pressure and possibly enhanced dissolution (e.g., in karstic limestones). The chosen aperture range has for example been reported by Bahat (1987) and Weisbrod et al. (1998), who studied fractured chalk systems in the northern Negev desert of Israel. For apertures on the order of tens of millimeters, the assumption of capillary-driven flow in the horizontal fracture is likely violated in natural systems, and low-capillary free surface flows may develop. The static contact angle θ 0 is chosen to vary from 5° to 85° corresponding to a wetting regime that has been observed in various studies (Sobolev et al., 2000;Su et al., 1999Su et al., , 2001Tokunaga & Wan, 2001). Even lower or higher angles are possible, but this changes the conceptual framework. Near-zero contact angles may indicate spreading of a fluid, for example, because of micro-roughness or interaction with adsorbed films in the presence of a porous matrix (Shigorina et al., 2017;Tokunaga & Wan, 1997). Contact angles above 90° indicate hydrophobic conditions and fluids would require external forces beyond the capillary drag to enter a void space such as a fracture. The latter conditions are rarely met and have been observed in soil systems (Bachmann et al., 2000). Finally, the vertical fracture spacing Δz ranges from 0.1 to 25 m. For the sake of simplicity, we have chosen an equidistant setup for the vertical spacing, while natural fracture systems are often classified with respect to their regularity and follow a distribution (Hooker et al., 2013;Marrett et al., 2018). As we nondimensionalize the bulk properties with respect to the constant fracture spacing Δz, it is difficult to include a distribution in the analytical expressions derived here. However, in principle, the model can employ an irregular fracture spacing via suitable distributions in the transfer function φ and computation of the outflow (and dispersion) via Equation 32.
For the given parameter ranges, τ m can take values between 8.1 × 10 −6 and 1.4 × 10 5 and τ 0 between 3.2 × 10 −15 and 6.46. While the above chosen parameters are within feasible ranges, we further limit the relevant range of τ m and τ 0 by constraining the Reynolds numbers within the horizontal fracture. Here we calculate the critical Reynolds numbers as where the characteristic velocity v c is computed from the critical length and time We chose a maximum value of Re c = 150 to stay within the steady nonlinear laminar flow regime as that studied by Dybbs and Edwards (1984). values of τ 0 < 1 are close to reaching this constant maximum, while for τ 0 > 1, the dispersion increases for the considered range of τ m . The smaller initial gradient of   Δ / Δ m D (e.g., τ 0 = 7.5) is caused by the nonlinear Washburn dynamics within the horizontal fracture. For smaller values of τ m , the initial rapid (potentially plug-flow type when t max < t c ) infiltration dominates the bulk flow, while for higher values of τ m , the classical t scaling comes into effect and causes stronger dispersion of the breakthrough signal. Furthermore this example demonstrates how the ratio of τ m and τ 0 affects the nondimensional dispersion coefficient. Increasing the ratio of τ 0 /τ m strengthens the dominance of the vertical flow paths and therefore decreases the overall dispersion, which in our model entirely stems from the horizontal fracture imbibition. However, it should be noted that this effect is negligible for values of τ 0 < 10 −5 (Re c restricted) and already vanishes for τ 0 < 0.1. This is similar to the behavior displayed by the dimensional example (Figures 9 and 10 . Note that τ m ≥ 1 for the left plot because the strict analytical solution for the pure plug-flow regime (t max < t c , t max /t c < 1) does not cause any type of dispersion.  number of fractures and, therefore, the magnitude of horizontal imbibition (inversely related to the fracture spacing Δz) is positively correlated with the dispersion and for plug-flow-regime dynamics, no dispersion occurs (t max > t c , equivalent to a very high ratio of τ 0 /τ m or low values of τ m ). Figure 11 (right) demonstrates the dependence of the dimensionless flow velocity  v on the horizontal fracture timescale τ m for a range of τ 0 between 0 and 7.5. Two regimes can be observed. For low values of τ m , the velocity converges toward a constant value, while for higher τ m , the nondimensional velocity scales as in accordance with Equation 53. This transition occurs at the time τ m that increases with τ 0 /τ m because of the increased impact of vertical film flow dynamics and plug-flow type dynamics in the horizontal fracture. For τ 0 ≤ 0.01 the velocity scales as over nearly the entire range of feasible τ m values, and the average breakthrough velocities decline with increasing magnitude of the horizontal fracture imbibition (e.g., deeper or wider fractures, higher static contact angles). For even lower values of τ 0 < 10 −5 in the Re c -restricted range, a perfect scaling governs the functional relation between nondimensional velocity and horizontal fracture imbibition timescale, and no regime transition occurs.
Next, we discuss the dependence of  D and  v on the vertical fracture timescale τ 0 . Figure 12 (left) demonstrates the limited influence of the vertical fracture timescale τ 0 on the nondimensional dispersion. Only for extreme end-members of the parameter range beyond τ 0 ≈ 1 is the effect of vertical fracture flow strong enough to counteract the dispersive action of the horizontal fracture and introduce a reduction in dispersion. However, within critical Re c ranges,  D is only a function of τ m . For values of τ m > 10 3 the relationship converges toward the constant value of   0.08 D .
Similarly, the nondimensional velocity is independent of the vertical fracture flow timescale τ 0 within critical Re c ranges and is only dependent on the flow dynamics within the horizontal fracture encoded by τ m with a scaling behavior. It should be noted that this scaling holds for values of τ m < 1 (see Figure 11, right); however, here the nondimensional dispersion is   0 D and flow is entirely governed by plug flow in the horizontal fracture and film flow on vertical surfaces.

Conclusion and Outlook
In this work, we developed an analytical solution for partially saturated flow through an arbitrarily large sugar-cube type fracture network consisting of wide-aperture horizontal fractures intersected by a vertical fracture. Based on numerical observations using an SPH code and former laboratory studies, we treat the partitioning dynamics at the fracture intersection as a sequential process whereby the fluid is channeled from the upper vertical surface into the horizontal fracture and finally onto the lower vertical fracture surface. Flow within the horizontal fracture is shown to follow plug-flow theory until critical pressure thresholds are exceeded. After the breakthrough, horizontal infiltration is governed by a Washburn-type scaling until the maximum horizontal fracture depth is reached and the outflow at the bottom of the system equals the inflow rate at the top. To model flow through arbitrarily large networks of the same internal structure, we capture this process with an analytical transfer function and carry out a convolution of the constant input signal following linear response theory. Given the complex parameter space of fluid and geometric properties, we analyze the outflow dynamics in terms of nondimensional values of τ m and τ 0 that encode the timescales of flow in the horizontal and vertical fractures and relate them to the nondimensional dispersion coefficient  D and velocity  v. It is shown that within the feasible Reynolds number range, the dimensionless dispersion coefficient converges to the value of   0.08 D with increasing τ m and is nearly independent of τ 0 , that is, the flow in the vertical fracture does not have an impact on the dispersion coefficient. Furthermore, the bulk flow velocities are characterized by a scaling that holds for all relevant values of τ m and is independent of τ 0 within critical Re c ranges.
Our work demonstrates the importance of horizontal fractures as drivers for the (lateral) dispersive action within a mainly vertically oriented flow field. This conclusion clearly deviates from the classical piston-flow dynamics that are often assumed in the field(continuum)-scale flow models in fractured-porous systems Lange et al., 2010). As demonstrated in this work, neglecting the effects of Washburn-type dynamics (i.e., assuming pure piston flow) and/or neglecting the impact of horizontal fractures heavily alters the dispersive bulk properties and decreases (mean) breakthrough times. Furthermore, our work sheds KORDILLA ET AL.

10.1029/2020WR028775
18 of 27 light on the relation between integral signals at outlet boundaries (e.g., water table fluctuations within boreholes) and the internal system geometry that transforms input signals (precipitation, recharge) and mainly contributes to the dispersion and bulk velocity within the vadose zone.
In our analysis, we simplify the infiltration process in terms of the fracture-network geometry as well as the partitioning process and flow mode occurrence. Our study assumes film flow on all vertical surfaces. This assumption is often made in studies related to preferential flow in soil systems and the respective macropore structure (Bogner & Germann, 2019;Germann et al., 2007;Nimmo, 2010Nimmo, , 2012Nimmo & Perkins, 2018). However, other flow modes such as flow-rate-dependent droplets (slugs) and rivulets are likely to occur on fracture surfaces (Dippenaar & Van Rooy, 2016;Dragila & Weisbrod, 2003, 2004Ghezzehei & Or, 2005;Jones et al., 2017) and are known to affect partitioning at intersections T. R. Wood et al., 2005;Xue et al., 2020). While droplets are more likely to bypass intersections due to their extended height (as compared to films) and hence gravitational impact , we have also demonstrated that consecutive routing of droplet flows through arrays of horizontal fractures will nearly always facilitate the formation of film (rivulet) flows on the vertical wide-aperture surfaces after the first partitioning (Noffz et al., 2019). Subsequently, flow is mostly channeled into horizontal fractures without bypass, supporting the assumption of sequential flow dynamics made in this study.
Because our study is limited to a two-dimensional fracture network, the observed dependence of the nondimensional dispersion  D and velocity  v on the vertical and horizontal fracture flow timescales must be interpreted with care. For stable infiltration fronts, the infiltration dynamics of three-dimensional systems can be accurately recovered with two-dimensional models (e.g., Kordilla et al., 2017) using homogenization over the third dimension. However, flow on vertical fracture surfaces tends to develop instabilities even when these surfaces are perfectly smooth (Shigorina et al., 2019). Such front instabilities can contribute to fracture-specific channeling and additional dispersion. Front instabilities can develop in horizontal fractures as well, even though the formation of instabilities here is not caused by gravitational pull but is mainly a result of viscous forces and velocity variations due to changes in fracture aperture (roughness) and variations of the capillary radius . Apart from the dimensionality of the system, the geometrical properties such as the inclinations of the horizontal and vertical fractures are clearly simplified because natural systems often deviate from the classical "sugar-cube" geometry. Nonhorizontal fractures will decrease the breakthrough time t c when inclination angles are positive (relative to the horizontal) due to additional gravitational drag on the fluid entering the fracture. Negative inclination angles may not only increase the breakthrough times (due to increasing drag on the fluid front into the fracture) but may even lead to a complete change in the conceptual model if fluid leaves the horizontal fractures when reaching the maximum depth l max . For horizontal fractures and under wetting conditions, this effect does not occur even when the fracture is open at l max . The conceptual base for such processes would require adding a sink term to the transfer function and therefore increasing the bulk dispersion. Such effects are important when assessing flow convergence and integrative properties of fracture networks (Glass & LaViolette, 2004;LaViolette et al., 2003) within the vadose zone and add another level of complexity to the problem at hand.
In contrast to other studies, we focus on the case of vertical fractures with wide apertures. Here, the term "wide" should be interpreted with respect to the probability of fluid wetting opposing sites of the vertical fracture. For contact angles in the range of 25° to 75° droplet heights (neglecting the dynamic flattening due to movement) would be on the order of 0.75-2.4 mm, therefore setting a lower limit where one-sided flow would persist. Studies focusing on "narrow" vertical fractures often observe a slightly wider range of partitioning patterns that stem from the erratic uptake and emittance of potentially chaotic droplet patterns from T-type (Xue et al., 2020;Yang et al., 2019) and X-type intersections (e.g., Glass et al., 2003;T. R. Wood et al., 2005;T. Wood & Huang, 2015). While the majority of fractures under common geological conditions will belong to the "narrow" category, wide-aperture fractures are more likely to be found in the vadose zone where overburden pressure is limited and especially in karstic environments where fractures can be affected by dissolution (Dahan et al., 1999(Dahan et al., , 2000. Because most studies are still focusing on individual intersections dynamics, a unified theory for a broad range of apertures and partitioning dynamics is still to be developed.
Upscaling of individual processes, such as the intersection uptake and partitioning dynamics, remains one of the most challenging aspects in the current state of infiltration dynamics in fractured-porous media. In this work, we demonstrated how to bridge the gap between small-scale process and larger-scale KORDILLA ET AL.

10.1029/2020WR028775
19 of 27 bulk application using a simple convolution and the analytical derivation of (nondimensional) dispersion and velocity parameters. In its current form, the model assumes that convolution occurs over an arbitrary number of equally structured intersections. In principle, this could be extended to sequences of intersections with dynamic properties by introducing parameter distributions that reflect changes in the transfer function φ and therefore outflow Q n (t) over the range of encountered fractures n (Equation 32). However, while this would enhance the applicability to natural geological systems, analytical forms of  D and  v would be more difficult to derive. Inclusion of more complex partitioning of different flow mode dynamics, (e.g., droplets; Xue et al., 2020;Yang et al., 2019) would be an interesting, yet highly challenging extension because the uptake and release of such flows at intersections introduces a highly erratic and chaotic component.
As our model assumes impermeable fracture walls, we cannot model effects of porous matrix storage. This is justified for low-permeable porous systems, such as granites, or impermeable limestone surfaces. When considering the porous matrix, both the vertical fracture walls as well as the horizontal walls will retard the movement within the fracture (e.g., Buscheck et al., 1991) and could be introduced via suitable storage (sink) terms into the transfer function. The potential influence of the matrix imbibition depends on the consolidated material. Sandstones are typically at the upper end of the porosity and permeability (sorptivity) range and are considered for the following brief example. Assuming a porosity of θ = 0.2 and sorptivity of S = 0.15 cm s −0.5 (e.g., Berea sandstone, Kang et al., 2013), a first-order approximation for the potential fluid uptake in the studied examples can be made. For the sake of simplicity, we assume that fluid imbibes into the porous matrix in the horizontal direction across the entire length of the upper vertical fractures l L m 2 . The volume stored within the fracture until t c is A f = al c = 5 × 10 −4 m 2 . Therefore the imbibition into the matrix could potentially reduce the inflow volume and hence velocity during this time period by about 50%.
To derive consistent and process-based infiltration functions for fractured-porous media, it is crucial to unify the various observed patterns for partitioning dynamics across the scientific community. This will require further studies at laboratory and field scales to elucidate the shortcomings of each approach and obtain a suitable array of methods adjusted for the respective study setting and data availability. Given the strong impacts of climate-change-induced transformation of precipitation patterns (Black, 2009), water resources management, specifically in arid and semi-arid regions, requires enhanced models for recharge prediction that take into account the rapid preferential flow component that may substantially contribute to groundwater replenishment under high evapotranspiration and short but extreme rainfall conditions (Pachauri et al., 2014).

Appendix A: Smoothed Particle Hydrodynamics Model
We employ the SPH method to model free surface flow of water described by the Navier-Stokes (NS) equations, including the momentum and conservation equations respectively, subject to the Young boundary condition at the liquid-gas-solid interface, the Young-Laplace boundary condition at the water-air interface and the no-slip boundary condition at the liquid-solid interface Tartakovsky & Meakin, 2005a;Tartakovsky & Panchenko, 2016). Here, v is the velocity, P is the pressure, ρ is the density, μ is the viscosity, and g is the gravitational acceleration.
To simplify the solution of the incompressible NS Equation A1, we employ the weakly compressible formulation where the continuity equation is replaced with its compressible form dρ/dt = −ρ∇ ⋅ v, and the equation of state is used to close the resulting compressible NS equations: KORDILLA ET AL. where ρ 0 is the reference water density and P 0 is a background pressure. The speed of sound c 0 is chosen such that |δρ|/ρ ≤ 0.03, where |δρ| is the maximum absolute change in density. This condition is sufficient for fluid to behave as an incompressible fluid and to obtain an accurate pressure field (Morris et al., 1997).
The SPH discretization of the weakly compressible NS equation is: where e r  ij ij ij r  / is the unit vector pointing from particle i to particle j, summations are over all fluid (f) and/or solid (s) particles, and W is a two-dimensional Wendland kernel (Wendland, 1995) that establishes a smoothed interaction over the range h between particles. To simulate surface tension an additional pairwise interaction term (Tartakovsky & Meakin, 2005a) is employed that consists of two overlapping cubic spline functions W 1 and W 2 with a short-range repulsive and long-range attractive component controlled by coefficients A and B (Kordilla et al., 2013: where n i is the normal unit vector (see the definition in Pan et al. (2014), β = μ/λ is the friction coefficient, and λ is the artificial slip length. For most fluids, the real slip length is on the order of several nanometers, and using such a small λ would result in a prohibitively small time step in the SPH method. We demonstrate in Section 2.1 that the no-slip boundary condition can be accurately modeled by setting λ to be 100 times smaller than the domain size. We note that in Equation A3, the summation of the viscosity term is only over fluid particles, while the no-slip condition is entirely enforced via Equation A6.
The highly nonlinear form of the equation A2 generates sufficiently high pressure (in addition to the repulsive part of the interaction force) to prevent fluid particles from penetrating solid surfaces.
To integrate Equation A3, a modified Velocity Verlet time-stepping scheme is employed and time steps are constrained as follows Pan et al., 2014;Tartakovsky & Meakin, 2005c):

Appendix B: Spontaneous Imbibition
The volumetric flow rate through a fracture conceptualized as a parallel plate is governed by: where l is the penetration depth (the length over which the pressure gradient ΔP acts), a is the aperture and W is the fracture (unit) width. The change in volume over time can be rewritten in terms of the penetration depth into the fracture Separating the variables and integration from l c to l(t) on the left and from t c to t on the right side yields Note that this Washburn-type flow behavior is only valid at times much larger than t c . We do not model the transitional flow between the plug and Washburn modes, but represent the change between the two modes as abrupt and match the linear and square root behaviors at time t c .

Water Resources Research Appendix C: Moments
To calculate the moments of φ(t), it is advantageous to write the equation as follows,