Statistical theory of probabilistic hazard maps: a probability distribution for the hazard boundary location

The study of volcanic flow hazards in a probabilistic framework centers around systematic experimental numerical modeling of the hazardous phenomenon and the subsequent generation and interpretation of a probabilistic hazard map (PHM). For a given volcanic flow (e.g., lava flow, lahar, pyroclastic flow, ash cloud), the PHM is typically interpreted as the point-wise probability of inundation by flow material.

In the current work, we present new methods for calculating spatial representations of the mean, standard deviation, median, and modal locations of the hazard's boundary as ensembles of many deterministic runs of a physical model. By formalizing its generation and properties, we show that a PHM may be used to construct these statistical measures of the hazard boundary which have been unrecognized in previous probabilistic hazard analyses. Our formalism shows that a typical PHM for a volcanic flow not only gives the point-wise inundation probability, but also represents a set of cumulative distribution functions for the location of the inundation boundary with a corresponding set of probability density functions. These distributions run over curves of steepest probability gradient ascent on the PHM. Consequently, 2-D space curves can be constructed on the map which represents the mean, median, and modal locations of the likely inundation boundary. These curves give well-defined answers to the question of the likely boundary location of the area impacted by the hazard. Additionally, methods of calculation for higher moments including the standard deviation are presented, which take the form of map regions surrounding the mean boundary location. These measures of central tendency and variance add significant value to spatial probabilistic hazard analyses, giving a new statistical description of the probability distributions underlying PHMs.

The theory presented here may be used to aid construction of improved hazard maps, which could prove useful for planning and emergency management purposes. This formalism also allows for application to simplified processes describable by analytic solutions. In that context, the connection between the PHM, its moments, and the underlying parameter variation is explicit, allowing for better source parameter estimation from natural data, yielding insights about natural controls on those parameters.

1 Introduction

The probabilistic study of volcanic hazards is part of an emerging paradigm in physical volcanology, one that openly admits and seeks to calculate and characterize the uncertainty inherent to physical modeling of hazardous volcanic processes (Bursik et al., 2012; Connor et al., 2015; Bevilacqua, 2016; Sandri et al., 2016). Critically, this approach is founded on the notion that physical parameters related to these processes can be represented by a probability distribution. The distribution is propagated by a physical model to yield a space of hazard distributions which put probabilistic constraints on quantities that would be fixed in a deterministic, scenario-based modeling framework. In practice this has typically been accomplished by accumulating the statistics of the processes through deterministic modeling over an ensemble of model inputs (Neri et al., 2015; Biass et al., 2016; Tierz et al., 2016; Bevilacqua et al., 2017; Gallant et al., 2018; Patra et al., 2018b). Polynomial chaos approximations and Gaussian model surrogates can provide fast algorithms to obtain these statistics (Bursik et al., 2012; Madankan et al., 2012; Stefanescu et al., 2012; Spiller et al., 2014; Bayarri et al., 2015; Tierz et al., 2018). Owing to the nature of the studied processes, the results are typically displayed in map form and contain a large amount of spatial information, the interpretation of which has become a source of extensive study (Calder et al., 2015; Thompson et al., 2015). Because the information in these maps is potentially of great interest to a wide community of stakeholders including governments, planners, emergency managers, local residents, and others spanning a wide range of familiarity with probabilistic reasoning, finding the best approach to condensing these data-rich maps is nontrivial. The raw probabilities are considered less useful to many end users despite their significant information content (Thompson et al., 2015). Often the information in these maps is presented qualitatively as high, medium, and low hazard or simply high and low hazard (Calder et al., 2015; Wagner et al., 2015). The reduction of a complex probability distribution to a relative categorization as in a high–low hazard scheme generally involves significant nuance in blending scientific information, decision-making, and administrative considerations. In the scientific component of the analysis, such a categorization could be based on a probability threshold, in many cases a certain percentile of the distribution, which corresponds to a probability contour on the map. While this is the most obvious method it is just one of many options for declaring this categorical boundary.

At present, a probabilistic hazard map (PHM) for a given hazard is easily constructed from numerical modeling; however, its most commonly used statistical product in hazard communication is just the set of probability contours. Consequently if a binary volcanic hazard map were to be constructed from a PHM showing the region of no likely impact versus the region of likely impact, the analyst would have to decide where to draw the line: possibly at the 5 % probability contour (a safe, conservative choice), the median contour (50 %) (a good estimate for the most likely impact–no-impact boundary), or some other choice of percentile. Often the choice of the contour level is based on agreements or official protocols made by civil protection authorities. Such potentially high-consequence decisions would benefit from additional statistical analyses which generate other tools for this task. For example, it may be useful to know the average geographical limits of a particular hazard in addition to the point-wise probability of hazard impact. For any sample data set or theoretical probability distribution, the mean and variance as well as one or more modal values are calculable in addition to the median and other percentiles. An important deficiency in the analysis of the PHM is that, previously, scientific estimates of the likely hazard boundary (a single curve on the map) have lacked the formal framework required for unique determination of statistical estimates such as the mean, variance, and higher moments of the distribution of hazard boundary locations.

For any given volcanic hazard analysis which maps probabilities of a given event, the probability values may represent a wide array of concepts in probability theory including cumulative probabilities, probability densities, conditional probabilities, and others. Here, the point-wise probability is interpreted as the relative frequency of hazard impact as predicted by repeated simulation using a deterministic physical model or ensemble of models equipped with an ensemble of inputs spanning the realistic range of underlying uncertain parameters. In the following theory and analysis, we only consider spatial probability distributions which map the probability of a given hazard impacting a given location where there is at least one point in space with a modeled probability of unity. Although this excludes the study of many possible PHMs, it covers a sufficiently wide range of models which are regular or regularized and span a realistically limited subspace of the general parameter space; that is, they depend continuously on the inputs or, in this case, maintain some level of spatial correlation across a realistically narrow, continuous range of inputs. In general, the methods detailed here are agnostic with respect to the type of uncertainty which is modeled. For any probabilistic geophysical model, there exist epistemic and aleatory uncertainties. The division between aleatory uncertainty and epistemic uncertainty corresponds to the distinction between intrinsic randomness of the system and the additional uncertainty that affects its representation, the latter originating from incomplete information about the underlying physical processes. Clearly PHMs constructed purely from aleatory uncertainties in the inputs have the most concrete meanings. However, PHMs constructed from a mix of aleatory and epistemic uncertainties can be considered as a single member of a time-dependent sequence of such assessments, each of which captures the types of uncertainty encoded in the generation of the PHM that is understood at that time. In general, as epistemic uncertainty decreases over time, the later assessments will become more accurate. As long as the uncertainties can be described as a probability space, these analyses hold up to current geophysical knowledge. Although varying the parameters in a single model gives the most concrete meaning, these methods could be applied to PHMs constructed from an ensemble of models provided that the ensemble shows regular behavior. Here, we focus on analyses of single hazards, although the boundary-finding methods could be justified to analyze PHMs constructed for multiple hazards.

Throughout this work, we realize these methods with examples using geophysical mass flow (GMF) models of inundation hazards. This is done primarily because the flow boundary is easily definable and conceptualized in these models, giving a consistent intuition to the complicated mathematical quantities discussed here. However, these concepts extend to many other types of hazard models where a particular threshold is of interest in defining the region impacted by a hazard in volcanology and other fields including concentration thresholds in volcanic clouds (Bursik et al., 2012), thickness thresholds in tephra fallout (Sandri et al., 2016), ballistic ejecta impact count thresholds (Biass et al., 2016), pyroclastic density current invasion (Bevilacqua et al., 2017), lava flow inundation (Gallant et al., 2018), tsunami run-up thresholds (Geist and Parsons, 2006; Grezio et al., 2017), ground-shaking thresholds in seismic hazards (Kvaerna et al., 1999), hurricane wind speed thresholds (Splitt et al., 2014), and many others. All of these examples can be cast in the same mathematical formulation that we detail below. A full evaluation of suitable and unsuitable types of probabilistic hazard assessments is beyond the scope of this work.

The work of Stefanescu et al. (2012) provided a rigorous framework for calculating the PHM from a Gaussian emulator of GMF simulations. Within that framework, the PHM is the point-wise mean probability of hazard impact, and the uncertainty bounds for that probability are well-defined. Although this method is robust, it has not been applied in many hazard studies. Furthermore, as with most calculations of PHMs, the interval 0 < Prob(Hazard Impact) < 1 occupies a very large region of space. As in the above hypothetical binary hazard-mapping scenario, the question remains – from the perspective of a probabilistic volcanic hazard analysis, where is the likely boundary of the impacted area?

The goal of the present work is to extend the probability theory of PHMs to construct statistically meaningful answers to this question. To this end, we give an explicit mathematical definition of the PHM, defining explicitly the connection between the cumulative exceedance probabilities and probability densities in a PHM and estimates of the likely hazard boundary. More specifically, we calculate spatial representations of the PHM's moments including the mean and variance of the hazard boundary location as well as its modality. We stipulate that the methods detailed here apply only to the scientific analysis component of any hazard-mapping project and must be considered as just one piece of the hazard-map-making process. This effort combines mathematical tools and concepts from probability theory, the study of dynamical systems, and differential geometry.

2 Hazard map distribution theory

Solutions to flow models which involve depth-integrated partial differential equations (e.g., Saint-Venant-type shallow water equations) typically include solving for the distribution of flow thicknesses over position and parameter space:

$$\begin{array}{}\text{(1)}& h=h(\mathit{x};\mathit{\beta }),\end{array}$$

where x is the position vector and β∈ℬ is an n vector in parameter space ℬ; that is, if there are n input parameters to a model, ℬ⊂ℝⁿ. The parameter space ℬ is taken as a probability space; however, we discuss it here only in a simplified sense without reference to sigma-algebras. For the present purposes, we note only that ℬ contains all values of the input parameters that may be considered reasonable or feasible for the application. The notation h(x;β) has been used to indicate that the model output is a function of position x that depends upon the input parameters in β. For any choice of model, the vector β is simply a collection of the input parameters to the model. From the solutions h(x;β), an indicator function can be constructed to denote the region in space, Ω(β), where the solution space has been inundated:

$$\begin{array}{}\text{(2)}& {\mathbf{1}}_{\mathrm{\Omega }\left(\mathit{\beta }\right)}\left(\mathit{x}\right):=\left\{\begin{array}{ll}\mathrm{1}& \left\{\mathit{x}\left|h\right(\mathit{x};\mathit{\beta })>h_{\mathrm{0}}\right\}\\ \mathrm{0}& \left\{\mathit{x}\left|h\right(\mathit{x};\mathit{\beta })\le h_{\mathrm{0}}\right\}\end{array}\right.,\end{array}$$

where h₀ is the inundation thickness threshold of interest. For many GMFs including lahars, lava flows, and pyroclastic density currents, the threshold of interest is h₀=0 since any amount of inundation is a significant hazard; however, in the typical practice of PHM construction from numerical modeling, a nonzero thickness tolerance is typically set due to inaccuracies in the solver near free boundaries (h(x)=0) or other considerations (e.g., Patra et al., 2018b). Although doing so is inaccurate for the flows considered in this work, a sufficiently small threshold captures the essential features of the flow boundary. In assessments of other types of hazards with nonzero thresholds of interest (e.g. tephra fallout required for collapse of a particular roofing material), the true threshold may be used without introducing problems due to these solver inaccuracies.

2.1 Probabilistic hazard map

The construction of any PHM requires the definition of a joint probability distribution of the inputs given by a probability density function (PDF) $f:\mathcal{B}\to {\mathbb{R}}^{+}$ which measures the parameter space ℬ. In general, the choice of probability distribution for β will represent a combination of aleatory and epistemic uncertainties. Consequently, calculating the PHM amounts to propagating f(β) from parameter space (ℬ) through a physical model into Euclidean space. The PHM is then calculated by integration through all of the parameter space:

$$\begin{array}{}\text{(3)}& \mathit{\varphi }\left(\mathit{x}\right):=\underset{\mathcal{B}}{\int }f\left(\mathit{\beta }\right){\mathbf{1}}_{\mathrm{\Omega }\left(\mathit{\beta }\right)}\left(\mathit{x}\right)\mathrm{d}\mathit{\beta }.\end{array}$$

In the integral, the indicator function acts as a filter, removing parameter space points from the integral which do not yield solutions exceeding the threshold at a given space point. In spatial subsets where all values of β in the support of f give solutions which exceed the threshold, the indicator is identically unity and ϕ=1 by the scaling condition for a valid PDF f(β). Consequently, $\mathrm{0}\le \mathit{\varphi }\le \mathrm{1}$ globally and measures the probability of finding material in the simulated flow at position x. An alternative interpretation is the following: any given location exists in a binary state, either impacted or not impacted by the flow, which is represented spatially by an indicator function with the model acting as a mapping from a set of inputs β to the indicator 1_Ω(β)(x). Invoking the law of the unconscious statistician (LOTUS) (DeGroot and Schervish, 2012), the PHM is then a formula for computing the probability of being in either member of the binary state at each point, since the indicator function (model output) is integrated against the PDF for the model inputs. This second interpretation (using the LOTUS) is the result of encoding the binary state as a statement of certainty (probability zero or unity) that a given point is inundated in each model output.

If the range of realistic parameter values were confined to a single point β₀ in parameter space, then the parameter PDF f is the Dirac delta function δ(β−β₀), resulting in a PHM which is merely the indicator function: $\mathit{\varphi }\left(\mathit{x}\right)={\mathbf{1}}_{\mathrm{\Omega }\left({\mathit{\beta }}_{\mathrm{0}}\right)}\left(\mathit{x}\right)$. This is also the result for a PHM constructed for a single deterministic run of the model. This corresponds to the intuition that complete certainty of input parameters corresponds with complete certainty of the spatial extent of the hazard (the region where ϕ=1).

If the parameter space is measured by a uniform distribution (maximal entropy distribution), the probabilistic hazard map (PHM) is constructed by mean value integration of the indicator function through the parameter space:

$$\begin{array}{}\text{(4)}& \mathit{\varphi }\left(\mathit{x}\right)=\underset{\mathcal{B}}{\int }{\mathbf{1}}_{\mathrm{\Omega }\left(\mathit{\beta }\right)}\left(\mathit{x}\right)\mathrm{d}\mathit{\beta },\end{array}$$

where the symbol ∫_ℬ is the mean value integral. This approach is a good assumption in hazard analysis if there is no a priori knowledge of the uncertain parameters except for reasonable bounds.

In the context of a typical probabilistic hazard modeling effort, the solution space and parameter space are discrete with a uniform distribution in parameter space, and these steps are accomplished merely by cell-wise averaging of the indicator functions derived from each model run to generate a discretized 2-D hazard map ϕ=ϕ_ij. In that context, ϕ_ij is typically contoured and presented in a hazard assessment as probabilistic zones defined by the specific contours. However, because the underlying physics of these flows is not well constrained, the parameter space must be sampled widely, resulting in large probability dispersion in some cases. In the present work we globally analyze the PHM as a type of probability distribution and generate its associated measures of central tendency which form curves in 2-D space. To analyze the PHM, we must first note several of its features and state some assumptions critical to the analysis. For a PHM generated from a finite, simply connected parameter space and physically realistic governing equations,

• i.

  ϕ(x) has compact support in a region Ω⁰⊂ℝ², with ϕ=0 on ∂Ω⁰;

• ii.

  $\mathrm{0}\le \mathit{\varphi }\le \mathrm{1}$; and

• iii.

  there exists a (possibly unconnected) set Ω¹⊂Ω⁰ in which ϕ(x)=1.

We will further assume that ϕ(x) is continuous and at least twice differentiable for $\forall \mathit{x}\in {\mathbb{R}}^{\mathrm{2}}\backslash \mathit{\{}\partial {\mathrm{\Omega }}^{\mathrm{0}}\cup \partial {\mathrm{\Omega }}^{\mathrm{1}}\mathit{\}}$, that is, in the regions away from the ϕ=0 and ϕ=1 contours. We focus the analysis to follow in regions away from local maxima with ϕ<1. This restriction can be mitigated by segmentation of the spatial domain as detailed below.

Formally, ϕ(x) is a cumulative distribution function (CDF) for the location of the edge of the simulated flow. This is represented schematically in Fig. 1. To see this, consider parametric curves along level set unit normals of ϕ:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{\mathit{x}}_n}{\mathrm{d}s}}=\widehat{\mathit{n}}:={\displaystyle \frac{\mathrm{\nabla }\mathit{\varphi }}{\left|\mathrm{\nabla }\mathit{\varphi }\right|}},\end{array}$$

parameterized by arc length s with $\mathrm{0}\le s\le L$, where ${\mathit{x}}_n\left(\mathrm{0}\right)\in \partial {\mathrm{\Omega }}^{\mathrm{0}}$ and ${\mathit{x}}_n\left(L\right)\in \partial {\mathrm{\Omega }}^{\mathrm{1}}$ and L is the finite total Euclidean path length of x_n. Note that $\widehat{\mathit{n}}$ is only calculable where $\left|\mathrm{\nabla }\mathit{\varphi }\right|\ne \mathbf{0}$. To ameliorate this, we define $\widehat{\mathit{n}}$ piecewise as

$$\begin{array}{}\text{(6)}& \widehat{\mathit{n}}\left(\mathit{x}\right):=\left\{\begin{array}{ll}\mathrm{\nabla }\mathit{\varphi }/\left|\mathrm{\nabla }\mathit{\varphi }\right|& \mathrm{where}\left|\mathrm{\nabla }\mathit{\varphi }\right|\ne \mathbf{0}\\ \mathrm{0}& \mathrm{elsewhere}\end{array}\right..\end{array}$$

Figure 1 Schematic representation of a hypothetical probabilistic hazard map (PHM) where contours (blue) are the percent of deterministic models inundating a given point in space. Steepest ascent integral curves on the PHM (x_n(s)) are shown in red. Profiles for the cumulative distribution ϕ(x_n(s)) (blue) and the probability density dϕ(x_n(s))∕ds (black dashes) are represented schematically. The dark grey transparent region is the region of influence for the local maximum (ℛ_m) in which all integral curves terminate at the local maximum. Two integral curves (one inside of ℛ_m, one outside of ℛ_m) show the typical behavior of x_n(s) near ∂ℛ_m. The set ${\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}$ is the entire white region, Ω⁰ is the region bounded by the zero level set (the white region plus the dark grey region), and Ω¹ is the subset of that region which is marked with 100 %.

Along each integral curve, profiles of the PHM are denoted as

$$\begin{array}{}\text{(7)}& \mathit{\varphi }\left(s\right):=\mathit{\varphi }\left({\mathit{x}}_n\left(s\right)\right)=\mathbb{P}\left(s_{\mathrm{F}}\le s\right),\end{array}$$

where ℙ(s_F≤s) is the cumulative probability that position x_n(s) is farther from x_n(0) than the true inundation front x_n(s_F) is from x_n(0) as measured along the integral curve. For example, where ϕ(s)=1, the inundation front is almost surely found somewhere further back on the curve towards x_n(0). Similarly, where $\mathit{\varphi }\left(s\right)=\frac{\mathrm{1}}{\mathrm{2}}$, there is a 50 % chance of finding the true inundation front between that point (x_n(s)) and either of the two endpoints of the integral curve (x_n(0) and x_n(L)). Because the integral curves follow the gradient, ϕ(s) is nondecreasing and along with $\mathrm{0}\le \mathit{\varphi }\left(s\right)\le \mathrm{1}$ satisfies the description of a CDF.

If local maxima less than unity are present, then the integral curves that ascend these maxima will not be valid CDFs (ϕ(s)) because they never reach the value ϕ=1. Despite this, as long as the PHM contains a finite number of local maxima less than unity (ϕ(x_m)<1 where the index “m” runs over the M local maxima), then for each local maximum we may define a region ℛ_m which surrounds x_m and is bounded by two integral curves as well as a segment of ∂Ω⁰ and ∂Ω¹. Because local maxima are sinks for gradient-ascending integral curves, all of the integral curves that begin on $\partial {\mathcal{R}}_{\mathrm{m}}\cap \partial {\mathrm{\Omega }}^{\mathrm{0}}$ will terminate at x_m by definition. The two bounding integral curves are therefore defined extrinsically as those integral curves producing valid CDFs originating closest to the segment of ∂Ω⁰ that originates the integral curves terminating at x_m. Although these regions are only defined in a weak sense, these regions are definable for practical purposes (Appendix B). If there is a saddle point between a local maximum and ∂Ω¹, then the two integral curves will meet at the saddle point and no part of ∂Ω¹ will make up the boundary of ℛ_m. In a sense, ℛ_m can be thought of as the shadow of trajectories that are influenced most by the local maximum ϕ(x_m) (Fig. 1). We remove these regions by constructing a new set on which the remainder of the analysis is valid:

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}={\mathrm{\Omega }}^{\mathrm{0}}\backslash \bigcup _{m=\mathrm{1}}^M{\mathcal{R}}_{\mathrm{m}}.\end{array}$$

2.2 Probabilistic hazard density map

To find the PDF dϕ∕ds associated with this distribution, consider that by the chain rule

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\mathit{\varphi }\left(s\right)}{\mathrm{d}s}}=\left(\mathrm{\nabla }\mathit{\varphi }\cdot {\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}s}}\right){|}_{{\mathit{x}}_n\left(s\right)}=\left|\mathrm{\nabla }\mathit{\varphi }\left({\mathit{x}}_n\left(s\right)\right)\right|,\end{array}$$

and consequently

$$\begin{array}{}\text{(10)}& \underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}s}}\mathrm{d}s=\underset{\mathrm{0}}{\overset{L}{\int }}\left|\mathrm{\nabla }\mathit{\varphi }\left({\mathit{x}}_n\right(s\left)\right)\right|\mathrm{d}s=\mathrm{1}.\end{array}$$

Although this description satisfies the requirements of a PDF along each individual integral curve, it cannot be applied directly to measure the probability density in two dimensions. To measure this probability density, we generate a normalized PDF over 2-D space which we refer to as the probabilistic hazard density map (PHDM):

$$\begin{array}{}\text{(11)}& \mathit{\psi }\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{Q}}\left|\mathrm{\nabla }\mathit{\varphi }\right|,\end{array}$$

where

$$\begin{array}{}\text{(12)}& Q=\underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\left|\mathrm{\nabla }\mathit{\varphi }\right|\mathrm{d}A=-\underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\mathit{\varphi }\mathrm{\nabla }\cdot \widehat{\mathit{n}}\mathrm{d}A,\end{array}$$

which is finite and positive in the sense of distributions. With this normalization, the integral over ${\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}$ is unity, giving ψ(x) the properties of a PDF for the inundation front location. Consequently, the value ψ(x)dA represents the probability of finding the inundation edge in an increment of area dA. As noted, this distribution is only valid in ${\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}$: the region of space where local maxima in the PHM do not influence the integral curves. A proof which extends Eq. (10) to this normalization of the PHDM is given in Appendix A.

In the notional example given above for a PHM generated by a parameter space Dirac delta, given a PHM that is identically the indicator function, the PHDM is defined only in the sense of distributions because of the apparently infinite gradient at ∂dΩ(β₀). However, because derivatives can be defined in the sense of distributions, the PHDM corresponding to certain input parameters is $\mathit{\psi }\left(\mathit{x}\right)=\mathit{\delta }(\mathit{x}-{\mathit{x}}^{*})$, where ${\mathit{x}}^{*}\in \partial \mathrm{\Omega }\left({\mathit{\beta }}_{\mathrm{0}}\right)$ form the locus of points representing the certain boundary location. In this case, Ω⁰=Ω(β₀). This scenario corresponds to the propagation of complete certainty in parameter space to complete certainty of the flow boundary location in Euclidean space.

2.3 Measures of central tendency and variance

Using these definitions, three measures of central tendency (mode, mean, median) can be generated for each integral curve, generating a locus of points (curve) for each measure parameterized by a new parameter (l) which increases along the PHM contours according to

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}{\mathit{x}}_t}{\mathrm{d}l}}=\mathbf{R}\widehat{\mathit{n}},\end{array}$$

where $\mathbf{R}=\left(\begin{array}{c}\mathrm{0}-\mathrm{1}\\ \mathrm{1}\mathrm{0}\end{array}\right)$ is a rotation matrix and x_t(l) are integral curves everywhere tangent to the PHM contours (schematized in Fig. 1).

By defining a density function for the inundation edge probability distribution as in Eq. (11), a notional modal set can be constructed by connecting the maximal value of ψ on each integral curve x_n. The modal set (parameterized by l) is

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{\mathit{x}}_{\mathcal{M}}\left(l\right):=\left\{\mathit{x}\in {\mathbb{R}}^{\mathrm{2}}\left|\mathit{\psi }\right(\mathit{x})={\max }_{s\in [\mathrm{0},L]}\mathit{\psi }\left({\mathit{x}}_n(s;l)\right)\right\}.\end{array}$$

For the inundation hazard map, the curve x_ℳ(l) represents the maximum likelihood curve along which to find the inundation front. Although each distribution may not be unimodal and individual distributions may have unimodal and multimodal regions owing to the complexity of the underlying physics of the problem, this estimate gives a starting point for understanding the central tendency of the distribution.

Furthermore, we may find the mean location of the inundation by considering the first moment of each univariate PDF $\frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}s}(s;l)$. For each integral curve x_n(s;l) connecting a point on ∂Ω⁰ to a point on ∂Ω¹, the first moment (mean) of dϕ∕ds may be written as

$$\begin{array}{}\text{(15)}& \mathit{\mu }\left(l\right)=\underset{\mathrm{0}}{\overset{L}{\int }}s{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}s}}(s;l)\mathrm{d}s.\end{array}$$

We compute this integral by parts but make a critical substitution before doing so. Note that for a new function 1−ϕ, differentiating yields d$\mathit{\varphi }=-$d(1−ϕ). We may substitute this into Eq. (15) and then integrate by parts:

$$\begin{array}{}\text{(16)}& \mathit{\mu }\left(l\right)=-s\left[\mathrm{1}-\mathit{\varphi }(s;l)\right]{|}_{s=\mathrm{0}}^{s=L}+\underset{\mathrm{0}}{\overset{L}{\int }}\mathrm{1}-\mathit{\varphi }(s;l)\mathrm{d}s.\end{array}$$

Because ϕ(L)=1, the endpoint-evaluated terms in the limit will vanish. This allows for a simpler representation for the first moment:

$$\begin{array}{}\text{(17)}& \mathit{\mu }\left(l\right)=\underset{\mathrm{0}}{\overset{L}{\int }}\mathrm{1}-\mathit{\varphi }\left({\mathit{x}}_n(s;l)\right)\mathrm{d}s.\end{array}$$

This formula gives the expected value of the arc length parameter, which can be substituted to give an expected boundary location. Once μ(l) is calculated for every integral curve x_n(s;l), we may connect these points to form a curve defined by

$$\begin{array}{}\text{(18)}& {\mathit{x}}_{\mathit{\mu }}\left(l\right):=\left\{\mathit{x}\in {\mathbb{R}}^{\mathrm{2}}\left|\mathit{x}\right.={\mathit{x}}_n\left(\mathit{\mu }\right(l\left)\right)\right\}.\end{array}$$

Finally, we may define a curve representing the median value of the probability distribution for inundation front location. This is simply the 50 % probability contour:

$$\begin{array}{}\text{(19)}& {\mathit{x}}_{\mathrm{med}}\left(l\right):=\left\{\mathit{x}\in {\mathbb{R}}^{\mathrm{2}}\left|\mathit{\varphi }\left(\mathit{x}\right)\right.={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right\},\end{array}$$

which is a level set of the hazard map. The three curves x_ℳ, x_μ, and x_med represent different measures of the central tendency of the hazard map probability distribution.

Additionally, a measure of the variance may also be constructed for each PDF $\frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}s}(s;l)$, although this requires integration of the density function over s:

$$\begin{array}{}\text{(20)}& {\mathit{\sigma }}^{\mathrm{2}}\left(l\right)=\underset{\mathrm{0}}{\overset{L}{\int }}{s}^{\mathrm{2}}{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}s}}(s;l)\mathrm{d}s-\mathit{\mu }(l{)}^{\mathrm{2}}.\end{array}$$

The standard deviation (σ) defines a region around the mean value in that for each integral curve, the points $s=\mathit{\mu }\left(l\right)\pm \mathit{\sigma }\left(l\right)$ map to points in space along the curves x_μ±σ(l), forming a region about  x_μ(l) between x_μ+σ(l) and x_μ−σ(l). Higher moments of the distribution may be calculated similarly along each curve and subsequently connected.

3 Results

To highlight the features of the new theory, we present three examples of its application in order of increasing complexity and realism.

3.1 Example 1: analytic model of time-dependent solidification of a viscous gravity current

A simple application of this theory is to a flow which can be calculated analytically. Here, we show the application of the above theory to the problem of time-dependent solidification of a viscous gravity current. This model used by Quick et al. (2017) is a modification of the famous model due to Huppert (1982) of the axisymmetric (radial) spread of viscous gravity currents, allowing for time-dependent viscosity where the height of the current is governed by

$$\begin{array}{}\text{(21)}& {\partial }_{t}h-{\displaystyle \frac{g}{\mathit{\nu }\left(t\right)}}{\displaystyle \frac{\mathrm{1}}{r}}{\partial }_r\left[r{h}^{\mathrm{3}}{\partial }_{r}h\right]=\mathrm{0},\end{array}$$

where g is gravitational acceleration. This model incorporates a simple model of time-dependent kinematic viscosity increase representing the solidification of the current,

$$\begin{array}{}\text{(22)}& \mathit{\nu }\left(t\right)={\mathit{\nu }}_{\mathrm{0}}{e}^{t/\mathrm{\Gamma }},\end{array}$$

where ν₀ is the kinematic viscosity (ν) at time t=0 and Γ is the time it takes for the viscosity to increase by a factor of e. The solution for the release of a constant volume V₀ with initial radius r₀ may be stated as

$$\begin{array}{}\text{(23)}& h(r,t)={\displaystyle \frac{\mathrm{4}{V}_{\mathrm{0}}}{\mathrm{3}\mathit{\pi }{r}_{\mathrm{0}}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+\mathrm{\Theta }(t)/\mathit{\tau }{)}^{\mathrm{1}/\mathrm{4}}}}{\left[\mathrm{1}-{\displaystyle \frac{r^{\mathrm{2}}}{r_{\mathrm{0}}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+\mathrm{\Theta }(t)/\mathit{\tau }{)}^{\mathrm{1}/\mathrm{4}}}}\right]}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where

$$\begin{array}{}\text{(24)}& \mathrm{\Theta }\left(t\right)=\mathrm{\Gamma }\left(\mathrm{1}-e^{-t/\mathrm{\Gamma }}\right)\end{array}$$

and

$$\begin{array}{}\text{(25)}& \mathit{\tau }={\left({\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\right)}^{\mathrm{5}}{\left({\displaystyle \frac{\mathit{\pi }}{V_{\mathrm{0}}}}\right)}^{\mathrm{3}}{\displaystyle \frac{{\mathit{\nu }}_{\mathrm{0}}{r}_{\mathrm{0}}^{\mathrm{8}}}{g}}.\end{array}$$

Notably, Θ→Γ as t→∞, implying that eventually the current will solidify and grow to a maximum size. Consequently, when the flow dynamics have ceased, the radius of the edge of the current (r_N), where the current thickness vanishes, will have reached a maximum at

$$\begin{array}{}\text{(26)}& r_N=r_{\mathrm{0}}(\mathrm{1}+\mathrm{\Gamma }/\mathit{\tau }{)}^{\mathrm{1}/\mathrm{8}}.\end{array}$$

To use this simple model to illustrate the probability distribution properties above, consider the case where all parameters are known precisely except for the viscosity freezing time coefficient Γ which can be constrained within a range ${\mathrm{\Gamma }}_{\min }<\mathrm{\Gamma }<{\mathrm{\Gamma }}_{\max }$. Because Eqs. (23) and (26) could be nondimensionalized so as to be rewritten in terms of one fundamental parameter $\mathit{\eta }:=\mathrm{\Gamma }/\mathit{\tau }$, which is proportional to Γ, this single random parameter choice captures the underlying uncertainty in this example. Under the assumption of no prior knowledge of the likelihood of any particular value of Γ, a uniform distribution may be used to construct an inundation indicator function for this model:

$$\begin{array}{}\text{(27)}& {\mathbf{1}}_{\mathrm{\Omega }\left(\mathrm{\Gamma }\right)}\left(r\right)=\left\{\begin{array}{ll}\mathrm{1}& r<r_N\left(\mathrm{\Gamma }\right)\\ \mathrm{0}& r\ge r_N\left(\mathrm{\Gamma }\right)\end{array}\right..\end{array}$$

To generate a PHM from this indicator function as above, the indicator function must be integrated through the parameter space. The only nonconstant segment of this integral is on the interval $r_N\left({\mathrm{\Gamma }}_{\min }\right)<r<r_N\left({\mathrm{\Gamma }}_{\max }\right)$ where

$$\begin{array}{}\text{(28)}& \mathit{\varphi }\left(r\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }\mathrm{\Gamma }}}\underset{{\mathrm{\Gamma }}_N\left(r\right)}{\overset{{\mathrm{\Gamma }}_{\max }}{\int }}\mathrm{d}\mathrm{\Gamma },\end{array}$$

with ${\mathrm{\Gamma }}_N\left(r\right)=\mathit{\tau }\left[\right(r/r_{\mathrm{0}}{)}^{\mathrm{8}}-\mathrm{1}]$ and $\mathrm{\Delta }\mathrm{\Gamma }={\mathrm{\Gamma }}_{\max }-{\mathrm{\Gamma }}_{\min }$. This results in the PHM (Fig. 2a and c):

$$\begin{array}{}\text{(29)}& \mathit{\varphi }\left(r\right)=\left\{\begin{array}{ll}\mathrm{1}& r\le r_N\left({\mathrm{\Gamma }}_{\min }\right)\\ \frac{{\mathrm{\Gamma }}_{\max }+\mathit{\tau }-\mathit{\tau }{\left(r/r_{\mathrm{0}}\right)}^{\mathrm{8}}}{{\mathrm{\Gamma }}_{\max }-{\mathrm{\Gamma }}_{\min }}& r_N\left({\mathrm{\Gamma }}_{\min }\right)<r<r_N\left({\mathrm{\Gamma }}_{\max }\right)\\ \mathrm{0}& r\ge r_N\left({\mathrm{\Gamma }}_{\max }\right)\end{array}\right..\end{array}$$

From the PHM, the PHDM (Fig. 2b and c) can be calculated as

$$\begin{array}{}\text{(30)}& \mathit{\psi }\left(r\right)={\displaystyle \frac{\mathrm{1}}{Q}}\left|\mathrm{\nabla }\mathit{\varphi }\right|=\left\{\begin{array}{ll}\frac{\mathrm{9}}{\mathrm{2}\mathit{\pi }}\frac{r^{\mathrm{7}}}{r_N{\left({\mathrm{\Gamma }}_{\max }\right)}^{\mathrm{9}}-r_N{\left({\mathrm{\Gamma }}_{\min }\right)}^{\mathrm{9}}}& r_N\left({\mathrm{\Gamma }}_{\min }\right)<r<r_N\left({\mathrm{\Gamma }}_{\max }\right)\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.,\end{array}$$

and the integral curves can be calculated as

$$\begin{array}{}\text{(31)}& r\left(s\right)=r_N\left({\mathrm{\Gamma }}_{\max }\right)-s,\end{array}$$

with $\mathrm{0}\le s\le L$ and $L=r_N\left({\mathrm{\Gamma }}_{\max }\right)-r_N\left({\mathrm{\Gamma }}_{\min }\right)$.

Figure 2 (a) PHM and (b) PHDM examples for the freezing viscous gravity current model under a uniform sampling of the freezing rate parameter Γ. (c) One-dimensional profile of the PHM and PHDM with the distribution centers indicated. This example was calculated with $r_N\left({\mathrm{\Gamma }}_{\min }\right)/r_{\mathrm{0}}=\mathrm{2}$ and $r_N\left({\mathrm{\Gamma }}_{\max }\right)/r_{\mathrm{0}}=\mathrm{3}$ or ${\mathrm{\Gamma }}_{\min }/\mathit{\tau }=\mathrm{255}$ and ${\mathrm{\Gamma }}_{\max }/\mathit{\tau }=\mathrm{6560}$.

These relationships lead naturally to the estimates of the center of the distribution (Fig. 2) as defined above,

$$\begin{array}{}\text{(32)}& \left\{\begin{array}{l}{r}_{\mathcal{M}}=r_N\left({\mathrm{\Gamma }}_{\max }\right)\\ r_{\mathrm{med}}=r_{\mathrm{0}}{\left[\mathrm{1}+{\displaystyle \frac{{\mathrm{\Gamma }}_{\max }+{\mathrm{\Gamma }}_{\min }}{\mathrm{2}\mathit{\tau }}}\right]}^{\mathrm{1}/\mathrm{8}}\\ r_{\mathit{\mu }}=r_N\left({\mathrm{\Gamma }}_{\max }\right)-\mathit{\mu }\end{array}\right.,\end{array}$$

and the standard deviation,

$$\begin{array}{}\text{(33)}& r_{\mathit{\mu }\pm \mathit{\sigma }}=r_N\left({\mathrm{\Gamma }}_{\max }\right)-\mathit{\mu }\mp \mathit{\sigma },\end{array}$$

where

$$\begin{array}{}\text{(34)}& \left\{\begin{array}{l}\mathit{\mu }=r_N\left({\mathrm{\Gamma }}_{\min }\right)\frac{\mathit{\tau }}{\mathrm{\Delta }\mathrm{\Gamma }}\left[\frac{r_N{\left({\mathrm{\Gamma }}_{\min }\right)}^{\mathrm{8}}}{\mathrm{9}{r}_{\mathrm{0}}^{\mathrm{8}}}\left({\mathit{\rho }}^{\mathrm{9}}-\mathrm{1}\right)-\frac{{\mathrm{\Gamma }}_{\min }+\mathit{\tau }}{\mathit{\tau }}(\mathit{\rho }-\mathrm{1})\right]\\ \mathit{\sigma }=\left\{r_N{\left({\mathrm{\Gamma }}_{\min }\right)}^{\mathrm{2}}\frac{\mathit{\tau }}{\mathrm{\Delta }\mathrm{\Gamma }}\left[\frac{r_N{\left({\mathrm{\Gamma }}_{\min }\right)}^{\mathrm{8}}}{\mathrm{45}{r}_{\mathrm{0}}^{\mathrm{8}}}\left({\mathit{\rho }}^{\mathrm{10}}-\mathrm{1}\right)-\frac{\mathrm{2}{r}_N{\left({\mathrm{\Gamma }}_{\min }\right)}^{\mathrm{8}}}{\mathrm{9}{r}_{\mathrm{0}}^{\mathrm{8}}}\left(\mathit{\rho }-\mathrm{1}\right)\right.\right.\\ {\left.\left.-\frac{{\mathrm{\Gamma }}_{\min }+\mathit{\tau }}{\mathit{\tau }}(\mathit{\rho }-\mathrm{1}{)}^{\mathrm{2}}\right]-{\mathit{\mu }}^{\mathrm{2}}\right\}}^{\mathrm{1}/\mathrm{2}}\\ \mathrm{and}:\\ \mathit{\rho }=\frac{r_N\left({\mathrm{\Gamma }}_{\max }\right)}{r_N\left({\mathrm{\Gamma }}_{\min }\right)}\end{array}\right.,\end{array}$$

where Eq. () was substituted into Eq. (), and Eqs. (15) and (20) were used to calculate Eq. (34). Higher moments of the distribution are theoretically calculable by repeated integration by parts, from which additional facts may be learned about the distribution.

These calculations can be made analytically because of the ability to write the equation of the boundary and because that equation is invertible for the unknown parameter; that is, it can be written as r_N(Γ) (a boundary in Euclidean space) and as Γ(r_N) (a boundary in parameter space). This example also serves to highlight the asymmetry of the PHM despite the underlying uniform distribution of physical parameters. This asymmetry is highlighted by the significant spread in the mean, median, and modal boundaries (Fig. 2).

This example was greatly simplified by the assumption of radial symmetry, allowing the distribution to be cast onto one random variable, the radius of the flow boundary, giving integral curves of constant azimuthal angle which are linear in map view. However, in most realistic cases, this is a very poor assumption and the distribution will vary between integral curves.

3.2 Example 2: Mohr–Coulomb flow on an inclined plane

To see the full realization of these concepts, a true two-dimensional problem is required. A simple example flow which had been well-studied is a lab-scale granular flow of sand down an inclined plane with a Mohr–Coulomb (MC) friction relation analogous to a natural-scale debris avalanche or landslide (Patra et al., 2005, 2018a, b; Dalbey et al., 2008; Yu et al., 2009; Webb and Bursik, 2016; Aghakhani et al., 2016). In the simulated version, a MC rheology variant of the Saint-Venant shallow water equations is solved numerically with the purpose-built solver TITAN2D (available from https://vhub.org/, Version 4.0.0, VHUB, 2016). In this setup (Fig. 3a), a 5 cm radius, 5 cm high cylindrical pile is placed 0.7 m up a 1 m long, 38.5^∘ inclined plane and is allowed to flow out onto a flat runout 0.5 m long by 0.7 m wide. In accordance with MC theory, the flows were stopped locally when the sum of the drag due to internal and bed friction exceeded the gravitational forces, and the simulations were stopped globally after 1.5 s, when the flow dynamics had typically ceased (Patra et al., 2018b). In this case, the solution variable of interest h is the maximum flow height over time, that is, the maximum flow height that each point experiences through the course of the flow.

Figure 3 (a) PHM contours for a Mohr–Coulomb (MC) flow on an inclined plane with integral curves x_n(s) (red) connecting points on ∂Ω⁰ to points on ∂Ω¹ as defined in the text. (b) Mean (x_μ, red curve) and median (x_med, blue curve) flow boundary sets as well as regions within 1 standard deviation from the mean (red) and containing 68 % of the possible flow boundaries around the median (blue).

The uncertain parameters in the problem for a given granular medium are the internal and basal friction angles of the medium, φ_int and φ_bed respectively. To construct the PHM for this flow, 512 TITAN2D simulations were performed over a Latin-hypercube-sampled input parameter space of

$$\begin{array}{}\text{(35)}& \mathcal{B}=\left\{\left({\mathit{\phi }}_{\mathrm{int}},{\mathit{\phi }}_{\mathrm{bed}}\right)\left|\mathrm{18}{}^{\circ }\le {\mathit{\phi }}_{\mathrm{bed}}\le \mathrm{30}{}^{\circ }\right.\cap \mathrm{2}{}^{\circ }\le {\mathit{\phi }}_{\mathrm{int}}-{\mathit{\phi }}_{\mathrm{bed}}\le \mathrm{10}{}^{\circ }\right\}.\end{array}$$

To construct the indicator functions for this experiment, a thickness threshold of h₀=1.0 mm was adopted. The PHM was then computed by averaging the indicator functions: a discrete analogy to Eq. (4). Owing to the discrete nature of the data, all continuous operations on the PHM were done by numerical finite differences and quadrature.

The PHM for this problem has many important characteristics, notably that the greatest dispersion in the probability contours occurs in the runout zone, owing principally to the uncertainty in φ_bed. Additionally, the spreading angle of the flow as a whole experiences a significant increase as the flow passes the break in slope (Fig. 3a). In order to calculate the mean flow boundary set and higher moments of this distribution, the gradient ascent curves of Eq. (5) were integrated numerically from a finite set of closely spaced points around an approximation of ∂Ω⁰ defined as

$$\begin{array}{}\text{(36)}& \widehat{\partial {\mathrm{\Omega }}^{\mathrm{0}}}:=\left\{{\mathit{x}}^j\left|\mathit{\varphi }\left({\mathit{x}}_j\right)\right.={\mathrm{10}}^{-\mathrm{4}},j\in [\mathrm{1},J]\right\},\end{array}$$

where J=2049 is the number of such points generated by numerical contouring. Equation (5) was integrated by Euler's method with a 1 mm step forward in s (Δs_k=1 mm), thus fixing $\left|\right|{\mathit{x}}_n^j\left(s_{k+\mathrm{1}}\right)-{\mathit{x}}_n^j\left(s_k\right)\left|\right|=\mathrm{1}$ mm. This yielded discrete integral curves ${\mathit{x}}_n^j\left(s_k\right)$ of varying lengths around the perimeter of the PHM (Fig. 3a). With each ${\mathit{x}}_n^j\left(s_k\right)$ calculated, ϕ(s_k;l^j) was constructed and the parameter μ^j calculated on each curve, yielding an estimate of x_μ(l) by interpolation of the discrete points ${\mathit{x}}_{\mathit{\mu }}^j$ (Fig. 3b). The standard deviations σ^j and the standard deviation space curves x_μ±σ(l) were calculated by similar means (Fig. 3b). An important result of these calculations is that for almost every part of the runout zone, the mean flow boundary set x_μ is dislocated from the median flow boundary set x_med ($\mathit{\varphi }=\frac{\mathrm{1}}{\mathrm{2}}$), indicating qualitatively the general asymmetry of this distribution. Additionally the region bounded by the curves x_μ±σ is remarkably different from the region which bounds the central 68 % of the data – the percent bounded by ±σ in the normal distribution (Fig. 3). All of this indicates that for PHMs of even simple realistic flows, the underlying statistics are not only nonnormal and nonsymmetric, they are variable along the boundary of the hazard.

3.3 Example 3: PHM and PHDM construction from simulations of the 1955 debris flow in Atenquique, Mexico

Following from the previous application of the PHM statistics to a simple flow, here we present a more advanced application: construction and analysis of a PHM for a more complex flow over natural topography. This example consists of numerical modeling detailed in Bevilacqua et al. (2019) of the 1955 volcaniclastic debris flow which destroyed the village of Atenquique at the base of Nevado de Colima, Mexico. On 16 October 1955, at 10:45 LT, residents of Atenquique experienced the arrival of an 8–9 m high debris flow wave front that subsequently destroyed much of the town's buildings and infrastructure and caused more than 23 deaths (Ponce Segura, 1983; Saucedo et al., 2008). The debris flow was initiated by landsliding in the upper reaches of the Atenquique, Arroyo Seco, Los Plátanos, and Tamazula fault ravines, with a minimum initial volume of 3.2×10⁶ m³ (Rupp, 2004; Saucedo et al., 2008). In the simulation of this flow, this volume was increased by 10 %–50 % to give volumes of $V=\mathrm{4.25}\pm \mathrm{0.75}\times {\mathrm{10}}^{\mathrm{6}}$ m³ (Bevilacqua et al., 2019).

We analyze an ensemble of TITAN2D simulations of this flow that were generated by Bevilacqua et al. (2019) using a physically realistic rheology for a debris flow, the Voellmy–Salm (VS) rheology. This rheology is parameterized by a bed friction coefficient μ_VS and a velocity-dependent friction parameter ξ which approximates turbulence-induced dissipation and is uncertain over several orders of magnitude. Consequently, ξ is sampled logarithmically (Patra et al., 2018a, b). The flow was initiated from five source paraboloidal piles of material with unitary aspect ratios and volumes represented as constant fractions of the total flow volume which were proportional to the drainage size (Fig. 4a). The input parameter space ℬ for this simulation was taken as a convex subset of {(arctan (μ_VS), log ₁₀(ξ), V)} with the bounds 1^∘ ≤ arctan (μ_VS)≤ 3^∘, $\mathrm{3}\le {\log }_{\mathrm{10}}\left(\mathit{\xi }\right)\le \mathrm{4}$, and $\mathrm{3.5}\times {\mathrm{10}}^{\mathrm{6}}\le V\le \mathrm{5.0}\times {\mathrm{10}}^{\mathrm{6}}$. This set was determined using information about the known flow runout, flow arrival time in the town, and exclusion of ravine overflow in conjunction with a round of preliminary simulations according to a new procedure of probabilistic modeling setup (Bevilacqua et al., 2019).

Figure 4 (a) PHM and (b) PHDM for a lahar inundating Atenquique, Mexico, using the numerical model TITAN2D incorporating the parameter space described in the text. (c) Profiles of the PHM and PHDM across sections X–X^′ and Y–Y^′ highlighting the fact that each side of the PHM profile acts as a CDF for the location of the inundation edge with the corresponding scaled PDF given by the same profile along the PHDM. Grey background (c) indicates regions where ϕ=0 or ϕ=1 (ψ=0).

As in the previous example, the solution variable of interest is the maximum flow height over time; however, here the height threshold h₀=5.0 cm is used. The 350 VS simulation indicator functions were averaged as in the previous example to obtain the PHM for this flow (Fig. 4a). For such a complex flow margin, the gradient-ascending integral curves could not be computed reliably at the resolution of the model; however, the PHDM was calculated using a discrete gradient operation, giving a sense of the maximum likelihood for the inundation boundary (Fig. 4b). Owing to the complexity of this flow and its underlying topography, profiles of the PHM and PHDM across the flow indicate that each margin of the flow is distributed nonsymmetrically, although the sense of asymmetry does not appear to be systematic (Fig. 4c). In numerous locations, the flow boundary distributions (profiles of the PHDM) are multimodal, which are interpreted to be controlled primarily by complexities in the underlying topography, relating to partial spillover of the steep bounding ravines. Of particular interest is the behavior of the PHM and PHDM at the distal margin of the flow in the town of Atenquique and as it leaves town downriver to the south and upriver to the north. To the south (downriver) the PHM decays over a very long distance (>5 km) down a constrained channel, giving small values for the PHDM along the direction of flow (Fig. 4b). Consequently in this region, the flow boundary distributions are multimodal with weaker modes more distally and only the lateral boundaries being well-defined. This portion of the flow highlights some important problems with PHMs in general; they may include regions where the distribution spans a wide range of distances and will therefore have a large standard deviation, yielding high uncertainty in the flow boundary location. The accuracy of the PDFs over the integral curves in this region is most sensitive to the number of samples affecting this region. In a finite sampling of the parameter space, these PDFs may contain spurious peaks in regions where very few models reach. For example if only a small number of models reach a region where the PHM decays over many grid cells, unless some smoothing or interpolation is applied before gradient calculation, the PHM will decay stepwise with many grid spaces between steps. When the numerical gradient is calculated over a set cell size, these steps will become peaks in the PDFs and PHDM. This problem is generally absent where the PHM step size is more comparable to the scale of the grid spacing used by the gradient operator.

This complex example highlights the need for consideration of additional probability measures than simply the probability contours. Where the flow is not narrowly confined or where spillover has occurred, the shape of the distribution becomes relevant and may yield significantly different results than estimating the likely flow boundary as one of the probability contours.

4 Discussion

As shown by the above three examples, our theory can be applied to give consistent results between continuous, analytic models and discrete, numerical models of flows. The following discussion highlights several consequences of the theory as applied to analytic models in research and more complex numerical models in hazard assessments.

4.1 Analytic models

Although an analytic model (Quick et al., 2017) was invoked in the first example to make an exact calculation of the mean, median, mode, and standard deviation of the flow boundary, the real value of these models is that they are generally invertible and can be used to extract information about the input parameters from the properties of a given flow. Specifically, the PDF measuring the space of input parameters f can be inverted from data collected from the natural phenomena that the analytic model represents. Using the example of the freezing viscous gravity current highlighted above (Quick et al., 2017), we show that an ensemble of natural examples of such flows can be used to invert the PDF of input parameters directly from the definition of the PHM. Suppose that remote sensing or geological field work yields data on the maximum radii (r_N) for an ensemble of these natural freezing gravity currents and that a natural PHM model is constructed called $\widehat{\mathit{\varphi }}\left(r\right)$, which represents the discrete CDF of radii that were measured. As before, we consider a simplified case wherein the most uncertain parameter from among the collection of natural flows is the freezing rate parameter Γ. If the researcher were sufficiently confident that these assumptions and the governing equation (Eq. 21) accurately describe these flows, then the only differences between $\widehat{\mathit{\varphi }}\left(r\right)$ and Eq. () would arise as a consequence of the PDF f(Γ) of the natural input space being nonuniform. Consequently, because the formula for the flow boundary ($r_N\left(\mathrm{\Gamma }\right)=r_{\mathrm{0}}(\mathrm{1}+\mathrm{\Gamma }/\mathit{\tau }{)}^{\mathrm{1}/\mathrm{8}}$) is known and is invertible, the formula for the PHM can also be inverted directly in this simple case as follows:

$$\begin{array}{}\text{(37)}& \widehat{\mathit{\varphi }}\left(r_N\left(\mathrm{\Gamma }\right)\right)=\underset{{\mathrm{\Gamma }}_{\min }}{\overset{{\mathrm{\Gamma }}_{\max }}{\int }}f\left(\mathrm{\Gamma }\right){\mathbf{1}}_{\mathrm{\Omega }\left(\mathrm{\Gamma }\right)}\left(r_N\left(\mathrm{\Gamma }\right)\right)\mathrm{d}\mathrm{\Gamma }=\underset{\mathrm{\Gamma }\left(r_N\right)}{\overset{{\mathrm{\Gamma }}_{\max }}{\int }}f\left({\mathrm{\Gamma }}^{\prime }\right)\mathrm{d}{\mathrm{\Gamma }}^{\prime }.\end{array}$$

Applying the fundamental theorem of calculus and the chain rule yields

$$\begin{array}{}\text{(38)}& f\left(\mathrm{\Gamma }\right)=-{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathrm{\Gamma }}}\widehat{\mathit{\varphi }}\left(r_N\left(\mathrm{\Gamma }\right)\right)=-{\displaystyle \frac{\mathrm{d}{r}_N}{\mathrm{d}\mathrm{\Gamma }}}\left(\mathrm{\Gamma }\right)\cdot {\displaystyle \frac{\mathrm{d}\widehat{\mathit{\varphi }}}{\mathrm{d}{r}_N}}\left(r_N\left(\mathrm{\Gamma }\right)\right).\end{array}$$

In practice, the collected data would be in the form of a normalized histogram $\widehat{p}\left(r_N\right)$ for the maximum flow radii. As a discrete approximation, $\widehat{\mathit{\varphi }}\left(r_N\right)$ is the associated complementary CDF for $\widehat{p}$. Consequently, the PDF over parameter space could be written approximately as

$$\begin{array}{}\text{(39)}& f\left({\mathrm{\Gamma }}_i\right)={\displaystyle \frac{\mathrm{d}{r}_N}{\mathrm{d}\mathrm{\Gamma }}}\left({\mathrm{\Gamma }}_i\right)\cdot \widehat{p}\left(r_N\left({\mathrm{\Gamma }}_i\right)\right),\end{array}$$

where Γ_i are sample points at which the calculation is performed. This illustrates the general fact that the distribution of flow boundaries observed reflects the underlying parameter space PDF weighted by the Jacobian (univariate derivative in this case) of the analytic mapping from parameter space to the flow boundary location in Euclidean space. Generally, unless the analytic model predicts that the boundary of a flow (Ω(β)) is proportional to a given parameter, the input space distribution f(β) undergoes a nonlinear transformation in the construction of a PHM and PHDM. More precisely, we remark that the calculated distribution is the conditional probability density f(Γ|β_k), where β_k represents all of the other parameters that were assumed to be known precisely. Inverting the PHM for the joint parameter space PDF is beyond the scope of this work.

4.2 Numerical models and probabilistic hazard assessment

If the parameter space of a given inundation hazard scenario and a realistic physical model of the process are known, a simple numerical experiment can be run using classical Monte Carlo methods or other Monte Carlo variants as was done for the granular flow on an inclined plane (Patra et al., 2018b) or the 1955 Atenquique lahar (Bevilacqua et al., 2019), as well as even more sophisticated sampling methods (e.g., Bursik et al., 2012, Bayarri et al., 2015). As in the examples in this study, the hazard map and hazard density map and the associated statistics may be generated from such experiments. This may then be used to give additional information to a probabilistic assessment of the hazard rather than merely the inundation probability contours. In this context, the method outlined here is capable of generating the typical measures that would appear in routine statistical analyses of distributions. The most obvious among these is the mean flow boundary and standard deviation region. These are fundamentally newly calculable measures for PHMs and should be included in probabilistic hazard assessments.

As illustrated in the examples, the mean boundary may differ significantly from the median boundary owing to the asymmetric statistics that are inherent to typical nonlinear flows. Furthermore, this asymmetry suggests that the modal or maximum likelihood boundary may be of interest to users of these assessments as well. For the scientific component of probabilistic hazard assessments to have maximum impact and accuracy, these analyses should incorporate the full statistics of the PHM.

5 Conclusions

With new tools to calculate the moments of any given hazard edge distribution from a PHM, hazard map makers and analysts become able to estimate the likely location of the hazard edge and the uncertainty in that estimate. As numerical modeling capabilities continue to grow, the increased ease and speed of this type of analysis will allow the scientific community to quickly construct the PHM, PHDM, the mean hazard boundary curve, standard deviation region, and higher moments of the hazard boundary distribution, giving improved estimates of the true hazard zone. Although this theory is useful for many applications where probabilities are given spatially, further work is needed to constrain which types of hazard assessments and models are most suitable for this analysis and how best to blend this scientific information within existing hazard-map-making procedures. Similarly, additional work is required to study the probabilistic inverse problem further than the simple example given and to evaluate the use of this theory in a time-dependent context. Overall, the analysis put forth here significantly enhances the study of spatial probabilistic hazards, yielding new estimates of the likely hazard boundary and the uncertainty in those estimates.

Appendix A: PHDM normalization

Using the definition of the integral curves in the text as steepest ascent integral curves x_n(s) such that

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{d}{\mathit{x}}_n}{\mathrm{d}s}}=\widehat{\mathit{n}}:={\displaystyle \frac{\mathrm{\nabla }\mathit{\varphi }}{\left|\mathrm{\nabla }\mathit{\varphi }\right|}},\end{array}$$

with $\mathrm{0}\le s<L$ where ${\mathit{x}}_n\left(\mathrm{0}\right)\in \partial {\mathrm{\Omega }}^{\mathrm{0}}$ and ${\mathit{x}}_n\left(L\right)\in \partial {\mathrm{\Omega }}^{\mathrm{1}}$, we prove the following lemma.

Lemma: PHDM normalization. Consider a probabilistic hazard map (PHM) ϕ(x) with the following properties (as given in the text):

• i.

  ϕ(x) has compact support in a region Ω⁰⊂ℝ², with ϕ=0 on ∂Ω⁰.

• ii.

  $\mathrm{0}\le \mathit{\varphi }\le \mathrm{1}$.

• iii.

  There exists a (possibly unconnected) set Ω¹⊂Ω⁰ in which ϕ(x)=1.

• iv.

  $\mathit{\varphi }\left(\mathit{x}\right)\in C^{\mathrm{2}}({\mathbb{R}}^{\mathrm{2}}\backslash \mathit{\{}\partial {\mathrm{\Omega }}^{\mathrm{0}}\cup \partial {\mathrm{\Omega }}^{\mathrm{1}}\mathit{\left\}}\right)$.

• v.

  If ϕ has local maxima ϕ(x_m)<1 with m∈1 … M, then each x_m can be surrounded by a region of influence ℛ_m bounded by integral curves as defined in the text.

Define

$$\begin{array}{}\text{(A2)}& {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}:={\mathrm{\Omega }}^{\mathrm{0}}\backslash \bigcup _{m=\mathrm{1}}^M{\mathcal{R}}_{\mathrm{m}}.\end{array}$$

Then

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\displaystyle }{\int }_{\mathrm{0}}^L\left|\mathrm{\nabla }\mathit{\varphi }\left(\mathit{x}{}_n\left(s\right)\right)\right|\mathrm{d}s& {\displaystyle }=\mathrm{1}\Rightarrow \underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\left|\mathrm{\nabla }\mathit{\varphi }\right(\mathit{x}\left)\right|\mathrm{d}A\\ {\displaystyle }& {\displaystyle }=-\underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\mathit{\varphi }\mathrm{\nabla }\cdot \widehat{\mathit{n}}\mathrm{d}A.\end{array}\end{array}$$

Proof. Consider the level sets of ϕ(x) as integral curves everywhere tangent to the level sets, given by

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}{\mathit{x}}_t}{\mathrm{d}l}}=\mathbf{R}\widehat{\mathit{n}},\end{array}$$

where $\mathbf{R}=\left(\begin{array}{c}\mathrm{0}-\mathrm{1}\\ \mathrm{1}\mathrm{0}\end{array}\right)$ is a rotation matrix. The curves x_t(l) are the level sets of ϕ. By invoking the tangent curves and parameterizing each normal curve x_n(s;l), each PDF $\frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}s}(s;l)$ may be extended away from a single integral curve:

$$\begin{array}{}\text{(A5)}& \mathrm{d}l\underset{\mathrm{0}}{\overset{L}{\int }}\left|\mathrm{\nabla }\mathit{\varphi }\left({\mathit{x}}_n(s;l)\right)\right|\mathrm{d}s=\mathrm{d}l.\end{array}$$

Integration through l yields an integral of the gradient magnitude over ${\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}$ with the area element parameterized in local orthonormal curvilinear coordinates by (s, l) which are derived from the PHM instead of ordinary space. Since both sets of integral curves (x_n, x_t) are parameterized by arc length, then dA(x)=dA(s, l). Because of these facts, the value of the integration will depend on the ϕ particular to the problem under consideration. Note that, in the coordinates (s, l), each of the domains of influence related to local maxima less than unity (ℛ_m) is a rectilinear patch of (s, l) space (segment of l∈[0, ∞)) which is removed, leaving the integration domain of l which contains gaps:

$$\begin{array}{}\text{(A6)}& \mathcal{L}={\mathbb{R}}_{+}\backslash \bigcup _{m=\mathrm{1}}^M{\mathcal{R}}_{\mathrm{m}}\left(l\right).\end{array}$$

The right hand side of Eq. (A3) can then be evaluated by invoking a special case of the divergence theorem. Introducing an outward unit normal ν on $\partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}$ and noting that

$$\begin{array}{}\text{(A7)}& \left\{\begin{array}{l}\mathit{\varphi }=\mathrm{0}\mathrm{on}\partial {\mathrm{\Omega }}^{\mathrm{0}}\cap \partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}\\ \widehat{\mathit{n}}\cdot \mathit{\nu }=\mathrm{0}\mathrm{on}\mathrm{each}\mathrm{segment}\partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}\cap \partial {\mathcal{R}}_{\mathrm{m}}\\ \partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}=\left(\partial {\mathrm{\Omega }}^{\mathrm{0}}\cap \partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}\right)\cup \bigcup _{m=\mathrm{1}}^M\left(\partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}\cap \partial {\mathcal{R}}_{\mathrm{m}}\right),\end{array}\right.\end{array}$$

the divergence theorem gives

$$\begin{array}{}\text{(A8)}& \underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\mathrm{\nabla }\cdot \left(\mathit{\varphi }\widehat{\mathit{n}}\right)\mathrm{d}A=\underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\left|\mathrm{\nabla }\mathit{\varphi }\right|+\mathit{\varphi }\mathrm{\nabla }\cdot \widehat{\mathit{n}}\mathrm{d}A=\underset{\partial {\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\mathit{\varphi }\widehat{\mathit{n}}\cdot \mathit{\nu }\mathrm{d}l\equiv \mathrm{0}.\end{array}$$

From this result, integration of Eq. () will yield

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{\displaystyle }\underset{\mathcal{L}}{\int }{\int }_{\mathrm{0}}^L\left|\mathrm{\nabla }\mathit{\varphi }\left({\mathit{x}}_n(s;l)\right)\right|\mathrm{d}s\mathrm{d}l& {\displaystyle }=\underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\left|\mathrm{\nabla }\mathit{\varphi }\right|\mathrm{d}A=-\underset{{\stackrel{\mathrm{̃}}{\mathrm{\Omega }}}^{\mathrm{0}}}{\int }\mathit{\varphi }\mathrm{\nabla }\cdot \widehat{\mathit{n}}\mathrm{d}A\\ {\displaystyle }& {\displaystyle }=\mathrm{constant}=:Q,\end{array}\end{array}$$

which is finite and positive due to the smoothness assumed in defining the hazard map ϕ(x) and the fact that $\left|\mathrm{\nabla }\mathit{\varphi }\right|$ is everywhere nonnegative. In cases where the gradient of the PHM tends to infinity, the smoothness of the PHM must be modified in the sense of distributions, allowing derivatives and integrals on the PHM in the sense of distributions. This defines the PHDM normalization Q and proves the lemma.

Appendix B: Recipe for generation and analysis of the PHM

The following procedure is designed to simplify and condense the many steps involved in the above analysis. This recipe has been written as a blueprint for carrying out the analysis using a numerical model of the process in question. This procedure assumes that a numerical model of the process has already been selected.

1. Identify the collection of uncertain input parameters (β) that are to be varied for the experiment as well as their ranges, correlation structure, and joint distribution. This may be done by the methods of Bevilacqua et al. (2019) or other methods chosen for the particular application. This process generates the joint PDF f(β) measuring the input space.

2. Generate random sample input vectors β distributed according to f and run the numerical model on these inputs. For example, Latin hypercube sampling is a well-established procedure for defining pseudorandom designs of samples in ℝⁿ for multivariate uniform probability distributions (McKay et al., 1979; Stein, 1987; Owen, 1992). If all of the parameters are uncorrelated but nonuniform, then each random variable may be sampled independently according to its marginal distribution alone. If f measures a correlated set of input random variables, then a sampling scheme must be used which respects the correlation structure and marginal distributions of the random variables. This sampling is very straightforward for the multivariate normal distribution which has implementations in most scientific computing languages; however, other distributions require different techniques such as the NORTA (NORmal To Anything) process (Cairo and Nelson, 1997). The number of samples generated is directly related to the smoothness of the output PHM to be created.

3. Identify the threshold of interest for a given model observable and generate the indicator function for each model run according to the threshold. In practice, this can be realized as a sequence of grids, each of which is an element-wise indicator. In a MATLAB environment, this could be accomplished for each grid with the following syntax in which the indicator and the variable of interest are arrays of double-precision floating-point:
   
   INDICATOR = double(VARIABLE >  THRESHOLD_VALUE).

4. The PHM ϕ is obtained by performing an element-wise average of all of these indicator grids. The random sampling according to f(β) ensures that the element-wise average represents the underlying distribution. Some smoothing of ϕ may be required for further processing to reduce the effects of finite input sampling.

5. Generate the level set unit normal vector field for ϕ using a numerical gradient operation g≈∇ϕ which can be implemented in most scientific computing language. Each component of g is a grid of the same size as ϕ. Calculate the unit normal vector field as

   $$\begin{array}{}\text{(B1)}& \widehat{\mathit{n}}=\left\{\begin{array}{ll}\mathit{g}/\left|\mathit{g}\right|& \mathrm{for}\left|\mathit{g}\right|\ne \mathrm{0}\\ \mathrm{0}& \mathrm{elsewhere}\end{array}\right.\end{array}$$

   or in a MATLAB script:
   
   [g_x, g_y] = gradient (PHI, dx, dy)
   
   g_mag = (g_x.^2 +  g_y.^2).^0.5

n_x=g_x./g_mag

n_y=g_y./g_mag

n_x(isnan(n_x))=0

n_y(isnan(n_y))=0.

6. Generate the gradient-ascending integral curves x_n. This is done numerically by any gradient ascent algorithm for solving

   $$\begin{array}{}\text{(B2)}& {\mathit{x}}_n\left(s_{k+\mathrm{1}}\right)={\mathit{x}}_n\left(s_k\right)+\mathrm{\Delta }{s}_k\widehat{\mathit{n}}\left({\mathit{x}}_n\left(s_k\right)\right).\end{array}$$

   $\widehat{\mathit{n}}\left({\mathit{x}}_n\right(s_k\left)\right)$ may be computed for each step by interpolating the vector field $\widehat{\mathit{n}}$ at the point x_n(s_k). The initial position for each curve is given by the set of points output by a numerical contouring algorithm to find points very close to the ϕ=0 level set. A large number output points should be specified in the contouring. In practice, a very small level set value is chosen for this contouring ($\mathit{\varphi }={\mathrm{10}}^{-\mathrm{4}}$ in our example 2). In general, small values of Δs_k give smoother curves. Iteration follows for each curve until g(x_n(s_k))→0 with either (a) ϕ(x_n(s_k))=1 or (b) ϕ(x_n(s_k))<1. These two conditions ensure that either the integral curve ends on ∂Ω¹ (condition a) or a local maximum ϕ(x_m) (condition b).

7. Generate the valid regions of the PHM by eliminating all regions in which integral curves terminate at local maxima. In practice, this could be accomplished by generating an indicator grid in which all array elements within a certain neighborhood of these curves are set to zero and all other elements are set to 1.

8. Generate the PHDM using the gradient magnitude $\left|\mathit{g}\right|$ that was already calculated. For the normalizing constant Q, choose a numerical integration scheme (trapezoidal rule, Simpson's rule, Monte Carlo, or others) and compute the integral over the entire grid of the gradient magnitude filtered by the valid region indicator created in the last step. The PHDM is then given as $\mathit{\psi }=\left|\mathit{g}\right|/Q$. Any integration of ψ must be computed with the same method.

9. To compute the mean on each curve, treat the sequence ϕ_k=ϕ(x_n(s_k)) as samples of a CDF over the sequence s_k, numerically approximate the formula

   $$\begin{array}{}\text{(B3)}& \mathit{\mu }=\underset{\mathrm{0}}{\overset{L}{\int }}\mathrm{1}-\mathit{\varphi }\left(s\right)\mathrm{d}s\end{array}$$

   with a suitable technique (trapezoidal rule, Simpson's rule, or others), and calculate the mean position on each curve as x_μ=x_n(μ) by interpolation. The variance on each curve can be found similarly by numerical approximation of

   $$\begin{array}{}\text{(B4)}& {\mathit{\sigma }}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{L}{\int }}\mathrm{2}s(\mathrm{1}-\mathit{\varphi }(s\left)\right)\mathrm{d}s-{\mathit{\mu }}^{\mathrm{2}}.\end{array}$$

   Higher moments can be calculated by similar formulas involving the term 1−ϕ(s) which comes from integration by parts of the classical probability distribution moment formulas. The median position x_med on each curve is the point at which $\mathit{\varphi }=\mathrm{1}/\mathrm{2}$ and the mode or maximum likelihood position x_ℳ is at the maximum of ψ(x_n(s)), which can be found by interpolation of the PHDM along the computed curves. Each x_μ, x_med, and x_ℳ is then connected for all of the integral curves, making a set of curves for each of the mean, median, and mode, and a region bounded by two curves for the standard deviation.