Application of Tempered-Stable Time Fractional-Derivative Model to Upscale Subdiffusion for Pollutant Transport in Field-Scale Discrete Fracture Networks

Fractional calculus provides efficient physical models to quantify non-Fickian dynamics broadly observed within the Earth system. The potential advantages of using fractional partial differential equations (fPDEs) for real-world problems are often limited by the current lack of understanding of how earth system properties influence observed non-Fickian dynamics. This study explores non-Fickian dynamics for pollutant transport in field-scale discrete fracture networks (DFNs), by investigating how fracture and rock matrix properties influence the leading and tailing edges of pollutant breakthrough curves (BTCs). Fractured reservoirs exhibit erratic internal structures and multi-scale heterogeneity, resulting in complex non-Fickian dynamics. A Monte Carlo approach is used to simulate pollutant transport through DFNs with a systematic variation of system properties, and the resultant non-Fickian transport is upscaled using a tempered-stable fractional in time advection–dispersion equation. Numerical results serve as a basis for determining both qualitative and quantitative relationships between BTC characteristics and model parameters, in addition to the impacts of fracture density, orientation, and rock matrix permeability on non-Fickian dynamics. The observed impacts of medium heterogeneity on tracer transport at late times tend to enhance the applicability of fPDEs that may be parameterized using measurable fracture–matrix characteristics.


Introduction
Fractional calculus, defined by non-integer order derivatives and integrals, has been applied to problems involving non-Fickian or anomalous dynamics for almost three decades [1,2].Despite their vast potential, both theoretical development and real-world applications of fractional partial differential equations (fPDEs) have been commonly constrained by the lack of understanding of how earth system properties influence non-Fickian transport dynamics, especially for the hydrologic sciences [3].This major challenge has historically reduced fPDEs to curve-fitting mathematical exercises, instead of routine hydrological modeling tools [4].Intensive efforts are needed to link medium properties and the fPDE model parameters, and the commonly used Monte Carlo approach is particularly suitable for providing target dynamics for heterogeneous systems with well-controlled and definable properties.
Applications of fPDEs in the hydrologic sciences have a relatively shorter history, compared to the other stochastic hydrologic approaches, such as the single or multiple rate mobile-immobile models used extensively in chemical engineering and hydrogeology [5,6].Anomalous transport has been observed in a variety of disordered systems in natural geological media, especially in heterogeneous porous media and fractured formations [6][7][8][9][10].After Benson et al. [11] first introduced the spatial fractional advection-dispersion equation (fADE) to capture super-diffusive transport in sand tanks and a relatively homogenous aquifer, the fADEs had been applied to model anomalous transport in saturated, heterogeneous porous media [12][13][14][15] and Earth surfaces, such as natural rivers [16][17][18].Applications of fADEs for flow and transport in fractured media are rather rare [19,20], although fractures are ubiquitous in geologic systems.More than 90% of natural aquifers are fractured.
Fractured media exhibit erratic heterogeneity and scale-dependent dynamics for flow and transport [21,22], challenging standard modeling tools and providing an ideal test case for fPDEs.Accurate and efficient simulation of contaminant transport within a fractured rock mass is practically important, since fractures play a large role in many natural and engineered processes, such as long-term disposal of high-level radioactive wastes in a geologic repository [21,[23][24][25].Fracture properties of natural rocks, such as fracture density, spatial location, length, aperture and orientation distributions, account for medium heterogeneity and result in non-Fickian contaminant transport [26][27][28][29], where the feasibility of fADEs and the potential correlation between fracture properties and fADE parameters have not been fully understood (see pioneer work in [19,20]).
This study selected the fPDE to quantify pollutant transport through field-scale fractured rock.Due to the prohibitive computational burden in mapping each individual fracture in standard grid-based models [20], different transport models have been developed to describe and predict complex transport behavior in fractured media at different scales [30][31][32][33].For example, the fractured medium was treated as an equivalent continuum and the advection-dispersion equation (ADE) with ensemble averaged parameters was used to model contaminant transport [22,34].However, the equivalent continuum assumption may only be valid for rocks with sufficiently high fracture densities that approximate an equivalent porous medium [35].This condition is rare in natural fractured rock, while sparsely-fractured rock is more common, and favored for geological disposal [35,36].Unlike the advection-dispersion equation (ADE), which characterizes Fickian diffusion at local scales, non-Fickian transport behaviors and asymmetric plumes, due to heterogeneity in fractures, have been observed in both the laboratory experiments and field tests [7,31,37].This provides our motivation to select a non-local fADE in this study, which can account for both spatial and temporal nonlocal processes [38][39][40][41] in describing non-Fickian transport in fractured media.
The rest of the paper is organized as follows.In Section 2, we present the Monte Carlo approach used to simulate multiple discrete fracture network (DFN) flow and transport scenarios.The DFN is an efficient and conceptually robust numerical approach to simulate the dynamics of flow and transport in fractured media with definable fracture properties.It assumes that the rock mass is dissected by a network of discrete fractures, and that fluid flow and contaminant transport in a low-permeability rock mass are controlled by geometric and hydraulic properties of interconnected fractures of a network [42,43].Results of the Monte Carlo simulations and applications of the fADE are then shown in Section 3. In Section 4, we explore the emergence of anomalous transport and its characterization, to gain insight into the effective solute dynamics and major mechanisms behind the observed anomalous behavior.We focus on the impact of fracture density and orientation and rock matrix permeability on non-Fickian dynamics and the fADE parameters.Section 5 draws the main conclusions.

Monte Carlo Simulations and the Fractional Advection-Dispersion Equation
The Monte Carlo approach utilized in this study contains three major steps.First, we generated multiple scenarios of stochastic fracture networks, with varying properties, using HydroGeoSphere (HGS) software (v.111,Aquanty Inc., Waterloo, ON, Canada) [44], which is a multi-dimensional, control-volume, finite element simulator designed to quantify the hydrologic cycle, including simulating flow and transport in fracture networks embedded in a porous medium with the discrete facture (DF) approach.Second, groundwater flow through the generated DFNs was modeled, which was expected to lead to strongly non-Fickian dynamics for conservative tracer transport (see further discussion below), providing the synthetic data to evaluate the influence of the fractured media properties on non-Fickian dynamics.Third, pollutant transport, characterized by breakthrough curves (BTCs), was then accounted for by a truncated power-law distribution memory function embodied in the fADE model.FracFit [41], a robust parameter estimation tool using a weighted nonlinear least squares (WNLS) algorithm, was used to parameterize our fractional calculus model, which could capture salient features of anomalous transport, such as skewness and power-law tails.Weights assigned by WNLS [45] were proportional to the reciprocal of the measured concentrations.Therefore, areas of lower concentration (representing low probability regions) received greater weight, which was essential for capturing early arrivals and late time-tail characteristics.

Generating the Random Discrete Fracture Networks
The model geometry, shown in Figure 1, is 50 × 25 m 2 , and represents a horizontal, two-dimensional (2D) map view of the fractures, with unit thickness and uniform blocks (100 × 50) along the x and z axes.One hundred 2D DFN realizations (for each scenario) were generated by superimposing two different sets of fractures onto the grid, with links between the fracture elements and matrix block elements, leading to realistic DFNs [46,47].Each scenario represented the ensemble average of 100 total individual DFN flow and transport realizations.Fracture locations, orientations, lengths, and hydraulic conductivities were generated from predefined distributions obtained in literature [32,48]: (1) Fracture location: distributed randomly over the entire domain, with a random seed selected for each realization; (2) Fracture orientation: two fracture sets followed mixed Gaussian distributions (each Gaussian distribution is equally weighted and has its own mean and variance), with a mean and standard deviation of 0 • ± 10 • for the first set, and 90 • ± 10 • for the second set; (3) Fracture length: log-normal distribution, with 10 length bins-the mean of the smallest length bin is 2.5 m (which is 1/10 of the minimum length of the x and z domain sizes), and the mean of the largest length bin is 25 m (which is the minimum length of the x and z domain sizes); (4) Fracture aperture: exponential distribution, with 10 aperture bins-the mean of the smallest aperture bin is 0.001 m, and the mean of the largest aperture bin is 0.0015 m.
The matrix permeability was assumed to be isotropic and in the range of an un-fractured, low-permeability rock, like a real silt medium ([1.15 × 10−16, 2.3 × 10−12] m 2 ).One realization of the DFN, with 100 fractures, is shown in Figure 1.A horizontal fracture with the length of 10 m and aperture of 0.001 m was inserted into the middle of the left boundary, in order to inject a conservative tracer.This zone provided a mobile source of pollutant, representing real-world contamination, where a pollutant usually enters firstly the mobile region in aquifers before spreading downstream.

Modeling Groundwater Flow and Pollutant Transport through the Discrete Fracture Networks
Steady-state groundwater flow through the confined aquifer generated above was then solved by HGS.Parameters for the flow model, including the specific storage coefficient and porosity, are defined in Table 1.Boundary conditions were chosen so that the main flow direction was from left to right, with a general hydraulic gradient of approximately 0.2.The high hydraulic gradient was selected for faster convergence, and decreased simulation times in this work (so that we could observe the late-time behavior of solute transport, which is critical to identify the dynamics of subdiffusive non-Fickian transport).The left and right boundaries were set to be constant head boundary conditions, and had values as 50 m and 40 m, respectively.

Modeling Groundwater Flow and Pollutant Transport through the Discrete Fracture Networks
Steady-state groundwater flow through the confined aquifer generated above was then solved by HGS.Parameters for the flow model, including the specific storage coefficient and porosity, are defined in Table 1.Boundary conditions were chosen so that the main flow direction was from left to right, with a general hydraulic gradient of approximately 0.2.The high hydraulic gradient was selected for faster convergence, and decreased simulation times in this work (so that we could observe the late-time behavior of solute transport, which is critical to identify the dynamics of sub-diffusive non-Fickian transport).The left and right boundaries were set to be constant head boundary conditions, and had values as 50 m and 40 m, respectively.
Table 1.Model parameters and initial and boundary conditions used in the simulations.
T tracer 10 s a ω f stands for aperture; µ stands for fluid viscosity -b α ω stands for fluid compressibility --Pollutant transport through the saturated DFN was then solved by HGS, using the same ILU (Incomplete lower upper factorization)-preconditioned ORTHOMIN (an iterative method using orthogonalization and minimization to achieve fast convergence) solver as is used for the flow problem.Transport parameters are shown in Table 1.The longitudinal, transverse, and vertical transverse dispersivities assigned to both the un-fractured rock matrix and rock fractures were 0.01, 0.001, and 0.001 m, respectively, where the longitudinal dispersivity was one order of magnitude smaller than the grid size, and the horizontal dispersivity was one order of magnitude smaller than that along the longitudinal direction.Tortuosity was assumed to be 1, so that we could focus on the impact of other measurable parameters (such as matrix permeability and porosity) on transport behaviors.A 10 s pulse of one conservative solute, representing an example of sodium chloride, with a nominal concentration of 1 kg/m 3 , was injected into the fracture, explicitly inserted at the middle of the left boundary.The free solution diffusion coefficient was 1 × 10 −9 m 2 /s (which is on the same order of the diffusion coefficient for hydrogen in water).The matrix porosity (0.30) was selected after referring to literature [52,53] where a value of 0.25-0.27was assumed for the matrix porosity.We also tested a much smaller matrix porosity (0.05), and results (not shown here) revealed a tracer BTC (shifted toward younger times, likely due to a higher diffusion rate) with the overall trend similar to that with a porosity of 0.30.
The HGS transport time weighting factor was set to a fully implicit numerical scheme.Fully implicit time weighting is less prone to exhibit oscillations than central and explicit time weighting.
The governing equation for subsurface flow in porous media used in the HGS was a modified form of Richards' equation [54]: where ω m [dimensionless] is the volumetric fraction of the total porosity occupied by the porous medium, Γ ex [L 3 L −3 T −1 ] represents the volumetric fluid exchange rate between the subsurface domain and discrete fractures, Q [L 3 L −3 T −1 ] represents a source/sink term to the porous medium system, and θ s [dimensionless] is the saturated water content which is assumed equal to the porosity.The fluid flux q [LT −1 ] in (1) is given by Darcy's law: where K [LT −1 ] is hydraulic conductivity tensor, k r [dimensionless] represents the relative permeability of the medium, ψ [L] is the pressure head, and z [L] is the elevation head.A fracture is idealized as the space between two-dimensional parallel surfaces, with the tacit assumption that the total head is uniform across the fracture width.The equation for saturated flow in a fracture of width/aperture ω f [L] can be written by using the analogy of Richards Equation (1) for the porous matrix.The governing two-dimensional flow equation in a fracture has the form [44]: where the fluid flux q f [LT −1 ] is given by: where ∇ is the two-dimensional gradient operator defined in the fracture plane, k rf [dimensionless] is the permeability of the fracture, ψ f and z f [L] are the pressure and the elevation heads within the fracture, and S sf [L −1 ] is the specific storage coefficient for the fractures.The saturated hydraulic conductivity of a fracture K f [LT −1 ] with a uniform aperture ω f is given by: K f = ρgω f 2 /12µ.Because it was assumed that the fractures were non-deformable and fluid-filled, there was no contribution to the storage term from fracture compressibility.Thus, the specific storage coefficient for a fracture under saturated conditions was only related to the compressibility of water α w [LT 2 M −1 ], according to: S sf = ρgα w .Three-dimensional transport of solutes in a saturated porous matrix is described by the following advection-dispersion equation [44]: where C is the solute concentration represents the mass exchange rate of solutes per unit volume between the subsurface domain and discrete fractures, and λ is a first-order decay constant [L −1 ].Solute exchange with the outside of the simulation domain is represented by Q c [ML −3 T −1 ] which represents a source or a sink to the porous medium system.The equation for solute transport in a saturated fracture is: where , and R f [dimensionless] represents retardation factor.

Applying the Fractional Advection-Dispersion Equations to Quantify Transport
If the solute is initially placed in the mobile zone (i.e., a fracture along the main flow direction), the corresponding fADE takes the form [55]: where C m and C im [ML −3 ] denote the solute concentration in the mobile and immobile phases, respectively; are the effective velocity and dispersion coefficient, respectively; and Γ(•) is the Gamma function.Meerschaert et al. [56] generalized the fADE ( 7) and ( 8) by introducing an exponentially-truncated power-law function, which is an incomplete gamma function, as the memory function: where λ > 0 [T −1 ] is the truncation parameter in time.This modification leads to the tempered time fADE (tt-fADE) [56,57]: where C 0 m = C m (x, t = 0) denotes the initial source, located in the mobile phase or fractures.A few parameters in the tt-fADE (10) and (11) including the effective velocity v and the dispersion coefficient D, may be determined by field experiments (i.e., monitoring wells) and laboratory experiments (transport through fractured media), while the other parameters should be approximated given the plumes observed in the field.The poor predictability of the fADE model motivated us to explore the relationship between DFN properties and model parameters in Section 4.
At the late time t << 1/λ (12) the tail of BTC for solutes in the mobile declines as a power-law function: The tail of the total-phase BTC also declines as a straight line in a log-log plot with a rate: while at a much later time t >> 1/λ, the slope of the mobile/total-phase BTC reaches infinity (i.e., the late-time BTC tail declines exponentially).Therefore, the value of λ controls the transition of the BTC late-time tail from a power-law function to exponential function [57].
In the following sections, we try to link the time index γ to characteristics of the fractured porous media.

Results of Monte Carlo Simulations
The ensemble average of BTCs for all 100 realizations for each DFN scenario are shown in Figure 2, along with the best-fit solutions using the tt-fADE model (10) and (11).Three scenarios of DFN were considered, by setting the total number of fractures to 20, 60, and 80, and holding mean fracture set orientations constant to 0 • and 90 • .An additional fourth scenario explored the impact of fracture orientation on transport dynamics, by changing mean fracture set orientation to −45 • and 45 • for a 60 fracture DFN.
The best-fit parameters of the tt-fADE model (10) and (11) for all DFNs are listed in Table 2.The effective velocity v, dispersion coefficient D, and fractional-order capacity coefficient β increase with an increasing fracture density.It is also noteworthy that β has the units of fractional-order in time, which can be converted to a dimensionless capacity coefficient β 0 = β γ λ γ−1 (where β 0 represents the long-term ratio of immobile versus mobile mass).The other two parameters, including the time index γ and the truncation parameter λ, first decrease and then increase as the number of fractures increases.This trend implies that the DFN domain containing 60 fractures might exhibit higher degree of heterogeneity, and therefore translate from non-Fickian to Fickian transport at a much later time, compared to the DFN with 20 or 80 fractures.When the matrix hydraulic conductivity is set to be 10 −7 m/s, the BTCs are more sensitive to fracture density.For example, when the number of fractures is 20, the peak concentration appears at day 365, while the DFN with 60 (or 80) fractures has the peak BTC at day 50 (or day 100).The truncation parameter λ can be calculated by the reciprocal of the cutoff time, representing the transition time from non-Fickian to Fickian transport; however, in this study, the power-law portion of BTCs persisted and the transition time could not be identified directly.10) and (11).The black line is a power law approximation to the late time decay in concentration.(e) shows the empirical linear relationship between the slope of the power-law decay in concentration and the time index γ in the tt-fADE model ( 10) and (11).
Table 2. Best-fit parameters for the tt-fADE model ( 10) and (11) for discrete fracture networks (DFNs) with different fracture densities and orientation scenarios.In the legend, v denotes the effective velocity, D is the dispersion coefficient, γ is the time/scale index, β is the fractional-order capacity coefficient, λ is the truncation parameter, and SBTC stands for the late-time BTC slope in a log-log plot.

Impact of Fracture Density on Non-Fickian Transport
The simulations indicated that fracture density dominates non-Fickian transport dynamics, and therefore significantly affects tt-fADE model parameters.For example, a denser DFN contains more connected flow paths, leading to faster plume movement (i.e., a larger effective velocity v), greater plume spreading in space (corresponding to a larger dispersion coefficient D), and more fracture/matrix interaction (i.e., tracer particles have a higher probability of entering the surrounding matrix, which correlates to a larger value of the capacity coefficient β).In addition, the complex relationship between DFN density and other tt-fADE model parameters (time index γ and truncation parameter λ) can be explained by the hypothesis that moderately dense DFNs (i.e., fracture scenario  10) and (11).The black line is a power law approximation to the late time decay in concentration.(e) shows the empirical linear relationship between the slope of the power-law decay in concentration and the time index γ in the tt-fADE model ( 10) and (11).
Table 2. Best-fit parameters for the tt-fADE model ( 10) and (11) for discrete fracture networks (DFNs) with different fracture densities and orientation scenarios.In the legend, v denotes the effective velocity, D is the dispersion coefficient, γ is the time/scale index, β is the fractional-order capacity coefficient, λ is the truncation parameter, and S BTC stands for the late-time BTC slope in a log-log plot.

Impact of Fracture Density on Non-Fickian Transport
The simulations indicated that fracture density dominates non-Fickian transport dynamics, and therefore significantly affects tt-fADE model parameters.For example, a denser DFN contains more connected flow paths, leading to faster plume movement (i.e., a larger effective velocity v), greater plume spreading in space (corresponding to a larger dispersion coefficient D), and more fracture/matrix interaction (i.e., tracer particles have a higher probability of entering the surrounding matrix, which correlates to a larger value of the capacity coefficient β).In addition, the complex relationship between DFN density and other tt-fADE model parameters (time index γ and truncation parameter λ) can be explained by the hypothesis that moderately dense DFNs (i.e., fracture scenario with 60 fractures) maximize variations in matrix-fracture exchange rates, which leads to the heaviest late-time BTC tail (i.e., the smallest time index γ) and the persistent power-law BTC (i.e., the smallest truncation parameter λ).DFNs with a broader range of densities are needed in a future study to further investigate the above hypothesis.
The dominant impact of fracture density on pollutant transport may be due to two reasons.First, the number of fractures significantly affects the architecture of DFNs, including both the interconnectivity of fractures and the thickness of surrounding matrix.It is well-known that the interconnected fracture network forms the preferred flow paths in the formation [22,58,59].These pathways are "paths of least resistance", where most of the flow and solute is concentrated [7,26,36,60], resulting in early arrivals.Matrix block size, which is typically characterized using fracture spacing, affects trapping time and influences late-time tailing in BTCs.Second, networks containing higher fracture densities may yield a broader distribution of apertures [61], resulting in a larger variation of advective velocities in the system.Flow in the fracture is constrained to the two-dimensional domain of the fracture aperture, where average flux at any location within a fracture is proportional to the cube of the aperture [35,50,62].Small variations in the aperture (which is 1 × 10 3 ~1.5 × 10 3 µm in this study), therefore, can cause large, non-linear variations in advective velocity.

Impact of Fracture Orientation on Non-Fickian Transport
Fracture orientation has a complex impact on pollutant transport behavior.The DFN with mean fracture set orientations of −45 • and 45 • (where the overall flow direction is along 0 • ) provides two major directions for pollutant particles to move, while the DFN with mean fracture set orientations of 0 • and 90 • offers horizontally dominant motion for pollutants.Hence, the former has a relatively larger effective velocity and less spatial spreading (i.e., a smaller dispersion coefficient).Mean fracture set orientations of −45 • and 45 • also enhance the movement of solute particles into matrix blocks, resulting in a greater immobile mass (and hence a larger capacity coefficient β) than that for the cases with mean fracture set orientations of 0 • and 90 • .It is noteworthy, however, that the BTC peak occurs at a time determined by both the effective velocity v and the capacity coefficient β.The larger for v and/or the smaller for β, the earlier for the BTC peak.Although the DFNs with mean fracture orientations of −45 • and 45 • have a relatively larger v, the ensemble BTC peak appears later (at day 100) than that for the DFNs with mean fracture orientations of 0 • and 90 • (BTC peak occurs at day 50), because the DFNs with mean fracture orientations of −45 • and 45 • have a much larger β, which acts as part of the retardation coefficient (1 + β) when γ = 1.The DFNs with mean fracture orientations of −45 • and 45 • might also produce more tortuous pathways, which lead to greater transverse dispersion than the DFNs with mean fracture orientations of 0 • and 90 • [19].Finally, regardless of fracture orientation, late-time BTCs exhibit similar slopes, and therefore the similar time index γ, if the DFNs have the same density and a constant distribution of fracture length.The latter impact (caused by the fracture length distribution), however, cannot be reliably investigated in this study, due to the field-scale DFNs (with a relatively small domain size).In a future study, we will increase the model domain size (to regional scale) and systematically explore the impact of fracture length distributions on non-Fickian dynamics.
The above analysis implies that the fracture density is more important than the fracture orientation in affecting the power-law late-time tailing of the BTC, while both properties can affect the plume mean displacement (captured mainly by the velocity in the tt-fADE model (10) and ( 11)), plume spreading (captured by the dispersion coefficient in model ( 10) and ( 11)), and fracture-matrix mass exchange ratio (captured by the capacity coefficient in (10) and ( 11)).

Impact of Matrix Permeability on Late-Time BTC
Rock matrix permeability may affect the late-time behavior of chemical transport.To explore this potential impact, here we conducted additional Monte Carlo simulations using different matrix permeability: (1) increasing matrix permeability by 10 times (from 1.1 × 10 −14 m 2 to 1.1 × 10 −13 m 2 ) and ( 2) decreasing matrix permeability by 2 times (from 1.1 × 10 −14 m 2 to 5.5 × 10 −15 m 2 ).Results are shown in Figure 3.We did not select the smaller value for matrix permeability, because preliminary numerical tests showed that if matrix permeability was smaller than 1.1 × 10 −15 m 2 , the numerical iteration (in the flow model) could not converge.10) and (11).The black/blue line is a power-law approximation to the late time decay in concentration.
The BTC declines faster at the late time with an increasing matrix permeability.This is because the matrix domain with higher permeability leads to shorter trapping times for chemical particles.Contrarily, the average velocity is smaller for a rock matrix with smaller permeability, resulting in heavier late-time tailing in the BTC.

How to Approximate the Time Index γ Given the Limited Late-Time BTC
Although the transport time considered by our Monte Carlo simulation is very long (about 5000 years), we could not observe the transition in the BTC from persistent non-Fickian to Fickian transport behaviors (in the Monte Carlo simulations performed in this study, the total modeling time could not be longer than 2 × 10 6 days, since the concentration at the outlet after 2 × 10 6 days is so small that HGS generates negative BTCs).It is common that the measured BTC cannot last long enough to cover the whole power-law portion and its transition in the BTC at late times, due to the usually limited sampling period, as well as the low concentration at the late time, which might be below the detection limit.The time index γ cannot be directly estimated when the power-law portion of the BTC at the late time cannot be identified.
One possible way to solve the above challenge is to approximate γ using the BTC's recession limb.In this study, the recession limb of each measured BTC exhibits a straight line in a log-log plot (Figure 2), which indicates prolonged non-Fickian or anomalous transport.There is a linear relationship between the slope of the BTC recession limb (whose absolute value is denoted as |SBTC|) and the time index γ in the tt-fADE model ( 10) and ( 11): where the corresponding coefficient of determination is 0.9933 (i.e., strong correlation).Since |SBTC| The BTC data (i.e., symbols) are simulated by HGS.The red/green line is the solution of the tt-fADE model (10) and (11).The black/blue line is a power-law approximation to the late time decay in concentration.
The BTC declines faster at the late time with an increasing matrix permeability.This is because the matrix domain with higher permeability leads to shorter trapping times for chemical particles.Contrarily, the average velocity is smaller for a rock matrix with smaller permeability, resulting in heavier late-time tailing in the BTC.

How to Approximate the Time Index γ Given the Limited Late-Time BTC
Although the transport time considered by our Monte Carlo simulation is very long (about 5000 years), we could not observe the transition in the BTC from persistent non-Fickian to Fickian transport behaviors (in the Monte Carlo simulations performed in this study, the total modeling time could not be longer than 2 × 10 6 days, since the concentration at the outlet after 2 × 10 6 days is so small that HGS generates negative BTCs).It is common that the measured BTC cannot last long enough to cover the whole power-law portion and its transition in the BTC at late times, due to the usually limited sampling period, as well as the low concentration at the late time, which might be below the detection limit.The time index γ cannot be directly estimated when the power-law portion of the BTC at the late time cannot be identified.
One possible way to solve the above challenge is to approximate γ using the BTC's recession limb.In this study, the recession limb of each measured BTC exhibits a straight line in a log-log plot (Figure 2), which indicates prolonged non-Fickian or anomalous transport.There is a linear relationship between the slope of the BTC recession limb (whose absolute value is denoted as |S BTC |) and the time index γ in the tt-fADE model (10) and (11): where the corresponding coefficient of determination is 0.9933 (i.e., strong correlation).Since |S BTC | can be measured easily (here "easily" means that one can use the observed BTC, shown for example in Figure 2, to directly obtain the slope of the BTC, which is the negative value of |S BTC |), the above linear relationship might be practically useful.However, this approximation will be ineffective if the power-law portion of the BTC is too short to observe (which represents a fast transfer from non-Fickian to Fickian diffusion), or the data contains apparent noises (which can be common for the measurement of low concentrations at the BTC tails in the field).

Subdiffusive Transport Shown by Plume Snapshots
We further explored the spatial distribution of a chemical in the DFNs, by depicting plume snapshots at different times for a single realization, and the ensemble average of the Monte Carlo results (Figures 4 and 5).power-law portion of the BTC is too short to observe (which represents a fast transfer from non-Fickian to Fickian diffusion), or the data contains apparent noises (which can be common for the measurement of low concentrations at the BTC tails in the field).

Subdiffusive Transport Shown by Plume Snapshots
We further explored the spatial distribution of a chemical in the DFNs, by depicting plume snapshots at different times for a single realization, and the ensemble average of the Monte Carlo results (Figures 4 and 5).We chose the first time step (t = 1 day) to show the early arrivals for the three different fracture densities with matrix permeability of 1.1 × 10 −14 m 2 .All the three cases show that chemical particles move along the preferential flow paths.With increasing fractures in the DFNs, the plume spreads to a wider area in the vertical direction with a more uniform moving front.For the 13th time step (t = 2000 years), the snapshots exhibited apparent retention and sequestration (meaning that the We chose the first time step (t = 1 day) to show the early arrivals for the three different fracture densities with matrix permeability of 1.1 × 10 −14 m 2 .All the three cases show that chemical particles move along the preferential flow paths.With increasing fractures in the DFNs, the plume spreads to a wider area in the vertical direction with a more uniform moving front.For the 13th time step (t = 2000 years), the snapshots exhibited apparent retention and sequestration (meaning that the surrounding low-permeable matrix trapped chemical particles, especially near the source [39,63]) in the least-dense fractured network.Competition between channeling and sequestration leads to persistent non-Fickian dynamics.Fast-and slow-moving chemical particles are counterbalanced, and the effective velocity becomes larger for the DFNs with more fractures.

Conclusions
This study aims to explore the relationship between real-world aquifer properties and non-Fickian transport dynamics, so that the fractional partial differential equations built upon fractional calculus can be reliably applied with appropriate hydrogeologic interpretations.We use the Monte Carlo approach to generate field-scale DFNs where the fracture properties change systematically, and then to simulate groundwater flow and pollutant transport through the complex DFNs.For a point source located initially in the mobile phase or fracture, the late-time behavior for the BTCs simulated by the Monte Carlo approach is then explained by the tempered-stable time fractional advectiondispersion equation.We build the relationship between medium heterogeneity and transport dynamics through the combination of numerical experiments and stochastic analysis.The following five primary conclusions were found.
First, the DFN density is the major factor affecting non-Fickian transport at the late time, since (1) the number of fractures significantly affects the internal structure of DFNs; and (2) DFNs containing a higher fracture density may yield broader distribution of apertures, resulting in greater

Conclusions
This study aims to explore the relationship between real-world aquifer properties and non-Fickian transport dynamics, so that the fractional partial differential equations built upon fractional calculus can be reliably applied with appropriate hydrogeologic interpretations.We use the Monte Carlo approach to generate field-scale DFNs where the fracture properties change systematically, and then to simulate groundwater flow and pollutant transport through the complex DFNs.For a point source located initially in the mobile phase or fracture, the late-time behavior for the BTCs simulated by the Monte Carlo approach is then explained by the tempered-stable time fractional advection-dispersion equation.We build the relationship between medium heterogeneity and transport dynamics through the combination of numerical experiments and stochastic analysis.The following five primary conclusions were found.
First, the DFN density is the major factor affecting non-Fickian transport at the late time, since (1) the number of fractures significantly affects the internal structure of DFNs; and (2) DFNs containing a higher fracture density may yield broader distribution of apertures, resulting in greater variation of advective velocities.A moderate dense DFN (filled with 60 fractures) results in the most heterogenous domain and the heaviest late-time tail in the BTCs, implying that there might exist a potential threshold of fracture density for non-Fickian dynamics when both the domain size and the fracture length distribution remain unchanged.
Second, for the DFNs with the same density and different mean orientations, the late time tailing behavior is affected by the ratio of fractures along different directions.More longitudinal fractures aligned with the overall flow direction will lead to an increase in earlier arrivals, while more transverse fractures, with flow deviating from the general flow direction, will enhance chemical particles to interact with the rock matrix, resulting in a greater immobile mass (a larger β in the tt-fADE) and a heavier late-time tailing.The BTC peak time is determined by both the effective velocity and capacity coefficient β (representing the immobile mass ratio).The time index and the slope of the late-time BTC, however, are similar for the DFNs with the similar density and different mean orientations, since the fracture density dominates the distribution of the mass exchange rate between fracture and matrix (and therefore defines the index of the time fractional derivative in the tt-fADE model).
Third, a rock matrix with higher permeability will lead to a relatively more homogeneous domain with a larger time index and a steeper slope in the late-time BTC, since the trapping times in rock matrix becomes relatively shorter.
Fourth, the DFN properties can be quantitatively linked to the slope of the recession limb of the BTC.Particularly, the time index γ in the tt-fADE model has a simple, linear empirical expression with the power-law slope of the BTC recession limb, which can help estimate the index γ in practical applications where the chemical transport time is usually not long enough to cover the full range of late-time transport behavior.
Fifth, plume snapshots show that all the three fracture densities built in this study exhibit early arrivals and sequestration effects.The denser DFNs can produce more early arrivals, due to more preferential flow paths, while the sparse DFNs tend to show heavier sequestration near the source.Therefore, competition between the channeling effect (in relative high-permeable fractures or the mobile phase) and the trapping effect (in the surrounding rock matrix or immobile domains) results in complex prolonged non-Fickian characteristics in chemical transport.With predictable parameters, the fractional partial differential equations can be used as an efficient upscaling tool to address real-world contamination problems.

Figure 1 .
Figure 1.Map of the fractured domain, showing flow geometry and boundary conditions.Note that the discrete fractures are curved in the figure, after they are incorporated into the grid by HydroGeoSphere (HGS).

Figure 1 .
Figure 1.Map of the fractured domain, showing flow geometry and boundary conditions.Note that the discrete fractures are curved in the figure, after they are incorporated into the grid by HydroGeoSphere (HGS).

Figure 2 .
Figure 2. Breakthrough curves (BTCs) at the control plane x = 50 m for three different fracture densities (ac); (b,d) have the same fracture density with different orientations.The BTC data (i.e., symbols) are simulated by HGS.The red line is the solution of the tempered time fADE (tt-fADE) model (10) and(11).The black line is a power law approximation to the late time decay in concentration.(e) shows the empirical linear relationship between the slope of the power-law decay in concentration and the time index γ in the tt-fADE model(10) and(11).

Figure 2 .
Figure 2. Breakthrough curves (BTCs) at the control plane x = 50 m for three different fracture densities (a-c); (b,d) have the same fracture density with different orientations.The BTC data (i.e., symbols) are simulated by HGS.The red line is the solution of the tempered time fADE (tt-fADE) model (10) and(11).The black line is a power law approximation to the late time decay in concentration.(e) shows the empirical linear relationship between the slope of the power-law decay in concentration and the time index γ in the tt-fADE model(10) and(11).
that if matrix permeability was smaller than 1.1 × 10 −15 m 2 , the numerical iteration (in the flow model) could not converge.

Figure 3 .
Figure 3. BTCs at the control plane x = 50 m for three different fracture densities (a-c) with matrix permeability of 1.1 × 10 −13 m 2 (red lines and circles) and 5.5 × 10 −15 m 2 (green lines and crosses).The BTC data (i.e., symbols) are simulated by HGS.The red/green line is the solution of the tt-fADE model (10) and(11).The black/blue line is a power-law approximation to the late time decay in concentration.

Figure 3 .
Figure 3. BTCs at the control plane x = 50 m for three different fracture densities (a-c) with matrix permeability of 1.1 × 10 −13 m 2 (red lines and circles) and 5.5 × 10 −15 m 2 (green lines and crosses).The BTC data (i.e., symbols) are simulated by HGS.The red/green line is the solution of the tt-fADE model(10) and(11).The black/blue line is a power-law approximation to the late time decay in concentration.

Figure 4 .
Figure 4.One realization of plume snapshots for the first time step (t = 1 day) for the three different fracture densities: (a) the DFNs with 20 fractures; (b) 60 fractures; and (c) 80 fractures.The 13th time step (t = 2000 years) for the three different fracture densities: (d) 20 fractures; (e) 60 fractures; and (f) 80 fractures.

Figure 4 .
Figure 4.One realization of plume snapshots for the first time step (t = 1 day) for the three different fracture densities: (a) the DFNs with 20 fractures; (b) 60 fractures; and (c) 80 fractures.The 13th time step (t = 2000 years) for the three different fracture densities: (d) 20 fractures; (e) 60 fractures; and (f) 80 fractures.

Figure 5 .
Figure 5. Monte Carlo results of plume snapshots for the first time step (t = 1 day) for the three different fracture densities: the DFNs with (a) 20 fractures; (b) 60 fractures and (c) 80 fractures.The 13th time step (t = 2000 years) for the three different fracture densities: (d) 20 fractures; (e) 60 fractures and (f) 80 fractures.

Figure 5 .
Figure 5. Monte Carlo results of plume snapshots for the first time step (t = 1 day) for the three different fracture densities: the DFNs with (a) 20 fractures; (b) 60 fractures and (c) 80 fractures.The 13th time step (t = 2000 years) for the three different fracture densities: (d) 20 fractures; (e) 60 fractures and (f) 80 fractures.

Table 1 .
Model parameters and initial and boundary conditions used in the simulations.