**Climate System Simulations: An Integrated, Multi-Scale Approach for Research and Decision-Making**

Ángel G. Muñoz, Alfredo Nuñez and Ramón J. Cova *Centro de Modelado Científico (CMC). Universidad del Zulia. Venezuela*

## **1. Introduction**

20 Numerical Simulations / Book 1

448 Computational Simulations and Applications

Strohbehn, J.W. & Clifford, S.F. (1967) Polarization and Angle-of-Arrival Fluctuations for a

Taylor, G.I. (1938). The Spectrum of Turbulence. *Proceeding of the Royal Society of London. Series*

Willis, G.E. & Deardorff, J. W. (1976). On the use of Taylor's hypothesis for diffusion in the

Zhu, X. & Kahn, J.M. (2002). Free-space optical communication through atmospheric

*and Propagation*, Vol. 15, No. 3 , (May 1967), pp. 416–421, ISSN 0018-926X. Tatarskii, V. I. (1971). *The Effects of the Turbulent Atmosphere on Wave Propagation*, McGraw-Hill,

ISBN 0 7065 06804, New York, USA.

2002), pp. 1293-1300, ISSN 0090-6778.

(October 1976), pp. 817 – 822, ISSN 0035-9009.

Plane Wave Propagated through a Turbulent Medium. *IEEE Transactions on Antennas*

*A, Mathematical and Physical Sciences,* Vol. 164, No. 919 (February 1938), pp. 476–490.

mixed layer. *Quarterly Journal of the Royal Meteorological Society,* Vol. 102, No. 434

turbulence channels. *IEEE Transactions on Communications*, Vol. 50, No.8, (August

Climate physicists define the Climate System in terms of the different interactions taking place between hydrosphere, atmosphere, cryosphere, biosphere and lithosphere. Being the associated general thermo-hydrodynamic partial differential equation system so complex for analytical solving, it is usually integrated numerically involving a good deal of computational resources.

Moreover, today different spatial scales involve different physical parametrisations, and each forecast horizon (few days, seasonal, annual, decadal and climate change periods) of interest deserves a special treatment, mostly defined by the predictability of the Climate System and the characteristic response time of the interacting components at the corresponding temporal scale. Thus, it is presently customary to distinguish between large-, meso- and micro-scale numerical models, their descriptions ranging from global climate state to basin availability of hydrological resources. An integrated, multi-scale approach considering the output of the various models is of vital to understand at the different levels the behaviour of the Climate System, and thus offer useful tools for decision makers and stake holders.

In this chapter, after introducing a few fundamental concepts and equations in Section 2, and a hierarchy flux of information among the models for different scales in Section 3, the methodologies involving downscaling executions are discussed in detail, regarding both scientific research and policy-making applications. In section 4 we explain the usefulness of conducting several simulations (realisations) to partially reduce the inherent uncertainties. Later, in Section 5 we offer a plausible way to put into operation all the components, presenting several examples based on the experience acquired through a regional initiative known as the *Latin-American Observatory for Climate Events*. Some concluding remarks and suggestions for future research are presented in Section 6.

## **2. Fundamentals**

In this section several fundamental aspects related to the definition of the Climate System, its temporal scales and the governing equations of the atmosphere and the oceans are discussed. While the formal definition in terms of the Climate Subsystems is standard, it is important to remark that different models tend to write the corresponding system of thermo-hydrodynamical equations in different ways, due to the employment of particular numerical methods, simplifications, coordinate systems or physical assumptions.

The altitude *z* (or the radial distance) is not the most convenient vertical coordinate for most applications (cf (Haltiner & Williams, 1980)). There is a general tendency to employ the pressure coordinates *p* or ln *p*/*p*0, being *p*<sup>0</sup> a reference level pressure (e.g. at sea level). Other common vertical coordinates are *σ* = *p*/*p*0, a close modification called *eta* coordinate, the potential temperature *θ* and the isoentropic vertical coordinate. For details the reader can

<sup>451</sup> Climate System Simulations: An Integrated,

It is possible to generalise the discussion considering a vertical coordinate *ξ*, which is a single-valued, monotonic function of the altitude. In general, *ξ* is a function *x*, *y*, *z* and *t*:

where *x* and *y* are horizontal cartesian coordinates. Thus, a scalar field *A* can be written as

Now, the partial derivative with respect *s* (a dummy parameter denoting *x*, *y* or *t*), reads

where the subscript *z* or *ξ* denotes the corresponding vertical coordinate. Since the vertical

This expression can iteratively be used with *s* = *x* and *s* = *y* to construct the gradient of *A*

*DtA* <sup>=</sup> (*∂tA*)*<sup>ξ</sup>* <sup>+</sup> **<sup>V</sup>** · ∇*<sup>ξ</sup> <sup>A</sup>* <sup>+</sup> ˙

Naturally, we also have to calculate the elements of Jacobian matrices, as we usually need to go from the rectangular coordinate system used for the computational domain, to the the real-world terrain-following coordinate system, and viceversa. The general curvilinear models are becoming more common today for they lead to better results at regional and

As mentioned before, here we follow the approach advanced in (Marshal et al., 2004). The equations representing a compressible and hydrostatic atmosphere (for the non-hydrostatic

*ξ* ≡ *ξ*(*x*, *y*, *z*, *t*) (4)

*A*(*x*, *y*, *ξ*, *t*) ≡ *A*(*x*, *y*, *z*(*x*, *y*, *ξ*, *t*), *t*) (5)

(*∂sA*)*<sup>ξ</sup>* = (*∂sA*)*<sup>z</sup>* + *∂zA* (*∂sz*)*<sup>ξ</sup>* (6)

(*∂sA*)*<sup>ξ</sup>* = (*∂sA*)*<sup>z</sup>* + *∂ξ A∂zξ* (*∂sz*)*<sup>ξ</sup>* (8)

∇*<sup>ξ</sup> A* = ∇*zA* + *∂ξ A∂zξ*∇*<sup>ξ</sup> z* (9)

∇*<sup>ξ</sup>* · **B** = ∇*<sup>z</sup>* · **B** + *∂ξ***B***∂zξ* · ∇*<sup>ξ</sup> z* (10)

(*∂tA*)*<sup>ξ</sup>* = (*∂tA*)*<sup>z</sup>* + *∂ξ A∂zξ* (*∂tz*)*<sup>ξ</sup>* . (11)

*ξ* = *dξ*/*dt* is the vertical velocity.

*ξ∂ξ A* (12)

*∂zA* = *∂ξ A∂zξ*, (7)

consult (Haltiner & Williams, 1980) and references therein.

Multi-Scale Approach for Research and Decision-Making

and the bidimensional divergence of the vector field **B**:

derivatives correspond to

For *s* = *t*,

microscale levels.

**2.3 Atmosphere equations**

so that equation (6) becomes

The total derivative in terms of *ξ* is

where **V** is the horizontal velocity and ˙

Nonetheless, for the sake of simplicity, we follow here an approach suggested by (Marshal et al., 2004), which enables a mathematical representation of both atmosphere and ocean dynamics via isomorphisms. This methodology has the advantage of providing an elegant way to better understand the underlying similarities between these subsystems by coupling physical quantities at the boundaries in a clear manner.

#### **2.1 The Climate System and its scales**

The Climate System (*S*) can be defined in terms of five subsystems, namely the atmosphere (*A*), hydrosphere (*H*), cryosphere (*C*), biosphere (*B*) and lithosphere (*L*) , and their complex, nonlinear interactions. These components can be described as non isolated, open and heterogeneous thermo-hydrodynamical systems, which are physically characterized by their corresponding mechanical (e.g. forces) and thermodynamical states (e.g. pressure, temperature, salinity). Actually *A*, *H*, *C* and *B* act as a cascading system involving mass, momentum and energy fluxes through their boundaries (Peixoto & Oort, 1992). There is a general consensus in defining *S* as

$$S \equiv A \cup H \cup \mathcal{C} \cup B \cup L \tag{1}$$

Due to the great complexity in the different interactions among the climate subsystems and taking into consideration the timescales of the involved phenomena, it is useful to define a hierarchy of internal systems in terms of the characteristic response time to stable state perturbations. First we consider the systems with shortest typical times, and the other components are considered to be part of the external systems. For example, the characteristic times associated to perturbations in the planetary boundary layer are in the order of minutes to hours, while for tectonic dynamic the timescale is in the order of tens of millons of years (Peixoto & Oort, 1992).

This hierarchy gives the possibility of defining the Climate System for a particular timescale in a simplified way without losing any important physical phenomena. If one is interested in weather forecast with timescales from a few hours to 1-2 weeks, then the only internal component is the atmosphere:

$$S \equiv A \tag{2}$$

considering in the boundary conditions the effects of the oceans, sea ice, land surface processess and other phenomena like being part of the external system. In a similar way, for studying the Climate System from months to decades, it is possible to write

$$S \equiv A \cup O \cup B \tag{3}$$

where the oceans *O* belong to the hydrosphere subsystem. Basically, equations 2 and 3 are normally used today for hindcast simulations, weather and seasonal forecasts and climate change scenarios.

#### **2.2 Coordinate systems**

It is customary to use latitude and longitude as rectangular coordinates on a tangent plane located at a certain point on Earth's surface or above it. However, the real world is better described through non-orthogonal, *terrain-following* coordinate systems, especially if we need high resolution domains and to consider the local orography. As one should expect from this discussion, the selection of the vertical coordinate is one of the most important issues in modern numerical models. In this section a general approach is presented, reviewing briefly some typical vertical coordinates.

The altitude *z* (or the radial distance) is not the most convenient vertical coordinate for most applications (cf (Haltiner & Williams, 1980)). There is a general tendency to employ the pressure coordinates *p* or ln *p*/*p*0, being *p*<sup>0</sup> a reference level pressure (e.g. at sea level). Other common vertical coordinates are *σ* = *p*/*p*0, a close modification called *eta* coordinate, the potential temperature *θ* and the isoentropic vertical coordinate. For details the reader can consult (Haltiner & Williams, 1980) and references therein.

It is possible to generalise the discussion considering a vertical coordinate *ξ*, which is a single-valued, monotonic function of the altitude. In general, *ξ* is a function *x*, *y*, *z* and *t*:

$$
\mathfrak{F} \equiv \mathfrak{F}(\mathfrak{x}, \mathfrak{y}, \mathfrak{z}, \mathfrak{t}) \tag{4}
$$

where *x* and *y* are horizontal cartesian coordinates. Thus, a scalar field *A* can be written as

$$A(\mathbf{x}, y, \xi, t) \equiv A(\mathbf{x}, y, z(\mathbf{x}, y, \xi, t), t) \tag{5}$$

Now, the partial derivative with respect *s* (a dummy parameter denoting *x*, *y* or *t*), reads

$$(\partial\_{\mathbf{s}}A)\_{\mathfrak{f}} = (\partial\_{\mathbf{s}}A)\_{\mathbf{z}} + \partial\_{\mathbf{z}}A \, (\partial\_{\mathbf{s}}z)\_{\mathfrak{f}} \tag{6}$$

where the subscript *z* or *ξ* denotes the corresponding vertical coordinate. Since the vertical derivatives correspond to

$$
\partial\_z A = \partial\_{\tilde{\zeta}} A \partial\_z \tilde{\zeta}\_{\prime} \tag{7}
$$

so that equation (6) becomes

$$(\partial\_s A)\_{\mathfrak{f}} = (\partial\_s A)\_z + \partial\_{\mathfrak{f}} A \partial\_z \mathfrak{f} \left(\partial\_s z\right)\_{\mathfrak{f}} \tag{8}$$

This expression can iteratively be used with *s* = *x* and *s* = *y* to construct the gradient of *A* and the bidimensional divergence of the vector field **B**:

$$
\nabla\_{\tilde{\xi}} A = \nabla\_{\tilde{z}} A + \partial\_{\tilde{\xi}} A \partial\_{\tilde{z}} \tilde{\xi} \nabla\_{\tilde{\xi}} z \tag{9}
$$

$$
\nabla\_{\tilde{\xi}} \cdot \mathbf{B} = \nabla\_z \cdot \mathbf{B} + \partial\_{\tilde{\xi}} \mathbf{B} \partial\_z \tilde{\xi} \cdot \nabla\_{\tilde{\xi}} z \tag{10}
$$

For *s* = *t*,

2 Numerical Simulations

Nonetheless, for the sake of simplicity, we follow here an approach suggested by (Marshal et al., 2004), which enables a mathematical representation of both atmosphere and ocean dynamics via isomorphisms. This methodology has the advantage of providing an elegant way to better understand the underlying similarities between these subsystems by coupling

The Climate System (*S*) can be defined in terms of five subsystems, namely the atmosphere (*A*), hydrosphere (*H*), cryosphere (*C*), biosphere (*B*) and lithosphere (*L*) , and their complex, nonlinear interactions. These components can be described as non isolated, open and heterogeneous thermo-hydrodynamical systems, which are physically characterized by their corresponding mechanical (e.g. forces) and thermodynamical states (e.g. pressure, temperature, salinity). Actually *A*, *H*, *C* and *B* act as a cascading system involving mass,

Due to the great complexity in the different interactions among the climate subsystems and taking into consideration the timescales of the involved phenomena, it is useful to define a hierarchy of internal systems in terms of the characteristic response time to stable state perturbations. First we consider the systems with shortest typical times, and the other components are considered to be part of the external systems. For example, the characteristic times associated to perturbations in the planetary boundary layer are in the order of minutes to hours, while for tectonic dynamic the timescale is in the order of tens of millons of years

This hierarchy gives the possibility of defining the Climate System for a particular timescale in a simplified way without losing any important physical phenomena. If one is interested in weather forecast with timescales from a few hours to 1-2 weeks, then the only internal

considering in the boundary conditions the effects of the oceans, sea ice, land surface processess and other phenomena like being part of the external system. In a similar way,

Basically, equations 2 and 3 are normally used today for hindcast simulations, weather and

It is customary to use latitude and longitude as rectangular coordinates on a tangent plane located at a certain point on Earth's surface or above it. However, the real world is better described through non-orthogonal, *terrain-following* coordinate systems, especially if we need high resolution domains and to consider the local orography. As one should expect from this discussion, the selection of the vertical coordinate is one of the most important issues in modern numerical models. In this section a general approach is presented, reviewing briefly

for studying the Climate System from months to decades, it is possible to write

where the oceans *O* belong to the hydrosphere subsystem.

seasonal forecasts and climate change scenarios.

*S* ≡ *A* ∪ *H* ∪ *C* ∪ *B* ∪ *L* (1)

*S* ≡ *A* (2)

*S* ≡ *A* ∪ *O* ∪ *B* (3)

momentum and energy fluxes through their boundaries (Peixoto & Oort, 1992).

physical quantities at the boundaries in a clear manner.

**2.1 The Climate System and its scales**

There is a general consensus in defining *S* as

(Peixoto & Oort, 1992).

**2.2 Coordinate systems**

some typical vertical coordinates.

component is the atmosphere:

$$(\partial\_l A)\_{\mathfrak{f}} = (\partial\_l A)\_z + \partial\_{\mathfrak{f}} A \partial\_{\mathfrak{z}} \mathfrak{f} \ (\partial\_l z)\_{\mathfrak{f}}.\tag{11}$$

The total derivative in terms of *ξ* is

$$D\_l A = (\partial\_l A)\_{\tilde{\xi}} + \mathbf{V} \cdot \nabla\_{\tilde{\xi}} A + \dot{\xi} \partial\_{\tilde{\xi}} A \tag{12}$$

where **V** is the horizontal velocity and ˙ *ξ* = *dξ*/*dt* is the vertical velocity.

Naturally, we also have to calculate the elements of Jacobian matrices, as we usually need to go from the rectangular coordinate system used for the computational domain, to the the real-world terrain-following coordinate system, and viceversa. The general curvilinear models are becoming more common today for they lead to better results at regional and microscale levels.

#### **2.3 Atmosphere equations**

As mentioned before, here we follow the approach advanced in (Marshal et al., 2004). The equations representing a compressible and hydrostatic atmosphere (for the non-hydrostatic

The boundary condition for integrating the hydrostatic equation (14) is

is the horizontal mean wind velocity. Figure 1 illustrates the case.

describing the terrain altitude in the lower boundary.

Multi-Scale Approach for Research and Decision-Making

The surface pressure evolves according to

and are parametrised in the numerical models.

*ω* =

at *z* = *η*, where we recall that *ps* is the sea level pressure.

The surface elevation evolves following

The boundary condition for the hydrostatic equation in this case reads

The boundary conditions are

where

**2.4 Ocean equations**

coordinates.

valid for *p* = *ps*, where *ps* is the sea surface level pressure and *H* stands for a scalar field

<sup>453</sup> Climate System Simulations: An Integrated,

*ps*

The hydrostatic equations describing the dynamics of an incompressible ocean in the

<sup>+</sup> *<sup>g</sup> <sup>p</sup> ρ*0

*Dt<sup>θ</sup>* <sup>=</sup> *<sup>Q</sup><sup>θ</sup>*

where *ω* = *Dtz* is the vertical velocity, *p* denotes pressure, *ρ* the density, *ρ*<sup>0</sup> a constant density reference, *S* the salinity and *Dt* = *∂<sup>t</sup>* + **V** · ∇*<sup>z</sup>* + *ω*(*∂z*) corresponds to (12) but in depth

Again, **F**, *Q<sup>θ</sup>* and *Qs* represent sources and sinks of momentum, heat and salinity, respectively,

<sup>−</sup>**<sup>V</sup>** · ∇*H*, *<sup>z</sup>* <sup>=</sup> <sup>−</sup>*<sup>H</sup>*

with *H* being a scalar field describing the bathymetry (−*H* is the sea bottom), *η* being the ocean surface and (*P* − *E*) being the net difference between precipitation and evaporation.

 *ps* 0

> *p ρ*0

<sup>&</sup>lt; **<sup>V</sup>** <sup>&</sup>gt;<sup>=</sup> <sup>1</sup>

Boussinesq approximation and using depth coordinates are (Marshal et al., 2004)

*∂z p ρ*0

*Dt***V** + *f***k** × **V** + ∇*<sup>z</sup>*

Φ = Φ*<sup>s</sup>* = *gH*, (20)

*∂<sup>t</sup> ps* = ∇ · (*ps* < **V** >) = 0 (21)

*dp***V** (22)

= **F** (23)

= 0 (24)

<sup>Π</sup> (27)

*DtS* = *Qs* (28)

*Dt<sup>η</sup>* <sup>=</sup> <sup>−</sup>(*<sup>P</sup>* <sup>−</sup> *<sup>E</sup>*), *<sup>z</sup>* <sup>=</sup> *<sup>η</sup>* (29)

*∂tη* + ∇ · (*H* + *η*) < **V** >= *P* − *E* (31)

*p* = *ps* (30)

∇*<sup>z</sup>* · **V** + *∂zω* = 0 (25) *ρ* = *ρ*(*θ*, *S*, *p*) (26)

Fig. 1. Mathematical description of the atmosphere boundary conditions. Source: (Marshal et al., 2004)

case see (Adcroft et al., 2011)) in pressure coordinates are

$$D\_t \mathbf{V} + f \mathbf{k} \times \mathbf{V} + \nabla\_p \Phi = \mathbf{F} \tag{13}$$

$$
\partial\_{\mathcal{P}} \Phi + \mathfrak{a} = 0 \tag{14}
$$

$$\nabla\_p \cdot \mathbf{V} + \partial\_p \omega = 0 \tag{15}$$

$$
\mathfrak{a} = \mathfrak{a}(\theta, p) = \theta \partial\_p \Pi \tag{16}
$$

$$D\_l \theta = \frac{Q\_\theta}{\Pi} \tag{17}$$

$$D\_t q = Q\_q \tag{18}$$

where *ω* = *Dt p* denotes the vertical velocity in pressure coordinates, *f* the Coriolis parameter, **k** a versor along the vertical coordinate, Φ = *gz* the geopotential, *α* the specific volume, *T* the temperature, *θ* = *cpT*/Π the potential temperature, Π = *cp*(*p*/*p*0)*<sup>κ</sup>* the Exner function and *q* the especific humidity. *cp* corresponds to the especific heat at constant pressure, *κ* = *R*/*cp*, with *R* the Rydberg's constant and *Dt* = *∂<sup>t</sup>* + **V** · ∇*<sup>p</sup>* + *ω*(*∂p*) is equivalent to (12) in pressure coordinates.

The terms **F**, *Q<sup>θ</sup>* and *Qq* represent sources and sinks of momentum, heat and humidity, respectively, and are parametrised in the numerical models. The boundary conditions are written as

$$\omega = \begin{cases} \ 0, \quad p = 0\\ \mathcal{D}\_t p\_{s\prime} \quad p = p\_s(\mathfrak{x}, y, t) \end{cases} \tag{19}$$

The boundary condition for integrating the hydrostatic equation (14) is

$$
\Phi = \Phi\_{\sf s} = \lg H\_{\sf s} \tag{20}
$$

valid for *p* = *ps*, where *ps* is the sea surface level pressure and *H* stands for a scalar field describing the terrain altitude in the lower boundary. The surface pressure evolves according to

$$
\partial\_t p\_s = \nabla \cdot (p\_s < \mathbf{V} >) = 0 \tag{21}
$$

where

4 Numerical Simulations

Fig. 1. Mathematical description of the atmosphere boundary conditions. Source: (Marshal et

*Dt<sup>θ</sup>* <sup>=</sup> *<sup>Q</sup><sup>θ</sup>*

where *ω* = *Dt p* denotes the vertical velocity in pressure coordinates, *f* the Coriolis parameter, **k** a versor along the vertical coordinate, Φ = *gz* the geopotential, *α* the specific volume, *T* the temperature, *θ* = *cpT*/Π the potential temperature, Π = *cp*(*p*/*p*0)*<sup>κ</sup>* the Exner function and *q* the especific humidity. *cp* corresponds to the especific heat at constant pressure, *κ* = *R*/*cp*, with *R* the Rydberg's constant and *Dt* = *∂<sup>t</sup>* + **V** · ∇*<sup>p</sup>* + *ω*(*∂p*) is equivalent to (12) in pressure

The terms **F**, *Q<sup>θ</sup>* and *Qq* represent sources and sinks of momentum, heat and humidity,

0, *p* = 0

*Dt***V** + *f***k** × **V** + ∇*p*Φ = **F** (13)

*∂p*Φ + *α* = 0 (14) ∇*<sup>p</sup>* · **V** + *∂pω* = 0 (15) *α* = *α*(*θ*, *p*) = *θ∂p*Π (16)

*Dtq* = *Qq* (18)

*Dt ps*, *<sup>p</sup>* <sup>=</sup> *ps*(*x*, *<sup>y</sup>*, *<sup>t</sup>*) (19)

<sup>Π</sup> (17)

case see (Adcroft et al., 2011)) in pressure coordinates are

respectively, and are parametrised in the numerical models.

*ω* =

The boundary conditions are written as

al., 2004)

coordinates.

$$<\mathbf{V}> = \frac{1}{p\_s} \int\_0^{p\_s} dp \mathbf{V} \tag{22}$$

is the horizontal mean wind velocity. Figure 1 illustrates the case.

#### **2.4 Ocean equations**

The hydrostatic equations describing the dynamics of an incompressible ocean in the Boussinesq approximation and using depth coordinates are (Marshal et al., 2004)

$$D\_t \mathbf{V} + f \mathbf{k} \times \mathbf{V} + \nabla\_z \frac{p}{\rho\_0} = \mathbf{F} \tag{23}$$

$$2\frac{p}{\rho\_0} + g\frac{p}{\rho\_0} = 0\tag{24}$$

$$\nabla\_z \cdot \mathbf{V} + \partial\_z \omega = 0 \tag{25}$$

$$
\rho = \rho(\theta, S, p) \tag{26}
$$

$$D\_l \theta = \frac{Q\_\theta}{\Pi} \tag{27}$$

$$D\_t S = Q\_s \tag{28}$$

where *ω* = *Dtz* is the vertical velocity, *p* denotes pressure, *ρ* the density, *ρ*<sup>0</sup> a constant density reference, *S* the salinity and *Dt* = *∂<sup>t</sup>* + **V** · ∇*<sup>z</sup>* + *ω*(*∂z*) corresponds to (12) but in depth coordinates.

Again, **F**, *Q<sup>θ</sup>* and *Qs* represent sources and sinks of momentum, heat and salinity, respectively, and are parametrised in the numerical models.

The boundary conditions are

$$\omega = \begin{cases} -\mathbf{V} \cdot \nabla H\_{\prime} & z = -H \\ D\_{t} \eta = -(P - E)\_{\prime} & z = \eta \end{cases} \tag{29}$$

with *H* being a scalar field describing the bathymetry (−*H* is the sea bottom), *η* being the ocean surface and (*P* − *E*) being the net difference between precipitation and evaporation. The boundary condition for the hydrostatic equation in this case reads

$$p = p\_s \tag{30}$$

at *z* = *η*, where we recall that *ps* is the sea level pressure. The surface elevation evolves following

$$
\partial\_l \eta + \nabla \cdot (H + \eta) < \mathbf{V} > = P - E \tag{31}
$$

Fig. 2. Mathematical description of the ocean boundary conditions. Source: (Marshal et al., 2004)

where

$$<\mathbf{V}> = \frac{1}{H+\eta} \int\_{-H}^{\eta} dz \mathbf{V} \tag{32}$$

Fig. 3. Isomorphic relations and coupling. Source: (Marshal et al., 2004)

and

where

atmosphere and ocean fluxes.

The corresponding boundary conditions are

*r*˙ =

*Rfixed* =

where *r* is the vertical, isomorphic coordinate. The total derivative takes the form

.

<sup>455</sup> Climate System Simulations: An Integrated,

Multi-Scale Approach for Research and Decision-Making

with ∇*<sup>h</sup>* being applied on the horizontal surfaces (constant *r*) and **k***∂<sup>r</sup>* acting along the vertical coordinate. *φ* is associated with pressure or geopotential, *b* is the buoyancy, *s* is the especific humidity for the atmosphere or the salinity for the ocean, **F** represents the external forcings and dissipation on **V**, *Q<sup>θ</sup>* the external forcings and dissipation on *θ* and *Qs* the ones on *s*. The terms **F**, *Q<sup>θ</sup>* and *Qs* are provided by physical parameterizations of the subgrid scale

<sup>−</sup>**<sup>V</sup>** · ∇*Rfixed*, *<sup>z</sup>* <sup>=</sup> *Rfixed*

<sup>−</sup>*H*, at sea bottom

*b* = *b*(*θ*,*s*,*r*) (37) *Dtθ* = *Q<sup>θ</sup>* (38) *Dts* = *Qs* (39)

*Dt* = *∂<sup>t</sup>* + **V** · ∇*A* + *r*˙*∂r*, (40)

∇ = ∇*<sup>h</sup>* + **k***∂r*, (41)

*Dtrs* <sup>−</sup> *Pr*, *<sup>r</sup>* <sup>=</sup> *Rs* (42)

0, at atmosphere top (43)

is the horizontal mean wind velocity. See figure 2 for illustration.

#### **2.5 Isomorphisms and physical coupling**

By taking a look at the previous subsections it is easy to recognize the change in variables for writing the isomorphism. Choosing *r* as the new isomorphic variable, such that (Marshal et al., 2004)


the boundary conditions (19,20,29,30) automatically remain isomorphic too. An illustration is depicted in Figure 3.

Then, the general equations are

$$D\_l \mathbf{V} + f \mathbf{k} \times \mathbf{V} + \nabla\_l \boldsymbol{\phi} = \mathbf{F} \tag{34}$$

$$
\partial\_r \phi - b = 0 \tag{35}
$$

$$\nabla\_r \cdot \mathbf{V} + \partial\_r \dot{r} = 0 \tag{36}$$

Fig. 3. Isomorphic relations and coupling. Source: (Marshal et al., 2004)

$$b = b(\theta, \mathbf{s}, r) \tag{37}$$

$$D\_t \theta = Q\_\theta \tag{38}$$

$$D\_l s = Q\_s \tag{39}$$

where *r* is the vertical, isomorphic coordinate. The total derivative takes the form

$$D\_t = \partial\_t + \mathbf{V} \cdot \nabla A + \dot{r} \partial\_{r\_\prime} \tag{40}$$

and

6 Numerical Simulations

Fig. 2. Mathematical description of the ocean boundary conditions. Source: (Marshal et al.,

*H* + *η*

By taking a look at the previous subsections it is easy to recognize the change in variables for writing the isomorphism. Choosing *r* as the new isomorphic variable, such that (Marshal et

> *ocean* ←→ *isomorphism* ←→ *atmosphere zr p ω r*˙ *ω*

*<sup>ρ</sup>*<sup>0</sup> *φ* Φ

�

*P* − *E Pr* 0

the boundary conditions (19,20,29,30) automatically remain isomorphic too. An illustration

*<sup>ρ</sup>*<sup>0</sup> *b* −*α θθ θ Ss q*

*<sup>s</sup> ps*

*Dt***V** + *f***k** × **V** + ∇*rφ* = **F** (34)

*∂rφ* − *b* = 0 (35) ∇*<sup>r</sup>* · **V** + *∂rr*˙ = 0 (36)

 *η* −*H*

*dz***V** (32)

(33)

<sup>&</sup>lt; **<sup>V</sup>** <sup>&</sup>gt;<sup>=</sup> <sup>1</sup>

is the horizontal mean wind velocity. See figure 2 for illustration.

*p*

<sup>−</sup>*<sup>g</sup> <sup>ρ</sup>*

*η r*

**2.5 Isomorphisms and physical coupling**

2004)

where

al., 2004)

is depicted in Figure 3.

Then, the general equations are

$$
\nabla = \nabla\_h + \mathbf{k} \partial\_{r\_l} \tag{41}
$$

with ∇*<sup>h</sup>* being applied on the horizontal surfaces (constant *r*) and **k***∂<sup>r</sup>* acting along the vertical coordinate. *φ* is associated with pressure or geopotential, *b* is the buoyancy, *s* is the especific humidity for the atmosphere or the salinity for the ocean, **F** represents the external forcings and dissipation on **V**, *Q<sup>θ</sup>* the external forcings and dissipation on *θ* and *Qs* the ones on *s*.

.

The terms **F**, *Q<sup>θ</sup>* and *Qs* are provided by physical parameterizations of the subgrid scale atmosphere and ocean fluxes.

The corresponding boundary conditions are

$$\dot{r} = \begin{cases} -\mathbf{V} \cdot \nabla R\_{f\text{ixed}} & z = R\_{f\text{ixed}} \\ \quad D\_l r\_s - P\_r & r = R\_s \end{cases} \tag{42}$$

where

$$R\_{fixed} = \begin{cases} -H\_{\prime} & \text{at sea bottom} \\ 0, & \text{at atmosphere top} \end{cases} \tag{43}$$

domain has been selected, the speed of the computations for all the methods is strongly dominated by the spatial resolution, even when in general the statistical methods will run many times faster than the dynamical ones. A coarse resolution implies less cells where to compute the models (physical, statistical or hybrid). The coarser the domain, the faster the computations, but the less the ability to accurately represent in a point the state of the Climate

<sup>457</sup> Climate System Simulations: An Integrated,

In this section we review some of these aspects and discuss how to construct a multi-scale

The planetary-scale, dynamical numerical models are called General Circulation Models (GCMs), and they can be coupled (CGCMs) or non coupled (e.g. atmospheric GCMs - AGCMs, and oceanic GCMs - OGCMs), with typical spatial resolutions in the order of 2.5*<sup>o</sup>* to 1*o*. The coupling makes reference to linking different models for mimetizing the various interactions of the Climate System components, following the ideas presented in section 2.1. Nonetheless, as it has been discussed previously, *S* can be simplified (see equation (2)) when considering certain temporal and spatial scales. Such cases are usually simulated by means of non-coupled GCMs, and due to their shorter execution times and some key advantages (see Mason (2008) and references therein), they have been the global dynamical models most

Coupled and non-coupled models may or may not coincide in their results, depending on several aspects. One of the main differences, as can be deduced from the previous discussion, is the inclusion or not of the various feedbacks. As an example, consider the simulation of the interactions between the atmosphere and the ocean on a seasonal (e.g. 3 months) scale. The coupled models are designed to guarantee the feedbacks, but the stand-alone atmospheric models consider the ocean effects via the boundary conditions (which in general evolve in time), without permiting any updates to the sea surface coming from the new atmospheric states. In other words, for these non-coupled atmospheric models to run it is first necessary to subscribe the evolution of the boundaries in a separate tier. This is why sometimes the stand-alone models are also called *two-tiered* models. One way to account for the differences among one-tiered and two-tiered simulations is to look at the energy balance between the surface and the top of the atmosphere. The radiative imbalance is usually less than 0.5 Wm−2, which is considered acceptable for seasonal timescales. However, the Climate System's slow component interactions can provide severe bias that must be considered in detail (Hazeleger

Amongst the main CGCMs used today are the CCSM-CESM (National Center for Atmospheric Research-NCAR Community Earth System Model (Blackmon et al., 2001)), COSMOS (Max Planck Institute-MPI Community Earth System Models (Roeckner et al., 2006)), HadGEM (Hadley Centre Global Environmental Model (Johns et al., 2006)) and CFS (NCEP Climate Forecast System (Saha et al., 2006)). Older and newer versions are available at

Both hindcast and forecast are important for climate simulations in research and decision-making. Retrospective simulations make use of grids filled with observations (e.g. gauge stations, satellite) for providing the GCMs' boundary conditions along the past period of interest while, obviously, there is no direct analog for a simulation of the future. The key idea here is that the GCM integrates exactly the same set of equations for the past or the future, the difference being on how the boundary evolves: in a hindcast, the observations are used for

approach for climate research and decision-making.

Multi-Scale Approach for Research and Decision-Making

**3.1 General Circulation Models: Coupled vs non-coupled**

a number of research and forecast centers around the world.

**3.2 Seasonal forecast methodologies for GCMs**

System.

widely utilised.

et al., 2010).

defines the location of the fixed boundary surface, while

$$R\_{\rm s} = R\_{\rm 0} + r\_{\rm s} \tag{44}$$

defines the location of the moving boundary surface, being

$$R\_0 = \begin{cases} 0, & \text{at sea surface} \\ R\_0(\mathbf{x}, y) = p\_\mathbf{s}(\mathbf{x}, y), & \text{at terrain surface} \end{cases} \tag{45}$$

The reference location of the moving boundary surface and its deviations are described by

$$r\_s = \begin{cases} \eta\_{\prime} & \text{at sea surface} \\ p\_{s\prime} & \text{at terrain surface} \end{cases} \tag{46}$$

Moreover,

$$P\_{\mathcal{I}} = \begin{cases} P - E\_{\prime} & \text{at sea surface} \\ 0, & \text{in the atmosphere} \end{cases} \tag{47}$$

is the volumetric flux through *rs*.

It is important to remark at this level how the system of equations presented are so complex that they need to be solved numerically in the general case. The last section shows a natural way to couple the needed quantities at the boundaries (see figure 3), which sketchs the coupling between ocean and atmosphere, the two most dynamic components of the Climate System. Bearing this review of the fundamental equations in mind, it will be easy to consider the particular expressions used in the different numerical models.

#### **3. An integrated, multi-scale approach**

The rationale behind climate modeling is to mathematically express the different interactions present in the Climate System. This is ussually accomplished by means of


The dynamical approach tries to describe the observed phenomenology using physics (Anderson, 2008), with sets of equations like the ones presented in the previous section. On the other hand, the statistical approach looks for the same goal, but it tries to recognize spatio-temporal patterns in the historical dataset of at least two different variables (usually known as the *predictor* and the *predictand*); it employs these patterns to build a mathematical model (Mason & Baddour, 2008), capable of reproducing the historical behaviour and to forecast the future, in analogous way to the dynamical approach. Finally, mixtures of both methodologies are also possible, entailing different kinds of hybrid models. The climate numerical simulations are the execution of these models. Executions reproducing the past are called *retrospective simulations* or *hindcasts*. Weather and seasonal forecasts, and climate change scenarios are examples of climate simulations looking at the future.

Each of the aforementioned approaches has its own pros and cons. For instance, the dynamical models need to solve the set of equations for each cell in the selected spatial domain, and thus they consume a lot of computational resources and are slow compared with the common statistical models. The latter, in turn, need homogeneous and quality datasets for the variables under study, an important issue in several regions of the planet. However, once the temporal 8 Numerical Simulations

0, at sea surface

The reference location of the moving boundary surface and its deviations are described by

*η*, at sea surface

*<sup>P</sup>* <sup>−</sup> *<sup>E</sup>*, at sea surface

It is important to remark at this level how the system of equations presented are so complex that they need to be solved numerically in the general case. The last section shows a natural way to couple the needed quantities at the boundaries (see figure 3), which sketchs the coupling between ocean and atmosphere, the two most dynamic components of the Climate System. Bearing this review of the fundamental equations in mind, it will be easy to consider

The rationale behind climate modeling is to mathematically express the different interactions

The dynamical approach tries to describe the observed phenomenology using physics (Anderson, 2008), with sets of equations like the ones presented in the previous section. On the other hand, the statistical approach looks for the same goal, but it tries to recognize spatio-temporal patterns in the historical dataset of at least two different variables (usually known as the *predictor* and the *predictand*); it employs these patterns to build a mathematical model (Mason & Baddour, 2008), capable of reproducing the historical behaviour and to forecast the future, in analogous way to the dynamical approach. Finally, mixtures of both methodologies are also possible, entailing different kinds of hybrid models. The climate numerical simulations are the execution of these models. Executions reproducing the past are called *retrospective simulations* or *hindcasts*. Weather and seasonal forecasts, and climate

Each of the aforementioned approaches has its own pros and cons. For instance, the dynamical models need to solve the set of equations for each cell in the selected spatial domain, and thus they consume a lot of computational resources and are slow compared with the common statistical models. The latter, in turn, need homogeneous and quality datasets for the variables under study, an important issue in several regions of the planet. However, once the temporal

*Rs* = *R*<sup>0</sup> + *rs* (44)

*ps*, at terrain surface (46)

0, in the atmosphere (47)

*<sup>R</sup>*0(*x*, *<sup>y</sup>*) = *ps*(*x*, *<sup>y</sup>*), at terrain surface (45)

defines the location of the fixed boundary surface, while

defines the location of the moving boundary surface, being

*rs* =

the particular expressions used in the different numerical models.

present in the Climate System. This is ussually accomplished by means of

change scenarios are examples of climate simulations looking at the future.

*Pr* =

*R*<sup>0</sup> =

Moreover,

is the volumetric flux through *rs*.

• Dynamical climate models • Statistical climate models. • Hybrid climate models.

**3. An integrated, multi-scale approach**

domain has been selected, the speed of the computations for all the methods is strongly dominated by the spatial resolution, even when in general the statistical methods will run many times faster than the dynamical ones. A coarse resolution implies less cells where to compute the models (physical, statistical or hybrid). The coarser the domain, the faster the computations, but the less the ability to accurately represent in a point the state of the Climate System.

In this section we review some of these aspects and discuss how to construct a multi-scale approach for climate research and decision-making.
