4.7. Surface Driven Convection

(in directory: verification/tutorial_deep_convection/)

deep convection setup

Figure 4.25 Schematic of simulation domain for the surface driven convection experiment. The domain is doubly periodic with an initially uniform temperature of 20 oC.

This experiment, Figure 4.25, showcasing MITgcm’s non-hydrostatic capability, was designed to explore the temporal and spatial characteristics of convection plumes as they might exist during a period of oceanic deep convection. It is

  • non-hydrostatic
  • doubly-periodic with cubic geometry
  • discretized with 50 m resolution in \(x, y, z\)
  • Cartesian
  • on an \(f\)-plane
  • using a linear equation of state

4.7.1. Overview

The model domain consists of an approximately 3 km square by 1 km deep box of initially unstratified, resting fluid. The domain is doubly periodic.

The experiment has 20 levels in the vertical, each of equal thickness \(\Delta z =\) 50 m (the horizontal resolution is also 50 m). The fluid is initially unstratified with a uniform reference potential temperature \(\theta =\) 20 oC. The equation of state used in this experiment is linear

(4.48)\[\rho = \rho_{0} ( 1 - \alpha_{\theta}\theta^{'} )\]

which is implemented in the model as a density anomaly equation

(4.49)\[\rho^{'} = -\rho_{0}\alpha_{\theta}\theta^{'}\]

with \(\rho_{0}=1000\,{\rm kg\,m}^{-3}\) and \(\alpha_{\theta}=2\times10^{-4}\,{\rm degrees}^{-1}\). Integrated forward in this configuration, the model state variable theta is equivalent to either in-situ temperature, \(T\), or potential temperature, \(\theta\). For consistency with other examples, in which the equation of state is non-linear, we use \(\theta\) to represent temperature here. This is the quantity that is carried in the model core equations.

As the fluid in the surface layer is cooled (at a mean rate of 800 Wm\(^2\)), it becomes convectively unstable and overturns, at first close to the grid-scale, but, as the flow matures, on larger scales (Figure 4.26 and Figure 4.27), under the influence of rotation (\(f_o = 10^{-4}\) s\(^{-1}\)).

vertical section deep cvct exp

Figure 4.26 Vertical section

surface section deep cvct exp

Figure 4.27 Surface section

Model parameters are specified in file input/data. The grid dimensions are prescribed in code/SIZE.h. The forcing (file input/Qsurf.bin) is specified in a binary data file generated using the Matlab script input/gendata.m.

4.7.2. Equations solved

The model is configured in non-hydrostatic form, that is, all terms in the Navier Stokes equations are retained and the pressure field is found, subject to appropriate boundary conditions, through inversion of a 3-D elliptic equation.

The implicit free surface form of the pressure equation described in Marshall et. al (1997) [marshall:97a] is employed. A horizontal Laplacian operator \(\nabla_{h}^2\) provides viscous dissipation. The thermodynamic forcing appears as a sink in the equation for potential temperature \(\theta\). This produces a set of equations solved in this configuration as follows:

\[\begin{split}\begin{aligned} \frac{Du}{Dt} - fv + \frac{1}{\rho}\frac{\partial p^{'}}{\partial x} - \nabla_{h}\cdot A_{h}\nabla_{h}u - \frac{\partial}{\partial z}A_{z}\frac{\partial u}{\partial z} & = & \begin{cases} 0 & \text{(surface)} \\ 0 & \text{(interior)} \end{cases} \\ \frac{Dv}{Dt} + fu + \frac{1}{\rho}\frac{\partial p^{'}}{\partial y} - \nabla_{h}\cdot A_{h}\nabla_{h}v - \frac{\partial}{\partial z}A_{z}\frac{\partial v}{\partial z} & = & \begin{cases} 0 & \text{(surface)} \\ 0 & \text{(interior)} \end{cases} \\ \frac{Dw}{Dt} + g \frac{\rho^{'}}{\rho} + \frac{1}{\rho}\frac{\partial p^{'}}{\partial z} - \nabla_{h}\cdot A_{h}\nabla_{h}w - \frac{\partial}{\partial z}A_{z}\frac{\partial w}{\partial z} & = & \begin{cases} 0 & \text{(surface)} \\ 0 & \text{(interior)} \end{cases} \\ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} + &=& 0 \\ \frac{D\theta}{Dt} - \nabla_{h}\cdot K_{h}\nabla_{h}\theta - \frac{\partial}{\partial z}K_{z}\frac{\partial\theta}{\partial z} & = & \begin{cases} {\cal F}_\theta & \text{(surface)} \\ 0 & \text{(interior)} \end{cases} \end{aligned}\end{split}\]

where \(u=\frac{Dx}{Dt}\), \(v=\frac{Dy}{Dt}\) and \(w=\frac{Dz}{Dt}\) are the components of the flow vector in directions \(x\), \(y\) and \(z\). The pressure is diagnosed through inversion (subject to appropriate boundary conditions) of a 3-D elliptic equation derived from the divergence of the momentum equations and continuity (see Section 1.3.6).

4.7.3. Discrete numerical configuration

The domain is discretized with a uniform grid spacing in each direction. There are 64 grid cells in directions \(x\) and \(y\) and 20 vertical levels thus the domain comprises a total of just over 80,000 gridpoints.

4.7.4. Numerical stability criteria and other considerations

For a heat flux of 800 Wm\(^2\) and a rotation rate of \(10^{-4}\) s\(^{-1}\) the plume-scale can be expected to be a few hundred meters guiding our choice of grid resolution. This in turn restricts the timestep we can take. It is also desirable to minimize the level of diffusion and viscosity we apply.

For this class of problem it is generally the advective time-scale which restricts the timestep.

For an extreme maximum flow speed of \(| \vec{u} | = 1 ms^{-1}\), at a resolution of 50 m, the implied maximum timestep for stability, \(\delta t_u\) is

\[\delta t_u = \frac{\Delta x}{| \vec{u} |} = 50 s\]

The choice of \(\delta t = 10\) s is a safe 20 percent of this maximum.

Interpreted in terms of a mixing-length hypothesis, a magnitude of Laplacian diffusion coefficient \(\kappa_h (= \kappa_v) = 0.1\) m\(^2\)s\(^{-1}\) is consistent with an eddy velocity of 2 mm s\(^{-1}\) correlated over 50 m.

4.7.5. Experiment configuration

The model configuration for this experiment resides under the directory verification/convection/. The experiment files

contain the code customizations and parameter settings for this experiment. Below we describe these experiment-specific customizations.

4.7.5.1. File code/CPP_OPTIONS.h

This file uses standard default values and does not contain customizations for this experiment.

4.7.5.2. File code/SIZE.h

Listing 4.22 verification/tutorial_deep_convection/code/SIZE.h
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
CBOP
C    !ROUTINE: SIZE.h
C    !INTERFACE:
C    include SIZE.h
C    !DESCRIPTION: \bv
C     *==========================================================*
C     | SIZE.h Declare size of underlying computational grid.
C     *==========================================================*
C     | The design here supports a three-dimensional model grid
C     | with indices I,J and K. The three-dimensional domain
C     | is comprised of nPx*nSx blocks (or tiles) of size sNx
C     | along the first (left-most index) axis, nPy*nSy blocks
C     | of size sNy along the second axis and one block of size
C     | Nr along the vertical (third) axis.
C     | Blocks/tiles have overlap regions of size OLx and OLy
C     | along the dimensions that are subdivided.
C     *==========================================================*
C     \ev
C
C     Voodoo numbers controlling data layout:
C     sNx :: Number of X points in tile.
C     sNy :: Number of Y points in tile.
C     OLx :: Tile overlap extent in X.
C     OLy :: Tile overlap extent in Y.
C     nSx :: Number of tiles per process in X.
C     nSy :: Number of tiles per process in Y.
C     nPx :: Number of processes to use in X.
C     nPy :: Number of processes to use in Y.
C     Nx  :: Number of points in X for the full domain.
C     Ny  :: Number of points in Y for the full domain.
C     Nr  :: Number of points in vertical direction.
CEOP
      INTEGER sNx
      INTEGER sNy
      INTEGER OLx
      INTEGER OLy
      INTEGER nSx
      INTEGER nSy
      INTEGER nPx
      INTEGER nPy
      INTEGER Nx
      INTEGER Ny
      INTEGER Nr
      PARAMETER (
     &           sNx =  50,
     &           sNy =  50,
     &           OLx =   2,
     &           OLy =   2,
     &           nSx =   2,
     &           nSy =   2,
     &           nPx =   1,
     &           nPy =   1,
     &           Nx  = sNx*nSx*nPx,
     &           Ny  = sNy*nSy*nPy,
     &           Nr  =  50)

C     MAX_OLX :: Set to the maximum overlap region size of any array
C     MAX_OLY    that will be exchanged. Controls the sizing of exch
C                routine buffers.
      INTEGER MAX_OLX
      INTEGER MAX_OLY
      PARAMETER ( MAX_OLX = OLx,
     &            MAX_OLY = OLy )

Three lines are customized in this file. These prescribe the domain grid dimensions.

  • Line 45,

    sNx=50,
    

    this line sets the lateral domain extent in grid points for the axis aligned with the \(x\)-coordinate.

  • Line 46,

    sNy=50,
    

    this line sets the lateral domain extent in grid points for the axis aligned with the \(y\)-coordinate.

  • Line 55,

    Nr=50,
    

    this line sets the vertical domain extent in grid points.

4.7.5.3. File input/data

Listing 4.23 verification/tutorial_deep_convection/input/data
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
 &PARM01
 tRef=20*20.,
 sRef=20*35.,
 viscAh=4.E-2,
 viscAz=4.E-2,
 no_slip_sides=.FALSE.,
 no_slip_bottom=.FALSE.,
 diffKhT=4.E-2,
 diffKzT=4.E-2,
 f0=1.E-4,
 beta=0.E-11,
 tAlpha=2.0E-4,
 sBeta =0.,
 gravity=10.,
 rhoConst=1000.,
 rhoNil=1000.,
 heatCapacity_Cp=4000.,
#rigidLid=.TRUE.,
 implicitFreeSurface=.TRUE.,
#exactConserv=.TRUE.,
 eosType='LINEAR',
 nonHydrostatic=.TRUE.,
 saltStepping=.FALSE.,
 &

# Elliptic solver parameters
 &PARM02
 cg2dMaxIters=1000,
 cg2dTargetResidual=1.E-9,
 cg3dMaxIters=100,
 cg3dTargetResidual=1.E-9,
 &

# Time stepping parameters
 &PARM03
 nIter0=0,
#endTime=43200.,
 nTimeSteps=3,
 deltaT=20.,
 abEps=0.1,
 pChkptFreq=43200.,
 chkptFreq=7200.,
 dumpFreq=1800.,
 monitorFreq=600.,
 monitorSelect=1,
 monitorFreq=1.,
 &

# Gridding parameters
 &PARM04
 usingCartesianGrid=.TRUE.,
 dXspacing=20.,
 dYspacing=20.,
 delZ=50*20.,
 &

# Input datasets
 &PARM05
 surfQnetFile=  'Qnet_p32.bin',
 hydrogThetaFile='T.120mn.bin',
 pSurfInitFile='Eta.120mn.bin',
 uVelInitFile =  'U.120mn.bin',
 vVelInitFile =  'V.120mn.bin',
 &

This file specifies the main parameters for the experiment. The parameters that are significant for this configuration are

  • Line 7,

    tRef=20*20.0,
    

    this line sets the initial and reference values of potential temperature at each model level in units of \(^{\circ}\mathrm{C}\). Here the value is arbitrary since, in this case, the flow evolves independently of the absolute magnitude of the reference temperature. For each depth level the initial and reference profiles will be uniform in \(x\) and \(y\).

  • Line 8,

    sRef=20*35.0,
    

    this line sets the initial and reference values of salinity at each model level in units of ppt. In this case salinity is set to an (arbitrary) uniform value of 35.0 ppt. However since, in this example, density is independent of salinity, an appropriately defined initial salinity could provide a useful passive tracer. For each depth level the initial and reference profiles will be uniform in \(x\) and \(y\).

  • Line 9,

    viscAh=0.1,
    

    this line sets the horizontal Laplacian dissipation coefficient to 0.1 \({\rm m^{2}s^{-1}}\). Boundary conditions for this operator are specified later.

  • Line 10,

    viscAz=0.1,
    

    this line sets the vertical Laplacian frictional dissipation coefficient to 0.1 \({\rm m^{2}s^{-1}}\). Boundary conditions for this operator are specified later.

  • Line 11,

    no_slip_sides=.FALSE.
    

    this line selects a free-slip lateral boundary condition for the horizontal Laplacian friction operator e.g. \(\frac{\partial u}{\partial y}\)=0 along boundaries in \(y\) and \(\frac{\partial v}{\partial x}\)=0 along boundaries in \(x\).

  • Lines 12,

    no_slip_bottom=.TRUE.
    

    this line selects a no-slip boundary condition for the bottom boundary condition in the vertical Laplacian friction operator e.g., \(u=v=0\) at \(z=-H\), where \(H\) is the local depth of the domain.

  • Line 13,

    diffKhT=0.1,
    

    this line sets the horizontal diffusion coefficient for temperature to 0.1 \(\rm m^{2}s^{-1}\). The boundary condition on this operator is \(\frac{\partial}{\partial x}=\frac{\partial}{\partial y}=0\) at all boundaries.

  • Line 14,

    diffKzT=0.1,
    

    this line sets the vertical diffusion coefficient for temperature to 0.1 \({\rm m^{2}s^{-1}}\). The boundary condition on this operator is \(\frac{\partial}{\partial z}\) = 0 on all boundaries.

  • Line 15,

    f0=1E-4,
    

    this line sets the Coriolis parameter to \(1 \times 10^{-4}\) s\(^{-1}\). Since \(\beta = 0.0\) this value is then adopted throughout the domain.

  • Line 16,

    beta=0.E-11,
    

    this line sets the the variation of Coriolis parameter with latitude to \(0\).

  • Line 17,

    tAlpha=2.E-4,
    

    This line sets the thermal expansion coefficient for the fluid to \(2 \times 10^{-4}\) oC\(^{-1}\).

  • Line 18,

    sBeta=0,
    

    This line sets the saline expansion coefficient for the fluid to \(0\), consistent with salt’s passive role in this example.

  • Line 23-24,

    rigidLid=.FALSE.,
    implicitFreeSurface=.TRUE.,
    

    Selects the barotropic pressure equation to be the implicit free surface formulation.

  • Line 26,

    eosType='LINEAR',
    

    Selects the linear form of the equation of state.

  • Line 27,

    nonHydrostatic=.TRUE.,
    

    Selects for non-hydrostatic code.

  • Line 33,

    cg2dMaxIters=1000,
    

    Inactive - the pressure field in a non-hydrostatic simulation is inverted through a 3-D elliptic equation.

  • Line 34,

    cg2dTargetResidual=1.E-9,
    

    Inactive - the pressure field in a non-hydrostatic simulation is inverted through a 3-D elliptic equation.

  • Line 35,

    cg3dMaxIters=40,
    

    This line sets the maximum number of iterations the 3-D conjugate gradient solver will use to 40, irrespective of the convergence criteria being met.

  • Line 36,

    cg3dTargetResidual=1.E-9,
    

    Sets the tolerance which the 3-D conjugate gradient solver will use to test for convergence in equation (2.60) to \(1 \times 10^{-9}\). The solver will iterate until the tolerance falls below this value or until the maximum number of solver iterations is reached.

  • Line43,

    nTimeSteps=8640.,
    

    Sets the number of timesteps at which this simulation will terminate (in this case 8640 timesteps or 1 day or \(\delta t = 10\) s). At the end of a simulation a checkpoint file is automatically written so that a numerical experiment can consist of multiple stages.

  • Line 44,

    deltaT=10,
    

    Sets the timestep \(\delta t\) to 10 s.

  • Line 57,

    dXspacing=50.0,
    

    Sets horizontal (\(x\)-direction) grid interval to 50 m.

  • Line 58,

    dYspacing=50.0,
    

    Sets horizontal (\(y\)-direction) grid interval to 50 m.

  • Line 59,

    delZ=20*50.0,
    

    Sets vertical grid spacing to 50 m. Must be consistent with code/SIZE.h. Here, 20 corresponds to the number of vertical levels.

  • Line64,

    surfQfile='Qsurf.bin'
    

    This line specifies the name of the file from which the surface heat flux is read. This file is a 2-D (\(x,y\)) map. It is assumed to contain 64-bit binary numbers giving the value of \(Q\) (W m\(^2\)) to be applied in each surface grid cell, ordered with the \(x\) coordinate varying fastest. The points are ordered from low coordinate to high coordinate for both axes. The matlab program input/gendata.m shows how to generate the surface heat flux file used in the example.

4.7.5.4. File input/data.pkg

This file uses standard default values and does not contain customizations for this experiment.

4.7.5.5. File input/eedata

This file uses standard default values and does not contain customizations for this experiment.

4.7.5.6. File input/Qsurf.bin

The file input/Qsurf.bin specifies a 2-D (\(x,y\)) map of heat flux values where \(Q = Q_o \times ( 0.5 + \mbox{random number between 0 and 1})\).

In the example \(Q_o = 800\) W m\(^{-2}\) so that values of \(Q\) lie in the range 400 to 1200 W m\(^{-2}\) with a mean of \(Q_o\). Although the flux models a loss, because it is directed upwards, according to the model’s sign convention, \(Q\) is positive.