Multiscale two-way embedding schemes for free-surface primitive equations in the “Multidisciplinary Simulation, Estimation and Assimilation System”

Abstract

We derive conservative time-dependent structured discretizations and two-way embedded (nested) schemes for multiscale ocean dynamics governed by primitive equations (PEs) with a nonlinear free surface. Our multiscale goal is to resolve tidal-to-mesoscale processes and interactions over large multiresolution telescoping domains with complex geometries including shallow seas with strong tides, steep shelfbreaks, and deep ocean interactions. We first provide an implicit time-stepping algorithm for the nonlinear free-surface PEs and then derive a consistent time-dependent spatial discretization with a generalized vertical grid. This leads to a novel time-dependent finite volume formulation for structured grids on spherical or Cartesian coordinates, second order in time and space, which preserves mass and tracers in the presence of a time-varying free surface. We then introduce the concept of two-way nesting, implicit in space and time, which exchanges all of the updated fields values across grids, as soon as they become available. A class of such powerful nesting schemes applicable to telescoping grids of PE models with a nonlinear free surface is derived. The schemes mainly differ in the fine-to-coarse scale transfers and in the interpolations and numerical filtering, specifically for the barotropic velocity and surface pressure components of the two-way exchanges. Our scheme comparisons show that for nesting with free surfaces, the most accurate scheme has the strongest implicit couplings among grids. We complete a theoretical truncation error analysis to confirm and mathematically explain findings. Results of our discretizations and two-way nesting are presented in realistic multiscale simulations with data assimilation for the middle Atlantic Bight shelfbreak region off the east coast of the USA, the Philippine archipelago, and the Taiwan–Kuroshio region. Multiscale modeling with two-way nesting enables an easy use of different sub-gridscale parameterizations in each nested domain. The new developments drastically enhance the predictive capability and robustness of our predictions, both qualitatively and quantitatively. Without them, our multiscale multiprocess simulations either were not possible or did not match ocean data.

Appendices

Appendix 1: Additional details of the conservative discrete equations

Our horizontal grids are structured Cartesian or spherical grids, either of which can be rotated from the standard geographic orientation (for more details, see Haley 1999, Section 2.2). To obtain the spherical representation of our Eqs. 37–44 and “Appendix  1.2”, identify x with longitude and y with latitude. Then multiply all Δx terms by the radius of the earth (R _earth) times the cosine of the latitude and all Δy terms by R _earth. Our model also includes options for atmospheric forcing, assimilation, tidal forcing, and a river forcing that employs relaxation time constants which can be tuned to reproduce the desired mass, salt, and internal energy (heat) transports. All of these options are compatible with the two-way nesting scheme.

1.1 Appendix 1.1: Vertical grid

In Section 2.4, we introduced our vertical discretization, defining first a set of terrain-following depths for the undisturbed mean sea level, $z^{\rm MSL}_{i,j,k}$. Here we present the details of $z^{\rm MSL}_{i,j,k}$. We can currently employ five different schemes for defining these vertical levels, two of which are new:

• (a) σ-Coordinates (Phillips 1957)

  $$ \label{Eq:SigmaCoords} z^{\rm MSL}_{i,j,k} = - \sigma_k H_{i,j} $$

  (56)

  where 0 ≤ σ _k  ≤ 1

• (b) Hybrid coordinates (Spall and Robinson 1989)

  $$ \label{Eq:HybridCoords} z^{\rm MSL}_{i,j,k} = \left\{ \begin{array}{lcc} \tilde{z}_k & {\rm if} & k \le k_c \\ -h_c - \sigma_k \left( H_{i,j} - h_c\right) & {\rm if} & k > k_c \end{array} \right. $$

  (57)

  where $\tilde{z}_k$ are a set of constant depths and h _c is the sum of the top k _c flat level depths

• (c) Double σ-coordinates (Lozano et al. 1994)

  $$ \label{Eq:DoubleSigmaCoords} z^{\rm MSL}_{i,j,k} = \left\{ \begin{array}{lcc} - \sigma_k \tilde{f}_{i,j} & {\rm if} & k \le k_c \\ -\tilde{f}_{i,j} - (\sigma_k - 1) \left( H_{i,j} - \tilde{f}_{i,j}\right) & {\rm if} & k > k_c \end{array} \right. $$

  (58)
  $$\begin{array}{rll} \tilde{f}_{i,j} &=& \frac{z_{c_1} + z_{c_2}}{2} + \frac{z_{c_2} - z_{c_1}}{2} \\&&\times \tanh\left[\frac{2 \alpha}{z_{c_2} - z_{c_1}} \left(H_{i,j} - h_{ref} \right)\right] \\ \sigma_k &\in& \left\{ \begin{array}{lcc} [0, 1] & {\rm if} & k \le k_c \\ {[1, 2]} & {\rm if} & k > k_c \end{array} \right. \end{array} $$

  where $\tilde{f}_{i,j}$ is the variable interface depth between the upper and lower σ-systems, $z_{c_1}$ and $z_{c_2}$ are the shallow and deep bounds for $\tilde{f}_{i,j}$, h _ref is the reference topographic depth at which the hyperbolic tangent term changes sign, and α is a nondimensional slope parameter ($|| \nabla \tilde{f} || \le \alpha || \nabla H ||$).

• (d) Multi-σ-coordinates This new system is a generalization of the double σ system in which, for Pσ-systems, we define P + 1 nonintersecting interface surfaces. Then the depths are found from

  $$ \begin{array}{rll}\label{Eq:MultiSigmaCoords} z^{\rm MSL}_{i,j,k} &=& -\tilde{f}^{p-1}_{i,j} - (\sigma_k - p + 1) \tilde{f}^p_{i,j} \quad \\&&\quad {\rm for}\quad k_{p-1} < k \le k_p \end{array} $$

  (59)
  $$ \begin{array}{lll}&&\tilde{f}^0_{i,j} = 0 ; \quad \tilde{f}^P_{i,j} = H_{i,j} \\ &&\sigma_k \in {[(p-1),\ p]} \quad {\rm for} \quad k_{p-1} < k \le k_p \end{array} $$

  The intermediate interfaces are free to be chosen from many criteria, including key σ _θ surfaces (e.g., top of mean thermocline) or large mean vertical gradients.

• (e) General coordinates For this new system, we provide a 3D field of level thicknesses, $\Delta z^{\rm MSL}_{i,j,k}$, under the constraint

  $$ \sum\limits^K_{k=1} \Delta z^{\rm MSL}_{i,j,k} = H_{i,j} \, . $$

  The unperturbed levels are then found from

  $$ \label{Eq:GeneralVerticalCoords} z^{\rm MSL}_{i,j,k} = \left\{ \begin{array}{lcc} \displaystyle\frac{-1}{2} \Delta z^{\rm MSL}_{i,j,1} & {\rm if} & k = 1 \\ z^{\rm MSL}_{i,j,k-1} - \displaystyle\frac{1}{2} \left( \Delta z^{\rm MSL}_{i,j,k-1} + \Delta z^{\rm MSL}_{i,j,k} \right) & {\rm if} & k > 1 \end{array} \right. $$

  (60)

  Note that our new general coordinate scheme contains schemes (a–d) as special cases. Hence, schemes (a–d) are now implemented by specifying $\Delta z^{\rm MSL}_{i,j,k}$ outside the model, according to their respective rules, and using the resulting $\Delta z^{\rm MSL}_{i,j,k}$ as input to the general coordinate scheme.

1.2 Appendix 1.2: Fluxes through boundaries of computational cells

To complete the conservative spatial discretizations of Section 2.4, we first establish some notation. Values taken at the centers of tracer volumes have integer indices, e.g., T _i,j,k, while values taken at the centers of velocity volumes have odd-half integer indices, e.g., $\mathbf u_{i+\frac{1}{2},j+\frac{1}{2},k}$. In the vertical, values taken at either the centers of tracer or velocity volumes have integer indices while those at the tops or bottoms of the computational volumes have odd-half integer indices, e.g., $\omega_{i,j,k+\frac{1}{2}}$. Using these rules, we define the following averaging and differencing operators:

$$ \begin{array}{rll} \label{Eq:avgop} \langle \phi \rangle^x_{i,j,k} &=& \frac{1}{2} \big(\phi_{i+\frac{1}{2},j,k} + \phi_{i-\frac{1}{2},j,k} \big) \\ \langle \phi \rangle^y_{i,j,k} &=& \frac{1}{2} \big(\phi_{i,j+\frac{1}{2},k} + \phi_{i,j-\frac{1}{2},k} \big)\\ \langle \phi \rangle^z_{i,j,k} &=& \frac{1}{2} \big(\phi_{i,j,k+\frac{1}{2}} + \phi_{i,j,k-\frac{1}{2}} \big) \\ \delta^x(\phi)_{i,j,k} &=& \phi_{i+\frac{1}{2},j,k} - \phi_{i-\frac{1}{2},j,k} \\ \delta^y(\phi)_{i,j,k} &=& \phi_{i,j+\frac{1}{2},k} - \phi_{i,j-\frac{1}{2},k}\\ \delta^z(\phi)_{i,j,k} &=& \phi_{i,j,k-\frac{1}{2}} - \phi_{i,j,k+\frac{1}{2}}. \end{array} $$

Note that in the above, i and j increase with increasing x and y while k increases with decreasing depth (negative below sea level).

Now we can define the fluxes through the sides of the computational cells. We start with the “flux velocities” evaluated at the centers of the sides. Following Dukowicz and Smith (1994, Appendix E), we define the integrated flows through the “east” and “north” lateral walls of the tracer volumes as

$$ \upsilon^n_{i+\frac{1}{2},j,k} = \Delta y_j \left\langle \Delta z^n {u^\prime} + \frac{1}{2} \Delta z^{n-1} \left(U^n + U^{n-1} \right) \right\rangle^y_{i+\frac{1}{2},j,k} \, , $$

$$ \nu^n_{i,j+\frac{1}{2},k} = \Delta x_i \left\langle \Delta z^n {v^\prime} + \frac{1}{2} \Delta z^{n-1} \left(V^n + V^{n-1} \right) \right\rangle^x_{i,j+\frac{1}{2},k} \, , $$

while at the velocity boxes, we define the integrated flows through the “east” and “north” lateral walls as

$$ \begin{array}{lll} &&{\kern-6pt} \upsilon^n_{i+1,j+\frac{1}{2},k} \\ && = \Delta y_{j+\frac{1}{2}} \Bigg\langle \Bigg\langle \Bigg\langle \Delta z^n {u^\prime} + \frac{1}{2} \Delta z^{n-1} \left(U^n + U^{n-1} \right) \Bigg\rangle^y \Bigg\rangle^x \Bigg\rangle^y_{i+1,j+\frac{1}{2},k} \, , \end{array} $$

$$\begin{array}{lll} &&{\kern-6pt} \nu^n_{i+\frac{1}{2},j+1,k} \\ &&= \Delta x_{i+\frac{1}{2}} \Bigg\langle \Bigg\langle \Bigg\langle \Delta z^n {v^\prime} + \frac{1}{2} \Delta z^{n-1} \left(V^n + V^{n-1} \right) \Bigg\rangle^x \Bigg\rangle^y \Bigg\rangle^x_{i+\frac{1}{2},j+1,k} \, . \end{array} $$

These particular spatial averagings are chosen to match the discrete transport constraint (Eq. 65 in “Appendix 1.4”). The new aspect here is the temporal evaluations. The baroclinic velocity components are evaluated at time n. However, the timings for the barotropic components are, again, chosen to match the transport constraint (Eq. 65). Also note that these timings assume θ = 1. To get the corresponding flows through the “west”(“south”) lateral walls, simply decrement i(j ) by one.

To evaluate the fluxes through the tops of the computational volumes, we use the above definitions in Eq. 37. At tracer volumes this yields

$$ \begin{array}{lll} \label{Eq:VertFluxTbox} &&{\kern-6pt} \delta^z (\omega^n)_{i,j,k} \Delta x_i \Delta y_j + \delta^x (\upsilon^n)_{i,j,k} + \delta^y(\nu^n)_{i,j,k}\\ &&{\kern2pt} + \frac{\Delta {\cal V}^n_{i,j,k}}{H_{i,j}+\eta^n_{i,j}} \frac{\delta (\eta_{i,j})^{n,n-2}}{\tau} = 0 \end{array} $$

(61)

while at velocity volumes, we get

$$ \begin{array}{lll} &&{\kern-6pt}\delta^z(\omega^n)_{i+\frac{1}{2},j+\frac{1}{2},k} \Delta x_{i+\frac{1}{2}} \Delta y_{j+\frac{1}{2}}\\ &&{\kern2pt} +\, \delta^x(\upsilon^n)_{i+\frac{1}{2},j+\frac{1}{2},k} + \delta^y(\nu^n)_{i+\frac{1}{2},j+\frac{1}{2},k}\\ \label{Eq:VertFluxVbox} &&{\kern2pt} +\, \frac{\Delta {\cal V}^n_{i+\frac{1}{2},j+\frac{1}{2},k}} {H_{i,j}+\langle\langle\eta^n\rangle^x\rangle^y_{i+\frac{1}{2},j+\frac{1}{2}}} \frac{\delta \left(\langle\langle\eta\rangle^x\rangle^y_{i+\frac{1}{2},j+\frac{1}{2}}\right)^{n,n-2}} {\tau}= 0 \, .\\ \end{array} $$

(62)

Using these definitions of the fluxes through the boundaries of the computational volumes, we can now simply write the discrete advection operator as

$$ \begin{array}{rll} \breve{\Gamma}(\phi)^n_{i,j,k} &=& \delta^x\left(\langle\phi^n\rangle^x\upsilon^n\right)_{i,j,k} + \delta^y\left(\langle\phi^n\rangle^y\nu^n\right)_{i,j,k}\\ &&+\; \delta^z\left(\langle\phi^n\rangle^z\omega^n\right)_{i,j,k} \Delta x_i \Delta y_j \, . \end{array} $$

This formulation is valid for both tracer and velocity computational volumes, with the understanding that for velocity volumes the i,j indices are shifted by one half.

We have evaluated the pressure force term, $-{1 \over \rho_0} \int_{{\cal S}^n} p^n_h \, \hat n_h \cdot d{\boldsymbol{\cal A}}$, both by directly discretizing the integrals of pressure along the cell walls (including the horizontal contributions from the sloping cell tops and bottoms) and by interpolating the pressure to the corresponding velocity depths and evaluating the differential gradient. Both give similar results, but the integral evaluation is conservative and produces less noise in the resulting velocities (especially near sloping bottoms).

1.3 Appendix 1.3: Open boundary conditions

For u ^′, T, S, and η, the application of boundary conditions is straightforward. Our options (see Haley et al. 2009; Lermusiaux 1997) include using values based on data, applying radiation conditions (Orlanski 1976; Spall and Robinson 1989), or, following Perkins et al. (1997), using radiation conditions to correct the provided values. For nested subdomains, we have first used the interpolated values directly or with Perkins et al. (1997) corrections. Some other promising options we have explored with nested subdomains include using the coarse grid values in a narrow buffer zone around the fine domain, which reduces discontinuities. Another important multiscale conservative boundary condition option is to feedback the averages of the fluxes across the boundary walls shared with the large domain (Fig. 3). These include the advective fluxes of momentum and tracers, the pressure force, and the diffusive fluxes of momentum and tracers.

We still need an additional boundary condition for $\overline{\hat{\cal F}^{n,n-1}}$ since we are unable to directly evaluate Eq. 42 at the boundaries. To derive this boundary condition, we recast Eq. 42 in the form of Eq. 24 and solve for $\overline{\hat{\cal F}^{n,n-1}}$:

$$ \label{Eq:SetUpUbarBC} \overline{\hat{\cal F}^{n,n-1}} = \frac{\delta(\mathbf{U})}{\tau} + f \hat k \times {\mathbf{U}}^\alpha + g \nabla \eta^\alpha_{i,j} \, . $$

(63)

Now, the right-hand side of Eq. 63 is made up entirely of quantities that can be directly evaluated at the boundary of the velocity grid. For the free surface, we have found that it is more stable to rewrite Eq. 63 in terms of transports:

$$ \begin{array}{rll} \label{Eq:UbarotropicBC} \overline{\hat{\cal F}^{n,n-1}} &=& \frac{1}{H+\eta^n} \left\{\frac{\delta\left[\left(H+\eta\right)\mathbf{U}\right]}{\tau} + f \hat k \times \left[\left(H+\eta\right){\mathbf{U}}\right]^\alpha\right\}\\ &&+\, g \nabla \eta^\alpha_{i,j} \, . \end{array} $$

(64)

Note: when evaluating Eq. 64, only values at time t _n + 1 are taken from the provided fields (or nesting interpolations). The fields at times t _n and t _n − 1 are both already in memory and in primitive equation balance. They are combined with the t _n + 1 fields to evaluate (Eq. 64).

Applying Perkins et al. ( 1997 ) boundary conditions:

Following the algorithm of Perkins et al. (1997), corrections to the provided values (and nesting interpolation values) are obtained by applying the Orlanski radiation algorithm to the difference between the PE model values and these provided values and using these differences to correct the boundary values.

For the barotropic transport, however, this is only done for the tangential component to the boundary. The correction to the normal component is derived from the correction obtained for the surface elevation, Δη, and the barotropic continuity equation

$$ \frac{\partial \Delta \eta}{\partial t} + \nabla \, \cdot \, \left[ (H+\eta) \Delta {\mathbf U} \right] = 0 \, . $$

1.4 Appendix 1.4: Maintaining the vertically integrated conservation of mass

To see how the free-surface algorithm maintains the vertically integrated conservation of mass, start from Eq. 44, multiply by θ(H + η ⁿ), and take the divergence of the result to get

$$ \begin{array}{lll} \nabla &\cdot& \left[ \left(H+\eta^n\right) \theta \mathbf U^{n+1}\right] = \nabla \cdot \left[ \left(H+\eta^n\right) \theta {\hat {\bf U}}^{n+1}\right]\\ &-& \alpha \theta g \tau \nabla \cdot \left[ \left(H+\eta^n\right) \nabla \delta \eta \right] + \theta \nabla \cdot \left(\left.\mathbf{u^\prime}^n\right|_\eta \frac{\delta \eta}{\tau}\right) \, . \end{array} $$

Substitute for the right-hand side of the above equation from Eq. 43 and rearrange to obtain

$$ \label{Eq:FreeSurfaceVertIntCont} \frac{\delta \eta}{\tau} {\kern-.5pt} +{\kern-.5pt} \frac{1}{2} \nabla {\kern-.5pt}\cdot{\kern-.5pt} \left[ \left(H{\kern-.5pt}+{\kern-.5pt}\eta^n\right) \left( \theta \mathbf U^{n+1} {\kern-.5pt}+{\kern-.5pt} \mathbf U^n {\kern-.5pt}+ (1 {\kern-.5pt}-{\kern-.5pt} \theta) \mathbf U^{n-1} \right) \right] {\kern-.5pt}={\kern-.5pt} 0 . $$

(65)

Equation 65 represents the discrete form of the barotropic continuity enforced by the free-surface algorithm. Imbalances in Eq. 65 produce unrealistic vertical velocities via Eqs. 61 and 62.

However, as illustrated by the above derivation, Eq. 65 is only satisfied to the same degree that Eqs. 43 and 44 are satisfied. This places restrictions on the valid avenues for nesting. For example, we can safely replace the coarse domain estimates of $(H+\eta){\hat {\bf U}}^{n+1}$ with averages from the fine domain without disturbing (Eq. 65). Moving this exchange one step later in the algorithm and trying to average $(H+\eta)[\theta{\hat {\bf U}}^{n+1} + \mathbf U^n + (1 - \theta) \mathbf U^{n-1}]$ would violate Eq. 43, in the sense that we would not be able to make the last substitution leading to Eq. 65 and hence we would violate Eq. 65. Similar to the bias argument of Section 4.7, these fields and their errors will persist over many time steps, leading to unsustainable vertical velocities.

Appendix 2: Additional details on running our fully implicit two-way nested free-surface PE model

1.1 Appendix 2.1: Setting up domains

1.1.1 Appendix 2.1.1: Topography

There are two main issues when defining topographies for nested simulations. The first is that the finer resolution grid can support finer topography scales, including sharper gradients. The bathymetry on the finer grid is not an interpolation of the coarser grid bathymetry, but the coarser grid bathymetry is a coarse-control-volume average of the finer grid bathymetry. The refinement in topographic scales can lead to abrupt artificial discontinuities in the topography where the fine and coarse domains meet. This can be exacerbated by conditioning the topography (Haley and Lozano 2001) to control the hydrostatic consistency condition (Haney 1991). For a given value of the hydrostatic consistency factor (roughly proportional to $\frac{d\mathbf{x} \times \nabla h}{h}$), the finer resolution domain can support steeper bathymetric features (e.g., shelfbreak). To ensure a smooth transition, we define a band of points around the outer edge of a fine domain (e.g., a band from the boundary to 6 points inside the boundary, see also Penven et al. 2006). In this band, we replace the fine grid topography with a blend of the coarse and fine grid topographies:

$$ \label{Eq:TopoBlending} h_{\rm blend} = \alpha h_{\rm f\/ine} + (1-\alpha) h_{\rm coarse} $$

(66)

where α varies from zero at the boundary to one at the inner edge of the band (e.g., 6 points).

The second issue comes about from the nesting algorithm itself. As mentioned in Section 3, we force the undisturbed vertical grid , $z^{\rm MSL}_{i,j,k}$, to satisfy the nesting rules of Eqs. 46 and 48. To ensure that the topographies in nested domains satisfy Eq. 46 and 48 and the blending Eq. 66 ,we usually follows these steps:

• 1. Apply the nesting constraints on the unconditioned topographies. Starting from the smallest domain, average the fine grid topographies on the successively larger topographies according to Eq. 46. Then starting from the coarsest domain, interpolate the topographies to the boundaries of the successively smaller domains according to Eq. 48.

• 2. Starting from the largest domain, apply the conditioning. After the largest domain is conditioned, apply the blending Eq. 66 to the second largest. Condition that domain and repeat the blending-conditioning cycle with the successively smaller domains.

• 3. Reapply the nesting constraints on the conditioned topography. Repeat step 1.

1.1.2 Appendix 2.1.2: Land masking

The first constraint for masking occurs at the boundaries of the finer domains. Considering any two nested domains, we want continuity of the masks across the domain boundary. In other words, a coastline that crosses the boundary of the fine domain should not have a jump or jog at the boundary of the fine domain. Enforcing this consistency, along with boundary constraints on the topography (“Appendix 2.1.1”), enforces consistent estimates of the areas of the lateral boundaries of the fine domain as measured in both the coarse and fine grids.

The second constraint is to have a certain degree of consistency in defining land and sea in the interior of the fine domain. This is a less exact statement because the fine domain supports a more detailed resolution of the land/sea boundary than the coarse domain. Because of the superior resolution, we take the view that the land mask in the interior of the fine domain is “more correct” than the coarse domain mask. Since we use collocated grids, this provides us with a simple algorithm for resetting the coarse mask. For each coarse grid point fully supported by fine grid points, we count how many of the supporting fine grid points are land and how many are sea. If at least one half the fine grid points are sea, the coarse grid point is marked as sea; otherwise, it is masked as land.

Our general procedure is to first define the land mask for the largest (coarsest) domain. Then use that mask to define a crude first guess for the mask in the fine domain. We then reset the interior nodes of the fine mask to better resolve the coasts (leaving a narrow band around the exterior untouched to ensure continuity through the boundary). If we have more then two domains, we use the current domain to initialize the mask for the next finest domain and repeat. When we finish the mask in the smallest (finest) domain, we use that mask to reset the mask in the next coarser domain, using the above sea/land counting procedure. We then examine the modified mask in that next coarser domain to eliminate any spurious artifacts that may have been created (e.g., a narrow mouth of a bay may have been closed leaving an isolated “lake” that we do not need to maintain). We repeat with the next coarser domain and so on until we get back to the coarsest domain.

1.2 Appendix 2.2: Initialization

Our common situation is to estimate the best initial synoptic state from temperature and salinity data (in situ, climatologies, satellite, etc.) but little or no direct velocity data. Our initialization scheme for this situation is described next, focusing mainly on the nesting considerations, first briefly for the rigid-lid procedures and then the extensions for initializations with a free surface.

1.2.1 Appendix 2.2.1: Rigid lid

Our procedures for rigid lid initializations in nested grids are based on (e.g., Haley et al. 2009). Starting from temperature and salinity data, climatologies, etc., we create 3D estimates of temperature and salinity, often using objective analyses (Carter and Robinson 1987; Agarwal and Lermusiaux 2010). From these 3D temperature and salinity estimates, we construct density (Eq. 6) and the hydrostatic pressure (Eq. 9). We then estimate the total velocity using the rigid lid geostrophic relation

$$ f \hat k \times (\mathbf u - \mathbf u_{\rm ref}) = \frac{1}{\rho_0} \nabla p_h = \frac{g}{\rho_0} \int^z_{Z_{\rm ref}} \nabla \rho \, d\zeta$$

where Z _ref is a suitably chosen reference level, which can be a “level of no motion”, u _ref is the absolute velocity at that depth, and we have interchanged the horizontal gradient with the vertical integral. When evaluating $\nabla \rho$ at a particular depth, if any of the ρ values used for the gradient would be below topography, we set $\nabla \rho$ to zero. To enforce no penetration of land, we find a stream function, ψ, which satisfies $\nabla^2 \psi = \nabla \times \mathbf u$ with ψ set to be constant along coasts. From this ψ, we recompute the velocity. We decompose this velocity into barotropic and baroclinic parts (Eq. 8). The baroclinic portion is fine as is, but barotropic velocities at this stage generally do not satisfy the nondivergence of transport. To enforce this, we define a transport stream function $\hat k \times \nabla \Psi = H \mathbf U$ and fit it to our estimated barotropic velocities via the Poisson equation

$$ \nabla \times \left( \frac{\hat k}{H} \times \nabla \Psi \right) = \nabla \times \mathbf U \, . $$

We derive Dirichlet boundary conditions for the above by first noting that the tangential derivative of Ψ to the boundary equals the normal component of transport, H U, through the boundary. We then integrate this relation around the boundary to obtain the Dirichlet values. For domains with islands, we also need to provide constant values for Ψ along the island coasts. We do this in a two-stage process in which we first compute Ψ assuming all the islands are open ocean. We then use that initial guess to derive constant island values that minimize the relative interisland transports using fast marching methods (Agarwal 2009).

Nesting considerations

For nesting the initial temperature, salinity, other tracers, and baroclinic velocity, we can directly enforce some conservation constraints by averaging estimates from finer to coarser grids. For the transport stream function, we go to the additional step of generating the Dirichlet boundary values for the Poisson equation in the fine domain by interpolating the stream function values from the coarse domain. This ensures that the same constant of integration is used for both domains and that the net flows through the fine domain are consistent in both the coarse and fine grids. For island values, if the island is represented in both the coarse and fine domains, the coarse domain value is used. If the island is only in the fine domain, then the procedure of the preceding paragraph is used.

1.2.2 Appendix 2.2.2: Free surface

The starting point for the free-surface initialization scheme is the above rigid-lid initialization. We start by explicitly computing the final, rigid-lid barotropic velocities from

$$ \mathbf U = \frac{\hat k}{H} \times \nabla \Psi \, . $$

We next create an equation for the initial surface elevation. We start from the geostrophic approximation with the full pressure

$$ f \hat k \times \mathbf u = g \nabla \eta + \frac{1}{\rho_0} \nabla p_h \, .$$

Integrating this equation in the vertical from − H to 0 and isolating η results in

$$ \label{Eq:FreeSurfaceGradientIC} g H \nabla \eta = f \hat k \times H \mathbf U - \frac{1}{\rho_0} \int^0_{-H} \nabla p_h dz \, . $$

(67)

Finally, we take the divergence of Eq. 67 to get

$$ \label{Eq:FreeSurfaceDivergenceIC} \nabla \cdot \left( g H \nabla \eta \right ) = \nabla \cdot \left( f H \mathbf U \right) - \frac{1}{\rho_0} \nabla \cdot \left( \int^0_{-H} \nabla p_h dz \right) \, . $$

(68)

To generate Dirichlet boundary values for Eq. 68, we integrate the tangential components of Eq. 67 around the boundary. Because the coastal boundary condition is zero normal derivative, no special action needs to be taken for islands.

Once an initial value for η is constructed, then, by Eq. 35, the initial depths are recomputed. The tracers (temperature, salinity, etc.) and baroclinic velocity are re-interpolated to these new initial depths. Finally, the barotropic velocities from the rigid-lid calculation are rescaled to preserve the transports:

$$ \mathbf U_{\rm free\ surface} = \frac{H}{H+\eta} \mathbf U_{\rm rigid\ lid} \, . $$

Nesting considerations

These are the same as for the rigid-lid case. The additional detail is that now we also interpolate the coarse grid estimate of η to generate Dirichlet boundary values for solving Eq. 68 in the fine domain.

1.3 Appendix 2.3: Tidal forcing

1.3.1 Appendix 2.3.1: Constructing the tidal forcing

When adding tidal forcing to our simulations, our underlying assumption is that our regional domains are small enough so that the tidal forcing through the lateral boundaries completely dominates the local body force effects. To model these lateral forcings, we employ the linearized barotropic tidal model (Logutov 2008; Logutov and Lermusiaux 2008). We use a shallow water spectral model and generate 2D fields for the amplitude and phase of tidal surface elevation and the barotropic tidal velocity. We dynamically balance these barotropic tidal fields with our best available topographic and coastal data along with the best exterior barotropic tidal fields (e.g., Egbert and Erofeeva 2002). Once we have constructed our tidal fields for the desired modes, we can simply evaluate them for any time.

The above procedures can provide tidal fields on different grids than used by our PEs. For example, Logutov (2008) and Logutov and Lermusiaux (2008) are formulated on a C-grid, instead of the B-grid being used here. In particular, this means that tidal fields interpolated from these grids will not, in general, exactly satisfy the same discrete continuity as in our grid. Our experience shows that satisfying the same discrete continuity leads to more robust solutions. To enforce this constraint, we solve the constrained minimization problem

$$ \begin{array}{rll} {\cal J} & = & \int \left\{ \alpha \eta^*_1 \eta_1 + \theta \beta \mathbf U^*_1 \cdot \mathbf U_1 + \phi \beta \mathbf U^* \cdot \mathbf U \right.\\ & & \left. + \lambda \Re\left[i\omega \eta + \nabla \cdot \left(H \mathbf U \right) \right] +\, \gamma \Im \left[i\omega \eta + \nabla \cdot \left(H \mathbf U \right) \right] \vphantom{\eta^*_1}\right\} d {\cal V} \nonumber \end{array} $$

where η ₀, U ₀ as the complex tidal surface elevation and barotropic tidal velocity interpolated from the original grid and η ₁, U ₁ are the additive “correction” complex tidal surface elevation and barotropic tidal velocity that minimize $\cal J$,

$$ \eta = \eta_0 + \eta_1 \; ; \quad \mathbf U = \mathbf U_0 + \mathbf U_1 $$

α and β are the weights (including nondimensionalizing factors), λ and γ are the Lagrange multipliers, the superscript ^* indicate complex conjugation, $\Re$ and $\Im$ refer to the real and imaginary parts, and θ, φ are penalty parameters to inhibit unreasonably large total velocities. Using the calculus of variations, the above minimization is equivalent to solving the following system of equations

$$ \begin{array}{rll} \omega \eta_1 - \nabla \cdot \left[ \frac{H^2}{(\theta + \phi) \omega \beta} \nabla \left( \alpha \eta_1 \right) \right] & = & - \omega \eta_0 + i \nabla \notag\\&&\cdot \left( \frac{\theta}{\theta + \phi} H \mathbf U_0 \right) \\ \left. \alpha \right|_{\rm open\ boundary} & = & 0 \\ \left.\left( i \omega \eta_1 + \sqrt{g H} \frac{\partial \eta_1}{\partial n}\right) \right|_{\rm open\ boundary} & = & 0 \\ \mathbf U_1 & = & - \frac{\phi}{\theta + \phi} \mathbf U_0 \nonumber\\&&- i \frac{H}{(\theta + \phi) \omega \beta} \nabla \left(\alpha \eta_1\right) \end{array} $$

Note that the radiation boundary condition does not come from the variations but is a useful addition we are free to make after obtaining α = 0 from the variations.

1.3.2 Appendix 2.3.2: Applying the tidal forcing

We use the barotropic tides both for initialization and for boundary forcing. For the surface elevation, we simply superimpose the tidal surface elevation with the subtidal elevation estimated in “Appendix 2.2.2”. For initialization, this superposition is done over the entire area before the final vertical interpolation of tracers and u ^′. For lateral forcing, this is done at run time in the PE model at the exterior boundaries (and also along two bands inside these boundaries for Perkins et al. (1997) boundary conditions). The resulting boundary values are used to generate Dirichlet boundary conditions for Eq. 43. A similar procedure is used for the barotropic velocities with two notable differences. First, the superposition is performed to preserve transport:

$$ \begin{array}{rll} \left(H + \eta_{\rm superimposed}\right) \mathbf U_{\rm superimposed} &=& \left(H + \eta_{\rm subtidal}\right) \mathbf U_{\rm subtidal}\\ &&+\, H \mathbf U_{\rm tidal} \, . \end{array} $$

Note that the tidal velocity is only multiplied by the undisturbed water depth. This is because we are using a linearized tidal model. The second difference is that the run-time boundary values of the barotropic velocity are used for Eq. 64, not directly applied to the final barotropic velocities.

Nesting considerations

For initialization, the process is as for the unnested case. The superpositions described above are done for the initial conditions of each domain. For the lateral forcing, however, the barotropic tidal fields are only applied at the boundaries of the coarsest domain. The reason being that applying the barotropic tidal fields at the boundary of the coarsest domain can produce the full tidal response in the interior, and hence, the barotropic tidal fields are unnecessary for the nested subdomains.

1.4 Appendix 2.4: Solving the equation for the surface elevation

Equation 43, with Dirichlet boundary conditions, represents an elliptic system of equations for the surface elevation, η. To numerically solve this system, we use a preconditioned conjugate gradient solver for sparse matrices (e.g., SPARSKIT; Saad 2009). A typical convergence test for such an iterative solver is an integrated measure of the reduction in the norm of the residual over all points. Specifically, if r is the residual of the current solver iteration and r ₀ is the residual of the initial guess, the convergence test is

$$\Vert r \Vert \le \tau_r \Vert r_0 \Vert + \tau_a$$

where τ _r is the relative tolerance and τ _a is the absolute tolerance. In practice, we tend to use very small values (10^− 12 and 10^− 25, respectively) to ensure a tight convergence. Using the results of “Appendix 3”, we also supplement this global constraint with the following point-wise constraint:

$$ \begin{array}{rll} \left\vert \frac{\partial \delta\eta^k}{\partial x} - \frac{\partial \delta\eta^{k-1}}{\partial x} \right\vert &\le& \tau^g_r \left\vert \frac{\partial \delta\eta^k}{\partial x} \right\vert + \tau_a \; ; \quad \\ \left\vert \frac{\partial \delta\eta^k}{\partial y} - \frac{\partial \delta\eta^{k-1}}{\partial y} \right\vert &\le& \tau^g_r \left\vert \frac{\partial \delta\eta^k}{\partial y} \right\vert + \tau_a\\ \left\vert \delta\eta^k - \delta\eta^{k-1} \right\vert \le \tau^g_r \left\vert \delta\eta^k \right\vert + \tau_a \end{array} $$

where the superscript k refers to the iteration number and $\tau^g_r$ is the relative tolerance for the gradient test (typically around 10^− 8). Here we test on both δη and its gradients to ensure the relative convergence of the barotropic velocities (Eq. 44).

Since we have discretized our equations on the B-grid, both Eq. 43 and, especially, Eq. 68 possess a well-known checkerboard mode in their null spaces (Deleersnijder and Campin 1995; le Roux et al. 2005; Wubs et al. 2006). For realistic geometries, we found that applying a Shapiro filter (Shapiro 1970) to the solution was sufficient to suppress the noise while maintaining the correct physical features. The one case where this approach failed was in creating an initialization for an idealized flow in a periodic channel. The lack of Dirichlet boundary values in that case and corresponding lack of structure they would have imposed allowed the checkerboard mode to suppress all other structures. To control this error, the matrix in Eq. 68 was augmented with a Laplacian filter (Deleersnijder and Campin 1995; Wickett 1999) to prevent the appearance of this mode. Again, this filter was only needed for the idealized periodic channel flow.

Appendix 3: Review of rigid-lid nesting algorithm

One of our first nesting schemes (see Section 4.1) for the nonlinear free-surface version of the PEs was built by analogy with our two-way nesting scheme (e.g., Spall and Holland 1991; Fox and Maskell 1995; Sloan 1996) for the HOPS rigid-lid formulation of the PEs. We have used this scheme in a variety of dynamical situations (Robinson et al. 2002; Leslie et al. 2008; Haley et al. 2009) and we present it next.

We start with the rigid-lid PEs:

$$ \nabla \cdot \mathbf u + \frac{\partial w}{\partial z} = 0 \, , $$

$$\frac{\partial \mathbf{u^\prime}}{\partial t} + f \hat k \times \mathbf{u^\prime} = \widehat{\cal F} - \overline{\widehat{\cal F}} \, , $$

$$\nabla \times \left( \frac{\hat k}{H} \times \nabla \delta \Psi \right) + \alpha_{\rm RL} \nabla \times \left[ f \, \hat k \times \left( \frac{\hat k}{H} \times \nabla \delta \Psi \right) \right] = \nabla \times \overline{\widehat{\cal F}} \,\, , $$

$$p_h = - \int^0_z \rho g \, d\zeta \, , $$

$$\frac{D T}{D t} = F^T \, , $$

$$\frac{D S}{D t} = F^S \, , $$

$$\rho = \rho (z,T,S) $$

where

$$ \mathbf u = \mathbf{u^\prime} + \mathbf U \; ; \quad \mathbf U = \frac{1}{H} \int^0_{-H} \mathbf u \, dz \; ; \quad H \mathbf U = \hat k \times \nabla \Psi \, , $$

$$ \widehat{\cal F} = -{1 \over \rho_0} \nabla p_h - {\boldsymbol{\Gamma}}(\mathbf u) + \mathbf F \; ; \quad \overline{\widehat{\cal F}} = \frac{1}{H} \int^0_{-H} \widehat{\cal F} \, dz $$

$$ {\boldsymbol \Gamma}(\mathbf u) = \left( \begin{array}{c} \Gamma(u) \\ \Gamma(v)\end{array}\right) \; ; \quad \Gamma(\phi) = u \frac{\partial \phi}{\partial x} + v \frac{\partial \phi}{\partial y} + w \frac{\partial \phi}{\partial z} . $$

Here we have introduced the transport stream function, Ψ, to replace the barotropic velocity as a state variable. Using this notation, we write our rigid-lid nesting algorithm as

• 1. Solve the rigid-lid PEs simultaneously in each domain for (u ^{′ n + 1}, T ^n + 1, S ^n + 1)

• 2. Replace values in the coarse domain at overlap nodes with the following averages from the fine domain values

  $$ \phi^{n+1}_{i_c,j_c,k} = \frac{1}{\Delta {\cal V}_{i_c,j_c}} \sum\limits^{j_{fc}+r_h}_{j=j_{fc}-r_h} \sum\limits^{i_{fc}+r_h}_{i=i_{fc}-r_h} \phi^{n+1}_{i,j,k} \Delta {\cal V}_{i,j,k} \, , $$

  (69)
  $$ \left(\nabla \times \overline{\widehat{\cal F}}\right)_{i_c,j_c} = \frac{1}{\Delta {\cal A}_{i_c,j_c}} \sum\limits^{j_{fc}+r_h}_{j=j_{fc}-r_h} \sum\limits^{i_{fc}+r_h}_{i=i_{fc}-r_h} \left(\nabla \times \overline{\widehat{\cal F}}\right)_{i,j} \Delta {\cal A}_{i,j}\nonumber\\ $$

  (70)

  where r is the ratio of the resolution of the coarse grid to fine grid,

  $$ \begin{array}{lll} &&{\kern-6pt} \phi = \mathbf{u^\prime}, T, S \; ; \quad \Delta {\cal V}_{i,j,k} = \Delta x_{i,j}\Delta y_{i,j} \Delta z_{i,j,k} \; ; \\ &&{\kern-6pt} \Delta {\cal A}_{i,j} = \Delta x_{i,j}\Delta y_{i,j} \; ; \quad r_h = \lfloor r/2\rfloor \; . \end{array} $$

• 3. In the coarse domain, solve the rigid-lid barotropic momentum equation for Ψ^n + 1.

• 4. Using piece-wise bi-cubic Bessel interpolation, $\cal B$, replace values in the fine grid boundary with values interpolated from the coarse grid

  $$ \begin{array}{rll} \phi^{n+1}_{i_{fb},j_{fb},k} &=& {\cal B}\left(\phi^{n+1}_{i_c,j_c,k}\right) \, , \\ \mathbf{u^\prime}^{n+1}_{i_{fb},j_{fb},k} \Delta z_{i_{fb},j_{fb},k} &=& {\cal B}\left(\mathbf{u^\prime}^{n+1}_{i_c,j_c,k} \Delta z_{i_c,j_c,k}\right) \, , \end{array} $$

  where

  $$ \phi = T,\,S,\,\Psi \, . $$

  Note that the interpolation of baroclinic velocity is written in terms of transport rather than velocity. This is done to preserve the baroclinicity of u ^′ in the fine domain.

• 5. In the fine domain, solve the barotropic momentum equation for Ψ^n + 1.

We found that this scheme maintains consistency between the estimates on the coarse and fine grids provided that a sufficiently stringent convergence criterion is used when solving for Ψ^n + 1. Only using a test on the integral of the residuals did not always maintain the consistency. Such a test is global in nature and can give different convergence results in the coarse and fine domains. Instead, we found the best results occurred when we supplemented the residual tests with the following point-wise test:

$$ \left| \frac{\frac{\partial \Psi^{k+1}}{\partial x} - \frac{\partial \Psi^k}{\partial x}} {\frac{\partial \Psi^k}{\partial x}} \right| \le \epsilon \; ; \quad \left| \frac{\frac{\partial \Psi^{k+1}}{\partial y} - \frac{\partial \Psi^k}{\partial y}} {\frac{\partial \Psi^k}{\partial y}} \right| \le \epsilon $$

where here the superscripts k refer to the iteration count in the solver and ϵ is the relative tolerance. We test the derivative of Ψ instead of Ψ so that the convergence is on the relative change in the transport, a physically more meaningful quantity than Ψ itself. This test ensures that in both domains the iterative solution proceeds at least until the specified tolerance is reached at every point, thereby ensuring at least that level of consistency between the solutions in the coarse and fine domains.