Solute Trapping and the Mechanisms of Non-Fickian Transport in Partially Saturated Porous Media

We study the upscaling of pore-scale solute transport in partially saturated porous media at different saturation degrees. The interaction between structural heterogeneity, phases distribution and small-scale flow dynamics induces complex flow patterns and broad probability distributions of flow, which control key aspects of transport, such as residence and arrival times, dispersion, and spatial solute distributions, as well as chemical reactions. A continuous-time random walk (CTRW) framework that integrates the processes of advection, diffusion

Recent numerical simulations showed that the simultaneous occurrence of preferential flow paths (of high fluid flow velocities) and stagnation zones (of low fluid flow velocities or immobile) give rise to broader probability distributions of fluid flow velocities than in more homogeneous cases (Triadis et al., 2019;Velásquez-Parra et al., 2022), which, in turn, enhance solute dispersion and increase or decrease chemical reactivity (Jiménez-Martínez et al., 2020;Nissan & Berkowitz, 2019;Sole-Mari & Fernàndez-Garcia, 2018).From a point of view of the applications, this may impact on soil and groundwater remediation actions (De Marsily, 1986;Moreno & Paster, 2017).
In this work, we analyze isothermal non-reactive transport of dissolved substances during horizontal steady water flow in variably saturated porous media using the continuous-time random walk (CTRW) framework.Under these conditions, water and air can be regarded as immiscible and incompressible.Steady water flow through partially saturated media occurs naturally when air bodies are entrapped within the media by capillary forces or in the case of constant low flow rates driven by gravity or suction (e.g., transpiration or evaporation), or at large enough distances from a periodical source (Ben-Noah & Friedman, 2019).In this case, the water and air distribution is stagnant, that is, the media properties vary spatially but are constant in time.
Travel time based approaches such as the CTRW framework account for the fact that the residence time of particles in low-velocity regions is longer than in the high-velocity regions and that the velocity varies in space on a characteristic length scale but not in time (De Josselin de Jong, 1958;Saffman, 1959;Scher & Lax, 1973).The CTRW approach models transport through spatio-temporal particle transition whose transition time distribution reflects the (broad) distribution of mass transfer times.This approach has been extensively used in different implementations for transport modeling in saturated media (Berkowitz et al., 2006;Delay et al., 2005;Noetinger et al., 2016).For transport in unsaturated media, CTRW models have been used as fitting approaches to adjust solute BTC (Bromly & Hinz, 2004;Bücker-Gittel et al., 2003;Cortis & Berkowitz, 2004;Zoia et al., 2010).Recently, the CTRW framework has been used for the interpretation and modeling of purely advective solute transport in partially saturated media (Aquino & Velásquez-Parra, 2022;Cortis & Berkowitz, 2004;Velásquez-Parra et al., 2022) following the approach of Dentz et al. (2016) for the quantification of the evolution of Lagrangian flow and transport properties in spatial random flows.Dentz et al. (2018) suggested a predictive CTRW framework integrating advection, diffusion, and trapping (in very low velocity and immobile regions) mechanisms for pore-scale transport in a saturated porous medium.Here, we adopt and adapt this framework to analyze and upscale solute transport in variably saturated porous media.
The paper is organized as follows.First, we introduce the experimental data and numerical models that are set as the basis for this work.Second, we present a step-by-step evolution of the integrated CTRW framework to predict macro-scale transport problems.Finally, we discuss on the relation between the physical model parameters for different saturation degrees and Péclet numbers (Pe).

Materials and Methods
In this section, we first describe the extraction of the saturation distributions corresponding to the experimental images from Jiménez-Martínez et al. (2017).Then, we pose the equations and boundary conditions governing pore-scale flow and transport.Finally, we describe the integrated CTRW framework for the interpretation of solute transport in partially saturated media.

Experimental Saturation Distribution and Image Processing
Media and phase distributions for different saturation degrees are obtained from the experimental work of Jiménez-Martínez et al. (2017).These authors conducted a set of multi-phase experiments of steady liquid (water and glycerol mixture) flow and solute transport in a hydrophilic millifluidic device also containing a stagnant gas (air) phase.Flow in the millifluidic device was controlled so that the phase distributions were not altered during the flow and transport experiments, that is, steady-state imbibing conditions were attained.The flow rates in the experimental setups correspond to a capillary number of about 3 ⋅ 10 −5 , which is considered a capillary dominated flow.Experimental images are processed using MATLAB R2021b to identify the grains, gas and liquid 10.1029/2022WR033613 3 of 15 regions and determine the effective flow domain and the effective saturation degree S e used for the numerical simulations and for the analysis below.The effective saturation degree S e is defined as the fraction of the porosity filled with the continuous and percolated liquid phase.The domain porosity in this study is fixed (ϕ = 0.72).
In this work we compare four effective saturation degrees, of S e = 1.00, 0.80, 0.73, and 0.65.Images of the medium at the different saturation degrees are presented in Figures 1a-1d in accordance.The connected liquid phase is in white, the unconnected liquid phase in red color.The gas phase and solid cylinders are in black.Here we treat the gas (air), solid and unconnected liquid phases similarly, as stagnant and inert media.The actual saturation degrees corresponding to the images in Figure 1 are S w = 1.00, 0.83, 0.77, and 0.71.
The image resolution is processed with a square pixel of 0.032 mm side length.The images do not include the entire original flow device and depict a flow domain of L = 105 mm (3,292 pixels) along the main flowing axis (x-axis) and width W = 70 mm (2,182 pixels).

Pore-Scale Flow and Solute Transport
Pore-scale flow and transport in the geometry within the phase distributions reported in the previous section are solved numerically using COMSOL Multiphysics (COMSO®, 2018).In the following, we describe the governing equations of pore-scale flow and transport and the corresponding boundary conditions.

Stokes Flow
The two-dimensional flow field is determined by solving the steady-state Stokes equation assuming incompressible laminar flow (1) where u [L/T] is the flow velocity, μ [M/(L•T)] the liquid viscosity, b [L] the device height.The term −12μu/b 2 in Equation 1 accounts for the viscous forces induced by the millifluidic device's top and bottom plates (Homsy, 1987).Constant inflow volumetric flux (q m ) was imposed on the left boundary, constant pressure at the outlet and no-slip boundary conditions at the liquid-solid and liquidgas interfaces.The value of the parameters for the numerical simulations are from Jiménez-Martínez et al. ( 2017) and summarized in Table 1.The linearity of the flow Equation 1 implies that viscosity has no impact on the flow velocities and transport but affects only the magnitude of fluid pressure.Also, both pressure and velocity (normalized by their mean values) are independent of the fluid properties.However, these parameters, and also the wettability (surface tension and contact angle), would affect the phase distribution in the experimental setup (Avraam & Payatakes, 1995), which is why they are reported here.
We determine the probability distribution of Eulerian flow speeds u e = |u|, denoted by p e (u), through areal sampling from the direct numerical flow simulations as where A i is the area of the ith mesh element, and u i is the flow speed at the Gaussian point of the ith mesh element.The indicator function   is one if the velocity u i is within the u k ≤ u i < u k + Δu k interval or zero otherwise.
where Δu k id the kth velocity bin.
The distributions    () and    () are defined accordingly.Furthermore, we consider the advective tortuosity χ, which measures the ratio of averaged trajectory length to the linear distance.χ accounts to the divergence of the streamline from the macroscopic flow direction.Thus, It can be calculated by the ratio of the mean flow speed   [L/T] and mean flow velocity (in the flow direction)   = ∕ (Koponen et al., 1996), where ϕ is the media's porosity Note that in this work, we treat the effect of the stagnant gas (air) phase on the liquid phase as equivalent to the solid phase.However, some theoretical works (Avraam & Payatakes, 1995;Kalaydjian, 1990;Whitaker, 1985) suggested that the gas-liquid interface may have a different effect than the solid-liquid, namely due to significant slip conditions at the gas-liquid interface.To analyze this effect, we compare two simulations with full-slip and no-slip conditions on the gas-liquid interface for the case S e = 0.65 and compute the respective probability density functions (pdfs) of velocity (shown in Figure S1a in Supporting Information S1) and BTC (shown in Figure S1b in Supporting Information S1).Both pdfs and BTCs are very similar, which is in accordance with results from others for example, (Guédon et al., 2019;Jiménez-Martínez et al., 2020;Markale et al., 2022;Triadis et al., 2019).This can be explained by the fact that the gas-liquid interface is restricted mainly to the junctions and dead-ends, where the velocities are very low.Therefore, and in the following, stagnant air will be treated in the same manner as the solid matrix (i.e., no-slip condition).It should be noted that for a more viscous non-wetting invading fluid (not the case of our evaluation), the assumption of no-slip condition may not hold (Lasseux & Valdés-Parada, 2022).

Solute Transport
The transient transport of a dissolved substance is simulated using the advection-diffusion equation where c [M/L 3 ] is the solute concentration, t [T] is time, and D [L 2 /T] is the molecular diffusion coefficient (in a free phase).Constant concentration c 0 is prescribed at the inlet, and D∇c ⋅ n = 0 at the outlet, where n is the unit vector normal to the outflow boundary.Zero flux boundaries are specified at the liquid-solid and liquid-gas The ratio between τ D and τ u defines the Péclet number Transport scenarios are uniquely defined in terms of Pe, that is, increasing D has the same effect on transport as decreasing flow rate.In the numerical simulations, we use three values for D: 1.5 ⋅ 10 −11 , 1.5 ⋅ 10 −10 , and 1.5 ⋅ 10 −9 m 2 /s; which are 0.1, 1, and 10 times the value used in Jiménez-Martínez et al. (2017).
Transport was analyzed in terms of the concentration BTCs at different control planes (x-direction).They were evaluated at each time step as the complementary normalized mean mass flux at the control plane In order to analyze the breakthrough behavior, we also consider the classical advection-dispersion approach.In this framework, the evolution of the average concentration   follows where   is the mean flow velocity, which is aligned with the x-direction of the coordinate system, D h [L 2 /T] is the hydrodynamic dispersion coefficient, θ e = S e ϕ is the effective water content.For a constant solute injection over the full width of the medium  ( (0, , ) = 0 ) the complementary cumulative BTC at a distance x from the source is given by (Lapidus & Amundson, 1952;Van Genuchten, 1982)

Integrated Lagrangian CTRW Framework
In this paper, we adopt the integrated CTRW framework of Dentz et al. (2018) to interpret and quantify solute BTC in partially saturated media.The integrated 1D Lagrangian CTRW framework models the transport of solute particles that move along a streamline with a random transition time (duration, τ [T]) over a characteristic length λ [L].In this framework the position x n and time t n of a solute particle (or parcel) along the mean flow direction after n + 1 steps are The spatial step ξ n [L] may be positive or negative, that is, down-or up-stream, depending on the advection and diffusion conditions.While particles transitions along their trajectories are made on a characteristic length (λ), the transition length |ξ n | = λ/χ accounts for the streamline (advective) tortuosity (χ).The transition time τ n depends on the pore-scale advection and diffusion mechanisms and is defined below.For example, for purely advective transport, it is given by the correlation length and the particle speed.
The breakthrough time τ L at a control plane located at a streamwise position x = L is given by 10.1029/2022WR033613 6 of 15 where N L is the number of random walk steps to reach the control plane.The BTC is given by where the angular brackets denote the average over all particles, and H(t) is the Heaviside step function.Details on the model derivation are given in Dentz et al. (2018).The model was implemented here in MATLAB.In the following, we summarize the key elements of the model.

Purely Advective Transport-CTRW a
We first consider the representation of purely advective transport in this framework.In this case, the transition time τ n is equal to the advection time over the distance λ, where u n is the particle speed at the nth random walk step.The particle speeds are modeled as independent identically distributed random variables whose distribution is denoted by p s (u).The transition length λ denotes the length scale at which subsequent speeds decorrelate, or in other words, the speed correlation length.The distribution p s (u) is related to the probability distribution p e (u) of Eulerian flow speeds through the speed-weighting relationship (Dentz et al., 2016) This relation can be rigorously derived based on volume conservation using the Reynolds theorem.It can be qualitatively understood as follows.The Eulerian speed pdf is obtained by volumetric sampling of the Eulerian speed, that is, it gives more weight to low than to high speeds because low speeds have a wider streamtube.The equidistant sampling of the Lagrangian speed weights all velocities equally, which is compensated for by the flux-weighting relationship.Note that at the extreme case of a zero velocity the weighted probability is also zero, which can be explained by the fact that a particle can not advectively enter an area where there is no velocity.Relation 16 has been confirmed numerically for pore and Darcy scale flows (Hakoun et al., 2019;Puyguiraud et al., 2019).
The transition length λ is related to the medium geometry, that is, the non-flowing matrix and the flow organization, and is related to the scale of heterogeneity.For example, for saturated homogeneous media, λ is of the order of the average pore size (Puyguiraud et al., 2019).For unsaturated media, λ is expected to increase with the reduction of the effective phase saturation (S e ) because of the emergence of irregular air clusters and preferential flow.
The arrival time due to advection only is given by where the superscript a denotes advective.The number of steps ) .

Advective-Diffusive Transport-CTRW m
For advective-diffusive transport, the transition time τ n at step n is defined as the harmonic sum of the advective ) and diffusive (λ 2 /2D) transition times (Dentz et al., 2018) If u n is much larger than D/λ, the transition time is approximately equal to the advective transition time.If u n ≪ D/λ it is approximately equal to the diffusion time τ D .Note that as an extension of the purely advective CTRW, this approach quantifies the impact of diffusion on transport along a streamline.It does not account for diffusion into stagnant regions perpendicular to the streamlines.This effect is quantified in the next section.
The arrival time due to advective-diffusive transport is given by where the superscript m denotes mobile, in contrast to trapped, which we consider in the next section.Also for advective-diffusive transport, the mean transition time  ⟨   ⟩ = ∕ because in average, diffusive transitions cancel out.However, the number of steps to reach the control plane     increases with upstream steps due to diffusion.Using Equations 18 and 19, we obtain where   ( ′ ) is the pdf of   ′ = ∕ .For the flow and transport scenarios under consideration here, Pe ≥ 10 2 , the relative error of using     =    approximation is less than 4%.If one considers a soil with water content of 0.3 and a characteristic length of 1 cm, than a Pe of 100, for Chloride (D = 10 −9 m 2 /s), means a flux of about 1 cm/hr which is not very uncommon in terms of rain intensity, infiltration in recharge basins, or even during irrigation, locally close enough to the source.

Advective-Diffusive Transport and Trapping-CTRW i
In partially saturated media, the flow domain is characterized by relatively wide regions of low velocity such as dead-ends or narrow water films near the grains of an air invaded pore.During transport in complex media, particles may get entrapped in these regions, causing long tailed BTC characteristic of non-Fickian transport (Gjetvaj et al., 2015;Liu & Kitanidis, 2012).Following Dentz et al. (2018), the trapping events are assumed to occur at a constant frequency γ [1/T].This implies that the number n t of trapping events that occur during a time t is Poisson distributed The average number of trapping events during the mobile arrival time     is then given by  ⟨    ⟩ =    .The trapping time τ f , that is, the amount of time a particle gets entrapped in the immobile regions can be determined as the escape time distribution of a diffusing particle.It is approximated here by the exponential distribution (Dentz et al., 2018) The ratio β of the immobile and mobile volume fractions is given by (Dentz et al., 2018) The volume fraction of immobile regions is expected to decrease the saturation degree (S w and S e ).
With these preparations, we can now write the arrival time τ L due to advective-diffusive transport and diffusion into dead-end zones can be written as The BTC is then given by expression 14.The long-time behavior of the BTC is given by the exponential (Dentz et al., 2018)

Results and Discussion
As discussed in the previous section, we use the micro-scale flow characteristics to construct a predictive macro-scale transport model.In addition to the velocity pdfs, the media are described by the characteristic length λ, the trapping frequency γ and the mean trapping time 〈τ f 〉.Notice that γ and 〈τ f 〉 are functions of both S e and D while λ is a function of S e only.In the following, we first study the impact of the saturation degree on the velocity distribution and advective tortuosity.Then, we discuss the solute BTC at different saturation degrees in the light of the upscaled integrate CTRW i model for different Péclet numbers.The length scale λ is determined by adjusting the purely advective CTRW a model to the intermediate time behaviors of the BTCs at high Pe.The trapping rate and trapping time are determined from the tail behavior using relation 25.Finally, we discuss the dependence of the model parameters on saturation and Péclet number.

Velocity Statistics
In this section, we discuss the results of steady liquid flow for different saturation degrees.The effect of the stagnant gas (air) phase on the spatial distribution of liquid flow velocity is not trivial.In fully saturated granular media the liquid phase flows in most of porous domain and preferential pathways prevails mainly in cases of correlated large pores (e.g., cracks).However, in partially saturated media, there are several counteracting mechanisms.On the one hand, capillary forces drive the water into the smaller pores, and in turn, air is driven to the larger pores, which forces water to flow in narrower paths (of greater resistance).This effect is expected to increase χ and to decrease water flow velocity variability.On the other hand, stagnant air generates dead-ends, which are areas of very low water velocity, thus significantly reducing the effective cross-section for the liquid flow, causing high velocity in the remaining paths, thus decreasing χ and increasing water velocities variability.

Velocity Distribution
The different variability of liquid flow velocities for two saturation degrees is illustrated in Figure 2. Results showing the relatively uniform water flow velocity distribution in saturated (Figure 2a), and the velocities variability and preferential pathways of the partially saturated (Figure 2b) media, suggest that the formation of deadends has a much more significant effect than the closing of preferential pathways, even for low Ca (3 ⋅ 10 −5 ).In this context, it should be mentioned that the average Ca may not be a good quantifier of the micro-scale interplay between the viscous and capillary forces, that is, that viscous forces may locally dominate in the pore-scale preferential pathways (Tang et al., 2019).The fact that the pore-scale velocities are distributed implies that the local capillary number is distributed itself.Its distribution is in fact similar to the water velocity distribution.The flow pattern and therefore the water velocity distribution is a result of the interaction between viscous and capillary forces, and therefore depends on the driving force.These effects will be studied elsewhere.
A noticeable velocity variability can be observed at the intra-pore scale (zoomed insets in Figure 2) for both the fully and partially saturated cases.This flow velocity variability for different saturation degrees are manifested  2022), for these same images.They separately demonstrate the fundamental effect of the trapped air on the flow behaviors.Both the fraction of small velocities and the magnitude of the maximum velocity increase with the reduction of S e (Figure 3a).For the saturated case, the pdf of u x has a clear bias toward positive values (Figure 3b).In contrast, the fraction of the counter direction flow is much larger in the partially saturated media.
Interestingly, the magnitude of the maximum velocity in x-direction is of the same order of magnitude as in y-direction (Figure 3c).The symmetry of the flow in the perpendicular (y-axis) velocity increases with the increase in saturation degree, that is, the absolute value of the sum of the velocities trajectories on the y-axis (|Σu y |) decreases with increasing S e (not presented).In these cases, the sum of u y is very small (two order of magnitude smaller) compared to   .However, this effect may suggests that a larger device (representative elementary volume [REV]) may be needed to account for lower saturation degrees.This statement is also supported by the effect of S e on λ as discussed in the next subsection.

Tortuosity
The dependence of the advective tortuosity χ, given by Equation 4, on the effective saturation degree S e is depicted in Figure 4. We observe that χ decreases with increasing S e .This dependence can be approximated by the power-law   = 1.1 −0.6  .As mentioned above, the effect of entrapped air on χ is dichotomous, on one hand the stagnant air phase decreases the effective cross-section, while on the other hand it increases the structural tortuosity.The latter, also termed diffusive tortuosity (χ D ) (Ghanbarian et al., 2013), is defined as the ratio between the product of the molecular diffusion D and effective water content θ e and the effective molecular diffusion coefficient D e (i.e., by χ D = θ e D/D e ).D e can be evaluated by Fick's law as the steady diffusive flux obtained under a unit concentration gradient (for no flow conditions in Equation 5, where D h reduces to D e ).The effect of S e on the structural tortuosity, is commonly described with a power-law (i.e., Archie's law (Archie, 1942)).However, the decrease of χ D with increase of S e is usually much steeper than the  ∼  −0.6  -dependence found here for the advective tortuosity χ (Friedman, 2005), especially in two-dimensional flow domain, in which the phase connectivity is more sensitive to the reduction in the saturation degree.Here, the effective diffusion coefficient D e is calculated from the steady diffusive flux using the numerical simulations with u m = 0. We find that   ∼  −2.8  (R 2 = 0.973, see Figure S2 in Supporting Information S1), that is, it decreases much faster with increasing S e than χ.This stronger dependence of χ D on S e compared to χ can be explained by the smoothness of advective particle trajectories, which align with streamlines of the flow field, in contrast to the irregular trajectories of a particle that diffuses through the pore space.The large structural tortuosity compared to χ and the relatively moderate effect of S e on χ compared to its effect on χ D (see Figure S3 in Supporting Information S1) are indicators of the significance of preferential flow paths for advective tortuosity.

Concentration Distribution and Breakthrough Curves
The effect of the saturation degree on the concentration distribution is depicted in Figure 5.The transport through the saturated medium demonstrates a piston type flow with a sharp concentration gradient in the main flow direction.The partially saturated cases, on the other hand, show sharp concentration gradients perpendicular to the flow path and organization of the displacement front into a lamellar pattern (Jiménez-Martínez et al., 2015).The relatively uniform intrapore concentrations in contrast to the strong velocity gradients (zoomed insets in Figure 2), suggest strong mixing induced by diffusion.For the saturated case (Figure 5a), the concentration near the device vertical walls (y = 0 and y = W) are significantly lower.This effect is because of the correlation of low velocities along and in close proximity to the walls.The wall effect is significant only for the saturated case (Figure 5a) due to its relatively uniform velocities (Figure 2a) so that the correlation of low velocities along the walls has a noticeable effect.The effect of the walls reduces with increasing saturation degree.This is because of the less uniform flow in the entire domain and also because the air clusters near the walls force the water to bypass through the inner sections of the device.The impact of the walls on the BTC can be eliminated by calculating solute transport on a frame with a distance Δy [L] from the walls.The value of Δy can be determined for example, from the distribution of the longitudinally averaged speed, denoted by 〈u〉 x , across the domain width.Figure S5a in Supporting Information S1 shows 〈u〉 x for the different S e , while Figure S5b in Supporting Information S1 depicts the BTC of the saturated case with different Δy.These results suggests that for the S e = 1, using a Δy = 6 mm is sufficient to avoid tailing caused by the walls, and that for the unsaturated cases there is no need of using a frame.
After these preliminary remarks, we now analyze the solute BTCs and their representation by the CTRW models discussed in Section 2.3.Figure 6 presents a comparison of the results for the S e = 0.73, D = 1.5 ⋅ 10 −10 m 2 /s case to the fitted analytical advection-dispersion solution (Equation 11), the purely advective CTRW a framework (Equation 13) and the advection-diffusion (mobile) CTRW m (Equation 18).The value of D h /θ e (1.5 ⋅ 10 −7 m 2 /s) in the analytical solution of the advection-dispersion equation was fitted to a least square difference.The advection-dispersion model (Figure 6 dash-dotted line) cannot correctly describe the BTC.It underestimates the long-time breakthrough and overestimates concentration at intermediate times.
The purely advective CTRW a accurately estimates the early and intermediate times breakthrough.However, it displays a heavier tailing than the data from the direct numerical simulations.This is expected because in the purely advective CTRW a , particles on very low velocity streamlines can experience transition times much larger   1. Simulations were conducted as detailed in Section 2.2.2.than the limit imposed by diffusive motion (Dentz et al., 2004).As outlined above, the correlation length λ is estimated by comparing the CTRW a BTCs to the numerical data of the advection-dominated scenario with the largesst Pe (Pe = 395, 1.17 ⋅ 10 4 , 1.60 ⋅ 10 4 , and 7.16 ⋅ 10 4 , for the S e = 1.00, 0.80, 0.73, and 0.65, in accordance, corresponding to D = 1.5 ⋅ 10 −11 m 2 /s in the different scenarios).
On the other hand, the advection-diffusion CTRW m (Figure 6 dashed line) underestimates the tailing.This is because the CTRW m accounts for diffusion along trajectories but not for the trapping of particles in dead-ends and immobile zones.The fact that the actual BTC tailing is bounded between the CTRW a and the CTRW m predictions may have practical implications on risk assessment, providing a relatively narrow range of expected residual times even without evaluating the immobile parameters.
The integrated CTRW i framework can be used for accurately modeling the arrival times.This framework is able to model the solute arrival times across different effective saturation degrees S e and diffusion coefficients D over more than six orders of magnitude as shown in Figure 7.The trapping rates and times are adjusted from the exponential tails of the BTCs from the direct numerical simulations as detailed in Section 2.3.3.In addition to the effect on the tailing, trapping in immobile zones, slows down transport in the mean flow direction.For times larger than the mean trapping time, this effect can be quantified by a retardation factor.
Interestingly, the solute tailing and early arrival of the case with saturation degree S e = 0.65 (Figure 7c) is much more pronounced than that of the S e = 0.80 (Figure 7a) or S e = 0.73 (Figure 7b) cases, which are only slightly different.This point may reflect on the significant effect of the stagnant air-phase connectivity on the formation and extent of dead-ends regions in the flowing phase.

Scalability of the CTRW
Figure 8 presents a comparison between the direct numerical simulations and the CTRW i models BTC, at different flow domain cross-sections (location of the control plane).The different CTRW i accounts for the slightly different velocity pdfs and χ of the different cross-sections but uses the same model parameters.The excellent fit, even at half the device length (Figure 8a), which stands for about 4.3λ, indicates that even the complex unsaturated flow domain decorrelates on a characteristic length.Using shorter domains, for example, 0.25 L (2.2λ), showed deviation between the numerical and the CTRW i models BTCs (not presented), suggesting that at this short lengths, the domain characteristics varies significantly from the full device and cannot stand as a REV.These results also reflects on the robustness of the CTRW i framework and the physical nature of the model parameters, not being a fully case sensitive, strictly empirical values.The parameters values as a function of saturation and diffusion coefficients are discussed in the next section.

Correlation Length and Trapping Parameters
In this section, we discuss the physical meaning of the parameters of the upscaled CTRW i model and their dependence on the system parameters in terms of saturation and Pe.For the advection-diffusion CTRW (broken line, Equation 18), the value of the molecular diffusion coefficient (D = 1.5•10 −10 m 2 /s) was taken from Jiménez-Martínez et al. (2017).We used 10 8 number of particles (N p ) for the CTRW models.

Correlation Length
The transition length λ of the system was estimated by fitting the BTC from the numerical simulations with the highest Pe (lowest D) to the results of the purely advective CTRW a model Equation 13 using the pdfs of flow speed and χ values reported in Section 3.1.Figure 9 shows that λ decreases as saturation increases.The results can be approximated by the exponential relation λ = 259 exp(−13S e ).This strong dependence can be partially explained by the shape and size of the stagnant air clusters which are mostly in form of large ganglia with a very few number of small bubbles.This suggests that correlations are case sensitive, and that for different conditions, for example, different fluid properties (e.g., wettability), different flow dynamics, a vertical orientation, or in a 3D flow domain, the relation would be significantly different.However, the general perception of strongly negative correlation will remain.Moreover, the strong increase in λ with the decrease in S e , along with the increasing asymmetry in the velocities on the y-axis (discussed above), suggests that the minimal REV increases with the decrease in S e .

Trapping Parameters
The Pe is by definition (Equation 8) inversely proportional to D and proportional to the mean flow velocity.In turn, the latter (for a given discharge rate, q m ) is related to τ and u, which depends on S e .Moreover, the Pe is proportional to λ, which is exponentially related to S e (Figure 9).We evaluate flow and transport, for different Pe values in the range of 10 2 -10 5 corresponding to different D and S e .
The effect of Pe on the immobile characteristics is depicted in Figure 10.The trapping frequency γ decreases with increasing Pe (Figure 10a).This is expected, because the probability per time for a particle to diffuse into an immobile zone at low Pe (diffusion dominated) is larger than for large Pe (advection-dominated).The single effect of D on γ is best fitted with a power-law exponent of 0.87 (see Figure S4 in Supporting Information S1).
However, the theoretical linear relation (exponent of 1), inferred from the explanation above, can also provide a close approximation (see Figure S4 in Supporting Information S1).Speculatively, this slight decrease from linearity may be explained by the effect of diffusion on the transition time in low velocity increments.
The normalized mean trapping time (〈τ f 〉/τ m ) increases approximately linearly with Pe (Figure 10b).If we assume that the trapping time is related to D with some characteristic trapping length (l 0 ), that is, that  ⟨ ⟩ =  2 0 ∕ (Dentz al., 2018) then, the ratio 〈τ f 〉/τ D is related to  (0∕) 2 , that is, independent of D. The fact that 〈τ f 〉/τ m is proportional to Pe (Figure 10b), implies that l 0 is proportional to (and about 5.5 times larger than) λ.This correlation may be explained by the fact that both the transition length and the fraction of stagnant areas are governed by the phase saturation degree and by its configuration, or in other words, the characteristic length scale of fluid clusters dominates both the velocity correlation and size of immobile regions.This  merit suggests that knowledge about l 0 can be used to evaluate both λ and 〈τ f 〉/τ m .Thus relation between l 0 and the air clusters shape and topology may prove to be a valuable future research path.
Figure 10c depicts the effect of Pe on β.As expected, the volume fraction β of immobile regions is strictly a parameter of the medium, that is, independent of D and increasing with S e .Interestingly, the β trends of the different D (black lines, Figure 10c) demonstrate a logarithmic relation with a similar slope.This is in accordance with other works (Markale et al., 2022) that found a non-linear variation in the stagnant volume fraction with the saturation degree.The single factor effect of D on τ f /τ m is well approximated with a power-law with a slope between −0.91 and −0.83 (for the different S e , not presented).

Conclusions
Water flow and solute transport in variably saturated media is of major concern in agricultural application (e.g., "fertigation") and environmental risk assessments from contaminants in the critical zone.In this paper, we propose a new framework to quantify and upscale anomalous transport in partially saturated media based on the approach of Dentz et al. (2018).Entrapped air significantly affects the water flow velocities.It dramatically increases the fraction of very low velocities in dead-end areas.Moreover, the trapped air enhances the velocity in directions perpendicular to and opposite to the main flow.Air clusters increase both velocities variability and advective tortuosity.As a result, the presence of trapped air clusters leads to non-Fickian solute transport that the classical advection-dispersion equation cannot explain.CTRW is a powerful tool to upscale heterogeneous micro-scale flow and transport features to macro-scale, accurately predicting the arrival times of contaminants and at trace amounts.Comparing the physics-based CTRW approach to direct numerical simulations and experimental results is used to derive up-scaling rules and correlations.Trapped air clusters increase the characteristic velocity correlation length (λ), the advective tortuosity, and the fraction of immobile area (both l 0 and β).The Pe was found to be strongly correlated to the trapping frequency (γ) and to the normalized mean trapping time (〈τ f 〉/τ m ).These findings provide new insights into the understanding of non-Fickian transport in partially saturated media, and its relation to medium geometry, saturation degree, and flow conditions.
The present analysis presents a scalable framework that is based on the pore-scale physics of particle motion and trapping.The model parameters are physics-based, and depend on the distribution and shape of the mobile and immobile regions which in turn depend on the flow dynamics, fluid properties (e.g., surface tension and viscosity ratio), and pore geometry and topology.Regarding the correlations and effects of S e and Pe on the transport and trapping mechanisms, we expect them to show qualitatively the same behavior in two-and three-dimensional flow domains.
In conclusion, we have shown that solute transport under different saturation degrees can be upscaled by the proposed integrated CTRW model, which connects Eulerian flow properties and medium geometry to large scale dispersion.A future challenge is how relate fluid properties, flow conditions and medium structure to the Eulerian flow properties and distribution of low flow regions.

Figure 1 .
Figure 1.Images of the porous medium at saturation (S w ) and effective saturation (S e ) degrees of (a) S w = S e = 1.00,(b) S w = 0.83, S e = 0.80, (c) S w = 0.77, S e = 0.73, and (d) S w = 0.71, S e = 0.65.In black is the solid and gas phases, in white the connected liquid phase, and in red is the unconnected liquid phase that set the difference between S w and S e .
Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022WR033613by Paul Scherrer Institut PSI, Wiley Online Library on [09/03/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License interfaces.The initial concentration in the domain is equal to zero.Next, we define two characteristic transport time scales.The advection time τ u [T] measures the mean advection time over the characteristic length scale λ [L] /T] is the mean Eulerian speed.The length scale λ is discussed below (in Section 2.3.1).Furthermore, we define the characteristic diffusion time τ D [T] over the distance λ as,

2023, 2 ,
Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022WR033613by Paul Scherrer Institut PSI, Wiley Online Library on [09/03/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License This relation allows to estimate the trapping rate γ and the mean trapping time 〈τ f 〉 from the long-time tails of the BTC.The dependence of γ and 〈τ f 〉 on S e and D are discussed in the next section.

Figure 2 .
Figure 2. Flow velocity field for (a) fully saturated and (b) effective saturation degree of 0.73.The velocities were calculated using COMSOL Multiphysics with the parameters described in the numerical simulation section andTable 1.

Figure 3 .
Figure 3. (a) Probability density function of flow velocity p(u), (b) in the main flow direction p(u x ), and (c) in the transverse direction p(u y ), for different effective saturation degrees S e : fully saturated (blue), 0.80 (red), 0.73 (orange), and 0.65 (purple).

Figure 4 .
Figure 4. Effect of the effective saturation degree (S e ) on the advective tortuosity (χ).

Figure 6 .
Figure 6.Breakthrough curves of the different models for S e = 0.73.The direct numerical simulation is indicated in circles.For the analytical solution (line-dot), the hydrodynamic dispersion coefficient (D h /θ e = 1.5 10 −7 m 2 /s) was fitted (least squares).The purely advective CTRW a (dotted line) is presented in Equation13.For the continuous-time random walk (CTRW) frameworks, the characteristic length value (λ = 0.014 m) was fitted against simulations with lower molecular diffusion coefficient (0.1 D).For the advection-diffusion CTRW (broken line, Equation18), the value of the molecular diffusion coefficient (D = 1.5•10 −10 m 2 /s) was taken from Jiménez-Martínez et al. (2017).We used 10 8 number of particles (N p ) for the CTRW models.

Figure 9 .
Figure 9.Effect of the effective saturation degree (S e ) on the correlation length (λ).