Multifidelity Monte Carlo estimation for efficient uncertainty quantification in climate-related modeling

Uncertainties in an output of interest that depends on the solution of a complex system (e.g., of partial differential equations with random inputs) are often, if not nearly ubiquitously, determined in practice using Monte Carlo (MC) estimation. While simple to implement, MC estimation fails to provide reliable information about statistical quantities (such as the expected value of the output of interest) in application settings such as climate modeling, for which obtaining a single realization of the output of interest is a costly endeavor. Specifically, the dilemma encountered is that many samples of the output of interest have to be collected in order to obtain an MC estimator that has sufficient accuracy – so many, in fact, that the available computational budget is not large enough to effect the number of samples needed. To circumvent this dilemma, we consider using multifidelity Monte Carlo (MFMC) estimation which leverages the use of less costly and less accurate surrogate models (such as coarser grids, reduced-order models, simplified physics, and/or interpolants) to achieve, for the same computational budget, higher accuracy compared to that obtained by an MC estimator – or, looking at it another way, an MFMC estimator obtains the same accuracy as the MC estimator at lower computational cost. The key to the efficacy of MFMC estimation is the fact that most of the required computational budget is loaded onto the less costly surrogate models so that very few samples are taken of the more expensive model of interest. We first provide a more detailed discussion about the need to consider an alternative to MC estimation for uncertainty quantification. Subsequently, we present a review, in an abstract setting, of the MFMC approach along with its application to three climate-related benchmark problems as a proof-of-concept exercise.

1 Introduction

In many application settings – climate modeling being a prominent one – large computational costs are incurred when solutions to a given model are approximated to within an acceptable accuracy tolerance. In fact, this cost can be prohibitively large when one has to obtain the results of multiple simulations, as is the case for, e.g., uncertainty quantification, control, and optimization, to name a few. Thus, there is often a need for compromise between the accuracy of simulation algorithms and the number of simulations needed to obtain, say, in the uncertainty quantification setting, accurate statistical information.

For example, consider the following case, which represents the focus of this paper. Suppose one has a complex system, say, a system of discretized partial differential equations, for which the input data depends on a vector of randomly distributed parameters z. Letting u(z) denote the solution of the discretized partial differential equations, we define an output of interest F(u(z))=F(z) that depends on u(z) so that, of course, it also depends on the choice of z. Next, suppose we wish to use Monte Carlo sampling to estimate the expected value Q=𝔼(F(z)) of the output of interest F(z), i.e., we have the Monte Carlo estimator of Q given by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{Q}^{\mathrm{MC}}& ={\displaystyle \frac{\mathrm{1}}{M}}\sum _{m=\mathrm{1}}^{M}F\left({\mathit{z}}_m\right)={\displaystyle \frac{\mathrm{1}}{M}}\sum _{m=\mathrm{1}}^{M}F\left(\mathit{u}\right({\mathit{z}}_m\left)\right)\approx \mathbb{E}\left(F\left(\mathit{z}\right)\right)\\ & =Q,\end{array}\end{array}$$

where $\mathit{\{}{\mathit{z}}_m{\mathit{\}}}_{m=\mathrm{1}}^M$ denotes a set of M independent and identically distributed samples of the random vector z. This is a common task in uncertainty quantification, which provides useful statistical information in both predictive and inferential applications.

On the other hand, Monte Carlo estimation is not without its drawbacks. Consider the following scenario. Let δ<1 denote a measure of the spatial grid size used to discretize the partial differential equation (PDE) system in question, and suppose that δ is normalized by the diameter of the domain in which that system is posed. Then, we have it that the error in the approximate solution is of $\mathcal{O}\left({\mathit{\delta }}^{{\mathit{\nu }}_{\mathrm{d}}}\right)$ for some ν_d>0 and that the cost (e.g., measured in seconds, days, or weeks) of a simulation which obtains the approximate solution is of $\mathcal{O}\left({\mathit{\delta }}^{-{\mathit{\nu }}_{\mathrm{c}}}\right)$ for some ν_c>0. Neglecting any additional costs connected with the evaluation of F(z)=F(u(z)) (for any given z) once u(z) is obtained, it is clear that the cost C incurred when obtaining the Monte Carlo estimator Q^MC is of $\mathcal{O}\left(M{\mathit{\delta }}^{-{\mathit{\nu }}_{\mathrm{c}}}\right)$. Of course, it is well known that the accuracy of a Monte Carlo estimator is of $\mathcal{O}(\mathrm{1}/\sqrt{M})$ so that in order for the accuracy of the estimator to be commensurate with the discretization error, we must have that $\mathrm{1}/\sqrt{M}\propto \mathcal{O}\left({\mathit{\delta }}^{{\mathit{\nu }}_{\mathrm{d}}}\right)$. Therefore, the needed number of samples is proportional to $M\propto {\mathit{\delta }}^{-\mathrm{2}{\mathit{\nu }}_{\mathrm{d}}}$, and the total cost CM incurred when determining the Monte Carlo estimator Q^MC is of $\mathcal{O}\left({\mathit{\delta }}^{-\mathrm{2}{\mathit{\nu }}_{\mathrm{d}}-{\mathit{\nu }}_{\mathrm{c}}}\right)$.

The cost CM can quickly get out of control when dealing with large-scale problems. For example, suppose that the approximate solution of the PDE system is second-order accurate, i.e., ν_d=2, and that a single simulation incurs a cost of 𝒪(δ⁻³) (i.e., ν_c=3), which is the best-case scenario in three dimensions. We then have that the number of samples needed is $M\propto {\mathit{\delta }}^{-\mathrm{4}}$, and the total cost CM is $\propto {\mathit{\delta }}^{-\mathrm{7}}$. Thus, with even a modest grid size of δ=0.01, we have it that the number of samples needed is M∝10⁸ and that the total cost is CM ∝10¹⁴, which is too high for practical use. Turning things around, suppose instead that the available computational budget allows for at most 10 thousand simulations, i.e., M=10⁴, so that the resulting accuracy of the Monte Carlo simulator is of 𝒪(10⁻²). Here, the estimator error and approximation error are already commensurate when δ=0.1, making the choice of a smaller delta redundant if not infeasible. In fact, in this case, there is no way to achieve the four digits of accuracy sought by choosing δ=0.01, since the best we can do with the available budget is on the order of 𝒪(10⁻²). Note that the situation becomes even worse when the PDE system in question is time dependent, as a single simulation incurs an even larger cost when one accounts for the number of time steps used in the simulation.

Given that the cost of MC estimation is at times prohibitively expensive, it comes as no surprise that many alternatives or run-arounds to such estimation have been proposed. One approach in this direction has led to the development of many different random-parameter sampling schemes (e.g., quasi-Monte Carlo sampling, sparse-grid sampling, importance sampling, Latin hypercube sampling, lattice sampling, and compressed sensing, to name just a few) for which the estimation error is guaranteed to be smaller than its Monte Carlo equivalent; see, e.g., Addcock et al. (2022), Evans and Swartz (2000), Gunzburger et al. (2014), Nierderreiter (1992), Sloan (1994), and Smith (2013). On the other hand, some of these alternate approaches require smoothness of solutions to achieve better accuracy. Furthermore, most, if not all, of these methods are superior to Monte Carlo sampling only for the moderate dimension of the parameter vector z; see the references just cited.

A second approach towards reducing the cost of Monte Carlo (and, for that matter, for any type of) uncertainty estimation is to use approximate solutions of the PDE system that are less costly to obtain compared to the cost of obtaining the approximation of actual interest. For example, using simulations obtained using coarser grids or using reduced-order models such as reduced-basis or proper-orthogonal-decomposition methods are less costly, as are interpolation and support vector machine approximations; see, e.g., Cristianini and Shawe-Taylor (2000), Fritzen and Ryckelynck (2019), Keiper et al. (2018), Quarteroni and Rozza (2014), Quarteroni et al. (2016), and Steinwart and Christmann (2008). However, such approaches, by definition, result in less accurate approximations compared to the accuracy that one wants to achieve.

In this paper, we do not consider any of the possible alternate sampling schemes, nor do we exclusively consider using less costly and less accurate approximate solutions of the PDE system. Instead, because of the near ubiquity of its use in practice, our goal is to outperform traditional Monte Carlo estimation by using a nontraditional Monte Carlo sampling strategy and, in so doing, to refrain from incurring any loss of accuracy. To meet this goal, we invoke multifidelity Monte Carlo estimation which, in addition to the expensive and accurate PDE system approximation of interest (hereafter referred to as the “truth” approximation), also uses cheaper-to-obtain and less accurate approximations (which are referred to as the “surrogates”). The bottom line is that multifidelity Monte Carlo estimation meets our goal by leveraging increased sampling of the less accurate and less costly approximations alongside low sampling of the more expensive and more accurate truth approximation. The multifidelity Monte Carlo algorithm systematically determines the number of samples taken from each surrogate (i.e., there is no guess-work involved) and systematically (i.e., again, there is no guess-work involved) combines the samples of the surrogates to obtain the desired estimator. We note that multifidelity Monte Carlo estimation has already been shown to outperform Monte Carlo estimation in a variety of application settings; see, e.g., Clare et al. (2022), Dimarco et al. (2023), Konrad (2019), Law et al. (2022), Modderman (2021), Quick et al. (2019), Rezaeiravesh et al. (2020), Romer et al. (2020), and Valero et al. (2022) for some examples.

In Sect. 2, we review how the Monte Carlo and multifidelity Monte Carlo estimators are constructed in an abstract setting; here we follow the expositions in Peherstorfer et al. (2016) and also in Gruber et al. (2022a). Then, in Sects. 3 and 4, we demonstrate the effectiveness of multifidelity Monte Carlo estimation using three well-known climate-related benchmarks. Section 3.1 and 3.2 are respectively devoted to single-layer shallow-water equation models based on the SOMA test case of Wolfram et al. (2015) and Test Case 5 of Williamson et al. (1992). Section 4 is devoted to a benchmark case for the first-order model for ice sheets; see Blatter (1995), Pattyn (2003), and Perego et al. (2012), Tezaur et al. (2015). We close by providing some concluding remarks in Sect. 5.

2 Monte Carlo and multifidelity Monte Carlo estimators

An abstraction of the specific settings considered in Sects. 3 and 4 is given as follows:

• It involves having in hand a (discretized) partial differential equation (PDE) system for which the solution u(z) depends on a random vector of parameters z∈Γ, where Γ denotes a parameter domain.

  Note that the input data to this PDE system (e.g., forcing terms, initial conditions, and coefficients) could depend on one or more of the components of the random vector z.

• We are interested in situations whereby, for any z∈Γ, obtaining u(z) is a costly endeavor.

• We define a scalar output of interest (OoI) F¹(z)=F¹(u(z)) that depends on the solution u(z) of the (discretized) PDE system; of course, if obtaining u(z) is costly, then so is obtaining F¹(z).

  While they are not considered in this work, vector-valued outputs of interest can also be treated with multifidelity Monte Carlo techniques. OoIs could be, e.g., averages or extremal values of the energy associated with the solution u(z).

• Having defined an OoI F¹(z), we let Q₁=𝔼(F¹) denote the quantity of interest (QoI) corresponding to the F¹(z), where 𝔼( ⋅ ) denotes the expected value with respect to z.

• Because the estimation of Q₁=𝔼(F¹) is the central goal of this paper, we refer to F¹(z) as the “truth” output of interest.

• Commonly, even ubiquitously, a Monte Carlo (MC) sampling method is used to (approximately) quantify the uncertainty in the chosen OoI F¹(z).

• Specifically, MC sampling is used to estimate the quantity of interest Q₁=𝔼(F¹) corresponding to F¹(z), i.e., we have it that the MC estimator of the QoI Q₁=𝔼(F¹) is given by

  $$\begin{array}{}\text{(2)}& Q_{\mathrm{1}}^{\mathrm{MC}}={\displaystyle \frac{\mathrm{1}}{M_{\mathrm{1}}}}\sum _{m=\mathrm{1}}^{M_{\mathrm{1}}}{F}^{\mathrm{1}}\left({\mathit{z}}_m\right)\approx Q_{\mathrm{1}}=\mathbb{E}\left(F^{\mathrm{1}}\right),\end{array}$$

  where $\mathit{\{}{\mathit{z}}_m{\mathit{\}}}_{m=\mathrm{1}}^{M_{\mathrm{1}}}$ denotes a set of M₁ randomly selected points in the parameter domain Γ.

• We are then faced with the following dilemma: on the one hand, obtaining an acceptably accurate MC estimator $Q_{\mathrm{1}}^{\mathrm{MC}}$ of the QoI Q₁ requires obtaining the solution u(z) of the discretized PDE system at many randomly selected points z∈Γ; on the other hand, each of those approximate solutions u(z) are computationally costly to obtain.

  Quantifying uncertainties in climate system settings are victimized by this two-headed dilemma to the extent that, e.g., accurate long-time integrations often cannot be realized in practice.

Due to the issue of prohibitive computational cost, we turn to multifidelity Monte Carlo (MFMC) methods for uncertainty quantification. MFMC methods leverage the availability of surrogate outputs of interest F^k(z), $k=\mathrm{2},\mathrm{\dots },K$, which have smaller computational complexity compared to that of the truth OoI F¹(z). As mentioned in Sect. 1, there are many types of surrogates that can be used for this purpose, e.g., discretized PDEs with coarser spatial grids and larger time steps, reduced-basis and proper-orthogonal-decomposition-based reduced-order models, and interpolants of F¹(z), to name just a few.

For the truth and the K−1 surrogate OoIs, Monte Carlo (MC) sampling yields the K unbiased estimators

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{Q}_k^{\mathrm{MC}}& ={\displaystyle \frac{\mathrm{1}}{M_k}}\sum _{m=\mathrm{1}}^{M_k}{F}^k\left({\mathit{z}}_m\right)\approx \mathbb{E}\left(F^k\right)\\ & =Q_k\text{for}k=\mathrm{1},\mathrm{\dots },K,\end{array}\end{array}$$

where ${\mathit{z}}_{\mathrm{1}},\mathrm{\dots },{\mathit{z}}_{M_k}$ denote M_k independent and identically distributed (i.i.d.) samples of z∈Γ. The goal of an MFMC method is to use an appropriate linear combination of the K MC estimators $Q_k^{\mathrm{MC}}$ in Eq. (3) to estimate the QoI Q₁=𝔼(F¹). Specifically, we define the MFMC estimator as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{Q}^{\mathrm{MFMC}}& =Q_{\mathrm{1}}^{\mathrm{MC}}+\sum _{k=\mathrm{2}}^K{\mathit{\alpha }}_k\left(Q_k^{\mathrm{MC}}-Q_{k-\mathrm{1}}^{\mathrm{MC}}\right)\approx Q_{\mathrm{1}}\\ & =\mathbb{E}\left(F^{\mathrm{1}}\right),\end{array}\end{array}$$

where $\mathit{\{}{\mathit{\alpha }}_k{\mathit{\}}}_{k=\mathrm{2}}^K$ denotes a sequence of scalar weights, and $\mathit{\{}{M}_k{\mathit{\}}}_{k=\mathrm{1}}^K$ denotes a nondecreasing sequence of integers defining the number of samples.

Letting C_k denote the cost of evaluating the kth output of interest F^k, the costs of computing the respective MC and MFMC estimators are given as

$$\begin{array}{ll}{\displaystyle }{C}_k^{\mathrm{MC}}& {\displaystyle }=C_k{M}_k\text{for}k=\mathrm{1},\mathrm{\dots },K\text{and}\\ {\displaystyle }& {\displaystyle }{C}^{\mathrm{MFMC}}=\sum _{k=\mathrm{1}}^K{C}_k^{\mathrm{MC}}=\mathbf{C}\cdot \mathbf{M},\end{array}$$

where $\mathbf{C}=\mathit{\{}{C}_k{\mathit{\}}}_{k=\mathrm{1}}^K$ and $\mathbf{M}=\mathit{\{}{M}_k{\mathit{\}}}_{k=\mathrm{1}}^K$ denote K vectors formed from the previously introduced sequences.

The variance ${\mathit{\sigma }}_k^{\mathrm{2}}$ of the kth output of interest F^k and the correlation ${\mathit{\zeta }}_{k,k^{\prime }}$ between the outputs of interest F^k and $F^{k^{\prime }}$ are given by, for $k^{\prime },k=\mathrm{1},\mathrm{\dots },K$,

$${\mathit{\sigma }}_k^{\mathrm{2}}=\mathrm{Var}\left[F^k\left(\mathit{z}\right)\right]\text{and}{\mathit{\zeta }}_{k,k^{\prime }}={\displaystyle \frac{\mathrm{Cov}\left[F^k\left(\mathit{z}\right),F^{k^{\prime }}\left(\mathit{z}\right)\right]}{{\mathit{\sigma }}_k{\mathit{\sigma }}_{k^{\prime }}}},$$

respectively, where $\mathrm{Cov}\left[\cdot ,\cdot \right]$ denotes the statistical covariance. We then have that the mean-squared error (MSE) incurred by the MC estimator $Q_k^{\mathrm{MC}}$ of the QoI Q_k is given by

$$\begin{array}{}\text{(5)}& \begin{array}{rl}e\left(Q_k^{\mathrm{MC}}\right)& =\mathbb{E}\left({\left(Q_k-Q_k^{\mathrm{MC}}\right)}^{\mathrm{2}}\right)\\ & =\mathbb{E}\left({\left(\mathbb{E}\left(F^k\left(\mathit{z}\right)\right)-Q_k^{\mathrm{MC}}\right)}^{\mathrm{2}}\right)\\ & =\mathbb{E}{\left(F^k\left(\mathit{z}\right)\right)}^{\mathrm{2}}-\mathbb{E}\left(F^k(\mathit{z}{)}^{\mathrm{2}}\right)\\ & ={\displaystyle \frac{{\mathit{\sigma }}_k^{\mathrm{2}}}{M_k}}\text{for }k=\mathrm{1},\mathrm{\dots },K,\end{array}\end{array}$$

whereas the MSE incurred by the Q^MFMC estimator of the QoI, Q₁ is given by (see Peherstorfer et al., 2016)

$$\begin{array}{}\text{(6)}& \begin{array}{rl}e\left(Q^{\mathrm{MFMC}}\right)& ={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}{M_{\mathrm{1}}}}+\sum _{k=\mathrm{2}}^K\left({\displaystyle \frac{\mathrm{1}}{M_{k-\mathrm{1}}}}-{\displaystyle \frac{\mathrm{1}}{M_k}}\right)\\ & \left({\mathit{\alpha }}_k^{\mathrm{2}}{\mathit{\sigma }}_k^{\mathrm{2}}-\mathrm{2}{\mathit{\alpha }}_k{\mathit{\sigma }}_k{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\zeta }}_{\mathrm{1},k}\right).\end{array}\end{array}$$

Note that it can be shown (see, e.g., Peherstorfer et al., 2016) that the MSE for Q^MFMC is lower than that for $Q_{\mathrm{1}}^{\mathrm{MC}}$ if and only if

$$\begin{array}{}\text{(7)}& \sqrt{{\displaystyle \frac{e\left(Q^{\mathrm{MFMC}}\right)}{e\left(Q_{\mathrm{1}}^{\mathrm{MC}}\right)}}}=\sum _{k=\mathrm{1}}^{K-\mathrm{1}}\sqrt{{\displaystyle \frac{C_k}{C_{\mathrm{1}}}}\left({\mathit{\zeta }}_{\mathrm{1},k}^{\mathrm{2}}-{\mathit{\zeta }}_{\mathrm{1},k+\mathrm{1}}^{\mathrm{2}}\right)}<\mathrm{1}.\end{array}$$

Given a fixed computational budget B, MFMC aims to construct an optimal sampling strategy $\mathbf{M}=\mathit{\{}{M}_k{\mathit{\}}}_{k=\mathrm{1}}^K$ along with an optimal set of weights $\mathit{\alpha }=\mathit{\{}{\mathit{\alpha }}_k{\mathit{\}}}_{k=\mathrm{2}}^K$ so that the MSE (Eq. 6) of the multifidelity estimator Q^MFMC is lower than the MSE (Eq. 5) of the Monte Carlo estimator $Q_{\mathrm{1}}^{\mathrm{MC}}$. Viewed differently, this means that an appropriate estimator Q^MFMC can achieve a fixed MSE ε>0 at a smaller computational cost compared to that incurred for the Monte Carlo estimator $Q_{\mathrm{1}}^{\mathrm{MC}}$ achieving MSE ε.

In Peherstorfer et al. (2016, 2018), unique optimal values of $\mathbf{M}=\mathit{\{}{M}_k{\mathit{\}}}_{k=\mathrm{1}}^K$ and $\mathit{\alpha }=\mathit{\{}{\mathit{\alpha }}_k{\mathit{\}}}_{k=\mathrm{2}}^K$ are analytically obtained by minimizing the MSE (Eq. 6) of the Q^MFMC estimator (Eq. 4) over the real numbers. However, the values M must certainly be integers for practical use, and the heuristic of Peherstorfer et al. (2016) which determines the optimal values can result in either a biased estimator of the expectation or (with naïve modification) a violation of the given budget B. For “small” computational budgets, the consequences of this can be quite severe, as illustrated in Gruber et al. (2022a).

Because “small” computational budgets are of high interest for climate modeling, here, instead of using the MFMC method of Peherstorfer et al. (2016, 2018), we use the modified MFMC method of Gruber et al. (2022a), which guarantees that the optimal sampling numbers are integers and that the computational budget is not exceeded. In that method, instead of simply minimizing the MSE e(Q^MFMC) defined in Eq. (6), the modified MFMC estimator is determined through the use of at least some of the sequential minimization problems given by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \text{for }k=\mathrm{1},\mathrm{2},\mathrm{\dots },K\text{and, if }k>\mathrm{1},\\ & \text{for given }{M}_{\mathrm{1},\mathrm{1}},\mathrm{\dots },M_{k-\mathrm{1},k-\mathrm{1}},\\ & \text{minimize the functional}\\ & {\mathcal{L}}^k({\mathbf{M}}_k,{\mathit{\alpha }}_k;{\mathit{\lambda }}_k,{\mathit{\mu }}_k,{\mathit{\xi }}_k)=\\ & e\left(Q^{\mathrm{MFMC}}\right)+{\mathit{\lambda }}_k(\sum _{k^{\prime }=k}^K{C}_{k^{\prime }}{M}_{k^{\prime },k}-(B-\sum _{k^{\prime }=\mathrm{1}}^{k-\mathrm{1}}{C}_{k^{\prime }}\left)\right)\\ & +\sum _{k^{\prime }=k+\mathrm{1}}^K{\mathit{\mu }}_{k^{\prime },k}(M_{k^{\prime },k}-M_{k^{\prime }-\mathrm{1},k})-{\mathit{\xi }}_k{M}_{k,k},\end{array}\end{array}$$

where ${\mathit{\mu }}_{k^{\prime },k}$ for $k^{\prime }=k+\mathrm{1},\mathrm{\dots },K$, λ_k, and ξ_k>0 are Lagrange multipliers.

In Eq. (8), the first term added to the MSE e(Q^MFMC) enforces the budget constraint – or, more precisely, for each k, enforces that the remaining available budget suffices to not exceed the given budget B after that budget has been depleted by the sampling already effected for the OoIs $F^{k^{\prime }}$, $k^{\prime }=\mathrm{1},\mathrm{\dots },k-\mathrm{1}$, where, of course, for k=1 the whole budget B is available. The second term added in Eq. (8) enforces the monotone non-decrease in the number of samples $M_{k^{\prime },k}$ for $k^{\prime }>k$. The third term added to the MSE enforces the positivity of M_k,k.

Repeating the arguments made in Peherstorfer et al. (2016) concerning optimal choices for the sample numbers and weights, the results given in Box 1 below are proved in Gruber et al. (2022a).

Box 1 Optimal non-integer sampling numbers for the modified MFMC method; see Gruber et al. (2022a).

It is remarkable that, due to Eq. (11), the weights ${\mathit{\alpha }}_{k^{\prime },k}^{*}$ depend only on input data, specifically on the variances and correlations of the OoIs $\mathit{\{}{F}^{k^{\prime }}{\mathit{\}}}_{k^{\prime }=\mathrm{1}}^K$. As such, the value of ${\mathit{\alpha }}_{k^{\prime },k}^{*}$ is independent of k so that

$$\begin{array}{rl}{\mathit{\alpha }}_{k^{\prime },k}^{*}={\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{1},k^{\prime }}{\mathit{\sigma }}_{\mathrm{1}}}{{\mathit{\sigma }}_{k^{\prime }}}}={\mathit{\alpha }}_{k^{\prime },\mathrm{1}}^{*}& \text{for all }k=\mathrm{1},\mathrm{\dots },K\\ & \text{ and all }{k}^{\prime }=k+\mathrm{1},\mathrm{\dots },K.\end{array}$$

Thus, the set of weights $\mathit{\{}{\mathit{\alpha }}_{k^{\prime },\mathrm{1}}^{*}{\mathit{\}}}_{k^{\prime }=\mathrm{2}}^K$ determined from the minimization of ℒ¹ suffices to determine ${\mathit{\alpha }}_{k^{\prime },k}^{*}$ for all such k^′ and k; i.e., all ${\mathit{\alpha }}_{k^{\prime },k}^{*}$ for all k^′ and k are determined once and for all by the minimization of ℒ¹.

Unfortunately, as was the case in Peherstorfer et al. (2016, 2018), the results given in Box 1 do not immediately lead to a practical MFMC method because the sample numbers $M_{k^{\prime },k}^{*}$ given in Eqs. (10) and (11) are not, in general, integers. Moreover, as was also the case in Peherstorfer et al. (2016, 2018), the first and perhaps even the first few sampling numbers at stage $k^{\prime }=\mathrm{1}$ may be <1 in addition to being non-integer. In those papers, a rounding procedure is implemented; i.e., the sample numbers $M_{k^{\prime },\mathrm{1}}$ are replaced by integers, though this is unsuitable for scenarios where $M_{\mathrm{1},\mathrm{1}}<\mathrm{1}$ as the choice $M_{\mathrm{1},\mathrm{1}}=\mathrm{0}$ leads to a biased estimator while the choice $M_{\mathrm{1},\mathrm{1}}=\mathrm{1}$ exceeds the computational budget. On the other hand, the modified MFMC method in Box 2 below is constructed to avoid these issues while preserving the optimality of the original solution (up to some amount of rounding). Note that in that box and elsewhere, ⌊ ⋅ ⌋ denotes rounding downwards to the nearest integer.

Box 2 Practical near-optimal integer sampling numbers for the modified MFMC method; see Gruber et al. (2022a).

3 Tests for the single-layer rotating shallow-water equations

Consider the single-layer rotating shallow-water equations (RSWEs) posed on the domain $\mathrm{\Gamma }\times [\mathrm{0},T]$ and given by (see Vallis, 2012)

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\displaystyle \frac{\partial h}{\partial t}}+\mathrm{\nabla }\cdot \left(h\mathit{u}\right)=\mathrm{0},\\ & {\displaystyle \frac{\partial \mathit{u}}{\partial t}}+(\mathit{k}\cdot \mathrm{\nabla }\times \mathit{u}+f)(\mathbf{k}\times \mathit{u})\\ & +\mathrm{\nabla }({\displaystyle \frac{|\mathit{u}{|}^{\mathrm{2}}}{\mathrm{2}}}+g(h+h_{\mathrm{b}}\left)\right)=G(h,\mathit{u}),\end{array}\end{array}$$

where Γ denotes the surface of a sphere or a subset of that surface, [0,T] denotes a time interval; k denotes a unit vector perpendicular to the surface of the sphere; h(x,t) denotes the fluid thickness; u(x,t) denotes the vector velocity field tangential to the surface of the sphere; G(h,u) denotes the a forcing term that depends on the specific setting; f denotes the Coriolis parameter; h_b(x,t) denotes the bottom topography; g denotes the (constant) acceleration due to gravity; and ∇ denotes the tangential gradient – i.e., $\mathrm{\nabla }f=Df-(Df\cdot \mathit{k})\mathit{k}$ where D is the derivative operator of ℝ³.

Note that the expression $(\mathbf{k}\cdot \mathrm{\nabla }\times \mathit{u}+f)/h$ that appears in the velocity equation is known as the potential vorticity; see Vallis (2012). Supplementing Eq. (12) are the initial conditions h=h₀ and u=u₀ at t=0; and, if Γ is a strict subset of the surface of the sphere, the boundary condition $\mathit{u}\cdot \mathit{n}=\mathrm{0}$ on the boundary ∂Γ of Γ, where n denotes the unit vector tangent to the surface of the sphere and also (outwardly) perpendicular to ∂Ω.

The RSWEs represent a useful simplification of the primitive equations (Vallis, 2012) which are commonly used in oceanic and atmospheric modeling and which are obtained by assuming a small ratio between the vertical and horizontal length scales. In this way, the single-layer RSWEs describe the motion of a thin layer of fluid which lies on a rigid surface, yielding conditions used extensively in the modeling of oceanic and atmospheric flows.

Spatial discretization of the system (Eq. 12) is effected using the TRiSK scheme (Ringler et al., 2010; Thuburn et al., 2009), which is a staggered C-grid mimetic finite difference–finite volume scheme that preserves desirable physical properties, including conservation laws for mass, energy, and potential vorticity. TRiSK discretization involves the approximation h_ℓ of the thickness h at the center of a grid cell P_ℓ, the approximation u_e of the normal to the edge component u⋅n of velocity at the centers of the edges of P_ℓ, and the approximation of the potential vorticity $(\mathbf{k}\cdot \mathrm{\nabla }\times \mathit{u}+f)/h$ at the vertices of P_ℓ. The necessary meshing is done using spherical centroidal Voronoi tessellation (SCVT) grids (Jacobsen et al., 2013; Ringler et al., 2008; Yang et al., 2018); examples of such grids are given in Sect. 3.1 and 3.2 for the specific settings of those sections.

Temporal discretization of the system (Eq. 12) is effected using an explicit fourth-order Runge–Kutta method, although many alternative time-stepping schemes have also been proposed for this purpose; see, e.g., Leng et al. (2019), Meng et al. (2020), and Trahan and Dawson (2012). We denote by $\mathit{\{}{t}_n{\mathit{\}}}_{n=\mathrm{0}}^{N_t}$ with t₀=0 and $t_{N_t}=T$ the time instants used for temporal discretization.

3.1 A wind-driven gyre test case for the RSWE system

The specific setting considered here is a modification of the benchmark test case referred to as “simulating ocean mesoscale activity” (SOMA; Wolfram et al., 2015), which involves a geodesic basin with radius 1250 km centered at latitude–longitude $(\mathrm{35}{}^{\circ },\mathrm{0}{}^{\circ })$ on the surface of the Earth. Inside the basin, the depth of fluid varies from 2500 m at the center to 100 m on the coastal shelf, creating a realistic topography which yields interesting dynamical behaviors which are useful for studying the propagation of fronts and eddies. The present experiment considers this same setting with a constant density ρ, leading to a computational model for a wind-driven barotropic gyre.

In Eq. (12),

$$\begin{array}{rl}G(h,\mathit{u})& =d_{\text{bottom}}(h,\mathit{u})+d_{\text{wind}}\left(h\right)\\ & \text{with}\left\{\begin{array}{ll}& d_{\text{bottom}}(h,\mathit{u})=-c_{\text{bottom}}\frac{\left|\mathit{u}\right|\mathit{u}}{\mathit{\rho }h}\\ & d_{\text{wind}}\left(h\right)=\frac{{\mathit{\tau }}_{\text{wind}}}{\mathit{\rho }h},\end{array}\right.\end{array}$$

with f_bottom(h,u) denoting a forcing term due to the bottom drag and f_wind(h) denoting a forcing term due to wind drag. Here, c_bottom denotes the bottom-drag coefficient chosen to be 10⁻³, ρ denotes the density, and τ_wind denotes the surface wind stress. See, e.g., Wolfram et al. (2015) for a discussion of this choice for G(h,u).

Table 1 Data corresponding to the three SCVT meshes used in the wind-driven gyre experiment.

For constructing the MC and MFMC estimators, we use three SCVT grids of the SOMA domain Γ listed in Table 1. A representative SCVT meshing of Γ is illustrated in Fig. 1. Note that for the SCVT grids used by the TRiSK scheme, the cells are almost all hexagonal, with a few pentagons and heptagons thrown in.

Figure 1 The 32 km grid used for the wind-driven gyre test case and a zoom-in on a portion of that mesh.

The output of interest we consider is the maximum fluid layer thickness, which is a common benchmark in oceanic RSWE simulations and which is relevant to, e.g., the detection of phenomena such as flooding. In particular, we simulate the RSWE until the final time of T=3 d using the finest grid (8 km resolution) of 120 953 cells, and then we choose the OoI $F_{\mathrm{gyre}}^{\mathrm{1}}$ given by

$$\begin{array}{}\text{(13)}& F_{\text{gyre}}^{\mathrm{1}}=\underset{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },\mathrm{120}\mathrm{953}}{max}{h}_{\mathrm{\ell }},\end{array}$$

where h_ℓ denotes the fluid thickness at the center of the cell P_ℓ at the final time T=3.

The quantity of interest we choose considers the effect that perturbations of the initial velocity u₀=u(0) have on $F_{\mathrm{gyre}}^{\mathrm{1}}$ and how that effect can be quantified using MC and MFMC estimation. To this end, consider random dilations (1+z)u₀ of the initial velocity depending on the i.i.d. random variable z that is uniformly distributed over the interval $[-\mathrm{0.5},\mathrm{0.5}]$. Then, the OoI defined in Eq. (13) depends on the choice of z in that interval; i.e., we have it that $F_{\mathrm{gyre}}^{\mathrm{1}}=F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$. In particular, we choose the QoI to be the expectation $Q_{\mathrm{1},\mathrm{gyre}}=\mathbb{E}\left(F_{\mathrm{gyre}}^{\mathrm{1}}\right)$ of the output of interest $F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$ defined in Eq. (13).

Note that in practical RSWE simulations, as is the case for more sophisticated models such as the primitive equations, an approximation of the initial data u₀ is often obtained from a pre-processing procedure in which the RSWE system is spun-up from rest, i.e., from a zero initial condition for u, up to some specified time; see Anderson et al. (1975) and Bleck and Boudra (1986). This procedure is invoked so as to eliminate transient artifacts which are not present in the current ocean or atmosphere and produces an initial configuration which is closer to observed oceanic and atmospheric data. The outcome of the spin-up calculation over the spin-up time frame of 15 d is illustrated in Fig. 2. The initial conditions u₀ and h₀ that supplement Eq. (12) are then simply set to the outcome of this pre-processing step, i.e., to the fluid thickness and velocity obtained at the end of the spin-up calculation.

Figure 2 For the setting of Sect. 3.1, the RSWE truth model solution (8 km resolution) with thickness h (a) and velocity field u (b) after integration of the system from rest for T=15 d.

3.1.1 MC and MFMC estimators

The MC estimator $Q_{\mathrm{1},\mathrm{gyre}}^{\mathrm{MC}}$ of $Q_{\mathrm{1},\mathrm{gyre}}=\mathbb{E}\left(F_{\mathrm{gyre}}^{\mathrm{1}}\right)$ is given by Eq. (2). Unfortunately, obtaining acceptable accuracy using an MC estimator suffers from a double shortcoming. First, for any given z, obtaining $F_{\mathrm{gyre}}^{\mathrm{1}}$ is a costly endeavor because it requires the solution of the discretized RSWE system to obtain the necessary 120 953 values of h_ℓ(z). Second, to obtain an MC estimator that is of acceptable accuracy requires obtaining $F_{\mathrm{gyre}}^{\mathrm{1}}$ for many sample values of z.

Naturally, to mitigate this double shortcoming, we turn to the MFMC estimator described in Box 2 of Sect. 2. To do so, we define three surrogate OoIs, $F_{\mathrm{gyre}}^{\mathrm{2}}\left(z\right)$, $F_{\mathrm{gyre}}^{\mathrm{3}}\left(z\right)$, and $F_{\mathrm{gyre}}^{\mathrm{4}}\left(z\right)$, all three of which are less costly to obtain compared to that of $F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$. Here, $F_{\mathrm{gyre}}^{\mathrm{2}}\left(z\right)$ and $F_{\mathrm{gyre}}^{\mathrm{3}}\left(z\right)$ are simply based on solving the discretized RSWE system using the coarser 16 and 32 km grids, respectively. The third surrogate, $F_{\mathrm{gyre}}^{\mathrm{4}}\left(z\right)$, is the piecewise-linear interpolant based on the values of $F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$, which are exact at the three points $z=\mathit{\{}-\mathrm{0.5},\mathrm{0},\mathrm{0.5}\mathit{\}}$. These four OoIs constitute a reasonable multifidelity ensemble of models that are available during computational ocean studies, even more realistic ones such as primitive equation models. Moreover, all four OoIs can be leveraged for MFMC estimation, whereas more traditional MC and multilevel MC estimators only make use of the first or the first three OoIs, respectively. Approximations of the costs and correlations for the four OoIs $\mathit{\left\{}{F}_{\mathrm{gyre}}^k\right(z){\mathit{\}}}_{k=\mathrm{1}}^{\mathrm{4}}$ are obtained by considering 100 uniform i.i.d. samples of $z\in [-\mathrm{0.5},\mathrm{0.5}]$. These are computed to be

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathbf{C}& =\left(\begin{array}{c}{C}_{\mathrm{1}}\\ C_{\mathrm{2}}\\ C_{\mathrm{3}}\\ C_{\mathrm{4}}\end{array}\right)=\left(\begin{array}{c}\mathrm{101.1}\\ \mathrm{12.83}\\ \mathrm{1.714}\\ \mathrm{0.05}\end{array}\right)\\ & \text{and}{\mathit{\zeta }}_{\mathrm{1}}=\left(\begin{array}{c}{\mathit{\zeta }}_{\mathrm{1},\mathrm{1}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{2}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{3}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{4}}\end{array}\right)=\left(\begin{array}{c}\mathrm{1.00000000}\\ \mathrm{0.99975045}\\ \mathrm{0.99975920}\\ \mathrm{0.99974835}\end{array}\right),\end{array}\end{array}$$

where C_k for $k=\mathrm{1},\mathrm{2},\mathrm{3}$ denotes the average computation time (in wall-clock seconds) necessary to advance the relevant discretized RSWE system by one time step, computed using 500 time steps for the simulation parameter z=0.001. Note that the cost C₄ is assigned arbitrarily because the cost of evaluating the interpolant is negligible. These cost-correlation pairs are ideal for MFMC estimation; i.e., the surrogates are very well correlated with $F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$ and are much less costly to obtain compared to $F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$. Note that the 100 samples of z used to determine the data in Eq. (14) can be reused when determining the MC and MFMC estimators.

From Eq. (14), we observe that the more expensive surrogate $F_{\mathrm{gyre}}^{\mathrm{2}}\left(z\right)$ is slightly less correlated to $F_{\mathrm{gyre}}^{\mathrm{1}}\left(z\right)$ than the cheaper surrogate $F_{\mathrm{gyre}}^{\mathrm{3}}\left(z\right)$ is. So, to satisfy the first criteria in Eq. (9), the correlations should be ordered in decreasing order: ζ_1,1, ζ_1,3, ζ_1,2, and ζ_1,4. However, it turns out that the second criteria in Eq. (9) for k=2 is then not satisfied; thus, $F_{\mathrm{gyre}}^{\mathrm{2}}\left(z\right)$ is removed from the list of surrogates when computing the MFMC estimator. That estimator is given by Eq. (4) with K=4 and where k=2 is excluded from the sum. This is an example which illustrates, as discussed in Sect. 2, that not all surrogates one chooses contribute to the efficiency of MFMC estimation and that, on the other hand, their omission does no harm to the accuracy of that estimator.

As already alluded to, the goal is to compare, for the same budget, the MSEs of the approximations $Q_{\mathrm{gyre}}^{\mathrm{MC}}$ and $Q_{\mathrm{gyre}}^{\mathrm{MFMC}}$ to the exact QoI $Q_{\mathrm{1},\mathrm{gyre}}=\mathbb{E}\left(F^{\mathrm{1}}\right)$ produced by the MC and MFMC estimators. However, there is still an additional challenge to be dealt with, namely that it is unclear how to choose a reference QoI that one can use to determine these MSEs, as only limited data (200 high-fidelity samples) are available, and an MC approximation $Q_{\mathrm{gyre}}^{\mathrm{MC}}$ with 200 samples is not sufficient for an accurate estimate of the expectation 𝔼(F¹). Therefore, we choose $Q_{\mathrm{gyre}}^{\mathrm{ref}}$ to denote the average of the MFMC approximation $Q_{\mathrm{gyre}}^{\mathrm{MFMC}}$ at the highest budget B=128C₁ taken over 250 runs using samples z that are independently drawn from the interval $[-\mathrm{0},\mathrm{5},\mathrm{0.5}]$, except for the fact that the first M₃ (recall that M₂ is omitted) samples of each run are drawn from the pre-collected sampling set of size 200. The use of MFMC in defining $Q_{\mathrm{gyre}}^{\mathrm{ref}}$ makes use of the fact that the MFMC method is at least as accurate as MC estimation, which can be verified through the inequality (Eq. 7). With this choice, the MSEs in the MFMC and MC estimators can be measured with respect to the “exact quantity” $Q_{\mathrm{gyre}}^{\mathrm{ref}}$.

The results of this procedure are given in Table 2 and Fig. 3, from which it is evident that MFMC produces a much more precise estimate compared to that of MC. In addition to MSE $e\left(Q_{\mathrm{gyre}}^{\mathrm{MFMC}}\right)$, we report the relative MSE defined as $e_{\mathrm{rel}}\left(Q_{\mathrm{gyre}}^{\mathrm{MFMC}}\right)=e\left(Q_{\mathrm{gyre}}^{\mathrm{MFMC}}\right)/(Q_{\mathrm{gyre}}^{\mathrm{ref}}{)}^{\mathrm{2}}$. It is interesting that most of the computational budget is loaded onto the very crude piecewise-linear interpolant approximation $F_{\mathrm{gyre}}^{\mathrm{4}}\left(z\right)$, allowing MFMC to achieve not only a lower MSE but also an estimate with much smaller variance from run to run. Especially telling is the bottom plot in Fig. 3, from which it is obvious that, for the same budget, MFMC estimation results in greater accuracy of 2 or 3 orders of magnitude compared to that MC estimation – or, looking at it another way, for the same relative MSE, MFMC estimation requires a much smaller budget compared to MC estimation.

Table 2 Results of the wind-driven barotropic gyre test with perturbed initial velocities for budgets B=2^kC₁ equivalent to 2^k highest-fidelity runs.

Figure 3 Results for the wind-driven barotropic gyre test with output of interest (Eq. 13) averaged over 250 applications of MC and MFMC estimation. (a) The quantity of interest Q as a function of the budget for the MC and MFMC estimators. (b) A zoom-in of the top figure. (c) The relative MSE of the MC and MFMC estimators as a function of the budget; shaded regions represent the standard deviations of the MC and MFMC predictions when compared to their averages over the 250 runs. Note that the melon region is contained in the blue region, indicating that the MFMC estimator has a smaller standard deviation than the MC estimator.

3.2 Test Case 5 for the RSWE system

We again consider the RSWE system (Eq. 12), but now Ω denotes the whole surface of the sphere. Specifically, we consider the configuration of Test Case 5 defined in Williamson et al. (1992), which is widely used as a benchmark and employed as a stepping stone towards more realistic atmospheric models; see, e.g., Leng et al. (2019), Meng et al. (2020), Ringler et al. (2010), and Williamson et al. (1992). Note that, for this example, there are no additional forces acting on the system so that $G(h,\mathit{u})=\mathrm{0}$ in Eq. (12).

Test Case 5 considers the flow over an isolated mountain centered at longitude ${\mathit{\lambda }}_{\mathrm{c}}=\frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}$ and latitude ${\mathit{\theta }}_{\mathrm{c}}=\frac{\mathit{\pi }}{\mathrm{6}}$ with height

$$\begin{array}{}\text{(15)}& h_{\mathrm{b}}\left(z_{\mathrm{1}}\right)=z_{\mathrm{1}}\left(\mathrm{1}-{\displaystyle \frac{r}{a}}\right),\end{array}$$

where $a=\frac{\mathit{\pi }}{\mathrm{9}}$, $r^{\mathrm{2}}=min\mathit{\{}{a}^{\mathrm{2}},(\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{c}}{)}^{\mathrm{2}}+(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{c}}{)}^{\mathrm{2}}\mathit{\}}$, and λ,θ denote longitude and latitude, respectively. In Eq. (15), z₁ denotes a random variable that is uniformly distributed over the interval [1 km, 3 km].

The initial tangential (to the sphere) velocity in the longitudinal and latitudinal directions is chosen to be ${\mathit{u}}_{\mathrm{0}}\left(z_{\mathrm{2}}\right)=(z_{\mathrm{2}}\cos \mathit{\theta },\mathrm{0})$, where z₂ denotes a random variable that is uniformly distributed over the interval [15 m s⁻¹, 25 m s⁻¹]. The initial fluid thickness h is chosen as

$$h_{\mathrm{0}}\left(z_{\mathrm{2}}\right)=\widehat{h}-{\displaystyle \frac{\mathrm{1}}{g}}\left(R\mathit{\omega }{z}_{\mathrm{2}}+{\displaystyle \frac{z_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}}}\right){\sin }^{\mathrm{2}}\mathit{\theta },$$

where $\widehat{h}=\mathrm{5.96}$ km, R=6371.22 km, and $\mathit{\omega }=\mathrm{7.292}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹. With this, the solutions u(z) and h(z) of the RSWE system (Eq. 12) depend on the random vector $\mathit{z}=(z_{\mathrm{1}},z_{\mathrm{2}})\in$ [1 km, 3 km] × [15 m s⁻¹, 25 m s⁻¹]. In Fig. 4, we provide an example of the initial thickness h₀(z₂) and the thickness h(z₁,z₂) after 10 d for specific values of z₁ and z₂.

Figure 4 The initial thickness h₀(z₂) with z₂=23.82865 m (a) and the thickness h after 10 d (b) and the mountain height h_b(z₁) with z₁=2023.78 m.

For the simulation results given in Fig. 4 and for other results in this subsection, we use the TRiSK scheme for spatial discretization and a fourth-order explicit Runge–Kutta method for temporal discretization. SCVT uniform gridding is employed as in Sect. 3.1, although now the grid covers the whole surface of the sphere. A comparison of two such SCVT grids with 480 and 240 km resolutions is provided in Fig. 5. In constructing the MC and MFMC estimators, we use the three globally refined SCVT meshes of the whole sphere according to the data in Table 3.

Table 3 Data corresponding to the three SCVT meshes used in the Test Case 5 experiment.

Figure 5 Two global SCVT meshes of the sphere surface with different grid resolutions: 480 km (a) and 240 km (b).

For the output of interest, we choose

$$F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\left(\mathit{z}\right)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{\ell }}^{\mathrm{1}}}}\underset{\mathrm{\ell }=\mathrm{1},\mathrm{\dots },N_{\mathrm{\ell }}^{\mathrm{1}}}{max}\left|u_{e,\mathrm{\ell }}\right(\mathit{z}\left)\right|,$$

where $N_{\mathrm{\ell }}^{\mathrm{1}}$ denotes the number of cell edges for the finest resolution case of 120 km, and u_e,ℓ(z) denotes the value of the normal component of velocity u_e at the ℓth cell edge for any choice of the random vector $\mathit{z}=(z_{\mathrm{1}},z_{\mathrm{2}})\in$ [1 km, 3 km] × [15 m s⁻¹, 25 m s⁻¹]. We then have it that the quantity of interest is given by

$$\begin{array}{}\text{(16)}& Q_{\text{1,test5}}=\mathbb{E}\left(F_{\text{test5}}^{\mathrm{1}}\left(\mathit{z}\right)\right).\end{array}$$

3.2.1 MC and MFMC estimators

Similarly to before, the goal is now to construct and compare, for the same computational budget, MC and MFMC estimators of the QoI defined in Eq. (16). The MC estimator $Q_{\text{1,test5}}^{\mathrm{MC}}$ of the QoI $Q_{\mathrm{1},\mathrm{test}\mathrm{5}}=\mathbb{E}\left(F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\right)$ is given by Eq. (2). The MFMC estimator $Q_{\mathrm{1},\mathrm{test}\mathrm{5}}^{\mathrm{MFMC}}$ of $Q_{\mathrm{1},\mathrm{test}\mathrm{5}}=\mathbb{E}\left(F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\right)$ makes use not only of the MC estimator $Q_{\mathrm{1},\mathrm{test}\mathrm{5}}^{\mathrm{MC}}$ for the finest 120 km grid but also of the MC estimators $Q_{\mathrm{2},\mathrm{test}\mathrm{5}}^{\mathrm{MC}}\approx \mathbb{E}\left(F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\right)$ and $Q_{\text{3,test5}}^{\mathrm{MC}}\approx \mathbb{E}\left(F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\right)$, corresponding the coarser 240 and 480 km grids, respectively.

In this case, 4500 uniform i.i.d. realizations of z are drawn for pre-computation, of which all models share 1200 and the two lower-fidelity surrogates share the entire 4500. Using 100 samples of the random variable z, the costs and correlation coefficients of these models are respectively approximated by

$$\begin{array}{rl}\mathbf{C}& =\left(\begin{array}{c}{C}_{\mathrm{1}}\\ C_{\mathrm{2}}\\ C_{\mathrm{3}}\end{array}\right)=\left(\begin{array}{c}\mathrm{434.8}\\ \mathrm{126.9}\\ \mathrm{58.04}\end{array}\right)\\ & \text{and}{\mathit{\zeta }}_{\mathrm{1}}=\left(\begin{array}{c}{\mathit{\zeta }}_{\mathrm{1},\mathrm{1}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{2}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{3}}\end{array}\right)=\left(\begin{array}{c}\mathrm{1.00000000}\\ \mathrm{0.99986604}\\ \mathrm{0.99925882}\end{array}\right),\end{array}$$

where the costs C₁, C₂, and C₃ denote the computation time (in wall-clock seconds) necessary to advance the relevant discretized RSWE system by one time step. As was the case for the wind-driven gyre experiment, these cost-correlation pairs are ideal for MFMC estimation; i.e., the surrogates OoIs $F_{\mathrm{test}\mathrm{5}}^{\mathrm{2}}\left(\mathit{z}\right)$ and $F_{\mathrm{test}\mathrm{5}}^{\mathrm{3}}\left(\mathit{z}\right)$ are very well correlated with $F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\left(\mathit{z}\right)$ and are much less costly to obtain compared to $F_{\mathrm{test}\mathrm{5}}^{\mathrm{1}}\left(\mathit{z}\right)$. Note that the 100 samples of z used to determine the data in Eq. (14) can be reused when determining the MC and MFMC estimators. From this point on, all other details about the construction and use of the MC and MFMC estimators are the same as for the wind-driven barotropic gyre test case of Sect. 3.1.

The results provided in Table 4 and Fig. 6 show that MFMC estimation produces a much more precise estimate compared to MC estimation. As was true for the wind-driven gyre test case, especially telling is the bottom plot in Fig. 6, from which it is obvious that, for the same budget, MFMC estimation results in an order of magnitude greater accuracy compared to that MC estimation – or, looking at it another way, for the same relative MSE, MFMC estimation requires a much smaller budget compared to MC estimation.

Table 4 Results of Test Case 5 with perturbed initial velocities for budgets B=2^kC₁ equivalent to 2^k high-fidelity runs.

Figure 6 Results over 250 runs of the RSWE system for the Test Case 5 experiment with quantity of interest given by Eq.(16). Shaded regions represent the variance in the MC and MFMC predictions over the 250 runs, which are independent except for the fact that all random sampling employs the same pre-collected (random) set of parameters. Note that the melon region is contained in the blue region, indicating that the MFMC estimator has a smaller standard deviation than the MC estimator.

4 First-order ice sheet model

The next experiment we consider illustrates the effectiveness of MFMC estimation on a QoI important for the realistic modeling of ice sheets such as those found near, e.g., Greenland, Antarctica, and various glaciers.

4.1 The first-order model for ice sheets

The dynamical behavior of ice sheets is commonly modeled by what is referred to as the first-order model or the Blatter–Pattyn model. Here, we provide a short review of that model; detailed descriptions are given in, e.g., Blatter (1995), Pattyn (2003), Perego et al. (2012), and Tezaur et al. (2015).

Let Ω denote the three-dimensional domain occupied by the ice sheet with the boundary $\mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{s}}\cup {\mathrm{\Gamma }}_{\mathrm{b}}\cup {\mathrm{\Gamma }}_{\mathrm{\ell }}$, given by

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathrm{\Gamma }}_{\mathrm{s}}& \Leftarrow x_{\mathrm{3}}=s(x_{\mathrm{1}},x_{\mathrm{2}})\Rightarrow \text{top surface of the ice sheet}\\ {\mathrm{\Gamma }}_{\mathrm{b}}& \Leftarrow x_{\mathrm{3}}=b(x_{\mathrm{1}},x_{\mathrm{2}})\Rightarrow \text{bottom (or basal) surface}\\ & \mathrm{\text{ of the ice sheet}}\\ {\mathrm{\Gamma }}_{\mathrm{\ell }}& \Leftarrow \mathrm{\ell }(x_{\mathrm{1}},x_{\mathrm{2}})=\mathrm{0}\Rightarrow \text{lateral boundary.}\end{array}\end{array}$$

Figure 7 provides an illustration of the boundary segments defined in Eq. (17).

Figure 7 An (x₁,x₃) cross section of the three-dimensional domain Ω occupied by an ice sheet and the boundary segments defined in Eq. (17).

Also, let u₁ and u₂ denote the x₁ and x₂ components of the velocity vector $\mathit{u}=(u_{\mathrm{1}}{u}_{\mathrm{2}}{u}_{\mathrm{3}}{)}^{⊺}$. Then, the first-order model equations for ice sheets is given by the following partial differential equations:

$$\begin{array}{}\text{(18)}& \left\{\begin{array}{ll}-\mathrm{\nabla }\cdot \left(\mathrm{2}\mathit{\mu }{\mathit{ϵ}}_{\mathrm{1}}\right)+\mathit{\rho }g\frac{\partial s}{\partial x_{\mathrm{1}}}& =\mathrm{0}\\ -\mathrm{\nabla }\cdot \left(\mathrm{2}\mathit{\mu }{\mathit{ϵ}}_{\mathrm{2}}\right)+\mathit{\rho }g\frac{\partial s}{\partial x_{\mathrm{2}}}& =\mathrm{0}\end{array}\right.\text{for}\mathit{x}=(x_{\mathrm{1}}{x}_{\mathrm{2}}{x}_{\mathrm{3}}{)}^{⊺}\in \mathrm{\Omega },\end{array}$$

where μ denotes the viscosity coefficient, g denotes the gravitational acceleration, and ρ denotes the density. In Eq. (18), the strain-rate tensor (ϵ₁,ϵ₂) is given by

$${\mathit{ϵ}}_{\mathrm{1}}=\left(\begin{array}{c}\mathrm{2}\frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{1}}}+\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}\\ \frac{\mathrm{1}}{\mathrm{2}}(\frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{2}}}+\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}})\\ \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{3}}}\end{array}\right)\text{and}{\mathit{ϵ}}_{\mathrm{2}}=\left(\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}}(\frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{2}}}+\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}})\\ \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{1}}}+\mathrm{2}\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}\\ \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}\end{array}\right)$$

and the nonlinear viscosity coefficient μ is given by the Glen flow law

$$\mathit{\mu }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\frac{\mathrm{1}}{n}}{\mathit{ϵ}}_{\mathrm{e}}^{\frac{\mathrm{1}}{n}-\mathrm{1}},$$

with n=3 being the usual choice. The effective strain rate ϵ_e is given by

$$\begin{array}{rl}{\mathit{ϵ}}_{\mathrm{e}}^{\mathrm{2}}& =({\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{1}}}}{)}^{\mathrm{2}}+({\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}}{)}^{\mathrm{2}}+{\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{1}}}}{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}({\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{2}}}}+{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{1}}}}{)}^{\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}({\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{3}}}}{)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}({\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{3}}}}{)}^{\mathrm{2}},\end{array}$$

and A is often chosen to obey the Arrhenius relation

$$A=A\left(T\right)=a\exp (-Q/RT),$$

where T denotes the absolute temperature measured in degrees kelvin, R denotes the universal gas constant, Q denotes the activation energy for creep, and a is an empirical flow constant often used as a tuning parameter.

The system (Eq. 18) is supplemented by the boundary conditions

$$\begin{array}{}\text{(19)}& \begin{array}{llll}{\mathit{ϵ}}_{\mathrm{1}}\cdot \mathit{n}=\mathrm{0}& \text{and}& {\mathit{ϵ}}_{\mathrm{2}}\cdot \mathit{n}=\mathrm{0}& \text{on }{\mathrm{\Gamma }}_{\mathrm{s}},\\ \mathrm{2}\mathit{\mu }{\mathit{ϵ}}_{\mathrm{1}}\cdot \mathit{n}+\mathit{\beta }u=\mathrm{0}& \text{and }& \mathrm{2}\mathit{\mu }{\mathit{ϵ}}_{\mathrm{2}}\cdot \mathit{n}+\mathit{\beta }v=\mathrm{0}& \text{on }{\mathrm{\Gamma }}_{\mathrm{b}},\\ u=\mathrm{0}& \text{and}& v=\mathrm{0}& \text{on}{\mathrm{\Gamma }}_{\mathrm{\ell }},\end{array}\end{array}$$

where β denotes a basal-friction parameter.

Once the horizontal components u₁ and u₂ of the velocity are determined, the vertical-velocity component w is determined by enforcing incompressibility; i.e., we have it that

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\partial w}{\partial x_{\mathrm{3}}}}=-{\displaystyle \frac{\partial u_{\mathrm{1}}}{\partial x_{\mathrm{1}}}}-{\displaystyle \frac{\partial u_{\mathrm{2}}}{\partial x_{\mathrm{2}}}}\text{for }\mathit{x}\in \mathrm{\Omega }.\end{array}$$

Because the right-hand side is known, this is an ordinary differential equation for w.

Of paramount interest in the modeling of ice sheets is the monitoring of the temporal evolution of the ice sheet domain Ω. However, one notices that there are no time derivatives of u₁ and u₂ appearing in Eq. (18); i.e., that system is a static one. The reasoning behind this curiosity is twofold. Firstly, we have it that the system (Eq. 18) is coupled to an equation for the temporal evolution of the temperature T within the ice sheet and is also coupled to an equation for the temporal evolution of the top surface of the ice sheet, which engenders changes in the ice sheet domain Ω. Secondly, the timescale of changes in the temperature is much shorter (e.g., hourly) when compared to the timescale (e.g., at least many days) of changes in the ice sheet domain Ω.

Thus, in a computational model of ice sheet dynamics, determination of the domain Ω is coupled to that of the velocity u. Precisely, the temperature and top surface are first advanced from a given domain–velocity pair (Ω,u) over several time steps, producing an evolution that prescribes a new domain–temperature pair (Ω,T). Following this, the pair (Ω,T) are used to generate a new velocity u by solving Eqs. (18) and (20). While both stages of this process are important, the present experiment focuses on the computation of u from Ω and T.

4.2 MC and MFMC estimation

The specific application setting we consider is based on Experiment C of the benchmark examples in Leng et al. (2012) and Pattyn et al. (2008). Here, the ice domain is a rectangular parallelepiped that has a square base with side length L and thickness H and which lies on a slanted bed with slope parameter θ. The upper surface boundary is given as

$$z_{\mathrm{s}}(x_{\mathrm{1}},x_{\mathrm{2}})=-x_{\mathrm{1}}\tan \mathit{\theta },$$

and the basal topography is given as

$$z_{\mathrm{b}}(x_{\mathrm{1}},x_{\mathrm{2}})=z_{\mathrm{s}}(x_{\mathrm{1}},x_{\mathrm{2}})-H.$$

For constructing the MC and MFMC estimators, we set the length L=80 km and height H=1 km and then use uniform tetrahedral grids of the ice domain Ω_t with 120, 60, and 30 grid intervals in each horizontal direction and, correspondingly, 20, 10, and 5 intervals in the vertical direction. As a result, we have that number of vertices, horizontal triangles, and tetrahedra are determined according to the data in Table 5.

Table 5 Data corresponding to the three tetrahedral meshes used in the ice sheet experiment.

The first-order ice sheet model is discretized using the stabilized P1–P1 finite elements given in Zhang et al. (2011). To solve the resulting nonlinear system of discrete equations, 15 Picard iteration steps are carried out, after which a switch is made to a Newton iteration, up to a maximum of a total of 40 iterations. However, if the residual error decrease is less than 75 % relative to the residual error of the previous step, the nonlinear iteration is switched back to a Picard iteration. An example illustration of the discrete solution of this ice sheet model for specific values of the parameters θ and β is provided in Fig. 8 using the highest-resolution grid.

Figure 8 A simulation of the ice sheet model with θ=0.7387 and β=901.5232, discretized at the highest grid resolution. The top two plots display contours of the velocity components u₁ and u₂, respectively, on the top surface Γ_s. The third plot displays a contour plot of velocity component u₃. The fourth plot displays a cutaway plot of the velocity magnitude $\left|\mathit{u}\right|$ inside the sheet.

Because the cracking and melting of ice sheets is an important indicator of climate change (Clark et al., 1999; Hanna et al., 2013), we consider the the output of interest to be

$$F_{\mathrm{ice}}^{\mathrm{1}}\left(\mathit{z}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{N}_{\mathrm{1}}}}\sum _{\mathrm{\ell }=\mathrm{1}}^{N_{\mathrm{1}}}\left|{\mathit{u}}_n\right(\mathit{z}){|}^{\mathrm{2}},$$

where N₁=307 461 is the number of vertices of the finest $\mathrm{120}\times \mathrm{120}\times \mathrm{20}$ grid, and u_n(z) denotes the value of the discrete velocity u at the nth vertex for any choice of the i.i.d. random vector $\mathit{z}=(\mathit{\theta },\mathit{\beta })\in [\mathrm{0.2},\mathrm{0.8}]\times [\mathrm{800}\mathrm{m},\mathrm{1200}\mathrm{m}]$. This OoI provides a measurement of how energetic the ice sheet is depending on its slope and friction coefficient and can be loosely related to how vigorously the ice will deform given a configuration specified by z. We then choose the quantity of interest, given by

$$\begin{array}{}\text{(21)}& Q_{\mathrm{1},\mathrm{ice}}=\mathbb{E}\left(F_{\mathrm{ice}}^{\mathrm{1}}\left(\mathit{z}\right)\right).\end{array}$$

Two surrogates for the ice sheet model,

$$\begin{array}{rl}{F}_{\mathrm{ice}}^{\mathrm{2}}\left(\mathit{z}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{N}_{\mathrm{2}}}}\sum _{n=\mathrm{1}}^{N_{\mathrm{2}}}\left|{\mathit{u}}_n\right(\mathit{z}){|}^{\mathrm{2}}\text{and}\\ F_{\mathrm{ice}}^{\mathrm{3}}\left(\mathit{z}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{N}_{\mathrm{3}}}}\sum _{n=\mathrm{1}}^{N_{\mathrm{3}}}\left|{\mathit{u}}_n\right(\mathit{z}){|}^{\mathrm{2}},\end{array}$$

are defined respectively using the coarser $\mathrm{60}\times \mathrm{60}\times \mathrm{10}$ grid with N₂=40 931 and the even coarser $\mathrm{30}\times \mathrm{30}\times \mathrm{5}$ grid with N₃=5766. Then, from Eq. (3), we have the corresponding three Monte Carlo estimators $\mathit{\{}{Q}_{\mathrm{1},\mathrm{ice}}^{\mathrm{MC}},Q_{\mathrm{2},\mathrm{ice}}^{\mathrm{MC}},Q_{\mathrm{3},\mathrm{ice}}^{\mathrm{MC}}\mathit{\}}$, and from Eq. (4), we have the MFMC estimator $Q_{\mathrm{ice}}^{\mathrm{MFMC}}$, which makes use of these three MC estimators.

At this point, we proceed as was done in Sect. 3.2. For example, we now have that the approximate costs (measured in wall-clock seconds and averaged over 100 random samples) for computing $F_{\mathrm{ice}}^{\mathrm{1}}\left(\mathit{z}\right)$, $F_{\mathrm{ice}}^{\mathrm{2}}\left(\mathit{z}\right)$, and $F_{\mathrm{ice}}^{\mathrm{3}}\left(\mathit{z}\right)$, as well as the approximate correlations for the three OoIs, are given by

$$\begin{array}{rl}\mathbf{C}& =\left(\begin{array}{c}{C}_{\mathrm{1}}\\ C_{\mathrm{2}}\\ C_{\mathrm{3}}\end{array}\right)=\left(\begin{array}{c}\mathrm{285.5}\\ \mathrm{23.05}\\ \mathrm{2.690}\end{array}\right)\\ & \text{and}{\mathit{\zeta }}_{\mathrm{1}}=\left(\begin{array}{c}{\mathit{\zeta }}_{\mathrm{1},\mathrm{1}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{2}}\\ {\mathit{\zeta }}_{\mathrm{1},\mathrm{3}}\end{array}\right)=\left(\begin{array}{c}\mathrm{1.00000000}\\ \mathrm{0.99999796}\\ \mathrm{0.99996691}\end{array}\right).\end{array}$$

The remaining experimental details are nearly identical to those of Sect. 3.2. As was the case for the oceanic simulations, the computational expense of these models precludes the collection of unlimited data, so a total of 8000 uniform i.i.d. realizations of z are drawn, of which all models share 200 and the lower-fidelity surrogates share 1000. The results of carrying out MFMC and MC estimation over 250 “independent” runs using this common data set are provided in Table 6 and Fig. 9. Here, it is especially apparent that the budget-preserving modifications to MFMC that led to the algorithm in Box 2 were necessary, as $B\ll \mathbf{C}\cdot \mathbf{R}$ and because only one high-fidelity run is ever selected by the algorithm.

Despite this, it is clear from Table 6 and the plots in Fig. 9 that, within these small budgets, MFMC estimation produces an estimator which is much more accurate and precise than that of MC estimation. This is not surprising because MFMC estimation has the flexibility to load most of its budget onto the cheaper lower-fidelity surrogate models, evaluating them many times in order to bring down the overall variance of their estimates. Again, this provides empirical validation for the use of MFMC estimation over MC estimation when estimating model statistics, even in practical cases where data availability is low.

Table 6 Results of the ice sheet experiment for budgets equivalent to 2^k high-fidelity runs.

Figure 9 Results over 250 runs of the ice sheet experiment with quantity of interest (Eq. 21). Shaded regions represent variance in the MC and MFMC predictions over the 250 runs. Note that the melon region is contained in the blue region, indicating that the MFMC estimator has a smaller standard deviation than the MC estimator.

5 Concluding remarks

This paper serves to introduce multifidelity Monte Carlo estimation as an alternative to standard Monte Carlo estimation for quantifying uncertainties in the outputs of climate system models, albeit in very simplified settings. Specifically, we consider benchmark problems for the single-layer shallow-water equations relevant to ocean and atmosphere dynamics, and we also consider a benchmark problem for the first-order model of ice sheet dynamics. The computational results presented here are promising in that they amply demonstrate the superiority of MFMC estimation when compared to MC estimation on these examples. Furthermore, the use of MFMC as an estimation method will surely be even more efficacious when quantifying uncertainties in more realistic climate-modeling settings for which the simulation costs are prohibitively large, e.g., for long-term climate simulations. Thus, our next goal is to apply MFMC estimation to more useful models of climate dynamics (such as the primitive equations for ocean and atmosphere and the Stokes model for ice sheets) that are also coupled to the dynamics of other climate system components and also to passive and active tracer equations.