Hierarchical sedimentary architecture governs basin-scale solute dispersion: from pre-asymptotic dynamics to uncertainty propagation

Real aquifers are structured as hierarchical sedimentary systems, where multiscale heterogeneity and geometric connectivity jointly govern groundwater flow and solute migration. Although the general influence of heterogeneity has been extensively investigated, the scale-dependent effects of hierarchical organization, particularly under basin-scale flow conditions, remain inadequately quantified. In this study, we reconstructed a series of three-dimensional heterogeneous sedimentary architectures at the basin scale and performed numerical simulations to explore solute dispersion behaviour. The results reveal that the geometry and connectivity of dominant lithofacies at larger scale control macrodispersivity, while finer scale heterogeneity has only a secondary influence on plume evolution. Furthermore, the evolution of macrodispersivity is characterized by a prolonged pre-asymptotic phase, approaching a quasi-steady state after around 5000 d, with an asymptotic stability value of 170 m. This timescale is nearly 10 times longer than that inferred from the Borden site, where macrodispersivity stabilizes at around 0.4 m after approximately 400 d. This indicates that basin-scale solute transport remains non-ergodic over extended times and distances. Uncertainty analysis further identifies a distinct buffering effect inherent to basin systems, in which the aggregation of numerous flow pathways dampens realization-to-realization variability caused by local heterogeneity. When integrated with previously reported laboratory- and sandbox-scale results from the same site, these findings establish a transferable framework linking hierarchical sedimentary architecture to multiscale dispersion and uncertainty. This framework advances theoretical understanding of non-Fickian transport and provides practical guidance for large-scale modeling and groundwater management in data-limited regions.

1 Introduction

The study of solute transport in heterogeneous porous media remains a central challenge in hydrogeology, with significant implications for groundwater science and engineering (Faroughi et al., 2023; Hansen and Berkowitz, 2020; Lee et al., 2018). Decades of field investigations, including the Borden, Cape Cod, MADE and Hanford tracer tests, and laboratory experiments have consistently demonstrated that spatial heterogeneity strongly governs plume migration, leading to non-Fickian dispersion, anomalous breakthrough behaviour, and scale-dependent transport parameters (Sudicky, 1986; Garabedian et al., 1991; Boggs et al., 1992; Berkowitz et al., 2006; Botter et al., 2008; Chen et al., 2010; Agbotui et al., 2025; Tellam and Barker, 2006).

Scheibe and Freyberg (1995) systematically introduced the concept of hierarchical organization from sedimentology into hydrogeology to characterize aquifer heterogeneity. They pointed out that aquifer heterogeneity is not a single-scale problem but rather results from the nested and superimposed geological units at different hierarchical levels (e.g., micro-, meso-, and macro-facies). Conceptually, this hierarchical-architecture view is consistent with sequence stratigraphy, which also organizes depositional heterogeneity in a nested manner through stratigraphic surfaces and stacking patterns. In practice, although stochastic theories and numerical upscaling frameworks have provided valuable tools to connect flow and mixing processes across spatial and temporal domains, many studies continue to compress heterogeneity into effective parameters (Dagan, 1984; Gelhar et al., 1992; Dentz et al., 2011). Such simplification obscures the mechanistic role of multiscale sedimentary architecture in shaping flow pathways and controlling plume dynamics, and can bias spreading and mixing-relevant predictions (Fitts, 1996; Neuman and Tartakovsky, 2009; Dentz et al., 2023; Lester et al. 2016; Yin et al., 2023).

To elucidate the mechanisms underlying scale dependent transport in the heterogeneous porous media, a series of Lagrangian-based models has been developed within the hierarchical-architecture framework. By characterizing the spatial organization of lithofacies across multiple scales, these models systematically integrate the multiscale heterogeneities into a unified representation, quantitatively linking sedimentary attributes (such as lithofacies volume proportion, mean length, and statistics of hydraulic conductivity) to transport metrics (such as dispersion and mixing). They have been successfully tested in laboratory and site-scale experiments (Soltanian et al., 2015b; Dai et al., 2020; Jia et al., 2023; Ma et al., 2025). Recent work further emphasizes that the role of lithofacies geometric attributes, particularly mean length and connectivity, control solute dispersion and mixing in alluvial systems (Ramanathan et al., 2008; Soltanian et al., 2020; Ershadnia et al., 2021; Soltanian et al., 2017). High-resolution numerical simulations within this framework provide more intuitive and detailed information on flow trajectories and interfacial reaction processes, thereby validating and supplementing the assumptions of the theoretical models (Ren et al., 2022).

To date, the validation of most Lagrangian-based models, along with their associated simulation efforts, has primarily been conducted at the laboratory or site scale. For regional to basin scale solute transport processes involving long migration distances and laterally extensive flow paths, there remains a lack of quantitative characterization explicitly constrained by hierarchical sedimentary structures. Although recent high-resolution, structure-resolving analytical approaches have shown that lithofacies geometry, connectivity, and conductivity contrasts can organize macroscopic dispersion and non-Fickian behavior at kilometer scale, (Carle et al., 2006; Pauloo et al., 2021; Guo et al., 2019), it remains unclear which hierarchical level of heterogeneity remains transport-dominant under basin-scale flow conditions, how long pre-asymptotic behavior persists, and whether the superposition of numerous large-scale pathways amplifies or buffers predictive uncertainty.

Against this background, this study uses the Qiqihar aquifer system as a basin-scale testbed and combines hierarchical sedimentary architecture reconstruction with groundwater flow and solute transport simulations to investigate the following key questions: (i) Is basin-scale dispersion behavior primarily controlled by dominant large-scale lithofacies structures, or does it still require the explicit characterization of heterogeneity at finer scales? (ii) How do the geometry of solute sources, sedimentary architecture and hydraulic contrasts between lithofacies influence the evolution of macrodispersion and its long-term pre-asymptotic behavior? (iii) Under field-representative basin-scale flow conditions, how do structural and hydraulic factors shape transport uncertainty and how uncertainties propagate and evolve? By comparing the results of this study with findings from previous laboratory- and field-scale studies conducted within the same conceptual framework (Dai et al., 2005; Ramanathan et al., 2010; Soltanian et al., 2015a; Ren et al., 2023), this study clearly distinguishes which controlling factors possess scale-universality and which factors manifest specifically at the basin scale.

The paper is organized as follows: Sect. 2 introduces the geographic background of the study area, borehole data, and sedimentary architecture analysis methods; Sect. 3 describes the construction of the multiscale heterogeneous structural model and the simulation process of solute transport; Sect. 4 shows the simulation results and conducts uncertainty analysis to explore the influence of sedimentary architecture and permeability parameters on solute dispersion; Sect. 5 summarizes the main conclusions and outlines future research directions.

2 Methods

2.1 Overview of the study area

The study area is located on the east bank of the Nen River in Qiqihar City, Heilongjiang Province, China (Fig. 1a). The Nen River defines the western boundary of the study area, and flows from the northeast to the southwest, with a total length of about 22 km and a width ranging from 400 to 900 m. The regional terrain is gently inclined, sloping from north to south at less than 1 ‰, which exerts a primary control on groundwater flow direction and results in a general alignment of groundwater flow with that of the river. The aquifer system consists of a confined aquifer and a phreatic aquifer. The confined aquifer generally contains groundwater of relatively high quality, whereas groundwater contamination is largely concentrated in the phreatic aquifer. The present study therefore investigates solute transport processes within the phreatic aquifer to clarify the mechanisms of solute migration under the influence of sedimentary architecture.

Figure 1 (a) Geographic location map of the study area; (b) Boreholes and cross-sections.

A total of 57 boreholes (see in Fig. 1b) were collected for this study. According to the borehole data, the phreatic aquifer is composed of thick Quaternary unconsolidated deposits. The lithofacies of the aquifer is dominated by sandy gravel, gravelly medium- to coarse-grained sand, and medium- to fine-grained sand. The aquifer thickness varies from 22.1 to 56.9 m, and the depth to the water table varies between 1.75 and 5.86 m.

The phreatic aquifer is primarily recharged by infiltration from the Nen River, as river stage is consistently higher than the groundwater level. Other recharge includes precipitation infiltration and by lateral inflow from the northern and northeastern boundaries. A relatively continuous aquitard, mainly composed of silty clay and clay, is located below this aquifer. However, it is locally absent near the Nen River and in the southern part of the study area, forming hydraulic windows that allow vertical connectivity between the two aquifer systems.

At the regional scale (kilometers), the phreatic aquifer exhibits a strongly layer-controlled, multiscale heterogeneity. Along the principal groundwater flow direction from upstream to downstream, the dominant lithofacies progressively transition from gravelly/medium-coarse sand to medium-fine sand. Away from the Nen River and into the floodplain, sediments become finer overall, with local development of low-permeability clay–silt lenses Vertically, a pervasive “coarse-over-fine” organization is observed: the upper section is dominated by gravel and gravelly coarse/medium sand with moderate sorting, whereas the lower section consists mainly of medium-fine to fine sand with better sorting.

2.2 Parameter determination

2.2.1 Sediments heterogeneity parameters

Eight cross-sections (Fig. 1b), oriented parallel to (defined as x-direction, sections 1–4) and perpendicular to (defined as y-direction, sections 5–8) the regional groundwater flow, were used to constrain the sedimentary structure characteristics. Following the multiscale hierarchical framework proposed by Dai et al. (2004) and Ritzi and Allen-King (2007), this study divided the lithofacies types of the aquifer into two hierarchical scales (Fig. 2). Under this classification, Scale I and Scale II are used in a relative sense to denote two levels of heterogeneity representation within the same aquifer: Scale I resolves finer lithofacies variability, whereas Scale II represents a coarsened description in which Scale I facies are aggregated into composite units that preserve the dominant architectural organization. At Scale I, eight mutually exclusive lithofacies were defined based on permeability contrasts: gravel (G), coarse sand (CS), medium-coarse sand (MCS), medium sand (MS), medium-fine sand (MFS), fine sand (FS), sub-sandy loam (SL), and clayey loam (C). At Scale II, these were aggregated into three composite units: gravel coarse sand (GCS), medium sand (MFS) and sandy clay (SC). Detailed delineations of the lithofacies at Scale I and Scale II are provided in Figs. S1–S8 and S9–S16 in the Supplement, respectively. From a sedimentological perspective, these sediments are interpreted as a river-alluvial-floodplain system associated with the Nen River. Coarse gravel/sand bodies represent channel zone (or paleochannel) deposits, while fine sand containing clay-silt lenses represents floodplain and overflow deposits. Statistical analysis of borehole lithofacies characteristics also revealed a lateral grain size reduction towards downstream from the river, and a prevalent vertical structure of coarse sand overlaying fine sand in the aquifer.

Figure 2 Schematic diagram of lithofacies composition of the layered structural model of Qiqihar phreatic aquifer.

The heterogeneous architecture models were produced via conditional stochastic simulation method under a transition-probability Markov-chain scheme, constrained by 57 boreholes and parameterized by lithofacies volume proportions and mean lengths. The volume proportions (P) define the stationary occurrence probabilities, while the directional mean lengths (L_X, L_Y, L_Z) quantify correlation scales of each lithofacies along the principal directions. To bias-correct parameter estimates affected by incomplete exposure of sections (see Fig. S17 in the Supplement), a Bayesian updating scheme (Dai et al., 2005; White and Willis, 2000) and pixel-based image analysis were employed to calculate lithofacies lengths and volume proportions at both hierarchical scales. Comparison of the volume proportions computed from whole-sections (P_S) with those calculated from boreholes along the sections (P_D) revealed a close agreement (see Table S1 in the Supplement), indicating that these sections are sufficiently representative for estimating structural parameters.

The lithofacies volume proportions and mean vertical thicknesses (L_Z) were derived from discrete borehole data, whereas mean horizontal lengths along x- and y- directions (L_X, L_Y) were obtained from the cross-sections. Statistics of heterogeneous parameters are listed in Table 1. Define the horizontal anisotropy coefficient as $\mathit{\alpha }=L_Y/L_X$, and the vertical anisotropy coefficient as $\mathit{\epsilon }=L_Z/$ [($L_X+L_Y)/\mathrm{2}$]. According to Table 1, most lithofacies are horizontally isotropic at both scales, however, the MCS and SL lithofacies extend longer in the x- direction and the C lithofacies extend longer in the y- direction at Scale I. The mean anisotropy ratio decreases from ∼ 0.005 (Scale I) to ∼ 0.003 (Scale II), evidencing reduced apparent heterogeneity with model coarsening.

Table 1 Statistics of heterogeneous architecture parameters at Scale I and Scale II.

In this study, the 57 boreholes provide hard conditioning data for lithofacies occurrence and aquifer thickness, whereas the eight cross-sections supply additional structural constraints on lateral continuity and stratigraphic organization along and across the principal directions. This study does not aim to deterministically reproduce specific in-situ plumes, but rather quantifies how hierarchical sedimentary architecture and associated parameter uncertainties govern basin-scale dispersion under field-representative flow conditions. The resulting heterogeneous model was intended to be statistically representative, while local connectivity in data-poor areas was treated as uncertain and quantified through a set of conditional realizations. Although boreholes and corresponding cross-sections are more densely packed in the central and western parts of the study area and relatively sparse in the eastern area, this is sufficient to serve the objectives of this study.

2.2.2 Hydraulic conductivity

Hydraulic conductivity (K) was constrained using a combination of pumping-test analyses and grain-size–based empirical equations. A total of 45 pumping experiments were collected in the study area (including steady flow pumping tests and unsteady flow pumping tests) and the corresponding K values were obtained using analytical solutions and type-curve matching. However, the resulting K values are most representative for coarse-grained media. For fine-grained sediments, K was estimated from empirical equations. The soil samples collected from the shallow part of the study area were first analyzed for grain-size distribution, and K values were then estimated using the empirical equations summarized by Vuković and Soro (1992):

$$\begin{array}{}\text{(1)}& K={\displaystyle \frac{g}{v}}C\mathit{\phi }\left(n\right)d_e^{\mathrm{2}}\end{array}$$

where: K is the conductivity (m d⁻¹); v is the kinematic viscosity (m² s⁻¹); g is the gravity, taking the value of 9.81m² s⁻¹; C is the empirical coefficient (provided in Table S2 in the Supplement); φ(n) is the dimensionless porosity function, in which $n=\mathrm{0.225}\times (\mathrm{1}+\mathrm{0.83}{)}^{\mathit{\eta }}$, $\mathit{\eta }=d_{\mathrm{60}}/d_{\mathrm{10}}$, usually referred to as coefficient uniformity; d_e is the particle diameter corresponding to e % finer on the cumulative grain-size distribution curve.

Grain-size analyses showed that all samples had characteristic diameters of less than 0.2 mm (data from Dai et al., 2022). By comparing the K results calculated by various empirical formulas, the USBR method was finally adopted for the lithofacies SL. For the very low-permeability C unit, where grain-size methods are less reliable, K was inferred from a regional Plasticity Index-K (PI-K) empirical relationship, where the PI data was collected from clayey soils in the Qiqihar area. Through the above methods, the statistical results of the mean and variance of the logarithmic conductivity (Ξ=ln (K)) of each lithofacies on two scales were obtained, as shown in Table 2.

Table 2 Univariate statistics of conductivity for different hierarchical lithofacies at different scales.

2.3 Construction of a heterogeneous structural model

Taking into account the aquifer structure and the topographic and geomorphological features of the study area, the simulation domain is set to be about 20 km in the x-direction, 22 km in the y-direction, and 50 m in the z-direction. Note that the purpose of this study is not aimed at accurately simulating groundwater dynamics or pollution plumes in the study area, but rather to discuss the influence of sedimentary structure and related parameters on solute dispersion. Therefore, for the convenience of modeling and parameter setting, the simulation area was simplified into an approximate square. As listed in Table 1, the SL facies type has the smallest mean horizontal extension (640.32 m), while the MS facies type has the smallest mean vertical extension (0.98 m). To capture lithological heterogeneity in detail, the grid cells were sized at 100 m × 100 m × 0.25 m, resulting in a total of 200 × 220 × 200 cells.

Indicated kriging method based on Markov chain transfer probabilities was used to complete the modeling of heterogeneous architectures. Following Proce et al. (2004) and Ren et al. (2022), the sedimentary architecture was modeled in two stages. First, 50 realizations at Scale II and 50 realizations at Scale I were generated based on proportion and length statistics, respectively. This number of realizations is sufficient to obtain stable ensemble statistics (Zhou et al., 2018; Henri et al., 2020). Second, Scale I lithofacies were then mapped onto their corresponding Scale II domains to yield hierarchical models (hereafter referred to as the multiscale models). The resulting multiscale models thus embed fine-scale lithofacies heterogeneity within a framework that retains large-scale sedimentary architecture, enabling hierarchical coupling across scales. As an example, the three-dimensional and two-dimensional cross-sections of a particular heterogeneous model are shown in Fig. 3. In all heterogeneous-architecture models, x=0 is adjacent to the Nen River; larger x indicates greater distance from the river. Larger y denotes upstream, and smaller y denotes downstream.

Figure 3 An example of multiscale model: (a) 3-D facies, (b) 2-D sections; and Scale II model: (c) 3-D facies, (d) 2-D sections. Section lines are located at x=600, 3000, 7600, 11 000 m and y=8000, 12 000, 16 000, 18 000 m.

Facies structures in the simulated 2-D sections closely match field patterns and exhibit pronounced stratification. Along the groundwater flow direction, dominant lithofacies transition from G and MCS to MS and MFS. Vertically, the succession fines downward, with coarser textures near the top and finer ones at the bottom.

2.4 Flow and solute transport simulation

2.4.1 Hydrogeologic conceptual model

Long-term groundwater level observations indicated that water levels had remained relatively stable over the years, suggesting a stable flow regime. According to multi-year groundwater level records, specified-head boundaries at y=0 m andy=22 km (corresponding to the upstream and downstream boundaries) with multi-year mean heads of 146 and 140 m were set, respectively (Fig. 4). The Nen River in the western part of the study area was also generalized as a steady head boundary (144 m). The eastern boundary was specified as a flux boundary, with flux values calculated from the natural hydraulic gradient. The bottom of the phreatic aquifer is an aquitard with overflow discharge, so it was set as the flux boundary. The top boundary is the phreatic surface where recharge is primarily from precipitation, while discharge occurs through evaporation and artificial pumping.

Figure 4 Conceptual hydrogeologic model of the study area and schematic diagram of solute plume planar source release conditions. The left side corresponds to the Nen River, and the pollution source is oriented approximately perpendicular to the regional groundwater flow direction.

2.4.2 Groundwater flow simulation

In this study, MODFLOW-2005 (Harbaugh, 2005) was used to calculate the groundwater transport process in saturated porous media. After constructing the heterogeneous sedimentary architecture, stochastic K fields were generated by assigning facies-conditioned K values to each cell based on the spatial distribution characteristics and statistics ($\overline{\mathrm{\Xi }}$ and ${\mathit{\sigma }}_{\mathrm{\Xi }}^{\mathrm{2}}$). Groundwater-flow simulations were run on the same grid as the architecture, ensuring a one-to-one correspondence between lithofacies and K. Based on the recharge conditions of phreatic water, the lithology and thickness of the vadose-zone, city land and farmland, four water-balance subregions were delineated (shown in Fig. 1). Detailed water-balance equations and parameter values are provided in the source and sink calculation section of the Supplement.

2.4.3 Solute transport simulation

The random walk particle tracer model program RWHet (LaBolle et al., 1996) was used to simulate the solute transport process. This study focuses on the transport characteristics of nonreactive solutes, and the governing equation for solute transport in saturated porous media defined as:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}& \left[\mathrm{\Theta }\right(\mathit{x},t\left)c\right(\mathit{x},t\left)\right]=-{\sum }_{i=\mathrm{1}}^{\mathrm{3}}{\displaystyle \frac{\partial }{\partial x_i}}\left[v_i\right(\mathit{x},t\left)\mathrm{\Theta }\right(\mathit{x},t\left)c\right(\mathit{x},t\left)\right]\\ & +{\sum }_{i,j=\mathrm{1}}^{\mathrm{3}}{\displaystyle \frac{\partial }{\partial x_i}}\left[D_{ij}\right(\mathit{x},t\left)\mathrm{\Theta }\right(\mathit{x},t\left){\displaystyle \frac{\partial c(\mathit{x},t)}{\partial x_j}}\right]\\ & +{\sum }_k{q}_k(\mathit{x},t)c_k(\mathit{x},t){\mathit{\delta }}_k(\mathit{x}-{\mathit{x}}_k)\end{array}\end{array}$$

where c [M L⁻³] is the dissolved resident concentration; v [L T⁻¹] is the velocity; Θ [L³ L⁻³] is the porosity; x_i,j [L] is the distance along the respective Cartesian coordinate axis; c_k [M L⁻³] is the aqueous phase concentration in the flux q_k [L³ T⁻¹] of water at x_k; and δ_k is a Dirac function. D_ij [L² T⁻¹] is the local-scale dispersion tensor.

In this study, Θ was assumed to be stably isotropic and was set to the value of 0.35. Although setting a constant porosity may lead to deviations in the time required to reach a certain stage and the absolute value of the dispersion index plotted on the time axis, this study, however, emphasized the influence of aquifer structure and K-statistics under consistent settings, where the spatial heterogeneity and connectivity of the corresponding velocity field were not significantly determined by subtle spatial variations in porosity. Another advantage of this choice is to avoid introducing other poorly constrained parameter fields into the model. For the same reason, the influence of local scale dispersion and molecular diffusion coefficient was not considered in this study, and therefore the corresponding dispersion and diffusion coefficients were taken as zero.

Solute transport was simulated under two distinct source scenarios: a point source and a planar source. The point source represents localized releases (e.g., spills, leaks or point discharges), whereas the planar source represents an extended source zone that intersects multiple flow paths. These two scenarios are used to quantify how source dimensions control early-time sampling of heterogeneity and uncertainty in plume metrics. Using the bottom of the aquifer near the Nen River upstream as the origin of the coordinate system, the planar source (shown in Fig. 4), oriented perpendicular to the groundwater flow, was centered at (10 000, 3000, 25) m and extended 14 000 m in the x-direction, 100 m in the y-direction, and 50 m in the z-direction. The point source was centered at (10 000, 3000, 44) m, and extended 200 m in the x-direction, 200 m in the y-direction, and 2.5 m in the z-direction. A continuous NaCl source with a concentration of 800 mg L⁻¹ was imposed based on groundwater samples, with a background concentration of zero. Absorption type boundaries were defined at y=0 m and y=22 km to allow particles exit the simulation domain at the inlet and outlet boundaries, while reflection boundaries were used for all other boundaries to ensure particles remained within the domain. The solute transport simulation was performed over a period of 10 000 d, utilizing the same spatial discretization as the flow model.

The solute transport was measured by the solute concentration moments. According to the definitions of Freyberg (1986), the first moment, normalized by the mass in the solution, represents the position of the solute plume, expressed as the centroid coordinates (x_c, y_c, z_c). The second spatial moment quantifies solute spreading around the centroid, given by the variances (${\mathit{\sigma }}_{xx}^{\mathrm{2}}$, ${\mathit{\sigma }}_{yy}^{\mathrm{2}}$, ${\mathit{\sigma }}_{zz}^{\mathrm{2}}$). In the subsequent analysis, only the longitudinal component along the groundwater flow direction (y-direction) was retained.

3 Results

3.1 Numerical simulation results

3.1.1 Water flow simulation results

Model validation was performed using groundwater levels measured at 15 observation wells in 2020, whose spatial distribution is shown in Fig. 5a. The average simulated heads from 50 realizations were compared with the observed values to evaluate model performance.

Figure 5 (a) Distribution of groundwater level measurement points in 2020; (b) Fitting diagram between simulated water level and measured data.

The simulated water levels show good agreement with the observed values, closely following the 1:1 line. This visual consistency is supported by a relatively small error (RMSE = 0.507 m), indicating that the water flow model reproduced the groundwater dynamics of the study area. This agreement further indicates that the hydraulic conductivity values derived from different methods are reasonable and consistent.

3.1.2 Solute transport simulation results

In stochastic transport modeling, aquifer heterogeneity is represented as independent random field with prescribed statistics. Consequently, transport characteristics are described using ensemble statistics over multiple realizations. Figure 6 shows the temporal evolution of plume center of mass y(t) and macrodispersivity α(t) for multiscale (red) and Scale II (blue) models under point- and planar-source release conditions. Lines denote ensemble means, and shaded envelopes depict 10 %–90 % confidence intervals from 50 stochastic realizations. Convergence tests indicate that 50 realizations are sufficient to achieve stable statistics. Two complementary metrics are employed to quantify macrodispersion (Dentz et al., 2000a, b), which differ by their averaging procedure: the effective dispersion coefficient (D^eff) and the ensemble dispersion coefficient (D^ens). The former is computed by first calculating the moment-based dispersion for each realization and then averaging across realizations. In contrast, D^ens is derived directly from the spatial moments of the ensemble-averaged solute concentration field. Under ergodic conditions or at sufficiently long times, these two measures are expected to converge, but under finite-time and non-ergodic conditions D^ens may include an artificial contribution from realization-to-realization wandering of the plume centroid. In this research, we report the corresponding macrodispersivities, α_eff and α_ens, obtained by normalizing dispersion by the average groundwater velocity.

Figure 6 Temporal evolution of longitudinal (a, c) solute transport distance y_c (t) and (b, d) macrodispersivity α(t) for multiscale (red) and Scale II (blue) models under point-source release (a, b) and planar-source release (c, d) conditions. Solid lines denote ensemble means, and shaded envelopes represent the 10th–90th percentile ranges.

The time derivative of the mass transport distance y_c(t) provides a simple consistency check on the simulated transport behavior. The average transport velocity of the solute plume is approximately 0.058 m d⁻¹ for a point-source release and 0.027 m d⁻¹ for a planar-source release (Fig. 6a and c). The latter is close to the Darcy-based estimate of the mean groundwater velocity, about 0.03 m d⁻¹, calculated from the mean hydraulic conductivity of 12.05 m d⁻¹ (Table 2), regional hydraulic gradient of 1 ‰, and effective porosity of ∼ 0.35 used in this study. This consistency indicates that the flow and transport simulations reproduce the macroscopic solute migration behavior of the study system.

A first-order result is the close agreement between the multiscale and Scale II simulations, particularly under planar-source release. This indicates that basin-scale plume evolution can be captured primarily by the geometry and connectivity of the dominant large-scale lithofacies, whereas finer-scale heterogeneity mainly modulates the details of transport. The source geometry mainly affects early-time sampling of heterogeneity and thus the magnitude of realization-to-realization variability. Under point-source release, the plume initially samples only a limited portion of the K field, resulting in larger fluctuations in transport distance and solute dispersion. In contrast, the planar source intersects a broader set of flow paths from the outset, which narrows the uncertainty envelopes and yields a more stable ensemble-like response. This result is consistent with previous theoretical and numerical studies showing that larger source areas enhance early sampling of heterogeneity, stabilize large-scale transport statistics, and reduce uncertainty in plume behavior (Dagan, 2017; de Barros, 2018). This also confirms that source dimensions remain important even at the basin scale.

Figure 6 also shows that α_eff(t) and α_ens(t) do not fully converge within the 10 000 d simulation period, although the degree of separation varies with source geometry and transport stage. In general, the difference between the two metrics remains more evident under point-source release condition. This incomplete convergence indicates that plume evolution remains in a prolonged pre-asymptotic state under the simulated basin-scale transport conditions. At this stage, the plume's evolution has not yet transitioned into the dispersion-dominated phase, which is characterized by the dominance of local-scale dispersion processes. Compared with previously reported site-scale results at the Borden site, where macrodispersivity approached a quasi-stationary value over a much shorter interval, the present basin-scale system exhibits a substantially slower approach to convergence.

3.2 Uncertainty analysis of solute dispersion

A recent global sensitivity analysis indicated that a small set of geologically interpretable factors, most notably lithofacies volume proportions and in-facies mean hydraulic conductivity, exerts first-order influence on non-reactive solute dispersion across regional to basin scales (Ren et al., 2023). Motivated by these findings, we next explore the individual contribution of each factor to solute dispersion. To ensure a stable and consistent comparison while minimizing noise between realizations, all subsequent simulations were conducted under planar-source releases and with Scale II heterogeneous models. This configuration focuses on the dominant lithofacies architecture, yielding an ensemble-like depiction of field-scale migration. This choice also provides a good chance to effectively narrow uncertainty bands and establish a clear baseline for systematically evaluating the impact of key parameters.

3.2.1 Effect of volume proportions

At Scale II, lithofacies are grouped into three types in decreasing order of permeability: GCS, MFS, and SC, with volumetric proportions of 0.503, 0.385 and 0.112, respectively. Two scenarios are designed by changing the lithofacies volume proportions, while holding all other statistics fixed. In Group A, the proportions were set to 0.33, 0.34, and 0.33, approximating an aquifer with equal volumetric fractions of all lithofacies. In Group B, the mixture was skewed toward the low-permeability unit, with volume ratios of 0.2, 0.3, and 0.5, respectively. From a sedimentological perspective, in fluvial–alluvial systems the areal proportion of coarse deposits (e.g., gravel/sand bodies produced in paleochannel zones) versus floodplain fine deposits can vary substantially at the basin scale, reflecting the coupled effects of stream power and sediment supply, channel migration, floodplain aggradation and development (Bridge, 2009). Accordingly, Group A (near-equal proportions) represents a more mixed and interbedded architecture consistent with frequent channel migration and lithofacies switching, whereas Group B (fine-dominated mixtures) represents a low-energy and/or distal floodplain setting where fine deposits are more prevalent and coarse bodies are more isolated.

Figure 7 shows the simulation results for different scenarios. With the increase of time, the solute transport distance increases proportionally, while the α_eff(t) shows a power function growth trend and approaches a quasi-stable value after roughly 5000 d. Across the tested scenarios, the asymptotic-scale α_eff(t) remains on the order of 170 m, indicating that the basin-scale system remains strongly dispersive over long travel times. This behaviour is consistent with the well-established theory of large-scale dispersion in heterogeneous media, where transport is characterized by an initial non-Fickian regime followed by a transition towards Fickian behaviour at late times (Dagan, 1989; Neuman and Tartakovsky, 2009). Increasing the proportion of GCS enhances plume migration and leads to larger y_c(t) and α_eff(t) values, whereas increasing the proportion of SC suppresses both plume migration and dispersivity. These trends indicate that modifying lithofacies proportions changes the balance between connected preferential pathways and retarding low-K domains, even though the overall basin-scale response remains relatively smooth. Previous theoretical and numerical studies (Fiori et al., 2010; Amooie et al., 2017; Soltanian and Ritzi, 2014; Puyguiraud et al., 2020), have also confirmed that large-scale solute dispersion is governed not by a single conductivity class, but by the combined influence of extreme permeability contrasts.

Figure 7 Effect of lithofacies proportions on longitudinal (a) solute transport distance y_c(t) and (b) macrodispersivity α(t) under planar-source release condition. Group A and Group B represent more mixed and fine-dominated architectures, respectively. Solid lines denote ensemble means, and shaded envelopes represent the 10th–90th percentile ranges.

The uncertainty envelopes in Fig. 7 also show a systematic but moderate response to lithofacies proportion. Increasing the fraction of high-K facies generally amplifies the variability among realizations, whereas increasing the fraction of low-K facies produces more uniform outcomes. However, the differences among scenarios remain much smaller than would typically be expected at the site scale. This suggests that, under basin-scale flow conditions, the effect of composition changes is present but moderated by large-scale spatial averaging.

3.2.2 Influence of the mean value of the conductivity

The influence of conductivity on solute dispersion is primarily based on the difference in the mean K values among lithofacies. In the Qiqihar site, the mean K values are 46.02, 10.34, and 0.12 m d⁻¹ for of GCS, MFS, and SC, respectively, indicating a pronounced permeability contrast. As is well known, K varies widely and is subject to considerable estimation and upscaling uncertainty at field scales. To isolate the role of individual lithofacies, three model groups were designed in which only the mean K of a single lithofacies was increased threefold, while the other two remained unchanged. The choice to expand by three times also takes into account the uncertainty of K at a medium to high level. In Group 1, the mean K of GCS was raised to 138.06 m d⁻¹. In Group 2, the mean K of MFS was increased to 31.02 m d⁻¹ and in Group 3, the mean K of SC was increased to 0.36 m d⁻¹. In all cases, the variance of K was preserved, and the underlying heterogeneous sedimentary architecture remained unchanged. Thus, any changes in dispersion can be attributed to altered inter-facies K contrast and the resulting redistribution of flow among lithofacies. Solute transport simulations for these scenarios, conducted under planar-source release conditions, are presented in Fig. 8.

Figure 8 Effect of the mean K of individual Scale II lithofacies on longitudinal (a) solute transport distance y_c (t) and (b) macrodispersivity α(t) under planar-source release condition. In each scenario, the mean Kis increased threefold while the lithofacies architecture and inter-facies variance remain unchanged. Solid lines denote ensemble means, and shaded envelopes represent the 10th–90th percentile ranges.

Figure 8a–c show the results for solute transport distance, and Fig. 8d–f are the α_eff(t) results. They all show a clear, asymmetric response when the mean K is tripled for a single lithofacies while keeping variance and architecture fixed. Increasing the mean K of the most permeable lithofacies (GCS) accelerates plume migration, raises the asymptotic α_eff(t), and slightly widens the uncertainty bands. This indicates that the most permeable lithofacies has a strong influence on basin-scale plume migration. By contrast, increasing the mean K of MFS does not accelerate plume migration, even though the domain-scale mean conductivity increases (from 12.05 m d⁻¹ of Qiqihar site to 19.49 m d⁻¹). Instead, the transport distance is reduced relative to the baseline case, and the corresponding α_eff curve shows a different evolution toward quasi-stationary behavior. Increasing the mean K of SC produces only a very limited change in both transport distance and α_eff(t). The early and late dynamics of α_eff(t) are also consistent with observations of non-Fickian transport, where temporary suppression and delayed convergence are characteristic of strongly heterogeneous aquifers (Fiori and Dagan, 2000; Dentz et al., 2011). However, more importantly, the spatial averaging effect at the basin scale ultimately restores ensemble-like behaviour.

The contrasting responses are likely related to the spatial roles of the three lithofacies in the modelled architecture. GCS is the most conductive lithofacies and is more influential near the release area, so increasing its mean K reinforces pre-existing fast routes and directly promotes plume migration. In contrast, MFS is less dominant at the release location and becomes more important farther along the transport domain. The mean K perturbation therefore reduces the contrast with GCS and redistributes part of the groundwater flux away from the fastest paths, which suppresses early-time plume spreading. Even after a threefold increase, SC remains much less conductive than the other two lithofacies and therefore cannot contribute substantially to the connected high-K transport network. These results show that basin-scale transport cannot be interpreted solely from the change in domain-average conductivity, because the effects of K perturbation depend strongly on which lithofacies is modified and where that lithofacies occurs within the connected architecture. These findings align with previous theoretical and numerical studies, which indicate that solute dispersion at large scales emerges from the interplay between lithofacies contrasts and connectivity, rather than from the mean conductivity values of single lithofacies (Zinn and Harvey, 2003; Soltanian and Ritzi, 2014).

4 Discussion

The present results clarify a central issue raised in the Introduction: under basin-scale flow conditions, the dominant controls on conservative solute dispersion are the geometry and connectivity of the larger-scale lithofacies framework, whereas finer-scale heterogeneity mainly modulates the details of plume evolution. Viewed together, the scenario analyses show that lithofacies proportions and mean K are not separate or competing explanations for this behavior, but two complementary ways in which the same architectural framework governs transport. Changing lithofacies proportions alters the statistical balance between preferential pathways and retarding domains, whereas changing the mean K of individual facies modifies how groundwater flux is partitioned among those pathways. What emerges from these results is that basin-scale dispersion is controlled less by domain-average hydraulic properties than by the way conductivity contrasts are embedded within the connected lithofacies architecture. This interpretation is consistent with previous regional and basin scale studies showing that aquifer-system heterogeneity, lithofacies architecture, and pathway organization remain influential controls on large-scale plume evolution (Carle et al., 2006; Pauloo et al., 2021). A recent global sensitivity analysis across multiscale heterogeneous media shows a robust ranking for non-reactive solute dispersion (Ren et al., 2023): the lithofacies mean K is typically the most influential factor, followed by volume proportions and mean lengths, whereas variance and some correlation scales contribute less. These basin-scale results are also broadly consistent with this ranking, especially in showing that lithofacies mean K and proportions exert the clearest influence on basin-scale dispersive behavior.

This study also indicates a markedly delayed approach to quasi-steady behavior at the basin scale. α_eff(t) approached a quasi-steady state only after roughly 5000 d, and α_eff(t) and α_ens(t) did not fully converge within the entire 10 000 d simulation period. This persistent separation indicates that plume evolution remains strongly pre-asymptotic over long times and distances. Compared with the Borden tracer test study, where dispersivity approached quasi-stationary behavior after around 400 d (Ren et al., 2022), the present basin-scale system shows a substantially slower approach to convergence. Moreover, site-scale studies have also shown stronger sensitivity of plume behavior and predictive uncertainty to realization-specific pathways and source sampling (Zheng et al., 2011; de Barros and Dentz, 2016), whereas the basin-scale response here is more muted. These differences reflect not only domain-scale spatial averaging, but also the gradual way in which the plume samples the heterogeneous flow field and the delaying effect of multiple overlapping pathways on the emergence of ensemble-like transport. Within this setting, source geometry mainly influences early-time uncertainty by determining how much heterogeneity is sampled at the outset, while the long-term structure of plume evolution remains organized by the dominant large-scale architecture. These interpretations are also consistent with previous regional and basin scale studies showing that heterogeneous flow systems can sustain persistent non-Fickian transport over long distances and times, with a delayed approach to simplified or asymptotic transport behavior (Guo et al., 2019; Pauloo et al., 2021).

Combined with earlier laboratory sandbox experiments and site-scale studies developed within the same hierarchical framework (Ma et al., 2022, 2025), the present basin-scale results help refine how the role of heterogeneity should be understood across sedimentary scales. At local and site scales, conductivity contrasts, lithofacies geometry, and connectivity are often expressed more directly through realization-sensitive plume spreading, stronger source-size dependence, and more immediate pathway control on uncertainty (Dai et al., 2004; Ramanathan et al., 2010; Yin et al., 2020). At the basin scale, however, these same controls remain operative, but their expression is increasingly filtered by long travel distances, pathway superposition, and spatial averaging, so that plume evolution becomes more slowly convergent and less sensitive to fine-scale structural differences. In this sense, the controlling mechanisms are not replaced at larger scales, rather, they are reorganized through scale-dependent averaging. In short, sedimentary heterogeneity exerts a scale-persistent but scale-reorganized control on solute transport.

This broader view also holds practical implications for model construction and interpretation. For problems focusing on bulk plume transport at regional to basin scales, it is not necessarily required to retain every detail of the fine-scale stratigraphic structure, provided that the dominant geological contrasts, lithofacies proportions, and directional organization of the connected flow framework are retained. This does not imply that finer-scale heterogeneity is irrelevant, on the contrary, it may remain important for early-time transport, near-source prediction, and problems sensitive to local concentration gradients. It must be acknowledged that neglecting porosity variations and molecular diffusion processes in this study may lead to an underestimation of early plume smoothing and lateral mixing, potentially delaying a significant convergence to Fick behaviour. However, at the basin scale and in the long-distance travel considered in this paper, structure-controlled velocity variations are expected to dominate the dispersion response, therefore, the main conclusions regarding the roles of lithofacies proportions, conductivity contrasts, and connected architecture remain unchanged. In that sense, the findings of this study support a scale-aware modeling strategy, in which the level of model complexity should be matched to the transport behavior of interest, rather than being increased indiscriminately across all scales. Such a scale-aware perspective provides a more practical theoretical basis for model simplification, monitoring design, and uncertainty assessment in data-limited regional aquifer systems.

5 Conclusions

This study provides a comprehensive analysis of multiscale heterogeneity and its influence on non-reactive solute dispersion and modeling uncertainty within the Nen River Basin aquifer system. By integrating a hierarchical architectural model with numerical simulations and cross-scale validation, we have drawn several key conclusions that advance the understanding of transport in multiscale heterogeneous media.

First, basin-scale solute dispersion is governed primarily by the geometry and connectivity of the dominant large-scale lithofacies framework, whereas finer-scale heterogeneity mainly modulates the details of transport. The evolution of dispersivity in our simulations is characterized by a long pre-asymptotic phase, which contrasts sharply with the rapid stabilization observed at sites like Borden. This unique behaviour highlights that at the basin scale, transport remains in a non-ergodic state for quite a long period of time, requiring a new perspective on long-term plume dynamics that accounts for this prolonged pre-asymptotic regime.

Second, the study identifies a distinct buffering effect in basin-scale transport. where long travel distances and numerous superimposed pathways dampen the variability caused by local heterogeneity. Although source geometry, lithofacies proportions, and inter-facies hydraulic contrasts all influence plume evolution, their effects on uncertainty are moderated by long travel distances, pathway superposition, and spatial averaging. Source dimensions remain important for early-time uncertainty, but their influence becomes less pronounced as transport integrates progressively larger parts of the heterogeneous domain. Third, the results support a scale-aware interpretation of solute transport in heterogeneous aquifers. Sedimentary heterogeneity exerts a control on transport that is persistent across scales but reorganized in its expression at the basin scale. For regional- to basin-scale prediction, it is therefore most important to preserve dominant geological contrasts, lithofacies proportions, and connected architectural organization. This provides a more practical basis for model simplification, monitoring design, and uncertainty assessment in data-limited regional aquifer systems.