Addressing challenges in uncertainty quantification: the case of geohazard assessments

We analyse some of the challenges in quantifying uncertainty when using geohazard models. Despite the availability of recently developed, sophisticated ways to parameterise models, a major remaining challenge is constraining the many model parameters involved. Additionally, there are challenges related to the credibility of predictions required in the assessments, the uncertainty of input quantities, and the conditional nature of the quantification, making it dependent on the choices and assumptions analysts make. Addressing these challenges calls for more insightful approaches yet to be developed. However, as discussed in this paper, clarifications and reinterpretations of some fundamental concepts and practical simplifications may be required first. The research thus aims to strengthen the foundation and practice of geohazard risk assessments.

1 Introduction

Uncertainty quantification (UQ) helps determine the uncertainty of a system's responses when some quantities and events in such a system are unknown. Using models, the system's responses can be calculated analytically, numerically, or by random sampling (including the Monte Carlo method, rejection sampling, Monte Carlo sampling using Markov chains, importance sampling, and subset simulation) (Metropolis and Ulam, 1949; Brown, 1956; Ulam, 1961; Hastings, 1970). Sampling methods are frequently used because of the high-dimensional nature of hazard events and associated quantities. Sampling methods result in less expensive and more tractable uncertainty quantification than analytical and numerical methods. In the sampling procedure, specified distributions of the input quantities and parameters are sampled, and respective outputs of the model are recorded. This process is repeated as many times as required to achieve the desired accuracy (Vanmarcke, 1984). Eventually, the distribution of the outputs can be used to calculate probability-based metrics, such as expectations or probabilities of critical events. Model-based uncertainty quantification using sampling is now more often used in geohazard assessments, e.g. Uzielli and Lacasse (2007), Wellmann and Regenauer-Lieb (2012), Rodríguez-Ochoa et al. (2015), Pakyuz-Charrier et al. (2018), Huang et al. (2021), Luo et al. (2021), and Sun et al. (2021a).

This paper considers recent advances in UQ and analyses some remaining challenges. For instance, we note that a major problem persists, namely constraining the many parameters involved. Only some parameters can be constrained in practice based solely on historical data (e.g. Albert et al., 2022). Another challenge is that model outputs are conditional on the choice of model parameters and the specified input quantities, including initial and boundary conditions. For example, a geological system model could be specified to include some geological boundary conditions (Juang et al., 2019). Such systems are usually time-dependent and spatial in nature and may involve, e.g. changing conditions (e.g. Chow et al., 2019). Incorporating uncertainties related to such conditions complicates the modelling and demands further data acquisition. Next, models could accurately reproduce data from past events but may be inadequate for unobserved outputs or predictions. This might be the case when predicting, e.g. extreme velocities in marine turbidity currents, which are driven by emerging and little-understood soil and fluid interactions (Vanneste et al., 2019). Overlooking these challenges implies that the quantification will only reflect some aspects of the uncertainty involved. These challenges are, unfortunately, neither exhaustively nor clearly discussed in the geohazard literature. Options and clarifications addressing these challenges are underreported in the field. Analysing these challenges can be useful in treating uncertainties consistently and providing meaningful results in an assessment. This paper's objective is to bridge the gap in the literature by providing an analysis and clarifications enabling a useful quantification of uncertainty.

It should be emphasised that, in this paper, we consider uncertainty quantification in terms of probabilities. Other approaches to measure or represent uncertainty have been studied by, for example, Zadeh (1968), Shafer (1976), Ferson and Ginzburg (1996), Helton and Oberkampf (2004), Dubois (2006), Aven (2010), Flage et al. (2013), Shortridge et al. (2017), Flage et al. (2018), and Gray et al. (2022a, b). These approaches will not be discussed here. The discussion about the complications in UQ related to computational issues generated by sampling procedures is also beyond the scope of the current work.

The remainder of the paper is as follows. In Sect. 2, based on recent advances, we describe how uncertainty quantification using geohazard models can be conducted. Next, some remaining challenges in UQ are identified and illustrated. Options to address the challenges in UQ are discussed in Sect. 3. A simplified example, further illustrating the discussion, is found in Sect. 4, while the final section provides some conclusions.

2 Quantifying uncertainty using geohazard models

In this section, we make explicit critical steps in uncertainty quantification (UQ). We describe a general approach to UQ that considers uncertainty as the analysts' incomplete knowledge about quantities or events. The UQ approach described is restricted to probabilistic analysis. Emphasis is made on the choices and assumptions usually made by analysts.

A geohazard model can be described as follows. We consider a system (e.g. debris flow) with a set of specified input quantities X (e.g. sediment concentration, entrainment rate) whose relationships to the model output Y (e.g. runout volume, velocity, or height of flow) can be expressed by a set of models 𝓜. Analysts identify or specify X, Y, and 𝓜. A vector Θ_m (including, e.g. friction, viscosity, turbulence coefficients) parameterises a model m in 𝓜. The parameters Θ_m determine specific functions among a family of potential functions modelling the system. Accordingly, a model m can be described as a multi-output function with, e.g. Y={runout volume, velocity, height of flow}. Based on Lu and Lermusiaux (2021), we can write

$$\begin{array}{}\text{(1)}& {\displaystyle }m:{\mathit{X}}_{s,t\times }{\mathbf{\Theta }}_m\to {\mathit{Y}}_{s,t}\text{(2)}& {\displaystyle }m\equiv ({\mathit{E}}_m,\mathit{S}{\mathit{G}}_m,\mathit{B}{\mathit{C}}_m,\mathit{I}{\mathit{C}}_m).\end{array}$$

Realisations of Y are the model responses y when elements in X take the values x at a spatial location s∈S and a specific time t∈T, and parameters θ_m∈Θ_m are used. In expression (1), $\mathit{X}\subset {\mathbb{R}}^{d_X}$ is the set of specified input quantities, $\mathit{T}\subset {\mathbb{R}}^{d_{\mathrm{T}}}$ is the time domain, $\mathit{S}\subset {\mathbb{R}}^{d_{\mathrm{S}}}$ is the spatial domain, ${\mathbf{\Theta }}_m\subset {\mathbb{R}}^{d_{{\mathrm{\Theta }}_m}}$ corresponds to a parameter vector, and $\mathit{Y}\subset {\mathbb{R}}^{d_Y}$ is the set of model outputs. To consider different dimensions, $d=\mathit{\{}\mathrm{1},\mathrm{2},\text{ or }\mathrm{3}\mathit{\}}$. The system is fully described if m is specified in terms of a set of equations E_m (e.g. conservation equations), the spatial domain geometry SG_m (e.g. extension, soil structure), the boundary conditions BC_m (e.g. downstream flow), and the initial conditions IC_m (e.g. flow at t=t₀); see Eq. (2).

Probabilities reflecting analysts' uncertainty about input quantities are specified in uncertainty quantification. Such distributions are then sampled many times, and the distribution of the produced outputs can be calculated. The output probability distribution for a model m can be denoted as $f\left(y\right|x,{\mathit{\theta }}_m,m)$, for realisations y, x, θ_m, m of Y, X, Θ_m, and 𝓜, respectively.

Betz (2017) has suggested that the parameter set is fully described by a parameter vector Θ; Eq. (3) is as follows:

$$\begin{array}{}\text{(3)}& \mathbf{\Theta }=\mathit{\{}{\mathbf{\Theta }}_m,{\mathbf{\Theta }}_X,{\mathbf{\Theta }}_{\mathit{\epsilon }},{\mathbf{\Theta }}_{\mathrm{o}}\mathit{\}},\end{array}$$

in which Θ_m refers to parameters of the model m, Θ_X are parameters linked to the input X, Θ_ε is the vector of the output-prediction error ε, and Θ_o is the vector associated with observation/measurement errors. More explicitly, to compute an overall joint probability distribution, we may have the following distributions:

• $f\left(y\right|x,{\mathit{\theta }}_m,m)$ is the distribution of Y when X takes the values x, and parameters θ_m∈Θ_m and a model m∈𝓜 are used to compute y;

• f(x|θ_X,m) is the conditional distribution of X given the parameters θ_X∈Θ_X and the model m. Note that each m defines which elements in X are to be considered in the analysis;

• $f\left(x\right|\widehat{x},{\mathit{\theta }}_{\mathrm{o}})$ is a distribution of X given the observed values $\widehat{\mathit{X}}=\widehat{x}$ and the observation/measurement error parameters θ_o∈Θ_o;

• additionally, one can consider $f\left(y^{*}\right|y,{\mathit{\theta }}_{\mathit{\epsilon }}$,m), which is a distribution of Y^*, the future system's response, conditioned on the model output y and the output-prediction error vector θ_ε∈Θ_ε. The output-prediction error ε is the mismatch between the model predictions and non-observed system responses y^*. ε is used to correct the imperfect model output y (Betz, 2017; Juang et al., 2019).

If, for example, the parameters Θ_m are poorly known, a prior distribution π(θ_m|m) weighing each parameter value θ_m for a model m is usually specified. A prior is a subjective probability distribution quantified by expert judgement representing uncertainty about the quantities prior to considering data (Raices-Cruz et al., 2022). When some measurements $\textcolor{red}{}=\mathit{\{}\widehat{\mathit{Y}}=\widehat{y},\widehat{\mathit{X}}=\widehat{x}\mathit{\}}$ are available, such parameter values θ_m, or their distributions π(θ_m|m), can be constrained by back-analysis methods. Note that measurements form part of different sources of data 𝓓, i.e. $\textcolor{red}{}\in \mathcal{D}$. Back-analysis methods include matching experimental measurements $\widehat{y}$ and calculated model outputs y using different assumed values ${\mathit{\theta }}_m^{\prime }$. Values for θ_m can be calculated as follows (based on Liu et al., 2022):

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}_m=\text{argmin}[\widehat{y}-y(\widehat{x},{\mathit{\theta }}_m^{\prime }\left)\right].\end{array}$$

The revision or updating of the prior π(θ_m|m) with measurements to obtain a posterior distribution denoted $\mathit{\pi }\left({\mathit{\theta }}_m\right|\textcolor{red}{},m)$ is also an option in back analysis. The updating can be calculated as follows (based on Juang et al., 2019; Liu et al., 2022):

$$\begin{array}{}\text{(5)}& \mathit{\pi }\left({\mathit{\theta }}_m|\textcolor{red}{},m\right)={\displaystyle \frac{\mathcal{L}\left({\mathit{\theta }}_m\right|\textcolor{red}{}\left)\mathit{\pi }\right({\mathit{\theta }}_m\mathrm{|}m)}{\int \mathcal{L}\left({\mathit{\theta }}_m\mathrm{|}\textcolor{red}{}\right)\mathit{\pi }\left({\mathit{\theta }}_m\mathrm{|}m\right)\mathrm{d}{\mathit{\theta }}_m}},\end{array}$$

where $\mathcal{L}\left({\mathit{\theta }}_m\right|\textcolor{red}{})=f(\textcolor{red}{}\left|{\mathit{\theta }}_m\right)$ is a likelihood function, i.e. a distribution that weighs given θ_m.

Similarly, we can constrain any of the distributions above, e.g. $f\left(y\right|x,{\mathit{\theta }}_m,m)$, or $f\left(x\right|{\mathit{\theta }}_X,m)$ to obtain $f\left(y\right|x,{\mathit{\theta }}_m,\textcolor{red}{},m)$ and $f\left(x\right|{\mathit{\theta }}_X,\textcolor{red}{},m)$, respectively.

For a geohazard problem, it is often possible to specify several competing models, e.g. distinct geological models with diverse boundary conditions; see expression (2). If the available knowledge is insufficient to determine the best model, different models m can be considered. The respective overall output probability distribution is computed as (Betz, 2017; Juang et al., 2019)

$$\begin{array}{}\text{(6)}& {\displaystyle }f\left(y\right|x,\mathbf{\Theta },\mathcal{D},\mathcal{M})=\sum f(y|x,\mathit{\theta },\textcolor{red}{},m)\mathit{\omega }\left(m\right|\mathcal{D},\mathcal{M})\text{(7)}& {\displaystyle }f\left(y\right|x,\mathit{\theta },\textcolor{red}{},m)=\int f(y|x,\mathit{\theta },m)\mathit{\pi }\left(\mathit{\theta }\right|\textcolor{red}{},m)\mathrm{d}\mathit{\theta }.\end{array}$$

In Eq. (6), $\mathit{\omega }\left(m\right|\mathcal{D},\mathcal{M})$ is a distribution weighing each model m in 𝓜.

The various models 𝓜, their inputs X, parameters Θ, outputs Y, and experimental data can be coupled all together through a Bayesian network, as has been suggested by Sankararaman and Mahadevan (2015) or Betz (2017). One possible configuration of a network coupling some elements in 𝓜, X, Θ, Y, and Y^* is illustrated in Fig. 1.

Figure 1 A configuration of a network coupling some elements in 𝓜, X, Θ, Y, and Y^*.

The previous description of a general approach to UQ considers uncertainty as that reflected in the analysts' incomplete knowledge about quantities or events. In UQ, to measure or describe uncertainty, subjective probabilities can be used and constrained using observations . It is also explicitly shown that model outputs are conditional on observations made available and models 𝓜 chosen by analysts. Analysts might also select several parameters Θ and initial and boundary conditions, BC_m and IC_m. Based on the above description, in the following, we analyse some of the challenges that arise when conducting UQ.

As mentioned, back-analysis methods help constrain some elements in Θ. However, given the considerable number of parameters (see expressions 1–3) and data scarcity, constraining Θ is often only achieved in a limited fashion. Back-analysis is further challenged by the potential dependency among Θ or 𝓜 and between Θ and SG_m, BC_m, and IC_m. We also note that back analysis, or, more specifically, inverse analysis, faces problems regarding non-identifiability, non-uniqueness, and instability. Non-identifiability occurs when some parameters do not drive changes in the inferred quantities. Non-uniqueness arises because more than one set of fitted or updated parameters may adequately reproduce observations. Instability in the solution arises from errors in observations and the non-linearity of models (Carrera and Neuman, 1986). Alternatively, in specifying a joint distribution f(x,θ) to be sampled, analysts may consider the use of e.g. Bayesian networks (Albert et al., 2022). However, under the usual circumstance of a lack of information, establishing such a joint distribution is challenging and requires that analysts encode many additional assumptions (e.g. prior distributions, likelihood functions, independence, linear relationships, normality, stationarity of the quantities and parameters considered); see e.g. Tang et al. (2020), Sun et al. (2021b), Albert et al. (2022), Pheulpin et al. (2022). A more conventional choice is that x or θ are specified using the maximum entropy principle (MEP), to specify the least biased distributions possible on the given information (Jaynes, 1957). Such distributions are subject to the system's physical constraints based on some available data. The information entropy of a probability distribution measures the amount of information contained in the distribution. The larger the entropy, the less information is provided by the distribution. Thus, by maximising the entropy over a suitable set of probability distributions, one finds the least informative distribution in the sense that it contains the least amount of information consistent with the system's constraints. Note that a distribution is sought over all the candidate distributions subject to a set of constraints. The MEP has been questioned since its validity and usefulness lie in the proper choice of physical constraints (Jaynes, 1957; Yano 2019). Doubts are also raised regarding the potential information loss when using the principle. Analysts usually strive to use all available knowledge and avoid unjustified information loss (Christakos, 1990; Flage et al., 2018).

Options to address the parametrisation challenge also include surrogate models, parameter reduction, and model learning (e.g. Lu and Lermusiaux, 2021; Sun et al., 2021b; Albert et al., 2022; Degen et al., 2022; Liu et al., 2022). Surrogate models are learnt to replace a complicated model with an inexpensive and fast approximation. Parameter reduction is achieved based on either principal component analysis or global sensitivity analysis to determine which parameters significantly impact model outputs and are essential to the analysis (Degen et al., 2022; Wagener et al., 2022). Remarkably, versions of the model learning option do not need any prior information about model equations E_m but require local verification of conservation laws in the data (Lu and Lermusiaux, 2021). These approaches still require large data sets sourced systematically, which is a frequent limitation in geohazard assessments. More importantly, however, is that, like many models, the credibility of unobserved surrogate model outputs can always be questioned, since, for instance, records may miss crucial events (Woo, 2019). Models may also fail to reproduce outputs caused by recorded abrupt changes (e.g. extreme velocities of turbidity currents) (Alley, 2004). An additional point is the issue of incomplete model response, which refers to a model not having a solution for some combinations of the specified input quantities (Cardenas, 2019; van den Eijnden et al., 2022).

In bypassing the described challenges when quantifying uncertainty, simplifications are usually enforced, sometimes unjustifiably, in the form of assumptions, denoted here by Ą. The set Ą can include one or more of the assumptions listed in Table 1. Note that the set of assumptions can be increased with those assumptions imposed by using specific models 𝓜 (e.g. conservation of energy, momentum, or mass, Mohr–Coulomb's failure criterion).

Table 1 Some enforced assumptions in UQ for geohazard assessments.

3 Addressing the challenges in uncertainty quantification

From the previous section, we saw that it is very difficult in geohazard assessments to meet data requirements for the ideal parameterisation of models. Further, we have noted that, although fully parameterised models could potentially be accurate at reproducing data from past events, these may turn out to be inadequate for unobserved outputs. We also made explicit that predictions are not only conditional on Θ but possibly also on SG_m, BC_m, and IC_m; see expressions (1)–(7). Ultimately, assumptions made also condition model outputs. More importantly, note that when only some model input quantities or parameters can be updated using data , the modelling will only reflect some aspects of the uncertainty involved. If the above challenges remain unaddressed, UQ lacks credibility. To address such challenges and provide increased credibility, clarifications and reinterpretation of some fundamental concepts and practical simplifications may be required, which are discussed in the following. Table 2 shows the major challenges found and how they are addressed in related literature, while in Table 3, some clarifications or considerations put forward by us are displayed. The discussion in this section builds on previous analysis by Aven and Pörn (1998), Apeland et al. (2002), Aven and Kvaløy (2002), Nilsen and Aven (2003), Aven and Zio (2013), Khorsandi and Aven (2017), and Aven (2019).

Table 2 Major challenges and options to address them in geohazard assessments.

Table 3 Some clarifications and considerations to address the challenges in UQ.

Among the clarifications, we consider a major conceptualisation suggested by the literature, which is the definition of uncertainty. Uncertainty refers to incomplete information or knowledge about a quantity or the occurrence of an event (Society for Risk Analysis, 2018). In Table 3, we denote this clarification as C1. Embracing this definition has some implications for uncertainty quantification using geohazard models. We use these implications to address the major complications and challenges. For instance, if uncertainty is measured in terms of probability, one such implication is that analysts are discouraged from using so-called frequentist probabilities. We note that frequentist probabilities do not measure uncertainty or lack of knowledge. Rather. such probabilities reflect frequency ratios representing fluctuation or variation in the outcomes of quantities. Frequentist probabilities are of limited use because these assume that quantities vary in large populations of identical settings, a condition which can be justified only for rather few geohazard quantities. The often one-off nature of many geohazard features and the impossibility of verifying or validating data by, e.g. a large number of repeated tests, make it difficult to develop such probabilities. Thus, a more meaningful and practical approach suggests to measure uncertainty by the use of knowledge-based (also referred to as judgemental or subjective) probabilities (Aven, 2019). A knowledge-based probability is an expression of the degree of belief in the occurrence of an event or quantity by a person assigning the probability conditional on the available knowledge 𝓚. Such knowledge 𝓚 includes not only data in the form of measurements made available but also other data sources in 𝓓. The models 𝓜 chosen for the prediction and the modelling assumptions Ą made by analysts are also part of 𝓚. Accordingly, to describe uncertainty about quantities, probabilities are assigned based on 𝓚, and, therefore, those probabilities are conditional on 𝓚. In the previous section, we have made evident the conditional nature of the uncertainty quantification (i.e. the probabilities) on measured data and models 𝓜 and wrote the expression f(y|x, Θ, 𝓓, 𝓜) for the overall output probability distribution (see Eq. 6). If assumptions Ą are also acknowledged as a conditional argument of the uncertainty quantification, we write more explicitly f(y|x, Θ, 𝓓, 𝓜, Ą) or equivalently $f\left(y\right|x,\mathbf{\Theta },\mathcal{K})$. We can therefore write

$$\begin{array}{}\text{(8)}& f\left(y\right|x,\mathbf{\Theta },\mathcal{K})=f(y|x,\mathbf{\Theta },\mathcal{D},\mathcal{M},\mathit{\text{Ą}}).\end{array}$$

The meaning of this expression is explained next. If, in a specific case, we would write $f\left(y\right|x,\mathbf{\Theta },\mathcal{K})=f(y|x,\mathit{\theta },\mathcal{D})$, it means that 𝓓 summarises all the knowledge that analysts have to calculate y given (realised or known) x and θ. Accordingly, the full expression in Eq. (8) implies that to calculate y, and given the knowledge of x and θ, the background knowledge includes 𝓓, 𝓜, and Ą. Note that 𝓚 can also be formed by observations, justifications, rationales, and arguments; thus, Eq. (8) can be further detailed to include these aspects of 𝓚. Structured methods exist to assign knowledge-based probabilities (see, e.g. Apeland et al., 2002; Aven, 2019). Here we should note, however, that since models form part of the available background knowledge 𝓚, models can also inform these knowledge-based probability assignments. It follows that, based on knowledge-based input probabilities, an overall output probability distribution calculated using models is also subjective or knowledge-based (Jaynes, 1957). Some of the implications of using knowledge-based probabilities are described throughout this section.

According to the left column in Table 2, the focus of the challenges relates to the model outputs, more specifically predictions (CH1 and CH2), input quantities (CH3–CH6), parameters (CH7–CH9), and models (CH10). We recall that uncertainty quantification helps determine the system's response uncertainty based on specified input quantities. Accordingly, an assessment focuses on the potential system's responses. The focus is often on uncertainty about future non-observed responses Y^*, which are approximated by the model output Y, considering some specified input quantities X. We recall that Y^* and X^* are quantities that are unknown at the time of the analysis but will take some value in the future and possibly become known. Thus, during an assessment, Y^* and X^* are the uncertain quantities of the system since we have incomplete knowledge about Y^* and X^*. Accordingly, the output-prediction error ε, the mismatch between the model prediction values y, and the non-observed system's response values y^* can only be specified based on the scrutiny of 𝓚.

There is another consequence of considering the definition of uncertainty put forward in C1, which links uncertainty solely to quantities or events. The consequence is that models, as such, are not to be linked to uncertainty. Models are merely mathematical artefacts. Models, per se, do not introduce uncertainty, but they are likely inaccurate. Accordingly, another major distinction is to be set in place. We recall that models, by definition, are simplifications, approximations of the system being analysed. They express or are part of the knowledge of the system. Models should therefore be solely used for understanding the performance of the system rather than for illusory perfect predictions. In Table 3, we denote the latter clarification as C2.

Regarding the challenges CH1 and CH2, we should note that geohazard analysts are often more interested in predictions rather than known system outputs. For instance, predictions are usually required to be calculated for input values not contained in the validation data. We consider that predictions are those model outputs not observed or recorded in the data, i.e. extrapolations out of the range of values covered by observations. Thus, the focus is on quantifying the uncertainty of the system's responses rather than on the accuracy of a model reproducing recorded data. This is the clarification C3 in Table 3. Considering this, models are yet to provide accuracy in reproducing observed outputs but, more importantly, afford credibility in predictions. Such credibility is to be assessed mainly in terms of judgements, since conventional validation cannot be conducted using non-observed outputs. Recall that model accuracy usually relates to comparing model outputs with experimental measurements (Roy and Oberkampf, 2011; Aven and Zio, 2013) and is the basis for validating models. Regarding the credibility of predictions, Wagener et al. (2022) have reported that such credibility can be mainly judged in terms of the physical consistency of the predictions. Such consistency is judged by checks rejecting physically impossible representations of the system. The credibility of predictions may also include the verification of the ability of models to accurately reproduce disruptive changes recorded in the data (Alley, 2004). However, as we have made explicit in the previous section, model predictions are conditional on a considerable number of critical assumptions and choices made by analysts (see Table 1 and clarification C4 in Table 3). Therefore, predictions can only be as good as the quality of the assumptions made. The assumptions could be wrong, and the examination of the impact of these deviations on the predictions must be assessed. To provide credibility of predictions, such assumptions and choices should be justified and scrutinised; see option O5 in Table 2. Option O5 addresses the challenge CH1; however, when conducting UQ, O5 has a major role when investigating input uncertainty, which is discussed next.

A critical task in UQ is the quantification of input uncertainty. Input uncertainty may originate when crucial historical events or disruptive changes are missing in the records (CH3). Some critical input quantities may also remain unidentified to analysts during an assessment (CH4). Analysts can unintendedly fail to identify relevant elements in X^* due to insufficiencies in data or limitations of existing models. For example, during many assessments, trigger factors that could bring a soil mass to failure could remain unknown to analysts (e.g. Hunt et al., 2013; Clare et al., 2016; Leynaud et al., 2017; Casalbore et al., 2020). UQ requires simulating sampled values from X, and elements in X can be mutually dependent. However, the joint distribution of X, namely f(x), is often also unknown. This is the challenge CH6. Considering the potential challenges CH3 to CH6, to specify f(x), we cannot solely rely on using the maximum entropy principle (MEP). The MEP may fail to advance an exhaustive uncertainty quantification in the input, e.g. by missing relevant values not recorded in the measured data. This would undermine the quality of predictions and, therefore, uncertainty quantification. Recall that the MEP suggests using the least informative distribution among candidate distributions constrained solely on measurements. Using counterfactual analysis, as described in Table 2, is an option. However, the counterfactual analysis will also fail to provide quality predictions, since this analysis focuses on counterfactuals (alternative events to observed facts , $\textcolor{red}{}\in \mathcal{D}$) and little on the overall knowledge available 𝓚. Note that the knowledge 𝓚 about the system includes, e.g. the assumptions made in the UQ, such as those shown in Table 1. Further note that such assumptions relate not only to data but also to input quantities, modelling, and predictions. Thus, it appears that the examination of these assumptions should be at the core of UQ in geohazard assessments, as suggested in Table 2, option O5. The risk assessment of deviations from assumptions was originally suggested by Aven (2013) and exemplified by Khorsandi and Aven (2017). An assumption deviation risk assessment evaluates different deviations, their associated probabilities of occurrence, and the effect of the deviations. A major distinctive feature of the assumption deviation risk assessment approach is the evaluation of the credibility of the knowledge 𝓚 supporting the assumptions made. Another feature of this approach is questioning the justifications supporting the potential for deviations. The examination of 𝓚 can be achieved by assessing the justifications for the assumptions made, the amount and relevance of data or information, the degree of agreement among experts, and the extent to which the phenomena involved are understood and can be modelled accurately. Justifications might be in the form of direct evidence becoming available, indirect evidence from other observable quantities, supported by modelling results, or possibly inferred by assessments of deviations of assumptions. This approach is succinctly demonstrated in the following section. Accordingly, we suggest specifying f(x) in terms of knowledge-based probabilities in conjunction with investigating input uncertainty using the assumptions deviation approach. This is identified as consideration C5 in Table 3.

Another point to consider is that when uncertainty is measured in terms of knowledge-based probabilities, analysts should be aware of what conditionality means. If, for example, a quantity X₂ is conditional on a quantity X₁, this implies that increased knowledge about X₁ will change the uncertainty about X₂. The expression that denotes this is conventionally written as X₂|X₁. Analysts may exploit this interpretation when specifying, e.g. the joint distribution f(x,θ). For example, when increased knowledge about a quantity X₁ will not result in increased knowledge about another quantity X₂, analysts may simplify the analysis according to the scrutiny of 𝓚, meaning that a distribution $f\left(y\right|x_{\mathrm{1},}{x}_{\mathrm{2}})$ to be specified may reduce to $f\left(y\right|x_{\mathrm{1}}\left)f\right(y\left|x_{\mathrm{2}}\right)$ according to probability theory. Apeland et al. (2002) have illustrated how conditionality in the setting of knowledge-based probabilities can inform the specification of a joint distribution.

The parameterisation problem, which involves the challenges CH7 to CH9 in Table 2, warrants exhaustive consideration. Addressing these challenges also requires some reinterpretation. To start, note that parameters are coefficients determining specific functions among a family of potential functions modelling the system. Those parameters constrain a model's output. Recall that y, as realisations of Y, are the model output when X takes the values x, and some parameters θ∈Θ, and models m∈𝓜 are used. Thus, as shown in the previous section, any output y is conditional on θ, and so is the uncertainty attached to y^*. We may also distinguish two types of parameters. We may have parameters associated with a property of the system. Other parameters exist that are merely artefacts in the models and are not properties of the system. As suggested, if uncertainty can solely be attached to events or quantities, we may say that parameters that are not properties of the system are not to be linked to any uncertainty. This is identified as clarification C6 in Table 3. For example, analysts may consider that the parameters not being part of the system as such are those linked to the output-prediction error ε, the vector associated with observation/measurement errors Θ_o, and the overall attached hyperparameters linked to probability distributions (including priors, likelihood functions). Analysts may consider the latter parameters as modelling artefacts, so it is questionable to attach uncertainty to them. Thus, focused on the uncertainty of the system responses rather than model inaccuracies, uncertainty is to be assigned to those parameters that represent physical quantities. Fixed single values can be assigned to those parameters that are not properties of the system. To help identify those parameters to which some uncertainty can be linked, we can scrutinise, e.g. the physical nature of these. In fixing parameters to a single value, we can still make use of back-analysis procedures, as mentioned previously. Analysts may have some additional basis to specify parameter values when the background knowledge available 𝓚 is scrutinised. 𝓚 can be examined to verify that not only data measurements but other sources of data, models, and assumptions made strongly support a specific parameter value. Based on this interpretation, setting the values of the parameters that are not properties of the system to a single value reduces the complications in quantifying uncertainty considerably. It also follows that analysts are encouraged to make explicit that model outputs are conditional on these fixed parameters, and on the model or models chosen, as we have shown in the previous section. The latter also leads us to argue that the focus of UQ is on the uncertainty of the system response rather than the inaccuracies of the models. This implies in a practical sense that in geohazard assessments, when parameters are clearly differentiated from specified input quantities, and models providing the most credible predictions are chosen, uncertainty quantification can then proceed. This parsimonious modelling approach is identified as consideration C7 in Table 3. This latter consideration addresses, to an extent, the challenge CH10.

In the following section, we further illustrate the above discussion by analysing a documented case in which UQ in a geohazard assessment was informed by modelling using sampling procedures.

4 Case analysis

To further describe the proposed considerations, we analyse a case reported in the specialised literature. The case deals with the quantification of uncertainty of geological structures, namely uncertainty about the subsurface stratigraphic configuration. Conditions in the subsurface are highly variable, whereas site investigations only provide sparse measurements. Consequently, subsurface models are usually inaccurate. At a given location, subsurface conditions are unknown until accurately measured. Soil investigation at all locations is usually impractical and uneconomical, and point-to-point condition variation cannot be known (Vanmarcke, 1984). Such uncertainty means significant engineering and environmental risk to, e.g. infrastructure built on the surface. One way to quantify this uncertainty is by calculating the probability of every possible configuration of the geological structures (Tacher et al., 2006; Thiele et al., 2016; Pakyuz-Charrier et al., 2018). Sampling procedures for UQ are helpful in this undertaking. We use an analysis and information from Zhao et al. (2021), which refer to a site located in the Central Business District, Perth, Western Australia, where six boreholes were executed. The case has been selected taking into account its simplicity to illustrate the points of this paper, but at the same time, it provides details to allow some discussion. Figure 2 displays the system being analysed.

Figure 2 Borehole logs in colours and longitudinal section reported by Zhao et al. (2021) located in the Central Business District, Perth, Western Australia. The records correspond to information on six boreholes. Three types of materials are revealed by the boreholes, including sand (yellow), clay (magenta), and gravel (blue).

In the system under consideration, a particular material type to be found in a non-bored point, a portion of terrain not penetrated during soil investigation, is unknown and thus uncertain. The goal is to compute the probability of encountering a given type of soil at these points. Zhao et al. (2021) focus on calculating the probabilities of encountering clay in the subsurface. The approach advocated was a sampling procedure to generate many plausible configurations of the geological structures and evaluate their probabilities. In a non-penetrated point in the ground, to calculate the probability of encountering a given type of soil c, p(y=c), Zhao et al. (2021) used a function that depends on two correlation parameters, namely the horizontal and vertical scale of fluctuation θ_h and θ_v. Note that spatial processes and their properties are conventionally assumed as spatially correlated. Such spatial variation may presumably be characterised by correlation functions, which depend on a scale of fluctuation parameter. The scale of fluctuation measures the distance within which points are significantly correlated (Vanmarcke, 1984). Equation (9) describes the basic components of the model chosen by Zhao et al. (2021) (specific details are given in the Appendix to this paper) as follows:

$$\begin{array}{}\text{(9)}& m:{\mathit{X}}_{s^{\times }}{\mathbf{\Theta }}_m\to {\mathit{Y}}_s\to p(\mathit{y}=c),\end{array}$$

where X is the collection of all specified quantities at borehole points s_x which can take values x from the set {sand, clay, gravel}, according to the setting in Fig. 2. Y is the collection of all model outputs with values y at non-borehole points s_y. Probabilities p(y=c) are computed based on the sampling of the values y and x, and a chosen model using the parameters θ_h=11.1 and θ_v=4.1 m, θ_h, θ_v∈Θ_m. Using the maximum likelihood method, the parameters were determined based on the borehole data revealed at the site. In determining parameters, the sampling from uniform and mutually independent distributions of θ_h and θ_v was the procedure advocated. The system is further described by a set of equations E_m (a correlation function and a probability function), the spatial domain geometry sg_m (a terrain block of 30×80 m), and the boundary conditions bc_m (the conditions at the borders). More details are given in the Appendix to this paper. Since this system is not considered time-dependent, the initial conditions IC_m were not specified.

Figure 3 Zhao et al. (2021) findings shown in their Fig. 9. (a) Most probable stratigraphic configuration. (b) Spatial distribution of the probability of the existence of clay. Reprinted from Engineering Geology, 238, Zhao et al. (2021), Probabilistic characterisation of subsurface stratigraphic configuration with modified random field approach, p. 106138, © 2021, with permission from the COPYRIGHT OWNER: Elsevier. Distances in metres.

The summary results reported by Zhao et al. (2021) are shown in Fig. 3. In Fig. 3, the most probable stratigraphic configuration, along with the spatial distribution of the probability of the existence of clay, is displayed. The authors focused on this sensitive material, which likely represents a risk to the infrastructure built on the surface.

Zhao et al. (2021) stated that “characterisation results of the stratigraphic configuration and its uncertainty are consistent with the intuition and the state of knowledge on site characterisation”. Next, throughout Zhao et al.'s (2021) analysis, the following assumptions were enforced (Table 4), although these were not explicitly disclosed by the authors.

Table 4 Assumptions enforced by Zhao et al. (2021).

Unfortunately, the authors did not report enough details on how the majority of these assumptions are justified. We should note, however, that providing these justifications was not the objective of their research. Yet, here we analyse how assumptions can be justified by scrutinising 𝓚 and using some elements of the assumption deviation approach described in the previous section. Table 5 summarises the analysis conducted and only reflects the most relevant observations and reservations we identified. Accordingly, the information in Table 5 may not be exhaustive but is still useful for the desired illustration. Table 5 displays some of our observations related to the credibility of the knowledge 𝓚. The examination of 𝓚 is achieved by assessing the amount and relevance of data or information, the extent to which the phenomena involved are understood and can be modelled accurately, the degree of agreement among experts, and the justifications for the assumptions made. Observations regarding the justifications for potential deviations from assumptions also form part of the analysis.

Not surprisingly, the observations in our analysis concentrate on the predictions' credibility. Recall that UQ focuses on the system's response, approximated by model predictions (considerations C2 and C3 in Table 3). For example, although using correlations is an accepted practice and a practical simplification, correlation functions appear counterintuitive to model geological structures or domains. Further, correlation functions do not help much in understanding the system (consideration C2 in Table 3). Recall that such structures are mainly disjoint domains linked to a finite set of possible categorical quantities (masses of soil or rock) rather than continuous quantities. Next, the variation of such structures can occur by abrupt changes in materials; thus, the use of smoothed correlation functions to represent them requires additional consideration. Moreover, the physical basis of the correlation functions is not clear, and physical models based on deposition processes may be suggested (e.g. Catuneanu et al., 2009). We should note a potential justification for the deviation from the assumption regarding the credibility of predictions. This is because knowledge from additional sources such as surface geology, sedimentology, local geomorphic setting, and structural geology was not explicitly taken into account in quantifying uncertainty. The revision of this knowledge can contribute to reducing the probability of deviation in predictions. Based on the observations in Table 5, we can conclude that there is potential to improve the credibility of predictions.

The choices made by Zhao et al. (2021) regarding the use of parameters with fixed values together with the choice for a single best model can be highlighted. These choices illustrate the points raised in considerations C6 and C7 (Table 3). The maximum likelihood method supported these choices; a back-analysis method focused on matching measurements and calculated model outputs using different assumed values for θ_h and θ_v. We highlight that a model judged to be the best model was chosen. This includes the specification of a particular spatial domain geometry in SG_m. Investigating the impact of the variation of SG_m was considered unnecessary. There was no need to specify several competing models, which is in line with our consideration labelled as C7 in this paper.

Zhao et al. (2021) investigated the joint distribution f(x), which was sampled to calculate probabilities. However, someone can suggest that the joint distribution f(x,θ,sg_m,bc_m) could have been produced. Nevertheless, we can argue that establishing such a joint distribution is challenging and requires, in many instances, that analysts encode many additional assumptions (e.g. prior distributions, likelihood functions, independence, linear relationships, normality, stationarity of the quantities and parameters considered).

Table 5 Examination of supporting knowledge 𝓚 and justifications for the potential deviation of assumptions.

A more crucial observation derived from the analysis of potential deviations of assumptions might considerably impact the credibility of predictions. This observation comes from revisiting the knowledge sources of Zhao et al.'s (2021) analysis, available from https://australiangeomechanics.org/downloads/ (last access: 29 June 2022). Another type of sensitive material was revealed by other soundings in the area, more specifically, silt. Depending on the revision of 𝓚, this fourth suspected material could be analysed in an extended uncertainty quantification of the system. Note that the specified input quantities X were originally assumed to take values x from the set {sand, clay, gravel}. Such an assumption was based on the records of six boreholes which were believed to be accurate. The latter illustrates the relevance of consideration C5 in Table 3.

Another choice by Zhao et al. (2021) is that they disregarded the possibility of incorporating measurement errors of the borehole data into the UQ, probably because these data were judged to be accurate. We recall in this respect that these errors reflect the inaccuracy of the measurements rather than the uncertainty about the system. As stated for consideration C6 (Table 3), we can hardly justify attaching uncertainty to measurement error parameters, since measurement errors are not a property of the system. The same can be said for the parameters θ_h and θ_v, which are not properties of the system. Note that their physical basis is questioned. We should note, however, that assuming global coefficients for the parameters θ_h and θ_v is an established practice (Vanmarcke, 1984; Lloret-Cabot et al., 2014; Juang et al., 2019). It can be pointed out that uncertainty quantification in this kind of system is, to an extent, sensitive to the choice of scale of fluctuation values (Vanmarcke, 1984). It can also be argued that using a global rather than local correlation between spatial quantities can misrepresent geological structure variation. Accordingly, further examination of the existing knowledge 𝓚 justifies some assessment of the impact of assuming a local rather than global scale of fluctuation.

Overall, the Zhao et al. (2021) analysis is, to an extent, based on the previously suggested definition of uncertainty; see the consideration C1 in Table 3.

We should stress that Zhao et al.'s (2021) uncertainty quantification refers specifically to the ground model described at the beginning of this section. In other words, the probabilities displayed in Fig. 3b are conditional on the parameters chosen (θ_h=11.1 and θ_v=4.1 m), the model selected (described by Eqs. 9, A1 and A2 in the Appendix to this paper), the specified spatial domain geometry sg_m (a terrain block of 30×80 m), and ultimately the assumptions made (listed in Table 4). This information is to be reported explicitly to the users of the results. This reflects the clarification C4 in Table 3.

Regarding the consideration of subjective probabilities, there has been some agreement on their use in this kind of UQ since Vanmarcke (1984). However, the use of knowledge-based probabilities in the extension described here is recommended, given the illustrated implications to advance UQ (as discussed in the previous section and stated in consideration C5). For example, increased examination of 𝓚 might have resulted in using a more informative distribution f(θ_h,θ_v) than the uniform distribution. The increased examination of 𝓚 might have led to different values for θ_h and θ_v, and a different model. Recall that the selection of the model and determination of parameters were based on the maximum likelihood method, which only uses measured data .

In our analysis of Zhao's et al. (2021) assessment, the examination of supporting knowledge 𝓚 resulted essentially in

1. judging the credibility of predictions;

2. providing justifications for assessing assumption deviations by considering the modelling of a fourth material;

3. considering additional data other than the borehole records, such as surface geology, sedimentology, local geomorphic setting, and structural geology;

4. analysing the possibility of distinct geological models with diverse spatial domain geometry and local correlations; and

5. ultimately, further examining the existing 𝓚.

5 Conclusions

In this paper, we have discussed challenges in uncertainty quantification (UQ) for geohazard assessments. Beyond the parameterisation problem, the challenges include assessing the quality of predictions required in the assessments, quantifying uncertainty in the input quantities, and considering the impact of choices and assumptions made by analysts. Such challenges arise from the commonplace situation of limited data and the one-off nature of geohazard features. If these challenges are kept unaddressed, UQ lacks credibility. Here, we have formulated seven considerations that may contribute to providing increased credibility in the quantifications. For example, we proposed understanding uncertainty as lack of knowledge, a condition that can only be attributed to quantities or events. Another consideration is that the focus of the quantification should be more on the uncertainty of the system response rather than the accuracy of the models used in the quantification. We drew attention to the clarification that models, in geohazard assessments, are simplifications used for predictions approximating the system's responses. We have also considered that since uncertainty is only to be linked to the properties of the system, models do not introduce uncertainty. Inaccurate models can, however, produce poor predictions and such models should be rejected. Then, an increased examination of background knowledge will be required to quantify uncertainty credibly. We also put forward that there could not be uncertainty about those elements in the parameter set that are not properties of the system. The latter also has pragmatic implications, including how the many parameters in a geohazard system could be constrained in a geohazard assessment.

We went into detail to show that predictions, and in turn UQ, are conditional on the model(s) chosen together with the assumptions made by analysts. We identified limitations of measured data to support the assessment of the quality of predictions. Accordingly, we have proposed that the quality of UQ needs to be judged based also on some additional crucial tasks. Such tasks include the exhaustive scrutiny of the knowledge coupled with the assessment of deviations of those assumptions made in the analysis.

Key to enacting the proposed clarifications and simplifications is the full consideration of knowledge-based probability. Considering this type of probability will help overcome the identified limitations of the maximum entropy principle or counterfactual analysis to quantify uncertainty in input quantities. We have exposed that the latter approaches are prone to produce unexhausted uncertainty quantification due to their reliance on measured data, which can miss crucial events or overlook relevant input quantities.

Appendix A

In this Appendix, the necessary details of the original analysis made by Zhao et al. (2021) are given. The following are the basic equations E_m used by these authors:

$$\begin{array}{}\text{(A1)}& {\displaystyle }p\left(\mathit{y}=c\right)\sim {\displaystyle \frac{{\sum }_{x_s\times \mathit{Y}}{\mathit{\rho }}_{x=c,y=c}}{{\sum }_{c=\mathrm{1}}^C\sum _{x_s\times \mathit{Y}}{\mathit{\rho }}_{x=c,y=c}}}\text{(A2)}& {\displaystyle }{\mathit{\rho }}_{xy}=\exp \left(-\mathit{\pi }{\displaystyle \frac{\stackrel{\mathrm{‾}}{s_x{s}_y}}{{\mathit{\theta }}_{\mathrm{h}}}}-\mathit{\pi }{\displaystyle \frac{\mathrm{|}{s}_x{s}_y}{{\mathit{\theta }}_{\mathrm{v}}}}\right),\end{array}$$

where X is the collection of all specified quantities at borehole points, which take values x. Y is the collection of all outputs at non-borehole points with values y. ρ_xy is the value of correlation between a quantity value x at a penetrated point s_x∈S_x and the value y at a non-penetrated point s_y∈S_y. $\stackrel{\mathrm{‾}}{s_x{s}_y}$ is the horizontal distance between points s_x and s_y, while |s_xs_y is the vertical one. θ_h and θ_v are the horizontal and vertical scales of fluctuation, respectively. Each material class considered is associated exclusively with an element in the set of integers $\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },C\mathit{\}}$. p(y=c) is the probability of encountering a type of material c in a point s_y. Such probability is initially approximated using Eq. (A1). More accurate probabilities are computed based on the repeated sampling of the joint distribution f(x,y), which was approximated using Eq. (A1). Equation (A1), described in short, approximates probabilities as the ratio of the sum of correlation values, calculated for a penetrated point in the set S_x and the set of non-penetrated points S_y for a given material c, to the sum of correlation values for all points and all materials.

Based on data collected at borehole locations, the selection of the type of correlation function and the scales of fluctuation took place using the maximum likelihood method. The authors considered three types of correlation functions, namely squared exponential, single exponential, and second-order Markov. In this case, the likelihood function $\mathcal{L}\left({\mathit{\theta }}_m\right|\textcolor{red}{})=f(\textcolor{red}{}\left|{\mathit{\theta }}_m\right)$ represents the likelihood of observing at borehole locations, given the spatial correlation structure θ_m. The squared exponential function yielded the maximum likelihood when the horizontal and vertical scales of fluctuation were set to 11.1 and 4.1 m, respectively. Hence, the squared exponential function correlation, whose expression is Eq. (A2) in this Appendix, was selected. Equations (A3) and (A4) correspond to the single exponential and the second-order Markov functions, respectively.

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\mathit{\rho }}_{xy}=\exp \left(-\mathrm{2}{\displaystyle \frac{\stackrel{\mathrm{‾}}{s_x{s}_y}}{{\mathit{\theta }}_{\mathrm{h}}}}-\mathrm{2}{\displaystyle \frac{\mathrm{|}{s}_x{s}_y}{{\mathit{\theta }}_{\mathrm{v}}}}\right)\text{(A4)}& {\displaystyle }\begin{array}{rl}{\mathit{\rho }}_{xy}& =\left(\mathrm{1}+\mathrm{4}{\displaystyle \frac{\stackrel{\mathrm{‾}}{s_x{s}_y}}{{\mathit{\theta }}_{\mathrm{h}}}}\right)\left(\mathrm{1}+\mathrm{4}{\displaystyle \frac{\mathrm{|}{s}_x{s}_y}{{\mathit{\theta }}_{\mathrm{v}}}}\right)\\ & \cdot \exp \left(-\mathrm{4}{\displaystyle \frac{\stackrel{\mathrm{‾}}{s_x{s}_y}}{{\mathit{\theta }}_{\mathrm{h}}}}-\mathrm{4}{\displaystyle \frac{\mathrm{|}{s}_x{s}_y}{{\mathit{\theta }}_{\mathrm{v}}}}\right)\end{array}\end{array}$$