On the Effects of Mixed and Deep Ocean Layers on Climate Change and Variability

Abstract

The ocean, one of the five major components of the Earth’s climate system, plays a key role in climate-forming processes, affecting its change and variability. The ocean influences climate over a wide range of time–space scales. To explore the climate, its components, interactions between them and, in particular, the effect of the ocean on weather and climate, researchers commonly use extremely complex mathematical models of the climate system that describe the atmospheric and ocean general circulations. However, this class of climate models requires enormous human and computing resources to simulate the climate system itself and to analyze the output results. For simple climate models, such as energy balance and similar models, the computational cost is insignificant, which is why these models represent a test tool to mimic a complex climate system and obtaining preliminary estimates of the influence of various internal and external factors on climate, its change and variability. The global mean surface temperature (GMST) and its fluctuations in time serve as critical indicators of changes in the climate system state. We apply a simple two-box ocean model to explore the effect of mixed and deep ocean layers on climate-forming processes and especially on climate change and variability. The effect of mixed and deep ocean layers on GMST is parameterized via the layers’ effective heat capacities and heat exchange between layers. For the listed parameters, the sensitivity functions were derived numerically and analytically, allowing one to obtain an idea of how the mixed and deep ocean layers affect climate change and variability. To study climate change, a deterministic version of the model was used with radiative forcing parameterized by both stepwise and linear functions. In climate variability experiments, a stochastic version of the model was applied in which the radiative forcing is considered as a delta-correlated random process.

1. Introduction

Exploring contemporary climate change and assessing its impact on people, communities, infrastructure, economic activity and natural systems represents one of the most important challenges of our time [1,2]. This extremely important and topical issue is in the focus of ongoing discussions not only among the professional communities, but also among politicians, public figures and the majority of ordinary people. The problem of climate formation and change under the influence of natural and anthropogenic factors is undoubtedly interdisciplinary and multifaceted, since the Earth’s climate system (ECS) is a large-scale hierarchical physical system of extreme complexity, and the main tools for its study are mathematical climate models numerically implemented on supercomputers (e.g., atmosphere-ocean general circulation models or AOGCM’s [3,4,5,6]).

In climate theory, a whole range of problems is associated with the assessment of the functional influence of the ocean on climate formation, its change and variability (e.g., [7,8,9,10]). The ocean, as we know, is one of the five major components of the ECS influencing climate over a wide range of time–space scales [11]. The ocean, having a huge heat capacity and being a store of the vast majority of excess heat, largely provides the property of thermal inertia in the climate system. Due to this property, the climate system acquires resilience or slowness to changes in significant (external) factors affecting the climate, such as, for example, the content of greenhouse gases (GHGs) in the atmosphere. In the context of the climate change problem, this means that the ocean’s thermal inertia delays global warming to some extent, and even if GHG emissions remain at current levels, the surface temperature on our planet will continue to rise (at least for decades or even centuries) as long as the climate system reaches a new equilibrium state [12,13,14,15]. In climate change and variability studies, the overwhelming majority of coupled atmospheric–ocean models take into account only the upper mixed layer of the ocean (with different depth) since this layer modulates the atmosphere–ocean interaction [16,17]. However, the importance of the deep ocean layers for the global climate system should not be underestimated as well [18,19], since the deep ocean holds excess heat, therefore affecting long-term climate formation processes. Generally, the atmosphere and ocean represent a large-scale interconnected physical system in which the ocean is responsible for the formation of the low-frequency part of the spectrum of variability and significant time lag between “action and effect”.

This article is aimed at studying the effect of mixed and deep ocean layers on climate-forming processes, climate change and variability. To achieve this objective, we use the conceptual two-box ocean model that describes the influence of mixed and deep ocean layers of the global mean surface temperature (GMST) [20,21,22]. This model is essentially a two-layer energy balance model (EBM) that does not consider the atmosphere. This assumption is physically justified, since, as can be easily shown [23], the effective heat capacity of the upper quasi-homogeneous (mixed) layer of the ocean is more than 20 times greater than the effective heat capacity of the entire atmosphere. In passing, we note that the original EBM, introduced by Budyko [24] and Sellers [25], belongs to the class of the simplest climate models. Although complex coupled atmosphere–ocean–land–ice models that primarily describe the atmospheric and ocean general circulation are the main tools used in climate research, simplified models such as EBM remain popular among researchers because they make it possible to quickly and efficiently obtain preliminary estimates of the influence of the main climate-forming factors on the essential features of climate change and variability. It has a certain logic behind it, since the response of extremely complex climate models to external forcing can be described with a few simple parameters, similar to those of a simple model. Moreover, complex climate models require enormous human and computing resources to simulate the climate system itself and analyze the outputs. Meanwhile, for simple climate models, such as being used in this study, the computational cost is insignificant, which is why these models can be successfully used as test tools to mimic a complex climate system. In this regard, a recently published paper [26] in which the performance of a two-box EBM (which is similar to the one used in our study) in simulating historical and future surface temperature projections from Coupled Model Intercomparison Project Phase 6 (CMIP6) models was analyzed in detail. To assess the EBM performance, the authors of this paper used the CMIP6 historical and future shared socio-economic pathway (SSP) projections for AOGSMs [27]. Results obtained showed that the EBM prediction errors for future global surface temperature projections differ significantly between AOGSM, radiative forcings produced by both greenhouse gases and aerosols, time periods and methods of EBM calibration. According to the authors of [26], the model performance can be improved by introducing an efficacy factor into the model equations, as well as by incorporating time variations in climate feedbacks. Nevertheless, the authors of the article concluded that the two-layer EBM allows for obtaining quite realistic results for pre-defined climate change scenarios and periods.

The two-box ocean model considered here predicts the evolution of two dependent variables: the GMST anomaly (deviation form a certain norm) that characterizes the state of the ocean mixed layer, and the global mean deep ocean temperature (GMDOT) anomaly that describes the deep ocean state. The change in GMST and its fluctuations serve as a key indicator of climate change and variability. For this reason, we will explore the sensitivity of GMST and its fluctuations with respect to the main parameters characterizing the influence of the ocean on the Earth’s climate within the framework of a two-box model.

This paper is organized as follows. In Section 2, we present the essential features of the two-box ocean model used in calculations, its analytic solution and sensitivity study approach. Results obtained are considered in Section 3. Discussions and concluding remarks are given in Section 4.

2. Materials and Methods

In this paper, we apply mathematical modeling and dynamical systems sensitivity analyses to study the effects of mixed and deep ocean layers on GMST and its fluctuations. All calculations were performed using the two-box ocean model. This model makes it possible to explore the dependence of climate system response to external radiative forcing on changes in the model parameters characterizing the thermal inertia of both the mixed and deep ocean layers. To quantify the influence of variations in the model parameters on the GMST fluctuations, sensitivity functions derived analytically are used.

2.1. The Model

The mathematical model used consists of two boxes. The first one (upper box) corresponds to the mixed ocean layer, and the second one (lower box) represents the deep ocean. We consider the intermediate layer waters that are several hundred meters deep and that affect the climate on time scales of years to decades. However, to be consistent with [20,21,22], we use the term “deep ocean” for these layers. It is usually presumed that only the upper few hundred meters of the ocean are involved in exchange of heat with the atmosphere on decadal time scales. In this case, the solutions of model equations reach an equilibrium state after several hundred years of integration. On longer time scales, the role of the deeper layers of the ocean increases.

The state of the model boxes is characterized by the temperature anomalies $T_a$ (upper box) and $T_o$ (lower box). In doing so, the change in $T_a$ is identified with change in the GMST. The model equations in stochastic formulation can be written as follows [20,21,22,28,29]:

$$C_a\left(d{T}_a/dt\right)=-\lambda T_a-\gamma \left(T_a-T_o\right)+F_d+F_r,$$

(1)

$$C_o\left(d{T}_o/dt\right)=\gamma \left(T_a-T_o\right),$$

(2)

where $C_a$ and $C_o$ (both in ${\mathrm{W}\mathrm{yr}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$) are effective heat capacities for the upper and lower boxes (note that $C_o\gg C_a$); $\lambda \left({\mathrm{W}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}\right)$ is the climate feedback parameter characterizing the change in downward top-of-the-atmosphere (TOA) radiative flux for a given change in surface temperature; $\gamma$ (${\mathrm{W}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$) is a heat exchange coefficient describing the deep ocean heat uptake; $F_d$ and $F_r$ are, respectively, top of the atmosphere deterministic and random radiative forcing terms (${\mathrm{W}\mathrm{m}}^{-2}$). Note that the EBM becomes completely deterministic if we set the random term $F_r$ to zero.

Using a specially designed fitting procedure, Geoffroy et al. [22] adjusted the two-box model parameters so that the GMST change calculated by this model matched the GMST change obtained in CMIP5 (Coupled Model Intercomparison Project, Phase 5) models [16] for step and linear forcing. The values of model parameters that are consistent with the CMIP5 multi-model mean are as follows: $C_a=7.3$ ${\mathrm{W}\mathrm{yr}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$, $C_o=105.5 {\mathrm{W}\mathrm{yr}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$, $\lambda =1.13$ ${\mathrm{W}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$, $\gamma =0.7 {\mathrm{W}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$. It can be easily shown that the mentioned values of ocean heat capacities $C_a$ and $C_o$ are equivalent, respectively, to a 75-meter-thick mixed layer and a deep-ocean layer depth equal to about 1100 m [22].

The model (1)–(2) is essentially a forced-damping oscillator, and its response to an external forcing is determined by fast ${\tau }_f$ and slow ${\tau }_s$ relaxation times defined by the following expressions, respectively [30]:

$${\tau }_f=\frac{\beta \left(1-\sqrt{1-{\omega }_0^2/{\beta }^2}\right)}{{\omega }_0^2}, {\tau }_s=\frac{\beta \left(1+\sqrt{1-{\omega }_0^2/{\beta }^2}\right)}{{\omega }_0^2},$$

(3)

where $\beta =\left[\left(\lambda +\gamma \right)C_D+\gamma C\right]/\left(2C{C}_D\right)$ is the damping constant, and ${\omega }_0=\sqrt{\lambda \gamma /C{C}_D}$ is the (angular) natural frequency of free oscillations. For the used values of model parameters (see above), we have ${\tau }_f\approx 3.9 \mathrm{yr}$ and ${\tau }_s\approx 240 \mathrm{yr}$.

In the general case, the response of climate system, represented by linear time invariant model Equations (1) and (2), to an arbitrary but sufficiently small time dependent deterministic radiative forcing $F_d\left(t\right)$, can be estimated via the so-called impulse response function $h\left(t\right)$ (IRF) as a convolution of two functions, $h\left(t\right)$ and $F_d\left(t\right)$:

$$T_a\left(t\right)=\left[h\ast F_d\right]\left(t\right)=\underset{0}{\overset{t}{{\displaystyle \int }}}h\left(\tau \right)F_d\left(t-\tau \right)d\tau ,$$

(4)

where the IRF corresponding to the two-box ocean model is given by [30,31].

$$h\left(t\right)=\frac{{\tau }_f{\tau }_f}{C\left({\tau }_f-{\tau }_f\right)}\left[\left(\frac{1}{{\tau }_f}-\frac{\gamma }{C_D}\right)e^{-t/{\tau }_f}-\left(\frac{1}{{\tau }_s}-\frac{\gamma }{C_D}\right)e^{-t/{\tau }_s}\right].$$

(5)

However, for forcing represented by a step or liner function, the solution for the GMST anomaly $T_a$ as a function of time can be found analytically. For forcing specified by a step function:

$$F_d\left(t\right)=\{\begin{array}{c}0 \mathrm{at} t<0 \\ F_A \mathrm{at} t \ge 0\end{array},$$

(6)

the corresponding solution is as follows [31]:

$$T_a\left(t\right)=\frac{F_A}{\lambda }\left[1-{\alpha }_1{e}^{-t/{\tau }_f}-{\alpha }_2{e}^{-t/{\tau }_s}\right],$$

(7)

where $F_A$ is a given constant.

For a linear forcing $F\left(t\right)=\eta t$, where $\eta$ is a parameter that determines the rate of forcing growth, the solution takes the form [31]:

$$T_a\left(t\right)={\alpha }_3{e}^{-t/{\tau }_f}+{\alpha }_4{e}^{-t/{\tau }_{fs}}+\left(\eta /\lambda \right)t-\eta \left(C_a+C_o\right)/{\lambda }^2.$$

(8)

In Expressions (7) and (8), the coefficients ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$ and ${\alpha }_4$ depend on the model parameters [31].

Deterministic version of the two-box ocean model ($F_d\ne 0$, $F_r=0$) is used to study the influence of non-random radiative forcing on a climate system, the state of which is characterized by the GMST anomaly $T_a$. Meanwhile, the stochastic version of the model ($F_d=0$, $F_r\ne 0$) is a tool for assessing the impact of random radiative forcing on climate variability, characterized by the variance of the GMST anomaly ${\sigma }_T^2$.

In climate research [32], random radiative forcing is usually considered additive and is parameterized by Gaussian $\delta -$ correlated in time process with zero mean $F_s\left(t\right)=0$ and correlation function given by $R_r\left(t_2-t_1\right)=\langle F_r\left(t_1\right)F_r\left(t_2\right)\rangle =2D_r\delta \left(t_2-t_1\right)$, where $\delta$ is the Dirac delta function and $D_r$ is the diffusion coefficient defined by the variance of Gaussian random process ${\sigma }_r^2$ and its correlation time ${\tau }_r$: $D_r={\sigma }_r^2{\tau }_r$ [33,34].

The variance of the GMST anomaly can be found via the power spectral density (PSD) of GMST fluctuations $S_{TT}\left(\omega \right)$, where $\omega$ is an angular frequency of them [28,29]. In turn, the PSD $S_{TT}\left(\omega \right)$ is expressed through the transfer function $H\left(\omega \right)$ of the two-box model and the power spectrum $S_{SS}\left(\omega \right)$ of the input signal as follows: $S_{TT}\left(\omega \right)={\left|H\left(\omega \right)\right|}^2{S}_{SS}\left(\omega \right)$. The model used represents the so-called linear time-invariant (LTI) dynamical system, since the latter is described by linear ordinary differential equations with constant coefficients [35]. LTI systems can be examined in the frequency domain using the system’s transfer function $H\left(s\right)$, which is the Laplace transform of its IRF $h\left(t\right)$ with zero initial conditions: $H\left(s\right)=ℒ\left\{h\left(t\right)\right\}$, where $ℒ$ denotes the Laplace transform operator and $s=\sigma +i\omega$ is the Laplace variable also known as a complex frequency. In our case, the input signal is a random Gaussian process for which the PSD is the Fourier transform of its autocorrelation function [36]: $S_{SS}\left(\omega \right)={{\displaystyle \int }}_{-\infty }^{\infty }{e}^{-i\omega \tau }{R}_r\left(\tau \right)d\tau ={\sigma }_r^2{\tau }_r/\pi$. It can be shown that the transfer function of the two-box ocean model $H\left(\omega \right)$ is as follows [29].

$$H\left(\omega \right)=\frac{\gamma +i\omega C_o}{C{C}_o\left[\left({\omega }_0^2-{\omega }^2\right)+i2\omega \beta \right]}.$$

(9)

Thus, we have the following formula for PSD of GMST fluctuations:

$$S_{TT}\left(\omega \right)={\left|H\left(\omega \right)\right|}^2{S}_{SS}\left(\omega \right)=\frac{q_s^2}{\pi C{C}_o}\frac{{\gamma }^2+{\omega }^2{C}_o^2}{{\left({\omega }_0^2-{\omega }^2\right)}^2+4{\omega }^2{\beta }^2},$$

(10)

where $q_r^2=2D_r=2{\sigma }_r^2{\tau }_r$.

Integrating (10) over all positive frequencies, we obtain an expression for the variance of GMST fluctuations [37]:

$${\sigma }_T^2=\frac{q_r^2}{2\lambda C}\frac{\gamma C+\lambda C_o}{\gamma C+\left(\gamma +\lambda \right)C_o}.$$

(11)

2.2. Sensitivity Analysis

To study the effect of mixed and deep ocean layers on the GMST and its fluctuations, we use sensitivity functions that are partial derivatives of the GMST and its variance with respect to the parameters $C_a$, $C_o$ and $\gamma$ [38]. Sensitivity functions characterizing the effect of the above-mentioned parameters on the GMST are defined as follows:

$$S_p={\partial T_a/\partial p|}_{p=p_0}, p=\left(C_a, C_0, \gamma \right),$$

(12)

where $p_0$ is the base parameter value around which the sensitivity is estimated.

Using sensitivity functions, one can easily evaluate the effect of model parameter variations on the GMST anomaly change: $\mathsf{\Delta }{T}_a=\delta p·{S_p|}_{p=p_0}$, where $\delta p$ is the variation in the parameter $p$ (sufficiently small departure of the parameter from its base value $p_0$). For simplicity, the sensitivity functions (12) are estimated numerically.

In contrast, the sensitivity functions that show the effect of model parameter variations on GMST fluctuations (i.e., on climate variability) are estimated analytically by differentiating the expression for the variance of GMST fluctuations (11) with respect the corresponding parameter:

$$S_{PSD,p}={\partial S_{TT}/\partial p|}_{p=p_0}, p=\left(C_a, C_0, \gamma \right),$$

(13)

For example, expression for sensitivity function $S_{PSD,C}$ is obtained by successively differentiating Equation (11) with respect to the parameter $C$:

$$S_{PSD,C}=\frac{1}{\pi }\frac{2q_r^2}{{C_a}^2{C}_D^2}\frac{{\gamma }^2+{\omega }^2{C}_D^2}{{\left[{\left({\omega }_0^2-{\omega }^2\right)}^2+4{\beta }^2{\omega }^2\right]}^2}\frac{{\omega }^2}{{C_a}^2}\left[2\beta \left(2\beta C_a-\lambda -\gamma \right)-\left({\omega }_0^2-{\omega }^2\right)C_a\right].$$

(14)

Similarly, one can derive expressions for sensitivity functions $S_{PSD,C_o}$ and $S_{PSD,\gamma }$. Differentiation of Expression (12) with respect to parameters $C_o$ and $\gamma$ leads to too bulky formulas, which are not presented here for the sake of brevity.

3. Results

It is assumed that in the two-box ocean model, climate change is associated with a change in GMST $T_a$ due to deterministic radiative forcing, while climate variability, in turn, is identified with the variance of GMST fluctuations ${\sigma }_T^2$ generated by a random radiative perturbation. The two-box model, both deterministic and stochastic, predicts the evolution of two dependent variables, $T_a$ and $T_o$, that characterize the state of the ocean mixed layer and the deep ocean state, respectively. The evolution of these variables depends on the model parameters, including those describing the climate inertia of the upper mixed ocean layer (parameter $C_a$), the climate inertia of the deep ocean layers (parameter $C_0$), and the heat exchange between these layers (parameter $\gamma$). Since, as already noted, the behavior of GMST serves as a key indicator of climate change and variability, we explore the effect of model parameters characterizing the influence of the ocean on the Earth’s climate in the two-box model, on GMST and its fluctuations using sensitivity functions.

As mentioned earlier, the model parameters used in calculations were obtained in [22] via a specifically designed fitting approach for step and linear radiative perturbations, namely for an instantaneous quadrupling of carbon dioxide (${\mathrm{CO}}_2$) content ($F_{4\times {\mathrm{CO}}_2}\approx 7.4 {\mathrm{W}\mathrm{m}}^2$ [39]) and for linear forcing, In order to be consistent with [22], we first consider the results obtained for abrupt radiative experiments, assuming $F=F_{4\times {\mathrm{CO}}_2}$ if $t\ge 0$, and second for logarithm relationship between radiative perturbation and the ${\mathrm{CO}}_2$ content [39], which gives the following approximate relation: $F=\eta t$, where $\eta =5.92·{10}^{-2} {\mathrm{W}\mathrm{m}}^{-2} {\mathrm{yr}}^{-1}$ [28].

3.1. Step Forcing (Deterministic Model)

The evolution of temperature anomalies of the upper and lower boxes, calculated in 1000-year intervals, is shown in Figure 1. Both $T_a$ and $T_0$ increase exponentially toward the equilibrium $T_{eq}$ determined by the radiative forcing $F_{4\times {\mathrm{CO}}_2}$ and climate feedback $\lambda$: $T_{eq}=F_{4\times {\mathrm{CO}}_2}/\lambda$. In the previous Section, we draw attention to the fact that the model response to an external radiative perturbation is characterized by fast ${\tau }_f$ and slow ${\tau }_s$ relaxation times defined by (3). An increase in temperature anomaly of the upper box (mixed ocean layer) is determined by the fast relaxation time (positive effect) and the deep ocean heat uptake (negative effect), while the increase in temperature anomaly of the lower box (deep ocean layers) is affected by the slow relaxation time and the deep ocean heat uptake, both effects having a positive sign.

Sensitivity functions of GMST and GMDOT anomalies to parameter $C_a$, denoted by $S_{C_a,T_a}$ and $S_{C_a,T_o}$, respectively, are illustrated in Figure 2. These functions provide an idea of the nature of upper mixed ocean layer effects on temporal change in $T_a$ and, for completeness, the change in $T_0$. As Figure 2 shows, the absolute value of $S_{C_a,T_a}$ reaches a maximum absolute value of $0.23 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}{\mathrm{yr}}^{-1}$ at $t={\tau }_f$, while the sensitivity function of $T_o$ with respect to parameter $C_a$ has a maximum absolute value of $0.014 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}{\mathrm{yr}}^{-1}$ at $t\approx 26 \mathrm{yr}$. The sensitivity of both $T_a$ and $T_o$ with respect to parameter $C_a$ weakens over time after they reach their maximum values.

Using sensitivity functions, we can estimate the influence of parametric uncertainty on changes in the GMST anomaly. Let the uncertainty in the parameter $C_a$ be 1% of its base value, i.e., $\delta C_a=0.073$ ${\mathrm{W}\mathrm{yr}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$. This uncertainty produces the change in the GMST anomaly of $\mathsf{\Delta }{T}_a=\pm 0.073·0.23=\pm 0.017 \mathrm{K}$ (this estimation is obtained for the maximum absolute value of $S_{C_a,T_a}$). The effect of variation in the parameter $C_a$ on the GMDOT anomaly is much less. For the maximum absolute value of the sensitivity function $S_{C_a,T_o}$ (see above), the uncertainty in GMDOT generated by 1% uncertainty in the parameter $C_a$ is $\mathsf{\Delta }{T}_o=\pm 0.073·0.014=\pm 1.02·{10}^{-3} \mathrm{K}$, which is 16 times less than $\mathsf{\Delta }{T}_a$. Thus, the thermal inertia of the upper ocean mixed layer, characterized by the effective heat capacity $C_a$, significantly affects the change in GMST anomaly, while its effect on the deep ocean temperature anomaly $T_o$ is insignificant, which is physically obvious.

A somewhat different picture emerges when estimating the sensitivity of GMST and GMDOT anomalies with respect to the thermal inertia of the deep ocean layers $C_o$. Figure 3 illustrates the sensitivity functions of GMST and GMDOT anomalies, denoted by $S_{C_o,T_a}$ and $S_{C_o,T_o}$, respectively, to parameter $C_o$. Both functions are convex downward and reach their absolute maximum values of $8.84·{10}^{-3}$ and $0.024 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}{\mathrm{yr}}^{-1}$ at $t={\tau }_f$. For the 1% of uncertainty in the parameter $C_o$ ($\delta C_D=1.055 {\mathrm{W}\mathrm{yr}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$), the changes in GMST and GMDOT anomalies are, respectively, as follows: $\mathsf{\Delta }{T}_a=\pm 9.34·{10}^{-3} \mathrm{K}$ and $\mathsf{\Delta }{T}_o=\pm 23.77 ·{10}^{-3} \mathrm{K}$. Thus, the value of $\mathsf{\Delta }{T}_a$ is quite comparable with the value of $\mathsf{\Delta }{T}_o$. After the functions $S_{C_o,T_a}$ and $S_{C_o,T_o}$ have reached their maximum modulo values, they decrease toward zero at $t\to \infty$.

Deep ocean heat uptake, or, in other words, the heat exchange between the mixed and deep ocean layers, also affects GMST change. As we mentioned earlier, deep ocean heat uptake contributes to a certain decrease in GMST and, at the same time, some increase in GMDOT. Sensitivity functions of $T_a$ and $T_o$ with respect to the parameter $\gamma$ characterizing the deep ocean heat uptake, presented in Figure 4, to a certain extent confirm this conclusion. Thus, for example, the sensitivity function $S_{\gamma , T_a}$ reaches an extremum of $-1.84 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}$ at $t\approx {\tau }_f$ and retains a negative sign on the time interval $t\in \left[0, {\tau }_s \right]$. At $t\approx {\tau }_s$, this sensitivity function transits into the domain of positive numbers, but at the same time, its values remain relatively small. In turn, the sensitivity function $S_{\gamma , T_o}$, being positive on the considered time interval $t\in \left[0, 1000 \right] \mathrm{yr}$, increases over time, reaching a maximum of $1.97 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}$ at $t\approx {\tau }_s$, and then smoothly decreases to the value of $0.34 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}$. at $t=1000$ yr.

3.2. Linear Forcing (Deterministic Model)

Since the climate system trajectory is largely determined by external radiative forcing, we also considered a radiative perturbation described by a linear function of time that corresponds to a 1% increase in the atmospheric ${\mathrm{CO}}_2$ content. The graphs of the functions $T_a\left(t\right)$ and $T_0\left(t\right)$ calculated for this case are shown in Figure 5. As expected, a linear increase in radiative forcing leads to a continuous increase in both $T_a$ and $T_0$.

Sensitivity function $S_{C_a,T_a}$, characterizing the effect of the upper mixed ocean layer on GMST, shows a rather rapid decrease during about the first 15 years ($t\approx 4{\tau }_f$), which continues further, but at a slower rate (see Figure 6). At the same time, the sensitivity function $S_{C_a,T_o}$ describing the influence of the upper mixed ocean layer on the deep ocean temperature decreases monotonically with time from the beginning. After a long time ($t\to \infty$), the two curves, $S_{C_a,T_a}$ and $S_{C_a,T_o}$, that are functions of time, de facto merge, reaching a value of $-4.1\times {10}^{-2} {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}{\mathrm{yr}}^{-1}$.

Figure 7 displays the graphs of sensitivity functions $S_{C_o, T_a}$ and $S_{C_o, T_o}$. The shape of the curves of these functions undoubtedly differs from similar graphs of functions $S_{C_o, T_a}$ and $S_{C_o, T_o}$ constructed for step-forcing and are shown in Figure 3. However, regardless of the type of function describing the radiative forcing, the effect of the upper mixed ocean layer is more pronounced on the temperature regime of the deep ocean than on the near-surface temperature. On the time interval $tϵ\left[0, 1000\right] \mathrm{yr}$, both functions gradually decrease, reaching the stationary mode.

Finally, let us consider the behavior of sensitivity functions $S_{\gamma , T_a}$ and $S_{\gamma , T_o}$ calculated for radiative forcing parameterized by a linear function of time (see Figure 8). In contrast to the sensitivity function $S_{\gamma , T_a}$ calculated for step-forcing (see Figure 4), in the case with linear forcing, $S_{\gamma , T_a}$ reaches its maximum absolute value of $1.32 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}$ at $t=t_f$. Then, the modulus function $\left|S_{\gamma , T_a}\right|$ decreases monotonically, tending to zero at $t\to \infty$. In turn, $S_{\gamma , T_o}$ is a monotonically increasing function tending to the value of $9.15 {\mathrm{K}}^2{\mathrm{m}}^2{\mathrm{W}}^{-1}$ in the limit of $t\to \infty$. Plots of functions $S_{\gamma , T_a}$ and $S_{\gamma , T_o}$ clearly show a significantly greater sensitivity of $T_0$ to the parameter $\gamma$ than sensitivity of $T_a$ (see note below).

Note. The sensitivity functions presented in each of the figures have the same dimension; we can easily determine which variable, $T_a$ or $T_0$, is more or less sensitive to one or another parameter.

3.3. Random Forcing

A stochastic version of the two-box ocean model is used to estimate the effect of ocean thermal inertia and deep ocean heat uptake on climate variability, which is identified with the variance of GMST fluctuations ${\sigma }_T^2$ induced by random radiative forcing. In turn, random radiative forcing, as we mentioned above, is parameterized by the delta-correlated in time random process and is asymptotically estimated as follows [34]: $q_r^2={\tilde{\sigma }}_r^2{\tau }_r$, where ${\tilde{\sigma }}_r^2$ is the radiative forcing variance averaged over the time interval $\left[0, {\tau }_r\right]$. Assuming ${\tau }_r=1 \mathrm{yr}$, we obtain that $q_r^2={\tilde{\sigma }}_r^2$. In this study, the standard deviation value of ${\tilde{\sigma }}_r\approx 0.26 {\mathrm{W}\mathrm{m}}^{-2}$ is used [28,29]. It is quite obvious that the parameterization of random radiative forcing by additive delta-correlated in the time process does not fully describe the entire spectrum of random perturbations both external and internal on the climate system. However, as shown in a number of previous publications (e.g., [28,29,30,40,41,42] and references therein), models similar to the one used in this study have proven to be a useful climate research tool able to reproduce climate variability on time scales from years to decades.

Sensitivity functions of power spectral density of GMST fluctuations with respect to the mixed ocean layer thermal inertia and the deep ocean thermal inertia determined, respectively, by the effective heat capacities $C_a$ and $C_o$ are presented in Figure 9 and Figure 10. In these figures, the x-axis is plotted on a logarithmic scale, where $\nu$ is a regular (linear) frequency. Sensitivity functions, both $S_{PSD, C_a}$ and $S_{PSD, C_o}$, are highly nonlinear, having pronounced extrema. For example, the function $S_{PSD, C_a}$ has a geometric shape that resembles an upside-down bell, the top of which corresponds to the temperature fluctuations with about a 41-year period $1/\nu ={\left(\lambda /2\pi C_a\right)}^{-1}$. This means that interdecadal surface temperature fluctuations are most sensitive to the thermal inertia of the upper mixed ocean layer. Therefore, when studying interannual and interdecadal climate variability, models of the atmosphere–ocean system are mainly used, in which only the upper ocean layer is described. In turn, the extremum of sensitivity function $S_{PSD, C_o}$, the geometric shape of which also resembles an inverted bell, is reached at a frequency of surface temperature fluctuations, the period of which is approximately 1150 years. Thus, when studying climate variability on intercentury time scales, coupled atmosphere–ocean models should take into consideration the climate inertia of the deep ocean layers.

Sensitivity function of power spectral density of GMST fluctuations with respect to the deep ocean heat uptake $S_{PSD, \gamma }$ determined by the heat exchange coefficient $\gamma$ is illustrated in Figure 11. This figure shows that climate variability on interannual, decadal, and interdecadal time scales is insignificant. The most sensitive are surface temperature fluctuations with about a 170-year period. The sensitivity function $S_{PSD, \gamma }$ has another extremum corresponding to surface temperature fluctuations, the period of which is about 2250 years. However, in absolute value, this extremum is about four times less than the first one.

4. Discussions and Concluding Remarks

In this paper, we have outlined the ability of sensitivity analysis for assessing the effects of upper and deep ocean layers on climate change and variability, which are identified with trends in GMST and variance of its fluctuations, respectively. As the main research tool, a conceptual climate model was used, namely a two-box ocean model in deterministic and stochastic formulation. To estimate the effects of thermal inertia of the upper and deep ocean layers and heat exchange between them on GMST trends (climate change), a deterministic version of the model was used, in which the climate system is affected by deterministic radiative forcing, approximated by stepwise and linearly increasing functions. In turn, the stochastic version of the model was applied to study the effects of mixed and deep ocean layers and deep ocean heat uptake on the temporal power spectrum of GMST fluctuations generated by random radiative forcing.

Sensitivity functions calculated via deterministic models give a clear idea of the GMST sensitivity with respect to model parameters $C_a, C_o$ and $\gamma$ in the time domain. Figure 2, Figure 3 and Figure 4, Figure 6, Figure 7 and Figure 8 are de facto time-domain graphs showing how GMST sensitivities with respect to model parameters change with time. In other words, time-domain analyses allow one to trace temporal dynamics of sensitivity functions and at any instant of time estimate the effect of variations in model parameters on GMST change. At the same time, using stochastic model one can calculate sensitivity functions of the power spectrum of GMST fluctuations with respect to model parameters. Graphs of these functions (Figure 9, Figure 10 and Figure 11) show the frequency distribution of sensitivity of GMST variance to model parameters (recall that climate variability is identified with the variance of GMST fluctuations).

The abrupt $-4\times {\mathrm{CO}}_2$ scenario used in this study corresponds to a constant radiative forcing, which means unaltered concentrations of greenhouse gases in the atmosphere. The sensitivity functions of GMST to variations in model parameters characterizing the thermal inertia of the upper and deep ocean layers were calculated for this scenario. Sensitivity functions allow for quantifying the thermal damping effect of the ocean on climate change and estimating the influence of the uncertainties in model parameters characterizing the thermal inertia of the upper and deep ocean layers on the change in GMST at constant radiative forcing. Analyzing the sensitivity functions, one can conclude that the thickness of the mixed ocean layer, which, along with the specific heat of sea water, determines the effective heat capacity of the upper ocean layer, is a critically important parameter that significantly affects the realistic simulation of climate system response to a constant external radiative forcing.

It is known that there are rather significant uncertainties in the climate change projections and long-term predictions of surface temperature trends obtained for given greenhouse gas emission scenarios from different complex climate models [43,44]. These uncertainties arise due to inter-model differences in the description (parameterization) of major climate-formation physical processes. Some of these uncertainties can be apparently explained by the differences in the parameterization of the upper layer of the ocean.

The GMST response to abrupt $-4\times {\mathrm{CO}}_2$ forcing used in this study and the corresponding sensitivity functions of GMST with respect to thermal inertia of upper and deep ocean layers and heat exchange between them can be used to estimate the climate system response to an arbitrary forcing scenario. This can be explained by the fact that an arbitrary radiative forcing can be approximated by a certain sequence of small stepwise perturbations, and the climate system response to these perturbations in the first approximation is a linear combination of individual responses [45]. As an example, we considered a linearly increasing radiative forcing [22]. The sensitivity functions of GMST to variations in the parameters $C_a$ and $C_o$ calculated for this scenario also confirm the thermal damping effect of the ocean on climate change and the importance of correct parameterization of the upper ocean layer including its thickness when projecting climate using complex climate models. The deep ocean heat uptake, characterizing by the parameter $\gamma$, is also the thermal damping mechanism that contributes to a decrease in GMST. However, the significance of this mechanism is noticeable on centennial time scales.

Thus, in a deterministic approach, sensitivity analysis of climate trends to model parameters shows how model parameters affect the GMST response on external radiative forcing over time: the greater the absolute value of the sensitivity function, the greater the response of the climate system. If on annual and decadal time scales the effect of thermal inertia of the upper ocean layer prevails over the thermal inertia of the deep ocean layers and the heat exchange between the layers, then over time, the picture changes significantly: these physical mechanisms begin to exert an ever greater influence on the surface temperature response to an external forcing, which is physically quite obvious.

As is known, random radiative forcing generates surface temperature fluctuations that contribute to climate variability. Temporal fluctuations in surface temperature are commonly analyzed by examining its power spectrum, which can be interpreted in terms of surface temperature variance at the respective frequency (period) of fluctuations. The effects of ocean thermal inertia and deep ocean heat uptake on climate variability can be seen from the change in PSD caused by variations in model parameters $C_a$, $C_o$ and $\gamma$. Analysis of sensitivity functions of the power spectrum of GMST fluctuations to variations in the listed above model parameters calculated under the influence of random forcing on the climate system confirms the well-known fact that in order to reproduce the essential features and patterns of climate variability on inter-annual, decadal and inter-decadal time scales, it is sufficient to take into consideration the thermal inertia of the mixed ocean layer in climate models (e.g., [46]), while the influence of deep ocean thermal inertia and deep ocean heat uptake on climate variability is manifested on much larger time scales. All calculated sensitivity functions of the power spectrum to model parameters are characterized by nonlinearity and significant changes in the considered frequency range. Extrema of sensitivity functions correspond to the frequency (period) of GMST fluctuations, which are the most responsive to change in the parameter under consideration. Although the two-box ocean model is an extremely simplified representation of the global climate system, the results obtained may be of some interest, at least from a qualitative point of view, for a better understanding of the effects of ocean thermal inertia and heat exchange between the upper and deep ocean layers on climate formation processes under the influence of both deterministic and stochastic radiative forcings. Our confidence in the validity and plausibility of the results obtained comes from the fact that, as shown in a number of research papers, the two-box ocean model and similar models have proven to be useful for studying climate change (e.g., [20,21,22,23,24,25,26,32,45]) and variability (e.g., [22,28,29,32,40,41,42,45,47,48]) on time scales from years to decades, since simulation results provide evidence that these models produce meaningful projections of GMST and its variability.

The application of other reasonable radiative forcing scenarios is unlikely to lead to qualitatively new results in the estimation of sensitivity functions, since deterministic forcing does not affect the power spectrum of GMST fluctuations (see Equation (10)), which is quite correctly reproduced by the two-box ocean model (e.g., [29,49]). It is clear that the results obtained from deterministic and stochastic models complement each other.

In conclusion, it is advisable to outline the problems that we intend to solve in the future. A large portion of GMST variability on interannual and decadal timescales is stimulated by random external forcings, also referred to as climate noise. This noise, being random in time, can spatially contribute to the formation of some persistent recurring features in the climate system state. As a result, under certain conditions, positive feedbacks that exist in multicomponent systems, such as the climate system, can be triggered, which will lead to the emergence of resonance phenomena. The model used in this study is two-component. In its current version, the thermocline fluctuations at the interface of the two layers of the model are not considered; therefore, the occurrence of resonance phenomena does not seem to be possible. Thus, the study of resonance, hysteresis, bistability and bifurcations is the area of our future research. At the same time, we intend to incorporate in the model time-varying feedback mechanisms, which will require further research. To study all of these problems in depth, we intend to use more complex climate models.