History matching for exploring and reducing climate model parameter space using observations and a large perturbed physics ensemble

Abstract

We apply an established statistical methodology called history matching to constrain the parameter space of a coupled non-flux-adjusted climate model (the third Hadley Centre Climate Model; HadCM3) by using a 10,000-member perturbed physics ensemble and observational metrics. History matching uses emulators (fast statistical representations of climate models that include a measure of uncertainty in the prediction of climate model output) to rule out regions of the parameter space of the climate model that are inconsistent with physical observations given the relevant uncertainties. Our methods rule out about half of the parameter space of the climate model even though we only use a small number of historical observations. We explore 2 dimensional projections of the remaining space and observe a region whose shape mainly depends on parameters controlling cloud processes and one ocean mixing parameter. We find that global mean surface air temperature (SAT) is the dominant constraint of those used, and that the others provide little further constraint after matching to SAT. The Atlantic meridional overturning circulation (AMOC) has a non linear relationship with SAT and is not a good proxy for the meridional heat transport in the unconstrained parameter space, but these relationships are linear in our reduced space. We find that the transient response of the AMOC to idealised CO₂ forcing at 1 and 2 % per year shows a greater average reduction in strength in the constrained parameter space than in the unconstrained space. We test extended ranges of a number of parameters of HadCM3 and discover that no part of the extended ranges can by ruled out using any of our constraints. Constraining parameter space using easy to emulate observational metrics prior to analysis of more complex processes is an important and powerful tool. It can remove complex and irrelevant behaviour in unrealistic parts of parameter space, allowing the processes in question to be more easily studied or emulated, perhaps as a precursor to the application of further relevant constraints.

Appendices

Appendix A: The statistical model derivation and justification for model discrepancy

We begin this technical Appendix A by deriving Eq. (7) using the statistical model in Sect. 4.1 given by Eqs. (5) and (6). Under our statistical model, supposing f ^h(x) is our climate model and is part of the MME and given the relation between the observations and the true climate to be as in (3), if x ₀ were $x^*$ we would have

$$\begin{aligned} z-\hbox{E}\left[f^{h}(x_{0})\right] &= (z-y)+(y-{{\mathcal{M}}}) + ({{\mathcal{M}}}-f^{h}(x_{0})) \\ &\quad + (f^{h}(x_{0})-\hbox{E}\left[f^{h}(x_{0})\right])\\ &= e + U - R_{h} + \xi(x_{0}),\\ \end{aligned}$$

The observation error, e, is uncorrelated with the other terms by definition, and U, the discrepancy between y and the expectation of the MME, and R _h , the discrepancy between f ^h and the expectation of the MME, are uncorrelated as a consequence of the representation theorem. We judge ξ(x ₀) to be uncorrelated with the other terms as this represents the error in our emulator based on a large ensemble at x ₀, which should have no relationship to these other errors. Hence, taking variances, Eq. (7) follows.

Our statistical model makes the explicit assumption that $f^{h}(x^{*}_{[h]})$, our climate model run at its best input is second order exchangeable with the MME. If the MME members are second order exchangeable, this would be fine so long as they represented different climate models run at their best inputs. However, this is not the case. Each MME member is a ‘tuned run’ of a different climate model run at some setting of its parameters x ^t_[j] where this setting is highly unlikely to be the best input, $x^{*}_{[j]}$, because, if it were, we would not be interested in the parametric uncertainty of the climate model. The usually adopted statistical model assumes y = f ^j $(x^{*}_{[h]})$ + η ^j, where we cannot learn about η ^j using our climate model (Rougier 2007). If x ^t_[j] were equivalent to $x^{*}_{[j]}$, under this model, a PPE tells you nothing more about climate than you already knew from $f^{j}(x^{*}_{[h]})$. Hence, usually, we would expect $f^{h}(x^{*}_{[h]})$ to be closer to y than $f^{h}(x^{*}_{[h]})$. The implication of this assumption is to have derived a larger discrepancy than is, perhaps, required.

To illustrate, suppose that we did view the climate models that contribute to the MME at their best inputs as second order exchangeable so that ${f^{j}(x^{*}_{[j]})={\mathcal{M}}\oplus R_{j}(f^{j}(x^{*}_{[j]}))={\mathcal{M}}\oplus R^{*}_{j}}$. Further, suppose that we have no a priori reason to believe that any of the MME members have been tuned closer to their best inputs than any other, so that

$$f^{j}(x^{t}_{[j]}) = f^{j}(x^{*}_{[j]})\oplus R^{t}_{j},$$

where R ^t _j is a mean-zero residual and Cov[R ^t _j ,R ^t _k ] = 0 for j ≠ k.

Then

$$f^{j}(x^{t}_{[j]}) = {{\mathcal{M}}} \oplus R^{*}_{j} \oplus R^{t}_{j},$$

which, by ignoring tuning as we have in (5), implies that R _j  = R ^* _j  + R ^t _j . This means that our discrepancy term has additional error due to tuning built in.

In practice it would be very difficult to assess the individual contribution of tuning error R ^t _j to the discrepancy without expert judgement from the climate modellers, even if the illustrative framework above were adopted.

However, by having a larger discrepancy we rule out less parameter space. This cautious approach is consistent with a first-wave history match, with the aim of removing those regions of parameter space that correspond to extremely unphysical climates. Though some of the remaining error can be attributed to tuning, this statistical model does not specifically account for the fact that the MME values we get from CMIP3 are not ten year means from climate models with a short spin up phase as our ensemble members are. This factor leads us to choose a statistical model that leads to larger discrepancy variance and less parameter space ruled out than might be obtained using expert judgement.

Appendix B: Constraint emulators

In order to fit each emulator mean function, we used a stepwise selection method to add or subtract functions of the model parameters to our vector g(x). The allowed functions were linear, quadratic and cubic terms with up to third order interactions between all parameters. Switch parameters were treated as factors (variables with a small number of distinct possible “levels”) and factor interactions with all continuous parameters were allowed.

First we perform a forward selection procedure where we permit each variable to be added to g(x) in its lowest available form. So if ct is not yet in g(x), ct can be added, but ct² cannot be added. If ct is already in g(x) then ct² can be added, but we simultaneously add all first order interactions with the other variables in g(x) so that the resulting statistical model will be robust to changes of scale (see Draper and Smith 1998, for discussion). We do similar for third order interactions. The term(s) added at each stage is the one that reduces the residual sum of squares from the regression the most.

When it becomes clear that adding more terms is not improving the predictive power of the emulator (a judgement made by the statistician based on looking at the proportion of variability explained by the regression and at plots of the residuals from the fit) we begin a backwards elimination algorithm. This removes terms, one at a time, with the least contribution to the sum of squares explained, without compromising the model. Lower order terms are not permitted to be removed from g(x) whilst higher order terms remain. We stop when the removing the next term chosen by the algorithm leads to a poorer statistical model. For more details on forwards selection methods and backwards elimination, see Draper and Smith (1998).

The emulators contained between 76 terms in g(x) (in the case of SAT) and 44 terms (in the case of PRECIP). For the SAT emulator, it was found that the response was easier to model as a function of log(entcoef), though this transformation produced inferior mean functions for each of the other 3 outputs. The terms for the SAT emulator can be seen in Table 3. Each header in the table refers to one of the parameters from Tables 5 and 6, shortened in an obvious way to save space. Numbers on the diagonal refer to power terms in g(x) in each of the relevant parameters. The number 1 on the diagonal implies only a linear term was included in g(x). The number 2 implies that both quadratic and linear terms were in g(x). The number 3 implies cubic, quadratic and linear terms.

Table 3 A table indicating which terms are in g(x) for our emulator of SAT in Eq. (1)

Table 4 Parameter descriptions and model section

Table 5 Ranges for each of the continuous parameters varied in our 10,000 member ensemble

Table 6 Switch parameters and their settings in our 10,000 member ensemble

Numbers on the upper triangle refer to the inclusion or not of interactions between the two relevant variables. For example, reading from the table, the term (vf1*dyndiff) is included in g(x), but the term (vf1*g0) is not. Variables indicated in bold on the lower triangle refer to three way interactions that are present in g(x). For example, the terms (ent*vf1²) and (vf1*ent²) are both included in g(x). Note that dynd and r_lay are both handled by the model as factors with pre-specified levels. Our emulator does contain a constant term, so the vector g(x) includes the element 1.

Having fitted a mean function, the next step in emulation is to provide a statistical model for the residual $\epsilon(x)$. Upon investigation, the residuals from the fitted mean functions appeared to behave like white noise. Hence we decided that it was not worth modelling $\epsilon(x)$ in Eq. (1) as a weakly stationary Gaussian process (where $\epsilon(x)$ and $\epsilon(x')$ are correlated via a covariance function R(|x − x′|) that depends on the distance between points in parameter space (See, for example, Williamson et al. 2012). We therefore opted to model $\epsilon(x)$ as mean zero uncorrelated error as is done by Sexton et al. (2011) and Rougier et al. (2009). We let the variance of $\epsilon(x)$ be equal to the variance of the residuals from the four regressions.

To validate the emulators, we kept 100 randomly chosen ensemble members back as a validation set. The emulators were not trained using these members. We then predict the value of the constraints in each ensemble member using the emulator and examine the residuals (the difference between our prediction and the truth). The residuals are plotted against the predicted values in Fig. 9.

Fig. 9

Validation plots for each of the 4 emulators fitted to the observational constraints. For each variable stated in the plot title, the predicted value of the emulator for that variable (x-axis) is plotted against the difference between the truth and the emulator prediction for the data in the validation set. The red dashed lines are ±3 standard deviations of the predictions. The validation set was not used when building the emulators

The emulator’s performance is acceptable as Fig. 9 shows that our emulator predictions are consistent with our uncertainty specification. If anything, the validation suggests that we have over specified our variance (at least in the case of SCYC and PRECIP), because no points fall outside 3 standard deviations. However, this will only lead to less parameter space being ruled out, and we should be cautious in terms of what parts of parameter space are ruled out and what is retained. We make the four emulators available to download along with our implausibility code so that our results may be replicated, and their sensitivity to any of our judgements explored.

Appendix C: Tables

Tables 4, 5 and 6 give descriptions and ranges for the parameters and the settings of switches used in our ensemble. Some parameters have relationships with other model parameters that were given to us by the Met Office so that a change in one leads to a derivable value for the other. CWland also determines CWsea, the cloud droplet to rain threshold over sea (kg/m³), MinSIA also determines dtice (the ocean ice diffusion coefficient) and k_gwd also determines kay_lee_gwave (the trapped lee wave constant for surface gravity waves m^3/2).