2.1 Installing the software

The Matilda software is available on GitHub. To install from our GitHub repository:

1)

2)

2.2 Software documentation

Full descriptions of package functions can be accessed in the package’s help documentation. Furthermore, detailed documentation and vignettes are available from our GitHub repository (github.com/jgcri/matilda).

2.3 Configuring a model core

An analysis in Matilda begins by setting up a Hector model instance, termed a “core”. A Hector core is an object that contains information about model inputs, current state, and outputs for a specific Hector run. The information contained in an initiated core comes from an INI file holding metadata, emissions scenarios, and model parameters needed to run Hector.

We call newcore() (a hector function) to initiate a core containing information to conduct model runs using the SSP 2–4.5 emission scenario:

3)

2.4 Parameter estimation and establishing parameter sets

The basis of running Hector in a probabilistic setup relies on establishing a set of parameter configurations that are used to run the model iteratively. Matilda uses parameter information gathered from the literature to inform prior distributions (Table 1). To build parameter sets, we draw parameter values from their prior distributions using Monte Carlo sampling. Each parameter is sampled from its marginal prior distribution independently of the rest of the parameters. The prior distribution for each parameter θ_i is defined using mean and standard deviation estimates as in: (1) where θ_i is a given Hector parameter and N(μ, σ) is the normal distribution of parameter θ_i using hyperparameters μ (mean) and σ (standard deviation). Some parameters have marginal distributions best represented using lognormal distribution, in such cases, (μ, σ) is substituted by log(μ, σ) (Table 1). Using informed prior marginal distributions from the literature as a starting point for building perturbed parameter sets enables the exploration of a range of possible parameter values from the multi-dimensional parameter space judged viable on the basis of existing knowledge [13,25]. Once each parameter set of draws is performed, the full parameter vector values are combined using a uniform multivariate distribution. This process assumes independence of parameters, meaning for example that any value from the univariate draws of parameter θ₁ is equally likely to be paired with any value from the univariate draws of parameter θ_2…n. This parameter estimation process ultimately establishes parameter sets that account for parameter uncertainty and can be used to build an ensemble of model runs that will be evaluated against observational evidence for some salient output of those runs. In other words, we use prior information about individual parameters to build parameter sets, remaining agnostic about which sets will result in the most skilled model results until confronting Hector’s output with observed data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Hector parameters used in Matilda.

Hector parameters used to generate parameter sets in this work. The distributions are indicated as mean ± standard deviation. References from where distributions are derived are included.

https://doi.org/10.1371/journal.pclm.0000295.t001

In Matilda, we build parameter sets by calling generate_params(). Parameter distributions are independent of the SSP scenario, however, to run this function the user must still provide an established Hector core. Additionally, the user must specify the number of parameter sets desired (draws). Using generate_params() will produce randomized draws each time it is run. Therefore, the user should either save the resulting data frame or use set.seed() if replication of parameter sets is critical to the analysis. In this example we use our previously established core to produce a set of 25 parameter configurations and display a subset of samples from the result:

4)

Parameters can be easily added or omitted from the new parameter set data frame. For example, to run the model with a subset of parameters, undesired columns can be omitted from the data frame. This will result in a data frame that only contains parameter distributions that the user wishes to perturb. Similarly, the user can characterize additional parameter distributions and add them as new columns to the data frame, as long as the parameter is described in Hector.

Once established, the parameter sets are used as inputs for independent Hector model runs. Thus, each model run represents a multivariate parameter combination as follows: (2) where m_i is an individual ensemble member and , are parameters sampled to build an independent configuration. Using different parameter sets to run Hector allows us to build PPEs and determine how different parameter combinations from a presumably viable parameter space interact to affect climate variable projections. This method effectively propagates parameter uncertainty to model ensemble uncertainty, a process described as forward uncertainty propagation [32].

2.5 Forward uncertainty propagation and running the model

We run Hector for each of the parameter sets by calling iterate_model(), which runs the model for each parameter set and combines the results into a data frame object representing the new PPE. To run iterate_model(), the same core object is used as in previous steps and we also must supply the object where parameter sets are stored:

5)

The resulting data frame returns 25 separate runs, as indicated by the run_number column; in this case, the total number of rows is 55600 (25 runs × 4 output variables × 556 years). Each run includes values for the major climate variables of a Hector default output (CO₂ concentration, total radiative forcing, CO₂ forcing, and global mean air temperature) for the years 1745–2300 (the time range defined by the SSP INI file we chose above).

While a core object and a data frame of parameter sets are the only required arguments to run iterate_model(), additional arguments can be supplied to reduce the variables and year range returned for each run using save_vars and save_years, respectively. This reduces the size of the data stored in memory, which may be important when running the model to build large ensembles (e.g., 15,000 members as in [13]). Any output variable from Hector can be returned using save_vars() for any year range subset from 1745–2300. In the following example, we supply these arguments to return values only for CO₂ concentration and global mean surface temperature anomaly for the year range 1745–2100:

6)

The resulting data frame has only 17800 rows, a 68% savings over the full example above.

2.6 Model evaluation approach and scoring model runs

Evaluating ensemble members is important in climate model assessment because we do not know a priori which parameter sets yield realistic simulations. Evaluation procedures allow us to gain insight into the uncertainty of model outputs, thereby enhancing the fidelity of results. The concept of weighting ensemble members is intuitive; members that are more skilled (i.e., agree better with the historical record) should receive a higher weight than members that are less skilled (i.e., present larger deviations from the historical record). Ensemble members closely aligned with historical climate data will contribute more information to our probabilistic projections than members with outputs deviating significantly from the historical record. Thus, the uncertainty from varying the model parameters a priori is “filtered” through an evaluation of the parameters performance against observations. Model weights are analogous to the likelihood factor in Bayes Theorem (defined in Sec. 2.6.1) and could potentially be used to update prior parameter distributions to approximate the posterior parameter distributions. While this is a capability of Matilda, we do not discuss it at length in this paper. It is also important to note the limitations to this approach. While the weighting approaches presented here can be effective for updating marginal distributions of parameters independently, they may not fully capture the complex dependencies present in multidimensional distributions. This limitation becomes noticeable as the dimensionality of parameter space increases. In these cases, formal Bayesian approaches (e.g., Markov Chain Monte Carlo) are better suited for exploring a multidimensional parameter space.

Scoring PPE members in Matilda is conducted using score_runs() which requires (1) a results data frame, (2) a scoring criterion, and (3) a scoring function/algorithm. The results data frame typically comes from calling iterate_model(), as above.

Scoring criteria define information used to compare ensemble members against observational data. A scoring criterion can be built by the user by calling new_criterion() and simply requires the climate variable to be used in the comparison, the years of comparison, and observed data values for the years specified. For example, a new criterion can be created based on global mean surface temperature from a dataset containing observed warming values from 1990–2023:

7)

This defines a custom criterion: a time series of 34 (1990–2023) values that will be compared against Hector’s “gmst” (global mean surface temperature anomaly) output variable.

The Matilda package has internally available scoring criteria for easy use, including criterion_co2_obs() and criterion_gmst_obs(). Data contained in criterion_co2_obs() is pulled from the Mauna Loa record of observed annual mean atmospheric CO₂ concentration [33], while criterion_gmst_obs() uses observed annual mean global surface temperature anomaly data retrieved from the HadCRUT5 data set [34].

Scoring functions in Matilda apply different mathematical algorithms to compute model weights based on the results and scoring criterion. We provide multiple mechanisms to weight model outputs against observations, and users can define their own custom functions as well. There are currently two internally available scoring functions called score_bayesian() and score_ramp(), that differ in functionality and computational complexity.

2.6.1 Scoring function: Score_bayesian().

Bayesian probability theory provides a rigorous framework combining prior information, observational data, and model simulations for parameter estimation, model evaluation, and uncertainty quantification [35,36]. The Bayesian weighting scheme presented here was influenced by Bayesian model averaging (BMA) methodology. In BMA, multiple model structures, each with specific parameterizations, are used to make projections [35,37–40]. In the BMA method each model is a priori assumed equiprobable but once observations are introduced, model probabilities change. Models that are more consistent with observations are considered more likely, while the likelihood of models that are inconsistent with observations is reduced [41]. During this process no model is completely excluded from contributing information, but models are down weighted if they do not accurately represent observations. In this way, score_bayesian() in Matilda behaves like BMA by computing weights for PPE members based on consistency with historical observations. Unlike BMA however, Matilda focuses solely on perturbing parameters of a single model (Hector) to quantify uncertainty and thus does not capture uncertainty associated with model structure. In this section, we detail the Bayesian weighting scheme performed in Matilda.

We approximate the formalism of Bayesian inference to analyze PPE members and assign weights according to the probability of ensemble member m_i conditional upon observed data Y. This is described as the posterior probability and, consistent with Bayes’ theorem, is proportional to the product of a chosen likelihood function of m_i given observed data and prior information of m_i (prior probability) [37]. We can express this in equation form as: (3) where P(m_i|Y) is the posterior probability of m_i conditional upon observed data Y, L(m_i|Y) is the chosen likelihood function, and P(m_i) is the prior information of ensemble member m_i.

As demonstrated in Eq 3, posterior probabilities are dependent on prior knowledge about ensemble member m_i and a likelihood function that quantifies the agreement between ensemble member m_i and observed data Y. Here, we use a normal likelihood distribution function based on root-mean-square error (RMSE) which is commonly used as a statistical evaluation of model performance in climate research and is optimal under the assumption that errors are normally distributed [42,43]. For a given time series of observed data _(t) (i.e., scoring criteria) and a corresponding time series of each ensemble member m_(t), RMSE is a quantification of the averaged difference between the observed and modeled data and is calculated using the following formula: (4)

here, RMSE_i is an independent RMSE value representing how well PPE member m_i agrees with observations where σ_Y(t) represents the standard deviation of the error associated with the observations Y at time t and N represents the total number of data points in the time series. We further use these RMSE_i values in our likelihood function. Assuming a normal distribution, our proportional likelihood can be calculated as: (5)

In this equation, a decay relationship exists between RMSE_i and (m_i|Y), indicating a gradual decrease in (m_i|Y) as RMSE_i increases. In other words, the likelihood of an ensemble member decreases as the disagreement with observations increases. The value of σ_RMSE in Eq 5 plays a role in establishing the sensitivity of the likelihood to increased RMSE values, controlling the rate of likelihood decay. The relationship is explained in detail below and can be visualized in Fig 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Decay rates from varying sensitivity values.

Root-mean-square error (RMSE) plotted against likelihood, conditional upon observed data. Different colors indicate decay rates for different sensitivity values in score_bayesian(). Setting higher values to sensitivity decreases the deviation penalty applied to ensemble members.

https://doi.org/10.1371/journal.pclm.0000295.g002

As described in Sec. 2.4, we enforce an equally distributed prior (m_i) across all ensemble members because although we use prior information about individual parameters of Hector to build parameter sets, we do not know which sets will result in the most skilled model results before considering observed data. Taking all this information together, weights are estimated in score_runs() as normalized probabilities of each ensemble member, taking into account both agreement with observed data and prior beliefs about ensemble members using the following formula: (6) where ω_i is a weight assigned to each ensemble member.

We provide an example of using score_bayesian() as the scoring function when calling score_runs(). For this example, we use the result produced in code block 5, and analyze the agreement between ensemble members and observed data with the criterion_co2_obs() scoring criterion:

8)

scored_hector_runs <- score_runs(results, criterion_co2_obs(), score_bayesian, sigma = 29.4, sensitivity = 2)

print(scored_hector_runs)

##    weights    run_number

##  1  0.093578592  1

##  2  0.005328194  2

##  3  0.138975159  3

##  4  0.008768552  4

##  5  0.040364027  5

The resulting data frame returns 25 weights assigned for each Hector run (indicated by run_number). Weights for each ensemble member will be nonzero positive values that sum to 1. Weights closer to 1 represent ensemble members in strong agreement with observed data and weights closer to 0 correspond to members that are not in strong agreement with observed data. It is important to note that weights are asymptotic and thus no one ensemble member can score exactly 1 or reach a value of exactly 0 due to the normalization process.

The sigma parameter in code block 8 corresponds to _(t) in Eq 7 and provides an option for the user to provide a time-varying penalty term in the RMSE_i calculation. This way, we address not only observational uncertainties but also changes in those uncertainties over time [44]. For example, uncertainties for global mean surface temperature data from ~1800s are several times higher in magnitude compared to temperature observations in more recent years. By providing an adjustable sigma parameter we account for such cases. By default, the sigma parameter assumes that the variability of residuals remains constant across all years (homoscedasticity) and uses the standard deviation of the observed data to penalize errors between Hector’s output and observations uniformly across all years. The user can indicate an alternative (constant) sigma value (as shown in code block 8). Alternatively, a vector of sigma values can be provided which will apply a time-varying penalty to the error term in the RMSE_i calculation.

The sensitivity parameter in code block 8 is a multiplier applied to σ_RMSE that sets the decay rate determining how quickly likelihood values decrease as RMSE values increase. Because sensitivity represents the unit of acceptable deviation, a lower sensitivity value leads to a faster decay rate, meaning that the likelihood decreases more rapidly with increasing RMSE values. Conversely, a higher sensitivity reduces the decay rate, and ensemble member likelihood decreases more slowly with respect to increasing RMSE values. This results in more weight being assigned to ensemble members that have a lower likelihood. In Fig 2, we show how setting different sensitivity values in score_bayesian() leads to different decay rates, meaning that weights are distributed differently depending on the magnitude of the sensitivity parameter. In effect, this determines how severely an ensemble member is penalized as it departs from a criterion.

This parameter thus gives users the ability to govern the sensitivity of the likelihood decay as RMSE values increase. By adjusting the sensitivity value, the user has control over the weight assigned to different ensemble members based on their RMSE values. Setting a lower sensitivity value will result in more weight placed on ensemble members with low RMSE values. In comparison, a higher sensitivity value will give relatively more weight to ensemble members with higher RMSE values. The default sensitivity value is set to `sensitivity = 1`which quantifies one-unit standard deviation of the RMSE values, this results in a relatively gradual tapering of the likelihood as RMSE increases. We also provide the option in score_bayesian() to set a value to sensitivity directly which will change the acceptable deviation from observed data.

We note that users should adjust sensitivity in line with the context and evaluation purpose specific to their analytical goals. The case study at the end of this paper provides an example of assessing ensemble member weights for different acceptable deviation limits and in Sec. 2.8 we offer a more general discussion of the caveats to model weighting.

2.6.2 Scoring function: Score_ramp().

The score_ramp() function is a simpler and more transparent scoring algorithm that computes the absolute difference between ensemble members m_i(t) and observed data Y(t) at each time step: (7)

Scores are then computed based on how far absolute differences (t) are from arbitrarily selected minimum (w₁) and maximum (w₂) divergence values. For example, D(t) ≤ w₁ indicates small differences between modeled and observed data at time t and will result in a score of 1, whereas D(t) ≥ w₂ indicates significant divergence of modeled data from observed data and will result in a score of 0 (Fig 3). In cases where D(t) falls between w₁ and w₂, scores are computed using a linear function that decreases as D(t) values get closer to w₂ and further from w₁ (Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Decay rate for score_ramp().

Example of decay method for score_ramp() where w₁ = 5 and w₂ = 10. Ensemble members with an average deviation from observation < 5 will score 1 and ensemble members with an average deviation > 10 will score 0. Scores of ensemble members with average deviation between w₁ and w₂ will decrease from 1 linearly as average deviation approaches w₂.

https://doi.org/10.1371/journal.pclm.0000295.g003

We can express this linear decay with the following formula: (8)

where S(t) is the score at time t. Once computed, scores are averaged across the entire time series, resulting in a single score for each ensemble member. Scores for ensemble members are normalized in score_runs() to assign a weight between 1–0 for skilled versus unskilled members, where more skilled ensemble members will be weighted closer to 1 while less skilled ensemble members will receive weights closer to 0. Similar to the normalization step above, weights are estimated using the following formula: (9) where ω_i is a weight assigned to each ensemble member from the normalized mean score of each member S_i(t).

Below we provide a code example using score_ramp() as the scoring function in a score_runs() call. As in code block 8, we use the results produced in code block 5 and assess the agreement between ensemble members and observed data with the criterion_co2_obs() scoring criterion:

9)

scores <- score_runs(result, criterion_co2_obs(), score_ramp, w1 = 0, w2 = 10)

print(scores)

##    weights    run_number

##  1  1.044760e-02  1

##  2  7.575543e-10  2

##  3  6.271913e-02  3

##  4  8.717546e-21  4

##  5  6.774904e-13  5

Similar to scoring ensemble members with score_bayesian(), the resulting data frame returns 25 weights assigned to each ensemble member (indicated by run_number).

Ensemble members weighted according to our scoring algorithms can be used to visualize the uncertainty ranges that different parameter combinations and/or parameter sampling distributions generate after observational evidence is brought to bear on their results. In Fig 4, we show all ensemble members weighted using our two scoring algorithms. The ensemble shading visually demonstrates how score_bayesian(), with a default sigma and sensitivity, distributes weights more evenly across likely ensemble members, whereas score_ramp() assigns higher weights to ensemble members falling closer to the minimum divergence range of w₁ (when w₁ = 0 and w₂ = 10 ppm).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Weighted ensemble members using different scoring algorithms.

Perturbed parameter ensemble (PPE) projections using 25 parameter sets plotted for atmospheric CO₂ concentration from 1960–2100 weighted using the A) score_ramp() and B) score_bayesian() algorithms. Ensemble member weights are indicated by color shading with the solid red line representing observed atmospheric CO₂ concentrations from 1959–2021.

https://doi.org/10.1371/journal.pclm.0000295.g004

A time test on an analysis weighting 1000 ensemble members shows that the difference in computation time between the two scoring algorithms is negligible, with both functions computing weights for N = 1000 runs in a fraction of a second.

2.7 Defining and calculating metrics and probabilities

Once the ensemble members are scored, we can use them to compute informative metrics from model projections. Calculating metrics from the final weighted PPE requires (1) a results data frame and (2) a metric object, which must first be defined by the user.

Metric objects determine what data the user is most interested in extracting and summarizing from the results data frame. For example, a metric object can identify information needed to estimate global mean air temperature anomaly (global_tas) for the 20-year average used by the IPCC to represent long-term temperature change (2081–2100). We complete this by calling the function new_metric():

10)

This defines a new custom metric object: obtain the mean global air temperature (global_tas) for the years 2081–2100. The argument op in code block 10 describes a statistical operation that can be performed on the model data to compute a descriptive statistic for each member of the ensemble. While we define a 20-year mean global_tas metric in this example, the user can also easily compute the median, min, max, standard deviation, etc. for each ensemble member. Additionally, users can specify a single year rather than quantifying statistics over a range of years (e.g., 2100 vs. 2081–2100).

Once this metric is defined, we call metric_calc() to compute metric values for each ensemble member using the results data frame:

11)

The resulting data frame returns 25 separate metric values (indicated by the metric_result column) representing the 2081–2100 mean warming of global air temperature.

When metrics of interest are calculated and weights are assigned to PPE members based on agreement with historical record, we have the necessary information to address questions such as, “What is the probability that long-term mean warming will be between 2.0–4.0°C relative to pre-industrial reference?”

We approach such a question by calling prob_calc(), a function that sums PPE weights as they are binned into metric ranges identified by the user. Running prob_calc() requires (1) a data frame where metric values can be identified, (2) bins defined by the user, and (3) a data frame where PPE weights can be identified. Here, we provide an example of prob_calc() usage:

12)

In this example we use PPE weights computed using the score_bayesian() algorithm. The prob_calc() result shows the total probability that long-term projections of mean warming will fall within each of the temperature ranges defined by our bins (1.5–2°C, 2–2.5°C, 2.5–3°C, 3–3.5°C, 3.5–4°C, and >4°C) for the SSP scenario represented in our core object (SSP 2–4.5). With the result above, for example, we can conclude that there is a ~90% probability that the long-term average global warming will be between 2.0–4.0°C relative to pre-industrial reference.

2.8 Caveats to model weighting for probabilistic analysis

While we chose to provide the user with the flexibility to tailor our model weighting framework to fit the specifics of individual analyses, we encourage users to exercise caution and consider some caveats when using Matilda for probabilistic analyses. First, we provide several methods for weighting PPE members with different scoring criteria. It is important to note that decisions regarding the choice of scoring algorithm and criteria can be open-ended, and different approaches will yield different results. For example, likelihood weights computed using observed temperature alone will differ from those computed with observed CO₂ alone and weights computed using multiple lines of evidence (can be completed in Matilda using multi_criteria_weighting()) will potentially differ from those using either temperature or CO₂ individually. One approach to ensure credible results may include a sensitivity test of different weighting procedures. For example, we conducted a test in which we weighted a 1000-member ensemble using multiple lines of evidence (historical temperature and CO₂ concentration). In this test we show that altering the influence (weight) of the scoring criterion does not significantly change resulting warming probabilities. This would provide evidence that our results are not highly sensitive to changes in the influence of the scoring criterion (S1 Fig). Additionally, Pressburger et al. [13] provide an example of using multiple lines of evidence to apply filters to a PPE in an informative way.

Despite the similarities between the two scoring algorithms, the differences in approach and customization options can lead to variations in the behavior and performance of each method (Fig 4). It is therefore important to consider the specific goals of an analysis and characteristics of the data when selecting which scoring algorithm to use. The benefits of ramp scoring (score_ramp) is its simplicity, making it particularly useful for educational purposes or in preliminary analyses. As Pressburger et al. [13] show, this weighting method can be useful for methodologies aiming to cull ensemble members that fall outside of hard bounds. Bayesian scoring (score_bayesian) is a more complex approach that is conceptually closer to formal Bayesian calibration methods like BMA. The approximation of Bayesian methods allows for a more rigorous approach to model weighting. Unlike our ramp scoring approach, score_bayesian() has the ability to incorporate time-varying error, making it useful in cases where the uncertainty of observed data (scoring criterion) is expected to vary temporally.

Finally, it is important to note that score_bayesian() assumes independence of model-observation residuals. We acknowledge this assumption does not always hold true. Therefore, while our framework provides valuable results, there are limitations by not addressing more complex error structures, such as autocorrelation of residuals [44]. In such cases, more sophisticated modeling techniques may be necessary to improve accuracy of inferences from results [29,44].