Mean flow direction modulates non-Fickian transport in a heterogeneous alluvial aquifer-aquitard system

.

represent non-Fickian (i.e., anomalous) transport under a range of geologic heterogeneity and transient flow boundary conditions encountered on large temporal and spatial scales.However, modeling regional nonpoint transport remains a challenging task due to ubiquitous non-Fickian spreading which is difficult to model at regional scales with existing approaches, especially when the boundary conditions shift from steady state to transient flow (Guo et al., 2020;LaBolle & Fogg, 2001).
Numerous studies show that anomalous transport in geologic media depends strongly on spatially heterogeneous hydraulic conductivity (Dagan & Neuman, 2005;De Marsily et al., 2005;Gelhar & Gelhar, 1993), which gives rise to the velocity field (Berkowitz et al., 2006).This study concerns Darcy scale processes in a typical alluvial aquifer-aquitard complex, wherein variability in the hydraulic conductivity (K) of sediments can range over 4-7 orders of magnitude (Fetter, 2001;Fogg & Zhang, 2016).The fluvial deposition processes that form alluvial aquifer-aquitard systems over millennia result in correlated fields (Galloway & Hobday, 1996), with interconnecting, high-K, coarse sand, and gravel lenses (Fogg, 1986) that can percolate in all spatial dimensions (Fogg et al., 2000), even at substantially low volume fractions (e.g., 13% as noted in 3D percolation studies by Harter (2005)).These high-K preferential flow paths facilitate anomalous early time arrival, and are embedded within low-K, fine-grained clays and silts, which typically comprise the majority of alluvial aquifer-aquitard volume fractions (Galloway & Hobday, 1996), and contribute to mass hold back and anomalous late-time solute tailing (Benson et al., 2000;Berkowitz et al., 2006;Haggerty et al., 2000;LaBolle & Fogg, 2001).Moreover, in typical alluvial heterogeneous material, where the horizontal length scales of K variability are often 2-3 orders of magnitude or more larger than vertical length scales of K variability, .Zhang et al. (2006) found that vertical transport is increasingly Fickian under more anisotropic conditions.
The interaction between spatial heterogeneity, hydraulics, and diffusion determines transport processes.However, anomalous transport studies commonly investigate simple hydraulic boundary conditions characterized by a unidirectional gradient (Benson et al., 2001;Edery et al., 2016Edery et al., , 2014;;Harvey & Gorelick, 2000;Henri & Fernàndez-Garcia, 2015;Julian et al., 2001;LeBlanc et al., 1991;Willmann et al., 2008;Zhang & Brusseau, 1999), which is typical of natural, ambient groundwater flow (e.g., a horizontal hydraulic gradient on the order of 10 −3 ; Fetter, 2001).Many of these studies hold hydraulic boundary conditions constant in order to investigate how changing the spatial heterogeneity of the domain impacts solute transport.However, because of horizontal stratification of coarse-grained and fine-grained facies typical of many clastic sedimentary environments and hence significant vertical anisotropy of effective hydraulic conductivity (K), the vertical hydraulic gradients commonly exceed horizontal gradients by a factor of 10 1 − 10 2 (e.g., Belitz et al., 1993;Fogg, 1986).In even moderately developed groundwater basins, and especially in agriculturally intensive ones, high rates of seasonal pumping and recharge from irrigation can further increase vertical hydraulic gradients (Gailey, 2018;Styles & Burt, 1999).Yet, only a few modeling studies investigate non-Fickian transport under such large shifts in regional hydraulic gradient, and these studies tend to emphasize local scale breakthrough measured at extraction wells or control planes.Matheron and De Marsily (1980) demonstrated in a 2D model with limited spatial heterogeneity that flow parallel to the bedding direction leads to preasymptotic anomalous transport, even at very long time scales because the solute travels along an increasingly nonstationary medium.Henri and Harter (2019) demonstrated that well design characteristics (e.g., well depth, pumping rate, and screen length) have a significant impact on the uncertainty of late arrival times generated by heterogeneity.Similarly, Libera et al. (2017) demonstrated that temporally PAULOO ET AL.

10.1029/2020WR028655
variable pumping rates strongly affected contaminant breakthrough curves at the pumping well.Guo et al. (2019) investigated anomalous transport arising from the transition from predominately horizontal to vertical flow, and found that internal hydraulic gradients between the high-K and low-K materials strongly affect mass transfer between hydrofacies.Furthermore (Guo et al., 2020), showed that an upscaled multirate mass transfer model failed to capture tailing effects when the mean flow direction changed.Thus, prior work suggests that hydraulic boundary condition transience can strongly impact non-Fickian transport yet we still lack studies that demonstrate how shifting hydraulic gradients (as opposed to the unidirectional gradient that most studies explore) and geologic heterogeneity may interact to influence anomalous transport observed in heterogeneous aquifer-aquitard complexes.
To the best of our knowledge, no studies have investigated the impact of nonhorizontal and nonvertical mean flow direction on non-Fickian transport within a 3D, kilometer-scale, alluvial aquifer-aquitard system.Such a question bears importantly on modeling macroscopic nonpoint source transport in large alluvial groundwater systems, and may shed light on how to build better upscaled transport models that perform well across a range of transient hydraulic gradient conditions.
In this study, we simulate groundwater flow and random-walk particle tracking in a detailed hydrofacies model of a typical alluvial aquifer-aquitard complex (with high vertical anisotropy in K) subject to an increasingly high vertical to horizontal hydraulic gradient ratio (VHGR), which is an indicator of the mean flow direction.In other words, if we increase the vertical hydraulic gradients in a system by a factor of 10, 100, or 1,000 (e.g., due to natural geologic variation in vertical anisotropy of K, or during the transition from a fallow winter season to a summer growing season characterized by pumping and recharge), one might expect some change in the non-Fickian transport, but by how much, and by what mechanisms?Changes in anomalous nonpoint source transport are compared between differing steady-state VHGR scenarios which collectively represent snapshots along a transition from predominantly horizontal (low VHGR) to vertical flow (high VHGR).We examine tailing, spatial variance of the nonpoint source plume, and the spatial distribution of final particle position to quantify and explain differences in non-Fickian spreading.Results are discussed in the context of the transition from advection to diffusion dominant systems, and related to the hydrogeology of the study site.Our findings support a conceptualization of nonpoint source contaminant migration in 3D, kilometer-scale, alluvial aquifer-aquitard systems characterized by oscillating patterns of non-Fickian transport driven by predominantly vertical flow alternating with predominantly horizontal flow, owing to seasonal changes in pumping and recharge.We explore how this perspective may inform the development of upscaled regional groundwater quality management models (Fogg & LaBolle, 2006).

Methods
First, a heterogeneous, hydrofacies-based geostatistical model of the Kings River alluvial fan in California's Central Valley was developed.Next, four VHGR scenarios of increasingly vertical groundwater flow were simulated, and used to solve the particle transport of a conservative solute.Lastly, the resulting 3D Lagrangian particle trajectories were analyzed in terms of their tailing, second spatial moment, and spatial distribution of final particle position.

Heterogeneous Hydrofacies Model and Study Site
The three-dimensional, spatially correlated hydraulic conductivity domain (Figure 1) was simulated with Transition Probability Geostatistical Software (T-PROGS), which uses the transition probability approach described by Carle and Fogg (1996), Carle andFogg (1997), andCarle (1999).Details of the facies model are provided in Weissmann et al. (1999), Weissmann and Fogg (1999), and Weissmann et al. (2002c, 2004), and a brief summary is provided here.
A subsection of the Central Valley aquifer (California, USA) in the Kings River fan, was chosen as a study site representative of a typical alluvial aquifer for its high vertical anisotropy (K h /K v = 238), and its validated performance in simulating both flow (hydraulic head) and transport (groundwater age distribution) (Weissmann et al., 2002a).The hydrofacies model was conditioned by analyzing more than 40,000 surface and shallow subsurface (<2 m) features (e.g., borings, excavations, road cuts, and natural exposures) (Huntington, 1980), ∼150 m of core and geophysical logs (Burow et al., 1997), ∼1,000 m of continuous core (Harter PAULOO ET AL. 10.1029/2020WR028655 et al., 2005), and numerous drillers' logs from water wells in the study area (Weissmann & Fogg, 1999).These analyses defined a conceptual geologic model and the T-PROGS model parameters.
The Kings River alluvial fan is composed of four sequences created during Pleistocene glacial cycles and a fifth basal depositional unit prior to major glaciation events, each separated by laterally extensive paleosols that formed during inter-glacial periods of pedogenesis (Weissmann et al., 2002c), which, owing to their relatively low hydraulic conductivity, act as aquitards.Notably, one unit represents an incised valley fill of coarse-grained sediment in the back-to-front y direction of the domain, created by downcutting rivers during high ice-stages (low sea-level) then filled with glacial outwash as glaciers retreated and sea-level rose.Key parameters derived from these analyses that constrained the T-PROGS model include the hydrofacies proportions, mean lengths in xyz, and the embedded transition probabilities in xyz, which reflect the cross-correlation between facies and thus the probability that one facies contacts another in a given direction (Carle & Fogg, 1996).Next, Markov chain models corresponding to each unit were developed and combined to represent the highly heterogeneous stratigraphy: overlapping alluvial fans, interspersed with laterally extensive paleosols, and a cross-cutting incised valley fill.Weissmann et al. (2004) provides a detailed overview of the development and combination of these models into a single conductivity field.
The model dimensions in xyz are 12.6 km × 15 km × 100.5 m, with discretization Δx = Δy = 200 m, and Δz = 0.5 m.Spatially variable K estimated from pumping tests, slug tests, core measurements, and literature values for similar lithologies (Burow, 1999;Weissmann et al., 1999) were assigned to each of the model cells by their simulated hydrofacies, which include gravel, sand, sandy mud, muddy sand, and paleosols (Table 1).
We only use one realization of the hydrofacies model (Figure 1) because previous studies demonstrate that the statistics of particle paths, and hence non-Fickian transport behavior do not appreciably differ between realizations due to good connectivity of aquifer facies (Guo et al., 2019;Weissmann et al., 2002a).PAULOO ET AL.

Groundwater Flow and Transport Simulation
Steady-state, groundwater flow in the geostatistical model grid (Section 2.1) is simulated in MODFLOW (Harbaugh et al., 2000;McDonald & Harbaugh, 1988).The 3D groundwater flow equation used to compute long-term, average confined flow is where K xx , K yy , and K zz [LT −1 ] are the hydraulic conductivities along the x, y, and z coordinate axes; h [L] is the potentiometric head; W [T −1 ] is the volumetric flux per unit volume of sources (W > 0) or sinks (W < 0); and t [T] is time.
The water table is represented on the top of the domain with a constant head boundary condition, which is the same across all scenarios.Vertical flow out of the bottom of the model and into deeper groundwater, due to a combination of pumping at depth and recharge from above, is represented by another constant head boundary condition on the bottom of the domain, and is altered depending on the VHGR scenario in question (discussed in Section 2.3).Two more constant head boundary conditions on the front and back faces of the model maintain an ambient regional horizontal groundwater flow gradient of 0.001 along the y direction (Figure 1) of the model.Head values are spatially distributed across all active faces such that the mean flow direction is uniform across the entire domain (Figure S1).No flow boundaries are specified along the right and left faces of the domain (along the x direction).
Solute transport is simulated by the random-walk particle tracking method (LaBolle et al., 1998(LaBolle et al., , 2000) ) applied in the code RW3D (Fernàndez-Garcia et al., 2005;Henri & Fernàndez-Garcia, 2014, 2015;Salamon et al., 2006) that solves the advection-dispersion equation: where c is the concentration of the dissolved solute [ML −3 ], θ is the effective porosity [−], and D is the dispersion tensor, defined as The source term is a discrete, instantaneous pulse injection of 10,000 particles at 10 m below the water table occurring across the entire model domain (Figure 1).To simulate the spatial variability of Contaminant mass arrival at the water table, flux-weighted injection is used to initialize particle motion.Simulations are run for a sufficient number of time steps such that all particles exit the domain (Table S2), and that the mean Lagrangian hydrofacies proportions (i.e., the proportion of particles in each hydrofacies) converge on the actual hydrofacies proportions over time (Figure S2).The Lagrangian particle tracking methods used in this study can be very computationally intensive, especially for the case of a regional-scale nonpoint source term, where the size of the domain (and hence the number of particles used to sample it) is large.Adding multiple injections of particles was computationally unrealistic, and the steady-state results here obtained PAULOO ET AL.

Hydraulic Gradient Scenarios
The aim of this study is to investigate the influence of mean flow direction on non-Fickian transport, thus we compare systems under different vertical to horizontal gradient ratio (VHGR) where i v is the vertical gradient ∂h/∂z and i h is the horizontal gradient ∂h/∂y (Figure 1).Note that these are mean gradients which will depart from mean values locally.Increasing VHGR corresponds to increasingly vertical flow.Four scenarios are simulated, in which VHGR is varied over four orders of magnitude by fixing i h = 0.001 (representing ambient horizontal groundwater flow from mountain front recharge), and adjusting i v : • VHGR = 0.1 (i v = 0.0001): primarily horizontal flow with minimal recharge • VHGR = 1 (i v = 0.001): equal horizontal and vertical gradients • VHGR = 10 (i v = 0.01): substantial component of vertical groundwater flow • VHGR = 100 (i v = 0.1): dominantly vertical flow Within-facies heterogeneity (e.g., stratification of coarse and fine facies) can create strong regional vertical anisotropy values on the order of 100:1 or larger (Freeze & Cherry, 1979).In our study site, permeameter tests indicate the ratio between regional-scale effective K h and K v is 238, which tends to produce natural values of VHGR > 1 in such systems (e.g., Fogg, 1986).Groundwater development by pumping can further increase VHGR, and the augmentation of recharge due to irrigation can still further increase it.For example, values of i v in portions of the Central Valley can vary from order 10 −3 in winter (nonpumping season) to 10 −1 in summer (pumping season).Hence, taken independently, the four steady-state VHGR flow fields represent stages along a transition from ambient, primarily horizontal groundwater flow to predominantly vertical flow induced by natural anisotropy in K and pumping and recharge, and permit an investigation of long-term plume evolution across these different hydraulic gradient conditions.
Across each VHGR, the Peclet number, Pe (i.e., the ratio of the characteristic time scales of advection to diffusion) was calculated for each hydrofacies (Table 1) where i is the vertical hydraulic gradient of the VHGR scenario and D* = 6.9E−10 m 2 /s (Weissmann et al., 2002b(Weissmann et al., , 2004)), and θ is the hydrofacies porosity.Due to the influence of vertical gradients, Pe of gravel, sand, and muddy sand were calculated using mean horizontal length scales (l x , l y ), whereas the Pe of mud and paleosol were calculated using vertical length scales (l z ).For gravel, sand, and muddy sand, θ = 0.3, and for mud and paleosol, θ = 0.5.Mud was used as a background category, thus its lateral mean lengths (l x and l y ) were already estimated in the process of Markov chain modeling of the transition probability (Weissmann et al., 1999).

10.1029/2020WR028655
6 of 20 Generally, system-wide Pe > 1 is considered advection dominant, Pe < 1 indicates diffusion dominoce, and slow advection (which may resemble a more diffusion-like process) may occur at relatively small Pe < 1 (Guswa & Freyberg, 2000).In heterogeneous aquifer-aquitard systems, which vary in terms of the hydraulic conductivity field, porosity, and hydraulic gradient, Pe varies spatially.Importantly, the Pe of low-K mud and paleosol hydrofacies in this study can transition between advection and diffusion dominoce due to changes in flow rates in these hydrofacies caused by different VHGR, with consequences for macroscopic transport behavior.

Metrics of Non-Fickian Transport
In order to compare the impact of various VHGR scenarios on anomalous transport, we measure tailing by fitting a tempered fractional advection-dispersion equation (T-fADE) (Meerschaert et al., 2008) to breakthrough curves and comparing the tempering parameter.Tailing is the late-time arrival of solutes at a control plane, and indicates mass holdback in fines.We also compute the time and space evolution of the mean squared displacement (MSD), or second spatial moment, which measures the variance in plume displacement and represents how "spread out" all of the particles are from one another due to mass holdback and preferential flow.Lastly, we investigate the spatial distribution of final particle positions with respect to the site's geologic heterogeneity, which illustrates preferential flow.
Between VHGR scenarios, travel times are rescaled by  / c c t l v, where t c is the characteristic travel time to travel the average facies length, and l c and v are the characteristic length and mean velocity along the mean flow direction y′ (details provided in Section S1 and Table S1).

Tailing
In order to quantify differences in tailing, a tempered fractional advection-dispersion equation (T-fADE) model (Meerschaert et al., 2008) was fit to breakthrough at a control plane along the bottom face of the model (5) where is the fractional capacity coefficient, which governs the immobile fraction; γ [−] is the time fractional order which controls non-Fickian plume variation; and λ [T −1 ] is the tempering parameter.Importantly, λ represents the characteristic time scale over which the waiting times in the immobile phase transfer from a power-law to exponential distribution, and is physically related to preasymptotic scaling due to mass retention times in large, low-conductivity immobile zones (Meerschaert et al., 2008;Zhang et al., 2007).Thus λ provides a quantitative measure of tailing due to mass holdback in fines: as λ approaches infinity, the T-fADE reverts to the classic advection-dispersion equation, and a decreasing λ indicates longer tailing.Because variations in fitted β and γ may also produce non-Fickian tailing, these parameters are kept relatively constant between VHGR scenarios (Table 2).
Tailing is quantified in the vertical dimension, because downward contaminant migration is of practical concern to water quality in groundwater wells, and under certain conditions may exacerbate site cleanup efforts.A T-fADE model was fit to breakthrough measured at a control plane along the bottom face of the model for each scenario, then compared to fitted values of the tempering parameter λ between scenarios (e.g., decreasing λ indicates more tailing).
We are principally interested in tailing along the bottom face, yet, as VHGR decreases, horizontal flow causes more particles to exit the down-gradient face rather than the bottom face.However, a substantial fraction of 10,000 initial particles exit the bottom face in each simulation: 78.break through the bottom of the domain after 1,000 characteristic times in VHGR = 100, 10, 1, and 0.1, respectively.Thus, to ensure the comparability of breakthrough, we assumed that the sample of particles that exit the bottom face of the domain represented macroscopic breakthrough along this plane if the domain was infinitely wide in the x and y dimensions, and if simulation time was unlimited.

Spatial Variance
Next, the time evolution of the plume spatial variance is calculated.The mean squared displacement along a given direction x is where  denotes the average position in the x direction over all particles, and every particle is first initialized to an identical reference location (i.e., centered spatial moment).We calculate the spatial variance in the longitudinal, transverse vertical, and transverse horizontal directions using 6.Fickian diffusion implies that plume MSD varies as t ϕ , where ϕ = 1 (Berkowitz & Braester, 1991;Neuman, 1990).By contrast, non-Fickian transport occurs when ϕ does not equal unity.Superdiffusion and subdiffusion refer to ϕ > 1 and ϕ < 1, respectively.In order to ensure the MSD represented an accurate measure of macroscopic spatial variance of the plume (and not boundary effects created by early exiting particles), we restrict our analysis of the MSD results to before 32 characteristic time scales (t/t c ≈ 10 1.5 ).

Final Particle Position
Lastly, we compare spatial distribution of final particle positions (i.e., the location at which particles exit the domain), and the geologic heterogeneity of the domain.The spatial patterns of particle exit points strongly illustrate preferential flow along high-K hydrofacies and lend insight into transport dynamics across VHGR scenarios.

Tailing in Breakthrough
The fitted T-fADE models agree well with the simulated breakthrough across all VHGR scenarios (Figures 2a-2d).Importantly, as VHGR decreases, the fitted tempering parameter λ also decreases even as β and γ remain relatively constant (Table 2), suggesting that between VHGR = 100 and VHGR = 0.1, the degree of non-Fickian tailing as measured by λ decreases by about a factor of 10 2 , roughly corresponding to the 10 3 change in vertical gradient across the scenarios (Table 1).
Late-time slopes of the breakthrough curves were also measured, and are consistent with the T-fADE results: as VHGR decreases and flow is increasingly horizontal, more tailing is observing, and consequently, late-time slopes are increasingly positive (Figure 2e).Furthermore, the change in λ between VHGR = 100 and 10 and between VHGR = 10 and 1 (about an order of magnitude difference) decreases to about a twofold difference of λ change between VHGR = 1 and 0.1 (Table 2).These differences in the change in λ are consistent with the change in the slope of the breakthrough curves (Figure 2e).Hence, results suggest that as VHGR decreases and the transport is increasingly dominated by diffusion limited trajectories, the latetime tails reflect this limit.

Temporal Evolution of the Plume Spatial Variance
Temporal evolution of the centered second spatial moments (Equation 6) in the longitudinal and transverse directions are superdiffusive and non-Fickian at early times across all VHGR scenarios; late-time transitions either remain persistently preasymptotic, or revert to Fickian scaling, depending on the hydraulic gradient conditions (Figure 3).Particles appear well mixed in all hydrofacies (Figure S2) and most particles still remain in the domain after 32 characteristic times, suggesting that plume ergodicity is reached by this time, and thus the interpretation of our results assumes the observed scaling regimes would persist if the domain was infinitely large and simulation time unlimited.We now summarize spreading in the three directions: PAULOO ET AL. 10.1029/2020WR028655

Longitudinal Spatial Variance
The magnitude of longitudinal spatial variance (Figure 3a) increases with decreasing VHGR, and is preasymptotic and non-Fickian after 32 characteristic times for VHGR = 10, 1, and 0.1.In contrast, when VHGR = 100, late-time longitudinal spreading is asymptotic and nearly Fickian (ϕ = 0.95).Spatial plume variance is around 4 orders of magnitude greater for VHGR = 0.1 compared to VHGR = 100 at all times.This large difference in the scale of spreading results in a wider-ranging residence time distribution between the scenarios (Figure 5a).A strong vertical hydraulic gradient when VHGR = 100 forces particles through nearly all facies in the downward, longitudinal direction (perpendicular to bedding layers), effectively generating an ergodic sample of the system heterogeneity much more rapidly than in scenarios with lower VHGR, thus late-time spreading becomes in more Fickian nature.By contrast, at low VHGR, some of the plume remains trapped in fines, and some spreads longitudinally along preferential pathways-both of which stretch the plume.
The smaller magnitude of transverse vertical spatial variance observed when VHGR = 100 and 0.1 is explained by plume interactions with aquifer heterogeneity.When flow is vertical and downward (VHGR = 100), transverse transport is horizontal, tends to occur along extensive preferential pathways, and is driven toward these pathways by a powerful vertical gradient and all advection dominant facies.When flow is horizontal (VHGR = 0.1), transverse vertical transport is limited by diffusion dominant, laterally extensive paleosols and muds, which prevent spreading along this direction.The middle scenarios (VHGR = 1 and 10) by contrast have stronger vertical gradients to move particles past the first laterally extensive low-K layers, and then a combination of advection, slow advection, and diffusion dominant facies, which creates in greater transverse vertical spreading.
The smaller magnitude of transverse vertical spatial variance observed when VHGR = 100 and 0.1 compared to VHGR = 10 and 1 is explained by plume interactions with aquifer heterogeneity.When flow is vertical and downward (VHGR = 100), transverse vertical transport occurs along the y direction, which due to relatively high K′ y compared to K′ z should facilitate rapid spreading.However, because all facies are advection dominant (Table 1) under the influence of powerful vertical gradients, particle trajectories tend to move relatively straight down along nontortuous particle trajectories.For a different reason, when flow is horizontal (VHGR = 0.1), transverse vertical transport is limited by diffusion dominant, laterally extensive paleosols and muds, which prevent spreading along this direction.The middle scenarios (VHGR = 1 and 10) have a combination of advection, slow advection, and diffusion dominant facies, which simultaneously hold some mass back, and also move particles through laterally extensive low-K facies and layers, all of which results in relatively greater transverse vertical spreading.

Transverse Horizontal Spatial Variance
The absence of a gradient in the transverse horizontal direction (Figure 3c) limits superdiffusive early time scaling (1.13 ≤ ϕ ≤ 2.09), and leads to late-time asymptotic Fickian scaling (ϕ = 1) after 32 characteristic times across all VHGR scenarios.We suspect that transverse horizontal early time superdiffusive scaling when VHGR = 100 is related to the same processes that create early time superdiffusive transverse vertical spreading: migration into preferential pathways.However, since the mean lengths of gravels (and the PAULOO ET AL.In all directions and across all VHGR scenarios tested in this study (except for the VHGR = 100 in the longitudinal direction), superdiffusive early time scaling becomes more Fickian in late times as particles explore more of the system space, and hence the transport is increasingly ergodic.When VHGR = 100, late-time longitudinal and transverse second moments approach Fickian scaling (rather than remaining persistently asymptotic), which suggests that an advection dominant system with powerful vertical hydraulic gradients can have a homogenizing impact on the anomalous transport created by system-scale heterogeneity.
PAULOO ET AL.Furthermore, it illustrates how the degree of non-Fickian effects created when the mean flow direction changes (e.g., persistent superdiffusive scaling) may challenge upscaled approaches that currently fail to capture non-Fickian effects under similar boundary condition transience (Guo et al., 2020).

Influence of Geology on Preferential Flow and Tailing
Non-Fickian transport is influenced by geologic heterogeneity in the hydrofacies model, including a high-K, interconnected incised valley fill, laterally extensive low-K paleosols, and the large volume fraction of low-K mud that plays an important role in diffusion and slow advection dominant transport.
PAULOO ET AL.The incised valley fill runs laterally along the y direction, and interconnects to an exit on the down-gradient side of the domain (Figure 4a and 4b).As VHGR decreases, particles take increasingly horizontal trajectories, and the incised valley fill and other interconnecting sand and gravel lenses act as a conduit for preferential flow that cause more particles exit the down-gradient side of the domain (Figure 4d,4f,4h,and 4j), where they form clusters around gravel and sand facies (Figure 4b).Across all VHGR scenarios, particle trajectories that exit via the incised valley fill on the down-gradient side tend to exit at earlier time scales compared to other particles that exit the down-gradient side (lighter colored dots in the black oval compared to darker colored dots without outlines in Figure 4c, 4e, 4g, and 4i), and compared to particles that take longer flow paths and exit the bottom of the model at late times (darker colored dots in Figure 5c, 5e, 5g, and 5i).
Notably, even at VHGR = 100 when most particle trajectories take vertical paths and exit the bottom face of the domain (yellow points in Figure 5d), some particles rapidly migrate in the transverse vertical direction along the incised valley fill to exit the down-gradient side (purple dots in the black oval in Figure 4d).Moreover, early time arrival through the incised valley fill when VHGR = 100 (Figure 4c) suggests that transverse horizontal preferential flow still occurs along interconnected coarse facies even under strong vertical gradients, and explains the ballistic early time transverse vertical scaling (ϕ = 2.84) in this scenario (Figure 3b).
Laterally extensive low-K paleosols limit (Figure 4b) the vertical migration of particles, and create segmented breakthrough along the down-gradient side of the domain where particles tend to exit above or below paleosol layers, but rarely through them (Figure 4d, 4f, 4h, and 4j).

Mean Flow Direction Modulates the Degree of Non-Fickian Transport
In clastic sedimentary environments, significant regional-scale vertical anisotropy in K caused by stratification of facies and seasonal pumping and recharge can create strong vertical hydraulic gradients that in turn cause oscillations in mean flow direction.This study aimed to assess impacts of these shifting flow directions on the regional-scale, non-Fickian transport of nonpoint source contaminant in a typical alluvial aquifer-aquitard complex.Results demonstrate that shifting VHGR can modify the degree of non-Fickian transport.Decreasing VHGR led to increasingly non-Fickian transport, measured by tailing at a control plane along the bottom face of the model, increased and persistent non-Fickian scaling of the second spatial moment, and increased preferential flow along interconnected high-K facies.
Our results are consistent with previous studies which demonstrate that non-Fickian spreading measured (at a target well) can accompany the change from plume migration along an ambient groundwater gradient, to transient pump-and-treat conditions (Guo et al., 2019;Guswa & Freyberg, 2000;LaBolle & Fogg, 2001), which facilitates the release of mass stored in fine-grained material.Results illustrate that non-Fickian spreading of a nonpoint source plume depends on the mean flow direction, modified in many aquifer systems worldwide by ubiquitous pumping wells (Gleeson et al., 2012;Pedretti et al., 2016;Wada et al., 2010) that increase vertical head gradients and drive greater vertical, commonly downward, flow.
Together, the increasingly non-Fickian effects (e.g., increased tailing, greater spatial spreading) that result as the mean flow direction shifts from horizontal to vertical suggest a conceptual model of nonpoint source transport in agriculturally intensive alluvial groundwater basins characterized by oscillating patterns of more Fickian vertical migration during periods of intense pumping and recharge, interspersed by periods of persistently non-Fickian lateral plume stretching under ambient groundwater flow when pumps are off and fields are not irrigated, except perhaps in areas with less vertical anisotropy in K.These insights into the dependence of non-Fickian transport on mean flow direction bear importantly on efforts to model and manage nonpoint source contamination like salts (Hansen et al., 2018) and nitrates (Pastén-Zapata et al., 2014) in large irrigated basins where the flow boundary conditions (and hence VHGR) are inherently transient.Thus, an unresolved challenge for existing methods that upscale the advection-dispersion equation (Neuman & Tartakovsky, 2009) is the ability to capture the varying non-Fickian behavior of regional, nonpoint source transport resulting from changes in the mean flow direction (Guo et al., 2020).This study advances a hydrogeologic perspective toward understanding the physical genesis of this challenge.

Hydrogeologic Heterogeneity Underpins Non-Fickian Transport
Differing degrees of non-Fickian transport and plume spreading arise from hydrogeologic heterogeneity in the study site, which can be understood in terms of the transition from diffusion to advection dominant transport and persistently non-Fickian transport caused by preferential flow and mass holdback.

Transitions Between Diffusion and Advection Dominant Transport
The results show that transience in the hydraulic gradients can create systems that range from fully advection dominant, to systems of mixed advection and diffusion dominoce, the latter of which will have greater non-Fickian effects.
As VHGR increases and flow is increasingly vertical, the hydrofacies Peclet number (Table 1) and hence advection dominant transport also increases.Importantly, in our study site, as VHGR increases from 0.1 to 100, the mud and paleosol facies shift from diffusion dominant (Pe < 1) to slow advection and advection dominant (Pe > 1).Consequently, high VHGR causes mass to migrate nearly straight through fines, consistent with the findings of Guswa and Freyberg (2000), who observed accelerated solute release from low-permeability zones under the influence of pumping in a slow advection mass transfer system.In contrast to a diffusion dominant system, pumping would have little to no effect on the mass release from fines.Slow advection is thus more related to the groundwater velocity, as opposed to a mass transfer model that strictly assumes diffusive controls on mass release from low-permeability lenses.In this study, the decrease in non-Fickian transport observed at increasing VHGR and hence Pe leads us to speculate that hydraulic forcings may exert a similar effect on transport as aquifer-aquitard spatial properties (e.g., Yin et al., 2020;Zhang et al., 2006Zhang et al., , 2007Zhang et al., , 2013)).For instance, Zhang et al. (2013) found that in domains containing hydrofacies with relatively short mean lengths, particles moved more rapidly between hydrofacies, which increased the likelihood for solutes to sample the full range of flow velocities, and thus more rapidly reach ergodicity.
Similarly, we find that strong vertical flow (e.g., VHGR = 100) can push particles through all facies alike, thus accelerating mass transport between hydrofacies and rapidly reaching ergodicity.Furthermore, Zhang et al. (2007Zhang et al. ( , 2013) ) observed that diffusion limited by vertical floodplain layer thicknesses controlled late time non-Fickian solute tailing.Similarly, our results suggest that higher VHGR can shift facies from diffusion to advection dominoce and diminish non-Fickian tailing.Hence, non-Fickian transport imparted by spatial heterogeneity is continually modified by hydraulic transience.

Persistent Non-Fickian Transport
Describing persistent, late time, non-Fickian transport observed in this study (Section 3.2) remains a modeling challenge.We now discuss the factors that control persistently preasymptotic non-Fickian transport, specifically how it relates to hydrofacies heterogeneity, preferential flow, and the transition from diffusion to advection dominant transport.
First, the super-linear growth of early time plume spatial variance observed in this study (i.e., superdiffusion, ϕ > 1) is consistent with other simulations in heterogeneous porous media (Hunt et al., 2011;Kang et al., 2014;Zhang et al., 2013).The conventional understanding of this phenomenon is that preferential flow in highly permeable facies stretches the plume front downstream, while the trailing plume tail remains trapped by low-permeability facies near the source (Dagan et al., 2003;Fiori et al., 2003;LaBolle & Fogg, 2001;Zhang et al., 2014).As ergodicity is reached in late times, plume spatial variance approaches Fickian scaling, which has been observed in numerous experiments (e.g., Cortis et al., 2004;Kang et al., 2014;Le Borgne et al., 2010) and is relatively easier to model.However, persistently preasymptotic MSD scaling observed in some VHGR scenarios reinforce the notion suggested by other researchers (e.g., LaBolle & Fogg, 2001;Matheron & De Marsily, 1980;Neuman, 1990;Zhang & Lv, 2007) that the late-time transition from non-Fickian to Fickian transport may never occur in natural geologic formations under certain gradient conditions, because the time needed to reach asymptotic behavior may be sufficiently large that the solute will encounter new spatial nonstationarities along the displacement distance before Fickian scaling is reached, and will thus always remain in a state of transition, never reaching asymptotic behavior.
We speculate that the multiple-scale geologic heterogeneity in our study area causes persistently non-Fickian superdiffusion at late times in longitudinal and transverse directions for all VHGR ≤ 10.This persistent late time non-Fickian transport is supported by LaBolle and Fogg (2001), who noted persistent preasymptotic PAULOO ET AL.
10.1029/2020WR028655 ambient transport in a similar well-connected, hydrofacies-based alluvial fan model.Moreover, Zhang and Lv (2007) found that transition times in simple, 2D uniform media were so large, that the authors speculated that in natural geologic media with multiple-scale heterogeneity, the transition from anomalous to Fickian transport may never complete.Complimentary to these observations, our results indicate that sufficiently high VHGR can reduce the impact of mass sequestration in fines, lead to a more compact plume in all directions, and hence, promote asymptotic late time Fickian scaling in all dimensions, thus lending itself toward being more appropriately described by an upscaled model.
The important role of preferential flow on persistent non-Fickian transport is illustrated by the impact of a lateral, connected, incised valley fill and other prominent interconnecting sand and gravel bodies (Figure 4a).When VHGR = 100, rapid particle migration into the incised valley fill transverse vertical to the mean flow direction resulted in early time ballistic scaling (ϕ = 2.84) of transverse vertical MSD that slowed to nearly Fickian scaling (ϕ = 0.95) at late times as ergodicity was reached.In contrast, when VHGR = 10, 1, and 0.1, preferential flow along high-K networks stretched some of the plume at the same time that fines sequestered other parts of it, leading to persistent and oftentimes a greater magnitude of late-time spreading (ϕ = 1.34-1.36).Across all scenarios, particles traveled laterally through the incised valley and exited the down-gradient side at relatively earlier times (Figure 4c).Observations of faster transport along highly permeable, interconnected sand and gravel lenses are consistent with numerous studies of preferential flow (Fogg, 1986;Heeren et al., 2010;Jussel et al., 1994;LaBolle & Fogg, 2001;Maxwell et al., 2008;Moreno & Tsang, 1994;Tsang & Neretnieks, 1998;Weissmann et al., 2004;Winograd & Pearson, 1976;Zheng & Gorelick, 2003).
Final particle positions and the hydrofacies they occupy at that time further illustrate the importance of advection and diffusion dominant transport and preferential flow.When VHGR = 100 and advection dominates transport, particles are nearly indiscriminate in terms of the facies they occupy upon exiting the bottom face of the domain (i.e., even spatial distribution of yellow dots in Figure 5d).By contrast, when VHGR = 0.1, spatial patterns of particles that exit the bottom face of the domain (yellow dots in Figure 5j) and the down-gradient side (purple dots in Figure 4j) correlate with more conductive, high-K facies (blue and dark blue sand and gravel in Figures 5b and 4b).This indicates preferential flow along interconnecting sand and gravel lenses, and partly explains the increased non-Fickian transport with decreasing VHGR as a result of relatively rapid plume stretching along these fast flow paths.

Implications for Upscaled Transport Models
The varying degrees of non-Fickian transport that arise from the interaction of hydraulics and geologic heterogeneity challenge the application of the advection-dispersion equation toward subsurface transport problems in natural geologic media.At regional scales (e.g., tens to hundreds of kilometers), subsurface characterization of the relevant heterogeneity may be intractable or cost-prohibitive, and even if reliable subsurface models were available, the computational complexity of solving transport over these scales would be enormous.Hence, regional upscaled transport models are needed.These upscaled models must be able to represent non-Fickian effects created by preferential flow in connected high-conductivity facies, mass transfer processes between high-conductivity and low-conductivity media, and as we show in this study, the different degrees of non-Fickian transport introduced by changes in the mean flow direction due to varying VHGR (e.g., in the case of significant pumping and recharge).
In a previous study, Guo et al. (2020) demonstrated an adaptive multirate mass transfer (aMMT) model which improved breakthrough curve estimation for diffusion dominant transient flow by incorporating time-dependent mass transfer coefficients based on the equivalent flow velocity.However, the aMMT approach failed to capture tailing in breakthrough curves when the mean flow direction changed, as in the case of shifting VHGR.In this study, we extend the work of Guo et al. (2020) and illustrate that as VHGR changes, the transition from diffusion to slow advection and advection dominoce in the discrete hydrofacies of a system can fundamentally change the mass transfer processes.Specifically, high VHGR (provided it creates system-wide advection dominoce) can cause solutes to move nearly straight through facies (e.g., Figures 1 and 5d) even though the vertical effective K is much less than the horizontal effective K, whereas low VHGR (diffusion to slow advection dominoce some facies) can lead to mass holdback in fines (Figure 2) and strong preferential flow (Figures 4h, 4j and 5h, 5j).Thus, a future promising direction toward the goal PAULOO ET AL.
10.1029/2020WR028655 of improving estimation of late time tailing under transient flow is to quantify changes in VHGR, Pe, mean flow direction, and characteristic length scales along the mean flow direction, and relate these to mass transfer processes parameters in a upscaled model (e.g., aMMT).

Study Limitations and Additional Considerations
We explore highly detailed nonpoint source plume evolution across a series of steady-state flow simulations in order to observe regional-scale, long-term transport under conditions that represent the transition from relatively ambient groundwater flow to predominantly vertical flow.However, these steady-state simulations do not the capture mass transfer effects that occur when hydraulic boundary conditions rapidly change, as in the case of a pumping well being turned on.For example, Bastani and Harter (2020) found that transport models based on steady-state flow did not capture short-term oscillations of nitrate concentrations in pumping wells at the local scale.Similarly Guo et al. (2020) showed that the aMMT model failed to reproduce transport when the flow direction changed, and (LaBolle & Fogg, 2001) observed the release of mass stored in fine-grained material when boundary conditions changed from plume migration along an ambient groundwater gradient, to transient pump-and-treat conditions.Thus, there is a need for future studies to explore nonpoint source plume evolution in a transient simulation that incorporates transitions between each VHGR scenario, and which represents the hydraulic conditions observed in natural irrigated basins with extensive pumping.
Heterogeneous porous media affects reactive solute transport by controlling solute mixing and residence times (Perez-Fodich & Derry, 2019), which may be more important for certain solutes and depth zones.In this study, we simulate a nonreactive solute, which may not be suitable for some nonpoint source contaminant in fully saturated groundwater.More research needed on how reaction rates may alter the non-Fickian effects explored in this study, for instance, by mineral dissolution of the domain which adds a source term (Schoups et al., 2005), or in the case of nitrates, by denitrification in the soil zone (Lee et al., 2006).
Exiting particles and a finite-size domain challenge the ability to measure long temporal scaling regimes.We assumed plume ergodicity was reached after 32 characteristic times based on observed mixing in all hydrofacies (Figure S2), but the evolution of plume spatial variance under varying VHGR in a domain large enough to track all particles across orders of magnitude of characteristic times would show how persistent the scaling regimes observed in this study are over longer time scales.
It is well known that hydrofacies geometry influences contaminant transport, and hence, our study is best interpreted with respect to the characteristic geology of the domain.For example, Yin et al. (2020) found that larger hydrofacies width increased high-permeability connectivity, and thus superdiffusion; conversely, larger mean hydrofacies thickness slowed solute migration and reduced plume scaling.Any extrapolation the results presented herein to other areas should be accompanied by a careful consideration of the hydrogeology of the new study site.

Conclusions
We simulated 3D steady-state groundwater flow and nonpoint source contaminant transport in a regional scale (15 km × 12.6 km × 100 m) geostatistical representation of a heterogeneous alluvial aquifer-aquitard system.In horizontally stratified clastic sedimentary aquifer systems, it is common for the strong vertical anisotropy of hydraulic conductivity to produce spatially and temporally varying vertical and horizontal components of the hydraulic gradient, in turn producing a complex interplay between the stratified heterogeneity and physical transport processes.These effects can be modified and made highly transient by the effects of pumping and recharge in, e.g., an agricultural system where oscillations between primarily horizontal and primarily vertical flow are caused by seasonal fluctuations in pumping as well as the augmented recharge commonly produced by irrigation.Thus, we tested the impact of these varying mean flow directions on non-Fickian transport, and observed increasingly non-Fickian transport as flow was increasingly horizontal (i.e., when the VHGR decreased).Importantly, when VHGR was large (e.g., 100), the plume was more compact: particle trajectories exhibited late-time asymptotic Fickian scaling of spatial variance, and moved perpendicular to the direction of aquifer bedding (even though the vertical effective K is much less than the horizontal effective K), and nearly straight through high-conductivity and low-conductivity facies PAULOO ET AL.
10.1029/2020WR028655 alike.As VHGR decreased (e.g., 10 and 1) flow became more horizontal and transport was more non-Fickian-particle trajectories increasingly moved along laterally extensive facies, took preferential flow paths, and exhibited greater (and preasymptotic) spatial spreading and tailing in breakthrough measured along a control plane at the bottom face of the model.Hence, the non-Fickian transport imparted by spatial heterogeneity is continually modified by hydraulic transience that may transition facies from diffusion to advection dominoce and thus alter non-Fickian mass transfer dynamics.These results illustrate why conventional methods to upscale anomalous transport under transient flow conditions remains a challenge, and indicates that future efforts to upscale regional-scale, transient nonpoint source transport may benefit by incorporating information on the time-dependent changes in VHGR, Pe, mean flow direction, and characteristic length scales along the direction of mean flow.

Figure 1 .
Figure 1.The Kings River alluvial fan hydrofacies model in California's Central Valley (USA) with a section cut away to reveal the geologic heterogeneity and four representative particle trajectories from the four VHGR scenarios in this study (Section 2.3).The location of the instantaneous nonpoint source term (10,000 particles, 10 m below the water table) is represented by a horizontal pink line.Particle trajectories are increasingly horizontal and non-Fickian with decreasing VHGR.Note that the vectors representing flow in z, y (and their geometric sum in the mean flow direction, y′) are illustrative, and not drawn to scale.
and α T are the longitudinal and transverse dispersivities [L], I is the identity matrix, and |q| is the magnitude of the velocity vector [LT −1 ].
for instantaneous pulse injection can provide an insight into a continuous injection case, which can be solved by superposition.LaBolle andFogg (2001) demonstrated for a similar hydrofacies-based domain that solute transport results were insensitive to local scale α, which was insignificant compared to the dispersion caused by the hydrofacies-scale heterogeneity of the conductivity domain.Moreover,Weissmann et al. (2002a)  investigated sensitivity analyses with α L and α T = 0.1 m and = 0.01 m in the same Kings River fan domain that we use in this study, and found that between geostatistical realizations, the simulation results did not differ appreciably in terms of the overall transport.Thus, we set longitudinal, transverse horizontal, and transverse vertical dispersivities to 0.01 times the cell discretization, α L = α TH = 2.0 m, and α TV = 0.005 m.In the transport simulations, most of the dispersion occurs due to the spatially variable velocity field and hence the dispersivity parameter is of little importance here.These local-scale dispersivities were maintained across the VHGR scenarios in order to ensure comparability.We assume a conservative (nonreactive, nonsorbing) solute and set the molecular diffusion coefficient D* = 6.9E−10 m 2 /s, based onWeissmann et al. (2002a).

Figure 2 .
Figure 2. (a-d) RW3D breakthrough curves (colored points) and fitted T-fADE models (black line) measured at the bottom face of the domain (Figure 1) for the four VHGR scenarios.Note the difference in x axis travel times.(e) Breakthrough curves are normalized by the characteristic travel time (see TableS1for details).Measured slopes of tails (denoted as m) increase as VHGR decreases, corresponding to increasingly anomalous transport.T-fADE, tempered fractional advection-dispersion equation; VHGR, vertical to horizontal hydraulic gradient ratio.

Figure 3 .
Figure 3. Temporal evolution of the centered second spatial moments (Equation6) calculated from the distribution of particle trajectories along the (a) longitudinal, (b) transverse vertical, and (c) transverse horizontal directions indicate superdiffusive early time scaling that reduces to less superdiffusive or Fickian scaling in late times, depending on the VHGR and direction of scaling.Triangles at early and late times show the slope of MSD temporal scaling, ϕ (Section 2.4.2).VHGR, vertical to horizontal hydraulic gradient ratio; MSD, mean squared displacement.

Figure 4 .
Figure 4. (a) Transparent 3D hydrofacies domain highlighting a lateral, interconnected sand and gravel incised valley fill; (b) the down-gradient side of the domain.A black oval marks the incised valley fill exit, and black outlines mark prominent interconnecting sands and gravel bodies that promote preferential flow.These same regions are marked in black in all other subplots.(c, e, g, i) The final z and x position of particles, colored by the time they exit the domain (residence time); (d, f, h, j) the final z and x position of particles, colored by the length (y) into the page, where yellow is the up-gradient side and purple is the down-gradient side.Particle clusters at z = 0 m indicate trajectories that exit the bottom face (Figure 5).

Figure 5 .
Figure 5. (a) The 3D hydrofacies domain; (b) the bottom of the domain with prominent sand and gravel bodies outlined in black; (c, e, g, i) the final x and y position of particles, colored by the time they exit the domain (residence time); (d, f, h, j) the final x and y position of particles, colored by the length (z) into the page, where yellow is the bottom face and purple is the top face.Particle clusters at y = 15,000 m indicate trajectories that exit the down-gradient side of the domain (Figure 4).

Table 1 Hyrdofacies
10.1029/2020WR028655 5 of 20 Hydrofacies l x (m) l y (m) l z (m) K (m/s) Properties Including Mean Lengths l in xyz Directions (Strike, Dip, and Vertical in 1), Hydraulic Conductivity K (m/s), and Average Peclet Number Across the VHGR Scenarios

Table 2
Fitted T-fADE Model Parameters Indicate an Increasingly Small Tempering Parameter as VHGR Decreases, Consistent With the Observation of Increased Tailing in Low VHGR Scenarios