Splitting Dynamics of Liquid Slugs at a T‐Junction

Understanding the mechanisms of liquid movement through fracture intersections is important for prediction of fluid flow and solute transport in unsaturated fractured media. Here we present a quasi‐static model to predict the dynamic splitting behavior of liquid slugs at a T‐junction, as a simplified representation of a fracture intersection and consisting of a main channel and a branch channel. The proposed model is validated against carefully controlled visualization experiments. We find that there exists a critical initial slug length at which the splitting behavior shifts from flow dominated by the main channel to that dominated by the branch. The influence of key parameters, including the inclination angle of the junction, the channel widths, and the dynamic contact angles, on the splitting dynamics is systematically investigated. We show that the splitting ratio depends nonmonotonically on the relative width of the branch channel to the main channel. Furthermore, it is demonstrated that the dynamic contact angles have a profound impact on the splitting ratios and meniscus velocities. It is shown that taking velocity‐dependent contact angle into account is essential to predict the meniscus velocities and the dynamic flow process.


Introduction
Water infiltration through fractured rock vadose zone has received extensive attention owing to its great influence on the application of hydrological cycle, contaminant treatment, geothermal energy production, and radioactive waste disposal (e.g., Berkowitz, 2002;Dobson et al., 2012;Essaid et al., 2015;Liu et al., 2003;Salve et al., 2008;Tsang et al., 2015). Field evidences show that focused, rapid, and deep infiltration along preferential pathways within fracture networks commonly occur (e.g., Dahan et al., 1999;Nimmo, 2012). However, owing to the random fracture aperture distribution and complex network structure, the preferential pathways of water infiltration show a significant spatial and temporal variability, making understanding and accurately predicting the complex behaviors of unsaturated flow a great challenge (e.g., Davidson et al., 1998;Nativ et al., 1995). In the past two decades, detailed investigations have been conducted to identify the fundamental processes of unsaturated flow and transport in fracture networks (e.g., Christensen et al., 2015;Glass et al., 2002Glass et al., , 2003Huang, Meakin, Liu, & McCreery, 2005). These studies have pointed to the critical role played by the fracture intersections in controlling unsaturated flow in fracture networks.
A number of experimental, theoretical, and numerical studies have investigated the fluid dynamics at fracture intersections of different structural types, such as X-, Y-, and T-shaped intersection (e.g., Dragila & Weisbrod, 2004;Kordilla et al., 2017;Wood et al., 2005). It has been well demonstrated by experiments that intersections can act as capillary barriers which cause local accumulation of the wetting fluid . The strength of this barrier effect depends on the combined influence of fracture apertures, intersection geometry, and wetting properties (e.g., Dragila & Weisbrod, 2004;Ji et al., 2006). Once the hydrostatic pressure of the accumulated fluid exceeds the capillary barrier, the liquid will quickly discharge and then the liquid pooling will be reestablished. The cycles of accumulation and discharge not only affect the size and frequency of the infiltrating liquid locally (Wood et al., 2005), but also cause flow pulsation in a fracture network . During liquid discharging and passing a bifurcating intersection, flow path selection and volume partitioning occur, which lead to flow pathway shifts and structural evolution in the networks (Glass & LaViolette, 2004). To better understand the liquid partitioning process, visualization experiments (Kordilla et al., 2017;Yang et al., 2019) have been performed to examine the fluid splitting behavior under different flow modes, from discrete droplet flows to continuous rivulets. The splitting process has been shown to exhibit significant dynamic characteristics in the meniscus velocities and the liquid branch lengths (Yang et al., 2019). The splitting dynamics is again influenced by multiple factors, including flow rate, intersection geometry, and hysteretic contact angles (e.g., Kordilla et al., 2017;Noffz et al., 2019). However, most previous experiments have not been designed to have carefully controlled conditions of the above influencing factors to allow for validation of mechanistic models of dynamic splitting at bifurcating junctions. Detailed mechanisms of the dynamic splitting process at fracture intersections remain not fully understood.
Wetting fluid flow through unsaturated factures and fracture intersections involves a complex key process, the motion of liquid-air contact lines. It is well understood that the apparent contact angle is not only influenced by the surface properties (e.g., Bonn et al., 2009;Hirasaki, 1991;Snoeijer & Andreotti, 2013) but also the contact line velocities (both direction and magnitude), giving rise to capillary hysteresis and dynamic contact angles (e.g., Hoffman, 1975;Lei et al., 2018;Petrov et al., 2003). To predict the relationship between dynamic contact angle and contact line velocity, theoretical frameworks have been proposed, including hydrodynamic models (e.g., Cox, 1986;Huh & Scriven, 1971;Voinov, 1976), molecular-kinetic models (Blake & Haynes, 1969), and combined molecular-hydrodynamic models (e.g., Brochard-Wyart & Gennes, 1992;Petrov & Petrov, 1992). Introducing the hydrodynamic theories in flow models has been shown to improve the prediction of wetting droplet flow velocity in a single vertical fracture (e.g., Bico & Quéré, 2001;Dragila & Weisbrod, 2003;Or & Ghezzehei, 2006) over the prediction based on static contact angles. However, the effect of velocity-dependent contact angles on the wetting liquid partitioning behavior at an intersection has not been explored.
A fracture intersection necessarily involves a nonhorizontal fracture. This means the wetting liquid flow through a fracture intersection or junction is controlled by the complex interplay between gravity, capillary, and viscous forces (Wood et al., 2005). For liquid drops/slugs in a fracture before entering an intersection, when the vertical length of droplets/liquid slugs is larger than the capillary length (l c = (σ/Δρg) 1/2 , where σ is the interfacial tension, Δρ is the density difference, and g is the gravitational constant), gravity or buoyancy becomes important in the droplet flow dynamics. This comparison can be quantified by the dimensionless Bond number, Bo = Δρg (Lsinα) 2 /σ, where L is the length of the droplet/liquid slug, α is the inclination angle of the fracture containing the droplet/slug. Simple static force balance analysis has been adopted in previous studies to qualitatively explain the flow accumulation and release behavior (e.g., Glass et al., 2002;Ji et al., 2006;Wood et al., 2005). However, the static force balance analysis is unable to describe the splitting dynamics, since the forces vary with time during fluid splitting. Controlling mechanisms of gravity-driven wetting liquid flow at an intersection remain to be elucidated. A comprehensive and quantitative description of dynamic fluid splitting at an intersection is still lacking.
Here we perform carefully controlled experiments of dynamic splitting of liquid slugs at a T-junction as a key representative process for wetting fluids passing a fracture intersection. We propose a quasi-static theory to describe the liquid splitting process at the T-junction. We develop a semianalytical solution approach to predict the dynamic splitting behavior which is quantified through the evolution of liquid branch lengths, meniscus velocities, and splitting volume ratios. A critical initial slug length is identified that separates the splitting behavior into two types of flow regimes. The key parameters affecting the splitting process are systematically investigated, including the initial liquid slug length, channel widths, inclination angle of the junction, and the dynamic contact angle.

Theory and Models
We consider a two-dimensional system of flow through a junction, a reasonable simplification for the water infiltration in fractured rocks where localized channel flows usually exist (Su et al., 2004). We present quantitative descriptions of liquid slug motion and splitting at a T-junction. The interplay between gravity, velocity-dependent capillary forces, and size of the slug during splitting is detailed. We highlight the important role of capillary hysteresis and dynamic contact angles in the splitting dynamics.

Dynamic Contact Angles and Liquid Slug Motion Within a Straight Channel
Consider a liquid slug of length L in a single idealized channel (see Figure 1a). The droplet motion is controlled by the capillary, gravity, and viscous forces. The slug's weight acts as a driving force to pull the droplet downward, while the capillary force difference between advancing and receding menisci acts to resist 10.1029/2020WR027730

Water Resources Research
motion. For a stationary liquid slug, the force balance is written as Lsinα = ψ a − ψ r , where α is the inclination angle of the channel relative to the horizontal plane, ψ a and ψ r are the capillary pressures (in units of hydraulic head) of the advancing and receding menisci, respectively. Ignoring the minor principal curvature, we calculate ψ a and ψ r by using the Young-Laplace equation, where θ a and θ r are the advancing and receding contact angles, respectively; σ is the interfacial tension; w is the channel width; ρ is the liquid density; g is the gravity acceleration. Assuming the contact angles satisfy the relationship: cosθ s = (cosθ a + cosθ r )/2 (e.g., Andrieu et al., 1994;Decker & Garoff, 1996), where θ s is the static contact angle which can be measured directly. By putting Equations 1a and 1b into the force balance equation, one can obtain the relationship between advancing and receding contact angles and the length of the stationary slug, given a static contact angle.
The receding contact angle decreases with increasing slug length, whereas the advancing contact angle increases. For water-air interfaces the critical contact angles θ * a and θ * r exist (e.g., Or & Ghezzehei, 2007;Su et al., 2004). When the advancing contact angle is greater than θ * a or the receding contact angle less than θ * r , the interface will slip, resulting in slug migration within the channel. Here we define a threshold length L * for triggering slug movement as When the liquid length L is larger than the threshold length L * , the slug maintains a steady velocity u, under the balance of capillary, gravity, and viscous forces (e.g., Or & Ghezzehei, 2007;Su et al., 2004): where μ is the dynamic viscosity. When the slug is in motion, the values of advancing and receding contact angles vary with the interfacial velocity. Here we adopt the hydrodynamic model for a liquid-air system proposed by Voinov (1976). This model (see Figure 1b) provides a simple scaling relationship between the contact angle and the capillary number Ca, where Ca = μu/σ, ξ 1 and ξ 2 are the shape parameters which can be obtained by fitting Equations 4a and 4b to experimentally measured data of interfacial velocities and advancing and receding contact angles.
In an idealized channel, the gravity can be eventually balanced by capillary and viscous forces. Given the droplets length and other known parameters, the unknowns, the velocity u and the dynamic contact angles θ a and θ r , can be obtained from Equations 3, 4a, and 4b by an iterative solution.

Quasi-Static Model of Liquid Slug Splitting at a T-Junction
Unlike the steady motion of a liquid slug in a single channel, the liquid slug splitting behavior at a T-junction is an unsteady process, involving transient slug lengths and interfacial velocities, as well as the associated velocity-dependent contact angles. Previous work (Yang et al., 2019) suggests that when the width of the branch channel is larger than that of the main channel, no splitting occurs and the wetting fluid bypasses the branch channel, which is a scenario same as that in a straight channel. Therefore, we focus on a geometry where the branch channel width is smaller than that of main channel. It is possible for the liquid slug to be completely diverted into the branch channel (e.g., when the inclination angle is small); we consider this scenario as a special case of splitting with a cumulative splitting ratio (as quantified later in section 3.1) of 1.
Here we propose a quasi-static model to describe liquid slug splitting at a T-junction. We assume that the motion of menisci during the splitting is completed in discrete steps of a small time interval Δt, within 10.1029/2020WR027730

Water Resources Research
which a force balance of the deformed (due to the bifurcation) liquid slug is maintained (Figure 1c). It is further assumed that during splitting, the flow velocities within the two individual channels C1 and C2 and the dynamic contact angles can be solved using the solution approach above for Equations 3, 4a, and 4b.
To describe the liquid motion, we denote the liquid length within channel C1 by L C1 (t), the length within the lower part of C1 by L C1,d (t), and the length within channel C2 by L C2 (t), see Figures 1c-1e. At the initial time (t = 0), the liquid is assumed to have just established interface menisci within the C2 and the lower part of C1, L C1 (0) = L, L C1,d (0) = w 2 /2, L C2 (0) = w 1 /2, where w 1 and w 2 are the channel width of C1 and C2, respectively. The splitting ends when the liquid completely passes the T-junction, L C1,d (t) = L C1 (t), i.e., no liquid is left above the junction. In an arbitrary time interval, the receding interface velocity u C1 r of channel C1 is given as where L C1 (t) = L − L C2 (t)w 2 /w 1 ; θ C1 r (t) and θ C1 a (t) are the corresponding transient receding and advancing contact angles.
Within channel C2, the fluid motion is driven by the pressure gradient and gravity. The interface velocity u C2 a (t) can be expressed by where ψ C2 a (t) is the capillary pressure in channel C2; P w is the fluid pressure (in unit of hydraulic head) at the junction, which is related to the capillary force ψ C1 a (t) and the liquid length L C1,d as

Water Resources Research
Fluid mass conservation provides a closure to the nonlinearly coupled equations to enable a solution. The supplied flow rate at the junction, u C1 r w 1 , is equal to the sum of the flow into channel C2 and the lower part of channel C1. The advancing interface velocity, u C1 a , can be expressed as There are six unknowns: u C1 r , u C1 a , u C2 a , θ C1 r , θ C1 a , and θ C2 a , and six equations: three equations 5-7 and three Voinov equations linking the dynamic contact angles and velocities, i.e., it is a determined system. However, a nested iterative solution strategy is required to solve the system. Here we describe a relatively simple procedure to achieve at a satisfactory solution (supporting information Figure S1). At each time step, we first obtain initial guesses of u C1 r , u C1 a , θ C1 a , and θ C1 a (assume u C1 r ¼ u C1 a ) using Equations 4a, 4b, and 5. Then, u C2 a and θ C2 a are computed within an internal iteration scheme to satisfy Equations 4a, 4b, and 6 with a preset tolerance of 0.1°. Next, we recalculate u C1 a according to Equation 7 and θ C1 a according to Equations 4a and 4b. Then we check convergence with respect to θ C1 a : if there is no convergence, we update θ C1 a and u C1 a and then redo the above steps of the time step. Upon convergence, we update the slug lengths, obtain the new meniscus locations, and continue the next time step. The procedure ends when the splitting is finished.

Laboratory Experiment
To test the quasi-static model described above, an experimental setup for liquid breakup at a T-junction is designed, as shown in Figure 2. The T-junction flow cell consists of three thick (3.0 mm), smooth fused-quartz plates with the polished surface and sharp edges. The quartz plates are aligned with controlled channel widths and then attached to a piece of acrylic plate by four spacers (3 mm thick) with UV curing glue. The acrylic plate is then fixed to a support system which is capable of allowing the T-junction to adjust its inclination angle α. The junction divides channel C1 and C2 into three parts, each being 7 cm in length. The channel widths are carefully controlled by feeler gauges. The accuracy of channel widths is double-checked with an optic microscope after the junction assemblage. We set the channel widths w 1 = 0.83 mm and w 2 = 0.33 mm to ensure the occurrence of liquid splitting and to encourage competition between capillary force and gravity.
At the top of the main channel C1, dyed water (0.075 g capsicum red per liter of pure water) is injected using a metal capillary tube which is connected to a syringe pump (Harvard Apparatus 703007, 0.25% injection accuracy) through soft silicon tubing. As soon as the injected water flows out of the capillary tube, the water spans the channel width and forms a liquid slug. We control the injection volume to acquire liquid slugs of various initial lengths.
Uniform and stable static contact angle is a key to ensure the experimental repeatability. Before each experiment, we clean the T-junction cell through plasma cleaning and subsequent repeated cycles of vacuum drying and ultrapure water immersion (Yang et al., 2019). We find that the static contact angle of the channel surface eventually reaches a stable state at 49.2°± 3.5°( Figure S2). Besides, we also measure the dynamic contact angles and find that with increasing interface velocity from 0.2 to 1.5 cm/s, the receding contact angle decreases from about 30°to 25°and the advancing contact angle increases from about 75°to 90°( Figure S3). These data allow us to obtain the shape parameters in Equations 4a and 4b.
Three different inclination angles (α = 90°, 70°, 50°) are controlled in the experiments to investigate the influence of gravity on slug splitting behavior. In each case, we also accurately control the injected volumes (V = 25, 30, 35, 40, 45, 50 μl) and the corresponding initial slug lengths (L = 1.05, 1.25, 1.45, 1.65, 1.85, 2.05 cm). The experimental processes are recorded by a camera at a frequency of 5 frames/s. Experiments for each slug length and each inclination angle are repeated five times. The Bond number (Bo = ΔρgL 2 /σ) ranges between 15 and 58.

Results
In this section, experimental results are presented together with the results from the quasi-static model described in section 2.1. The comparison between experiments and model also serves as a validation of the theory. Then, using the model, we systematically investigate the impact of slug length L, channel width ratio w 2 /w 1 , inclination angle α, and dynamics contact angles θ r and θ a on the splitting behavior.

Characteristic Behaviors of Liquid Splitting at the T-Junction 3.1.1. Splitting Characteristics
When the liquid front arrives at the junction, splitting occurs under the influence of capillary force in channel C2, resulting in part of water invading channel C2 and the rest entering the lower part of channel C1. After completely passing the T-junction, the liquid slug is divided into two daughter slugs. We experimentally find that the initial slug length L has a strong influence on the volume ratio of the daughter slugs. When L = 1.25 cm, most liquid enters C2 (Figures 3a-3d). When L = 1.65 cm, most water prefers to enter the lower part of C1 (Figures 4a-4d). According to the volume ratio of the two daughter slugs, the splitting behavior can be generally divided into two typical patterns: (1) type I for small L, flow dominated by channel C2; (2) type II large L, flow dominated by channel C1.
The detailed evolution of liquid lengths and interface velocities for the case of initial length L = 1.25 cm is shown in Figures 3e and 3f. We reproduce the experiments by the quasi-static model. The basic parameters used in the model are w 1 = 0.83 mm, w 2 = 0.33 mm, σ = 0.072 N/m, ρ = 1,000 kg/m 3 , g = 9.8 N/kg, μ = 1 × 10 −3 Pa·s, θ s = 50°, ξ 1 = 12,000, ξ 2 = 500, L * = 7 mm. As can be seen, the modeling results of liquid lengths well match the experimental observations in terms of the overall trends. The maximum deviation of slug lengths from experiments for L = 1.25 cm is about 0.2 cm, which is an acceptable error. For L = 1.05 cm, the maximum deviation is comparable.
The modeling results also show that velocities u C1 a , u C1 r , and u C2 a a all rapidly decline due to fluid partitioning into the side branch. The ratio between u C2 a and u C1 r continuously increases from 1.25 at t = 0 to 2.52 (which equals to w 1 /w 2 ) at t = 1.36 s, indicating that the capillary force of channel C2 dominates the splitting process. According to the evolution of the velocity u C1 a , the splitting process can be divided into two stages. During the

10.1029/2020WR027730
Water Resources Research first stage, the fluid arriving at the junction flows through both exits with velocities u C1 a and u C2 a , since the flow rate in channel C2 is less than the supplied flow rate above the junction, u C1 r w 1 > u C2 a w 2 . The separation between the first and the second stage is marked by the velocity u C1 a ¼ 0. In the second stage, the liquid no longer enters the lower part of C1, and the receding interfacial velocity u C1 r is controlled by the advancing interfacial velocity u C2 a in C2. The evolution of the type II splitting for large L = 1.64 cm is presented in Figure 5. In this scenario, most of the liquid arriving at the junction enters the lower part of C1. The splitting process is excellently captured by the model (Figure 4e). The deviation between predictions and experimental observations is smaller than that in the cases of L = 1.25 cm. During the splitting process, the interface velocities also decrease continuously until the fluid completely passes the junction. The advancing interface velocity u C2 a is always smaller than the receding velocity u C1 r , indicating that the gravity of the liquid in C1 dominates the splitting behavior. Different from the type I flow pattern, velocities u C1 a and u C2 r slowly converge with time, which is another indication of the flow dominance in C1. This result implies that the more water entering the lower part of channel C1, the more unfavorable the capillary-dominated flow in channel C2 is.

Splitting Ratio
We use the dimensionless parameters t * and η to quantify the dynamic splitting process. The dimensionless time t * is defined as the volume ratio of water passing the junction to the initial volume. The splitting ratio η(t * ) is a transient variable defined as the flow rate of water entering channel C2 to the total flow rate (as measured at the receding meniscus) at time t * , The evolution of splitting ratio against time for initial liquid slug lengths from 1.05 to 2.05 cm is shown in Figure 5a. The splitting processes show a stark contrast between type I (L ≤ 1.25 cm) and type II flow According to the quasi-static calculations, we find that there is a critical initial length (L crit = 1.41 cm) at which the splitting ratio maintains a stable value. The critical length can be used as a criterion to

Water Resources Research
determine the flow regimes of dynamic splitting at the T-junction. However, a fully analytical expression for L crit is difficult to obtain due to the strong nonlinearity of system.
The cumulative splitting ratio η * is also analyzed, η * = V C2 /V; i.e., η * is defined as the volume ratio of the daughter liquid slug within channel C2 to the total volume. The relationship between η * and the initial length L is presented in Figure 5b. Generally, the cumulative splitting ratio is smaller for longer liquid slugs. Consistent with η-t * data in Figure 5a, the critical liquid slug length is at the inflection point of the obtained η * -L relationship, separating the two flow regimes (see Figure 5b). Within the type I flow regime, η * is strongly sensitive to the variation of L, indicating that the key controlling force of the splitting behavior is gravity. When L is much larger than L crit , η * is insensitive to the change of L. This is because in the type II splitting regime, the flow in the branch channel is limited by the capillary force in C2. The dynamic contact angle acts as a flow regulator: the higher the flow rate, the smaller the capillary force in channel C2.

Effect of Inclination Angle
Previous studies have shown that the inclination angle is an important parameter in the wetting fluid flow in a conduit (Su et al., 2004). When the T-junction is inclined, gravitational effects are present in the liquid flow of both channels C1 and C2. Note here the inclination angle of 90°means channel C1 is vertical. Our experimental results of liquid splitting with different inclination angles (α = 70°, 50°) for L = 1.65 cm are shown in Figure 6. The model calculations are again in excellent agreement with the experimental data ( Figure 6a). As expected, reducing the inclination angle α leads to more water invading channel C2. As the inclination angle α lowers from 70°to 50°, the initial velocity u C1 r drops from 1.3 to 0.8 cm/s (Figure 6b), giving more time for the water to enter channel C2. When α = 50°, the velocity u C2 a in C2 first remains constant before the meniscus in C1 stops advancing, and then rises as the liquid in the upper C1 continue to feed C2 (Figure 6b).

Water Resources Research
The evolution of splitting ratio η with time t * for the T-junction with different inclination angles is presented in Figures 6c-6d. Decreasing the inclination angle α leads to a larger critical liquid slug length: when α = 70°, L crit = 1.75 cm (compare with 1.41 cm at α = 90°). When α further decreases to 50°, the splitting behavior completely shifts to the type I flow regime for L between 1.05 and 2.05 cm. As a result, given an initial liquid slug length L, the splitting ratio is higher for a smaller inclination angle.

Effect of Channel Widths
In addition to the inclination angle, the channel widths w 2 and w 1 also have a strong impact on the splitting behavior. The flow capacity in channel C2 is determined by the meniscus advancing velocity u C2 a and the conductivity or permeability. On the one hand, when w 2 decreases, u C2 a is expected to rise as the capillary force in C2 becomes stronger, which is favorable to a higher flow rate in C2. On the other hand, the permeability of the channel C2 decreases with w 2 as k C2 = w 2 2 /12, resulting in a reduction of flow. In other words, there exist two competing effects of the channel width. Thus, one would expect a nonmonotonic dependence of the splitting ratio on the channel width w 2 . Here, because no experimental data are available, we only use the developed model to quantitatively explore the influence of the channel width ratio w 2 /w 1 on the splitting ratio η and the cumulative ratio η * by varying w 2 . The inclination angle α is set to 90°. Figures 7a and 7b, when w 2 /w 1 changes from 0.2 to 0.4, the transient splitting ratio generally becomes higher for a given initial liquid slug length, especially for the Type I splitting regime (small L).

Water Resources Research
However, as w 2 /w 1 increases from 0.4 to 0.6 (compare Figures 7b and 7c), η is noticeably smaller for L between 1.25 and 1.45 cm. When w 2 /w 1 is large enough (Figure 7d), the transient splitting ratio η even drops to 0 (i.e., the fluid no longer invades C2) before the slug completely penetrates the T-junction.
Furthermore, we find, by the quasi-static model, that the critical length L crit depends nonmonotonically on w 2 /w 1 (Figure 7f). At about w 2 /w 1 = 0.4, L crit reaches the maximum, about 1.43 cm, for α = 90°(vertical C1). The relationship between w 2 /w 1 and L crit can be fitted by a parabola (Figure 7f). When w 2 /w 1 < 0.4, L crit increases with w 2 /w 1 , indicating that the hydraulic conductivity of C2 is the main limiting factor for water partitioning into C2. When w 2 /w 1 > 0.4, L crit declines with the width ratio, which is a sign that the capillary force in C2 is taking control of the dynamic splitting. When w 2 /w 1 > 0.7, the splitting behaviors for L between 1.05 and 2.05 cm all belong to the type II regime. The cumulative splitting ratio η * for different w 2 /w 1 values is presented in Figure 7e. A marked nonmonotonic dependence of η * on w 2 /w 1 can be observed for all initial lengths except for L = 2.05 cm (in which case the splitting is dominated by gravity-driven flow in C1). The shorter the initial length, the stronger the nonmonotonicity. This also manifests the salient nonlinearity caused by the complex interplay between gravity and velocity-dependent capillarity during liquid splitting.

Effect of Dynamic Contact Angles
The influence of dynamic contact angles on the splitting behavior is studied by varying the critical contact angles, θ * a , θ * r , and the shape parameters, ξ 1 , ξ 2 (see Equations 4a and 4b). Since these parameters cannot be controlled independently in experiments, here we only present a modeling analysis. Here we set α to 90°, and the channel widths to w 1 = 0.83 mm, w 2 = 0.33 mm.

The Receding Contact Angle
To investigate the effect of the receding contact angle, we set the Voinov theory-based advancing contact angle as θ * a ¼ 60°and ξ 1 = 12,000 (see dashed line in the inset of Figure 8a). Different values of θ * r ranging from 10°to 40°and of ξ 2 between 50 and 1,000 are tested (Figures 8a-8b, insets). The results show that

Water Resources Research
with the decrease of the critical contact angle θ * r (with ξ 2 fixed at 250), the cumulative splitting ratio η * increases obviously, especially for small slugs (Figure 8a). This is attributed to the fact that a smaller θ * r leads to a significant decrease of the velocity u C1 r compared with u C2 r , which favors liquid partitioning into the branch C2 (Figure 8c). This conclusion can be also supported by the variation of the initial splitting ratio η(t → 0) at a small L = 1.25 cm (Figure 8c). The higher the initial splitting ratio η is, the higher the cumulative splitting ratio for type I regime becomes. However, as the initial slug length L increases to above 1.7 cm, the effect of θ * r gradually diminishes, overshadowed by the gravity effect. The shape parameter ξ 2 dictates the relationship between the contact angle and the capillary number (Equation 4b). It can be seen from Figure 8b that η * appears to be insensitive to changes in ξ 2 . This is because as the parameter ξ 2 increases from 50 to 1,000, the contact angle shows a significant change only for interfacial velocity u > 1 cm/s, while the interfacial velocities, given the splitting scenario and the parameter space considered in this study, are mostly <1 cm/s. Besides, the influence of ξ 2 can also be explained by the evolution of the initial splitting ratio η(t → 0). With ξ 2 increasing from 50 to 1,000, η only increase from 0.46 to 0.56 (Figure 8d).

The Advancing Contact Angle
The impact of the advancing contact angle on the cumulative splitting ratio η * is investigated by using different values of θ * a (= 40°, 50°, 60°, 70°) and the shape parameter ξ 1 (= 1,000, 11,000, 21,000, 31,000). With increasing advancing contact angle, both the velocities u C1 r and u C2 a decrease. On the one hand, it is favorable for liquid splitting into channel C2 with u C1 r decreases. On the other hand, the decrease of u C2 a is unfavorable for the splitting behavior. Thus, the cumulative splitting ratio depends on the relative change of u C1 r and u C2 a . Interestingly, we find that η * presents the opposite evolution for type I regime (flow dominated by channel C2) and type II regime (flow dominated by channel C1) (Figure 9a). When L <~1.45 cm, η * increases with θ * a ;

10.1029/2020WR027730
Water Resources Research otherwise, η * decreases with θ * a . Qualitatively, this can also be explained by the evolution of initial interface velocities, u C1 r and u C2 a , and the corresponding initial splitting ratio η(t → 0). When L is small (type I regime), with increasing θ * a , u C1 r will decrease significantly compared with u C2 a (Figure 9c), resulting a larger η(t → 0), which is in favor of liquid entering channel C2. This is owing to the fact that the variation of advancing contact angle is more sensitive to the smaller velocity (Figure 9a inset). Conversely, when L is large (type II regime, u C1 r > u C2 a , η(t → 0) decreases with the increase of θ * a . The shape parameter ξ 1 regulates the relationship between the advancing contact angle and the interfacial velocity. We study the influence of ξ 1 on the splitting with θ * a ¼ 60°and θ * r ¼ 30°. The results show that ξ 1 , compared with θ * a , has a much less appreciable impact on η * (Figure 9b), even though the meniscus velocities show a strong dependence on ξ 1 (Figure 9b inset). Furthermore, we find that when ξ 1 = 1,000, the advancing contact angle approaches a constant for interface velocity between 0 and 5 cm/s, resulting in increase of u C1 r and u C2 a by near an order of magnitude. This result implies that the shape parameters ξ 1 and ξ 2 play an important role in controlling the interface velocities.

Slugs Splitting at the Junction
When a liquid slug passes a dry fracture intersection, there exist two different flow behaviors: splitting and bypass. Here bypass means the liquid flows only in the main channel and ignores the branch channel, i.e., the cumulative η * = 0. The criterion for separating the two behaviors can be determined using the interface velocity u C2 a . When the interface velocity u C2 a > 0, slug splitting occurs; otherwise, bypass occurs. According to Equation 6, the critical condition, u C2 a > 0, for the splitting behavior can be expressed as This equation means that the flow behavior is controlled by the competition between channel C2 and the lower part of C1. Consider an initial status, L C1,d = 0 and L C2 = 0, and assuming θ C1 a ¼ θ C2 a , we find that the wetting liquid invades channel C2 only when the widths (apertures) satisfy w 2 < w 1 . This is consistent with previous experimental studies (e.g., Ji et al., 2006;Yang et al., 2019). However, it should be noted that the competition between the two exits at the intersection is also mediated by the velocity-dependent contact angles. In addition, the lengths L C1,d and L C2 may not be initially zero due to preexisting liquids. Thus, the condition w 2 < w 1 for the occurrence of splitting may not always be applicable.
When splitting occurs, the transient splitting ratio is determined by the relative magnitude of interface velocities u C1 r and u C2 a , which are dependent on and, at the same time, alter the capillary and gravitational forces. As shown in Equation 5, the interface velocity u C1 r is related to the channel width, inclination angle, and slug length in channel C1. Longer slugs tend to have higher velocity u C1 r , which is in favor of flow dominated by channel C1. On the other hand, Equation 6 indicates that the interface velocity u C2 a is not only affected by channel C2 but also by the pressure at the junction, P w . During the splitting process, the more liquid enters the lower part of the main channel, the lower the pressure P w is, which leads to a smaller velocity u C2 a in C2. This explains the observation that a preexisting hanging column of water in C1 inhibits the invasion of the channel C2 . Additionally, according to Equation 7, the liquid is completely diverted into channel C2 when u C1 a t ð Þ ¼ 0. This special behavior with a cumulative splitting ratio of 1 tends to occur for small initial slug lengths and small inclination angles (see Table S1).
Two types of splitting flow regimes are distinguished: the branch channel-dominated flow and the main channel-dominated flow, which are found to be separated by a critical initial slug length L crit . This length L crit is important as it gives us a criterion to determine the main flow direction at junctions/intersections and even in simple networks. This length is a function of the inclination angle, the wetting parameters, and the channel widths/apertures. In Figure 7f, we present a fitted function of the dependence of L crit on w 2 /w 1 . However, it is currently difficult to derive an analytical expression of the dependence considering all influencing factors. Further investigations are thus needed in this regard.

Water Resources Research
It is important to further clarify the role of viscous force in liquid splitting at an intersection. The viscous pressure drop ΔP can be estimated as ΔP = 12 μuL/w 2 ; given u = 2.0 cm/s, w = 0.33 mm, and L = 1 cm, ΔP = 22 Pa. This is about one order of magnitude smaller than the capillary pressure 2σcosθ s /w (= 280 Pa for θ s = 50°and w = 0.33 mm). In this sense, the viscous force seems to play an insignificant role for the widths considered in this study; this is also suggested by previous analyses (e.g., Ji et al., 2006;Wood et al., 2005). However, our work here also draws attention to a previously overlooked issue: the effect of visco-capillary balance at the contact region, as embedded in the Voinov (1976) theory, on wetting liquid flow through intersections. Through the model calculations, we find that ignoring the velocity dependence of the apparent contact angle can overpredict the interface velocities by over one order of magnitude (see Figure 9d).

Experimental and Conceptual Limitations
To reveal the essential mechanisms of fluid splitting at an intersection, we make some necessary simplifications in this study. Our experimental and model scenarios correspond to the droplet flow mode of unsaturated flow (Kordilla et al., 2017). However, in natural fractures, other flow patterns, such as rivulet and film flow, can also occur, which can lead to alternative splitting behaviors at fracture intersections (e.g., Dragila & Weisbrod, 2004;Kordilla et al., 2017;Yang et al., 2019). In addition, we use a two-dimensional system in this study to investigate wetting liquid flow through fracture intersections, while in reality water infiltration in fractured media presents a three-dimensional flow problem. Nevertheless, this simplification can be considered reasonable as localized, channeling behavior has been recognized to be commonplace in unsaturated flows (Su et al., 2004).
During the calculation of the interface velocity, we have assumed the receding meniscus to be in a stable state without slip. In fact, when the droplet velocity is high, the droplet slip will occur, leaving a macroscopic film separated from the receding interface (Zhao et al., 2018). This film will increase the droplet/slug velocity and should be considered in the calculation of capillary pressures (Bico & Quéré, 2001). Additionally, only the scenario of initially dry intersections is considered here, while in natural conditions preexisting fluid either as pooling water or as residual droplets/liquid bridges can make the problem more complex. Future work is needed to address the above limitations.

Extensions and Applications of the Model
In this study, the idealized case of liquid slug motion in smooth channels is considered, while in nature channels can have rough surfaces and be nonstraight (Shigorina et al., 2019). Rough surfaces can significantly enhance the effect of contact angle hysteresis (Bonn et al., 2009), resulting in a smaller interface velocities of the liquid slug (Su et al., 2004). In this respect, the effect of roughness on slug motion and splitting can be possibly accounted for by adjusting the parameters (θ* a, θ* r, ξ 1 , ξ 2 ) in Equations 4a and 4b. To fully address this effect, more experimental and modeling work is needed. For nonstraight channels, the slug length in the gravity force terms can be replaced by its projection length in the gravity direction. Additionally, by changing the angle in the pressure gradient term (cosα) in Equation 6, the model can be straightforwardly extended to predict slug splitting at nonperpendicular junctions, i.e., the model is not limited to T-shaped geometries.
In subsurface environmental and engineering applications, such as underground waste isolation, it is important to be able to predict seepage into or flow diversion around natural cavities and excavated openings (Finsterle et al., 2003;Ghezzehei, 2005). The unsaturated flow mechanism in the fractured rocks surrounding the cavities is the key in this prediction. Under low flow conditions, the droplet/slug flow mode commonly prevails, giving rise to spatially and temporally discontinuous flows (Wood et al., 2005). The resulting complex flow diversion and fluctuation behaviors may not be well captured by continuum unsaturated flow models with effective parameters (e.g., Finsterle et al., 2003). Under this circumstance, our mechanistic model of fluid splitting can be particularly useful in predicting the flow diversion through fractures around the cavities. The liquid splitting ratio as a function of initial slug length (see, e.g., Figure 5b) can also be incorporated into network models to study macroscopic flow features (e.g., flow convergence/divergence) in fracture networks, which may provide insights for understanding field-scale infiltration and contaminant transport problems in the vadose zone.

Conclusions
In this work, the splitting dynamics of liquid slugs at a T-junction has been studied. A quasi-static model and a semianalytical solution approach are proposed to predict the dynamic splitting behavior. The model is shown to well reproduce the liquid splitting process as observed in visualization experiments. Combining model predictions and experimental data, we find that there is a critical initial slug length dividing the splitting process into two different regimes: flow dominated by the main channel and flow dominated by the branch channel. We systematically investigate the influence of key parameters on the splitting dynamics, including the inclination angle of the junction, the channel widths/apertures, and the dynamic contact angles. It is shown that the splitting ratio depends nonmonotonically on the relative width of the branch channel to the main channel. Furthermore, it is demonstrated that the dynamic contact angles have a profound impact on the splitting ratios and meniscus velocities. We highlight the important, previously overlooked role of velocity-dependent contact angle in the splitting dynamics. Ignoring the velocity dependence of contact angles can lead to overprediction of interface velocities by over an order of magnitude.
This work sheds light on the mechanisms of two-phase flow driven by gravity/buoyancy in fractured media. The improved understanding on the liquid splitting dynamics at a junction is also of direct relevance to applications of controlling bifurcation of liquid slugs/plugs in microfluidics. This work has important implications in hydrological and environmental applications, where the mechanisms of liquid splitting at channel intersections in a complex network need to be elucidated in order to mechanistically predict unsaturated flow and contaminant transport behavior in the subsurface. Splitting dynamics of liquid flow in channel networks composed of a series of junctions/intersections remains to be explored in future work.

Data Availability Statement
Experimental data and movies are available in the supporting information. All data can be accessed from a public repository (https://www.hydroshare.org/resource/36326719eabd4d39864d17c5a64c25ad/).