the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Statistical characterization of environmental hot spots and hot moments and applications in groundwater hydrology
Jiancong Chen
Bhavna Arora
Alberto Bellin
Environmental hot spots and hot moments (HSHMs) represent rare locations and events that exert disproportionate influence over the environment. While several mechanistic models have been used to characterize HSHM behavior at specific sites, a critical missing component of research on HSHMs has been the development of clear, conventional statistical models. In this paper, we introduced a novel stochastic framework for analyzing HSHMs and the uncertainties. This framework can easily incorporate heterogeneous features into the spatiotemporal domain and can offer inexpensive solutions for testing future scenarios. The proposed approach utilizes indicator random variables (RVs) to construct a statistical model for HSHMs. The HSHM indicator RVs are comprised of spatial and temporal components, which can be used to represent the unique characteristics of HSHMs. We identified three categories of HSHMs and demonstrated how our statistical framework is adjusted for each category. The three categories are (1) HSHMs defined only by spatial (static) components, (2) HSHMs defined by both spatial and temporal (dynamic) components, and (3) HSHMs defined by multiple dynamic components. The representation of an HSHM through its spatial and temporal components allows researchers to relate the HSHM's uncertainty to the uncertainty of its components. We illustrated the proposed statistical framework through several HSHM case studies covering a variety of surface, subsurface, and coupled systems.
Environmental hot spots and hot moments (HSHMs) were originally defined as rare locations or events that support or induce disproportionately high activity levels (e.g., chemical reaction rates) compared to surrounding areas or preceding times (McClain et al., 2003). Vidon et al. (2010) further classified HSHMs into either transportdriven or biogeochemically driven HSHMs, based on the mechanisms causing the HSHMs. Bernhardt et al. (2017) derived the concept of ecological control points (CPs) related to HSHMs, defining CPs as areas of the landscape that exert a disproportionate influence on the biogeochemical behavior of an ecosystem under study. These definitions have mainly focused on HSHMs related to elevated biogeochemical activities triggered by hydrological or biogeochemical processes or a confluence of both processes. The concept of HSHMs is also used in climate science, where it is related to elevated greenhouse gas emissions or specific locations that are subject to extreme natural hazards (e.g., sealevel rise, floods, hurricanes, or earthquakes) caused by climate change (Arora et al., 2021; Shrestha and Wang, 2018). Further, Henri et al. (2015) related HSHMs to locations experiencing elevated environmental risks and developed the incremental lifetime cancer risk (ILCR) model to quantify the effects of hot spots on human health. In the present study, we provide a unified treatment of both positive and negative impacts of HSHMs, which allows us to present an integrative analytical framework for understanding and modeling HSHMs in various fields.
Various approaches have been developed to better quantify HSHM dynamics, including numerical modeling, empirical modeling and databased approaches with statistics. For example, Dwivedi (2017) developed a threedimensional highresolution numerical model to investigate whether organiccarbonrich and chemically reduced sediments located within the riparian zone act as denitrification hot spots. Their study demonstrated a significantly higher potential (∼70 %) of the naturally reduced zones (NRZs) to remove nitrate than the nonNRZ locations. Arora et al. (2016b) used a twodimensional transect model and showed that temperature fluctuations constituted carbon hot moments in a contaminated floodplain aquifer that resulted in a 170 % increase in annual groundwater carbon fluxes. Gu et al. (2012) developed a Monte Carlo reactive transport approach and discovered how denitrification HSHMs are triggered by river stage fluctuations. Abbott et al. (2016) developed the HotDam framework that combines the HSHM concept and Darmköhler number with multiple tracers to advance our understanding of ecohydrology. Statistical concepts have also been used to identify HSHMs through simple comparison to the average, substantial percentage of total flux, outlier in distribution of data, statistically significant difference between or among landscape elements and contribution to flux/total area or time (Bernhardt et al., 2017, and references therein). Wavelet and entropybased approaches have also been used to identify nonuniform regions and times and consequently HSHMs (Arora et al., 2013, 2019a). However, most of these quantitative methods are derived based on sitespecific data, which severely limits the transferability of these approaches. In contrast, a unified HSHM approach offers multiple advantages. First, a unified strategy based on commonly used parameters for a given HSHM would allow modelers to create probability priors that could be used for prediction of said HSHMs at unsampled or poorly sampled sites (Li et al., 2018). Second, such a standardized approach for modeling HSHMs could be beneficial to developing and implementing monitoring standards and regulations for environmentally sensitive HSHMs. Last, but not least, a unified approach can be used together with mechanistic models to capture uncertainty and heterogeneity for HSHMs in environmentally relevant applications.
Successful characterization of HSHMs through deterministic physically based models or purely statistical approaches relies on experts' knowledge of a site, intensive field characterization, and possibly continuous field sampling to provide the data to develop and validate these approaches. Understandably, intensive site characterization and longterm sampling can be quite challenging due to the associated costs and efforts. In this regard, having access to a stochastic approach that could improve predictions through builtin model updating (i.e., Bayesian) capabilities could prove to be an advantage.
Stochastic concepts and models have been widely applied in hydrology and hydrogeology for addressing situations subject to uncertainty, including but not limited to modeling flow and contaminant transport, quantifying subsurface heterogeneity and the associated uncertainties, developing strategies for site characterization, and providing informative priors for ungauged watersheds. Bayesian approaches were found to be particularly useful, especially through concepts such as conditioning and updating. In this paper, we aim to bring the experience gained in hydrology and hydrogeology to HSHM modeling.
An important characteristic of HSHMs is that they occupy a limited portion of the investigated domain and may be active for a limited amount of time since they are activated when the control variable exceeds a given threshold. Physical and geochemical heterogeneities and the impossibility of fully characterizing them render the deterministic identification of HSHMs a vanishing objective. To address this hurdle, we propose to cast the problem into a probabilistic framework by seeking the probability of HSHM occurrence at a given position and time. For a given time and/or space interval and for a priori specified criteria, an HSHM occurrence could be viewed as a binary variable where the ensemble mean is the probability of occurrence. Indicator statistics have previously been applied to model flow and transport phenomena in groundwater (Rubin and Journel, 1991), where indicators were used to model the spatial distribution in a sandshale formation. Wilson and Rubin (2002) and Bellin and Rubin (2004) used indicator statistics to characterize aquifer heterogeneity. These studies suggest that representation of a system's structure through indicator formulation holds the potential to make informed decisions, for example concerning remediation actions, under incomplete site characterization.
Based on the mechanisms that trigger HSHMs, we identified three categories of HSHMs: (1) those triggered only by spatial (static) contributors, (2) those triggered by both spatial (static) and temporal (dynamic) contributors, and (3) those triggered by multiple dynamic contributors. Applications of the proposed indicator formulation to a diverse range of HSHM situations are presented to illustrate the generality of our proposed approach. The remainder of the paper is structured as follows. Section 2 outlines the proposed statistical framework for predicting HSHMs. In Sect. 3, various reported cases from previous HSHM studies are presented using the framework of our proposed approach, intended to demonstrate its generality. In Sect. 4, we present an HSHM application in groundwater hydrology and show how the HSHM uncertainty relates to the spatial variability of the hydraulic conductivity. Advantages and limitations of our approach are discussed in Sect. 5.
Herein, we present a probabilistic formulation of hot spots and hot moments, which considers the HSHM occurrence as a binary event, expressed through indicator statistics embedded with the HSHM underlying physics. Section 2.1 summarizes the indicator formulation of HSHMs. Based on the contributors to HSHMs, we classified HSHMs into three different types, and we demonstrate how indicators are constructed for each type of HSHM in Sect. 2.2–2.4, respectively. Section 2.5 focuses on the linkages between indicators and Bayesian concepts. Case studies for each class of HSHM are provided in Sect. 3.
2.1 Indicator formulation of HSHMs
HSHMs represent intervals in space and/or time characterized by hydrobiogeochemical activity rates or fluxes that differ significantly from the background conditions, thus exerting disproportionate influences over an ecosystem's dynamics. We define Ω^{*} as the volume (subdomain) within which the hot spot may be located and t^{*} the time at which the hot moment occurs. An indicator random variable, ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$, is used to identify whether an HSHM occurs at (${\mathbf{\Omega}}^{*},{t}^{*}$) or not. If userdefined critical conditions needed to trigger an HSHM are met at (${\mathbf{\Omega}}^{*},{t}^{*}$), then ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})=\mathrm{1}$, and it is equal to zero otherwise. These userdefined critical conditions may require concentration or reaction rates above certain thresholds or require both concentration and rate limits depending on the type of HSHM. What makes the indicator a random variable is the uncertainty in the spatial and temporal distribution of the quantities triggering the HSHM event in reallife applications. Following the original definition by McClain et al. (2003), in our method, ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ can take a value of 0 or 1, depending on whether suitable thresholds are exceeded or not as follows:
where $C(\mathit{x},{t}^{*})$ and $R(\mathit{x},{t}^{*})$ are the concentration and reaction rate at position x and time t^{*}, respectively. C_{th} and R_{th} represent the concentration and reaction rate thresholds, respectively, which identify whether an HSHM is triggered or not. Defining indicators with concentration or reaction rate depends on the target of an HSHM. For example, indicators defined by concentration are preferred with transportdriven HSHMs, whereas biogeochemically driven HSHMs may require indicators defined by both concentration and reaction rates. The threshold values can also be based on regulatory limits or defined by the user.
The critical values of C_{th} and/or R_{th} are keys to an effective application of the above framework and should be determined based on the specific scenario under investigation. For example, in the case of contaminants that are associated with significant environmental or health risks (e.g., nuclear waste or a cancerous substance), C_{th}=0 or R_{th}=0 can be used so that the HSHM will be triggered as soon as there is the presence of such contaminants and relevant chemical reactions. As an alternative, a limit in the total accumulated mass or fluxes within hot spots may also be set, such as suggested by the EPA (USEPA, 2001), but in this case the definition (1) of the indicators should be modified. For water quality parameters, C_{th}=MCL or ${R}_{\mathrm{th}}={R}^{*}$ can be assigned, where MCL represents the maximum concentration limit for a specific solute, whereas R^{*} could represent a critical reaction rate. The critical thresholds can be determined based on statistics, such as percentiles and extremes as defined by regulations or analytical studies. Alternatively, C_{th} and R_{th} could also be chosen based on the experts' domain knowledge or from welldocumented studies in similar environments. Through the flexibility to adopt different choices for activation thresholds, our approach could allow users to compare relevant indicator models and assess their applicability by testing how different thresholds would influence the probability of the HSHM occurring and assessing said probabilities against risk tolerance and/or regulations.
Following the definition of ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ as a binary random variable (Eq. 1), we propose to model it with a Bernoulli distribution, such as ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})\sim \mathrm{Bernoulli}(<{I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})>)$, where $<.>$ is the ensembleaveraging operator. An important characteristic of the Bernoulli distribution is that all the statistical moments of ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ can be expressed as a function of the ensemble mean $<{I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})>$. For example, the variance is given by $\mathrm{var}\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},{t}^{*}\left)\right)=<{I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})>\cdot (\mathrm{1}<{I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})>)$. Thus, being able to fold the HSHM physics into an indicator formulation, a simplified approach is presented through Eq. (1).
Casebased formulation of the Bernoulli distribution of ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ requires the incorporation of the mechanisms that govern the development and occurrence of HSHMs into the indicator model. To facilitate this undertaking, we propose to decompose ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ into a typeA (static) indicator random variable, I_{s}(Ω^{*}), and a typeB (dynamic) indicator random variable, ${I}_{\mathrm{d}}({\mathbf{\Omega}}^{*},{t}^{*})$. Definitions of the typeA and typeB indicators are provided herein.

TypeA (static) contributors. This category covers discrete spatial elements (and their associated critical states) that could trigger an HSHM once they come into contact with typeB contributors (see discussion below). Critical states are the range of values needed to trigger an HSHM (either in standalone mode or when coupled with typeB contributors).

TypeB (dynamic) contributors. This category covers dynamic variables (and their associated critical states) that could trigger an HSHM once they come into contact with typeA contributors. This category includes, for example, mass transport variables. It also includes changes in local hydrological and environmental conditions (e.g., water table fluctuations). The displacements of solutes in the subsurface (trajectories and travel times) from below and aboveground processes are prime examples of typeB contributors.
As an example, naturally reduced sediments (typeA contributor) occurring next to the river corridor at the Rifle site were identified as carbon export hot spots (Arora et al., 2016a; Wainwright et al., 2015). Studies showed that these hot spots were triggered when temperature conditions (typeB contributor) varied in the subsurface, resulting in a 170 % increase in groundwater carbon export from the floodplain site to the river (Arora et al., 2016b). In another example, topographic features, such as the backslope of the lower montane hillslope (typeA contributor) within the East River Watershed (Hubbard et al., 2018), were considered denitrification hot spots, which can have a significant impact on the watershedscale nitrogen loss pathway. These hot spots were often triggered by spring snowmelt and storm events (typeB contributor).
Both indicators of the typeA and typeB contributors assume a value of either 0 or 1. If one of these indicators takes a value of 1, it can be viewed as an HSHM contributor. However, for an HSHM to occur, both indicators must have a value of 1 at the same location and time. This idea can be expressed as follows:
Based on the mechanisms of HSHMs, we can distinguish three different HSHM categories as discussed below. These categories can be used to guide the application of the above statistical framework in a variety of complex HSHM scenarios, and they can also be used to develop analytical or numerical solutions for both typeA and typeB contributors.
2.2 HSHMs induced by typeA (static) indicators
In this section, we consider HSHMs that are defined by static indicators only (Fig. 1a). This list can include zones of high, persistent concentration and reactivity that are due to the subsurface or the ecosystem's unique hydrological and biogeochemical properties. For example, the accumulation of contaminants in the subsurface (e.g., the high nuclide concentration in the subsurface at the Hanford site) could lead to the evolution of persistent, highreactivity zones. An aquifer's reactivity is another example that could distinguish certain regions with high reactivity compared to surrounding areas (Loschko et al., 2016). Such highreactivity spots (hereafter denoted as Ω^{*}) can be characterized by static indicator RVs due to the persistence of high concentration or reactivity. The static indicators are defined as follows:
where ${Z}_{\mathrm{s}}^{*}$ represents the conditions needed to trigger a hot spot at Ω^{*}, and Z(Ω^{*}) represents the corresponding local conditions at Ω^{*}. Notice that Ω^{*} is a volume centered at a selected position of the domain where the probability of developing an HSHM is evaluated.
2.3 HSHMs induced by typeA (static) and typeB (dynamic) indicators
HSHMs can also result from dynamic processes encountering specific local conditions at Ω^{*} (Fig. 1b). This is the situation described by Eq. (2), where the typeA indicators are determined first and then used jointly with the typeB indicators for complete HSHM characterization. For example, Bundt et al. (2001) concluded that locations (Ω^{*}) intersected by preferential flow paths are possible biological hot spots for soil microbial activities. Meanwhile, dynamic factors, such as snowmelt or rainfall infiltration control contaminant transport via the preferential flow paths, thus constitute the hot moment component. Additional case studies are presented in Sect. 3.
For HSHMs induced by both typeA and typeB indicators, the static locations are selected first, based on their HSHMrelated properties. After this, we can focus on characterizing the HSHM dynamics as they relate to the relevant locations. A selected location, Ω^{*}, could become an HSHM site based on characteristics defined through the following typeA and typeB indicators, respectively:
where ${Z}_{\mathrm{d}}^{*}$ represents the critical conditions needed to trigger the hot moment, whereas ${Z}_{\mathrm{d}}({\mathbf{\Omega}}^{*},{t}^{*})$ represents the corresponding, criticalstate local conditions at t^{*} and Ω^{*}. The statistical model of ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ can be expressed using the statistical models of I_{s} and I_{d} as shown in Eq. (2).
2.4 HSHMs induced by multiple typeB (dynamic) indicators
A confluence of dynamic processes could result in the formation of an HSHM (Fig. 1c). Unlike the previous scenarios where static locations can be determined through known characteristics provided by geophysical or other types of data, HSHMs can also emerge due to a confluence of dynamic processes. For example, Gu et al. (2012) analyzed how streamflow fluctuations could trigger a nitrogen HSHM. In their example, the dynamics of the streamflow and groundwater controlled the transport and mixing of the chemical reactants, thus triggering the occurrences of HSHMs. For this case, the static locations of Ω^{*} are determined by the confluence of multiple dynamic processes, not being restricted by a set of local conditions. In this case, only typeB indicators need to be modeled.
We can consider the case where an HSHM is predicated on m dynamic processes, d_{j}, where ${I}_{\mathrm{d},j}({\mathbf{\Omega}}^{*},{t}^{*})$ is the dynamic for each dynamic d_{j} at Ω^{*} and time t^{*}. The hotspot location Ω^{*} is determined by the confluence of all dynamic processes at time t^{*}. These dynamic processes are not necessarily independent. Therefore, generally, the statistical model for the comprehensive dynamic indicator (which covers all dynamic contributors) assumes the following form:
In situations where the various dynamic contributors can be viewed as independent (e.g., Destouni and Cvetkovic, 1991) – i.e., where the reactants travel via different paths – then, assuming independence, we can state that
Here, the mean of the dynamic indicator becomes
If Ω^{*} is a hot spot, then Eq. (9) also defines $<{I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})>$. However, if Ω^{*} is not a hot spot, then we need to resort to coupled statistical modeling, as suggested by Eq. (2).
2.5 Additional advantages of stochastic formalism
The statistical framework provides several benefits. Unified formulations of HSHMs through indicators provide us with a platform to evaluate alternative HSHM models thoroughly and objectively. For example, the Akaike information criterion (AIC, Akaike, 1974) and Bayesian information criterion (Schwarz, 1978) can be used to rank between alternative indicator formulations and evaluate their ability to explain HSHM observations. Smaller AIC and BIC values indicate more information preserved in a given indicator HSHM model and imply better model quality than other indicator models. On the other hand, if larger AIC and BIC values are observed, important processes for HSHMs are likely missing, indicating the necessity of increasing site characterization and refinement of conceptual models.
In addition, informative priors constructed from similar HSHM sites (Cucchi et al., 2019; Li et al., 2018) could advance early stage planning for HSHM site investigation. Knowledge from studies at similar HSHM sites can be summarized into prior distributions, which can account for variabilities within and between sites. For example, Cucchi et al. (2019) demonstrated how the distribution of hydraulic parameters at unknown target sites can be predicted using information from hydrologically similar sites with existing tool packages such as exPrior. Goaloriented site characterization also becomes feasible with informative priors; for example, Li et al. (2018) demonstrated the usefulness of informative priors in reducing model uncertainty and potential risks for estimating groundwater drawdown at Mintang tunnel in China. Therefore, through integration with statistical concepts, unified formulations of HSHMs enable us to integrate Bayesian concepts to obtain combined and less risky estimations of HSHMs at new sites, which can help us gain better understanding of the underlying mechanism.
In this section, we selected numerous examples from published research to present how our approach can be used to derive statistical representations for the HSHMs investigated in these studies. We grouped these studies into three categories based on the similarities of their underlying HSHM mechanisms, as described in Sect. 2. Section 3.1 demonstrates the formulation of staticonly HSHM; Sect. 3.2 presents the case with static and dynamictriggered HSHMs, and Sect. 3.3 summarizes the steps to construct multiple dynamic indicators for HSHMs. Table 1 presents a summary of these cases, where environmental risk levels as well as impacts on the ecosystem were also included.
3.1 HSHMs triggered by static contributors only
In this section, we use Wainwright et al. (2015) as an example to illustrate our process to construct ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})$ following Eq. (3), where an HSHM is triggered by static contributors only (Sect. 2.1). NRZs within floodplain environments at the Rifle site are considered biogeochemical hot spots because they represent elevated concentrations of uranium, organic matter, and geochemically reduced minerals, and they have been found to contribute significant carbon fluxes to the atmosphere and to local rivers (Arora et al., 2016a). Due to its characteristics, we considered the spatial distribution of an NRZ to be a staticmechanismbased hot spot. Wainwright et al. (2015) used geophysical data (e.g., induced polarization) to map the distribution of an NRZ at the subsurface level. They found that the phase shift (ϕ) from the induced polarization data of the NRZ was within [4.5, 5] mrad compared to nonNRZ locations at ϕ⊆[1, 3.5] mrad. Thus, ϕ can be used to construct the static indicator with a critical condition of [4.5, 5] mrad. Therefore,
Other static attributes, including but not limited to elevation, hydraulic conductivity, and resistivity, can also be used to define the critical conditions to construct the static indicator for hot spots through Bayesian conditioning.
3.2 HSHMs occurring when dynamic contributors coincide at locations defined by static contributors
The second case we present here utilizes Eqs. (4)–(6), where HSHMs are triggered when dynamic contributors coincide at hot spots determined by static contributors. Here, we present the case investigated by Duncan et al. (2013), where riparian hollows representing less than 1 % of the total catchment area contributed more than 99 % of the total denitrification within the watershed. In their study, the denitrification rates peaked during the base flow (midsummer) period, when the riparian hollows were partially oxygenated and the hydrologic fluxes were at a minimum. The site was considered to have low inorganic N availability and, thus, nitrate was supplied via nitrification. The highest rates of denitrification were therefore tied to nitrification and the partially aerated conditions.
The static indicator needs to be constructed based on the microtopographical features within the riparian zone. Specifically, the topographic wetness index (TWI) (Beven and Kirkby, 1979; Sørensen et al., 2006) was used in Duncan et al. (2013) to delineate the riparian hollows from other riparian locations. Terrain analysis indicated TWI threshold values of 6.0 and 8.0 for riparian hollows under wet and dry conditions, respectively, whereas 4.8 and smaller TWI values corresponded to other riparian locations (e.g., hummocks). Thus, the static indicator can be constructed using the TWI values within the riparian zone to determine the hotspot locations – the hollows. Hence,
Multiple dynamic processes control the denitrification rate at the riparian hollows. As examined by Duncan et al. (2013), a partially aerated condition (${C}_{{\mathrm{O}}_{\mathrm{2}}}>\mathrm{5}$ %) is needed to support nitrification, which supplies the nitrate for denitrification. As quiescent, nonstorm periods during base flow favor the coupled nitrification–denitrification mechanism, this is another key process that needs to be represented by a dynamic indicator. Although Duncan et al. (2013) did not mention specific concentration ranges for nitrogen species, the major components, such as organic N, should be available. Therefore, we can construct the dynamic indicators as follows:
where ${I}_{\mathrm{d},\mathrm{Hydro}}({\mathbf{\Omega}}^{*},{t}^{*})$ is the dynamic indicator representing the streamflow stages; this will be 1 if the base flow conditions are met. Additionally, here, ${I}_{\mathrm{d},\mathrm{N}}({\mathbf{\Omega}}^{*},{t}^{*})$ is the dynamic indicator for the transport of the nitrogen species in the subsurface that support the coupled nitrification–denitrification mechanism.
It is noted that these dynamic processes are not statistically independent. Usually, when one condition is met (e.g., base flow conditions), other conditions may consistently be satisfied (e.g., the transport of nitrogen in riparian hollows). Alternatively, numerical modeling approaches can be used to construct the dynamic indicators based on the critical conditions at riparian hollows (Ω^{*}), where we could directly target N_{2} fluxes using a Monte Carlo approach. The statistical formulation used here is constructed specifically for the mechanisms described by Duncan et al. (2013). Thus, the detailed threshold limits could change under other denitrification HSHM cases, such as the case presented in Hill et al. (2000), who focus on desert landscapes, or the one by Harms and Grimm (2008), where the monsoon season is influential for the nitrogen transport. Nonetheless, the general formulation of HSHMs using indicators is still applicable.
3.3 HSHMs occurring when multiple dynamic processes converge in space
HSHMs can also be triggered by the confluence of multiple dynamic processes that lead to the convergence of complementary reactants at Ω^{*}. Complementary reactants can be mobilized and transported via different hydrologic flow paths. They can converge at hotspot locations and trigger hot moments during the mixing. Following the statistical framework developed in this study, Eqs. (7) to (9) are suitable for this condition. In order to illustrate how the dynamic indicators are constructed, we consider here the case reported by Gu et al. (2012), where high biogeochemical activity was observed at the interface of groundwater and surface water during the stream stage fluctuations, which resulted in significant instream denitrification and ${\mathrm{NO}}_{\mathrm{3}}^{}$ removal.
In their study, hot spots form around the nearstream riparian subsurface during river stage fluctuations, where active biogeochemical reaction (e.g., denitrification) requires both O_{2} depletion and the simultaneous presence of ${\mathrm{NO}}_{\mathrm{3}}^{}$ and the dissolved organic carbon (DOC). Specifically, the spatiotemporal distribution of denitrification hot spots coincides with an O_{2} depletion zone along the DOC infiltration flow paths. In order to determine the mixing of groundwater and surface water during stage fluctuations, Gu et al. (2012) defined bank storage volume V(t) and maximum bank storage volume V_{max}. The flood hydrograph was subdivided into the rising limbs, recession limbs and return flow, the latter representing the slow restitution of part of the water that infiltrated during the previous stages. Considering the different dynamics of these components, they observed that the largest infiltration rate occurred prior to the maximum stage rise, while V_{max}=5 m^{3} m^{−1} (critical condition) occurred in the recession limb of the flood event. Instead, maximum return flow occurred toward the end of the recession curve before the stream hydrograph stabilizes. Maximum ${\mathrm{NO}}_{\mathrm{3}}^{}$ rate removal occurred when the return flow phase was almost complete and then decreased until the depletion of ${\mathrm{NO}}_{\mathrm{3}}^{}$. Through statistical analysis, they found that V_{max}, viewed as an integrated index for hydrological exchange, could explain 64 % of the variation in the ${\mathrm{NO}}_{\mathrm{3}}^{}$ removal. Thus, V_{max} can be used as the critical state to determine whether or not the hyporheic dynamics is significant to enhance relevant biogeochemical processes. In order for the hot moments to be significant, the stream–riparian zone should also be microbially active. Based on these conditions, the dynamic indicators can be constructed as follows:
where ${I}_{\mathrm{d},\mathrm{Hydro}}({\mathbf{\Omega}}^{*},{t}^{*})$ represents the dynamic process induced by the hydrologic conditions (e.g., stage fluctuation), and ${I}_{\mathrm{d},\mathrm{Chem}}({\mathbf{\Omega}}^{*},{t}^{*})$ represents the dynamic process controlled by the transport and accumulation of chemical reactants. Based on the critical values or ranges, we formulate the indicators as follows:
Typically, because of the complexity of the processes, no analytical solutions are available for formulating the indicators. However, Monte Carlo simulations can be useful in constructing such indicators. For this case, an HSHM at any given location and tim (${\mathbf{\Omega}}^{*},{t}^{*}$) will only be triggered when all of the conditions are met and the ensemble mean of the indicator assumes the following form:
where ${I}_{\mathrm{d},i}({\mathbf{\Omega}}^{*},{t}^{*})$ is the value that the indicator assumes in the ith realization and N is the total number of realizations.
Overall, our choices of the three studies should not limit the generalizability of the indicator statistics approach for deriving statistical formulations for HSHM applications. The critical conditions chosen to construct the indicators are determined solely on the findings from these selected studies, and they will vary under different scenarios.
This section focuses on HSHMs in the subsurface for demonstration of linking HSHM models with the contributing physical processes, such as the migration of groundwater carrying reducing substrates, nuclear waste transport within the subsurface, the accumulation and transport of dense nonaqueous phase liquid (DNAPL) and other biogeochemical processes. Some current modeling approaches that focus on subsurface HSHMs assume simplified hydrologic structures (e.g., homogeneous and isotropic domains) in quantifying the fate and transport of solutes in the subsurface. However, such assumptions neglect the effect of the heterogeneity in the subsurface, potentially missing localized HSHMs arising as the combined effect of heterogeneity in physical and geochemical properties, and do not allow us to assess uncertainties in the HSHM occurrences.
This section therefore focuses on HSHMs taking place in the subsurface, with a particular emphasis on the role of spatial variability of the hydrologic parameters. Section 4.1 illustrates the potential of subsurface heterogeneity for triggering and timing of HSHMs. In Sect. 4.2, we develop closedform analytical solutions for HSHM probability. In doing this, we demonstrate the linking between our indicator model and the physics of the HSHMs in the subsurface. In Sect. 4.3 we demonstrate applications under various conditions of spatial variability.
4.1 Importance of spatial variability in the subsurface
The heterogeneous structure of hydraulic conductivity leads to significant variability in the transport of solutes in the subsurface, which couples with heterogeneous geochemical properties leading to a spatially varying reactivity (Arora et al., 2019b; Loschko et al., 2016; Sassen et al., 2012; Wainwright et al., 2015). Figure 2 demonstrates the uncertainty associated with HSHMs by looking at the flow fields in twodimensional loghydraulic conductivity (Y=ln (K)) fields with streamlines resulting from a uniform mean head gradient, left to right. The three panels differ in terms of the variance, ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$, of the log conductivity. The covariance function used for generating the fields is exponential and isotropic. ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$ is shown to have a profound impact upon the conductivity field. As the variance increases, regions of high and low log conductivity emerge, creating preferential flow paths bypassing the lowconductivity zones as shown by particle trajectories. At smaller variance (i.e., ${\mathit{\sigma}}_{Y}^{\mathrm{2}}=\mathrm{0.1}$), particles mainly travel along the mean flow direction with very limited departure from the mean trajectory, which are the straight lines connecting the left and right boundaries. In this situation, the arrival times of solute particles at a critical location (i.e., Ω^{*}, where for example geochemical conditions are favorable for certain types of reactions to occur) are predictable. With large variances (i.e., ${\mathit{\sigma}}_{Y}^{\mathrm{2}}=\mathrm{2}$), the streamlines assume a very irregular, hardtopredict geometry, and we can observe the emergence of flow channels, where particles can move quickly next to stagnant flow regions. Arrival times become more uncertain, because the exact geometry of the streamlines is hard to predict unless the Y field is known deterministically. However, since this is never the case, in another equally likely realization of the Y field, the situation may be different, resulting in significant uncertainties in predicting the particle travel times. Thus, spatial variability of log conductivity is a major uncertaintyinducing factor and, by extension, requires the need for stochastic modeling of HSHMs in situations where the associated processes and attributes are subject to uncertainty. In the following sections, we will present illustrative examples to analyze how subsurface spatial variability influences $<{I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},{t}^{*})>$, including variance and anisotropy ratio of the log conductivity.
4.2 Illustrative example and indicator formulation
In this section we illustrate the proposed indicator approach by means of synthetic case studies developed by using methods of stochastic hydrogeology. The choice of the synthetic case studies does not limit our approaches to broader applications where stochastic modeling with Monte Carlo simulations is applicable. Figure 3 displays the configuration of this case example. Consider the case of an instantaneous point source release of a target compound at location x_{0} and time t_{0}. HSHMs are triggered at any (${\mathbf{\Omega}}^{*},{t}^{*}$) if the solute is present. Consider the hot spot (Ω^{*}) to be confined within the following volume: ${w}_{\mathrm{1}}\le {x}_{\mathrm{1}}\le {w}_{\mathrm{1}}^{\prime}$; ${w}_{\mathrm{2}}\le {x}_{\mathrm{2}}\le {w}_{\mathrm{2}}^{\prime}$; ${w}_{\mathrm{3}}\le {x}_{\mathrm{3}}\le {w}_{\mathrm{3}}^{\prime}$; the dynamic indicator is therefore defined as follows:
The injected solute can be modeled in a Lagrangian framework as a particle moving according to the velocity field without changing its volume. The latter is the consequence of neglecting porescale dispersion. The expected value of this dynamic indicator at t^{*} is therefore
where ${f}_{\mathbf{X}\left({\mathit{t}}^{*}t\right)}(\mathit{a}\mathrm{}{\mathit{x}}_{\mathrm{0}},{t}_{\mathrm{0}})$ is the probability distribution function (pdf) of the particle's trajectory at time t^{*} (Dagan and Nguyen, 1989; Rubin, 2003). If we also assume steady, uniform in the average flow with mild heterogeneity of the log hydraulic conductivity field with the Gaussian displacement pdf, then we can compute $<{I}_{\mathrm{HSHM}}\left({\mathbf{\Omega}}^{*}{t}^{*}\right)>$ analytically using the following equation:
which can be integrated to yield
For simplicity, but without lack of generality, in Eq. (20) we assumed ${\mathit{x}}_{\mathrm{0}}=(\mathrm{0},\mathrm{0},\mathrm{0})$. The displacement variances X_{ii}, i=1, 2, 3 depend on the spatial distribution of the hydraulic conductivity in the subsurface. Equations (A4) to (A6) present the displacement variances for an axisymmetric exponential covariance function of the log conductivity (Eq. A3) given in Appendix A. In Eqs. (A7) to (A19), we have provided derivations of indicator formulations for other HSHM scenarios, including indicator formulation for complex concentration thresholds and indicator formulation for hot moment durations. Notice that in obtaining Eq. (20) we postulated ergodicity, which in practical terms reflects the actual situation of an instantaneous injection into a source zone with a transverse dimension much larger than the integral scale of the hydraulic conductivity (Dagan, 1990), such that the ensemble mean is representative of the effects of the actual, but unknown, distribution of hydraulic conductivity.
4.3 Probability of HSHM occurrence controlled by subsurface heterogeneity
In the following sections, we present the results from the case study described in Sect. 4.2. Specifically, in Sect. 4.3.1 and 4.3.2, we explore how heterogeneity of loghydraulic conductivity influences the probability of HSHM occurrences. To make results as general as possible, lengths are made dimensionless with respect to the integral scales (${I}_{{Y}_{\mathrm{h}}}$ in the two horizontal directions and ${I}_{{Y}_{\mathrm{v}}}$ in the vertical one) and time with respect to the following advective timescale: ${I}_{{Y}_{\mathrm{h}}}/U$, where U is the mean velocity. In the following, we explore the effect of the remaining parameters, i.e., the anisotropy ratio $e=\frac{{I}_{{Y}_{\mathrm{V}}}}{{I}_{{Y}_{\mathrm{H}}}}$ and the variance of the log conductivity ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$, on the emergence of an HSHM. We placed Ω^{*} along the mean trajectory at ($\mathrm{21}{I}_{{Y}_{\mathrm{H}}}$, 0, 0) with dimensions as ($\mathrm{2}{I}_{{Y}_{\mathrm{H}}}$, $\mathrm{2}{I}_{{Y}_{\mathrm{H}}}$, $\mathrm{2}{I}_{{Y}_{\mathrm{V}}}$). The dimensions of the hot spot are therefore of two integral scales in the three coordinate directions (x_{1}, x_{2}, x_{3}) and are placed at a dimensionless distance of 21 from the point source.
4.3.1 Dependence of $P\left({I}_{\mathrm{HSHM}}\right({\mathrm{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$ on variance in the spatial correlation structure of the log conductivity
Figure 4 shows $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$, which is the probability that an HSHM is triggered at Ω^{*} and the dimensionless time $\mathit{\tau}=\frac{tU}{{I}_{{Y}_{\mathrm{h}}}}$, for a few values of ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. Here, ${I}_{\mathrm{HSHM}}({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1}$ represents the situation in which an HSHM is triggered. At early time (e.g., τ<5), larger probability $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$ is observed with increase in ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. At intermediate time, i.e., at times comparable with the mean travel time τ=21, $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$ is inversely proportional to ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. At late time (e.g., τ>40), the largest $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$ occurs at intermediate ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. We observe that ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$ regulates the timing of the peak in $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$, which is located in the proximity of the mean travel time, τ=21, for weak heterogeneity and shifts towards earlier times as ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$ increases. From the practical perspective, Fig. 4 shows the probability of developing an HSHM at the identified position Ω^{*} at the given time τ.
These effects relate to the relationship between travel times (from the source to Ω^{*}) and ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. The key point to note is that ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$ controls the spread of the travel time around its mean value. A larger ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$ enhances channeling effects (Fiori and Jankovic, 2012; Moreno and Tsang, 1994, also in Fig. 2), which in turn enable earlier arrival times. However, at the same time, it also leads to the emergence of lowconductivity zones with low velocity or stagnant groundwater. The solute tends to bypass low hydraulic conductivity zones, as shown by the streamlines depicted in Fig. 2; however, the small amount of solute that actually penetrates these zones by slow advection and diffusion gets trapped for long time before being released, and this results in an extended tailing with low concentration, which increases the probability of observing an HSHM at later times. Thus, with an increase in ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$, we notice an increase in the probability of observing both increasingly earlier and increasingly delayed arrival times, which widens the probability distribution. By contrast at small variance, particles deviate little from the ensemble mean trajectory, because of the small contrast in conductivity between high and lowconductivity zones. This results in small particle spreading and travel times that differ only slightly from the mean travel time (τ=21) and a probability distribution less spread around the mean, where the peak is observed.
In summary, hydraulic conductivity contrast between low and highconductivity lithofacies increases with ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$, leading to the emergence of organized highconductivity pathways sneaking through surrounding lowconductivity zones, with the latter acting as “trapping” elements. This causes the emergence of both early and late arrival times, with the consequent larger probability of triggering HSHMs at early and later times, with respect to the case of low heterogeneity. Early arrival times are controlled by the connected highconductivity pathways, and the late arrival times are influenced by the lowconductivity zones, which act as lowrelease reservoirs for solutes.
4.3.2 Dependence of $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$ on anisotropy in the spatial correlation structure of the loghydraulic conductivity
The discussion here (accompanying Fig. 5) focuses on the impact of the anisotropy ratio in the correlation structure (e, defined above) on the probability of triggering HSHMs. The anisotropy ratio, e, provides an indication of the persistence of the log conductivity (Y) in the principal directions. The spatial correlation model used here for demonstration is that of axis symmetry, which is common to sedimentary formations (Dagan, 1989; Rubin, 2003), with e providing the ratio between the persistence of Y in the vertical (x_{3}) direction, represented by ${I}_{{Y}_{\mathrm{V}}}$, and the ones on the horizontal plane (x_{1}–x_{2}), represented by ${I}_{{Y}_{\mathrm{H}}}$. In unconsolidated sedimentary formations, ${I}_{{Y}_{\mathrm{V}}}$ is typically smaller than ${I}_{{Y}_{\mathrm{H}}}$ by as much as 1 order of magnitude, due to the different timescales of the depositional process taking place in the horizontal and vertical directions, which leads to thin and elongated lithofacies and consequently to a much smaller persistence of Y values in the normal to horizontal plane (Miall, 1985, 1988; Ritzi et al., 2004).
Figure 5 compares $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$ between formations defined by different anisotropy ratios and different ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. It shows that we have two factors to consider when explaining the differences in $P\left({I}_{\mathrm{HSHM}}\right({\mathbf{\Omega}}^{*},\mathit{\tau})=\mathrm{1})$. The first factor, as discussed earlier, is the widening of the probability distribution (direct consequence of the widening of the travel time distribution) due to increase in ${\mathit{\sigma}}_{Y}^{\mathrm{2}}$. With larger variance, we observe higher probabilities of departure of the travel times away from the average. The anisotropy ratio e adds a compounding factor. To understand its effect, we should recall the analyses of lateral displacement variances of solute particles moving in heterogeneous formations (cf. Dagan, 1989, and Eqs. A4 to A6 here), showing that smaller e leads to smaller lateral (both vertical and horizontal) displacement variances, implying smaller probabilities of lateral departures from the mean flow trajectory. Smaller e limits lateral spreads and increases the probability of particles entering Ω^{*}, sooner or later, and triggering HSHMs. The effect could also be viewed as a channeling effect of sorts: smaller e implies Y blocks of small aspect ratio (i.e., long but thin elements), which provide fast tracks for particles when defined by high Y values while blocking lateral spreads when defined by low Y values.
Additional notes: first, Ω^{*} known in the present analysis is located downstream from the source, along with the mean trajectory of the solute displacement. We expect different results in situations where Ω^{*} is positioned at an offset with respect to the mean flow direction or when its position is unknown. In both cases we expect a reduction of the probability of triggering an HSHM. Relevant to our discussion is that the proposed probabilistic framework can address the case of an unknown position for Ω^{*} as well. Second, we note that the analytical models used to compute the displacement statistics are formally limited to small variance of the log conductivity (${\mathit{\sigma}}_{Y}^{\mathrm{2}}<\mathrm{1}$), although they are shown to provide good approximations for large variances (Bellin et al., 1992; Salandin and Fiorotto, 1998). Third, the stochastic formulation provides the theoretical and computational formalism for conditioning the probabilities on in situ measurements (Ezzedine and Rubin, 1996; Rubin and Dagan, 1992) as well as on information borrowed from similar sites (Li et al., 2018; Cucchi et al., 2019).
In this study, we developed a general stochastic framework for characterizing the spatiotemporal distribution of environmental hot spots and hot moments (HSHMs). The stochastic formulation is built around the following principles.

The HSHMs are defined as random variables, and the goal is to derive their stochastic distribution in terms of the relevant processes and attributes.

The processes and attributes are modeled as stochastic processes and random variables, respectively, based on the underlying physics.

The static contributors are modeled stochastically using geostatistical space random functions.

The dynamic contributors are modeled stochastically using probability distribution functions derived from the underlying mathematical–physical models.

Several HSHM categories are defined, based on the contributing factors, as follows: HSHMs defined by dynamic contributors only, HSHMs defined by static contributors, and, most commonly, HSHMs requiring the coupling of static and dynamic contributors. The HSHM stochastic formulations are expressed in terms of the stochastic formulations of the relevant contributors.
We provided a detailed review of multiple HSHMs and showed how they relate to our definitions.
The framework we proposed in this study is advantageous in that it allows us to calculate the uncertainty associated with HSHMs based on the uncertainty associated with its contributors. Additionally, it provides a formalism, well established by Bayesian theory, for conditioning the HSHM probabilities on in situ measurements as well as on information borrowed from geologically and otherwise similar sites.
We demonstrated our proposed approach through applications in the area of subsurface transport and hydrogeology, focusing on the impacts of subsurface heterogeneity on HSHMs. We analyzed, quantitatively, how subsurface heterogeneity of the conductivity field controls the HSHM statistics, for example, the time expected for the probability of the HSHM to occur to reach a priori set thresholds or time to peak probability.
Lastly, as mentioned both here and in previous studies, statistical methods for quantifying the occurrences of HSHMs and the associated uncertainties are needed to advance our understanding of the mechanisms that cause HSHMs, as well as to enhance our ability to predict HSHMs and manage their consequences.
A1 Equations for the displacement pdf
Assuming steady uniformity in the average flow, with mild heterogeneity of the loghydraulic conductivity field with Gaussian displacement, the pdf of the longitudinal (x_{1}, Fig. 3) displacement of a solute particle starting at time t_{0}=0 at ${\mathit{x}}_{\mathrm{0}}=(\mathrm{0},\mathrm{0},\mathrm{0})$ is given by the following equation (Dagan and Nguyen, 1989; Dagan and Rubin, 1992):
Additionally, the displacement pdf in the transverse directions (x_{2} and x_{3}) is given by
A2 Equations for displacement variances under anisotropic conditions
Dagan (1984) developed a solution of the displacement variances X_{ii}, i=1, 2, 3 for an exponential and axisymmetric logconductivity covariance:
where r is the twopoint separation distance and 〈Y〉 the ensemble mean of the log conductivity Y=ln K. $\stackrel{\mathrm{\u203e}}{{J}_{\mathrm{0}}}$ and $\stackrel{\mathrm{\u203e}}{{J}_{\mathrm{1}}}$ are, respectively, the zero and first orders of the firstkind Bessel functions.
A3 Equations for displacement variances under isotropic conditions
Dagan (1984) provided analytical solutions for longitudinal and transverse displacement variances. This is a special case for the anisotropic case with e=1. The solutions are as follows:
A4 Indicator formulation for complex concentration thresholds
When considering local dispersion or in case of a reactive tracer, the condition that the particle is inside the volume Ω^{*} does not suffice to define the dynamic indicator, and a concentration threshold C_{th} should be introduced:
In the absence of local dispersion and for a reactive solute decaying at a (spatially) constant rate k, the ensemble mean assumes the following expression (Cvetkovic and Shapiro, 1990):
where C_{0} is the initial concentration and H[⋅] is the Heaviside step function. The ensemble mean (Eq. A10) is the product of the probability that the particle assumes a concentration larger than the threshold at t^{*} (given that reaction rate k is constant, this probability is either 0 or 1) and the probability that at the same time t^{*} the particle is within the hot spot Ω^{*}. In other words, Eq. (A10) expresses the fact that a particle inside Ω^{*} contributes to the hot moment only if its concentration is greater than the threshold, and this occurs for ${t}^{*}<\frac{\mathrm{1}}{k}\mathrm{ln}\left(\frac{{C}_{\mathrm{0}}}{{C}_{\mathrm{th}}}\right)$. Equation (A10) can be generalized to the cases of instantaneous injection into a source of volume V_{0}, as discussed before for the nonreactive case. For other complex situations, such as that in which k is spatially variable and complex reaction networks, the ensemble mean of the indicators can be addressed by Eq. (16) in a Monte Carlo framework.
A5 Indicator formulation for hot moment durations
As hot moments can persist over short time periods, estimating the corresponding probabilities for any given time interval becomes also very important. The probability that the hot moment persists over the interval [t_{1}, t_{2}] at Ω^{*} can be formally computed as follows:
where $P({t}_{\mathrm{1}},{\mathbf{\Omega}}^{*})$ is the probability that the particle is inside Ω^{*} at time ${t}^{*}={t}_{\mathrm{1}}$ and $P({t}_{\mathrm{2}}\mathrm{}{t}_{\mathrm{1}},{\mathbf{\Omega}}^{*})$ is the probability that the particle is still inside Ω^{*} at time ${t}^{*}={t}_{\mathrm{2}}$, provided that at time t_{1}, it was also inside Ω^{*}. If the particle exits Ω^{*} during interval [t_{1}, t_{2}], this time interval will not be qualified as a hot moment, and thus the probability computation needs to ensure the particle stays within Ω^{*} during the entire time interval.
Under the firstorder approximation (FOA) (see, e.g., Dagan, 1989; Gelhar, 1993; Rubin, 2003), the pdf of the particle displacement is normal with mean $<\mathbf{X}({t}^{*};{\mathit{x}}_{\mathrm{0}},{t}_{\mathrm{0}})>$ and autocovariance tensor of the residual displacements ${\mathbf{X}}^{\prime}\left({t}^{*}\right)=\mathbf{X}\left({t}^{*}\right)\langle \mathbf{X}\left({t}^{*}\right)\rangle $ defined by ${\mathbf{X}}_{ij}({t}^{*};{\mathit{x}}_{\mathrm{0}},{t}_{\mathrm{0}})=\langle {\mathbf{X}}_{i}^{\prime}({t}^{*};{\mathit{x}}_{\mathrm{0}},{t}_{\mathrm{0}}){\mathbf{X}}_{j}^{\prime}({t}^{*};{\mathit{x}}_{\mathrm{0}},{t}_{\mathrm{0}})\rangle $, $i,j=\mathrm{1}$, 2, 3. For simplicity in the following, we assume x_{0}=0 and t_{0}=0. Under these assumptions,
where the conditional pdf ${f}_{\mathbf{X}\left({t}_{\mathrm{2}}\right)}^{\mathrm{c}}\left(\mathit{b}\mathrm{}\mathbf{X}\right({t}_{\mathrm{1}})=\mathit{a})$ is multinormally distributed with conditional mean and variance tensor given by
and
respectively, which further yields the following:
where $\cdot $ indicates the determinant, exp is the exponential function and the exponent T indicates the transpose of the vector.
In Eqs. (A13) and (A14), ${\mathbf{X}}^{\prime}\left({t}^{*}\right)=\mathbf{X}=\left({t}^{*}\right)\langle \mathbf{X}\left({t}^{*}\right)\rangle $ stands for the departure of the particle's displacement with respect to the ensemble mean trajectory, and Var[X]^{−1} is the autocovariance tensor of the residual displacement whose elements are defined above. Similarly, Cov[X^{′}(t_{1}), X^{′}(t_{2})] is the covariance tensor of residual displacement whose elements are ${\mathbf{X}}_{ij}({t}_{\mathrm{1}},{t}_{\mathrm{2}};{\mathit{x}}_{\mathbf{0}},{t}_{\mathrm{0}})=\langle {X}_{i}^{\prime}\left({t}_{\mathrm{1}}\right){X}_{j}^{\prime}\left({t}_{\mathrm{2}}\right)\rangle $, $i,j=\mathrm{1}$, 2, 3. Note that in the general threedimensional case 〈X(t_{2})X(t_{1})=a〉 is a threedimensional vector and σ(t_{1},t_{2}) is a 3×3 secondorder tensor.
For t_{2}→t_{1}, ${f}_{X\left({t}_{\mathrm{2}}\right)}^{\mathrm{c}}\left[\mathit{b}\mathrm{}\mathbf{X}\right({t}_{\mathrm{1}})=\mathit{a}]\to \mathit{\delta}(\mathit{b}\mathit{a})$, where δ(⋅) is the Dirac delta, such that $P({t}_{\mathrm{2}}\mathrm{}{t}_{\mathrm{1}},{\mathbf{\Omega}}^{*})\to \mathrm{1}$. On the other hand, for t_{2}≫t_{1}, Cov[X^{′}(t_{1}), ${\mathbf{X}}^{\prime}\left({t}_{\mathrm{2}}\right)]\to \mathrm{0}$ and $P({t}_{\mathrm{2}}\mathrm{}{t}_{\mathrm{1}},{\mathbf{\Omega}}^{*})\to P({t}_{\mathrm{2}},{\mathbf{\Omega}}^{*})$, the marginal probability that the particle is within Ω^{*} at time ${t}^{*}={t}_{\mathrm{2}}$. Equations (A12) to (A15) are obtained under the FOA and assume that the particle can enter Ω^{*} only once. Such an assumption is needed to obtain analytical solutions and is reasonable for situations with small to mild subsurface heterogeneity (e.g., ${\mathit{\sigma}}_{Y}^{\mathrm{2}}\le \mathrm{1.6}$), such as the cases presented in Bellin et al. (1992, 1994) and Cvetkovic et al. (1992). In particular, FOA assumes small heterogeneity, and under this assumption the particle trajectory deviates slightly from its ensemble mean, which is directed along the regional hydraulic head gradient. For a regular volume Ω^{*}, this reduces the probability of the particle entering more than once the hot spot. This probability reduces further if in the horizontal and vertical transverse directions Ω^{*} is much larger than the respective integral scales, because the probability of observing negative longitudinal velocity components (i.e., along the mean flow field) is much smaller than in the transverse directions (Bellin et al., 1992) and vanishes as formation heterogeneity reduces.
If the hotspot Ω^{*} is the volume confined between two planes at ${x}_{\mathrm{1}}\frac{{l}_{\mathrm{1}}}{\mathrm{2}}$ and ${x}_{\mathrm{1}}+\frac{{l}_{\mathrm{1}}}{\mathrm{2}}$, with the other two dimensions much larger than the transverse horizontal and vertical integral scales, l_{2}≫I_{h}, l_{3}≫I_{v}, Eq. (A13), simplify to
where X_{1} is the longitudinal component of the particle's trajectory and ${f}_{{X}_{\mathrm{1}}\left({t}^{*}\right)}^{\mathrm{c}}$ is its conditional pdf, which is normal with conditional mean and variance given by
and
respectively. Consequently, ${f}_{{X}^{\mathrm{c}}\left({t}^{*}\right)}$ in Eq. (A16) assumes the following form:
Substituting Eq. (A16) into Eq. (A19) allows us to compute $<{I}_{\mathrm{d}}({\mathbf{\Omega}}^{*}$, t_{1}, ${t}_{\mathrm{2}})>.$ For situations where the FOA assumptions are not valid (e.g., large heterogeneity), a Monte Carlo simulation framework is still applicable as an alternative approach to construct the dynamics indicators (see Eq. 16).
No data sets were used in this article.
JC performed the analysis and prepared the initial manuscript. All the coauthors contributed to the discussion of the results and revision of the paper.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank the editor Nunzio Romano and the anonymous reviewers for their comments and advice.
This research has been supported by the Jane Lewis Fellowship from University of California, Berkeley, Office of Science, Office of Advanced Scientific Computing (“Deduce: Distributed Dynamics Data Analytics Infrastructure for Collaborative Environments” under contract no. DEAC0205CH11231) and the Italian Ministry of Education, University and Research (MIUR) in the framework of the Departments of Excellence Initiative 20182022 granted to DICAM of the University of Trento.
This paper was edited by Nunzio Romano and reviewed by James Stegen and two anonymous referees.
Abbott, B. W., Baranov, V., MendozaLera, C., Nikolakopoulou, M., Harjung, A., Kolbe, T., Balasubramanian, M. N., Vaessen, T. N., Ciocca, F., Campeau, A., Wallin, M. B., Romeijn, P., Antonelli, M., Gonçalves, J., Datry, T., Laverman, A. M., de Dreuzy, J. R., Hannah, D. M., Krause, S., Oldham, C., and Pinay, G.: Using multitracer inference to move beyond singlecatchment ecohydrology, EarthSci. Rev., 160, 19–42, https://doi.org/10.1016/j.earscirev.2016.06.014, 2016.
Akaike, H.: A New Look at the Statistical Model Identification, IEEE Trans. Automat. Contr., 19, 716–723, https://doi.org/10.1109/TAC.1974.1100705, 1974.
Andrews, D. M., Lin, H., Zhu, Q., Jin, L., and Brantley, S. L.: Hot spots and hot moments of dissolved organic carbon export and soil organic carbon storage in the shale hills catchment, Vadose Zone J., 10, 943–954, https://doi.org/10.2136/vzj2010.0149, 2011.
Arora, B., Mohanty, B. P., McGuire, J. T., and Cozzarelli, I. M.: Temporal dynamics of biogeochemical processes at the Norman Landfill site, Water Resour. Res., 49, 6909–6926, https://doi.org/10.1002/wrcr.20484, 2013.
Arora, B., Dwivedi, D., Hubbard, S. S., Steefel, C. I., and Williams, K. H.: Identifying geochemical hot moments and their controls on a contaminated river floodplain system using wavelet and entropy approaches, Environ. Model. Softw., 85, 27–41, https://doi.org/10.1016/j.envsoft.2016.08.005, 2016a.
Arora, B., Spycher, N. F., Steefel, C. I., Molins, S., Bill, M., Conrad, M. E., Dong, W., Faybishenko, B., Tokunaga, T. K., Wan, J., Williams, K. H., and Yabusaki, S. B.: Influence of hydrological, biogeochemical and temperature transients on subsurface carbon fluxes in a flood plain environment, Biogeochemistry, 127, 367–396, https://doi.org/10.1007/s1053301601868, 2016b.
Arora, B., Wainwright, H. M., Dwivedi, D., Vaughn, L. J. S., Curtis, J. B., Torn, M. S., Dafflon, B., and Hubbard, S. S.: Evaluating temporal controls on greenhouse gas (GHG) fluxes in an Arctic tundra environment: An entropybased approach, Sci. Total Environ., 649 284–299, https://doi.org/10.1016/j.scitotenv.2018.08.251, 2019a.
Arora, B., Dwivedi, D., Faybishenko, B., Jana, R. B., and Wainwright, H. M.: Understanding and Predicting Vadose Zone Processes, Rev. Mineral. Geochem., 85, 303–328, https://doi.org/10.2138/rmg.2019.85.10, 2019b.
Arora, B., Briggs, M. A., Zarnetske, J., Stegen, J., GomezValez, J., Dwivedi, D., and Steefel, C. I.: Hot Spots and Hot Moments in the Critical Zone: Identification of and Incorporation into Reactive Transport Models, In: Biogeochemistry of the Critical Zone, edited by: Wymore, A., Yang, W., Silver, W., McDowell, B., and Chorover, J., SpringerNature, in press, 2021.
Bellin, A. and Rubin, Y.: On the use of peak concentration arrival times for the inference of hydrogeological parameters, Water Resour. Res., 40, W07401, https://doi.org/10.1029/2003WR002179, 2004.
Bellin, A., Salandin, P., and Rinaldo, A.: Simulation of dispersion in heterogeneous porous formations: Statistics, firstorder theories, convergence of computations, Water Resour. Res., 28, 2211–2227, https://doi.org/10.1029/92WR00578, 1992.
Bellin, A., Rubin, Y., and Rinaldo, A.: EulerianLagrangian approach for modeling of flow and transport in heterogeneous geological formations, Water Resour. Res., 30, 2913–2924, https://doi.org/10.1029/94WR01489, 1994.
Bernhardt, E. S., Blaszczak, J. R., Ficken, C. D., Fork, M. L., Kaiser, K. E., and Seybold, E. C.: Control Points in Ecosystems: Moving Beyond the Hot Spot Hot Moment Concept, Ecosystems, 20, 665–682, https://doi.org/10.1007/s100210160103y, 2017.
Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology, Hydrol. Sci. Bull., 24, 43–69, https://doi.org/10.1080/02626667909491834, 1979.
Bundt, M., Widmer, F., Pesaro, M., Zeyer, J., and Blaser, P.: Preferential flow paths: Biological “hot spots” in soils, Soil Biol. Biochem., 33, 729–738, https://doi.org/10.1016/S00380717(00)002182, 2001.
Cucchi, K., Heße, F., Kawa, N., Wang, C., and Rubin, Y.: Exsitu priors: A Bayesian hierarchical framework for defining informative prior distributions in hydrogeology, Adv. Water Resour., 126, 65–78, https://doi.org/10.1016/j.advwatres.2019.02.003, 2019.
Cvetkovic, V., Shapiro, A. M., and Dagan, G.: A solute flux approach to transport in heterogeneous formations: 2. Uncertainty analysis, Water Resour. Res., 28, 1377–1388, https://doi.org/10.1029/91WR03085, 1992.
Cvetkovic, V. D. and Shapiro, A. M.: Mass arrival of sorptive solute in heterogeneous porous media, Water Resour. Res., 26, 2057–2067, https://doi.org/10.1029/WR026i009p02057, 1990.
Dagan, G.: Solute transport in heterogeneous porous formations, J. Fluid Mech., 145, 151–177, https://doi.org/10.1017/S0022112084002858, 1984.
Dagan, G.: Flow and Transport in Porous Formations, Springer Verlag, Berlin, 1989.
Dagan, G.: Transport in heterogeneous porous formations: Spatial moments, ergodicity, and effective dispersion, Water Resour. Res., 26, 1281–1290, https://doi.org/10.1029/WR026i006p01281, 1990.
Dagan, G. and Nguyen, V.: A comparison of travel time and concentration approaches to modeling transport by groundwater, J. Contam. Hydrol., 4, 79–91, https://doi.org/10.1016/01697722(89)900272, 1989.
Dagan, G. and Rubin, Y.: Conditional estimation of solute travel time in heterogeneous formations: Impact of transmissivity measurements, Water Resour. Res., 28, 1033–1040,, 1992.
Destouni, G. and Cvetkovic, V.: Field scale mass arrival of sorptive solute into the groundwater, Water Resour. Res., 27, 1315–1325, https://doi.org/10.1029/91WR00182, 1991.
Duncan, J. M., Groffman, P. M., and Band, L. E.: Towards closing the watershed nitrogen budget: Spatial and temporal scaling of denitrification, J. Geophys. Res.Biogeo., 118, 1105–1119, https://doi.org/10.1002/jgrg.20090, 2013,
Dwivedi, D.: Hot Spots and Hot Moments of Nitrogen in a Riparian Corridor, Water Resour. Res., 54, 205–222, https://doi.org/10.1002/2017WR022346, 2017.
Ezzedine, S. and Rubin, Y.: A geostatistical approach to the conditional estimation of spatially distributed solute concentration and notes on the use of tracer data in the inverse problem, Water Resour. Res., 32, 853–861, https://doi.org/10.1029/95WR02285, 1996.
Fiori, A. and Jankovic, I.: On Preferential Flow, Channeling and Connectivity in Heterogeneous Porous Formations, Math. Geosci., 44, 133–145, https://doi.org/10.1007/s1100401193652, 2012.
Frei, S., Knorr, K. H., Peiffer, S., and Fleckenstein, J. H.: Surface microtopography causes hot spots of biogeochemical activity in wetland systems: A virtual modeling experiment, J. Geophys. Res.Biogeo., 117, 1–18, https://doi.org/10.1029/2012JG002012, 2012.
Gelhar, L. W.: Stochastic Subsurface Hydrogeology, PrenticeHall, Upper Saddle River, N.J., 1993.
Gu, C., Anderson, W., and Maggi, F.: Riparian biogeochemical hot moments induced by stream fluctuations, Water Resour. Res., 48, 1–17, https://doi.org/10.1029/2011WR011720, 2012.
Harms, T. K. and Grimm, N. B.: Hot spots and hot moments of carbon and nitrogen dynamics in a semiarid riparian zone, J. Geophys. Res.Biogeo., 113, 1–14, https://doi.org/10.1029/2007JG000588, 2008.
Henri, C. V., FernàndezGarcia, D., and de Barros, F. P. J.: Probabilistic human health risk assessment of degradationrelated chemical mixtures in heterogeneous aquifers: Risk statistics, hot spots, and preferential channels, Water Resour. Res., 51, 4086–4108, https://doi.org/10.1002/2014WR016717, 2015.
Hill, A. R., Devito, K. J., and Campagnolo, S.: Subsurface Denitrification in a Forest Riparian Zone: Interactions between Hydrology and Supplies of Nitrate and Organic Carbon Author(s): Alan R. Hill, Kevin J. Devito, S. Campagnolo and K. Sanmugadas, Springer Stable, 51, 193–223, 2000.
Hubbard, S. S., Williams, K. H., Agarwal, D., Banfield, J., Beller, H., Bouskill, N., Brodie, E., Carroll, R., Dafflon, B., Dwivedi, D., Falco, N., Faybishenko, B., Maxwell, R., Nico, P., Steefel, C., Steltzer, H., Tokunaga, T., Tran, P. A., Wainwright, H., and Varadharajan, C.: The East River, Colorado, Watershed: A Mountainous Community Testbed for Improving Predictive Understanding of Multiscale Hydrological–Biogeochemical Dynamics, Vadose Zone J., 17, 180061, https://doi.org/10.2136/vzj2018.03.0061, 2018.
Li, X., Li, Y., Chang, C. F., Tan, B., Chen, Z., Sege, J., Wang, C., and Rubin, Y.: Stochastic, goaloriented rapid impact modeling of uncertainty and environmental impacts in poorlysampled sites using exsitu priors, Adv. Water Resour., 111, 174–191, https://doi.org/10.1016/j.advwatres.2017.11.008, 2018.
Loschko, M., Woehling, T., Rudolph, D. L., and Cirpka, O. A.: Cumulative relative reactivity: A concept for modeling aquiferscale reactive transport, Water Resour. Res., 52, 8117–8137, https://doi.org/10.1002/2016WR019080, 2016.
McClain, M. E., Boyer, E. W., Dent, C. L., Gergel, S. E., Grimm, N. B., Groffman, P. M., Hart, S. C., Harvey, J. W., Johnston, C. A., Mayorga, E., McDowell, W. H., and Pinay, G.: Biogeochemical Hot Spots and Hot Moments at the Interface of Terrestrial and Aquatic Ecosystems, Ecosystems, 6, 301–312, https://doi.org/10.1007/s1002100301619, 2003.
Miall, A. D.: Architecturalelement analysis: A new method of facies analysis applied to fluvial deposits, EarthSci. Rev., 22, 261–308, https://doi.org/10.1016/00128252(85)900017, 1985.
Miall A. D.: Facies Architecture in Clastic Sedimentary Basins, in: New Perspectives in Basin Analysis, Frontiers in Sedimentary Geology, edited by: Kleinspehn, K. L. and Paola, C. Springer, New York, NY, https://doi.org/10.1007/9781461237884_4, 1988.
Mitchell, C. P. J., Branfireun, B. A., and Kolka, R. K.: Spatial characteristics of net methylmercury production hot spots in peatlands, Environ. Sci. Technol., 42, 1010–1016, https://doi.org/10.1021/es0704986, 2008.
Moreno, L. and Tsang, C.F.: Flow channeling in strongly heterogeneous porous media: A numerical study, Water Resour. Res., 30, 1421–1430, https://doi.org/10.1029/93WR02978, 1994.
Ritzi, R. W., Dai, Z., Dominic, D. F., and Rubin, Y. N.: Spatial correlation of permeability in crossstratified sediment with hierarchical architecture, Water Resour. Res., 40, W03513, https://doi.org/10.1029/2003WR002420, 2004.
Rubin, Y.: Applied Stochastic Hydrogeology, Oxford University Press, Oxford, UK, 2003.
Rubin, Y. and Dagan, G.: Conditional estimation of solute travel time in heterogeneous formations: Impact of transmissivity measurements, Water Resour. Res., 28, 1033–1040, https://doi.org/10.1029/91WR02759, 1992.
Rubin, Y. and Journel, A. G.: Simulation of nonGaussian space random functions for modeling transport in groundwater, Water Resour. Res., 27, 1711–1721, https://doi.org/10.1029/91WR00838, 1991.
Salandin, P. and Fiorotto, V.: Solute transport in highly heterogeneous aquifers, Water Resour. Res., 34, 949–961, https://doi.org/10.1029/98WR00219, 1998.
Sassen, D. S., Hubbard, S. S., Bea, S. A., Chen, J., Spycher, N., and Denham, M. E.: Reactive facies: An approach for parameterizing fieldscale reactive transport models using geophysical methods, Water Resour. Res., 48, W10526, https://doi.org/10.1029/2011WR011047, 2012.
Schwarz, G.: Estimating the Dimension of a Model, Ann. Stat., 6, 461–464, https://doi.org/10.1214/aos/1176344136, 1978.
Shrestha, N. K. and Wang, J.: Current and future hotspots and hotmoments of nitrous oxide emission in a cold climate river basin, Environ. Pollut., 239, 648–660, https://doi.org/10.1016/j.envpol.2018.04.068, 2018.
Sørensen, R., Zinko, U., and Seibert, J.: On the calculation of the topographic wetness index: evaluation of different methods based on field observations, Hydrol. Earth Syst. Sci., 10, 101–112, https://doi.org/10.5194/hess101012006, 2006.
United States Environmental Protection Agency (USEPA): Risk assessment guidance for superfund, in Process for Conducting Probabilistic Risk Assessment Part A, Office of Emergency and Remedial Response, U.S. Environmental Protection Agency, Washington, D.C.,Vol. 3, 385 pp., 2001.
Vidon, P., Allan, C., Burns, D., Duval, T. P., Gurwick, N., Inamdar, S., Lowrance, R., Okay, J., Scott, D., and Sebestyen, S.: Hot spots and hot moments in riparian zones: Potential for improved water quality management, J. Am. Water Resour. Assoc., 46, 278–298, https://doi.org/10.1111/j.17521688.2010.00420.x, 2010.
Wainwright, H. M., Orozco, A. F., Bucker, M., Dafflon, B., Chen, J., Hubbard, S. S., and Williams, K. H.: Hierarchical Bayesian method for mapping biogeochemical hot spots using induced polarization imaging, Water Resour. Res., 51, 9127–9140, https://doi.org/10.1002/2014WR016259, 2015.
Wilson, A. and Rubin, Y.: Characterization of aquifer heterogeneity using indicator variables for solute concentrations, Water Resour. Res., 38, 191–1912, https://doi.org/10.1029/2000wr000116, 2002.
 Abstract
 Introduction
 Methodology and statistical formulation of HSHMs
 Examples of the statistical formulation of HSHMs with case studies
 HSHM applications in groundwater hydrology
 Discussion and summary
 Appendix A: Equations and derivations
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Methodology and statistical formulation of HSHMs
 Examples of the statistical formulation of HSHMs with case studies
 HSHM applications in groundwater hydrology
 Discussion and summary
 Appendix A: Equations and derivations
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References