Appendix B: Field of Dreams - The General Theory Behind DREAM Forcing and Solutions

B1. Quick guide

Whenever I try to explain the way DREAM is forced I seem to run into difficulties. Is it because it’s really complicated ? Is it because I’m bad at explaining things ? Is it because people are generally a bit thick ? I don’t really accept any of these reasons ! But I know it’s impossible to explain it in a short conference talk, and even in a longer presentation it is tempting to go into too much depth. But here we can relax and look at all the implications of a forcing method that I still believe to be simple and effective.

Let’s start with notation that is intuitive for people who were bought up in the geosciences. Consider the development of some variable, say the potential vorticity \(q\), in the real atmosphere:

\[ \begin{equation*} \frac{\partial q}{\partial t} + {\bf v}.\nabla q = {\cal F}(t) - {\cal D}(q) \qquad\qquad (B1) \end{equation*} \]

What we have here is tendency, advection, external forcing and state-dependent dissipation. We can simulate the advection with our model, and we can specify damping and diffusion and tune it as we please. The biggest challenge is to find a way of specifying the forcing \(\cal{F}\). We know what it is physically. It represents source terms like radiative heating, boundary fluxes or condensation. But it’s a tough job to specify it in a simple way.

What we know for sure is that if we ran the model without any forcing, from any realistic initial condition, it would quickly develop large systematic errors. So one prescription for the forcing might be a source term for q that corrects the average systematic error of an unforced model. Let’s assume that we are in a statistically steady climate state (i.e. we consider one season). And let’s also decide to impose an external forcing that time independent. This is not realistic, but it does lead to the great simplification that our external forcing is just the time mean of \(\cal{F}\). We can readily see this by taking the time mean of (B1)

\[ \begin{equation*} \overline{\cal F} = \overline{{\bf v}.\nabla q} + \overline{{\cal D}(q)} \qquad\qquad (B2) \end{equation*} \]

Evidently if we use this to force a model that develops with the same advection and dissipation, i.e.

\[ \begin{equation*} \frac{\partial q}{\partial t} + {\bf v}.\nabla q + {\cal D}(q)= \overline{{\cal F}} \qquad\qquad (B3) \end{equation*} \]

then any statistically steady integration of this model will also satisfy equation (B2). This does not mean that the model’s climatology will be identical to the observed climatology. There is more than one way to satisfy (B2), because at least the advection term is nonlinear. The time mean advection contains contributions from advection of the mean flow and transient fluxes. Let’s have a look:

\[ \begin{equation*} \overline{\cal F} = \overline{\bf v}.\nabla \overline{q} + \overline{{\bf v}'.\nabla q'} + {{\cal D}(\overline{q})} \qquad\qquad (B4) \end{equation*} \]

So you see the same forcing could balance different partitions between the first two terms on the right hand side, associated with different model climatologies. A model forced in this way may have systematic errors, just like any GCM. Only the time mean generalised total flux is constrained to conform to observations. Note that in (B4) we have assumed that the dissipation is linear.

It remains to find a way of calculating this time-independent forcing. If you’ve read Chapter 4, part 2 then you already know the answer. The time mean of \(\cal{F}\) is just the exact opposite of the time mean of all the tendencies that would be generated by the unforced model if it were initialised from a representative set of observed states. So the trick is to take the model and initialise it many times, from a range of initial conditions, without any forcing. Run it for one timestep for each initial condition. Take all those one-timestep tendencies and average them together. \(\cal{F}\) is just minus that.

B2. Forcing a simple GCM

I tried to keep the quick guide quick. If you’re still reading, it’s because you want more. Be careful what you wish for. We’re now going to switch to a more generalised notation in which the instantaneous state of the atmosphere is represented by a state vector \({\bf \Phi}=(\zeta,D,T,q,p_* …)\), which represents vorticity, divergence, temperature, specific humidity and surface pressure in the same basis as the model. The development of the observed atmosphere can then be written as

\[ \begin{equation*} \frac{d {\bf \Phi}}{dt} + ({\cal A} + {\cal D}){\bf \Phi} = {\cal F}(t) \qquad\qquad (B5) \end{equation*} \]

where \({\cal A}\) is a nonlinear advection operator. Let’s separate the state vector into time-mean and transient components.

\[ \begin{equation*} \frac{d {\bf \Phi}'}{dt} + ({\cal A} + {\cal D})(\overline{\bf \Phi} + {\bf \Phi}') = \overline{\cal F} + {\cal F}' \qquad\qquad (B6) \end{equation*} \]

As before, the time mean of this equation can be written

\[ \begin{equation*} ({\cal A} + {\cal D})\overline{\bf \Phi} + \overline{O({\bf \Phi}'^2)} = \overline{\cal F} \qquad\qquad (B7) \end{equation*} \]

This is the generalised form of (B4). We see that time-mean advection is balanced by transient eddy fluxes and forcing. We could move the second term to the right hand side and call it the “transient eddy forcing”. Subtracting (B7) from (B6) we obtain the transient budget with respect to the time mean

\[ \begin{equation*} \frac{d {\bf \Phi}'}{dt} + {\cal L}\Phi' + \left[O({\bf \Phi}'^2) - \overline{O({\bf \Phi}'^2)} \right] = {\cal F}' \qquad\qquad (B8) \end{equation*} \]

In (B8) all the terms have zero time mean, so there is no large cancellation. This means that the time-mean state is a realistic basis for perturbation experiments. The development of small perturbations is conditioned by the linear operator \({\cal L}\) , which itself depends on the time mean state. And the structure of these small perturbations may be relevant to observed transient systems. For deeper philosophical insight into these matters see Hall and Sardeshmukh (1998).

Now back to the practical problem of forcing a simple GCM. Let’s introduce a model, with a state vector \({\bf \Psi}\) in the same basis as \({\bf \Phi}\) . If it is forced with a constant source field, it will develop according to:

\[ \begin{equation*} \frac{d {\bf \Psi}}{dt} + ({\cal A} + {\cal D}){\bf \Psi} = {\cal G}_{cm} \qquad\qquad (B9). \end{equation*} \]

As before let’s set our simple GCM forcing \({\cal G}_{cm} = \overline{\cal F}\). To find it we run the model with no forcing for one timestep, providing us with a set of tendencies:

\[ \begin{equation*} \frac{d {\bf \Psi}}{dt} + ({\cal A} + {\cal D}){\bf \Psi} = 0 \:\:\: \Rightarrow \:\:\: ({\cal A} + {\cal D}){\bf \Psi} = -\frac{d {\bf \Psi}}{dt} \end{equation*} \]

Do this many times with a set of observed states \({\bf \Phi}_i\) as initial conditions, and take the average of all the tendencies:

\[ \begin{equation*} {\cal G}_{cm} = \frac{1}{n} \sum_{i=1}^n ({\cal A} + {\cal D}){\bf \Phi}_i \qquad\qquad. \end{equation*} \]

If we use this forcing to perform a long integration of the model we can compare our simulation with the dataset we used to generate the empirical forcing. And we find that this prescription for the forcing guarantees that the total generalised flux from the model simulation will be the same as in the observations, i.e.

\[ \begin{equation*} \overline{({\cal A} + {\cal D}){\bf \Psi}} = \overline{({\cal A} + {\cal D}){\bf \Phi}} \qquad\qquad. \end{equation*} \]

But for reasons already explained above, it does not guarantee that the simulated time-mean flow will be realistic, i.e. \(\overline{\Psi} \ne \overline{\bf \Phi}\).

Neither does it guarantee that the transient fluxes will be realistic. This is because the balance of terms in:

\[ \begin{equation*} ({\cal A} + {\cal D})\overline{\bf \Psi} + \overline{O({\bf \Psi}'^2)} = ({\cal A} + {\cal D})\overline{\bf \Phi} + \overline{O({\bf \Phi}'^2)} \qquad\qquad \end{equation*} \]

can be achieved differently on the two sides of the equation, as already discussed above directly in terms of mean flow and transient fluxes (equation B4).

So we have a forcing that is time-independent. It formally corrects the average initial systematic error that you’d get if you had no forcing, and it ensures that the generalised flux convergence in the model solution is identical to that in the observational dataset. But it does not guarantee a realistic climate, or realistic transient fluxes. The model will typically have systematic errors, like any other GCM, associated with compensating errors in mean and transient fluxes, and between transient fluxes on different timescales.

B3. Forcing a perturbation model with a fixed basic state

To maintain a climatological mean state against time-mean advection requires a combination of diabatic forcing and transient eddy fluxes. The sum of these two effects could be equated to the first term in (B7). So let’s do that, and define a new forcing:

\[ \begin{equation*} {\cal G}_{bs} = ({\cal A} + {\cal D})\overline{\bf{\Phi}} \end{equation*} \]

From (B7) we see that this forcing is the same forcing we had before (\({\cal G}_{cm}\)) plus the “transient eddy forcing” alluded to previously. And no, that’s not what “bs” stands for, it stands for “basic state”. What would happen if we forced the model with \({\cal G}_{bs}\) and initialised it with the climatology \({\bf \Psi} = \overline{\bf{\Phi}}\) ? Well, of course, nothing would happen ! The forcing would exactly cancel the tendency produced by the initial condition and the integration would proceed with no development. It’s actually much easier to find this forcing that to find the simple GCM forcing, because all you have to do is run the model for one timestep from a single initial condition. In this example that initial condition would be the time-mean state, but this can work for any basic state you want. What is the point of doing this ? Well, you can then study the development of of the solution due to perturbations, either in the initial condition or in the forcing.

If we add a perturbation \({\bf \Psi}_1\) to the initial condition, the development equation becomes

\[ \begin{equation*} \frac{d {\bf \Psi}_1}{dt} + {\cal L}{\bf \Psi}_1 + O({\bf \Psi}_1^2) = 0 \end{equation*} \]

We can make sure that \({\bf \Psi}_1\) is small (and remains small) we have a linear perturbation model

\[ \begin{equation*} \frac{d {\bf \Psi}_1}{dt} + {\cal L}{\bf \Psi}_1 = 0 \qquad\qquad (B10) \end{equation*} \]

The solution to this equation gives the normal mode structure associated with the climatology \(\overline{\bf \Phi}\) (or any other basic state, with appropriate \({\cal G}_{bs}\) ). If \({\cal L}{\bf e}_n(x,y,z) = \lambda_n {\bf e}_n (x,y,z)\) and \(\lambda_n = \sigma + i\omega\)

then \({\bf \Psi}_1 = {\bf e}_n(x,y,z)e^{(\sigma + i\omega)t}\) or \({\bf \Psi}_1 = \left[ A(x,y,z) \cos\omega t + B(x,y,z) \sin\omega t \right] e^{\sigma t}\) , an exponentially developing normal mode with a spatial structure that cycles periodically between A, B, -A, -B.

Instead of adding a perturbation to the initial condition we add also perturbation to the forcing (to obtain the linear response it is sufficient to multiply some realistic forcing perturbation by a small factor, and then scale the solution back up by the same factor for diagnsotics). If the forcing anomaly is added to the right hand side of (B10)

\[ \begin{equation*} \frac{d {\bf \Psi}_1}{dt} + {\cal L}{\bf \Psi}_1 = {\bf f}_1 \end{equation*} \]

the solution is:

\[ \begin{equation*} \Psi_{1j} (t) = \frac{f_{1j}}{\lambda_j} \left( e^{\lambda_j t} -1 \right) \end{equation*} \]

for jth projection of \({\bf \Psi}_1\) and \({\bf f}_1\) onto \({\cal L}^T\) and jth eigenvalue of \({\cal L}\). If the eigenvalues \lambda_j all have negative real parts \(\sigma\), the basic state is stable for the model’s current set of dissipation parameters. Under these circumstances the model can be used to integrate forwards in time and find asymptotic time-independent perturbation solution \({\bf \Psi}_1 = {\cal L}^{-1}{\bf f}_1\).

But basic states that have realistic zonal winds are usually unstable, and so the normal modes described above will usually emerge after about 20 days regardless of how the model is perturbed. Structures and growth rates of normal modes can be elegantly teased out of a long integration using the modefinder (Chapter 4, section 4a).

But by their nature such normal mode solutions have little to do with the details of a forcing perturbation, so unless they are themselves the object of study, they can be a nuisance. For stationary wave perturbation studies we often need a stable system, so that a time integration will converge to a steady response. An approximate way to get this is to damp the entire system for a selection of damping rates. The trick is to damp all degrees of freedom equally, so as not to interfere with the structures of the eigenmodes. The procedure is also outlined inChapter 4, section 4a, section 4i.

B4. A word on damping and restoration

As discussed in Chapter 1, section 2, a popular way to force simple advective models is through “restoration” to a radiative-convective equilibrium state or a reference climatology. We imagine this as a spring, which returns the atmospheric state to what it would be if there were no dynamical fluxes, or at least prevents a model from straying too far from a realistic state. It’s worth briefly exploring the connection between restoration forcing and our empirical forcing approach in which we specify the dissipation, but not the equilibrium state.

Let’s rewrite the model development equation (B9):

\[ \begin{equation*} \frac{d {\bf \Psi}}{dt} + {\cal A}{\bf \Psi} = {\cal G}_{cm} - {\cal D}{\bf \Psi} = R({\bf \Phi^*} - {\bf \Psi}) \end{equation*} \]

On the right we have restoration towards a specified fixed equilibrium state \({\bf \Phi}^*\) at rate \(R\). We can identify this state as \({\bf \Phi^*} = {\cal G}_{cm}/R\) and the rate \(R\) is simply the local damping rate \({\cal D}\) . So it appears that the two approaches are mathematically identical. Instead of specifying a damping rate and an ad-hoc restoration state, we specify a damping rate and an objective empirical forcing.

Strictly, for this equivalence to hold, our dissipation needs to be just a local damping, so \({\cal D}\) is linear and diagonal. The dissipation in DREAM actually has off-diagonal elements associated with diffusion. We could still calculate the associated restoration state if we wanted, it’s just hard to think of a scientific motivation for doing it.

B5. Nudging

Imagine you want to examine the influence of an observed sequence of events \({\bf \Phi}_i\) in a specified region on the large scale circulation. It is possible to artificially constrain the model within this specified region by nudging. This is another kind of restoration forcing but it does not correspond to any real physical process. Think of it as an additional term in (B9):

\[ \begin{equation*} \frac{d {\bf \Psi}}{dt} + ({\cal A} + {\cal D}){\bf \Psi} = {\cal G}_{cm} + \left( \frac{{\bf \Phi}_i - {\bf \Psi}}{\tau} \right) \end{equation*} \]

where \(\tau\) is the nudging timescale. The procedure for implementing nudging is outlined in Chapter 4,section 4b.

Note that you need to be careful if you try to use nudging as a linear perturbation. This is normally achieved by scaling down the forcing perturbation with some small factor \(\epsilon\). But in this case you can’t just multiply the nudging term by , because if the model state is always close to a specified basic state, then the nudging term above will just become an additional constant forcing, not a spring. So instead of nudging towards the state \({\bf \Phi}_i\) , we need to nudge towards a state \({\bf \Phi}_i^*\) whose departure from the climatology \(\overline{\bf \Phi}\) is scaled down \({\bf \Phi}_i^* = \overline{\bf \Phi} + \epsilon({\bf \Phi}_i - \overline{\bf \Phi})\).

If we do this, then the nudging term becomes a small perturbation forcing and a damping on the anomaly response. The nudged version of (B10) can be written as:

\[ \begin{equation*} \frac{d {\bf \Psi}_1}{dt} + {\cal L}{\bf \Psi}_1 = \epsilon \left( \frac{{\bf \Phi}_i - \overline{\bf \Phi}}{\tau} \right) -\frac{{\bf \Psi}_1}{\tau} \end{equation*} \]

Fig. B1: Say no more !

B6. Diagnosing nonlinear and transient forcing

There are many studies in which some kind of perturbation is added to a dynamical model, and the response is analysed in terms of the direct response to the perturbation, plus an additional response to the modified transient fluxes. The direct response might be linear or nonlinear (the time independent nonlinear component is sometimes called the “stationary nonlinearity’).

If the transient fluxes have been modified in a perturbation run, then this is surely a finite amplitude nonlinear solution. In this case the direct response to the implied change in “transient eddy forcing” can also be diagnosed, and used to explain how the total response is different to the direct response.

A lot of studies start to get a bit qualitative at this point. One can look at eddy fluxes, E-vectors, their divergences, the local tendencies due to the anomalies in these fluxes, etc. And one can reason that this is why, for example, the PNA pattern is stronger in the GCM than it is in the simple stationary wave model. This is all good, traditional, midlatitude dynamics.

The neat thing about DREAM is that it is relatively easy to be more quantitative. In fact, the direct response to any forcing anomaly is not just local. This is equally true for a transient eddy forcing anomaly as it is for the original forcing anomaly that triggered the direct response. With DREAM, you can diagnose every source, and separately calculate every response, both local and non-local.

Consider a simple GCM experiment with a forcing perturbation. First do a control run (B9) that yields a model climatology \(\overline{\bf \Psi}_c\) . Now do a perturbation run :

\[ \begin{equation*} \frac{d {\bf \Psi}_p}{dt} + ({\cal A} + {\cal D}){\bf \Psi}_p = {\cal G}_{cm} + {\bf f}_1 \end{equation*} \]

Take the time mean and subtract from the time mean of the climatology run, and you get:

\[ \begin{equation*} ({\cal A} + {\cal D})\overline{\bf \Psi}_p - ({\cal A} + {\cal D})\overline{\bf \Psi}_c = {\bf f_1} + {\bf TE} \end{equation*} \]

Where \({\bf TE}\) is time-mean difference between quadratic terms in the transients from the two runs. This is the transient eddy forcing anomaly, and it is easy to find for the entire model state vector. We know the time mean for both runs and we know the forcing perturbation. So we can work out \({\bf TE}\) using the same techniques as we used to find above. Now we know both \({\bf f}_1\) and \({\bf TE}\), we can use them together and separately in perturbation experiments about the fixed basic state \(\overline{\bf \Psi}_c\). You’ll probably find that the linear solution about this basic state to the forcing perturbation \({\bf f_1} + {\bf TE}\) is pretty close to the difference between the two GCM runs. I told you it was neat !

B7. Forcing the annual cycle

Some uses of DREAM require a direct correspondence to real-world events, like SSTs or nudging data. For these use cases the model simulation must have a simple way of changing with the seasons. The constant forcing described up to now needs to be extended to include a repeating annual cycle. If we use the notation tilde to represent such a cycle, which contains the annual frequency and its harmonics, then the forcing we seek can be expressed as \({\cal F} = \overline{\cal F} + \widetilde{\cal F}\), and we can consider a full spectrum of time variation in observed states \({\bf \Phi} = \overline{\bf \Phi} + \widetilde{\bf \Phi} + {\bf \Phi}'\).

Where prime now denotes any variation that does not repeat annually. If we substitute this into (B5) then our expression for the forcing (equivalent to B7) becomes

\[ \begin{equation*} \overline{\cal F} + \widetilde{\cal F} = \frac{d \widetilde{\bf \Phi}}{dt} + (\overline{\cal A} + \widetilde{\cal A})(\overline{\bf \Phi} + \widetilde{\bf \Phi} + {\bf \Phi}') + {\cal D}(\overline{\bf \Phi} + \widetilde{\bf \Phi}) \end{equation*} \]

The main change is the first term on the right hand side which must be calculated directly from the data. The rest is calculated from one-timestep unforced model runs as before. The complete procedure is described in Chapter 4, section 2d. The second term contains a large number of timescale interactions which are exhaustively explored in Hall, Leroux and Ambrizzi (2019). They use a dummy multiplier technique in the data to isolate each term, but this is beyond the scope of this user guide. If you just want to run with an annual cycle then it is sufficient to calculate \(\overline{\cal F} + \widetilde{\cal F}\) as a single diagnostic.