# Chapter 10Time-Dependent Response and the Ocean

Understanding and estimating the evolving temporal response of the system due to a time-dependent forcing are key problems in climate theory. In particular, the layers of air, land, and water have an effective heat capacity that can delay the response to a time-dependent stimulus. The column of air above land for forcing frequencies in the annual cycle range involves only a fraction of the atmospheric column's heat capacity. This effect can delay the warmest day of the year from the day of maximum heating by up to a month. Over open ocean, the same delay can be a whole season or mathematically a quarter of a cycle. A quarter cycle delay turns out to be the maximum when the effective heat capacity is large enough. The reason for the delay difference is that for the column of air, the relaxation time is small compared to the period of the periodic forcing. But the mixed-layer of the ocean has an effective heat capacity that leads to a radiative relaxation time of several years, depending on the depth chosen for its thickness. In a simple linear system of the seasonal cycle, this corresponds to a high frequency forcing (period of 1 year compared to a relaxation time of several years). Land surface response is probably confined to a meter or so into the soil and this, together with the air column, has a relaxation time of the order of a month when forced at roughly annual frequencies. A meter or two below the surface, say in a cave, we find the mean temperatures independent of season – the ground above filters out signals of period shorter than a few months. Lower frequencies tend to penetrate deeper with an effective depth inversely proportional to the square root of the frequency. We will explain this in Section 10.2 where a pure diffusive vertical heat transport is employed (valid in homogeneous soil and perhaps in the upper layers of the ocean).

A good opening example is the response to a sudden spike of forcing such as the negative spike from the dust veil following a volcanic eruption. This example is often used in linear analysis and is the so-called impulse/response function or the temporal Green's function.1 A second example is the response to an instantaneous doubling of . This forcing is proportional to the discontinuous Heavyside step function (0 for for ). In addition, we would like to see how the system responds to a periodic forcing. The latter has application to the seasonal cycle, the solar cycle, and the response to white noise forcing according to which amplitudes of the Fourier components of the forcing are spread with variance evenly across all frequencies (see Chapters 2 and 9). The damped-diffusion dynamics acts as a low-pass filter, yielding a response more concentrated in the lower frequency bands. The (simplified) forcing most directly bearing on the climate change problem is the response to a linear ramp forcing in time, beginning with a “cold start” at . As discussed in Chapter 2, this occurs when the concentration of carbon dioxide is increasing exponentially starting at . The forcing is linear in time because the forcing due to the greenhouse effect (Chapter 4) is nearly logarithmic in the concentration of (see Section 4.10). It is an observed fact that the concentration is increasing at a rate of 0.5% . But other greenhouse gases are also increasing at about the same rate, and it is conventional to take the forcing to be one that is increasing at a rate of 1.0% year, which yields a doubling time of about 70 years. Many general circulation climate model results are presented for this prototype experiment. In fact, the so-called transient climate sensitivity is the amount of global average temperature change from the onset to the time of doubling. All of the models treated in this chapter are governed by systems of equations that are linear and with time-independent coefficients. These exercises constitute an approach to linear systems through probing with some representative forcings followed by examination of the system's response characteristics.

In this chapter, we examine a variety of models, each treating thermal participation of different levels of the ocean in a hierarchy of complexity (a nice early review of transient models is given by Harvey and Schneider, 1985). Most of the models in the chapter are for a planet covered completely with ocean, the exception being one (see Section 10.6) where the full two-dimensional land–sea geographical surface distribution is considered. The simplest model is one which only considers the mixed layer of the ocean. This is the layer that is stirred by wind stresses at the ocean–air interface. Its depth depends on the wind stress and the static stability of ocean at a particular point. It tends to be deeper in the Southern Hemisphere. According to the same factors, the depth of the mixed layer depends on season and location, but we will mostly be concerned with global averages and we will use different depths from one exercise to another (mainly because of the variety of depths used by different authors in the literature). Using time averages of a month or so, this idealization should hold reasonably well. A similar scheme is often used in common atmospheric general circulation model (GCM) experiments involving the seasonal cycle.

As with virtually every chapter in this book, we remind the reader that the models studied are highly idealized. The ocean models considered here are far from realistic from the point of view of lateral heat transport, but their level of simplicity is roughly commensurate with the energy balance models (EBMs) introduced so far in the text. As with the other chapters, the models have the advantage that the reader can solve or at least understand the steps involved in detail. There are phenomenological coefficients (i.e., they are fitted to observations), but nothing is hidden.

## 10.1 Single-Slab Ocean

We begin with a couple of crude but surprisingly helpful exercises with the single-slab ocean. First is the response to a Dirac delta function , then a step-function forcing.

### 10.1.1 Examples with a Single Slab

Here we have in mind the global surface temperature response to a prototype volcanic eruption. Such an eruption has to blast material vertically with sufficient thrust to populate the stratosphere with debris (hygroscopic gases such as sulfur oxides or nitrous oxides can attract and adhere water molecules to form electrolyte solutions of sulfuric and sulfurous acid, also nitrous and nitric acids – the resulting aerosol particles reflect sunlight back to space). The light-reflecting debris can remain in the stratosphere for several years before coagulating to form larger aerosol particles and settling out (with some help from stratospheric circulation). The aerosol particles form a thin layer, inducing a negative forcing that cools the planet. The heating can be taken here as a negative delta function with its spike at time . Such a pulse causes a quick depression of the surface temperature (ideally globally distributed, but often the homogenization takes less than a year). The depressed temperature leads to a radiation imbalance and is followed by an exponential-like recharge that depends on the time constant of the Earth–atmosphere–ocean system. The recovery is long compared to the duration of the negative spike of forcing. The behavior of the temperature anomaly can be described by an energy balance equation:

10.1

where , is the heat capacity of a column of water whose horizontal cross section is 1 m, whose thickness is the depth of the mixed layer (typically 50–100 m), and W m K is Budyko's coefficient of the surface temperature for the outgoing radiation to space – it is sometimes referred to as the radiation damping coefficient. The portion of the heat capacity due to the atmospheric column is small compared to that of the mixed layer and is therefore neglected. This notation has been used often in previous chapters. To solve the ordinary differential equation, we multiply through by the integrating factor, and rearrange to obtain

10.2

Next we integrate each side from to , noting that , where means the value just infinitesimally below 0:

10.3

and finally, we have

10.4

As crosses the origin, the temperature abruptly falls by an amount . After the pulse at , the temperature “recharges” exponentially to zero as shown in Figure 10.1a. The mathematical steps for the adjustment from a flat forcing to a higher one of strength are almost identical to those just described for the impulse/response. The results are shown in Figure 10.1b. These experiments with GCMs are standard practice (e.g., Donahoe et al., 2015).

Next we take up the ramp forcing. As already mentioned, the dependence of outgoing terrestrial radiation on concentration is approximately logarithmic. Hence, if the amount of is increasing exponentially in time, the increase in radiative forcing is approximately linear in time. In this section, we imagine the system to be in steady state with the exponentially increasing being switched on at time . Atop the atmosphere of our all-ocean planet, we express this as

10.5

where is time in years, the Heaviside step function is defined by for , and W m/ with the equivalent doubling time for (conventionally taken to be 70 years).

This case and its response to different radiative-imbalance forcings can be found in many papers (e.g., Kim et al., 1992, and Watts et al., 1994). This model will serve as a convenient benchmark for us to compare with more complicated models later in this chapter. The slab world is defined by its energy balance equation:

10.6

where is the slab heat capacity per unit horizontal area (J m K) and is the net flux density of heat (W m) from the layers below and is an external time-dependent heat source applied at the surface (such as in the ramp case). The slab is assumed to have a high thermal conductivity such that its vertical temperature profile is homogenized in a short time compared to the radiative relaxation time (months to years). In the single-slab model, —there is no responding medium below the mixed layer. The time constant for relaxation is years, for a representative mixed-layer slab of about 80 m thickness.

The departure from steady state for is the solution to

where . We have a first-order linear nonhomogeneous equation to be solved for . The solution will consist of two parts: a homogeneous solution and a particular solution

10.8

A satisfactory particular solution (found by trial and error) is given by

10.9

which is the same as a solution except for the time lag . The homogeneous solution is given by

The integration constant is to be chosen so as to fit the initial condition, :

This partitioning of the solution is useful for physical insight (and for comparison later with more complex models where a similar trick will be used). Note that the particular solution as we have constructed it is the asymptotic form as becomes large compared to . We could think of the particular solution as the attractor, as the homogeneous solution always decays away because of the way we have partitioned it. The homogeneous solution which depends on the initial condition decays away in a characteristic time leaving only the particular or asymptotic solution. It is interesting that we have a familiar characteristic adjustment time for the solution to reach its asymptotic form (this adjustment structure and characteristic time will vary according to the model's complexity). The asymptotic form consists of a straight line corresponding to a slab with no heat capacity, but lagged by the characteristic time . These curves are illustrated in Figure 10.2. In the construction of the two terms, we have contrived to have the constant to be canceled at by the coefficient of the exponential in the homogeneous term. This decomposition involving the cancellation strategy holds in some of the more complex models that have a varying vertical structure below the surface.

Figure 10.3a shows the evolution of global and hemispherical surface temperatures based on an early version of the Geophysical Fluid Dynamics Laboratory (GFDL) coupled ocean–atmosphere GCM. After the transient dies out, one can trace a straight line to the time axis (shown as the heavy gray lines in the figure) to find lags of about 20 years for the global and Southern Hemisphere temperatures, while the Northern Hemisphere appears to have a lag closer to 15 years (Manabe et al., 1991). Other models, such as the Goddard Institute for Space Studies (GISS) coupled model of the same era, show similar behavior (Russell and Rind, 1999) as indicated by Figure 10.3b. In these experiments, the extrapolation is closer to 10 years for the lags. While these particular GCMs have evolved since these figures were produced, they indicate that the slab model concept of a decaying transient followed by the ramp of upward temperatures for large-scale surface temperature averages has some heuristic value.

### 10.1.2 Eventual Leveling of the Forcing

Next consider the increase of greenhouse gases to end at time followed by a leveling of the right-hand side of (10.7) to become where we have dropped the subscript on the heat capacity of the slab to simplify the notation. For times earlier than , the solution is the same as in the previous section if we continue to use the initial condition that . After time , the solution can be shown to be

10.12

For , we see that the solution finally settles down to the constant value of , which is just the value of the equilibrium-to-equilibrium change in global temperature for a forcing of . Note that the latter is the total accumulated forcing (proportional to the total log of concentration) in the atmosphere.

Figure 10.4 shows a numerical example of what happens after the forcing ramp is flattened, starting at 50 years. The forcing has been raised to a new constant level, but because of the lag as seen in the second term of (10.11), the temperature has not caught up with the forcing. The temperature continues to rise until it reaches the equilibrium-to-equilibrium value associated with a change in forcing from until . The adjustment after the leveling begins is exponentially damped to the new level with the same time constant .

The rise above the dashed line in Figure 10.4 is referred to as the commitment, because we would be committed to this much additional warming although we have zeroed the rate of rise of the forcing after time .

## 10.2 Penetration of a Periodic Heating at the Surface

To gain some insight into how heat is transmitted in a layered medium, we next consider a (vertically) thermally diffusive medium. It may seem odd to consider the more complicated continuous medium before looking at multiple slabs, but the analysis is not difficult and it helps to understand how the finite slab layers work. We take up first the case of periodic heating at the surface, because the frequency analysis of such a process gives an idea how far down the heat is conducted before the next round of cyclical forcing is applied. As usual, we allow the cooling of the surface to obey (see Chapter 2). Homogeneous terrestrial soil is the perfect candidate for this problem. We use it as a prototype for the ocean. Initially, we imagine a semi-infinite slab of soil (or ocean) (), hence is at the surface and is deep in the soil (or ocean). For simplicity, we take the thermal conductivity to be constant in space and time. For the temperature anomaly, (the departure of the temperature from its steady-state form for no external heating), we can write

10.13

where is the heat capacity of the medium (for seawater it is J m K; we ignore any density variation with depth), and , with units W K m, is the (macro)thermal conductivity. The term has units of s or 8000 m year. The vertical component of the heat flux density in J s m is given by

10.14

We will be imposing a sinusoidal forcing (heating) at the surface, , of where is the angular frequency and . We take the response temperature to be of the form .

At the top of the ocean, the boundary condition is

10.15

Moreover, the lower boundary condition is that must be finite as .

Since the problem is linear with real, constant coefficients, we may take the real part of the solution we find at the end. The ratio of the imaginary part of the solution to the real part will be the tangent of the phase lag in radians as a function of . Inserting this form, we find

where, in the last two equations, the prime indicates differentiation with respect to . Next we assume the form

10.18

where is a coefficient (likely complex) to be determined. Substitute this into (10.16) and (10.17). After noting that

10.19

we find

We dismiss the negative choice in the solution (10.20) because it would lead to unacceptable behavior as . Finally, we arrive at the complex amplitude:

10.21

Then, with and . The magnitude of the complex amplitude (after some algebra) is

10.22

where, in the last step, we substituted for . Note that the dependence is in the numerator, and it indicates that the warming signal is modulated by an exponentially decreasing function as . The characteristic depth of heat penetration is inversely proportional to the square root of the angular frequency of the driver.

10.23

Note that the penetration depth is also proportional to the square root of the thermal diffusivity . Low-frequency components can penetrate deeper than those with higher frequency. We can get an intuitive idea of how the heat is diffused from the surface toward the depths by looking at how a Gaussian distribution of diffusing substance spreads in a time . The distance, , of the one-standard-deviation width of a spreading Gaussian is (see Section 6.4). If we replace by , we find

10.24

At , we find that the amplitude squared (power or variance) can be written as follows:

10.25

This form has more power at low frequencies and for large , it monotonically decays as . Unlike the case of white noise, its integral over all frequencies is not bounded. In the statistics and engineering literature, it is known as “pink noise” ( with ). In reality, the forcing (white noise) is never really constant as , but rather it must cut off at some finite frequency. In this book, we usually think of the white noise drivers as atmospheric weather which has a timescale of a few days. Note that if there is no cooling to space (), the power spectrum diverges as . But note that the numerator, which is constant for , is proportional to and hence, for , causes convergence for .

Figure 10.5 shows the real part of the complex amplitude : Figure 10.5a is the amplitude without the damping factor and Figure 10.5b shows the same but including the damping factor. In this example, m year (i.e., below the surface) leading to a value of m. We can compute the phase lag by using W m (C).

In real soil, the medium is not usually homogeneous; hence, both and the local heat capacity are -dependent. Moreover, water in the column of soil can influence both of these parameters as well. Ignoring these complications for our purposes, we see that long timescales in the forcing can penetrate deeper in the soil. One might even posit that the effective heat capacity is proportional to the e-folding length scale, , which depends on the frequency . Heuristically, the penetration of the daily heating of the soil is much less than that of the annual cycle. Anyone who has toured a cave will find that the annual cycle of temperatures below ground is very small. We encounter the mean annual temperature instead. A similar effect can be noted with permafrost. Permafrost is mostly a few meters below the surface. For this reason, it is the mean annual temperature that matters at such depths. One could say the high-frequency components of the driving signal are filtered out by the vertical diffusion process.

An application of a more refined ground-heat-conduction model with a knowledge of the nonhomogeneous properties can be used as a proxy for past surface temperatures (see, e.g., Pollack and Huang, 2000; NRC, 2006). Proxy data generated from preexisting boreholes are very noisy, but when many of them are averaged together, the warming signal is clear.

While the vertical transport of heat in the ocean is hardly a diffusive medium, some of the intuition gained from the soil case is helpful. The case of upwelling of ocean waters combined with thermal diffusion gives the best2 vertical model as we will see in this chapter. Keep in mind that diffusion is a molecular process in which heat energy is transported by molecular collisions. Eddies in the fluid motion also transport heat energy and we imitate this transport mechanism by diffusion, but the real turbulent stirring process is more complicated than diffusion. However, in our ocean modeling slabs, diffusion and upwelling are appropriate models consistent with our overall strategy in this book.

## 10.3 Two-Slab Ocean

In this section, we consider a very simple two-slab model following the work of Gregory (2000)3 who provides background and motivation for the model. This was also the same class of model used by Held et al. (2010) to infer the long-term behavior of global temperatures after emissions are leveled or even stopped.

We denote the mixed-layer temperature as and the deeper-layer temperature as . The equation is for the anomaly from equilibrium response to the forcing at the top of the ocean . The flux density of heat between the two layers is proportional to their temperature differences with an exchange coefficient . We can introduce the two-component vector:

10.28

The homogeneous (differential equation) system can be expressed as follows:

10.29

with

10.30

Next, insert :

10.31

The eigenvalues can be found by solving the resulting quadratic equation for . The inverse are the decay times for the two eigenvectors. Figure 10.6 shows the two relaxation times corresponding to the inverse of the eigenvalues as a function of the exchange coefficient . When , the two slabs decouple and the shorter timescale becomes that of the mixed layer alone, 5.0 years, while that of the lower slab becomes indefinitely large. For large values of , the larger root tends toward the relaxation time for both slabs taken together, while the shorter time tends to zero. The solution corresponding to the lower timescale becomes a singular solution, because if (10.27) is divided through by , the time derivative term drops out, reducing the order of the system by one.

The actual solution for the homogeneous case can be written as a pair of exponentials that will decay according to the mix of the eigenvectors corresponding to the eigenvalues and .

where stands for the component of the unit eigenvector 1 (corresponding to ) in the direction of , and so on. The coefficients must be found from the initial conditions . We refer to the coefficients as the eigenmode amplitudes, as they indicate the strength of the eigen patterns represented by the unit eigenvectors . Once the coupling parameter is fixed, the eigenvector components are known constants. Note that the matrix is not symmetric, and therefore the eigenvectors are not orthogonal, but they still can be used as a basis set in two-dimensional space, provided they are not parallel. This means the temperature of each slab is a different linear combination of two exponential decays. For weak coupling, , this would be a short decay determined by the radiative decay time added to a component with long decay corresponding to the situation where there might be a lot of heat stored below that takes a long time to transfer to the upper layer where it can be radiated to space. This is the effect referred to as the recalcitrant climate effect (Held et al., 2010) that might occur after a long period of slow warming due to, say, greenhouse gases followed by a shutdown or even a leveling of the greenhouse gases. The recalcitrance is the long time required to release the heat from lower layers to the upper layer, where it can cool to space. If the coupling to lower layers is weak, the warming due to greenhouse gases does not reverse itself rapidly.

The eigenvalues are the roots of the quadratic equation derived from the null determinant .

10.33

The eigenvectors are given by complicated expressions (but easily handled by MATHEMATICA as analytical expressions, and numerically as well). It is convenient to use (normalized) unit vectors, . These unit vectors are not orthogonal. Hence, it is useful to introduce reciprocal vectors, , which4 are defined by

10.35

The eigenvectors have been normalized to have unit length, whereas the reciprocal vectors are not unit vectors. Next consider the relationship of decay times and and the thermal coupling . Figure 10.6 shows how the decay times (inverse of the eigenvalues) vary as a function of . The inverses of the eigenvalues are the decay time constants for the eigenmodes. As , one decay time is infinite and the other approaches 3.31 years, the decay times for the upper mixed-layer slab alone. In this limit, the corresponding eigenvectors are orthogonal as shown in Figure 10.7. The passage of is a singular limit because as the coupling goes to zero, the degree (dimension of the system) is reduced from two components to one. The result is that one of the timescales of the eigensystem tends to infinity.

But as is increased to very large values, we see that the decay time of becomes vanishingly small, while the other approaches 36.4 years, the decay time of a 550 m single, aggregated slab. In this limit, both slabs are tightly coupled and their entire mass is participating as a unit. As increases, the eigenvectors become more skewed at an angular separation of (see Figure 10.7). Figure 10.8 shows how the two unit vectors and rotate as increases from 0 to large values. Note that points in the negative direction in the plane.

For reference, we list the components of in terms of the :

10.37
10.38

and

10.39

Note that det is the cross product of and . It vanishes if these two unit vectors are orthogonal, which is the case as (see Figure 10.7).

### 10.3.1 Decay of an Anomaly with Two Slabs

Consider an anomaly or initial condition where K and K. (Note that this is equivalent to the upper slab being subjected to a delta function heating at time .) This anomaly will decay according to (10.32) after the constants are determined by the conditions

10.40

A numerical example is shown in Figure 10.9 and the late part of the decay is shown in Figure 10.10. In this numerical example, the upper slab has thickness 50 m and the lower slab 500 m, as in the previous subsection. We chose a fairly large value of the coupling parameter W m for illustrative purposes. Given this value of , the characteristic times (inverse of the eigenvalues) of the eigenmodes are 61.5 and 1.43 years. As the initial anomaly decays in the upper layer because of the fairly strong coupling, heat is transferred to the lower layer, which soon heats up to a warmer temperature than the upper layer. In this case, the cooling of the upper layer to the lower one is more efficient than the cooling to space. In weaker coupling (e.g., W m) scenarios, this does not happen until much later (20 years). Figure 10.10 shows a magnification where we see that both slabs are at nearly 0.02 K, even after many characteristic time lengths of the shorter e-folding time. The lower level is still warmer than the upper even at 50 years after the initial anomaly is launched.

Using as the heat capacity of the mixed layer, the decay time of an anomaly is related to sensitivity for a one-slab model, as the equilibrium sensitivity to doubling (see Chapter 2) is proportional to . Here we might be including in the feedbacks that are working at longer timescales, say, a few years. It is tempting to speculate that one could estimate the sensitivity by observing the decay time or equivalently the autocorrelation time. Aside from not having a good estimate of the slightly ambiguous term, , we have to worry about the coupling to the slab below the mixed layer. The coupling constant is similar to , in that it is not well known. Figure 10.11 shows several decays for unit anomaly and various values of the coupling parameter . In the figure, there are four decay curves that look similar to exponentials (but consist of a linear combination of two of them), whose time constants can be read from Figure 10.6. From Figure 10.11, we see that for very weak coupling the decay time is close to 3.31 years, the time for an uncoupled single upper slab. But as is increased, the decay time shrinks. The casual observer might think the sensitivity is less than that inferred by a value of . Lindzen and colleagues (Lindzen, 1994; Lindzen and Giannitsis, 1998) have used simple models similar to5 those in this chapter to show that, because of the unknown coupling parameter, we cannot get an unambiguous estimate of the sensitivity of climate. This seems to be borne out in our model as well.

### 10.3.2 Response to Ramp Forcing with Two Slabs

In this section, we consider the same model as in the previous section that has a mixed layer of thickness m (whose effective volumetric heat capacity is ).6 Below this thin layer lies another well-mixed slab whose bottom is at the thermocline. The thickness of the lower slab is m with effective volumetric heat capacity . This places the thermocline at about 550 m below the surface. In our model, there is no vertical heat transport below the thermocline. Several models of this genre were studied numerically by Harvey and Schneider (1985). The radiative decay time of the uncoupled mixed layer is taken to be 3.31 years and that of the lower slab to be 33.1 years. As before, let the two slabs transfer heat proportional to their temperature difference (warmer to cooler) with a coupling coefficient . The surface temperature is (taken to be the same as that of the mixed layer) with the usual heat flux densities and . The governing equations are (10.26) and (10.27).

Next consider obtaining a solution for ramp forcing applied to the upper layer only. The problem can be stated as follows:

with

10.42

We can now use

We can now insert (10.43) into (10.41) to find

10.44

where the overdot denotes time differentiation. Next, project out the equations for and by multiplying from the left by , then using (10.34)–(10.36). The result is two uncoupled ordinary differential equations, each of which is of the same form as the one-slab problem:

10.45
10.46

We can now solve each of these equations as in the solution of (10.7)–(10.11). The temperature in each slab can now be constructed for a given value of by inserting and into (10.43). Thus each slab's transient will consist of linear combinations of exponentials with the characteristic timescales for the eigenmodes. Without solving the system explicitly, we can conclude that after the transients have been exhausted, the simulated temperature in each slab gets onto the straight-line ramp.

The procedure just derived for two slabs can be used to solve for slabs by constructing the eigenvectors and their reciprocal vectors for the matrix . Further, we can use the technique to solve for any forcing on the right-hand side of the governing equation. For example, one might consider white noise forcing from above the ocean surface generated by weather. In this case, we insert for the forcing as we have done already in previous chapters.

## 10.4 Box-Diffusion Ocean Model

Lebedeff (1988) examined a box-diffusion model in which there is a deep layer below the mixed layer. The deep layer transports heat by the thermal diffusion mechanism. We follow his derivation but omit some detail because of the technical details making use of the Laplace transform.7

The model has a finite-thickness mixed layer (slab) on top with a continuum below:

10.47

where represents the departure from the steady state of the mixed-layer temperature, and is the thickness of the mixed layer (meters), is the volumetric heat capacity of seawater (4.18 J m K), is the flux density of heat entering the surface layer from below. In this last equation, is used in place of to include possibly other feedbacks such as ice-cap feedback.8

where is the (macro)thermal diffusivity of the ocean below the mixed layer (typically, in the range 1–10 m year)—here the macro-diffusion is considered to be driven by random eddies. is the departure of the ocean temperature below the mixed layer () from its steady-state profile. We have denoted the vertical coordinate in the lower layer as with the origin at the bottom of the mixed layer and having negative values below. The temperature in the lower levels is governed by

10.49

The boundary conditions are specified by (10.48) and that

10.50

We take the “bottom” of the ocean to be insulating. A null heat flux at the bottom of the layer at can be expressed as

10.51

Lebedeff's solution for step-function forcing is

10.52

where9

10.53

The formula for is not defined, but Lebedeff provides a solution for this value of the constant, . Plots for and 10 are shown in Figure 10.12. The formula for tells us that it is the ratio of two timescales: the time for a particle to diffuse across a distance of the mixed-layer depth to the radiative relaxation time for an uncoupled mixed layer: . In the numerator, the distance is diffused in accordance with . We discussed the spread of diffused heat in Section 6.4. Readers are referred to Lebedeff (1988) for more discussion and details of his solution. In his paper, he also derives simplified forms for long-term behavior of the solutions. Note that we have reversed the labels and in Lebedeff's paper to be consistent with our notation elsewhere in this chapter.

## 10.5 Steady State of Upwelling-Diffusion Ocean

An upwelling-diffusion (UD) ocean is yet another highly simplified model, but slightly more realistic than the purely diffusive one. It dates back to early works of Robinson and Stommel (1959) and Walter (1966). Many others have worked on this class of models as well (e.g., Hoffert et al., 1980; Harvey and Schneider, 1985; Morantine and Watts, 1994; Watts et al., 1994). One very attractive feature of this model is that a thermocline-like feature comes out naturally in the solution to the steady-state problem. The level of realism of this model is pretty low (for an up-to-date survey of ocean circulation theories, see, e.g., Huang, 2013). We adopt the UD model in the same spirit as other schematic approximations throughout this book.

The thinking behind the UD model is that throughout the ocean, eddies transport heat down gradient from the warm water above to the colder waters below. In addition cold water in the abyss is constantly supplied by downwelling of cold, dense waters in the polar regions. These volumes of colder and more saline volumes of water sink to the bottom or in some cases intermediate levels, spreading out horizontally mostly across the entire world ocean floor. This supply of cold water and its upward displacement of water above is the reason for the upwelling everywhere. It can also suggest us to make our lower boundary condition at be at a temperature of C or zero difference from the equilibrium solution when we deal with anomalies.

The local rate of change of thermal energy in the mixed layer is given by the equation:

10.54

where is the specific heat of seawater (in per volume units), is the effective heat capacity per unit horizontal area of the mixed layer, whose thickness is , is an external forcing, and as in Chapter 2, is the infra-red outgoing radiation to space, is the absorbed solar radiation flux density. Below the mixed layer is the ocean interior where the energy changes are governed by

where and are called the vertical diffusivity and the upwelling velocity, respectively. Note that is taken to be at the bottom of the mixed layer. In this formulation has dimensions /time (m year) with typical values 2000–6000 m year The upwelling velocity per unit heat capacity has dimensions length/time (m year) with typical values a few m year. These parameters suggest a vertical length scale ranging from 400 to 1000 m and a timescale years. The lower boundary condition is:

10.56

In obtaining the solution first consider the steady-state solution for the UD model with . Equation (10.55) yields after one integration

10.57

where we have set the integration constant to zero in order to satisfy the boundary condition at . Using we find

10.58

This very simple model of the world ocean leads us to a thermocline whose characteristic depth is below the mixed layer.

If the is suddenly doubled, the asymptotic solution (long times) is

10.59

where the subscript “p” stands for particular solution. The transient solution will be uniform over the all-ocean planet subjected to step-function forcing. Its solution is presented by Morantine and Watts (1994). They begin by finding the impulse/response function via a Laplace transform method. The steps in the derivation are more complicated than necessary for this book, but can be found in their paper.

### 10.5.1 All-Ocean Planetary Responses

One concern with the box-diffusion model raised by Hoffert and Flannery (1985) is that due to the response after long times, the entire column of water is heated during ramp forcing, whereas the UD model suggests that only the (effective) mass of water above the thermocline participates in the long-term (asymptotic) solution. The reason for this is that the upwelling cool water presses against the downward diffusing warm water and a match is found at the thermocline. We can find the so-called impulse/response function (also called Green's function) by finding the response to a forcing spiked at the origin, the Dirac delta function, . As mentioned in the last sections, this is equivalent to finding the decay of an initial anomaly in the mixed layer.

Watts et al. (1994) summarize several models including the pure diffusion model, the UD model, and a UD model that includes land masses in much the same way as the beach-ball model of Chapter 8. In this elegant survey, they used the Laplace transform method introduced by Lebedeff (1988) to provide analytical solutions. Here we show one example, the UD model for an all-ocean planet. The solution for this model is given by the following formula (Equation 9 of Morantine and Watts, 1994):

10.60

where is the departure of the mixed-layer temperature from its initial steady-state value, is time normalized by the upwelling timescale with defined below. The denominator of the left-hand side of the last equation consists of, the radiative forcing of the Earth–atmosphere system (e.g., the imbalance due to a sudden change [step function] in ), is the climate sensitivity parameter,10 typically approximately 2.0 W m K. The climate sensitivity parameter is the inverse of the change in global average temperature per unit of radiative forcing expressed in units, W m K. The special functions in the formula are the error function, erf(), and the complementary error function, erfc defined in Section 10.4. Other terms are defined below starting with and

10.61

An example of the response to step-function forcing is shown in Figure 10.13.

### 10.5.2 Ramp Forcing

Watts and colleagues also solve the problem with ramp forcing for the global ocean (no land as with the step-function forcing in the previous section). They approach the problem in much the same style as we have earlier in this chapter. The particular solution for the mixed-layer temperature responding to a forcing is

where , , and with in per unit volume units. In its asymptotic (in time) form, the lag behind the response with no ocean or atmosphere is the sum of the lag for the mixed layer and the lag associated with the effective heat capacity of the water between the thermocline and the surface. The transient solution is much more complicated, involving a mixture of exponentials, error functions, and so on. Note the very long adjustment time in Figure 10.13 of the UD ocean compared with the mixed-layer-only case. The transient or homogeneous solution contains two timescales, and . Morantine and Watts (1994) also point out in their discussion that naturally induced changes in (the response to which they give solutions) can cause changes in the global surface temperature that are as large as the global warming signal.

## 10.6 Upwelling Diffusion with (and without) Geography

The study by Kim et al. (1992) takes the past studies by Morantine and Watts (1990) and others (e.g., Wigley and Schlesinger, 1985; Lebedeff, 1988; Morantine and Watts, 1994) to include the land–sea geography in the sense of the previous chapters of this book, namely, the local surface is characterized by its effective heat capacity. The approach is similar with the starting point being the trial of the particular solution:

where is a lag that depends on the position on the globe. It turns out that , where can be thought of as the depth of the thermocline minus the depth of the mixed layer. By insertion of this last form into the governing equation, we find that satisfies

10.64

Over ocean, the terms on the right-hand side sum up to the radiative relaxation time for the entire thermocline. Over land, these terms add up to the relaxation time of an effective column of air. These quantities are known from the land–sea geography. This inhomogeneous equation has the form of the steady-state EBM, but with as dependent variable. The dependent variable is the response to the source terms on the RHS. It is damped by the term and smeared out by the diffusion term. Through this process of filtering the geography is entered. The solution for lag will be large over oceans, small over land, but the hard edges are removed over a length scale of the order of . Figure 10.14a shows the lag field for parameter values m year and m year. As expected, the lag is large over the oceans and small over land masses. The features are smoothed by the diffusion operator in accordance with the local length scale . Figure 10.14b shows a contour map of the asymptotic form for large values of time after the transients have died out, with the same parameters employing a deep thermocline of about 1100 m. Kim et al. (1992) used this deeper thermocline depth in these experiments to enhance the land–sea contrast of the signal. Figure 10.15 shows the global average temperature for a case where the more shallow thermocline is at 600 m. This figure illustrates that the scenario needs to start from rest in the middle 1700s because of the long timescales associated with the transient relaxation times. Figure 10.16 shows the anomaly response including both transient and particular solutions.

## 10.7 Influence of Initial Conditions

During the adjustment to the particular solution, there is a potentially long period during which the initial conditions are important because of the amount of ocean water that has to be heated. The common thread in all of these solutions is the separation of the homogeneous or transient solution from the particular or long-term solution. Recalling (10.11), we see that in order to make the solution consistent at , the coefficient in (10.10) has to contain a term that will cancel the lag term in the particular solution for the right-hand side to be equal to . This same cancellation trick comes up again in the UD ocean model for the step function as well as in the ramp-heating cases. In these cases, the initial condition must involve the profile in the deep ocean as well as just its value at the surface (see Figure 10.17).

In the case of the all-ocean planet, is a constant and from the terms in (10.63) one must cancel the exponential profile (of water mass) which has coefficient and the factor . Each of these terms will have a time constant of the order of the relaxation time of the whole thermocline (approximately few decades). As mentioned above, Morantine and Watts (1990) show that this has a complicated but closed-form analytical solution. The transient solution will also contain the timescale of the mixed layer alone, .

## 10.8 Response to Periodic Forcing with Upwelling Diffusion Ocean

In this section, we sketch the response of the UD Ocean to a periodic forcing at the surface. The approach is similar to that in Section 10.2, the difference being the inclusion of upwelling. The complete details (including land–sea geography) can be found in Kim and North (1992). This problem is important because it allows us to compute the spectral density of this model system under white noise forcing. We summarize the procedure here for the all-ocean planet. Toward the end, we quote some results from published papers that include geography.

The governing equations are the same as in the last section except that the forcing is now periodic in time. The procedure is similar to that in Section 10.2, but the upwelling term complicates the analysis. As in that chapter, we are interested in the steady-state periodic solution as opposed to the transient solutions which die out after several characteristic timescale intervals. We begin with the all-ocean planet. The details of this section are from Kim and North (1992). As in Section 10.2, we substitute for :

10.65

where is the angular frequency of the periodic forcing at the surface, and and are the real and imaginary parts of the depth dependence in the exponential. After inserting this form into the governing equations and the boundary conditions, we find a relation between the parameters , , and the angular frequency :

10.66

or equivalently:

10.67
10.68

At this point, it is convenient to normalize the variables to characteristic lengths (nominally, the depth of the thermocline):

10.69

Following this normalization, we drop the asterisks for simplicity. Then the previous equations become

10.70
10.71

Solving for , we have a quartic equation with four roots for a given :

10.72

Recall that in the case with no upwelling (), the penetration decreases as the driving frequency is increased (. Figure 10.18 shows us that when the frequency (squared) is zero, the penetration depth is at (the thermocline). As the frequency is increased (dotted line) from 0.5 to 3.0, the penetration depth changes from 1.25 to about 1.5 times the thermocline depth. In other words, the upwelling term induces quite a different penetration structure. This result that the effective depth is close to the thermocline might have been anticipated by the lag found in this model for ramp forcing in (10.62) where the thermocline relaxation time is added to the relaxation time of the mixed layer.

In the paper by Kim and North (1992), the model is developed with full geography as in the ramp-warming case of Section 10.6. In that paper, the seasonal cycle is updated with the new UD ocean and the model is subjected to white noise in space as well as time. The response to the noise forcing is shown in Figure 10.19 in the form of the spectral density. In constructing model-generated spectral densities, the strength of the forcing variance is arbitrary. We adjusted that strength to match the spectral density of the observations in the frequency range from a few decades to a period of a few months. Here we gain some insight especially when we consider low frequencies. At low frequencies, the models only differ by about a factor of 2, the UD model having more power as expected than the mixed-layer-only model. The UD model has a penetration of heat only just below the thermocline even at these very low frequencies. If upwelling is turned off in the infinitely deep ocean model, the power in this graph would be unbounded at lowest frequencies. This is a genuine difference between these two model structures.

## 10.9 Summary and Conclusions

In this chapter, we have considered a variety of oceanic models that are suitable for energy balance climate model (EBCM) studies. By moving through a list of models, each more complicated than the earlier one, we were able to uncover a number of features of the models that are important for determining whether they might be useful for a particular application. The simplest of these is the single-slab model. It satisfies the most elementary first-order linear ordinary differential equation with purely exponential-type solutions when it is unforced. There is a single timescale for this model, the radiative relaxation time of the slab. If found out of equilibrium, the solution returns to equilibrium with this time constant. If the model is forced at the top with white temporal noise, it reveals a spectral density identifiable as classical red noise.

Before going to multiple slabs, we treated the vertically continuous model with thermally diffusive transport of heat energy. If this is the only mechanism of vertical heat transfer and if the ocean is infinitely deep, we find some rather peculiar properties. The first is that the equilibrium solution is uniform in the vertical in disagreement with observations where a cooling of water from the surface to the thermocline is observed at a few hundred to a thousand meters. Another curious property of this model is that the depth of the penetration is strongly dependent on frequency of the surface heating. The penetration is inversely proportional to the square root of the frequency. Hence, low frequencies can penetrate very deep into the ocean, giving rise to a likelihood of very large power in the spectral density at low frequencies. In fact, if the radiation to space is made very small, the spectral density diverges at low frequencies. At high frequencies, there is another catastrophe, namely, the spectral density is of the pink-noise type (), which means it is not integrable – not as bad as white noise, but still not finite.

Numerous investigators have known of this problem for many years. The logical cure that does not defy the EBCM framework and spirit seems to be to introduce an upwelling term. The basis for such a radical idea is that cool saline waters in the polar regions downwell water that is cold compared to the waters above. This dense water slips under the bulk of the abyssal waters, gradually lifting the waters above at a rate of a few meters per year. This is the so-called UD model. But before taking it up, we showed how models with two layers work. In this case, the layers are coupled with flux densities that are proportional to the temperature difference. To be consistent with the discussion above, we introduced an upwelling term as well. The solutions to this class are somewhat like the single-slab model in that it is often convenient to separate the transient solutions from the particular solution, which can often be identified as an asymptotic form for large time. In this way, the problem has two parts. The transient solutions tend to die away, leaving the asymptotic particular solution to emerge after all the transients have died out. In the case of several slabs, there are time constants that can be related to each slab, but as they are coupled, the time constants emerge as the inverses of eigenvalues. Even if the coupling of the slabs is weak, the long timescales of the deep tend to linger in the adjustment times, but we will find that the upwelling term limits the deep eigenvalues to some extent.

The most commonly desired product of interest throughout the model hierarchy is the response to a ramp-like forcing that would result from an exponential increase of . In this case, up and down through the hierarchy, one finds that after the transients have died down, the temperatures in the response fall on a straight-line ramp. The response temperatures lag behind the situation where there is no thermal inertia at the surface. The lags can often be guessed by knowing the finite slabs or by solving for the continuum solutions. We note one difference right away: the UD model leads to lags in the asymptotic forms that suggest that the heat capacity of the thermocline might dominate. This is in contrast to the case where the upwelling term is omitted. In that latter case, the contribution of the very deep layers are more important. In the better-behaved UD models, the asymptotic solutions are not too hard to extract from the equations. The transient solutions are more difficult to obtain, especially if one seeks analytical forms. Often, as we have seen in the UD models, the adjustment process to new conditions involves warming waters well below the surface and this can take a very long time. Is the time taken for adjustment onto the ramp tens or hundreds of years? Modern coupled GCMs suggest that it may be only a few decades until we are well into the ramp regime.

Once we have learned how to generate the solutions for the UD model for the all-ocean planet, it is possible to extend it to the case where the land–sea geography is specified as in the earlier chapters of this book. Solutions are of course dependent on parameter choice and these are not uniquely determined (the nemesis of all climate models great and otherwise). One easily sees that in the ramp-like scenarios, the land masses lead the oceans and the lags are shorter on land as well, both results agreeing with data and intuitively satisfying. The spectral density of these UD models driven by white noise in space and time are interesting. They exhibit power at low frequencies, somewhat higher than comparably generated solutions with a mixed layer. But the power at very low frequencies is not much higher than one might expect. This interesting feature appears to be limited to the lack of deep penetration of power in these models (contrasting with a pure diffusion model). If true, this could mean that our current period of global warming might be forced (as opposed to natural variability driven by noise at the surface) as many investigators seem to believe.

## Exercises

1. 10.1 a. Determine the complete solution of a one-dimensional EBM with a sudden impulse forcing
10.73
2. where is the Dirac delta function.
3. b. What is the initial temperature right after the forcing is applied?
4. c. How long does it take for the initial temperature change to be reduced to an (1/e) level?
1. 10.2 a. Obtain the solution of a one-dimensional EBM with a step forcing
10.74
2. Here , and are constants, and the Heaviside step function is defined by
10.75
1. 10.3 a. Obtain the solution of a one-dimensional EBM with a ramp forcing
10.76
2. where , and are constants. Assume that .
3. b. Determine the lag between the solution and the linear forcing as time approaches infinity.
4. c. Let us further assert that the forcing levels at a constant value at time . Find the complete solution of the problem.
1. 10.4 Let us consider a one-dimensional EBM on top of a diffusive ocean with an infinite depth. Inside the diffusive ocean, temperature changes according to the heat diffusion equation

10.77

where is heat capacity and is thermal diffusivity. At the top of the ocean, energy balance is described as

10.78

where we essentially ignored the atmospheric column and assumed that the atmospheric temperature is identical with the surface temperature of the ocean.

1. a. Set up an equation together with the boundary condition for anomalous temperature driven by additional sinusoidal forcing .
2. b. Find the solution of the problem in Part (a).
2. 10.5 Let us consider a two-layer EBM in the form

10.79
10.80

where the subscript “m” and “d” denote mixed layer and deeper layer, respectively.

1. a. Set up the equation as a system of coupled linear equations.
2. b. Let us assume a homogeneous solution of the linear system in Part (a) in the form . Determine the parameters .
3. c. Plot the decay timescale (inverse of ) and the angle between two eigenvectors as a function of the coupling parameter . Use .
3. 10.6 Let us consider an energy balance system for a diffusive ocean with a mixed-layer at the top:

10.81
10.82

where

10.83

is the flux entering the mixed layer from below. Further, the boundary condition is written as

10.84

Derive the solution for a step forcing for .

4. 10.7 Let us consider an energy balance equation with a slab ocean:

where W year m C represents the heat capacity of the mixed-layer ocean which is approximately 70 m deep and is a ramp forcing. Find the complete solution.

5. 10.8 Let us consider an ocean whose interior temperature is governed by a diffusive process:

1. a. Explain the physical meaning of the boundary conditions.
2. b. Obtain the solution of the problem.
6. 10.9 Let us consider a UD ocean model in the form

where is the upwelling speed.

1. a. Calculate the vertical heat flux density for this model.
2. b. Boundary conditions for anomalous temperature are given by
3. where the ocean is assumed to be insulated from below. Non-dimensionalize the governing equation using and .
4. c. Find the complete solution of the problem.
7. 10.10 Let us consider an UD ocean coupled with a two-dimensional EBM as a top boundary condition:
where represents a point on the sphere, is local heat capacity per unit area, is local horizontal diffusion coefficient in the atmosphere and the mixed layer, is the slope of the best linear!it for the infrared emission to space, is upwelling speed, is vertical heat diffusion coefficient, is the heat capacity of seawater per unit volume, and is radiative forcing.
1. a. Show that thermal flux at the bottom of the mixed layer is given by
2. b. Let us consider a ramp forcing , where is the Heaviside step function. Let us consider the particular solution of the problem in the form
3. Show that the lag satisfies
4. where and
5. c. Let us consider the radiative forcing in the form
6. Let us consider the solution in the form
7. where is the position-dependent temporal lag of the solution, and are constants determining the vertical structure of the solution, and is a complex-valued amplitude. Determine and as a function of .