HESS Opinions: Chemical transport modeling in subsurface hydrological 1 systems – Space, time, and the “holy grail” of “upscaling”

. 13 Extensive efforts over decades have focused on quantifying chemical transport in subsurface 14 geological formations, from microfluidic laboratory cells to aquifer field scales. Outcomes of 15 resulting models have remained largely unsatisfactory, however, largely because domain 16 heterogeneity – characterized for example by porosity, hydraulic conductivity and geochemical 17 properties – is present over multiple length scales, and “unresolved”, practically unmeasurable 18 heterogeneities and preferential pathways arise at virtually every scale. While spatial averaging 19 approaches are effective when considering overall fluid flow, wherein pressure propagation is 20 essentially instantaneous, purely spatial averaging approaches are far less effective for chemical 21 transport essentially because well-mixed conditions do not prevail. We assert here that an explicit 22 accounting of temporal information, under uncertainty, is an additional but fundamental 23 component in an effective modeling formulation. As an outcome, we further assert that “upscaling” 24 of chemical transport equations – in the sense of attempting to develop and apply chemical 25 transport equations at large length scales, based on measurements and model parameter values 26 obtained at significantly smaller length scales – can be considered an unattainable “holy grail”. 27 Rather, we maintain that it is necessary to formulate, calibrate and apply models using 28 measurements at similar scales of interest. 29 31

There have been extensive efforts over the last ~60 years to model and otherwise quantify fluid 36 flow and chemical (contaminant) transport in soils and subsurface geological formations, from 37 millimeter-size, laboratory microfluidics cells to aquifer field scales extending to hundreds of 38 meters and even tens of kilometers.
(1) Basic structural information on "conducting elements" in a system representing a porous 126 and/or fractured geological domain can provide insight regarding overall fluid conduction 127 in the domain, as a function of "conducting element" density. We emphasize that without 128 direct simulation of fluid flow in such a system, this type of analysis does not delineate the 129 actual flow field and velocity distributions throughout the domain. 130 (2) Spatial information on the hydraulic conductivity distribution at a continuum scale, or solid 131 phase distribution at the pore scale, throughout the domain, can be used to determine the 132 flow field. We then show that this is insufficient to define chemical transport. 133 (3) Temporal information on chemical species migration, which quantifies distributions of 134 retention and release times (or rates) of chemicals by advective-dispersive-diffusive and/or 135 chemical mechanisms, can be used to determine the full spatial and temporal evolution of 136 a migrating chemical plume, either by solution of a transport equation or use of particle 137 tracking on the velocity field. 138 We comment, parenthetically, that in conceptual-philosophical terms, this hierarchy and the  Analysis of the geometry of structural elements in a domain can yield basic insights on fluid flow 152 patterns. This approach is used, for example, when examining fracture networks in essentially 153 impermeable host rock. As discussed below, however, full delineation of the underlying velocity 154 field ultimately requires solution of equations for fluid flow. 155 In this context, percolation theory (Stauffer and Aharony, 1994) is particularly useful in 156 determining, statistically, whether or not a domain with N "conducting elements" (e.g. fractures) 157 includes sufficient element density to form a connected pathway enabling fluid flow across the 158 domain. One can estimate, for example, the critical value, Nc, for which the domain is "just" 159 connected, as a function of fracture length distribution, or the critical average fracture length as a 160 function of N needed to reach domain connectivity (Berkowitz, 1995). Similarly, percolation 161 theory shows how the overall hydraulic conductivity of the domain scales as the number of 162 conducting elements, N, relative to the Nc critical number of conducting elements required for the 163 system to begin to conduct fluid. Percolation theory also addresses diffusivity scaling behavior of 164 chemical species. But, fundamentally, percolation is a statistical framework suitable for large 165 ("infinite") domains, and provides universal scaling behaviors with no coefficient of equality; see 166 e.g., Sahimi (2021) for detailed discussion. 167 Other approaches have been advanced to analyze domain connectivity, for example using 168 graph theory and concepts of identification of paths of least resistance in porous medium domains 169 (e.g. Rizzo and de Barros, 2017), or topological methods (e.g. Sanderson and Nixon, 2015). Like  Why is knowledge only of the geometrical "static" structure (spatial distribution of solid phase) Clearly, then, except in highly idealized and simplified geometries, use of a purely analytical 208 solution to identify the full velocity field and streamline patterns at the pore scale is not feasible.   Considering now continuum-scale domains, but in analogy to the example shown in Sect. 2.1, we 239 illustrate why knowledge only of the geometrical "static" structure is insufficient to know the flow 240 dynamics, without solution of the Darcy equation. Here, the geometrical structure refers to the 241 spatial distribution of the hydraulic conductivity, K. and boundary conditions. [Note, too, that critical path analysis from percolation theory (discussed 262 in Sect. 2) -again from purely "static" information without solution of the flow fieldyields an 263 incorrect interpretation, as shown in detail by Edery et al. (2014).] 264 We emphasize that the delineation of "preferential flow paths" is usually relevant only for 265 study of chemical transport; if water quantity, alone, is the focus, then specific "flow paths" 266 travelled by water moleculesand their advective and diffusive migration along and between 267 streamlines, and into/out of less mobile regionsare of little practical interest. On the other hand, 268 the movement of chemical species, which experience similar advective and diffusive transfers, 269 must be monitored closely to be able to quantify overall migration through a domain. We return to  are seen to contain lower K "bottlenecks". 308 309 It is significant, too, that fluid flow and chemical transport occur in preferential pathways that 310 contain low conductivity sections (indicated by circles in Figs. 2c,d). How do we explain passage 311 through low hydraulic conductivity "bottlenecks" within the preferential pathways, rather than 312 migration "only" through the highest conductivity patches? 313 To address this question, we first consider what happens in a 1d path. Consider two paths, each 314 containing a series of five porous medium elements or blocks, with distinct hydraulic conductivity 315 values, Ki. Consider Path 1, with a series hydraulic conductivity values of 3, 3, 3, 3, 3, and Path 2, 316 with values 6, 6, 1, 6, 6 (specific length/time units are irrelevant here). The value of K = 1 317 represents a clear "bottleneck" in an otherwise higher K path than that of Path 1. In a 1d series, 318 however, the overall hydraulic conductivity (Koverall) of the path is given by the harmonic mean of 319 the conductivities of the elements comprising the path: Koverall = 5 / (i=1,5 1/Ki); significantly, in 320 the two cases here, both paths have Koverall = 3. So a "bottleneck" (K=1) can be "overcome" and 321 does not cause necessarily a potential pathway to be less "desirable" than a pathway without such 322 "bottlenecks". In other words, flow through pathways containing some low K regions should be 323 expected. Of course, in 2d and 3d systems, patterns of heterogeneity and pathway "selection" by 324 water/chemicals are significantly more "complicated", but the principle discussed here for 1d 325 systems still holds, in the sense that lower hydraulic conductivity ("bottleneck") elements can (and    We point out, too, that for both pore-scale and continuum-level scenarios, one can solve 350 explicitly a governing equation for transport. Alternatively, one can obtain an "equivalent" 351 solution by solving for "particle tracking" of transport along the calculated streamlines, in a 352 Lagrangian framework. In other words, particle tracking methods essentially represent an 353 alternative means to solve an (integro-)partial differential equation for chemical transport; such 354 methods can be applied, too, when the precise partial differential equation is unknown or the 355 subject of debate. We also note that solution of the relevant equations for fluid flow and chemical 356 transport is sometimes achieved by (semi-)analytical methods, if the flow/transport system can be 357 treated sufficiently simply (e.g. as macroscopically, section-averaged 1d flow and transport in a 358 rectangular domain). 359 We first discuss principal features of pore-scale (Sect. 3.1) and continuum-scale (Sect. 3.2) 360 chemical transport, and in Sect. 3.3, we focus on effective model formulations. We focus on   Here, the macroscopic Pe is based on the mean particle velocity 380 and mean particle displacement distance per transition (or "step"). profiles, but rather asymmetrical, "heavy-tailed" profiles.

479
At this juncture, note that here and below we use the terms "non-Fickian", or "anomalous" -480 others sometimes use the terms "pre-asymptotic" or "pre-ergodic"to denote any chemical 481 transport behavior that differs from that described by the classical ADE or similar type of  It is in this context that the term "homogeneous" packing used above is placed in quotation 491 marks, to indicate that in natural geological media, "homogeneity" does not really exist. Any 492 natural geological sample of porous medium contains multiple scales of heterogeneity; and at each 493 particular scale of measurement, "unresolved" heterogeneities that are essentially unmeasurable 494 are present. And thus, as seen in Fig. 4 for example, the overall transport pattern even in an 495 "homogeneous" system can be non-Fickian (anomalous). We therefore emphasize that because 496 natural heterogeneity in geological formations occurs over a broad range of scales, "normal" 497 (Fickian) transport tends to be the "anomaly", whereas "anomalous" (non-Fickian) transport is 498 ubiquitous, and should be considered "normal".   modeling of chemical transport is more contentious, the reasons for which we expand upon below. 525 We argue here that modeling of chemical species transport requires us to think in terms of time, 526 not just space. To assist the reader to enter this frame of thinking, and to sharpen our Brownian motion (Einstein, 1905). Prior to his analysis, researchers appliedwith puzzlement -  (1) and (2).

585
The debate in the literature between "effective" and high-resolution hydrogeological modeling, 586 as well as various preconceptions and misconceptions discussed below and in Sect. 4, lead 587 naturally to consideration of the (often incorrectly invoked) argument that "fewer model 588 parameters is better". 589 We first discuss briefly aspects of high-resolution hydrogeological modeling in Sect. 3.3.1, 590 and then focus on "effective" transport equation modeling in Sect. 3.3.2. We emphasize that the 591 latter approach is applicable to both small-and large-scale domains. The former approach is 592 generally intended for large-(field-)scale systems, although it is in some sense often applied for 593 detailed pore-scale modeling; this approach is not particularly contentious, per se, but is hampered 594 by the complexity and cost associated with the demand for highly detailed hydrogeological 595 information. Therefore, research work remains heavily invested in "effective" transport equation  This approach is attractive in terms of the ability to "reproduce" detailed heterogeneous 631 hydraulic conductivity structures, and can provide useful "overall assessments" of fluid flow and 632 chemical transport pathways, and migration of a chemical plume. Moreover, solutions for fluid 633 flow and chemical transport can be considered "exact", at least at the scale at which the domain is 634 discretized; they can thus also capture at least some aspects of non-Fickian transport. But even at 635 this type of spatial resolution, the ability to effectively quantify actual chemical transport, even 636 relative to the limited available field measurements, remains a question of debate; the research 637 community, as well as practicing engineers, still often prefer to analyze chemical transport in a 638 domain by use of relatively simple (often 1d, section-averaged) model formulations.

639
Finally, we point out that in the context of efforts to obtain increasing amounts of structural 640 and hydrological information at a given field site, due consideration should also be given to the

685
To explain this approach, we refer to the continuous time random walk (CTRW) framework, 686 which is particularly broad and general (Berkowitz et al., 2006). Significantly, and conveniently, 687 it turns out that special, or limit, cases of a general CTRW formulation lead to other well-known 688 formulations that can also quantify various types of non-Fickian transport, as explained in, e.g., 689 Dentz and Berkowitz (2003) and Berkowitz et al. (2006). These "subsets" include mobile- inappropriate, for example, to apply a mobile-immobile (two domain) model to interpret chemical 697 transport in a "uniform, homogeneous" porous medium when it displays non-Fickian transport 698 behavior (recall Fig. 4).

699
Here, we describe only briefly the principle and basic aspects of the CTRW formulation; 700 detailed explanations and developments are available elsewhere (e.g. Berkowitz et al., 2006).

701
To introduce "temporal thinking" in the context of non-Fickian transport, we begin by 702 mentioning the analogy between a classical random walk (RW) -which leads to Fick's lawand 703 the CTRW. A classical random walk is given in Eq. 1: where (ℓ, ℓ ′ ) represents the probability of a random walker ("particle") advancing from location 708 ℓ′ to ℓ, (ℓ ′ ) denotes the probability of a particle being located at ℓ′ at (fixed) time step n, and 709 +1 (ℓ) denotes the probability of the particle then being located at ℓ at step n+1. With this formulation in mind, Einstein (1905) and Smoluchowski (1906a,b) demonstrated that for n 711 sufficiently large and a sufficient number of particles undergoing purely (statistically) random 712 movements in space, the spatial evolution of the particle distribution is equivalent to the solution 713 of the (Fickian) diffusion equation. This elegant discovery demonstrated that a partial differential 714 equation and its solution can be represented by following, numerically, the statistical movement 715 of particles (i.e. particle tracking) following a random walk. Remarkably, random walk 716 formulations are "generic" in the sense that they can be applied in a broad range of phenomena in 717 physics, chemistry, mathematics, and life sciences; here, they describe naturally migration of 718 chemical species (dissolved "particles" or "packets") in water-saturated porous media.

719
Generalizing the partial differential equation to include transport by advection, solution of the 720 ADE under various boundary conditions can then be determined by an appropriate random walk 721 method.

722
The simple random walk given in Eq. 1 can be generalized by accounting for time, replacing where +1 ( , ) is the probability per time for a particle to just arrive at site s at time t after n+1 730 steps and (s, t) is the probability rate for a displacement from location s' to time s with a difference 731 of arrival times of t-t'. It is clear that (s, t) is the generalization of (ℓ, ℓ ′ ) in Eq. 1, and that the 732 particle steps can each now take place at different times. Indeed, it is precisely this explicit 733 accounting of a distribution of temporal contributions to particle transport, not just spatial 734 contributions, that offers the ability to effectively quantify transport behaviors as expressed by, 735 e.g. heavy-tailed, non-Fickian particle arrival times.

736
To where does the generalization in Eq. 2 lead us? In a mindset similar to that of Brownian where v and D are generalized particle velocity and dispersion, respectively, and M(t) is a 766 temporal memory function based on (t).

767
The strength of this type of formulation is that it effectively quantifies (non-Fickian) early

776
Equation (4) is essentially an ADE weighted by a temporal memory. When (t) is an 777 exponential function (or power law but for β  2), M(t)  (t) and we recover Fickian transport 778 described by the ADE; thus, the ADE assumes, implicitly, that particle transition times are 779 distributed exponentially. But with a power law form ( ) ~ −1− for 0 <  < 2, the transport is The term "modeling" is used in many contexts and with differing intents. However, in the 841 literature dealing with chemical transport in subsurface hydrological systems, there are frequent 842 but often misguided "arguments" regarding "which model is better", with a major point of some 843 authors being the claim that "fewer parameters is always best". Not always. Indeed, some models 844 involve more parameters than others, but if these parameters have physical meaning and are needed 845 as factors to quantify key mechanisms, then "more parameters" is not a "weakness". We  It should be recognized that quantitative model information criteria, or model selection criteria, 863 can be used to assess and compare various model formulations that are applied to diverse scenarios 864 (such as fluid flow, chemical transport) in subsurface geological formations. These information 865 criteria include AIC (Akaike, 1974), AICc (Hurvich and Tsai, 1989), and KIC (Kashyap, 1982) 866 measures, as well as the Bayesian (or Schwarz) BIC (Schwarz, 1978). They are formulated to rank 867 models, or assign (probabilistic) posterior weights to various models in a multimodel comparative 868 framework, and therefore focus on model parameter estimates and the associated estimation 869 uncertainty. As such, these information criteria discriminate among various models according to 870 (i) the ability to reproduce system behavior, and (ii) the structural complexity and number of 871 parameters. Discussion of theoretical and applied features of these criteria is given elsewhere (e.g.  To conclude this section: Notwithstanding the above arguments, some readers might continue 880 to argue that the approach discussed hereviz., the need for time considerations as well as space 881 such as embodied in the CTRW framework and related formulations -is "inelegant" because it 882 requires more parameters relative to the classical ADE. In response, the reader is encouraged to 883 recall the words of Albert Einstein following criticism that his theory of gravitation was "far more 884 complex" than Newton's. His response was simply: "If you are out to describe the truth, leave 885 elegance to the tailor". We begin by defining the term "upscaling" in the context of the discussion here on chemical significantly smaller length scales. 895 We attempt "upscaling" in the hope of developing governing equations for chemical transport 896 at larger and larger scales, from pore, to core, to plot, and to field length scales. Clearly, then,

897
"upscaling" is relevant to the modeling approach discussed in Sect. claiming that "parameters are empirical because they are estimated by calibration (fitting) to 927 experiments"; additional "criticisms" follow, for example, that such as a model is therefore not 928 "universal", and/or "it therefore has no predictive capability". We address these latter "criticisms"  With regard to model "universality", recall that, for example, percolation theory (discussed at 941 the beginning of Sect. 2) offers "universal" exponents in scaling relationships. But even for this 942 type of convenient and useful, statistical model, such scaling relationships, too, can only advance 943 from "scaling" (e.g. A ~ B) to a full "equation" (e.g. A = kB) by calibration of a coefficient of 944 equality (k) against actual measurements. So even in "simple" models, model calibration cannot 945 be avoided.

946
To address "empiricism"here enters the question of whether parameters of a particular model 947 (in this case, equations for chemical transport) have a physical meaning. As discussed in Sect. Pe distribution, the impact on the overall transport pattern evolution cannot be determined without 985 actually solving for transport in the domain.

986
For chemically reactive species, the transport situation becomes even more complex, because 987 the local residence time, not just the local Pe, must be taken into consideration. Moreover, when 988