Improved Representation of Ocean Heat Content in Energy Balance Models

TL;DR

This paper improves energy balance models by incorporating ventilation processes, enabling better representation of surface warming and vertical heat distribution in the ocean, highlighting the importance of ocean circulation in climate modeling.

Contribution

It introduces a modified anomaly-diffusing energy balance model that includes ventilation effects, enhancing the accuracy of heat uptake representation.

Findings

The modified model accurately captures surface warming and heat distribution.

A simple four-layer model suffices for effective representation.

Ventilation processes are crucial for vertical heat uptake modeling.

Abstract

Anomaly-diffusing energy balance models (AD-EBM) are routinely employed to analyze and emulate the warming response of both observed and simulated Earth systems. We demonstrate a deficiency in common multi-layer as well as continuous-diffusion AD-EBM variants: They are unable to, simultaneously, properly represent surface warming and the vertical distribution of heat uptake. We show that this inability is due to the diffusion approximation. On the other hand, it is well understood that transport of water from the surface mixed layer into the ocean interior is achieved, in large part, by the process of ventilation---a process associated with outcropping isopycnals. We, therefore, start from a configuration of outcropping isopycnals and demonstrate how an AD-EBM can be modified to include the effect of ventilation on ocean uptake of anomalous radiative forcing. The resulting EBM is able to successfully represent both surface warming and the vertical distribution of heat uptake, and indeed, a simple four layer model suffices. The simplicity of the models notwithstanding, the analysis presented and the necessity of the modification highlight the role played by processes related to the down-welling branch of global ocean circulation in shaping the vertical distribution of ocean heat uptake.

Tables

Table 1. Table 1: A brief overview of the experiments conducted. Units for RMS error of SAT is degrees C, and units for RMS error of OHC is 10 22 superscript 10 22 10^{22} J.

Expt. #	SCM	Calibration Details			RMS Error: SAT				OHC
		Forcing	AOGCM QoI	Obs. Data	4x	1%	HIS	OBS	$\le$700	$>$700
1	2L	4x	SAT		0.13	0.11	0.16	0.18	7.5	3.2
2	2L	4x	SAT, RN		0.14	0.12	0.16	0.18	8.5	3.1
3	2L	4x	SAT, RN		0.21	0.14	0.16	0.18	5.8	1.5
		Historical		0-700m OHC	Degraded 4x SAT
				700-2000m OHC	Improved but Poor OHC
4	4L	4x	SAT, RN		0.17	0.13	0.16	0.18	5.7	1.4
		Historical		0-700m OHC	Reasonable SAT
				700-2000m OHC	Poor OHC
5	60L	4x	SAT, RN		0.18	0.13	0.16	0.18	6.5	1.7
	(many variants)	Historical		0-700m OHC	Reasonable SAT
				700-2000m OHC	Poor OHC
6	4LE	4x	SAT, RN		0.18	0.13	0.16	0.18	1.4	1.5
		Historical		0-700m OHC	Reasonable SAT
				700-2000m OHC	Reasonable OHC

Table 2. Table 2: Some of the continuous AD-EBM variants considered

Expt. #	SCM	Comment
5a	40L; MLD inferred
	${\kappa }_v^{TC}$ inferred	Degraded 4x SAT
	${\kappa }_v^{M{L}_{bot}}={\kappa }_v^{TC}$	Poor OHC
5b	40L	Degraded 4x SAT
	MLD = 70 m	Poor OHC
5c	40L	Degraded 4x SAT
	${\kappa }_v^{M{L}_{bot}}\ne {\kappa }_v^{TC}$
	${\kappa }_v^{M{L}_{bot}},{\kappa }_v^{TC}$ inferred	Poor OHC
5d	40L	Poor SAT
	${\kappa }_v^{TC}={10}^{-5}$ m2/s	Poor RN
	${\kappa }_v^{M{L}_{bot}}\ne {\kappa }_v^{TC}$	Poor OHC
5e	40L; w inferred	Same as in Expt. 5c
	${\kappa }_v^{M{L}_{bot}}\ne {\kappa }_v^{TC}$	Degraded 4x SAT
	${\kappa }_v^{M{L}_{bot}},{\kappa }_v^{TC}$ inferred	Poor OHC
5f	40L; ${\kappa }_v^{>650m}$ inferred	Same as in Expt. 5c
	${\kappa }_v^{M{L}_{bot}}\ne {\kappa }_v^{TC}$	Degraded 4x SAT
	${\kappa }_v^{M{L}_{bot}},{\kappa }_v^{TC}$ inferred	Poor OHC
5g	60L; Same as 5c	Same as in Expt. 5c
7a	40L; Same as 5c	Reasonable SAT
	Non-local interaction	Reasonable OHC
7b	60L; Same as 5c/g/7	Same as in Expt. 7a
	Non-local interaction

Table 3. Table 3: RMS error of SAT, with respect to CCSM4 experiments abrupt4x, 1%, historical forcing, and with respect to HadCRUT, and RMS error of OHC with respect to 0-700m and 700-2000m estimates from NODC in historical-forcing runs. Units for RMS error of SAT is degrees C, and units for RMS error of OHC is 10 22 superscript 10 22 10^{22} J.

Expt.	SAT RMS Error				OHC RMS Error
	4x	1%	Hi	HC	0-700m	700-2000m
5a	0.19	0.11	0.16	0.17	7.5	2.2
5b	0.26	0.12	0.16	0.16	8.2	2.3
5c	0.18	0.13	0.16	0.18	6.5	1.7
5d	0.39	0.14	0.15	0.17	7.0	3.1
5e	0.16	0.13	0.16	0.18	5.9	1.4
5f	0.17	0.13	0.16	0.18	6.5	1.7
5g	0.18	0.13	0.16	0.18	6.5	1.7
7a	0.15	0.12	0.16	0.18	2.0	1.5
7b	0.15	0.12	0.16	0.18	1.6	1.5

Equations

$$C_u\frac{d{T}_u}{dt}$$

$$C_d\frac{d{T}_d}{dt}$$

$$\frac{d\mathcal{H}}{dt}=\mathcal{F}-\lambda T_u;\mathcal{H}=c(h_u{T}_u+h_d{T}_d)$$

$$c\left({\alpha }_i\sum _{j=1}^{i-1}{h}_j+h_i\sum _{j=1}^i{\alpha }_j\right)\frac{d{T}_i}{dt}={\alpha }_i\left(\mathcal{F}-\lambda T_i\right)+\sum _{j=\max (1,i-1)}^{\min (N,i+1)}{\gamma }_{ij}(T_j-T_i)$$

$$\frac{dT}{dt}=C+AT$$

$$\sum _{j=i+1}^N{\alpha }_j{T}_j+T_i\sum _{j=1}^i{\alpha }_j={\tilde{T}}_i$$

$$BT=\tilde{T}$$

$$T=B^{-1}\tilde{T}.$$

$$\frac{d\tilde{T}}{dt}=\tilde{C}+\tilde{A}\tilde{T}.$$

$${\tilde{A}}_{ij}=\sum _{k=1}^n\sum _{l=1}^n{B}_{ik}{A}_{kl}{B}_{lj}^{-1}=\sum _{k=i}^n\sum _{l=1}^j{B}_{ik}{A}_{kl}{B}_{lj}^{-1}$$

$$T_u^{ESM}(t)=T_u^{EBM}(t)+{ϵ}_t,$$

$$\Sigma (t-s)={\sigma }^2{\rho }^{|t-s|}$$

$$\frac{\partial T}{\partial t}+w\frac{\partial T}{\partial z}=\frac{\partial }{\partial z}\left(\kappa (z)\frac{\partial T}{\partial z}\right)+{\gamma }_d\exp \left(-\frac{(z-z_{mx}{)}^2}{\Delta z^2}\right)(T_1-T).$$

Full text

††journal: Climatic Change††journal: Climatic Change

1 Introduction

Energy balance models (EBM) have long been used in climate science as conceptual tools to understand heat transfer within the Earth system (e.g., Budyko, 1969; Sellers, 1969; Hoffert et al, 1980; Schneider and Thompson, 1981; North et al, 1981; Murphy, 1995; Gregory, 2000), as computationally-efficient emulators of and diagnostic tools for more complex Earth system models111e.g., because of their computational intensity, ESMs themselves are never run long enough to directly assess climate sensitivity. (ESMs) (Wigley and Raper, 1992; Raper and Cubasch, 1996; Houghton et al, 1997; Raper et al, 2001; Meinshausen et al, 2011a; Geoffroy et al, 2012), and to explore and quantify uncertainties in climate change projections (Urban and Keller, 2010; Padilla et al, 2011; Aldrin et al, 2012; Johansson et al, 2015; Bodman and Jones, 2016). EBMs represent changes in climate as the transfer of heat between appropriately aggregated units, and where the units are loosely referred to as layers when they are distributed in the vertical dimension alone (cf. isopycnal layers) and as boxes if other dimensions are involved as well. As an example of the former, in two-layer models that are widely used (e.g., see Schneider and Thompson, 1981; Gregory, 2000; Held et al, 2010, and references therein) both for climate projection purposes and to emulate Atmosphere-Ocean General Circulation Models (AOGCM), the atmosphere, the land surface, and the oceanic mixed layer components of the earth system are aggregated into one surface layer with a smaller thermal inertia and it interacts with a second layer that is representative of the sub-surface ocean and which has a larger thermal inertia. On the other hand, EBMs may allow for a coninuous representation in the relevant dimension (e.g., Budyko, 1969; Sellers, 1969; Raper et al, 2001; Meinshausen et al, 2011a, and others). Thus, in the context of upwelling-diffusion or pure-diffusion in the vertical, such EBMs consist of many layers (typically $\sim$40) with the intent of representing the limit of a continuous vertical dimension. Further combinations are easily realized.

In its use as a physics-based emulator, we would ask that an EBM be able to reproduce both ESM-simulated and observed surface and subsurface warming in a range of scenarios. Diffusion-class models have shown skill at reproducing observed surface warming and overall ocean warming as well as warming behaviors of ESMs under historical and future CO2 emissions scenarios (e.g., see Gregory, 2000; Raper et al, 2001, 2002, and others). In particular, a popular form of a two-layer model has proven to be capable of capturing the typical two-timescale surface temperature response seen in abruptly forced AOGCM simulations (e.g., see Raper et al, 2001; Friend, 2011; Geoffroy et al, 2013, and others) across a range of AOGCMs. However, how well that model can represent the vertical structure of ocean heat uptake has not received much attention. We find that the two-layer model is unable to simultaneously properly represent surface warming responses and observational estimates of 0-700m and 700-2000m heat storage. Perhaps more surprising is that we find that more highly-resolved diffusion-class models are similarly deficient.222In order to exclude the possibility that our results are due to poor parameter tuning in the EBMs we consider, we apply a Bayesian statistical framework and Monte Carlo sampling to widely explore the space of parameter uncertainties.

On the other hand, from the point of view of dynamical processes underlying global ocean circulation, this last finding—that vertically-resolved EBMs are unable to simultaneously fit both the surface response and the vertical distribution of ocean heat uptake—is not surprising for the following reasons: (a) From the point of view of ocean dynamics, it is well understood that transport of water from the surface mixed layer into the ocean interior is achieved, in large part, by the process of ventilation—a process associated with outcropping isopycnals. That is, when isopycnals outcrop at high latitudes, water masses that are at the surface at high latitudes are able to move laterally back to middle and low latitudes at deeper levels (e.g., see Chapter 4 of Pedlosky, 2013). (b) The deep ocean loses heat to the atmosphere at high latitudes through convective instabilities (e.g., see Chapter 7 of Talley, 2011), and when the stability of the water column is incrementally enhanced in a warming scenario (that is, the water column continues to remain unstable, but its degree of instability is slightly reduced), the heat loss from the deep ocean is slightly reduced leading effectively to a warming of the deep. Indeed, the effects of such processes have been verified in various climate models (e.g., see Gregory, 2000; Raper et al, 2001, 2002, and others).

Clearly, while there are many directions in which one could improve the realism of an EBM’s warming response, in this work we present one simple modification that enables simultaneously fitting a wide range of scenarios: Outcropping of isopycnals/neutral surfaces is a common feature of global ocean circulation (e.g., see Talley, 2011; Pedlosky, 2013). Nevertheless, this configuration has not received much attention in the context of EBMs. Indeed, we find that a cure to the problem identified in AD-EBMs lies in considering such a configuration. We, therefore, start from a conceptual configuration of outcropping isopycnals. Next, we perform a transformation of coordinates in order to cast this configuration into a framework of horizontally-averaged layers, a framework that is more commonly encountered in the EBM context. We then show mathematically that when such a configuration of outcropping isopycnal/neutral surfaces is viewed in terms of horizontally-averaged layers, non-local interactions between near-surface and deeper ocean layers result. We then show that such a modification, a modification that represents the lowest order effects of outcropping isopycnal/neutral surfaces, is sufficient to simultaneously capture surface warming and deeper ocean heat uptake. Indeed, we find that it is not necessary to use a fully one-dimensional model; a four-layer model suffices.

2 Methodology and Results

2.1 The Two-Layer EBM

A particular form of simple climate model (SCM) that has been popular in summarizing integral thermal properties of AOGCMs and ESMs is the anomaly-diffusing energy-balance model (AD-EBM) (e.g, see Schneider and Thompson (1981); Gregory (2000); Held et al (2010)). In this horizontally-integrated model of the climate system, the upper (surface) layer is comprised of the oceanic mixed layer, the atmosphere and the land surface and the bottom layer represents the ocean beneath the mixed layer. Evolution of the upper surface heat content anomaly per unit area is given by

[TABLE]

where ${\cal F}$ represents the effective (anomalous) radiative forcing (e.g., as due to green-house gases recently), $\lambda$ represents the climate feedback parameter (e.g., see Knutti and Hegerl (2008)). Exchange of heat between the surface layer and the ocean beneath is parameterized by the difference in temperature between the two layers, with $\gamma$ being the constant of proportionality. Similarly, evolution of the subsurface ocean heat content (OHC) anomaly (again per unit area) is given by

[TABLE]

The heat capacities $C_{u}$ and $C_{d}$ can be equivalently expressed in terms of ocean layer depths, $h_{u}$ and $h_{d}$ as $C_{(\cdot)}=ch_{(\cdot)}$, where $c$ is the volumetric heat capacity with a nominal value of 4.1$\times$ 106 Jm*-3K-1*.

2.2 Skillful Representation of Surface Warming in the Two-Layer EBM

This two-layer model has been used extensively for estimating climate sensitivity of AOGCMs and ESMs (e.g., see Held et al (2010); Geoffroy et al (2012), and others). Figure 1 shows an example of its typical use. In this figure, we have calibrated the four EBM parameters ($\lambda$, $\gamma$, $h_{u}$, $h_{d}$) using the Community Climate System Model v. 4.0 (CCSM4) AOGCM’s annual global surface air temperature (SAT) output from an abrupt CO2 quadrupling (abrupt4x) experiment333We note that the intention of the CMIP abrupt 4$\times$CO2 forcing experiment is one of calibration and diagnosis of equilibrium climate sensitivity since post-industrial surface warming itself is compatible with a wide range of climate sensitivities as the sole calibration data. The parameters are inferred probabilistically from the CCSM4 SAT time series using Bayesian Monte Carlo sampling, assuming a first-order autoregressive noise model to capture the AOGCM natural variability not represented in the EBM (e.g., as in Urban and Keller (2010) or Nadiga and Urban (2016); see online supplement section titled “Further Methodological Details” for further details). Bayesian inference allows us to explore the full set of EBM parameters that are capable of fitting the AOGCM SAT data, to ensure that our results are not fragile with respect to one particular choice of parameters.

Next, we validate the calibration of the two-layer model by considering new forcing scenarios—1%/year CO2 increase and historical forcings. By propagating the, calibrated, joint posterior distribution of parameters, we produce posterior predictive distributions of the respective SAT responses. The SAT responses and their 95 percentile ranges are shown in Fig. 1 (see Sec. A of ES for further details). The good correspondence between the posterior predictive distribution of SAT in the 1% CO2 (RMS error of 0.11 C) and historical forcing experiments (RMS error of 0.16 C), as compared to the corresponding AOGCM results attest to both the reasonableness and the usefulness of the model. The SAT response of a handful of other AOGCMs in the abrupt4x and 1% CO2 forcing CMIP5 (Coupled Model Intercomparison Project Phase 5; Taylor et al (2012)) experiments (lighter lines) and the Hadley center estimate (HadCRUT4) are shown for future reference.

2.3 Ocean Heat Uptake Considerations

Use of SCMs like in the previous example is well established. However since in excess of 90% of the anomalous radiative forcing is sequestered in the world oceans (e.g., Domingues et al, 2008; Levitus et al, 2012; Balmaseda et al, 2013, and others) and since the nature of ocean heat uptake is important in determing surface response, we next examine the ability of the two-layer EBM to represent aspects of the ocean heat uptake. For this, the main data product we consider is the sequestration of heat in the top 700 m and between 700 m and 2 km in the global ocean, as provided by the National Oceanic Data Center (NODC) (Levitus et al, 2012). We compare these two observation-based timeseries against the corresponding timeseries in the SCMs.

We have considered a specific AOGCM, CCSM4, for our example here and it’s clear that we could have chosen any other AOGCM from the CMIP5 model ensemble. While the particular choice of the AOGCM whose response we use to calibrate the SCM will affect details of the results presented, we believe that the importance of this choice is secondary. Figure 1 shows in light lines the spread of SAT in the abrupt4x and 1% CO2 CMIP5 experiments from a handful of other modeling centers as well. While there is considerable spread in the response of the different models in these experiments, the spread in the averaged response of the different models in the historical forcing simulations is much less given the required tuning of AOGCMs to best match the current state of the Earth’s climate system (e.g., Mauritsen et al, 2012). Multi-model AOGCM differences are therefore not likely to account for the ability of an SCM, fit to a particular AOGCM, to reproduce historical climate trends.

In the top-left panel of Fig. 2, we compare the 0-700m (in the main plot) and 700-2000 m OHC (in the inset plot) in the previous experiment444All experiments are run with each of the three forcings under consideration (abrupt4x, 1%, and historical), and the quantities of interest (QoI) are the SAT with each forcing and the 0-700m and 700-2000m OHC with historical forcing. While we also considered RN, those results are left out for brevity. against the NODC (Levitus et al, 2012) observational estimates; again the posterior mean and 95% envelopes are shown.555In the two-layer EBM, warming of the ocean extends only to the depth of the two layers. That is, there is no warming below the depth of the second layer. Furthermore, the depth of the top layer should be thought of as roughly equivalent to the globally-averaged mixed-layer depth (O(100m); $\ll 700m$) The 0-700 m and 700-2000 m warming is obtained using linear interpolation. The SCM estimates are seen to fit the observational estimates very poorly.

This may not be surprising considering that in the previous experiment there was no constraint on heat storage in the climate system; as shown in Table. 1, calibration data comprised of only the surface warming in the abrupt4x AOGCM experiment. After all, the main motivation and the typical use of the two-layer EBM has been to capture the two timescales that are thought to characterize the surface warming response (e.g, see Gregory (2000) or Held et al (2010)).

2.4 Constraining the EBM’s Sub-surface Heat Content

Summing Eqs. 1 and 2 results in

[TABLE]

and Eq. 3 represents the evolution of the heat content anomaly per unit area. The right hand side of Eq. 3 is the radiative imbalance (RI) or equivalently radiative non-equilibrium (RN).

Given the poor fit of OHC when the two-layer EBM is calibrated using the SAT response of the EBM in the abrupt4x experiment, we repeat the previous experiment, but now calibrating against both the SAT and radiative imbalance responses of the AOGCM. (A list of the experiments considered here is shown in Table 1; a few other supporting experiments discussed in the online supplement are shown in a table there.) While the SAT behavior is similar to that in Fig. 1 (not shown, but see Sec. 2 of supplement for figure and details), the OHC comparison is shown in the top-right panel of Fig. 2. It is seen that constraining the total heat storage leads to only minor improvements in the two layer EBM’s representation of ocean heat uptake. However, from the point of view of calibration itself, the differences are more dramatic: Not only are the EBM parameters inferred to be significantly different in the two experiments, but the inferences themselves are sharper and correlations that exist between model parameters in experiment 1 are seen to be greatly reduced in experiment 2 (see Sec. 2 of supplement for details).

2.5 A Shortcoming of the Two-Layer EBM

Since constraining overall heat storage in the abrupt4x CO2 forcing experiment was seen to be insufficient to elicit an observationally-consistent heat storage response in the historical forcing experiment, we consider experiment 3, wherein we calibrate the two layer EBM against the AOGCM’s SAT and RN responses in the abrupt4x CO2 forcing experiment, and calibrate against the observational estimate (0-700m and 700-2000m) of OHC itself, when the SCM is forced with an estimate of the historical forcing followed by the RCP8.5 scenario (Meinshausen et al, 2011b) to 2050.

In experiment 3, the simple two layer EBM’s representation of OHC (bottom left panel of Fig. 2) is seen to be better with improvements in the representation of the 700-2000m heat content. However, the 0-700m warming continues to be poorly represented. Attendant with the partial improvement in the representation of heat storage, unsurprisingly, is a slight degradation in the ability of the EBM to fit the SAT response in the abrupt4x CO2 experiment (see Sec. 2 of supplement for further details and discussion).

We consider the inability of the two layer model to simultaneously properly represent the surface warming and the 0-700m and 700-2000m heat storage a shortcoming of the model and seek a simple extension of the model that will not suffer from this shortcoming. For brevity, when we refer to the problem or the shortcoming in the rest of the article, we are specifically referring to the shortcoming just described.

2.6 Persistence of the Shortcoming in Vertically-Resolved AD-EBMs

To investigate the possibility that it is the lack of resolution in the vertical that is the cause of the inability of the two layer EBM to represent heat storage in a manner that is consistent with observations, we considered a three layer extension and a four layer extension (experiment 4) of the two layer model. In each of these cases, calibration proceeded as in experiment 3. The close similarity in the SAT and heat storage responses of the three (not shown) and four layer models (see Table 1) to that in the two-layer suggests that a lack of resolution in the vertical may not be the direct cause of the problem.

Nevertheless, to probe the issue of vertical resolution further, we considered a number of variants of the continuously stratified version of the AD-EBM. A few of these experiments are detailed in the supplement under the section titled “Experiments with Continuous AD-EBM Variants”. (Bryan and Lewis, 1979; Nadiga and Urban, 2016). In these experiments, we typically discretized the subsurface thermocline region down to 2500 m into a series of layers of equal thickness since observational estimates of warming below 2000 m in the global ocean tend to be small (e.g., Cheng et al, 2017): $<$10% of the full ocean warming or $\lesssim$ 3$\times$ 1022 J currently. In these experiments, the main variations were in the form of vertical or diapycnal diffusivity assumed, and as to whether the coefficient(s) of diffusivity itself was (themselves were) fixed at realistic values or inferred in the calibration procedure. Other variations included fixing or inferring the depth of the upper layer and allowing for a (depth-independent) upwelling velocity or not (see supplement for details and results).666When a number of the variations mentioned are considered in combination, and with a large number of layers, calibration of the resulting model can sometimes lead to, e.g., the top layer depth being different from the a priori estimate of around 70m. In such cases, we performed a companion experiment in which the top layer depth is fixed at 70m, and so on. Like in the other panels of Fig. 2, the bottom-right panel shows the OHC comparison for one such model that fit AOGCM responses and observational estimates of OHC best (Experiment 5g in supplement).

The improvements in going from a two-layer model to a large number of layers (here 60) is seen to be rather small. The inability of the continuously stratified AD-EBM may be understood as follows: In the suite of experiments conducted with the continuous model (and described in Sec. 3 of supplement), in the cases which fit AOGCM responses and observational estimates of OHC best, the inferred value of thermocline diffusivity ranged between 0.7 and 0.9 cm2/s. These values are much higher than the accepted levels of 0.1 cm2/s. One reason for these large values of vertical diffusivity is that the model is attempting to get sufficient heat down to the 700-2000 m depths in order to satisfy the observational constraint. And in trying to do so, it is clear from these experiments that this is causing too large a heating in the 0-700 m depth range—a feature that is consistent across a wide suite of experiments (described in Sec. 3 of supplement). That is, assuming that vertical diffusivity is not varying with time (and there is no reason to expect it to have varied substantially with time in the post industrial period) getting observed levels of warming into the 700-2000 m layer while limiting the warming of the 0-700m layer to observed levels is difficult to achieve in AD-EBMs.

2.7 Moving Beyond Diffusion of Anomalies: Outcropping Surfaces

The inability of the continous variants of AD-EBM model to reasonably represent the vertical distribution of OHC leads us to next consider if it is the anomaly diffusion aspect of these models that prevents the SCM from being able to simultaneously properly represent surface warming and the 0-700m and 700-2000m heat storage. After all, the global ocean circulation tends to be highly dynamic in the sense that it is largely adiabatic outside of the surface layer and bottom boundary layers, and an anomaly diffusing model attempts to parameterize the horizontally-averaged effect of the myriad dynamical processes through diffusion across adjoining layers. This may be restrictive in the sense that the effect of processes such as deep-water formation, and subduction and mode water formation processes on ocean heat storage, when averaged horizontally, may no longer appear as local diffusion operators. For example, when water in the midlatitudes whose properties are set by its last interaction with the atmosphere slides down an appropriate isopycnal, in a horizontally averaged sense this process will appear as if two non-adjoining layers are interacting.

To demonstrate this effect, consider a schematic of the meridional-depth cross-section of the world oceans, as shown in Fig. 3. In this figure, the outcropping isopycnal/neutral surfaces in thin continuous lines are idealized into the horizontal and vertical surfaces shown in thick continuous lines. In Fig. 3, $\alpha_{i}$ is the meridional fraction of layer $i$ that outcrops and is exposed to the atmosphere. The EBM for such layers is:

[TABLE]

and where $\gamma_{ij}$ ($=\gamma_{ji}$) parameterizes the interactions between adjoining layers $i$ and $j$ and the summation terms drop to 0 when the starting lower index is greater than the ending upper index.

With trivial manipulation, Eq. 4 can be rewritten in a compact symbolic form as

[TABLE]

where $T$ is the vector of temperatures $T_{i}$ in Eq. 4, $C$ is a constant vector related to the anomalous radiative forcing and $A$ is a tridiagonal matrix that is independent of $T$. (Note that $A$ in Eq. 5 is not symmetric even though the right hand side of Eq. 4 is. This is because each row of the matrix form of Eq. 4 gets divided by a different constant (the coefficient of $\frac{d{T}}{dt}$) in Eq. 4.) Next, consider horizontal averaging of temperature in Fig. 3 to go from the idealized outcropping layers (shown in thick continuous lines) to horizontally-averaged layers (shown in dashed lines):

[TABLE]

where $\widetilde{T}_{i}$ is the temperature of the $i^{th}$ horizontally-averaged layer, and again the same rules hold for summation. Again, Eq. 6 can be written in a compact form as:

[TABLE]

where by virtue of the limits of the summation in the first term on the left hand side of Eq. 6, $B$ is seen to be upper triangular. In the generic case of $\sum_{j=1}^{i}\alpha_{j}\neq 0$, Eq. 7 can be inverted as

[TABLE]

$B^{-1}$ is upper triangular because $B$ is upper triangular. Inserting Eq. 8 in Eq. 5 and simplifying leads to

[TABLE]

While $A$ is tridiagonal (see Eq. 4), $\widetilde{A}$ is not: To see that $\widetilde{A}$ is not tridiagonal, consider that

[TABLE]

where the second equality is due to the upper triangular nature of $B$ and ${B}^{-1}$ ($B_{ij}=B_{ij}^{-1}=0$ for $i>j$). Therefore if $i>j+1$, the values of $k$ and $l$ in the above sum satisfy the inequality $k\geq i>j+1\geq l+1$, meaning that $k>l+1$. Therefore $A_{kl}=0$ for all the terms in the above sum for $i>j+1$, making ${\widetilde{A}}_{ij}=0$ for $i>j+1$.

That is while the outcropping isopycnal/neutral layers interact purely locally in isopycnal/neutral coordinates (tridiagonal $A$ in Eq. 5), the same dynamic when represented using horizontally averaged layers leads to nonlocal interactions (e.g., the top layer interacts with as many of the underlying layers as there are outcropping isopycnals.)

2.8 Results with the Modified EBM

Following the mathematical development above, we next evaluate if the modified EBM given by Eq. 9 is capable of overcoming the shortcoming identified in the AD-EBMs. As previously mentioned, $\widetilde{A}$ not being tridiagonal amounts to allowing for interactions between non-adjoining, horizontally-averaged layers. As for adjoining layers, we again parameterize such interactions in terms of the difference in their temperature anomalies. The simplest extension then consists of a three layer model with a new interaction between the top layer and the deepest layer. Within the range of experimentation that we conducted, we find that such a three layer extension is incapable of adequately resolving the shortcoming; for that we need at least four layers.

Figure. 4 shows the response of such an extended four layer model (Expt. 6 in Table. 1). Both SAT and OHC representations are reasonable enough to suggest that this is a simple extension of the two layer EBM that overcomes the shortcoming of the latter of not being able to simultaneously represent the surface responses and subsurface ocean heat uptake in a variety of settings.

3 Discussion and Conclusions

Simple climate models play a valuable role in helping interpret both observations and responses of comprehensive earth system models. To this end, we considered the anomaly-diffusion energy balance model and used it in a Bayesian framework to interpret a few CMIP5 experiments—the abrupt 4$\times$CO2 forcing, the 1% CO2 forcing and the historical forcing experiments. The intention of the CMIP abrupt 4$\times$CO2 forcing experiment is one of calibration and diagnosis of equilibrium climate sensitivity since post-industrial surface warming itself is compatible with a wide range of climate sensitivities. As such, we used the abrupt 4$\times$CO2 forcing CCSM4 experiment in conjunction with an observational estimate of OHC to infer the SCM parameters and used those parameters to compute the SCM responses to other forcings. We demonstrated that the commonly-used anomaly-diffusing energy balance SCM is deficient in being unable to simultaneously properly represent the SAT responses and observational estimates of OHC from NODC.

After verifying that this deficiency is not related to vertical resolution, we considered the possibility that it is due to the anomaly diffusion approximation. To this end, we demonstrated mathematically that when outcropping isopycnals are considered in a framework of horizontally-averaged layers, the outcropping of isopycnals leads to non-local interactions in the framework of horizontally-averaged layers. With this modification, we then showed that four layers were sufficient to successfully resolve the shortcoming of the AD-EBM class of SCMs. We also go on to show that that with such apparent non-local interactions—see supplement section “Further Experiments with Non Local Interactions” for a particular form—a continuous version of the model is similarly able to successfully resolve the shortcoming.

It is possible, however, that other kinds of augmentations to the commonly-used AD-EBM may fix the identified deficiency as well. On the other hand, since the simple extension that allows for simultaneously properly representing the surface warming and the 0-700m and 700-2000m heat storage is consistent with the effects of processes such as subduction and mode water formation and deepwater formation on ocean heat storage, our finding may be interpreted as indicative of the importance of direct sequestration of heat in the vertical interior of the world oceans by such processes. In particular, and as discussed in the introduction, since the (deep) ocean loses heat to the atmosphere at high latitudes through convective instabilities, when the stability of the water column is incrementally enhanced in a warming scenario (that is, the water column continues to remain unstable, but its degree of instability is slightly reduced), the heat loss from the deep ocean is slightly reduced leading effectively to a warming of the deep. At the next order, the resulting slowdown of the overturning circulation itself would further affect the heat uptake. Indeed, to examine such effects, we note that in ongoing work we consider extensions to a class of models that are used to study the overturning circulation (Tziperman, 1986; Gnanadesikan, 1999; Goodwin, 2012; Marshall and Zanna, 2014) to simultaneously represent surface warming and deeper ocean warming.

In both the four-layer and continuous variants of the EBM with non-local interactions, we also find that being able to reasonably represent the surface response and the 0-2000 m OHC is often associated with heat storage at depths below 2000 m that is greater than 10%. Clearly this finding needs to be investigated further in a hierarchical framework—a framework that will allow us to infer model parameters using the multi-model CMIP5 ensemble. In ongoing work, we are also examining (a) the implications of this study on estimates of uncertainty of climate sensitivity and (b) if these modified SCMs will help in better analyzing and inter-comparing ESMs.

Acknowledgements

BTN would like to thank W. Riley Casper for help in identifying the zero entries of $\widetilde{A}$ in Eq. 10. This research was supported by the U.S. Department of Energy (DOE) Office of Science (Biological and Environmental Research), Early Career Research program. All of the data used in this article has been previously archived and may be obtained as follows: The Levitus ocean heat content data may be obtained from the the National Oceanographic Data Center at https://www.nodc.noaa.gov. CMIP5 data may be obtained from one of the Earth System Grid Federation nodes, e.g., https://esgf-node.jpl.nasa.gov. The SAT data may be obtained from https://crudata.uea.ac.uk/cru/data/temp-erature The historical radiative forcing and the RCP8.5 scenario forcing may be obtained from (Meinshausen et al, 2011b). We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table XX of this paper) for producing and making available their model output. For CMIP the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals

Supplementary Material

∎

Improved Representation of Ocean Heat Content in Energy Balance Models. Supplementary Material

B.T. Nadiga and N.M. Urban

Appendix A Further Methodological Details

EBMs have been used extensively to obtain estimates of Equilibrium Climate Sensitivity (ECS) of AOGCMs and ESMs. However, unlike the comprehensive AOGCMs and ESMs, the simpler of the EBMs such as the ones considered in this article do not represent natural climate variability. For this reason, the Bayesian calibration procedure requires the specification of a structure for natural variability. For example, for the temperature of the upper layer, $T_{u}$,

[TABLE]

where $\epsilon_{t}$ is the discrepancy and similarly for other variables such as radiative imbalance and ocean heat content. Depending on the variable, the discrepancy structure can be reasonably approximated by simple stochastic processes such as independent and identically distributed (iid), auto-regressive processes (AR; e.g., AR(1)) or Gaussian Processes (GP). As mentioned in the main body, we used AR(1) for each of the variables considered in the Bayesian calibration procedure.

There is, however, a dependency of parameter estimates and their uncertainties on the assumed temporal structure of the discrepancy. Likewise, how the posterior distribution of calibrated responses match the original data (observational estimates or ESM responses) depends on the discrepancy model used. See Nadiga and Urban (2016) for an example of how ECS estimates and their uncertainty depends on the structure assumed for the discrepancy in the context of the two-layer EBM considered here.

Even after the structure of the discrepancy model is assumed, the results of the Bayesian inference procedure depend on whether the parameters of the discrepancy model are specified or inferred. A number of the experiments presented here were performed both ways: once while inferring, also, the parameters of the AR(1) model and a second time while specifying the parameters of the AR(1) model after first fitting that particular response alone using least-squares (i.e., after obtaining a point estimate of only the EBM parameters). In most of the cases, we found that attempting to infer the parameters of the AR(1) model simultaneously led to issues of identifiability. That is, for a significant number of the parameters related to the discrepancy model that we were, also, trying to infer, the calibration data were unable to move the prior distribution of that parameter. That is, for these parameters, the posterior distribution remained very close to the specified prior distribution. Consequently, the width of the assumed uninformative priors for these parameters led to further augmenting uncertainties in the EBM parameters. For this reason, it was deemed more appropriate to perform the inferences while specifying the parameters for the discrepancy model.

The covariance structure for the AR(1) process can be written as

[TABLE]

The specified values of $\sigma$ and $\rho$ were as follows abrupt4x SAT: (0.15 K, 0.5); abrupt4x RN: (0.27 W/m2, 0.1); 0-700 m OHC: (2 $\times 10^{22}$ J, 0.8); 700-2000 m OHC: (1 $\times 10^{22}$ J, 0.8). This is consistent with the anamolous forcing varying fastest, the SAT response which is an integrated response to the anamolous forcing varying on a slower time scale related to the lower heat capacity of the upper layer and the OHC varying slowest because of its higher heat capacity. There may be legitimate reasons to consider other values for these parameters. Further experimentation was conducted along these lines. Results of such experimentation suggest that there may be a tradeoff between fitting the various pieces of calibration data but that such tradeoff will not affect the overall nature of our results and conclusions.

Appendix B Further Two-Layer EBM Details

As mentioned in Sec. 2 of the article and seen in the top row of Fig. 2 of the article, while there are differences in representation of OHC in the two-layer EBM in experiments 1 and 2, these differences are not dramatic. For completeness, the SAT response is shown in the left panel of Fig. 5 and the differences in the SAT responses are seen to be minor as well.

However, constraining the total heat storage in the climate system in experiment 2 by calibrating against the radiative imbalance/nonequilibrium (RI/RN) response of the ESM results in significant changes to the inferred (posterior) distributions of the SCM parameters. As seen in Fig. 6, using the RN constraint leads to the inferred EBM parameters being significantly different. It is also seen in that figure that the posterior distributions themselves are much narrower.

Furthermore, correlations that exist between model parameters in experiment 1 are seen to be greatly reduced in experiment 2 (see Fig. 7).

For completeness, SAT responses in experiment 3 is shown in the right panel of Fig. 1 and the parameter distributions, and their correlations in experiment 3 are shown in Fig. 8.

Appendix C Experiments with Continuous AD-EBM Variants

Table. 3 lists some of the continuous variants of AD-EBM discussed here. We note that the calibration procedure in each of these experiments were identical to that in Experiment 3 (see Table 1). In the first version (Experiment 5a), we consider a variable thickness top layer and discretize the sub-surface thermocline region down to 2500 m into 39 equi-depth layers. Figure 9 shows that the continuously-stratified AD-EBM suffers from a similar shortcoming in not being able to simultaneously properly represent surface warming and the 0-700m and 700-2000m heat storage.

In Experiment 5a, the thickness of the top layer in the continuously stratified AD-EBM was seen to be very shallow ($\sim$ 10 m). For this reason, we conduct Experiment 5b in which the thickness of the top layer is fixed at 70 m. The SAT and OHC behavior is shown in Fig. 10. Fixing the top layer depth at a reasonable value seem to have little effect on the representation of heat storage, and the results are similar to that in Experiment 5a.

In Experiment 5c, the vertical diffusivity at the bottom of the top layer is allowed to be different from its thermocline value in order to see if this will permit a better representation of ocean heat content. Significant improvements are indeed seen in the responses of the SCM and qualitatively, these responses are similar to those in the two-layer model. Indeed, the results of this experiment seem to strongly suggest that the continuously-stratified AD-EBM suffers from a similar shortcoming in not being able to simultaneously properly represent surface warming and the 0-700m and 700-2000m heat storage.

Next, we consider vertical diffusivity to have realistic values of 10*-5* m2/s in the thermocline and 10*-4* m2/s in the abyssal ocean with a vertical profile following Bryan and Lewis, 1979. However, the transition from the thermocline value to the abyssal value occurs (and is conventionally thought to be so) at 2500 m. As discussed in the main body, observational estimates of present day heating below 2000 m is small and so, in Experiment 5d, the value of vertical diffusivity is fixed at the thermocline value of 10*-5* m2/s. Experiment 5d was among the worst in being able to match the ESM responses and the observed estimates of OHC.

Experiment 5e differs from Experiment 5c in allowing for a depth-independent upwelling velocity that is inferred in the Bayesian calibration procedure. The degree of mismatch between the SCM responses and the ESM responses and observational estimates of OHC are similar to that in Experiment 5c.

In experiment 5f, a variant of 5c, we allow the diffusivity to take on a different value below 650 m. This was motivated more by the nature of the observational estimate of OHC rather than physics. However, only minor differences are seen with respect to experiment 5c.

In Experiment 5g, we increase the number of layers in the thermocline to 60, and Fig. 15, demonstrates convergent behavior of the inability to simultaneously properly represent surface warming and the 0-700m and 700-2000m heat storage in the continuously stratified AD- EBM.

The inability of the continuously stratified AD-EBM may be understood as follows: Experiment 5c (equivalently 5g) fits ESM responses and observational estimates of OHC best. The inferred value of thermocline diffusivity in these cases range between 0.7 and 0.9 cm2/s. These values are much higher than the accepted levels of 0.1 cm2/s. One reason for these large values of vertical diffusivity is that the model is attempting to get sufficient heat down to the 700-2000 m depths in order to satisfy the observational constraint. And in trying to do so, it is clear from these experiments that this is causing too large a heating in the 0-700 m depths—a feature that is consistent across a number of the experiments described here. That is, assuming that vertical diffusivity is not varying with time (and there is no reason to expect it to have varied substantially with time in the post industrial period) getting observed levels of warming into the 700-2000 m layer while limiting the warming of the 0-700m layer to observed levels is difficult to achieve in the AD-EBM.

Appendix D Further Experiments with Non Local Interactions

We next consider a modification of the continuous version of the AD-EBM that addresses the deficiency. In the original approach of Hoffert, 1980 to the AD-EBM, the direct sequestration of heat by the polar downwelling branch in the cartoon of the overturning circulation in the world oceans is not considered; only the effect of the resultant weak and broad upwelling is considered. We hypothesize that the deficiency of the AD-EBM identified above can be remedied by introducing a parameterization of the direct sequestration of heat by processes related to polar downwelling in the AD-EBM. To test this hypothesis, we consider a parameterization of the direct sequestration of heat by polar downwelling, and in anology with the formation of the North Atlantic Deep Waters. The structure of a nominal overturning cell in which the NADW is a component has a maximum in the overturning streamfunction around a km in depth. Consequently, we consider a profile of exchange of heat between the first sub-surface layer and the interior that has a Gaussian profile in the vertical centered at $z_{mx}$ and a vertical scale of $\Delta z$, and where both $z_{mx}$ and $\Delta z$ are inferred. That is, the anomaly-diffusion nature of the thermocline is augmented as

[TABLE]

Here $\gamma_{d}$ is a scalar coefficient of heat exchange associated with the parameterization and is inferred in the calibration process. Figure 16 shows that the SAT and OHC responses are fit reasonably in the model. Furthermore, in this simplified parameterization, if half of the depth $z_{mx}$ is taken as the depth at which the maximum of the overturning streamfunction occurs, the maximum overturning depth is inferred to lie at about a km, which is somewhat realistic. Finally, in Fig. 17, it is seen that the changes are minimal when the number of vertical layers is increased from 40 to 60.