Physical processes within geoscientific models are sometimes described by simplified schemes known as parameterisations. The values of the parameters within these schemes can be poorly constrained by theory or observation. Uncertainty in the parameter values translates into uncertainty in the outputs of the models. Proper quantification of the uncertainty in model predictions therefore requires a systematic approach for sampling parameter space. In this study, we develop a simple and efficient approach to identify regions of multi-dimensional parameter space that are consistent with observations. Using the Parallel Ice Sheet Model to simulate the present-day state of the Antarctic Ice Sheet, we find that co-dependencies between parameters preclude any simple identification of a single optimal set of parameter values. Approaches such as large ensemble modelling are therefore required in order to generate model predictions that incorporate proper quantification of the uncertainty arising from the parameterisation of physical processes.

The aim of any geoscientific model is typically to replicate the state and behaviour of real-world systems as accurately as possible, or at least with sufficient accuracy to generate useful insights into the problem being studied. This requires the model to incorporate a sufficiently accurate description of the real world, as well as sufficiently accurate data to provide boundary conditions and an initial state. However, geoscientific modelling inevitably involves making compromises in model design and implementation. Observational data, which is typically used to provide both initial and boundary conditions and to evaluate the models, can also be limited in accuracy and spatial coverage. Model error can therefore result from a number of sources, including missing or incomplete physics, missing or incomplete boundary conditions, and missing or incomplete initial conditions. This study focuses upon the first of these three potential sources of error, aiming to explore the contribution to model prediction error that arises from the simplifications made in the representation of physical processes within geoscientific models.

Either through choice or through necessity, physical processes are sometimes
described by simplified schemes known as parameterisations

There are a number of reasons why parameterisations are required, or why a
decision might be made to use them. These reasons include the following:

In practice, because of the time and effort required for parameter optimisation, the typical approach is simply to rely on prior published values. This can be the case even if these values were derived using a different version of the model or a different model configuration (for example, a different model resolution). This approach has the potential to result in inappropriate parameter choices, with the potential to introduce bias into model experiments. Using fixed parameter values in model ensembles will also mean that any uncertainties derived in the model predictions will not incorporate the contribution arising from the uncertainty in parameter values. This neglection of parameter uncertainty will result in an underestimation of the uncertainty in the model predictions.

Despite the challenges involved, there is an increasing emphasis on
quantifying uncertainty in geoscientific modelling. This is apparent, for
example, within the reports of the Intergovernmental Panel on Climate Change

The process of parameter optimisation is commonly referred to as “tuning”

The techniques used to optimise geoscientific models can be grouped into four
broad categories:

The iterative parameter optimisation process that we develop in this study
consists of five key steps:

Identify the parameters to be optimised.

Select the initial ranges for each parameter.

Construct and run a perturbed-physics ensemble.

Evaluate each member of the ensemble against observations and determine which regions of parameter space, if any, can be rejected.

Repeat steps 3 and 4 until the process has converged, i.e. until no further changes are made to the ranges of any of the parameters.

To illustrate the application of this process using a real-world geoscientific model, we will now present a demonstration using the Parallel Ice Sheet Model to simulate the present-day state of the Antarctic Ice Sheet (AIS).

In this study, we use version 0.7.3 of the Parallel Ice Sheet Model (PISM), a
highly parallel, open-source model that is suitable for large ensemble
modelling and the simulation of large-scale marine ice sheets. Computational
efficiency is achieved by employing a hybrid stress balance model to calculate
ice velocities. PISM is described and evaluated for Antarctica, by

The model is written in C

The following is an outline summary of the key model physics. Specific aspects that form the target of the parameter optimisation are discussed in further detail later in this section. This description also reflects the specific configuration of the model used in this study and may not therefore accurately describe the configurations used in other studies.

The model is configured to simulate the present-day state of
the entire AIS, with a horizontal resolution of

The two prognostic model variables
are ice temperature and ice thickness. A hybrid approach is used to calculate
the stress balance within the ice sheet, with the ice velocity being
calculated by the superposition of two shallow stress balance approximations

PISM assumes that the ice sheet rests on a layer of till

The lateral boundaries of the ice sheet are free
to evolve

The topography follows Bedmap2

Through reading prior published work and the model documentation, 10
parameters are selected for optimisation on the basis that they describe key
physical processes and that their values are not well constrained by either
theory or observations. Out of the 10 parameters, 6 relate to the description of
basal sliding, which is highly parameterised and involves parameters that are
particularly poorly constrained

The selection of initial parameter ranges is based on prior knowledge. PISM is
distributed with a number of pre-configured experiments

The 10 physical parameters that are varied in this study: the default values for the pre-configured SeaRISE-Antarctica experiment in PISM and the values used in previous studies. M11:

The shallow ice enhancement factor is a flow enhancement factor for the
non-sliding Shallow Ice Approximation, which is used to model grounded ice

Previous studies have used values of the shallow ice enhancement factor in the
range 1.0–5.0 (Table

The shallow shelf enhancement factor is a flow enhancement factor for the
shallow shelf approximation

Most previous studies have used values of the shallow shelf enhancement factor
in the range 0.4–0.8; however,

PISM uses a pseudo-plastic power law model to describe basal resistance

Theoretically, the value of the exponent must lie in the range

The yield stress

The effective till pressure scaling factor sets the value of

We are only aware of one previous study that has varied the effective till
pressure scaling factor parameter:

PISM uses the calving scheme of

The units of

The calving scheme in PISM removes any ice at the calving front that is
thinner than a specified minimum thickness

Previous studies have used values in the range 200–225

PISM calculates the till friction angle

Previous studies have used values for

Previous studies have used values for

A 100-member perturbed-physics ensemble is constructed. Each of the 10
parameters is perturbed independently, using a Latin hypercube approach

The specific set of parameter values used is generated using the

The model is run for 1 year in bootstrapping mode

The model is run for 100 years, with basal sliding disabled, to smooth the surface of the ice sheet.

The model is run for 250 000 years, with fixed geometry (basal sliding and mass conservation disabled), in order to improve the enthalpy field.

Finally, the model is run for 100 000 years with full physics. This period is chosen to allow the simulated ice sheet to come into equilibrium with the boundary conditions, and thus to lose the memory of its initial state.

Stages 1–3 are fast and typically account for only

For this study, PISM is run on a Huawei E9000 cluster with each ensemble
member using 224 cores. The queueing system for the cluster imposes a time
limit of 48

PISM uses an adaptive time-stepping scheme controlled by both the maximum diffusivity and the solutions to the mass conservation and energy conservation equations

Ensemble members can fail to complete the simulation because of numerical errors.

The convergence of the iterative parameter optimisation process for the 10 physical parameters varied in this study. For each iteration, the values shown are the minimum and maximum of the range sampled for that parameter. Bold text indicates a change in the value from the previous iteration. Also shown are the number of ensemble members to complete the simulation successfully during each iteration, as well as the critical values of the error metrics

Each ensemble member is evaluated to determine whether the simulated state of
the AIS is realistic. Two error metrics are calculated, as follows.

To determine whether each individual ensemble member can be considered to be
“realistic”, it is necessary to define critical values of

To determine whether part of the range for each parameter can be rejected, it
is necessary to perform a statistical test. For parameter

If the probability

Similarly, if the minimum value of

If

The value of

In this study

Table

During the optimisation process, the ranges for all four of the parameters
used to determine the till friction angle remain unchanged. However, for the
other six parameters, the ranges are reduced in width by between
14.5

The final ranges for two of the parameters are noteworthy within the context of previous work. We find that values of the shallow shelf enhancement factor smaller than 0.68 are not consistent with observations, or at least not for the specific version and specific configuration of the model used here. This excludes both the default value used by PISM (0.6) and the final values used in all previous studies with the exception of

For the final iteration of the optimisation process, 91 of the 100 ensemble
members complete successfully (i.e. they do not fail because of numerical
errors or because of exceeding the time limit). Of these members, by
definition, one-third satisfy the criterion

The simulated ice thickness for each of the 14 ensemble members is shown in
Fig.

The simulated ice thickness (m) for the 14 ensemble members that are in best agreement with observations. The observed ice thickness

The error in the simulated ice thickness (m) relative to Bedmap2

The values of each physical parameter for the 14 ensemble members that are in best agreement with observations. For the name of each parameter, see Tables

The Pearson correlation coefficients between each pair of the 10 physical parameters, for the 14 ensemble members that are in best agreement with observations. Bold text indicates that the correlation coefficient is statistically significant at the

The parameter combinations used in each ensemble member, as well as the values
of

Scatter plots of the relationships between the three pairs of physical parameters identified in Table

In each case, there are plausible physical explanations for the relationships:

Nonetheless, the existence of such relationships indicates that there is no single configuration of the model that can be considered to be optimal, at least on the basis of the evaluation conducted in this study.

We have shown that the optimal ranges for each parameter can be dependent on other variables. While we have been able to substantially reduce the volume of plausible parameter space, limitations on our understanding of the underlying physical system ensure that the plausible ranges remain large for some parameters. Using additional observational and palaeoclimate datasets to evaluate the model, such as the surface ice velocity and vertical profiles of temperature and age from ice cores, might allow us to constrain the parameter ranges further.

Nonetheless, to sample equally from amongst all plausible model configurations
requires the systematic sampling of parameter space. Different sampling
strategies might be superior to others, with the application of these
strategies (for example, the size of Latin hypercube ensembles) being
potentially constrained by computational considerations. The size of the
ensemble presented in this study (100) is relatively small, particularly given
the large number of parameters being optimised (10).

Fundamentally, however, the systematic sampling of parameter space requires
the application of large ensemble modelling approaches. We also emphasise the
importance of such modelling approaches when generating projections. Proper
ensemble design is necessary not just to quantify uncertainty around the mean
or median response of the system, but also to correctly identify the mean or
median response itself. The importance of these points is demonstrated by

Finally, we note that, whereas the approach developed in this study allows for the rigorous quantification of uncertainty in model parameters, and therefore the quantification of uncertainty in model projections, the technique presented here does not allow these uncertainty ranges to be interpreted in probabilistic terms. Extending our approach to generate future projections with associated probability distributions would require larger ensembles and further understanding of the uncertainties inherent in the physical system. This would include uncertainties in our physical understanding of that system, in the numerical representation of that physical understanding within the model and in the boundary conditions applied to the model.

We have developed a simple and efficient iterative technique for optimising
parameters in geoscientific models. Specifically, our approach is able to
eliminate regions of parameter space that are inconsistent with
observations. While it is analogous to other techniques that use large
ensemble modelling to refine parameter ranges

We have demonstrated the application of our technique by using PISM to simulate the present-day state of the AIS. After five iterations, we were able to refine the ranges of 6 out of 10 parameters. Most significantly, we find that multiple different parameter combinations are able to generate equally skillful simulations. This suggests that, at least for the model and for the experiments used in this study, ice sheet models have no single optimal configuration and therefore cannot be meaningfully “tuned”. Using single model configurations to generate predictions is therefore likely to underestimate the magnitude of the uncertainty around the best estimates. The solution to this is to use large ensemble modelling approaches with perturbed parameters.

Given the existence of parameter uncertainty, exploring parameter space is essential. However, the behaviour of the model may depend upon the parameter values in non-trivial (and, in particular, non-linear) ways: this requires systematic exploration of parameter space. Identifying implausible regions of parameter space allows for more efficient exploration of those regions that are potentially consistent with observations.

Correct values for geoscientific model parameters may be unknown or may not
even exist given that parameterisations by their very nature represent
simplifications of real-world processes. The parameter uncertainties
identified in this study, and in other studies that have used analogous large
ensemble modelling approaches

All the source code and data associated with this study are deposited at Zenodo:

The supplement related to this article is available online at:

SJP and JLR conceived the experiments. SJP conducted the experiments and analysed the output. All authors wrote the paper.

The authors declare that they have no conflict of interest.

This work was supported by the Centre for Southern Hemisphere Oceans Research, a joint research centre between the Qingdao National Laboratory for Marine Science and Technology (China) and the Commonwealth Scientific and Industrial Research Organisation (CSIRO; Australia). It was also supported by the University of Tasmania and the Australian Government's NCRIS program, through the Tasmanian Partnership for Advanced Computing. Development of PISM is supported by NSF grants PLR-1603799 and PLR-1644277 and NASA grant NNX17AG65G. The authors acknowledge use of the PyFerret program to generate the graphics in this paper. We thank Xuebin Zhang, two anonymous reviewers and the editor, Heiko Goelzer, for their comments, which have greatly improved the paper.

This research has been supported by the Australian Research Council (grant nos. SR140300001 and SR200100008).

This paper was edited by Heiko Goelzer and reviewed by two anonymous referees.