# Bayesian Emulation of Geotechnical Deterioration Curves Using Quadratic and
B-Spline Hierarchical Models

Jordan L. Oakley and Aleksandra Svalova
School of Mathematics The authors are grateful for the financial support of
the Engineering and Physical Sciences Research Council (EPSRC) through the
programme Grant ACHILLES (EP/R034575/1).  Statistics and Physics  Newcastle
University  Newcastle upon Tyne  NE1 7RU  UK
Peter Helm Dr. Helm’s contribution to the manuscript preparation was supported
by the Research Hub for Decarbonised Adaptable and Resilient Transport
Infrastructures (DARe) funded by the UK Department for Transport and EPSRC
(EP/Y024257/1).  Mohamed Rouainia  and Stephanie Glendinning
School of Engineering  Newcastle University  Newcastle upon Tyne  NE1 7RU  UK
Dennis Prangle
School of Mathematics  University of Bristol  Bristol  BS8 1TW  UK
and Darren Wilkinson
Department of Mathematical Sciences  Durham University  Durham  DH1 3LE  UK

###### Abstract

The stability of geotechnical infrastructure assets, such as cuttings and
embankments, is crucial to the safe and efficient delivery of transport
services. The successful emulation of geotechnical models of deterioration of
infrastructure slopes has the potential to inform slope design, maintenance
and remediation by introducing the time dependency of deterioration into
geotechnical asset management. We have performed computer experiments of
deterioration, measured by the factor of safety (FoS), for a set of cutting
slope geometries and soil properties that are common in the southern UK.
Whilst computer experiments are an extremely useful and cost-effective method
of better understanding deterioration mechanisms, it would not be practical to
run enough experiments to understand relations between high-dimensional inputs
and outputs. Therefore, we trained a fully-Bayesian Gaussian process emulator
using an ensemble of 75 computer experiments to predict the FoS. We construct
two different emulator models, one approximating the FoS temporal evolution
with a quadratic model and one approximating the temporal evolution with a
B-spline model; and we emulated their parameters. We also compare the ability
of our models to predict failure time. The developed models could be used to
inform infrastructure cutting slope design and management, and extend
serviceable life.

##  1 Introduction

Engineered earthworks are a key part of national infrastructure in the United
Kingdom. Modern highways earthworks are designed so that their actual
stability exceeds the required minimum (typically by reducing the design slope
angle). The ratio between excess capacity and the minimum required to maintain
stability is known as factor of safety (FoS) against ultimate limit state
failure (i.e., states associated with collapse, see for example Burland et al.
(2023) and BSI (2023)). Factor of safety can be defined in a number of ways;
here the specific definition of FoS relating to slope stability is as
described by Duncan (1996), whereby FoS is the ratio of available shear
strength (the ability of a material to resist the component that acts parallel
to a plane) of the soil to the shear strength required to barely maintain
equilibrium/stability, where the slope is stable if FoS is greater than 1.
As the strength of slopes deteriorates due to short duration weather events,
seasonal cycles, longer term climate change and other cyclic loading (for
example, vehicular traffic), the FoS will reduce. This reduction in FoS due to
deterioration poses an additional problem for UK railway earthworks,
significant proportions of which are greater than 150 years old (Abbott,
2018). These earthworks tend to be overly steep and hence have a lower FoS at
construction (Perry et al., 2003) compared to the more modern highways
earthworks, and together with their longer exposure to processes that cause
deterioration (Briggs et al., 2023) mean they may be closer to ULS failure
than more recently constructed engineered slopes.
In addition to the ultimate limit state (ULS) as mentioned above, there are
additional limit state constraints that infrastructure assets tend to operate
within, these are known as serviceability limit states (SLS). While ULS is
concerned with the final collapse/failure of an asset and has significant
implications for safety if it occurs, SLS defines the minimum performance
required by an asset to satisfactorily fulfil its intended function. In the
context of slope stability, this could relate to deformations driven by
gradual deterioration. These deformations can affect rail track alignment,
where they begin to affect ride quality and in more extreme cases, if they
exceed a specified (SLS) threshold, would necessitate the imposition of
reduced speed limits for safety until remediation works could be undertaken.
These same drivers of deterioration can also, over long periods, lead to ULS
failure. An early identification of the assets likely to be affected by these
serviceability related issues in turn reduces maintenance, delay and
intervention costs and can reduce the risk of ULS failures occurring.
While the prediction of the time at which ULS may occur has been addressed by
Svalova et al. (2021) who modelled time to failure of transport infrastructure
cuttings as a function of geometry and soil strength, that work did not
address the rate of asset deterioration that can lead to exceedance of SLS
prior to ULS failure.
Prior geotechnical modelling work (Postill et al., 2021) has shown that the
change in FoS with time is not linear and so explicitly modelling the FoS
change as opposed to the time to failure allows an estimate to be made of the
rate of asset deterioration as a function of material strength and slope
geometry. Furthermore the magnitude of relative deterioration that has
occurred when interventions are undertaken has been shown to influence their
effectiveness (Armstrong et al., 2024). This shows that a model of changing
FoS with time prior to failure helps identify and prioritise the selection of
assets for monitoring and remediation and allows more effective management by
giving a clearer picture of the proportion of the asset portfolio which may be
affected at a given time.
As such, in this work, we address the deterioration towards shear failure of
cut slopes, where shear failure can be defined as the collapse of the slope by
sliding along a shear surface in the soil mass. This is performed using
computer experiments of deterioration and Gaussian process emulation (GPE).
Computer experiments can be used to simulate FoS evolution over time for a
wide range of geotechnical properties (see for example, Helm et al., 2024a)
and weather scenarios (Rouainia et al., 2020; Postill et al., 2021). The GPE-
computer experiment symbiosis has been a frequent occurrence in the
statistical machine learning literature over the last thirty years (see
Santner, 2018; Gramacy, 2020, and references therein).
A key body of related Bayesian literature has been developed by O’Hagan and
colleagues, who define the “emulator” to be a Gaussian process conditioned on
observations (i.e., training data from a computer simulator) that are assumed
to be normally distributed and whose parameters are inferred using a Bayesian
framework (O’Hagan, 2006; Bastos and O’Hagan, 2009). This work has been
extended to multiple-output problems (Fricker et al., 2013), dynamical
problems (Conti et al., 2009), and non-deterministic (stochastic) simulators
(Oyebamiji et al., 2019). These and other works (Kennedy and O’Hagan, 2001;
Oakley and Youngman, 2017) also outline how to perform Bayesian calibration,
whilst Saltelli-style sensitivity indices (Homma and Saltelli, 1996) have been
developed into a fully-Bayesian sensitivity analysis and applied to GPEs to
better understand input-output relations (Farah and Kottas, 2014). Regarding
its (very numerous) applications, GPE has been used in modelling mitochondrial
DNA deletions (Henderson et al., 2009), influenza epidemic models (Farah et
al., 2014), and microbial communities (Oyebamiji et al., 2019) to name just a
few.
Dynamical GPE may be used to model time-dependent outputs of computer
experiments (Conti et al., 2009; Farah and Kottas, 2014; Mohammadi et al.,
2019; Oyebamiji et al., 2019). Such methods exploit the availability of the
conditional form of the multivariate normal distribution to emulate the
temporal structure of the data. This requires updating the values of outputs
(state vector) at time $t_{i}$ using corresponding time-varying forcing inputs
and the value of the state vector at $t_{i-1}$, as well as static
parameters/initial conditions if any (Conti et al., 2009). We only have access
to the initial conditions (ICs) of the experiments and are unable to use
fully-dynamical emulation methods.
Instead, we approximate the FoS curves with a (1) single quadratic polynomial
and (2) two piecewise quadratic polynomials, and emulate the resulting
model/polynomial parameters by assigning Gaussian process prior distributions
to them. This gives rise to a multi-output hierarchical GP model whereby we
relate the parameters determining FoS curves with the ICs of the computer
experiments. Therefore, the emulator that we developed uses the ICs of the
slopes to predict the FoS behaviour over a time horizon.
We perform Bayesian inference on the unknown parameters using Markov chain
Monte Carlo (MCMC). The resulting posterior distributions allow us to make FoS
predictions and accurately quantify their uncertainty. This work builds on
Svalova et al. (2021) who modelled time to slope failure using the slope
geometry and soil strength using Gaussian process emulation (GPE). Bayesian
emulation of computer experiments is highly cost-effective and allows full
uncertainty analysis (O’Hagan, 2006), which motivates us to use it for FoS
modelling. Predictions for earthworks with specific geometries and materials
can be obtained on-line and used in design, maintenance, and management of
infrastructure.
The structure of this paper is as follows. Section 2 describes the computer
experiments and the simulator or geotechnical model (GM), that simulates
deterioration processes in cut slopes. In Section 3 we describe hierarchical
Bayesian modelling of FoS time series and introduce the quadratic model and
B-spline model of FoS. In Section 4 we propose the Gaussian process emulator
for the quadratic and B-spline models, followed by outlining how we elicited
the corresponding prior distributions. In Section 5 we use the quadratic and
B-spline emulators to model FoS. We obtain within-sample and out-of-sample
posterior distributions of FoS and compare the performance of both models
using the mean squared error $(MSE)$ and the continuous ranked probability
score $(crps)$. We also show how the models can be used to predict time-to-
failure of cut slopes. Section 6 concludes the paper and our findings and
proposes areas for future research. Supplementary Material provides additional
derivations and plots not included in the main manuscript.

##  2 Computer experiments of cut slope deterioration

In this section we describe the computer experiments and the GM (our
simulator) that simulates deterioration processes in cut slopes.

###  2.1 The geotechnical model

This subsection summarises the mechanisms that characterise the slope
deterioration behaviour in the GM (our simulator). The GM is implemented
within FLAC-TP (Fast Lagrangian Analysis of Continua with Two-Phase Flow
(Itasca, 2016)), which treats the soil as a porous medium with variable
saturation and depth-dependent saturated permeability (Postill et al., 2021).
The pore fluids are separated into air and water phases and treated as
immiscible fluids with differing density and viscosity. Water and airflow
velocity is a function of the respective pore fluid pressures and the relative
permeability of the soil. The latter is derived using the van Genuchten-Maulem
relation (van Genuchten, 1980).
The soil mechanical constitutive model uses a non-linear pressure dependent
stiffness to control the bulk and shear moduli. A nonlocal Mohr-Coulomb strain
softening model is used to describe the change in shear strength as a function
of plastic shear strains, whereby the initial peak strength (effective peak
cohesion $c^{\prime}_{p}$ and friction angle $\phi^{\prime}_{p}$) of the
material reduces towards a minimum or residual value with increasing plastic
strain.
The hydrological and mechanical properties were adopted from prior modelling
studies (Potts et al., 1997; Jurečič et al., 2013; Tsiampousi et al., 2017;
Summersgill et al., 2018) and from field and laboratory data (Bromhead and
Dixon, 1986). The effect of weather is included by application of a flux
boundary condition where water is added and removed from the model upper
boundary to replicate the surface water balance driven by precipitation and
evapotranspiration. More details about the GM, including its validation and
quality, can be found in Postill et al. (2020); Rouainia et al. (2020);
Postill et al. (2021), and Helm et al. (2024b).

###  2.2 Initial conditions

While over 40 model ICs can be varied (Itasca, 2016; Rouainia et al., 2020;
Postill et al., 2021), in this paper we only vary five of them. Namely, the
slope height, $h_{s}$, slope angle $\nu$, derived from light detection and
ranging (LiDAR) survey data provided by project stakeholders (Mott MacDonald
and Network Rail), the peak shear strength parameters (effective peak cohesion
$c^{\prime}_{p}$ and friction angle $\phi^{\prime}_{p}$) before the material
has undergone any strength reduction, derived from previous laboratory data
and modelling studies (Apted, 1977; Ellis and O’Brien, 2007), and the
reference coefficient of permeability of the soil at 1 metre depth with
respect to water, $k_{1}^{w}$, derived from field measurements (Dixon et al.,
2019). This simulation study specifically investigated high plasticity
overconsolidated clay, as these were identified by relevant infrastructure
asset owners as being amongst the most vulnerable to slope failures (Mellor et
al., 2017). The five ICs/variables were chosen as previous work has
demonstrated their importance in assessing the deterioration and stability of
infrastructure earthworks (Potts et al., 1997; Ellis and O’Brien, 2007; Mellor
et al., 2017; Rouainia et al., 2020; Postill et al., 2023). This allows us to
limit the experimental time while obtaining an informative training data set.
Additional parameters would increase computation time, which was around three
months for 75 computer runs using 10 machines. To ensure an optimal coverage
of the parameter space, a Latin hypercube experimental design is used to
obtain 75 IC vectors (height, angle, peak cohesion, friction angle, and
permeability). For more information on the adopted GM and the values adopted
for the remaining ICs, see Itasca (2016), Rouainia et al. (2020) and Postill
et al. (2021). Table 1 summarises ranges of the ICs used in the computer
simulations which were selected based on previous studies (Rouainia et al.,
2020; Postill et al., 2021) and expert opinion from partners in the ACHILLES
project. During emulation, permeability was scaled by $10^{8}$ to put the
explanatory variables on a common scale.

Property |  $h_{s}$ (m) |  $\nu$ (degrees) |  $c^{\prime}_{p}$ (kPa) |  $\phi^{\prime}_{p}$ (degrees) |  $\kappa^{w}_{1}$ (m/s)
---|---|---|---|---|---
Range | $[4,20]$ | $[7.6,63.4]$ | $[3,10]$ | $[18.5,25]$ | $[0.145,2.5]\times 10^{-8}$
Table 1: Material and geometry ICs used in the computer experiments. The
ranges were selected based on previous studies and expert opinion from our
project partners. Figure 1: Examples of FoS time series illustrating different
deterioration patterns.

###  2.3 Model run time

The duration of the geotechnical modelling runs varied from approximately 75
minutes to 10 days. The latter depends on a number of factors, including model
geometry, permeability, and the adopted strength parameters. The simulator
runs were implemented until failure or 184 years of model time. Therefore,
some computer runs provide time series of FoS from slope initiation until
slope failure and others are censored at 184 years.

###  2.4 Computer experiments of FoS

Figure 1 illustrates five computer runs of winter FoS which can be used to
understand an earthwork’s proximity to failure (a FoS of unity). Steeply-
declining curves indicate rapid deterioration, whereas curves which have a
plateau over a large number of years are indicative of earthworks undergoing
minimal deterioration. The different behaviours are mainly driven by slope
geometry, for example runs 1, 6, and 31 constitute short (in height) slopes
with low angles, whereas the others are higher and steeper. The geology also
influences deterioration behaviour, whereby cut slopes excavated within weaker
materials will deteriorate more quickly and may also have a lower starting
value of FoS. Finally, the dominant failure mechanism can change over time.
For example, an initially deep seated rotational mechanism (base failure) may
switch to a translational mechanism at a shallower depth within the slope
(landslide) due to near-surface weathering and resultant near-surface strength
reduction (Briggs et al., 2023). This would lead to a change in the FoS
gradient and shape. The time step in our geotechnical model can be between 10s
to 1000s of seconds (here, it is capped at less than or equal to 3600 seconds)
and the weather data is daily. However, the FoS analysis procedure is a time
consuming iterative process. We are emulating yearly measurements in this
work, as more frequent iterations would not be practical in terms of model run
time. Instead, we use an estimate of the annual wettest and driest days, which
broadly correspond to the highest and lowest values of annual FoS. In general,
within a given year, the slopes have the lowest FoS and so are most vulnerable
at the wettest time so the winter results are presented.
Despite the runs in Figure 1 appearing to have a stochastic component, the
computer simulator is deterministic. The FoS behaviour is a function of the
variability of the weather and the effects of antecedence and of the current
state of deterioration. A rainfall or drying event will have a different
effect on slope stability; this effect depends on the initial state of the
slope prior to the event and the adopted permeability. Although each model has
identical boundary conditions, due to the differences in permeability,
geometry, and adopted strength properties, their responses will vary. Once
models undergo significant strength reduction due to deterioration, their FoS
is more strongly affected by seasonal cycles of wetting and drying, which
increases the variability of annual FoS for slopes nearing ULS failure.
In what follows, FoS curves for a set of 75 computer experiment runs of
earthwork (cutting) deterioration (Svalova et al., 2021; Postill et al., 2021;
Helm et al., 2024a) will be modelled as a function of slope height and angle,
soil strength properties (peak cohesion, peak friction angle), and
permeability. This work builds on Svalova et al. (2021) who modelled time to
slope failure using Gaussian process emulation (GPE).

##  3 Hierarchical Bayesian modelling of FoS time series

Hierarchical (multilevel) Bayesian modelling has been referred to as one of
the most important statistical ideas over the past 50 years (Gelman and
Vehtari, 2021), with the key idea of combining/pooling information about
different categories of parameters for better inference. Here, hierarchical
Bayesian modelling can be used to relate the temporal evolution of FoS to the
static GM ICs. We will assume that the FoS time series can be approximated by
a function $g$, the parameters of which will be emulated/related to simulator
ICs.
To simplify some of the methodology, we will be modelling FoS measurements
shifted by a unit, $Y=\,\,$FoS$\,-1$, such that failure is defined by reaching
a (shifted) FoS of zero. For $i\in\\{1,\ldots,N\\}$ and
$j\in\\{1,\ldots,N_{i}\\}$, where $N$ is the number of computer runs and
$N_{i}$ is the number of FoS measurements in run $i$, our hierarchical
Bayesian FoS model is:

$\begin{cases}\text{FoS}_{i,j}=1+g(t_{i,j};\bm{\theta}_{g,i})+\varepsilon_{i,j},\quad\varepsilon_{i,j}\sim\text{N}(0,\sigma^{2}_{i}),\\\
\text{Y}_{i,j}=g(t_{i,j};\bm{\theta}_{g,i})+\varepsilon_{i,j},\\\
\bm{\theta}_{g,i}\sim\text{GP}(m(\bm{x};\bm{\phi}_{1}),V(\bm{x},\bm{x}^{\prime};\bm{\phi}_{2})),\\\
\bm{\phi}_{1}\sim\Phi_{1},\quad\bm{\phi}_{2}\sim\Phi_{2}.\end{cases}$ (1)

In the above, the top level $g+\varepsilon_{i,j}$ is the model for FoS
temporal evolution, and $t_{i,j}$ denotes time in years. We will compare two
$g$ parametrisations, a single and a piecewise quadratic polynomial in the
form of a B-spline. The parameters of $g$,
$\bm{\theta}_{g}=(\bm{\theta}_{g,1},\ldots,\bm{\theta}_{g,N})^{T}$, are
assigned GP prior distributions defined by mean and variance functions $m$ and
$V$. The mean function $m$ is a function of the ICs $\bm{x}$ and
hyperparameters $\bm{\phi}_{1}$; and the covariance function is a function of
a pair of ICs $(\bm{x},\bm{x}^{\prime})$, such that
$V_{i,j}=\text{Cov}(\bm{x}_{i},\bm{x}_{j})$, and hyperparameters
$\bm{\phi}_{2}$. Finally, the hyperparameters (vectors) $\bm{\phi}_{1}$ and
$\bm{\phi}_{2}$ will be assigned some appropriate prior distributions
$\Phi_{1}$ and $\Phi_{2}$. The rest of this Section details the
parameterisation of FoS time series using a quadratic polynomial curve and a
quadratic B-spline with two pieces.
We note that both models are trained on FoS time series with at least 4
measurements, otherwise the B-spline model is not identifiable. This limits
the use of the developed emulator for predicting very sudden failures; however
infrastructure earthworks are typically built to last at least several decades
(Perry et al., 2003), thus this should not be a problem in practice.
Furthermore, we define three types of time-to-failure (TTF). Namely, the model
TTF $(\omega)$, the predicted TTF $(\rho)$, and the true TTF. For computer run
$i$, model TTF $(\omega_{i})$ is defined as the minimum time such that
$g(t)\leq 0$, and predicted TTF $(\rho_{i})$ is defined as the minimum time
such that $g(t)+\varepsilon_{i}\leq 0$. The true TTF is taken from the data.
Some of our computer runs do not reach failure and, for those, the true TTF is
not available.

###  3.1 Quadratic model of FoS

For a computer run $i$ at time $j$, assume that $Y_{i,j}\in\mathbb{R}$ changes
with time $t_{i,j}\in\mathbb{R}^{+}$ through a quadratic relation

$Y_{i,j}=\alpha_{0,i}+\alpha_{1,i}t_{i,j}+\alpha_{2,i}t_{i,j}^{2}+\varepsilon_{i,j},\quad\varepsilon_{i,j}\sim\text{N}(0,\sigma_{i}^{2}),$
(2)

where $a_{0,i}+1$ is the FoS at time zero and $a_{1,i}$ is the FoS gradient at
time zero. Equation (2) can be re-parametrised to include the _model_ TTF for
computer run $i$, $\omega_{i}$:

$Y_{i,j}=\alpha_{0,i}+\alpha_{1,i}t_{i,j}-\frac{\alpha_{1,i}+\alpha_{0,i}/\omega_{i}}{\omega_{i}}t_{i,j}^{2}+\varepsilon_{i,j}.$
(3)

It is useful to rewrite Equation (3) as a quadratic B-spline (De Boor, 1978;
Knott, 2000) with no interior knots. In this form, the quadratic model can be
directly compared to the B-spline FoS model described in Section 3.2.
Therefore,

$\displaystyle Y_{i,j}$
$\displaystyle=g_{1}(t_{i,j};\bm{\theta}_{g_{1},i})+\varepsilon_{i,j},$ (4)
$\displaystyle g_{1}(t_{i,j};\bm{\theta}_{g_{1},i})$
$\displaystyle=\bigg{\\{}\gamma_{0,i}+t_{i,j}\bigg{(}\frac{2\gamma_{1,i}}{\omega_{i}}-\frac{2\gamma_{0,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{\gamma_{0,i}}{\omega_{i}^{2}}-\frac{2\gamma_{1,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}\text{I}_{t_{i,j}}(0;\omega_{i}),$

where $\bm{\theta}_{h_{1},i}=(\gamma_{0,i},\gamma_{1,i},\omega_{i})$,

$\text{I}_{t_{i,j}}(0;\omega_{i})=\begin{cases}1,\quad 0\leq
t_{i,j}<\omega_{i},\\\ 0,\quad\text{otherwise},\end{cases}$

and $\gamma_{0,i}=\alpha_{0,i}$ and
$\gamma_{1,i}=(\alpha_{1,i}\omega_{i}+2\alpha_{0,i})/2$. Here, $\gamma_{0,i}$
is the intercept of $Y_{i,j}$ at $t_{i,j}=0$, and we require $\gamma_{0,i}\geq
0$ to ensure non-negative $Y_{i,j}$ at $t_{i,j}=0$. Thus, we can write
$\gamma_{0,i}=\,$exp(A${}_{0,i})$. Similarly, we set the failure time
parameter $\omega_{i}=\,$exp$(\Omega_{i})$ and
$\sigma_{i}=\,$exp$(\Sigma_{i})$.

####  3.1.1 Constraining the quadratic model

Figure 2: Constraining the quadratic and B-spline models. Plot (left)
illustrates deterministic (noise-free) FoS curves when $\gamma_{1,i}<0$. Plot
(middle) presents B-splines with $\gamma_{1,i}<0,\gamma_{2,i}<0$ and plot
(right) illustrates an unrealistic B-spline with $\gamma_{1,i}\geq
0,\gamma_{2,i}\geq 0$.

The model TTF, $\omega_{i}$, is the smallest positive time that solves the
deterministic part of Equation (4). This can be achieved by constraining
$\gamma_{1,i}>0$. Figure 2 (left) illustrates deterministic (noise-free) FoS
curves when $\gamma_{1,i}$ is unconstrained. From Figure 2 (left) we can see
that the FoS curves can cross the failure-threshold prior to $\omega_{i}$ when
$\gamma_{1,i}<0$. Thus, we write $\gamma_{1,i}=\,$exp($\text{A}_{1,i})$. The
parameter set, to be emulated, of the quadratic model defined by Equation (4),
is
$\bm{\theta}_{g_{1},i}=(\text{A}_{0,i},\text{A}_{1,i},\Omega_{i},\Sigma_{i})$.

###  3.2 B-spline model of FoS

The geometry and location of the critical failure mechanism in a slope may
change over time as deterioration progresses (Postill et al., 2021). This may
lead to a sudden change in the gradient of the FoS deterioration curve. Thus,
a quadratic curve might not be sufficient in explaining such variation,
especially as the earthwork failure time increases. A piecewise polynomial
could be more advantageous over a simple polynomial, whereby local curves are
fitted within non-overlapping regions. One such model is the B-spline (De
Boor, 1978; Knott, 2000), which is a piecewise polynomial defined within a set
of knots/breakpoints that is forced to be smooth and continuous at these
knots. B-splines are commonly used in shape-optimisation methods (Talebitooti
et al., 2015), and can have a desired polynomial degree, which has to be the
same for all piecewise curves (and greater than zero). A quadratic B-spline
with no inner knots is equivalent to a simple quadratic (Equation (4)).
We use a B-spline to fit two piecewise quadratic polynomials to the FoS time
series with one interior knot at half-failure time, $\omega_{i}/2$. This is a
natural extension to the quadratic model in (4) as we are fitting two
polynomials corresponding to a potential change in the slope failure
mechanism, instead of one. What follows provides a brief definition of the
B-spline, the full analytical definition can be found in Supplementary
Material A.
The B-spline approximation of $Y_{i,j}$ is a linear combination of B-spline
basis functions,

$Y_{i,j}=\sum_{l=0}^{m+o-1}\gamma_{l,i}\phi_{l,o-1}(t_{i,j})+\varepsilon_{i,j}=\sum_{l=0}^{3}\gamma_{l,i}\phi_{l,2}(t_{i,j})+\varepsilon_{i,j},\quad\varepsilon_{i,j}\sim\text{N}(0,\sigma^{2}_{i}),$
(5)

where $\phi_{l,2}(t_{i,j})$ are basis functions and $\gamma_{l,i}$ are linear
coefficients. Following De Boor recursive relationships (De Boor, 1978,
Supplementary Material A), the spline model can be rewritten as a sum of two
quadratic functions over non-overlapping intervals,

$\displaystyle Y_{i,j}$
$\displaystyle=g_{2}(t_{i,j};\bm{\theta}_{g_{2},i})+\varepsilon_{i,j},$ (6)
$\displaystyle
g_{2}(t_{i,j};\bm{\theta}_{g_{2},i})=\bigg{\\{}\gamma_{0,i}+t_{i,j}\bigg{(}\frac{4\gamma_{1,i}}{\omega_{i}}-\frac{4\gamma_{0,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{4\gamma_{0,i}}{\omega_{i}^{2}}-\frac{6\gamma_{1,i}}{\omega_{i}^{2}}+\frac{2\gamma_{2,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}\text{I}_{t_{i,j}}(0;\omega_{i}/2)+$
$\displaystyle\bigg{\\{}2(\gamma_{1,i}-\gamma_{2,i}+\gamma_{3,i})+t_{i,j}\bigg{(}\frac{8\gamma_{2,i}}{\omega_{i}}-\frac{4\gamma_{1,i}}{\omega_{i}}-\frac{4\gamma_{3,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{2\gamma_{1,i}}{\omega_{i}^{2}}-\frac{6\gamma_{2,i}}{\omega_{i}^{2}}+\frac{4\gamma_{3,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}\text{I}_{t_{i,j}}(\omega_{i}/2;\omega_{i})$

where
$\bm{\theta}_{g_{2},i}=(\gamma_{0,i},\gamma_{1,i},\gamma_{2,i},\gamma_{3,i},\omega_{i})$,

$\text{I}_{t_{i,j}}(a;b)=\begin{cases}1,\quad\text{if }a\leq t_{i,j}<b,\\\
0,\quad\text{otherwise},\end{cases}$

and $\varepsilon_{i,j}\sim\text{N}(0,\sigma^{2}_{i})$. Removing $\gamma_{3,i}$
from Equation (LABEL:B-spline_expanded_full) forces the B-spline curve to
cross zero at the last knot, which is fitting for modelling the FoS. Then, our
B-spline model becomes

$\displaystyle Y_{i,j}$
$\displaystyle=g_{2}(t_{i,j};\bm{\theta}_{g_{2},i})+\varepsilon_{i,j},$ (7)
$\displaystyle
g_{2}(t_{i,j};\bm{\theta}_{g_{2},i})=\bigg{\\{}\gamma_{0,i}+t_{i,j}\bigg{(}\frac{4\gamma_{1,i}}{\omega_{i}}-\frac{4\gamma_{0,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{4\gamma_{0,i}}{\omega_{i}^{2}}-\frac{6\gamma_{1,i}}{\omega_{i}^{2}}+\frac{2\gamma_{2,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}\text{I}_{t_{i,j}}(0;\omega_{i}/2)+$
$\displaystyle\bigg{\\{}2(\gamma_{1,i}-\gamma_{2,i})+t_{i,j}\bigg{(}\frac{8\gamma_{2,i}}{\omega_{i}}-\frac{4\gamma_{1,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{2\gamma_{1,i}}{\omega_{i}^{2}}-\frac{6\gamma_{2,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}\text{I}_{t_{i,j}}(\omega_{i}/2;\omega_{i}).$

where
$\bm{\theta}_{g_{2},i}=(\gamma_{0,i},\gamma_{1,i},\gamma_{2,i},\omega_{i})$.
As before, $\gamma_{0,i}$ is the intercept of $Y_{i,j}$ at $t_{i,j}=0$ and
$\gamma_{0,i}\geq 0$ to ensure non-negative $Y_{i,j}$ at $t_{i,j}=0$. We
define $\gamma_{0,i}=\,$exp(A${}_{0,i})$ as with the quadratic model.
Similarly, $\omega_{i}=\,$exp$(\Omega_{i})$ and
$\sigma_{i}=\,$exp$(\Sigma_{i})$.
The B-spline model has an additional coefficient $\gamma_{2,i}=\,$exp(A2,i),
that arises due to the second polynomial piece. Its interpretation is
difficult, but we can get an intuition into its effect on the B-spline
graphically (Section 3.2.1).

####  3.2.1 Constraining the B-spline model

The model TTF $\omega_{i}$ is the smallest positive time that solves the
deterministic part of (7). Figure 2 (middle) illustrates deterministic (noise-
free) FoS curves when $\gamma_{1,i}$ and $\gamma_{2,i}$ are unconstrained.
From Figure 2 (middle) we can see that the FoS curves can cross zero (failure
threshold) prior to the model TTF. This behaviour can be avoided by
constraining $\gamma_{1,i},\gamma_{2,i}>0$, thus we set
$\gamma_{1,i}=\,$exp($\text{A}_{1,i})$ and $\gamma_{2,i}=\,$exp(A${}_{2,i})$.

Figure 2 (right), however, is not plausible for modelling FoS as the peak of
the second polynomial is higher than the initial FoS, $\gamma_{0,i}$. Indeed,
Figure 2 (right) would imply a very dramatic increase in soil strength after
an initial decrease. This would typically only occur if there were significant
changes made to the slope/soil such as geotechnical intervention measures,
which may take the form of the installation of soil nails, piles or retaining
structures or replacement of the soil with a higher strength material (see
Nowak, 2012). In the absence of such interventions this large increase in FoS
should not occur. Additionally, allowing this behaviour can cause difficulties
in estimating the the B-spline coefficients and using the resulting model for
prediction. Therefore, we impose an additional constraint, that the FoS time
series is nonincreasing after the internal knot, i.e. $Y_{i,j}(t_{2})\leq
Y_{i,j}(t_{1})$ for $t_{1}\leq\omega/2\leq t_{2}$. This constraint eliminates
FoS curves that are not plausible while still allowing us to adequately model
the FoS computer experiment data. The parameter set, to be emulated, of the
B-spline model, defined by Equation (7), is
$\bm{\theta}_{g_{2},i}=(\text{A}_{0,i},\text{A}_{1,i},\text{A}_{2,i},\Omega_{i},\Sigma_{i})$.

##  4 Gaussian process emulator

In this section, we define a Bayesian GPE for the B-spline model parameters,
$\bm{\theta}_{g_{2},i}=(\text{A}_{0,i},\text{A}_{1,i},\text{A}_{2,i},\Omega_{i},\Sigma_{i})$,
(see Equation (7)). Let
$\textbf{A}_{0}=(\text{A}_{0,1},\ldots,\text{A}_{0,N})$,
$\textbf{A}_{1}=(\text{A}_{1,1},\ldots,\text{A}_{1,N})$,
$\textbf{A}_{2}=(\text{A}_{2,1},\ldots,\text{A}_{2,N})$,
$\bm{\Omega}=(\Omega_{1,1},\ldots,\Omega_{1,N})$ and
$\bm{\Sigma}=(\Sigma_{1,1},\ldots,\Sigma_{1,N})$. The Gaussian process
emulator for the quadratic model is defined identically to the GPE for the
B-spline model defined in this section, excluding the parameter vector
$\text{{A}}_{2}$ and the associated hyperparameters $\bm{\beta}_{2}$ and
$\tau_{2}$. For a computer run $i$, $\bm{x}_{i}=(x_{i,1},\ldots,x_{i,5})^{T}$
is a set of ICs corresponding to height $x_{1}$ (m), angle $x_{2}$ (degrees),
cohesion $x_{3}$ (kPa), friction angle $x_{4}$ (degrees), and permeability
$x_{5}$ (m/s). We define $\bm{z}_{i}=(z_{i,1},\ldots,z_{i,5})^{T}$, where

$z_{i,k}=\frac{x_{i,k}-\bm{\bar{x}}_{\text{training},k}}{\sigma_{\text{training},k}}\text{
for }k\in\\{1,3,4,5\\},\text{ and
}z_{i,2}=\frac{x_{i,2}-\bm{\bar{x}}_{\text{training},2}}{1.5\times\sigma_{\text{training},2}},$
(8)

where $\bm{x}_{\text{training},k}=(x_{1,k},\dots x_{N,k})$, $N=75$ is the
number of training samples and $\bm{\bar{x}}_{\text{training},k}$ and
$\sigma_{\text{training},k}$ correspond to the mean and standard deviation of
$\bm{x}_{\text{training},k}$, respectively. This standardisation puts all ICs
on the same scale. For a computer run $i$, the FoS parameters are related to
$\bm{z}_{i}$ following a normal distribution

$\begin{pmatrix}\text{A}_{0,i}\\\ \text{A}_{1,i}\\\ \text{A}_{2,i}\\\
\Omega_{i}\\\ \Sigma_{i}\end{pmatrix}\sim
N\left(\begin{bmatrix}h(\bm{z}_{i})^{T}\bm{\beta}_{0}\\\
h(\bm{z}_{i})^{T}\bm{\beta}_{1}\\\ h(\bm{z}_{i})^{T}\bm{\beta}_{2}\\\
h(\bm{z}_{i})^{T}\bm{\beta}_{\Omega}\\\
h(\bm{z}_{i})^{T}\bm{\beta}_{\Sigma}\\\
\end{bmatrix},\begin{bmatrix}\tau_{0}&0&0&0&0\\\ 0&\tau_{1}&0&0&0\\\
0&0&\tau_{2}&0&0\\\ 0&0&0&\tau_{\Omega}&0\\\ 0&0&0&0&\tau_{\Sigma}\\\
\end{bmatrix}\right).$ (9)

The parameter vectors
$\bm{\beta}_{l}=(\beta_{l,0},\ldots,\beta_{l,5})^{T},l\in\\{0,1,2,\Omega,\Sigma\\}$
are regressor coefficients and $h(\bm{z}_{i})=(1,z_{i,1},\ldots,z_{i,5})^{T}$
is a regressor function of the ICs. The variance matrix in (9) is diagonal
with marginal variances $\tau_{k}$ and zero prior correlation between
$\text{A}_{0,i},\text{A}_{1,i},\text{A}_{2,i},\Omega_{i}$, and $\Sigma_{i}$,
which therefore can be expressed as individual normal distributions.
Transforming the data from $\bm{x}_{i}$ to $\bm{z}_{i}$ makes the B-spline
parameters ($\text{A}_{0,i},\text{A}_{1,i},\text{A}_{2,i},\Omega_{i}$ and
$\Sigma_{i}$) easier to interpret. For example, when $\bm{z}_{i}=\bm{0}$,
$\exp(\beta_{\Omega,0})$ can now be interpreted as the expected failure time
of a slope with “average” characteristics. On the other hand, low or zero
values for some of the (untransformed) ICs have a non-linear relationship with
deterioration, for example, zero cohesion would imply a slope buillt with “dry
sand” which becomes difficult to model and interpret. Jointly, A0, A1, A2,
$\bm{\Omega}$, and $\bm{\Sigma}$ follow a matrix-normal distribution,

$\displaystyle\text{vec}(\text{{A}}_{0},\text{{A}}_{1},\text{{A}}_{2},\bm{\Omega},\bm{\Sigma})\sim\text{N}\left(\text{vec}(H\bm{\beta}_{0},H\bm{\beta}_{1},H\bm{\beta}_{2},H\bm{\beta}_{\Omega},H\bm{\beta}_{\Sigma}),\bm{V}\otimes\bm{U}\right),$
(10)
$\displaystyle\bm{U}=C(\bm{x},\bm{x}^{\prime},\bm{\delta},\zeta),\quad\bm{V}=\text{diag}(\tau_{0},\tau_{1},\tau_{2},\tau_{\Omega},\tau_{\Sigma}).$

In the above, $H$ is the regressor matrix, whose $i$th row is $h(\bm{z}_{i})$,
$\bm{U}$ and $\bm{V}$ are the row- and column-wise covariance matrices,
respectively. We use a non-isotropic Gaussian correlation function
$C(\bm{z},\bm{z}^{\prime},\bm{\delta},\zeta)=\exp\left\\{-\sum_{i=1}^{5}\frac{(z_{i}-z_{i}^{\prime})^{2}}{\delta_{i}^{2}}\right\\}+\zeta\mathbb{I}(\bm{z},\bm{z}^{\prime})$,
where $\bm{\delta}=(\delta_{1},\ldots,\delta_{5})^{T}$ is a vector of
correlation lengths. The nugget (Andrianakis and Challenor, 2012) $\zeta$ is
added to the diagonal of $C$ to improve numerical stability of inversion
calculations, and the function $\mathbb{I}(\bm{z},\bm{z}^{\prime})$ is an
indicator for the event $\bm{z}=\bm{z}^{\prime}$. Note that as A0, A1, A2,
$\bm{\Omega}$ and $\bm{\Sigma}$ are uncorrelated a priori, their joint prior
distribution can be written as a product of the following marginal
distributions

$\displaystyle\text{{A}}_{0}\sim\text{N}\left(H\bm{\beta}_{0},\tau_{0}\bm{U}\right),\quad\text{{A}}_{1}\sim\text{N}\left(H\bm{\beta}_{1},\tau_{1}\bm{U}\right),\quad\text{{A}}_{2}\sim\text{N}\left(H\bm{\beta}_{2},\tau_{2}\bm{U}\right),$
(11)
$\displaystyle\bm{\Omega}\sim\text{N}\left(H\bm{\beta}_{\Omega},\tau_{\Omega}\bm{U}\right),\quad\bm{\Sigma}\sim\text{N}\left(H\bm{\beta}_{\Sigma},\tau_{\Sigma}\bm{U}\right).$

Concerning the hierarchical model given by Equation (1), we have defined two
models for FoS temporal evolution; the quadratic model, defined by Equation
(4), and the B-spline model, defined by Equation (7). Furthermore, Equation
(11) completes the second level of the hierarchical model, assigning GP priors
to the temporal model parameters. Removing $\text{{A}}_{2}$ from Equation (11)
provides the GPE for the quadratic model. In Sections 4.1 and 4.2 we complete
the hierarchical model defined by Equation (1) by assigning prior
distributions to the hyperparameters of the B-spline and quadratic model
parameters, respectively.

###  4.1 Prior distributions for the B-spline model

This section describes how we elicited prior distributions for the B-spline
models, using engineering judgement where applicable. The IC $\bm{z}=\bm{0}$
corresponds to a slope with $x_{1}=9.77$ (height), $x_{2}=20$ (angle),
$x_{3}=6.81$ (cohesion), $x_{4}=21.7$ (friction) and $x_{5}=0.77$
(permeability). Using Equation (9), $\bm{z}=\bm{0}$ corresponds to a slope
with “average” characteristics of our training data with the following prior
distributions

$\displaystyle\text{A}_{0}(\bm{0})\sim\text{N}(\beta_{0,0},\tau_{0}),\quad\text{A}_{1}(\bm{0})\sim\text{N}(\beta_{1,0},\tau_{1}),\quad\text{A}_{2}(\bm{0})\sim\text{N}(\beta_{2,0},\tau_{2}),$
(12)
$\displaystyle\Omega(\bm{0})\sim\text{N}(\beta_{\omega,0},\tau_{\omega}),\quad\Sigma(\bm{0})\sim\text{N}(\beta_{\Sigma,0},\tau_{\Sigma}).$

The means of distributions in Equation (12) reduce to intercept-only instead
of the full regressor function form. These $\beta_{l,0}$,
$l\in\\{0,1,2,\Omega,\Sigma\\}$, represent our uncertainty about the initial
FoS averaged across $\bm{x}$. By generating synthetic data and ensuring
reasonable relations between the ICs and the resulting FoS behaviour in the
prior predictive distributions, we elicited the following prior distributions

$\displaystyle\beta_{0,0}\sim\text{N}(\text{log}(1),0.5),\quad\beta_{1,0}\sim\text{N}(\text{log}(0.6),0.4),\quad\beta_{2,0}\sim\text{N}(-0.5,2.5),$
(13)
$\displaystyle\beta_{\omega,0}\sim\text{N}(5.25,1),\quad\beta_{\Sigma,0}\sim\text{N}(\text{log}(0.1),0.5).$

In addition,

$\beta_{0,i}\sim\text{N}(0,0.5),\quad\beta_{1,i}\sim\text{N}(0,0.5),\quad\beta_{2,i}\sim\text{N}(0,1),\quad\beta_{\omega,i}\sim\text{N}(0,1),\quad\beta_{\Sigma_{i}}\sim\text{N}(0,0.5),$
(14)

for $i=1,\dots,5$. Table 2 presents prior mean and central $95\%$ credible
intervals for the expected initial FoS, model failure time and noise using the
priors given by Equation (13). The shape parameters $\text{A}_{1}$ and
$\text{A}_{2}$ have no intuitive interpretation and are not included in Table
2. From Table 2, we see that a priori we expect the initial FoS for a slope
with $\bm{z}=\bm{0}$ to be between 1.38 and 3.66, $95\%$ of the time with a
mean of 2; $\omega_{i}$ is expected to be between 26.8 days and 1350 days,
$95\%$ of the time with a mean of 191 days; and we expect the noise to be
between 0.0375 and 0.266, $95\%$ of the time with a mean of 0.1. These priors
are vaguely informative and represent our uncertainty about the parameters of
the B-spline model for a slope with $\bm{z}=0$ and put prior weighting in
regions of the parameter space that produce plausible synthetic data.

Prior | Mean | Central $95\%$ Credible Interval
---|---|---
$\text{E}[\gamma_{0}]+1$ | 2.00 | $(1.38,3.66)$
$\text{E}[\omega]$ | 191.0 | $(26.8,1350.0)$
$\text{E}[\sigma]$ | 0.1000 | $(0.0375,0.2660)$
Table 2: Prior mean and central $95\%$ credible intervals for
$\text{E}[\gamma_{0}]+1$, $\text{E}[\omega]$ and $\text{E}[\sigma]$.

Similarly to Svalova et al. (2021), we selected vague prior distributions for
the covariance function parameters. Specifically,

$\tau_{l}\sim\text{IGa}(3,0.5),\quad\delta_{k}\sim\text{Exp}(0.2),$ (15)

for $l\in\\{0,1,2,\Omega,\Sigma\\}$ and $k=1,\dots,5$. Figure 3 illustrates
some example prior deterministic deterioration curves with $\gamma_{0}=1$ (the
prior mean for a slope with $\bm{z}=\bm{0}$) and $\Omega=5.25$ (the prior mean
for a slope with $\bm{z}=\bm{0}$). Figure 3 (left) illustrates the mean
deterministic deterioration curve (dashed line) and the central $95\%$
prediction interval (solid lines) when sampling from $\beta_{1,0}$ and
$\beta_{2,0}$ in Equation (13); the red line corresponds to
$\gamma_{1},\gamma_{2}=0$; and the evolution depicted by the blue line is a
B-spline curve we believe is not plausible. From Figure 3 (left) we can see
that the priors specified on $\beta_{1,0}$ and $\beta_{2,0}$ are vague. These
priors represent our uncertainty in the evolution of the FoS of the slope with
$\bm{z}=\bm{0}$. It is difficult to specify how the FoS of a slope with
$\bm{z}=\bm{0}$ will evolve over time. However, there is little weighting
between the 0.025 quantile and the red line, as we do not expect curves to
deteriorate rapidly (up to the knot at $\omega/2$, shown by the vertical black
line) and then asymptotically deteriorate to zero. Through simulation of
synthetic data we were able to elicit priors that put little weight on
implausible curves, but there is no published data on how the FoS of a slope
with $\bm{z}=\bm{0}$ will evolve over time.
Figure 3 (right) presents some example deterministic FoS time series sampled
from $\beta_{1,0}$ and $\beta_{2,0}$ in Equation (13). The priors specified on
$\beta_{1,0}$ and $\beta_{2,0}$ allow for a variety of time series of slope
deterioration that we see in reality. The priors for $\beta_{1,0}$ and
$\beta_{2,0}$ govern the evolution of deterioration for a slope with
$\bm{z}=\bm{0}$; and these priors put little weighting on deterioration curves
that we do not expect to see in reality (such as the red and blue curves in
Figure 3 (left)), but it is difficult to be more specific than that.

Figure 3: Some example prior deterioration curves with
$\gamma_{0}=1,\Omega=5.25$ sampling from $\beta_{1,0}$ and $\beta_{2,0}$ in
Equation (13).

The remaining parameters, i.e., $\beta_{l,i}$, for
$l\in\\{0,1,2,\Omega,\Sigma\\}$ and $i=1,\dots,5$, were given vaguely
informative priors centered around zero, as can be seen by Equation (14).

The posterior distributions for the GPE parameters are intractable and will be
sampled using Markov chain methods.

###  4.2 Prior distributions for the quadratic model

The prior distributions for the Gaussian process emulator for the quadratic
model are identical to those used under the B-spline model. Once again, the
posterior distributions for the GPE parameters are intractable and will be
sampled using Markov chain methods.
Figure 4 presents some example prior deterministic deterioration curves with
$\gamma_{0}=1$ (the prior mean for a slope with $\bm{z}=\bm{0}$) and
$\Omega=5.25$ (the prior mean for a slope with $\bm{z}=\bm{0}$) . Figure 4
(left) presents the mean deterministic deterioration curve (dashed line) and
the central $95\%$ prediction interval (solid lines) when sampling from
$\beta_{1,0}$ in Equation (13); the red line corresponds to $\gamma_{1}=0$;
and the evolution depicted by the blue line is a quadratic curve we believe is
not plausible. From Figure 4 (left) we can see that the prior specified on
$\beta_{1,0}$ is vague. This prior represents our uncertainty in the evolution
of the FoS of the slope with $\bm{z}=\bm{0}$. Similarly to Section 4.1,
through simulation of synthetic data, we were able to elicit a prior that put
little weight on implausible curves, but there is no published data on how the
FoS of a slope with $\bm{z}=\bm{0}$ will evolve over time. Figure 4 (right)
presents some example deterministic FoS time series sampled from $\beta_{1,0}$
in Equation (13). The prior specified on $\beta_{1,0}$ allows for a variety of
time series of slope deterioration that we see in reality. Similarly to
Section 4.1, the prior for $\beta_{1,0}$ governs the evolution of
deterioration for a slope with $\bm{z}=\bm{0}$; and this prior puts little
weighting on deterioration curves that we do not expect to see in reality
(such as the red and blue curves in Figure 4 (left)), but it is difficult to
be more specific than that.

Figure 4: Example prior deterioration curves with $\gamma_{0}=1,\Omega=5.25$
sampling from $\beta_{1,0}$ in Equation (13).

##  5 Results

We obtained posterior distributions of the two hierarchical models using
RStan, the R interface to Stan (Carpenter et al., 2017; Stan Development Team,
2022), with 4 chains. The quadratic and B-spline emulators took approximately
5 and 33 hours, respectively, to run using 4 cores with Intel CPU (Xeon,
E5-2699 v4, base frequency 2.2GHz) through the Slurm workload manager Jette
and Wickberg (2023) on Newcastle University’s high performance computing
service, Rocket.
Figure B1-B3 in the Supplementary Material B present some example trace plots,
for the upper-level parameters, under the quadratic and B-spline models,
respectively. In addition, Figures B5 and B6 present trace plots for the
emulator parameters under the quadratic and B-spline models, respectively. We
can see that all four chains have converged to the same values for the
emulator parameters under the quadratic and B-spline models, respectively.
Similar conclusions can be drawn from Figures B1 and B2, although FoS run 33
in Figure C5 may hint at multi-modality under the quadratic model or that the
chains have not converged to the stationary distribution. Some trace plots
show multi-modality under the B-spline model (see Figure B3). Figure B4
presents three within-sample computer runs with (at least one) multimodal
(upper level) posterior distribution; we overlay a sample posterior FoS curve
sampled from each mode. The chains in FoS run 34 switch between modes and the
chains in FoS runs 39 and 68 do not switch between modes. The three computer
runs in Figure B4 all follow a similar time series, with a relatively stable
initial period followed by a rapid decline in FoS. Both modes fit reasonably
well to the data. It may be possible to eliminate one mode by penalising FoS
curves that are non-decreasing (by penalising values of
$\gamma_{1}>\gamma_{0}$).

###  5.1 Posterior distributions of FoS

Figure 5, and Figures C1 - C8 in the Supplementary Material C, present
posterior distributions of FoS obtained under the quadratic model (blue) and
the B-spline model (red), for examples of computer experiment data. From
Figure 5, and Figures C1 - C8, we can see that the B-spline model captures the
FoS trends well. The quadratic model and the B-spline model perform
approximately the same for computer runs 12, 14 and 26, within the range of
the data (5). However, the quadratic model is unable to capture the trends of
computer runs 1, 18, 22, 30, and 64. Similar conclusions can be drawn from
Figures C1 - C8 in the Supplementary Material C, e.g. run 43 in Figure C1.
This is expected as some computer runs have trends that are not adequately
modelled by a quadratic curve.

Figure 5: Examples of posterior distributions of FoS using the quadratic model
(blue) and the B-spline model (red). The dashed lines represent the posterior
means and the solid lines represent the central posterior 95% prediction
intervals.

Figure 6 presents the difference in the mean squared error $(MSE)$ and the
continuous ranked probability score $(crps)$ (Gneiting et al., 2007) between
the quadratic model ($MSE^{Q}$ and $crps^{Q}$ respectively) and the B-spline
model ($MSE^{BS}$ and $crps^{BS}$ respectively) for the nine within-sample
computer runs presented in Figure 5. The $crps$ is calculated at each time
point for each computer run and is defined in terms of the predictive CDF,
$F_{t}$, and an observation, $x_{t}$,

$crps(F_{t},x_{t})=\int_{-\infty}^{\infty}\\{F_{t}(y)-\textbf{I}(y\geq
x_{t})\\}^{2}dy,$ (16)

where $F_{t}$ is the CDF of FoS at time $t$ and $x_{t}$ is the observed FoS at
time $t$. The $crps$ is a proper scoring rule that evaluates the performance
of a full predictive density.

Figure 6:  Difference between the quadratic model and the B-spline model for
nine within-sample computer runs in the in the $MSE$ (left) and $crps$
(right).

As expected for computer runs 12, 14 and 26 the quadratic model and the
B-spline model perform approximately the same, as their distributions are
narrowly concentrated around zero. On the other hand, the $MSE$ and the $crps$
are larger for the quadratic model for computer runs 1, 18, 22, 30, 43 and 64
(c.f. Figure 5). For the latter runs, the boxplots are centred in the positive
range and have a longer tail in the same direction. Figure D1 in Supplementary
Material D provides the difference in the $MSE$ and the $crps$ between the
quadratic model and the B-spline model for the remaining sixty-six within-
sample computer runs (Figures C1 - C8).

###  5.2 Posterior distributions of TTF

We can obtain posterior distributions of the _predicted_ TTF, $\rho_{i}$, for
$i=1,\dots,N$, by simulating FoS time series for every posterior draw of the
model parameters and obtaining a distribution of the first time when FoS
reaches failure, in most cases due to the error term. The distributions of
predicted TTF, $\rho_{i}$, are presented in Figure 7 and Figures E1 - E8 in
the Supplementary Material E. The solid red lines indicate true TTF for
computer runs that reached failure.

Figure 7: Examples of posterior distributions of predicted TTF, $\rho_{i}$,
using the quadratic model (blue) and the B-spline model (red). The dashed
lines represent the central posterior 95% prediction intervals under each
model. The solid red lines indicate true TTF for computer runs that reached
failure.

From Figure 7, and Figures E1 - E8, we can see that the TTF predictions are
noticeably different under the quadratic model compared to the B-spline model
for some computer runs. In general, the B-spline model leads to a posterior
uncertainty decrease in the predicted TTF, and in some cases, a dramatic one.
Figure 8 presents examples of posterior distributions of FoS using the
quadratic model (blue) and the B-spline model (red) for the nine computer runs
shown in Figure 7, extrapolating beyond the range of the data. The quadratic
model is unable to capture some of the trends in the FoS data, such as
computer runs 1, 18, 22, 30, and 64 in Figure 8 or run 43 in Figure E1. For
computer run 1, we can see there is a substantial difference in the TTF
prediction between the two models. Computer run 1 decays for around 50-75
years and then reaches a stable plateau. The quantiles of the quadratic model
fail to capture this behaviour, whereas the B-spline model captures the decay
and stability periods well. The B-spline model then predicts the slope will
begin to deteriorate to failure. This is a consequence of having the internal
knot fixed at $k=\omega/2$.

Figure 8: Examples of posterior distributions of FoS using the quadratic model
(blue) and the B-spline model (red). The dashed lines represent the posterior
means and the solid lines represent the central posterior 95% prediction
intervals.

Computer run 1 could remain stable for a long time after year 184, as
illustrated in Figure E9 (Supplementary Material E). Emulating the internal
knot, $k$, as well as the model failure time, $\omega$, would allow for time
series such as those shown in Figure E9. At present, we avoid modelling the
internal knot as our training data is of moderate size, but this should be
considered in future work.

###  5.3 Model validation

We also validated both models using 20 held-out computer experiment runs. As
GPs are closed under conditioning, it is possible to derive an analytical
expression for the GPE conditioned on a set of computer training runs. Assume
a (finite) collection of $n$ “observed” experimental outputs,
$\bm{y}=\big{(}\text{A}_{0}(\bm{z}_{1}),\text{A}_{0}(\bm{z}_{2}),\dots,\text{A}_{0}(\bm{z}_{n})\big{)}$
performed on the ICs $\bm{z}_{1},\bm{z}_{2},\dots,\bm{z}_{n}$, where $n=75$,
which comprise the training data. We assume no repeated inputs, so
$\bm{z}_{i}=\bm{z}_{j}$ if and only if $i=j$. The $n-$vector $\bm{y}$ follows
a multivariate normal distribution,

$\bm{y}\mid\bm{\beta}_{0},\tau_{0},\bm{\delta},\zeta\sim\text{N}(H_{z}\bm{\beta}_{0},\tau_{0}\Sigma_{z}),$
(17)

where $H_{z}$ is a matrix of regressors whose $i$th row is $h(\bm{z}_{i})$,
$\Sigma_{z}(i,j)=C(\bm{z}_{i},\bm{z}_{j})+\zeta\mathbb{I}(\bm{z},\bm{z}^{\prime})$,
$\zeta$ is the nugget and the function $\mathbb{I}(\bm{z},\bm{z}^{\prime})$ is
an indicator for the event $\bm{z}=\bm{z}^{\prime}$. Using standard rules for
conditioning on a subset of observations (Gramacy, 2020), we find

$\displaystyle\text{A}_{0}(\cdot)|\bm{y},\bm{\beta}_{0},\tau_{0},\bm{\delta},\zeta\sim\text{GP}(m^{*}(\cdot),V^{*}(\cdot,\cdot)),$
(18) $\displaystyle
m^{*}(\bm{z})=h(\bm{z})^{\top}\bm{\beta}_{0}+t(\bm{z})^{\top}\Sigma^{-1}_{z}(\bm{y}-H_{z}\bm{\beta}_{0}),\quad
V^{*}(\bm{z},\bm{z}^{\prime})=\tau_{0}\big{(}C(\bm{z},\bm{z}^{\prime},\bm{\delta})-t(\bm{z})^{\top}\Sigma^{-1}_{z}t(\bm{z}^{\prime})\big{)},$

where
$t(\bm{z})=\big{(}C(\bm{z},\bm{z}_{1},\bm{\delta}),C(\bm{z},\bm{z}_{2},\bm{\delta}),\dots,C(\bm{z},\bm{z}_{n},\bm{\delta})\big{)}$
is a column vector of correlations between the (generic) emulator IC $\bm{z}$
and training ICs $\bm{z}_{1},\bm{z}_{2},\dots,\bm{z}_{n}$. Similar expressions
can be derived for the conditional distributions of
$\text{A}_{1}(\cdot),\text{A}_{2}(\cdot),\Omega(\cdot)$ and $\Sigma(\cdot)$.
Using Equation (18) for a held-out computer run with $\bm{z}$, provides us
with a distribution for
$\text{A}_{0}(\bm{z})|\bm{y},\bm{\beta}_{0},\tau_{0},\bm{\delta},\zeta$, and
this procedure can be repeated for
$\text{A}_{1}(\bm{z}),\text{A}_{2}(\bm{z}),\Omega(\bm{z})$ and
$\Sigma(\bm{z})$. Using these distributions, we can sample from Equation (7),
and obtain FoS time series for a computer run with some IC $\bm{z}$.
Figure 9 and Figures F1 - F2 in Supplementary Material F present posterior
distributions of FoS under the quadratic (blue) and B-spline (red) models for
20 held-out computer runs. For both models, the posterior means capture the
FoS behaviour well. Validation runs 1 (Figure F1) and 2 (Figure 9) have the
widest 95% highest posterior density intervals (HPDI) under both models. These
computer runs are quite unusual, as they both begin with high initial FoS
measurements which then increase before plateauing at values larger than their
initial strength. Both models do not reach failure at the end of simulation
and at 184 years (the simulation end time) both computer experiements have
higher FoS measurements than their starting FoS. This behaviour leads to a
large uncertainty in the model and predicted TTF and in the posterior
predictions of FoS. Despite this uncertainty, the mean predicted FoS is
extremely close to the validation data. For the remaining runs in Figure 9 and
Figures F1 - F2 in Supplementary Material F, the 95% HPDIs are much narrower
(under both models), and the mean posterior behaviour is close to the true
values for most out-of-sample computer runs.

Figure 9: Examples of posterior distributions of FoS using the quadratic model
(blue) and the B-spline model (red) for out-of sample computer runs. The
dashed lines represent the posterior means and the solid lines represent the
central posterior 95% prediction intervals. Figure 10: Plots presenting the
difference in the $MSE$ and $crps$ between the quadratic model and the
B-spline model for 20 out-of-sample computer runs.

Figure 10 presents the difference in the $MSE$ and the $crps$ between the
quadratic model and the B-spline model for the 20 out-of-sample computer runs
presented in Figure 9, and Figures F1 - F2 in Supplementary Material F. The
differences in the $MSE$ are small for most out-of-sample runs. For computer
runs e.g. 5 (Figure F1), 11, and 12, the mean behaviour of the B-spline model
more closely captures the FoS behaviour compared to the quadratic model. The
opposite can be observed for computer runs 3 (Figure F1) and 6, although for
run 3 the mean FoS at failure time is closer to the true FoS than that of the
quadratic model. For runs with unoberved failure times (e.g. 7, 12 (Figure 9),
and 16 (Figure F1)), the B-spline model tends to provide a more conservative
estimate of the FoS at the last time point. An exception is run 14 where the
mean FoS is marginally larger.
The $95\%$ prediction intervals tend to be narrower under the B-spline model,
except for failure times of 10 years or less. Computer run 16 (Figure F1)
illustrates behaviour that is difficult to model, as it begins with a stable
period followed by a sudden, rapid deterioration before an increase in FoS.
Most of its measurements are outside of the quadratic model prediction
intervals from initiation until around 100 years. The corresponding $95\%$
interval is wider under the B-spline model, however it captures all of the FOS
measurements.

##  6 Conclusions

Factor of safety is an important diagnostic tool to predict and monitor the
safety of transport infrastructure. For a given slope, the change in FoS can
be estimated using computer experiments of deterioration, however they are
very time-consuming to run. Here, we built a GP emulator of FoS for cuttings
in high-plasticity soils using Bayesian hierarchical modelling to approximate
FoS measurements faster for any IC combination. We approximated the FoS curves
with a single quadratic model and a B-spline model of two quadratic
polynomials with one interior knot; and we emulated the model parameters by
assigning Gaussian process prior distributions to them. This produced a multi-
output hierarchical GP model that relates the parameters determining FoS
curves with the computer experiment ICs.
We obtained within-sample posterior distributions of FoS and demonstrated how
the models could be used to obtain posterior distributions of predicted TTF.
The quadratic model was able to emulate some computer experiments well, as
several FoS time series deteriorated with quadratic trends. However, in many
cases it did not have the flexibility required to emulate FoS curves that were
driven by more than one failure mechanism. For example, computer run 1 (Figure
5), decays until around year 40, and then stabilises for approximately 20
years, before deteriorating to failure; and computer run 64 remains stable for
approximately 75-100 years, before rapidly declining. These time series are
much better captured under the B-spline model. We compared model performance
using the $MSE$ and the $crps$ and found that the B-spline model performed as
well as or outperformed the quadratic model for most computer runs.
Furthermore, we illustrated how to use the models for out-of-sample prediction
by conditioning on the “observed” upper level model parameters. The successful
emulation of geotechnical models of deterioration of infrastructure slopes,
has the potential to inform slope design, maintenance and remediation by
introducing the time dependency of deterioration into geotechnical asset
management.
Future developments in our work include learning the interior knot of the
B-spline from the data rather than fixing it at half-model TTF, as setting
$k=\omega/2$ limits the B-spline’s predictive capability outside of the range
of the data. For example, for computer run 1, the B-spline model predicts the
FoS time series will deteriorate rapidly to failure and predicts the expected
failure time of computer run 1 to be 318 years with a central 95% prediction
interval of $(295,347)$ years. Computer run 1 could remain stable for a
relatively long time after year 184, as illustrated in Figure E9. Other
developments would include emulating the FoS using linear dynamical models and
Gaussian processes, for example, as in Oyebamiji et al. (2017), as well as
including wider ranges of soil properties and earthwork types in the computer
model.

## References

  * Abbott (2018) Abbott, S. (2018, June).  Earthworks technical strategy.  Technical report, Network Rail.
  * Andrianakis and Challenor (2012) Andrianakis, I. and P. G. Challenor (2012).  The effect of the nugget on Gaussian process emulators of computer models.  Computational Statistics & Data Analysis 56, 4215–4228.
  * Apted (1977) Apted, J. P. (1977).  Effects of weathering on some geotechnical properties of London clay.  PhD thesis, Imperial College (University of London).
  * Armstrong et al. (2024) Armstrong, J., P. Helm, J. Preston, and F. Loveridge (2024).  Economics of geotechnical asset deterioration, maintenance and renewal.  Transportation Geotechnics, 101185.
  * Bastos and O’Hagan (2009) Bastos, L. S. and A. O’Hagan (2009).  Diagnostics for Gaussian process emulators.  Technometrics 51, 425–438.
  * Briggs et al. (2023) Briggs, K. M., P. R. Helm, J. A. Smethurst, A. Smith, R. Stirling, A. Svalova, Y. Trinidad González, F. A. Loveridge, and S. Glendinning (2023).  Evidence for the weather-driven deterioration of ageing transportation earthworks in the uk.  Transportation Geotechnics 43, 101130.
  * Bromhead and Dixon (1986) Bromhead, E. N. and N. Dixon (1986).  The field residual strength of London Clay and its correlation with laboratory measurements, especially ring shear tests.  Géotechnique 36, 449–452.
  * BSI (2023) BSI (2023).  BS EN 1990: Eurocode. Basis of structural and geotechnical design.  OCLC: 1400083613.
  * Burland et al. (2023) Burland, J. B., T. Chapman, and M. Walsh (2023).  Foundations and other geotechnical elements in context–their role.  In ICE Manual of Geotechnical Engineering, Second edition, Volume I: Geotechnical engineering principles, problematic soils and site investigation, pp. 5–12. Emerald Publishing Limited.
  * Carpenter et al. (2017) Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. A. Brubaker, J. Guo, P. Li, and A. Riddell (2017).  Stan: A probabilistic programming language.  Journal of statistical software 76.
  * Conti et al. (2009) Conti, S., J. P. Gosling, J. E. Oakley, and A. O’Hagan (2009).  Gaussian process emulation of dynamic computer codes.  Biometrika 3, 663–676.
  * De Boor (1978) De Boor, C. (1978).  A practical guide to splines.  New York: Springer-Verlag.
  * Dixon et al. (2019) Dixon, N., C. J. Crosby, R. Stirling, P. N. Hughes, J. Smethurst, K. Briggs, D. Hughes, D. Gunn, P. Hobbs, F. Loveridge, S. Glendinning, T. Dijkstra, and A. Hudson (2019).  In situ measurements of near-surface hydraulic conductivity in engineered clay slopes.  Quarterly Journal of Engineering Geology and Hydrogeology 52, 123–135.
  * Duncan (1996) Duncan, J. M. (1996).  State of the art: Limit equilibrium and finite-element analysis of slopes.  Journal of Geotechnical Engineering 122, 577–596.
  * Ellis and O’Brien (2007) Ellis, E. A. and A. S. O’Brien (2007).  Effect of height on delayed collapse of cuttings in stiff clay.  Proceedings of the Institution of Civil Engineers- Geotechnical Engineering 160, 73–84.
  * Farah et al. (2014) Farah, M., P. Birrell, S. Conti, and D. De Angelis (2014).  Bayesian emulation and calibration of a dynamic epidemic model for A/H1N1 influenza.  Journal of the American Statistical Association 109, 1398–1411.
  * Farah and Kottas (2014) Farah, M. and A. Kottas (2014).  Bayesian inference for sensitivity analysis of computer simulations, with an application to radiative transfer models.  Technometrics 56, 159–173.
  * Fricker et al. (2013) Fricker, T. E., J. E. Oakley, and N. M. Urban (2013).  Multivariate Gaussian process emulators with nonseparable covariance structures.  Technometrics 55.
  * Gelman and Vehtari (2021) Gelman, A. and A. Vehtari (2021).  What are the most important statistical ideas of the past 50 years?  Journal of the American Statistical Association 116, 2087–2097.
  * Gneiting et al. (2007) Gneiting, T., F. Balabdaoui, and A. E. Raftery (2007).  Probabilistic forecasts, calibration and sharpness.  Journal of the Royal Statistical Society Series B: Statistical Methodology 69(2), 243–268.
  * Gramacy (2020) Gramacy, R. (2020).  Surrogates.  Boca Raton: CRC Press.
  * Helm et al. (2024a) Helm, P., A. Svalova, A. Morsy, M. Rouainia, A. Smith, A. El-Hamalawi, D. Wilkinson, H. Postill, and S. Glendinning (2024a).  Emulating long-term weather-driven transportation earthworks deterioration models to support asset management.  Transportation Geotechnics 44, 101155.
  * Helm et al. (2024b) Helm, P. R., A. Svalova, A. M. Morsy, M. Rouainia, A. Smith, A. El-Hamalawi, D. J. Wilkinson, H. Postill, and S. Glendinning (2024b).  Emulating long-term weather-driven transportation earthworks deterioration models to support asset management.  Transportation Geotechnics 44, 101155.
  * Henderson et al. (2009) Henderson, D. A., R. J. Boys, K. J. Krishnan, C. Lawless, and D. J. Wilkinson (2009).  Bayesian emulation and calibration of a stochastic computer model of mitochondrial DNA deletions in substantia nigra neurons.  Journal of the American Statistical Association 104, 76–87.
  * Homma and Saltelli (1996) Homma, T. and A. Saltelli (1996).  Importance measures in global sensitivity analysis of nonlinear models.  Reliability Engineering and System Safety 52, 1–17.
  * Itasca (2016) Itasca (2016).  Fast Lagrangian Analysis of Continua, Ver. 8.0.  Minneapolis: Itasca: Itasca Consulting Group, Inc.
  * Jette and Wickberg (2023) Jette, M. A. and T. Wickberg (2023).  Architecture of the slurm workload manager.  In D. Klusáček, J. Corbalán, and G. Rodrigo (Eds.), Job Scheduling Strategies for Parallel Processing. JSSPP 2023, Lecture Notes in Computer Science, Cham. Springer.
  * Jurečič et al. (2013) Jurečič, N., L. Zdravković, and V. Jovičić (2013).  Predicting ground movements in London Clay.  Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 166, 466–482.
  * Kennedy and O’Hagan (2001) Kennedy, M. C. and A. O’Hagan (2001).  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63(3), 425–464.
  * Knott (2000) Knott, G. D. (2000).  Interpolating cubic splines.  Boston: Springer.
  * Mellor et al. (2017) Mellor, R., L. Parry, and C. Power (2017, August).  CP6 earthworks asset policy developemnt. Task 63 - global stability and resilience appraisal. Interim report.  Technical report, Network Rail.
  * Mohammadi et al. (2019) Mohammadi, H., P. Challenor, and M. Goodfellow (2019).  Emulating dynamic non-linear simulators using Gaussian processes.  Computational Statistics and Data Analysis 139, 178–196.
  * Nowak (2012) Nowak, P. (2012).  Slope stabilisation methods.  In J. B. Burland, T. Chapman, H. Skinner, and M. Brown (Eds.), ICE Manual of geotechnical engineering. 2: Geotechnical design, construction and verification, Volume 2. London: ICE Publishing.
  * Oakley and Youngman (2017) Oakley, J. E. and B. D. Youngman (2017).  Calibration of stochastic computer simulators using likelihood emulation.  Technometrics 59(1), 80–92.
  * O’Hagan (2006) O’Hagan, A. (2006).  Bayesian analysis of computer code outputs: A tutorial.  Reliability Engineering & System Safety 91, 1290–1300.
  * Oyebamiji et al. (2017) Oyebamiji, O. K., D. J. Wilkinson, P. G. Jayathilake, T. P. Curtis, S. P. Rushton, B. Li, and P. Gupta (2017).  Gaussian process emulation of an individual-based model simulation of microbial communities.  Journal of Computational Science 22, 69–84.
  * Oyebamiji et al. (2019) Oyebamiji, O. K., D. J. Wilkinson, B. Li, P. G. Jayathilake, P. Zuliani, and C. T. P. (2019).  Bayesian emulation and calibration of an individual-based model of microbial communities.  Journal of Computational Science 30, 194–208.
  * Perry et al. (2003) Perry, J., M. Pedley, and K. Brady (2003).  Infrastructure cuttings - condition appraisal and remedial treatment.  Technical Report C591, CIRIA, London.
  * Postill et al. (2020) Postill, H., N. Dixon, G. Fowmes, A. El-Hamalawi, and W. A. Take (2020).  Modelling seasonal ratcheting and progressive failure in clay slopes: a validation.  Canadian Geotechnical Journal 57(9), 1265–1279.
  * Postill et al. (2023) Postill, H., P. R. Helm, N. Dixon, A. El-Hamalawi, S. Glendinning, and W. A. Take (2023).  Strength parameter selection framework for evaluating the design life of clay cut slopes.  Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 176(3), 254–273.
  * Postill et al. (2021) Postill, H., P. R. Helm, N. Dixon, S. Glendinning, J. A. Smethurst, M. Rouainia, K. M. Briggs, A. El-Hamalawi, and A. P. Blake (2021).  Forecasting the long-term deterioration of a cut slope in high-plasticity clay using a numerical model.  Engineering Geology 280, 1–19.
  * Potts et al. (1997) Potts, D. M., N. Kovacevic, and P. R. Vaughan (1997).  Delayed collapse of cut slopes in stiff clay.  Géotechnique 47, 953–982.
  * Rouainia et al. (2020) Rouainia, M., P. Helm, O. Davies, and S. Glendinning (2020).  Deterioration of an infrastructure cutting subjected to climate change.  Acta Geotechnica 15, 2997–3016.
  * Santner (2018) Santner, T. J. (2018).  The Design and Analysis of Computer Experiments (2 ed.).  Springer Series in Statistics.
  * Stan Development Team (2022) Stan Development Team (2022).  Stan Modeling Language Users Guide and Reference Manual (VERSION 2.29 ed.).  https://mc-stan.org.
  * Summersgill et al. (2018) Summersgill, F. C., S. Kontoe, and D. M. Potts (2018).  Stabilisation of excavated slopes in strain-softening materials with piles.  Géotechnique 68, 1–14.
  * Svalova et al. (2021) Svalova, A., P. Helm, D. Prangle, M. Rouainia, S. Glendinning, and D. J. Wilkinson (2021).  Emulating computer experiments of transport infrastructure slope stability using Gaussian processes and Bayesian inference.  Data-Centric Engineering 2, e12.
  * Talebitooti et al. (2015) Talebitooti, R., M. H. Shojaeefard, and S. Yarmohammadisatri (2015).  Shape design optimization of cylindrical tank using B-spline curves.  Computers & Fluids 109, 100–112.
  * Tsiampousi et al. (2017) Tsiampousi, A., L. Zdravkovic, and D. Potts (2017).  Numerical study of the effect of soil–atmosphere interaction on the stability and serviceability of cut slopes in London Clay.  Canadian Geotechnical Journal 54, 405–418.
  * van Genuchten (1980) van Genuchten, M. T. (1980).  A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils.  Soil Science Society of America Journal 44, 892–898.

Supplementary Materials

##  Appendix A B-spline definition

Consider a vector of knots $\bm{k}=(k_{0},\ldots,k_{m+1})$ over which the
B-spline is to be fitted. The boundary knots are equivalent to zero-time
($k_{0}=0$) and the failure time $k_{m+1}=\omega$, and there are $m$ interior
knots. As in the quadratic model, we assume that $\omega$ is unobserved and
will be estimated with the remaining parameters. To fit a B-spline of order
$o$ (degree $o-1$) over $\bm{k}$, define an augmented knot sequence
$\bm{k}^{*}$ with $m+2o$ elements

$k_{-o+1}=\cdots=k_{0}\leq\cdots\leq k_{m}\leq k_{m+1}=\cdots=k_{m+o}.$ (19)

It is noteworthy that in Equation (19) the first and the last $o$ knots are
repeated to fix the boundary constraints. For easier computation, the
subscripts of the elements of $\bm{k}^{*}$ can be re-labelled as
$k_{l}\subset\bm{k}^{*},l\in\\{0,\ldots,m+2o-1\\}$. We assume that two
piecewise polynomials could be sufficient to explain the two dominant failure
modes (deep and shallow), thus $m=1$ and $o=3$. This gives the augmented knot
sequence
$\bm{k}^{*}=(k_{0},k_{0},k_{0},k_{1},k_{2},k_{2},k_{2})=(0,0,0,k,\omega,\omega,\omega)$.
Frequently, the interior knots are equidistant, and here we assume that
$k=0.5\omega$.
The B-spline is defined using the following recursive relationships. Firstly,
we define a set of basis (indicator) functions $\phi_{l,j}$ for
$j\in\\{0,1,\ldots,o-1\\}$ and $l\in\\{0,1,\ldots,m+2o-2\\}$:

$\phi_{l,0}(x)=\begin{cases}1,&\text{if }k_{l}\leq x<k_{l+1},\\\
0,&\text{otherwise}.\end{cases}$ (20)

Then, for $j\in\\{0,1,\ldots,o-1\\}$ and $l\in\\{0,1,\ldots,m+2o-2-(j+1)\\}$
we define the set of functions:

$\phi_{l,j+1}(x)=\alpha_{l,j+1}(x)\phi_{l,j}(x)+[1-\alpha_{l+1,j+1}(x)]\phi_{l+1,j}(x),$
(21)

where

$\alpha_{l,j}(x)=\begin{cases}\frac{x-k_{l}}{k_{l+j}-k_{l}},&\text{if
}k_{l+j}\neq k_{l},\\\ 0,&\text{otherwise}.\end{cases}$ (22)

The B-spline approximation of $Y_{i,j}=\,$FoS${}_{i,j}-1$ can be written as

$Y_{i,j}=\sum_{l=0}^{m+o-1}\gamma_{l,i}\phi_{l,o-1}(t_{i,j})+\varepsilon_{i,j}=\sum_{l=0}^{3}\gamma_{l,i}\phi_{l,2}(t_{i,j})+\varepsilon_{i,j},\quad\varepsilon_{i,j}\sim\text{N}(0,\sigma^{2}_{i}).$
(23)

Expanding Equation (23) using Equations (20), (21) and (22) we find that the
spline model can be written as

$\displaystyle Y_{i,j}=$
$\displaystyle\bigg{\\{}\gamma_{0,i}+t_{i,j}\bigg{(}\frac{4\gamma_{1,i}}{\omega_{i}}-\frac{4\gamma_{0,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{4\gamma_{0,i}}{\omega_{i}^{2}}-\frac{6\gamma_{1,i}}{\omega_{i}^{2}}+\frac{2\gamma_{2,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}I(0\leq
t_{i,j}<\omega_{i}/2)+$ (24)
$\displaystyle\bigg{\\{}2(\gamma_{1,i}-\gamma_{2,i}+\gamma_{3,i})+t_{i,j}\bigg{(}\frac{8\gamma_{2,i}}{\omega_{i}}-\frac{4\gamma_{1,i}}{\omega_{i}}-\frac{4\gamma_{3,i}}{\omega_{i}}\bigg{)}+$
$\displaystyle
t_{i,j}^{2}\bigg{(}\frac{2\gamma_{1,i}}{\omega_{i}^{2}}-\frac{6\gamma_{2,i}}{\omega_{i}^{2}}+\frac{4\gamma_{3,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}I(\omega_{i}/2\leq
t_{i,j}<\omega_{i})+\varepsilon_{i,j},$

where $\varepsilon_{i,j}\sim\text{N}(0,\sigma^{2}_{i})$. Conveniently,
removing $\gamma_{3,i}$ forces the B-spline curve to cross zero at the last
knot, which is well-suited for modelling the FoS. Therefore, the B-spline that
will be used in this study is

$\displaystyle Y_{i,j}$
$\displaystyle=\bigg{\\{}\gamma_{0,i}+t_{i,j}\bigg{(}\frac{4\gamma_{1,i}}{\omega_{i}}-\frac{4\gamma_{0,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{4\gamma_{0,i}}{\omega_{i}^{2}}-\frac{6\gamma_{1,i}}{\omega_{i}^{2}}+\frac{2\gamma_{2,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}I(0\leq
t_{i,j}<\omega_{i}/2)+$ (25)
$\displaystyle\bigg{\\{}2(\gamma_{1,i}-\gamma_{2,i})+t_{i,j}\bigg{(}\frac{8\gamma_{2,i}}{\omega_{i}}-\frac{4\gamma_{1,i}}{\omega_{i}}\bigg{)}+t_{i,j}^{2}\bigg{(}\frac{2\gamma_{1,i}}{\omega_{i}^{2}}-\frac{6\gamma_{2,i}}{\omega_{i}^{2}}\bigg{)}\bigg{\\}}I(\omega_{i}/2\leq
t_{i,j}<\omega_{i})+\varepsilon_{i,j}.$

##  Appendix B Posterior diagnostic plots

Figure B1: Examples of posterior trace plots of the quadratic model
parameters. Figure B2: Examples of posterior trace plots of the B-spline model
parameters. Figure B3: Examples of multi-modal posterior trace plots of the
B-spline model parameters corresponding to the FoS plots shown in Figure B4.
Figure B4: Examples of multi-modal posterior plots of FoS. Figure B5:
Posterior trace plots of the emulator parameters using the quadratic model.
Figure B6: Posterior trace plots of the emulator parameters using the B-spline
model.

##  Appendix C Posterior plots of FoS

Figure C1: Examples of posterior distributions of FoS using the quadratic
model (blue) and the B-spline model (red). The dashed lines represent the
posterior means and the solid lines represent the central posterior 95%
prediction intervals. Figure C2: Examples of posterior distributions of FoS
using the quadratic model (blue) and the B-spline model (red). The dashed
lines represent the posterior means and the solid lines represent the central
posterior 95% prediction intervals. Figure C3: Examples of posterior
distributions of FoS using the quadratic model (blue) and the B-spline model
(red). The dashed lines represent the posterior means and the solid lines
represent the central posterior 95% prediction intervals. Figure C4: Examples
of posterior distributions of FoS using the quadratic model (blue) and the
B-spline model (red). The dashed lines represent the posterior means and the
solid lines represent the central posterior 95% prediction intervals. Figure
C5: Examples of posterior distributions of FoS using the quadratic model
(blue) and the B-spline model (red). The dashed lines represent the posterior
means and the solid lines represent the central posterior 95% prediction
intervals. Figure C6: Examples of posterior distributions of FoS using the
quadratic model (blue) and the B-spline model (red). The dashed lines
represent the posterior means and the solid lines represent the central
posterior 95% prediction intervals. Figure C7: Examples of posterior
distributions of FoS using the quadratic model (blue) and the B-spline model
(red). The dashed lines represent the posterior means and the solid lines
represent the central posterior 95% prediction intervals. Figure C8: Examples
of posterior distributions of FoS using the quadratic model (blue) and the
B-spline model (red). The dashed lines represent the posterior means and the
solid lines represent the central posterior 95% prediction intervals.

##  Appendix D Diagnostics of posterior fit

Figure D1: Plots (left) presents the difference in the $MSE$ between the
quadratic model and the B-spline model for sixty-six within-sample computer
runs. Plot (right) presents the difference in the $crps$ between the quadratic
model and the B-spline model for sixty six within-sample computer runs. The
numbers on the x-axis correspond to run numbers.

##  Appendix E Posterior distributions of predicted TTF

Figure E1: Examples of posterior distributions of predicted TTF, $\rho_{i}$,
using the quadratic model (blue) and the B-spline model (red). The dashed
lines represent the central posterior 95% prediction intervals under each
model. The solid red lines indicate true TTF for computer runs that reached
failure. Figure E2: Examples of posterior distributions of predicted TTF,
$\rho_{i}$, using the quadratic model (blue) and the B-spline model (red). The
dashed lines represent the central posterior 95% prediction intervals under
each model. The solid red lines indicate true TTF for computer runs that
reached failure. Figure E3: Examples of posterior distributions of predicted
TTF, $\rho_{i}$, using the quadratic model (blue) and the B-spline model
(red). The dashed lines represent the central posterior 95% prediction
intervals under each model. The solid red lines indicate true TTF for computer
runs that reached failure. Figure E4: Examples of posterior distributions of
predicted TTF, $\rho_{i}$, using the quadratic model (blue) and the B-spline
model (red). The dashed lines represent the central posterior 95% prediction
intervals under each model. The solid red lines indicate true TTF for computer
runs that reached failure. Figure E5: Examples of posterior distributions of
predicted TTF, $\rho_{i}$, using the quadratic model (blue) and the B-spline
model (red). The dashed lines represent the central posterior 95% prediction
intervals under each model. The solid red lines indicate true TTF for computer
runs that reached failure. Figure E6: Examples of posterior distributions of
predicted TTF, $\rho_{i}$, using the quadratic model (blue) and the B-spline
model (red). The dashed lines represent the central posterior 95% prediction
intervals under each model. The solid red lines indicate true TTF for computer
runs that reached failure. Figure E7: Examples of posterior distributions of
predicted TTF, $\rho_{i}$, using the quadratic model (blue) and the B-spline
model (red). The dashed lines represent the central posterior 95% prediction
intervals under each model. The solid red lines indicate true TTF for computer
runs that reached failure. Figure E8: Examples of posterior distributions of
predicted TTF, $\rho_{i}$, using the quadratic model (blue) and the B-spline
model (red). The dashed lines represent the central posterior 95% prediction
intervals under each model. The solid red lines indicate true TTF for computer
runs that reached failure. Figure E9: Examples of deterministic (noise-free)
FoS curves with different internal knot positions.

##  Appendix F Posterior distributions of FoS for the validation data sets

Figure F1: Examples of posterior distributions of FoS using the quadratic
model (blue) and the B-spline model (red) for four out-of sample computer
runs. The dashed lines represent the posterior means and the solid lines
represent the central posterior 95% prediction intervals. Figure F2: Examples
of posterior distributions of FoS using the quadratic model (blue) and the
B-spline model (red) for four out-of sample computer runs. The dashed lines
represent the posterior means and the solid lines represent the central
posterior 95% prediction intervals.