Numerical modeling of Atmosphere and Ocean

[vinu]

Many of us are familiar in using Re-analysis data sets from
various sources. A few may be aware of how these are modelled
in real time. A lot of effort has been made since 1940 (the early
stages of computer) to model the atmospher and ocean. Since then,
thoughts about 'Numerical Schemes with stable solutions'
were also involved.

Here you can download a technical report by Global
Atmospheric Research Programme- (GARP) published in early 1970's.
Any of todays models can not exclude these basic laws of modeling
structured by Messinger and Arakawa (1970's). Most of us may be aware
of NCEP/NCAR, ECMWF, GFDL models and many others. All these models
are designed according to these 'basic' rules
developed in early 1970's.


Numerical Methods used in Atmospheric Models (Mesinger and Arakawa,1976),
GARP publication series, No.17. (69 pages, 14 MB)

PS: If the file or link does not work properly, please feed back.

----------------------------------------------------------------------------------------------------

[vinu]

This is a short note on Numerical Ocean Model developement.
A Good and Short note to understand the basics about
Ocean Modeling like,

What is an Ocean Model
What are the governing equations of a simple ocean model
What are the various ocean models usually in practice
What is a grid and how they are designed
How to arrange these governing equations in a grid
How to integrate it in time for a forecast.

read more in Ocean Model (PDF file ) -74 pages (1.4 MB)

Source: http://csep1.phy.ornl.gov/om/om.html
Copyright (C) 1991, 1992, 1993, 1994, 1995 by the
Computational Science Education Project

-----------------------------------------------------------------------------------------------------

[Suprit]

I have browsed through the modelling book and its informative indeed, thanks.

A valuable book for the fellows interesting in modelling ( although not available free on the net and not so popular in oceanography student circles but certainly a must for mets ) is "Numerical prediction and dynamic meteorology" by Haltiner & Williams published by wiley. The concepts are dealt from the basics, and in straightforward manner.

[vinu]

Books which are worth refering to understand
about the basics of Ocean models-

Numerical Models of Oceans and Oceanic Processes
by L.H. Kantha, C.A. Clayson, International Geophysics Series,
Vol.66, 2000.
(-This book combines the fundementals of GFD and
a good start point for a beginner. It include examples for simple
models and expanded it in the finite difference form which will be
useful in understanding the development of a model-)

Numerical Modeling of Ocean Dynamics
by Z. Kowalik, T. S. Murty, Advanced Series on Ocean Engineering,
World Scientific, 1993
(-This book points to the essentials of advanced Numerical
methods in Ocean Modeling-).

Numerical Ocean Circulation Modeling
by Dale B. Haidvogel, Aike Beckmann, Series on Environmental
Science and Management. Vol.2, Imperial College Press, 1999.
(-This is a short book wich introduces Numerical methods in
Ocean Modeling and discuss the various OGCM usually under practice-).


Small Scale Processess in Geophysical Fluid Flows
by L. H. Kantha, C. A. Clayson, International Geophysics Series,
Vol.67, 2000.
(-This book is usefull to understand more about upper layer
(boundary layer) processes. It encorporate the modeling of Chemical
and Biological component in the Ocean-)

[vinu]

Here in this post, let us see a simple Hydro-dynamical Ocean Model.

What is a model?

Roughly speaking, it is a 'structure' to represent something. For example, a model of a building which is made up of 'card-board'- is a physical model. Now in our discipline, an ocean process (like its motion, heating, mixing etc.) can be modeled using some governing equations', under which these processes are happening. For example, Miss. A has 'X' income and 'Y' expenditure. Then her savings is, 'S = X-Y'. This is a simple 'model' for 'savings of Miss. A'. Likewise, every process in the ocean is also under some basic 'physical balance' and it can be reasonably represented using some sets of equations. [b]This is our key.

[b]Here let us see a simple hydrodynamical ocean model which is usually called, under practice, as 1&1/2 Layer model .

What is a 1&1/2 layer model?

Here the ocean is assumed into two layer- an upper layer and lower layer with slightly different density. However only the upper layer is 'under motion'. The 'lower layer' is very deep so that it is assumed to be at 'rest'. This is called 1&1/2 layer (means 'one' active layer and the rest 'quiescent layer'). These layers represent the 'well-mixed' upper layer of the real ocean and relatively slow and inactive 'deep-ocean'. This means that 1&1/2 layer version of ocean model can well represent a real ocean and that is its importance. (An upper layer thickness of 100m is a good practice).

The figure represents the structure of a 1&1/2 layer model with some necessary dimensional representations.



The upper layer density is 'rho-1' and lower layer density is'rho-2' the upper layer velocities are 'U'(zonal) and 'V'(meridional) and lower layer velocities are zero . The 'free-surface' elevation is 'eta' and 'inter-facial' elevation is 'h'. The upper layer thickness is 'H' and lower layer thickness is 'Hd'.

If the lower layer is under 'rest', why is it important? Why the entire bottom layer is not taken as an 'ocean floor?

The lower layer interacts with the upper layer through the density stratification at the interface. Eventhough the lower layer is under rest, the difference in density at the interface causes internal waves to develop. And thus it affects the 'upper-layer'.

What are the governing equations in such situation?

The main process to model is its motion. Let us use the shallow-water
wave equation for the upper layer


The equation (1) and (2) are the 'momentum equations'. The L.H.S represents the 'evolution of velocity ' in time. The first term in the R.H.S represents the 'coriolis acceleration'. The second term in the R.H.S is the 'pressure-gradient' term. Equation (3) is the 'continuity equation', which means that the 'evolution' of the layer thickness (L.H.S) is balanced by the 'horizontal divergence' (R.H.S) of the layer.

All the terms in the equation (1) and (2) are known for the upper layer. However in (3) we need 'h', which is the 'inter-facial' displacement due to the 'density stratification'.

How can we model 'h' (the inter-facial displacement)?

Since the lower layer is in rest, it implies that, there is no pressure gradient between two points in the 'lower layer' in order to drive the motion. Let us use this tool and try to find 'h' consider the pressure at 'A' (shown in Figure-1). Using the simple hydrostatic balance (p=h.rho.g), the pressure at A is the sum of the pressure of the water column above it. From a careful examination of the figure we can estimate it as shown in Figure-2 as 'Pa'. Now take another point Pb on the same level, and subtract Pa-Pb. According to above constraint (i.e. pressure gradient is zero) Pa-Pb=0. This will lead us to get the equation (4).
Now we can find all the terms in the R.H.S of equation 1 to 3.

Now look at the equation (4). Usually in the ocean, the lower layer is heavier(denser) than the upper layer, resulting a negative value of (rho-1 - rho-2)in equation (4). This implies that, the 'free surface' elevation (eta) is in opposite 'direction' with the interface-displacement (h).

How do these equations work (or how does it model the ocean)?

As mentioned earlier, the savings of Miss. A is 'S = X-Y'. Well, this is a useless equation, because it doesnt tells us what is the savings of Miss. A last month, or next month. Suppose we have the 'rate-of-savings' of Miss. A, (something like, dS/dt = (X-Y), where 't' is the time.) Let us 'integrate this equation', with the knowledge of todays 'X and Y'. PS: S(t2) denotes S_suffix_t2, means value of S at time t2 i.e. (S(t2) - S(t1) )/(t2-t1) = X-Y (in finite difference form).

i.e. S (t2) = S (t1) + (X-Y)*(t2-t1) , where 't1'= today, 't2' = tomorrow.

Now we end-up in S (t2) , that is Miss. A's tomorrows of savings,(provided todays rate is 'dS/dt'). This is how 'a forecast' is being done. The L.H.S of equation 1 to 3 is a 'time-derivative'. We know that if we 'integrate' a time derivative in the present time, we will get values for the future-time' (depending upon the time-interval we used in the integration).

So what?

Now we have to think about the way-of-integration . At first we need a 'space' to do the integration. Simply, we need to design our ocean in a regular 'grid' which represents some 'points (location)' in the real ocean. And then we have to express our equations in a 'finite difference form' for integrating'numerically'. (Just like what we have seen in Miss. A's example). (A good note about stable numerical schemes is here).

What is a Grid?

Roughly speaking, it is a group of points representing the ocean. We use these points to 'integrate' the above sets of equations. We can not take as many points as we want, because it needs laborious computational power. So let us restrict our 'points' necessary enough to represent the physical processes of the ocean in which we are interested in. Likewise we can-not restrict our points 'to-a-very-minimum' number to save the computational time because we may lose useful physical details of the ocean. There are well defined criteria, though they are not the concern of this post.

What type of grid we use?

In equation 1 to 3, we have three 'prognostic variables', i.e. variables to find out for 'future time'. They are 'eta or h, U and V. There are certain rules to define these variables in certain points of a grid in the ocean. Depending on various architecture, there are a wide variety of grids that are in practice. Let us see three of them.

A-Grid.


It is a simple grid. The variable 'eta or h, U, V' are defined at the 'same location' in the ocean. (see figure above)

B-Grid.


It is a 'staggered grid'. The 'U, V locations are at the corner points of the grid, whereas 'eta or h' is defined at the center of the grid.(see figure above)

C-Grid


It is also a 'staggered grid'. The 'U' locations are 'half-grid size away' to the 'west of 'eta'. The 'V' locations are 'half-grid size away' to the south of 'eta'. (See figure above).

Why these jargons?

The staggering of the grids is related with the 'computational stability of various numerical schemes'. For designing appropriate finite difference scheme in space we need to 'locate' these variables in suitable positions so that the resulting scheme yields maximum stability in calculation. A note about this is included here.



Which grid is appropriate?

This depends upon the problem we are interested in. For the above set of 3 equations, any of the above grid can be used. Each has its own advantages and disadvantages. For example, if we want to solve 'U' component of velocity at a 'grid-location', we need 'dP/dX and 'f.V'. If we set a C-grid (suppose), we have two P (or eta) points on the right and left side of 'U'. So it is easy to take dP/dx in this case. However, we have no 'V' values at the same position of 'U' (see figure C-grid for understanding the situation). So we have to average the 'V' ('f' as well) from the neighboring corners to find the 'f.V'. This is not precise, but theres no other choice. This example tells the various technical details of a particular grid. If we use A-grid instead, all the values are available at the same location. However the resulting finite difference scheme is not necessarily computationally stable. Concluding these technical side (not discussing much in this post), let us set for a C-grid which is more appropriate.

Now, we have a Grid-, then?

In the Miss.A's saving example, we have seen a step of 'integration'. If we use the time-step of integration as 'one-month', we will get the savings of Miss. A for one month at once. (In a single step calculation). However, from the common logic we can follow that it is not-accurate, because the saving rate 'dS/dt ' may be changing from 'day to day'. Likewise in numerical modeling of geophysical flows, the selection of 'time-step' for the integration needs much more attention.

CFL limit

In numerical integration of Geophysical Fluids, it is necessary to follow a constraint in the time-step of the integration. The C.F.L condition says that (Courant-Friedrichs-Lewy condition) the time step chosen, dt < dx/C, where 'dx' is the grid size and 'C' is the phase speed of 'fastest wave' in the domain. In Geophysical Fluid flow, the fastest wave is the 'gravity wave' with speed = square-root(g.H), where g=acceleration due to gravity and 'H' is the fluid depth. For a 'barotropic' mode (the entire fluid column is assumed to be in uniform speed) H is usually 4000m. This yields the 'gravity wave' speed c = 200m/sec, which is very large. In this case for example, for a grid-size that is of 1 degree, the CFL limit says dt < (dx/C) , which comes around 600 seconds. So this is the time-step we can grant for our integration. We have to forecast all the variable in a short-future of 600 sec (or less). From there we do the next integration to get the next future values and so on.

However in 1&1/2 Layer version of ocean models, there is no 'barotropic' mode because we are concerned only with the motion of the'upper-layer' (lower layer is in rest). As well as the depth of the fluid which will be in some order of 100's, so the time-step can be effectively increased.

What is next?

Now we have a set of equations, we have grid, we have finite-difference form, we have chosen a 'time-step' for integration. Now comes the initial-condition as well as the boundary-conditions. In simple words, initial conditions means, the state of the ocean we refer to the model initially. The consequent solutions in time (or evolution in time) will be depending on this initial state of the ocean. Thus an accurate initial condition is necessary to yield a reliable evolution of the model in time.

For example, we can assume, the ocean has a particular Temperature and Salinity values initially. This may be selected from the available climatological values (Data). For the momentum equations, we can assume, for example, Ocean is initially at rest. That means zero velocity. In our above 1&1/2 layer, we are concerned only with the momentum equations. Sea-Surface elevation (eta) can be assumed to be flat initially so that there is no pressure-gradient across the surface and hence no-motion.

Now let us see what is meant by a boundary condition. The boundary means the portion, where the model is separated from its surroundings. Simply speaking, we can expect boundaries at the surface, bottom and the lateral walls. Sometimes the lateral walls may be open to some other ocean basin. From our simple logic we can expect that the water is not allowed (or should not) to flow out of these boundaries (except for open boundaries where the ocean is connected to another ocean basins or channels).

Thus generally we can classify the boundary condition as Closed boundary conditions and Open boundary conditions. For example the lateral wall as well as the bottom of the sea-floor is the closed boundaries whereas the connection of a semi-enclosed bay to the open ocean is an open boundary. At the closed boundary (i.e. for lateral walls) the velocity that is normal to the wall is assumed as zero (flow cannot penetrate the wall). The velocity parallel to the wall can be either set to non-zero or set to zero. In the non-zero case it is called 'free-slip' boundary condition in-which the water does not feel any frictional force of the wall. The zero-velocity case is called as'no-slip means that the friction within the grid close to the lateral wall is important. These two are ideal cases whereas a good practice for general exercises in ocean modeling. (The description about open boundaries are not mentioning here).


(...Continued)

If any of the above statements found misleading PLEASE FEED BACK.
- - - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - -

[vinu]

The following is Praveen's commend over the above post, I quote it for him here,
-----------------------------------------------------------------------------------------------
under C grid picture, it's description

there "east" of eta is making problem, is it "west" ??

is it like this?
The 'U' locations are 'half grid size away'


g'=g.(rho2-rho1)/rho2

what is g', is it reduced gravity something like that?

the thread is much informative. try to include different cordinates
systems like sigma cordinates,......

with regards
-------------------------------------------------------------------------------------------------

Dear Praveen

Thanks for the commend and type mistake you noticed in the
above post. I have corrected those.

yes, the term g' is the 'reduced gravity'.
However, for the 1&1/2 layer version, g'
need not necessarily find during the calculation, because if you
substitute g' in 'eta' equation, you will reach equation (4) by
eliminating the 'reduced gravity'.

Now, If we apply a Rigid-Lid (i.e. ocean surface is believed to as rest,
or just like a solid slab, then we can eliminate 'eta' from the momentum
equation, and can be replaced by the 'interfacial displacement h'
from the equation (4). As well as, the 'eta' term in equation (3) can
be neglected comparing with the 'interfacial displacement 'h' '. Then
these equations will give the motion under the 'reduced gravity' .Actualy
this is what being done in a 1&1/2 layer model. This way, we can avoid the
fastest moving gravity wave with speed sqrt(g.H) and reduces it to sqrt(g'.H).
This will provide an effective way to increase the time-step for
integration as with the C.F.L limit.


thanks for your feed back

PS: this time i quote your commend, you can post your replies and notices
in the thread itslef. we meant for this 'Forum Structure' to open
discussions, If in future if the postings are overflowing, let us merge it
with appropriate ones.

[vinu]

Shall we begin the modeling?

Thus we have looked at the necessary aspects to begin the simple Numerical Ocean Model. Now a question remains, well, we have a set of equation in 1 to 3. And we mentioned the initial model velocity as zero and sea-surface is assumed to be flat so that there is no pressure gradient either. In such a situation, how does the Ocean begin to move?. Thus it is clear that somehow the motion should initiate. It is a good practice to give some particular 'sea-surface elevation' as the initial condition. This will start the motion in our Ocean Model. However there are other important 'forcing' in the atmosphere which immediately impart momentum to the Ocean model (perhaps most of the large scale ocean circulation are driven by this 'forcing', and that is the 'wind-forcing'. Thus surface wind stress can be added to the momentum equations (1) and (2).

How does these wind go into the ocean model?

According to the Navier-Stokes equation, the stress component in the momentum equation is represented as X/rho, where X is the body force. Here in the real Ocean case the body force is the wind stress and can be added to the equation 1 and 2 as Tx/rho and Ty/rho, where 'Tx and 'Ty' are wind stress in 'x' (zonal) and 'y'(meridional) direction.

'To what depth is this wind stress going to impart its momentum to the ocean?.

Well, we know that, wind stress causes the 'Ekman spiral' in the ocean applicable to an 'Ekman depth'. In our simple 1&1/2 layer version we have only one active layer which has a 'unique' velocity. So we are going to assume that the surface wind is going to affect the entire layer. This can be simply derived using the following concept. In the 'active layer' the velocity is unique. Thus if we integrate the equation 1 and 2 vertically over the depth of the layer 'H', we getH.du/dt = H.(fv) - H.dP/dx + (1/rho).(Ts-Tb), where Ts is the surface stress (of course due to the wind) andTb is the bottom stress- due to the bottom layer friction.

However we can safely neglect the term Tb, so that the above equation becomes du/dt = f.v - dP/dx + (1/rho).Ts/H. So the wind stress is ordered as (1/rho).Ts/H, where 'H' is the depth of the 'upper-layer'.

The above all points are some necessary notes for modeling the Ocean circulation. However, lets consider a situation. Suppose the ocean is forced with a wind stress for some period of time and suddenly the wind is switched off. Now the ocean will reach stability according to equation 1 to 3. However, from the common knowledge we know that after a long time the system eventually comes back to rest (in the real ocean). Why? . It is because the ocean has friction. So we have to account this amount of friction in the ocean so that it accounts the 'dissipation of the motion'. Of course the ocean water has molecular viscosity' as every other fluid. However the scale of such viscous force (frictional force) is very few and not appropriate for large scale motion. However we have the enlarged version of this friction for the ocean, and we refer it as 'Diffusion' or 'eddy diffusion'. This term is represented as Am.d^2u/dx^2., where Am. is called the 'eddy viscos
ity coefficient'. Perhaps this kind of frictional term is necessary for computational stability.

So finally we have to add two more term to the R.H.S of equation 1 and 2.

They are + Tx/(rho.H) + Am.d^2u/dx^2 . Reference

The pioneer works in ocean modeling started from 60's and 70's. Most of the fundamental works in reduced gravity layered model are initiated by several scientist. Among them, the team led by J.J. O'Brien is outstanding. He and his students made many contributions to the development of simple layered models to the state-of-the-art Ocean General Circulation models (OGCM). A number of classical papers can be obtained from here (scanned copy)
(J.J. O'Brien.



Jensen T. G., 1991: Modeling of the Seasonal Undercurrents in the Somali currentSystem, J. Geophysical Research, Vol. 96, No. C12, 22,151-22,167.

Jensen T. G., 1996: Artificial retardation of barotropic waves in layered ocean models. Monthly Weather Review, 124, 1272-1283.

Busalacchi, A. J., and J. J. O'Brien, 1980: The seasonal variability in a model of the tropical Pacific. J. Physical Oceanography, 10, 1929-1951.

McCreary, J. P., and P. K. Kundu, 1988: Numerical investigation of Somali current during southwest monsoon. J. Marine Research, 46, 25-58.

McCreary, J. P., P. K. Kundu and R. L. Molinary, 1993: A numerical investigationof dynamics and thermodynamics and mixed layer process in Indian Ocean, Progress in Oceanography, 31, 181-244




PS: If any of the above snetence found miss-leading PLEASE FEED BACK,
Excuse is requested for type-mistakes.
-------------------------------------------------------------------------------------------------------

[praveen]

hi vinu
The above thread is becoming more attractive. I have a doubt about open boundaries in models. If we set an open boundary in a semi enclosed region, like if we model Indian ocean basin. we have to set southern part as open boundary. then how the system become stable. Do we need any aditional information in that region?? . In some model I saw inclusion of a sponge in that region. Can you explain this in your thread?

[vinu]

Hi Praveen,
Yes, open boundaries need some more explanation. I hope to add it
with the above thread as time comes.

There are standard methods to deal with 'Open boundaries', In short,
some of them can be listed as (1) Dirichlet B.C. (2) Neumann B.C. (3)
Newtonian damping (4) Sponges (5) Radiation B.C (6) Cyclic B.C etc...

In the Dirichlet, prescribed values are specified at the boundary (can be
from observational data or some other model) and 'clamped' to the
model at each time-step.

In Neumann, it assumes that, there is 'no gradient' of flux across the
boundary, i.e. the flux are freely allowed to flow out/in to the domain.

In Newtonian damping, the variable (eg. Temp or Salinity )are damped to
a reference value at the open boundaries. This is like a 'relaxation to the
pre-defined values'.

Sponges are 'like a cushion' at the boundaries. This is accomplished by
making the fluid 'high viscous' (like oil) (increasing the values of Am) at
the boundaries so that high fluctiations 'impinging' on the boundary is damped quickly.

Radiation boundaries, as the name implies, it allows the 'waves to radiate' out of the domain.

Cyclic boundary let the waves exiting from one side of the domain to
enter through the another side. For example,the waves leaving the western
boundary will enter into the model domain through the easter boundary.
This is accompanied by restoring values of the 'boundaries' to the values
of its opposite side.

well, this is a short discription, each of the above method comes with
advantage and penalty.

The first book listed here explains these in detail.
------------------------------------------------------------------------------------------------------

[vinu]

Detailed Physical and Numerical aspects of the Ocean and Atmospheric
Modeling is described in the following text book

Fundamentals of Ocean Climate Models,

by Stephen Griffies, GFDL, Princeton University.

Here is a work book for Ocean modeling practices,

by Zygmunt Kowalik, which comes with necessary FORTRAN subroutines.

PS: These two books are contributed to our web-site by our member Ramon Garcia.
Thank you Ramon, your contribution is highly appreciated.


Thank you.
------------------------------------------------------------------------------------------------