Thesis Status Report: Behavior of a Simple Coupled Numerical Model of Mid-latitude Climate Variability

Jason Goodman

1. Introduction

This document describes the progress I've made in my thesis project between June 1998 and February 1999. I begin by recapitulating the proposed work for this time period, and then go on to describe what I've actually done, and where this work leads.

In pre-thesis work, Dr. Marshall and I proposed a simplified physical mechanism which produced interannual and decadal variability in the middle-latitude coupled atmosphere-ocean system, finding a coupled growing perturbation which arises through thermal forcing of the atmosphere by sea-surface temperature and feedback onto the ocean via windstress. Preparations for publication of this work have been completed; an article will appear in the February 1999 Journal of Climate - this work will be referred to as ``GM'' in what follows. A schematic form of the GM equations is

$$J(\psi_a,q_a) -\lambda ( \theta_a - \theta^{*}(\mbox{SST}))$$ = (1)

$$\frac{\partial}{\partial t} \mbox{SST}+ J(\psi_o,\mbox{SST}) \mathcal{F}(\psi_o,\psi_a)$$ = (2)

$$\frac{\partial}{\partial t}q_o + J(\psi_o,q_o) \tau(\psi_a)$$ = (3)

Here, q_a and $\psi_a$ represent atmospheric PV and streamfunction, and q_oand $\psi_o$ the same for the ocean. The atmosphere is forced to relax to a temperature set by the sea surface; SST is forced by fluxes across the air-sea and mixed-layer interfaces, and the ocean is forced by a surface windstress. GM uses 2 levels in the atmosphere and 1.5 levels in the ocean to represent these dynamics, and then fully specifies the physics on the right-hand side, forms linear perturbation equations about a spatially-uniform basic state, and solves to find exponentially growing waves in the coupled system.

The analytically-soluble ``chalkboard'' model of the atmosphere-ocean system discussed in GM has quite a few drastic simplifications. Among the most drastic (and easiest to address) are that the background atmospheric winds are purely zonal and of constant strength, and that the ocean domain is unbounded. I proposed to assess the validity of the basic physical mechanism described in GM in more realistic geometries by discretizing the equations and solving them numerically to find eigenvalues and eigenvectors of the coupled system. Eigenvalues with positive real part correspond to growing modes of the coupled system.

In GM and in the proposal, we noted that the atmospheric pressure patterns which made up the coupled mode were those which corresponded to nearly-stationary free planetary waves - that is, solutions to the unforced atmospheric evolution where wave speed nearly balances flow speed. For a realistic background flow in a 3-layer QG model, Marshall and Molteni (1993) (M&M hereafter) were able to find ``neutral vectors'', free modes which were nearly stationary. Several of these patterns beak a marked resemblance to the NAO and PNA patterns. Like GM's modes and like observed patterns, the M&M neutral vectors have a barotropic structure. We hypothesized that in the coupled model with a wavy background flow, it would be these ``neutral vectors'' which coupled to the ocean to produce interannual variability.

In my proposal, I claimed that the linearized perturbation form of the GM equations could be finite-differenced to form a large set of ODEs, as follows:

$$\frac{\partial}{\partial t}\mathbf{A}\left[\begin{array}{c}\m... ...\Psi}_a \\ \mbox{\bf SST} \\ \mathbf{\Psi}_o\end{array}\right]$$ (4)

where $\mathbf\Psi_a$, $\mathbf\Psi_o$, and SST are vectors holding the values of the model variables at each of N gridpoints, and A and B are large, sparse matrices which can be specified to permit coastlines and non-rectilinear background flows in atmosphere and ocean. The eigenvalues and eigenvectors of this system with largest positive real part could be identified by the Arnoldi method. I set out to develop Matlab code to generate these matrices.

2. Model Development

It turns out that (4) above is not the correct way to state the problem. The atmospheric state of the GM model is slaved to its SST forcing: on the long timescales under consideration here, its local time tendency is taken to be zero. $\mathbf\Psi_a$ is determined by a diagnostic equation at each timestep, e.g.:

$$\mathbf{M} \mathbf{\Psi}_a = \mathbf{F}\cdot\mbox{\bf SST}$$ (5)

The tendency operator as implemented is an algorithm rather than an explicit matrix: it first solves (5) to find $\Psi_a$, then uses this atmospheric state to force SST and ocean streamfunction. This linear operation on SST and $\Psi_o$ could be written as a matrix, but the matrix would be large and dense, and thus difficult to use and store. The current code handles only the ``entrainment-only'' SST parameterization discussed by GM: in this case, SST is also a diagnostic linear function of $\psi_o$, so the entire model system can be written as

$$\frac{\partial}{\partial t}\mathbf{\Psi}_o = \mathcal{L}\left(\mathbf{\Psi}_o\right)$$ (6)

The most time-consuming part of this research was correctly formulating the M matrix in (5) to correctly discretize the model equations while exactly conserving PV. Since solving directly for the steady state is the same as finding the solution after an infinite amount of time has elapsed, one can see that even a tiny ``PV leak'' will have disastrous consequences.

I have had many weeks worth of difficulty using the Arnoldi eigensolver supplied with Matlab on this system. While the method should be applicable to sparse matrices generally, and also accepts linear operator functions like $\mathcal{L}$ in (6), I have been unable to get convergence with this procedure. ¹

However, I have had success with a ``poor man's eigensolver'' - straight forward integration of the time-dependent coupled system (6). After sufficiently long integration, only the fastest-growing (or least-damped) mode will remain. This is horribly inefficient, but it works.

3. Model Behavior

Illustrations of the time-evolution of several model configurations are attached. Figure 1 shows the evolution of a simple channel model with zonal extent equal to the 50$^{\circ}$ latitude line, no meridional boundaries in the ocean, and a uniform zonal background flow of 14 and 5 m/s in the upper and lower layers of the atmosphere, respectively. The model configuration is essentially identical to the analytically-soluble case discussed in GM. The model is initialized with smoothed random values of unit variance at each gridpoint. A snapshot of the model state is shown at the top of the figure, superimposed with contours of the background flow. A Hovmoeller diagram taken along a zonal slice through the `x' in the top figure is shown at bottom left, and a time-series of the model state at the `x' is shown at bottom right. The fastest-growing mode has zonal wavelength 8000 km, meridional wavelength 6400 km, an oscillation period of 6 years, with a similar e-folding timescale. These results are almost indistinguishable from GM.

However, growth is extremely sensitive to the amplitude of the zonal winds. The second figure shows the same configuration with a reduction of 1 m/s in the upper-level winds. Coupled growth ceases; the slow decay and westward propagation of perturbations is consistent with viscous oceanic Rossby wave dynamics without atmospheric feedback. Why does this happen? In GM, we found that the wavelength of the fastest-growing mode was highly sensitive to U₁, U₂, $\beta$, and meridional wavenumber l. In figure 1, the channel length is an integer number of wavelengths of the fastest-growing mode. In figure 2, this wavelength has changed, so the fastest-growing mode cannot meet the periodic boundary conditions, and thus cannot exist. In a spherical geometry, the coupled wave's meridional wavenumber l could change with U to allow the zonal periodic constraint to be met: however, in this narrow channel l is also highly restricted.

At the outset of this research, we noted that the results of M&M (date?) demonstrate that a wavy background flow can often create ``potential wells'' which trap planetary wave energy, localizing atmospheric perturbation modes. We speculated that localization might make the periodic boundary constraint which restricts the growth rate in figure 2 less important.

Figure 3 shows the model behavior for the (rather unrealistic) case of a small-amplitude wavenumber-1 ``wavy background flow''. The contour lines in the upper panels show streamlines of the background flow. It can be seen that while shape of the mode is altered slightly, growth still occurs. Also, no significant wave-energy concentration of the sort described above can be seen. Reducing the mean upper-level winds by 1 m/s (figure 4) eliminates spontaneous growth as it did in figure 2, but a coupled wavenumber-3 pattern is still visible as the least-damped mode (unlike figure 2). Note in the Hovmoeller diagram that wave energy seems to be concentrated at particular longitudes, rather than propagating evenly around the globe as it did in figure 2.

These flows are not very wavy: however, the atmospheric model appears to be numerically unstable for background flows with closed PV contours, giving gridscale noise where $\frac{\partial}{\partial y}\bar{q}$ changes sign. The reason for this remains unclear.

The second issue to be addressed with this model is an ocean with meridional boundaries. Figure 5 shows the model run with the same parameters as figure 1 but with an ocean basin covering only 1/3 of the model domain. We see that the perturbations decay extremely rapidly, with an e-folding timescale of about 5 years. Ocean streamfunction perturbations propagate westward as Rossby waves; when they strike the western boundary, they are dissipated and destroyed in the boundary current. Perturbations thus have a finite lifetime. The atmospheric response to these perturbations excites a new perturbation of opposite sign eastward of the dying one; however, the excitation mechanism is insufficient to build a new anomaly as strong as the old one before the old one dies, and thus perturbation energy decays away with time. The ratio of coupled growth e-folding timescale to the time required for a Rossby wave to cross the basin is the crucial parameter here.

A model run (not shown) with two basins instead of one gives similar highly damped results; in addition, the global pattern of atmospheric response ensures that the wave patterns in both basins vary in perfect synchrony. This seems pretty unlikely to occur in real observations - although a quick check shows a slight anticorrelation between the NAO and PNA pattern winter timeseries.

Finally, we come to the most interesting case. In the run shown in figure 6, I have put a strong regional wavy pattern in the background flow (the wavelength is consistent with observations of the stationary wave in the upper troposphere above the Atlantic), and constrained the ocean to a 6000 km basin. The wavy flow has a confluence over the western shoreline and a diffluence on the eastern side of the basin.

The atmospheric response is still global, but it is now phase-locked over the ocean. The atmospheric Hovmoeller diagram shows a an alternating sequence of highs and lows; the plan view shows a dipole pattern meridionally. This addresses one criticism of the GM model: observed interannual atmospheric variability fluctuates in place rather than propagating westward. The mode is still highly damped (due to the presence of the western boundary) the characteristic wavelength and oscillation frequency of the basic coupled interaction are visible.

In summary, then, forward integrations of the model with a wavy background flow show that

1.

An ocean of finite extent causes rapid damping of the GM coupled mode. This is a robust result with a huge impact on the GM mechanism.

2.

The requirement that an integer number of waves fit within the re-entrant channel imposes a strong constraint on the coupled growing mode: this makes growth highly sensitive to the background flow strength.

3.

A wavy background atmospheric flow changes the shape of the coupled mode, and may cause phase-locking in the atmosphere. However, the pattern is still global, even when the background flow causes phase-locking. It is not yet clear whether the presence of a wavy background flow can relax the ingeger-$\lambda$ constraint in (2) above, as we had hoped.

4.

The model is numerically unstable when the atmospheric background flow is strongly wavy; this, along with the reentrant channel geometry makes it very difficult to compare this model's atmospheric response with the Marshall & Molteni model, or with observations.

   
4. Is the Coupled Mode's Atmospheric State a ``Neutral Vector''?

We discussed in section 1 the hypothesis that in a wavy background flow, it is the ``neutral vectors'' of the atmosphere which produce a coupled atmosphere-ocean interaction. While the characteristic patterns of the present model's atmosphere look nothing like those of Marshall and Molteni, and are rather unlike observations, we can still use this model to find out if the atmospheric part of the coupled mode in this model is a ``neutral vector'' response.

Marshall & Molteni begin with the unforced QG planetary wave equation

$$\frac{\partial}{\partial t}q + J(q,\psi) = \mathcal{D}(\psi)$$

and discretize and linearize it, writing it matrix form:

$$\frac{\partial}{\partial t}\mathbf{Q}\Psi + \mathbf{M}\Psi = 0$$ (7)

They define the neutral vectors as those values of $\mathbf{\Psi}$ which minimize the vector norm of the PV tendency²

$$\frac{\left\langle\frac{\partial}{\partial t}\mathbf Q \Psi,\... ...\right\rangle}{\left\langle\Psi,\Psi\right\rangle} = \lambda^2$$

Substituting from (7),

$$\frac{\left\langle\mathbf M \Psi,\mathbf M \Psi\right\rangle}{\left\langle\Psi,\Psi\right\rangle} = \lambda^2$$ (8)

$$\frac{\left\langle\Psi,\mathbf M^{\dagger}\mathbf{M} \Psi\right\rangle}{\left\langle\Psi,\Psi\right\rangle} = \lambda^2$$

This can be minimized by looking for the eigenvectors $\Psi_n$ of $\mathbf{M}^{\dagger} \mathbf{M}$ with the smallest eigenvalues $\lambda_n^2$. These will be the modes with the smallest time tendency.

In this work, I have been solving for the steady atmospheric response to surface forcing by using the relation

$$J(q,\psi) = \mathcal{D}(\psi) + \mathcal{F}(\mbox{SST})$$

This equation can be discretized:

$$\mathbf{M} \Psi = \mathbf{f}$$ (9)

(see 5), where M is the same in both (7) and (9). This means we can re-interpret (8) as

$$\frac{\left\langle\mathbf f,\mathbf f\right\rangle}{\left\langle\Psi,\Psi\right\rangle} = \lambda^2$$

The $\Psi_n$ which minimize $\lambda^2$ are also the same patterns which maximize the magnitude of streamfunction response $\left\langle\Psi,\Psi\right\rangle$ for a given magnitude of thermal forcing $\left\langle\mathbf f,\mathbf f\right\rangle$. Thus, the neutral vectors of the free planetary wave problem are the same as the maximal excitation vectors of the steady forced problem.

In this light, our hypothesis that the neutral vectors are the ones which couple to the ocean makes perfect sense, for it is the neutral patterns which will respond most strongly to thermal forcing.

To test this idea, I have computed the first few neutral vectors for the Mmatrix corresponding to figures 1 and 6, the zonal channel model and the ``localized wiggle'' model. These are shown as figures 7 and 8.

The first six neutral vectors of the zonal flow case have constant barotropic zonal flow with increasing meridional wavenumbers. Since they have no zonal structure, advection of PV is identically zero, and the waves are affected only by the dissipation in the model. The first dynamically active modes are the eigenpair 6 and 7, which have a structure identical to that found in the coupled run (compare figure 1). Note that these modes are equivalent barotropic, like those in both GM and M&M.

The first neutral vector of the localized-wiggle case is also uninteresting; the second and third, however, look very much like the behavior of the coupled run in figure 6. (The upper part of Figure 6 shows only a snapshot of the coupled model; as it evolves, it tends to oscillate between a state resembling eigenvectors 2 and 3.)

Thus, we see that the neutral vectors are indeed the prominent patterns acting in the coupled model.

We can go one step further and answer the very interesting question: ``What patterns of thermal forcing would excite these neutral vectors''? This would enable us to find out what characteristics of the ocean state are relevant in determing the atmospheric flow. The forcing patterns an be found quite simply by plugging the neutral vectors $\Psi_n$ into (9) to obtain a set of forcing vectors f_n.

Figure 9 shows, at top right and left, eigenvectors 2 and 3 from the ``localized-wiggle'' atmosphere. Below that, I have plotted the projection of the forward-run coupled model's atmospheric state onto the first five atmospheric neutral vectors. Since the model's total energy is dissipating, I have multiplied these projections by a factor of $e^{t/\mbox{\small 2.5 yrs}}$ to offset their exponential decay to illustrate long-term trends. Note that $\Psi_2$ and $\Psi_3$ oscillate together about 60$^{\circ}$ out of phase. All other modes decay rapidly to zero.

Why do they do this? I have computed the first few thermal forcing eigenvectors f_n as well, and computed their projection onto the model's SST. This is shown in the second line plot. As expected, the timeseries are identical. Thus, we can say that ``the atmosphere looks like $\Psi_2$ (top left) at t=32 years because SST looks like f₂ (bottom left); the atmosphere looks like $\Psi_3$ (top right) at t=30.5 years because SST looks like f₃ (bottom right).''

Figure 10 shows how the changing SST pattern determines the atmospheric response. Figure 10a shows the forward model's SST at 32 years, when the atmosphere displays a $\Psi_2$ pattern. In 10b, we multiply the SST state at each point in the model by the forcing eigenvector f₂: the projection is strongest in the blue areas at the northwest and southwest corners of the basin. 10c shows the SST at a slightly earlier time (30.5 years); note the SST blobs near the center of the basin. Correspondingly, in figure 10d we see that this projects onto forcing mode f₃ at a point nearer the eastern side of the basin. As SST blobs propagate from east to west, they ``light up'' f₂ and f₂ sequentially, exciting $\Psi_2$ and $\Psi_3$ in the atmosphere.

5. What next?

The simplicity and limitations of this 2-level atmospheric channel model model make it unable to reliably find the neutral vectors of a background flow based on observations, and the numerical problems which appear in some cases make it difficult to trust the output. However, the conclusion reached above -- that the ``neutral vectors'' of the are, in fact, the modes which couple to the ocean and provide windstress feedback -- is a general result of the linear algebra in section 4, and should be valid no matter the accuracy of the solution.

The atmospheric model used by Marshall & Molteni is more numerically sophisticated, being a 3-layer spectral model on the sphere rather than a 2-layer finite-difference channel. In addition, some of its neutral vectors are very strongly correlated with the PNA and NAO patterns in the observed atmosphere. It would be very easy and quite illuminating to solve the M&M version of (9) to find the thermal forcing function which excites these patterns. If these patterns project onto the observed surface heat-flux pattern associated with the NAO and PNA, we have evidence that the ocean state can play an active role in producing these patterns.

The $\mathbf{M}^{\dagger} \mathbf{M}$ operator can also be produced using a (linearized) GCM and its adjoint. The automatic adjoint compiler can produce code linearized about a specific forward model run - I'm not sure, but it should be possible to linearize about an arbitrarily climatology as well. One could thus look for neutral vectors and their associated forcing functions in a model with more complete physics.

This approach merges into the second topic in my proposal, ``Maximal Excitation of Variability Patterns'', in which I proposed to use an adjoint GCM to find the atmospheric forcing pattern which maximizes the projection of the resulting atmospheric state onto the NAO, PNA, and other observed patterns.

About this document ...

Thesis Status Report: Behavior of a Simple Coupled Numerical Model of Mid-latitude Climate Variability

This document was generated using the LaTeX2HTML translator Version 98.1p1 release (March 2nd, 1998)

Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

The command line arguments were:
latex2html -split 0 thesis_status.

The translation was initiated by Jason Goodman on 1999-03-31

Footnotes

... procedure.¹

Just before printing the final draft of this report today (Feb 24 1999), I figured out the solution to this problem. I have found the growing eigenmodes of the coupled model for several simple cases: they agree precisely with the analysis and figures below. I will show results when we meet tomorrow.

... tendency²

Actually, they try to minimize $\left\langle\frac{\partial}{\partial t}\Psi,\frac{\partial}{\partial t}\Psi\right\rangle / \left\langle\Psi,\Psi\right\rangle$; however, a pattern with very small streamfunction tendency must also have very small PV tendency, so the two patterns should look alike. There are conceptual advantages to M&M's approach and numerical advantages to the method described here.