The Meridionally Averaged Model of Eastern Boundary Upwelling Systems (MAMEBUSv1.0)

Eastern boundary upwelling systems (EBUSs) are physically and biologically active regions of the ocean with substantial impacts on ocean biogeochemistry, ecology, and global fish catch. Previous studies have used models of varying complexity to study EBUS dynamics, ranging from minimal two-dimensional (2-D) models to comprehensive regional and global models. An advantage of 2-D models is that they are more computationally efficient and easier to interpret than comprehensive regional models, but their key drawback is the lack of explicit representations of important three-dimensional processes that control biology in upwelling systems. These processes include eddy quenching of nutrients and meridional transport of nutrients and heat. The authors present the Meridionally Averaged Model of Eastern Boundary Upwelling Systems (MAMEBUS) that aims at combining the benefits of 2-D and 3-D approaches to modeling EBUSs by parameterizing the key 3-D processes in a 2-D framework. MAMEBUS couples the primitive equations for the physical state of the ocean with a nutrient–phytoplankton–zooplankton–detritus model of the ecosystem, solved in terrain-following coordinates. This article defines the equations that describe the tracer, momentum, and biological evolution, along with physical parameterizations of eddy advection, isopycnal mixing, and boundary layer mixing. It describes the details of the numerical schemes and their implementation in the model code, and provides a reference solution validated against observations from the California Current. The goal of MAMEBUS is to facilitate future studies to efficiently explore the wide space of physical and biogeochemical parameters that control the zonal variations in EBUSs.

1 Introduction

Eastern boundary upwelling systems (EBUSs) are among of the most biologically productive regions in the ocean, supporting diverse ecosystems and contributing to a significant portion of the global fish catch (Bakun and Parrish, 1982). The characteristic wind-driven upwelling dominant in EBUSs is forced by an equatorward meridional wind stress that decreases toward the shore, driving a zonal Ekman transport offshore. The resulting Ekman pumping brings cold, nutrient-rich water to the surface, fueling primary productivity (Jacox and Edwards, 2012; Chavez and Messié, 2009; Rykaczewski and Dunne, 2010).

The upwelling-favorable winds also drive baroclinic, equatorward geostrophic current, which sheds mesoscale eddies (Colas et al., 2013). Together with offshore Ekman transport, mesoscale eddies redistribute nutrients zonally and subduct nutrients and other tracers into the ocean subsurface (Capet et al., 2008; Gruber et al., 2011; Renault et al., 2016). The resulting cross-shore gradient of nutrients at the surface supports a zonal variation in the abundance of phytoplankton, with high biomass and chlorophyll nearshore, and low offshore (Chavez and Messié, 2009). The size structure of phytoplankton is similarly affected, with larger cells with higher nutrient demand onshore, and smaller cells offshore (Cabre et al., 2013).

While these qualitative patterns of productivity are common to upwelling systems, previous studies have shown that productivity varies substantially between EBUSs, but the causes of these inter-EBUS variations are not well understood. Possible physical drivers of these inter-EBUS variations include the shape and strength of the wind-stress curl, which set the upwelling strength and source depth (Bakun and Nelson, 1991; Jacox and Edwards, 2012). This in turn controls the energy transferred to the baroclinic eddy field, modulating surface nutrient availability via the “eddy quenching” mechanism (Gruber et al., 2011; Renault et al., 2016). Additionally, inter-EBUS variations may have biogeochemical origins, for example, due differing subsurface oxygen inventories (Chavez and Messié, 2009).

Our understanding of these drivers is hindered in part by the observational limitations and in part by the computational expense of regional models that can resolve the processes mentioned above. A range of models of varying complexity have been used to study EBUSs, from minimal two-dimensional (2-D) models (Jacox and Edwards, 2012; Jacox et al., 2014) to comprehensive regional models (Shchepetkin and McWilliams, 2005; Chenillat et al., 2018). While 2-D models require fewer computational resources than comprehensive regional model studies and thus allow a more comprehensive exploration of the relevant parameter space, they lack the explicit representation of important physical processes that affect biology in upwelling systems (i.e., eddy-quenching and meridional transport of nutrients).

Here, we aim to close the current gap in understanding by developing an idealized, quasi-2-D model of the physics and biogeochemistry of EBUSs. The model includes parameterizations of the key three-dimensional processes, while retaining the computational efficiency of a 2-D model. The model is cast in a residual-mean framework (Plumb and Ferrari, 2005a) in terrain-following coordinates (Song and Haidvogel, 1994) and is referred to as the Meridionally Averaged Model of Eastern Boundary Upwelling Systems (MAMEBUS). A schematic of all the important processes in MAMEBUS is shown in Fig. 1.

The rest of the paper is organized as follows. In Sect. 2, we describe the equations and physical parameterizations implemented in MAMEBUS, including general formulation of tracer advection and diffusion, the time-dependent turbulent thermal wind approximation of the momentum equations (T3W), eddy and boundary layer parameterizations, and our ecosystem formulation. In Sect. 3, we detail the algorithms and discretizations, including mesh specification, vertical coordinate transformation, and time integration. In Sect. 4, we describe the implementation of MAMEBUS including the various options available to the user, parameter choices, initialization, and output. In Sect. 5, we describe reference solutions for MAMEBUS, discussing model sensitivities to changes in bathymetry, wind forcing, and surface heat fluxes. Finally, in Sect. 6, we discuss further model development and future work.

Figure 1 A schematic of the essential components of the Meridionally Averaged Model of Eastern Boundary Upwelling Systems (MAMEBUS). This schematic highlights some components that the user is able to control including the offshore restoring conditions, the eddy mixing along isopycnals, the wind forcing, the surface mixed layer and bottom boundary layer parameterizations, and grid spacing.

2 MAMEBUS framework

MAMEBUS is comprised of a series of components that are necessary to capture physical–biogeochemical dynamics in EBUSs: (1) explicit momentum conservation in form of geostrophic, hydrostatic, and Ekman balances implemented as part of the T3W formulation; (2) eddies and their effect on material transport; (3) surface and bottom boundary layers; (4) nutrient and plankton cycles in form of a size-structured “NPZD”-type model (Banas, 2011).

With the exception of the velocity field, all tracers in MAMEBUS evolve according to the following conservation equation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{phys}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{bio}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{nct}},\end{array}$$

where the bar indicates a meridional average. The key physical tracer that follows Eq. (1) is temperature, θ, which serves as the thermodynamic variable in our model. We choose temperature as our thermodynamic variable because of its important effects on biogeochemistry (Sarmiento and Gruber, 2006). The biogeochemical tracers that are affected by the biogeochemical evolution term, ${\partial }_{\mathrm{t}}\stackrel{\mathrm{‾}}{c}{|}_{\text{bio}}$, are a limiting nutrient N (here expressed in nitrogen units, akin to nitrate); a phytoplankton tracer, P; a zooplankton tracer, Z; and a detrital pool, D. The non-conservative terms, ${\partial }_{\mathrm{t}}\stackrel{\mathrm{‾}}{c}{|}_{\text{nct}}$, represent physical sources and sinks of tracers, including surface fluxes, restoring at the offshore boundary, and optional restoring throughout the domain.

2.1 Tracer evolution

We first formulate an evolution equation for the meridionally averaged concentration of an arbitrary tracer c. We assume that c evolves according to a combination of advection by the three-dimensional ocean flow and diffusion by microscale mixing processes:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial c}{\partial t}}{|}_{\text{phys}}=\underset{\text{advection}}{\underbrace{-{\mathrm{\nabla }}_{\mathrm{3}}\cdot \left({\mathit{u}}_{\mathrm{3}}c\right)}}+\underset{\text{mixing}}{\underbrace{{\mathrm{\nabla }}_{\mathrm{3}}\cdot \left({\mathit{\kappa }}_{\text{dia}}{\mathrm{\nabla }}_{\mathrm{3}}c\right)}}.\end{array}$$

Here, u₃ is the three-dimensional velocity vector, ∇₃ is the three-dimensional gradient operator, and κ_dia the microscale diffusivity. In Eq. (2), we have assumed that the velocity field is non-divergent, i.e., ${\mathrm{\nabla }}_{\mathrm{3}}\cdot {\mathit{u}}_{\mathrm{3}}=\mathrm{0}$. We further assume that u₃ and c have already been averaged over a short timescale to exclude fluctuations associated with microscale eddies, whose effects are parameterized via the microscale mixing term (e.g., Aiki and Richards, 2008). We further simplify Eq. (2) by assuming that horizontal tracer gradients are small compared with vertical gradients, i.e., ${\partial }_{z}c\gg {\partial }_{x}c,{\partial }_{y}c$, as is typical for oceanic scales of evolution (e.g., Vallis, 2017). This implies that the microscale mixing acts primarily in the vertical, i.e.,

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial c}{\partial t}}{|}_{\text{phys}}\approx -{\mathrm{\nabla }}_{\mathrm{3}}\cdot \left({\mathit{u}}_{\mathrm{3}}c\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial c}{\partial z}}\right).\end{array}$$

We now reduce the dimensionality of Eq. (3) by taking a meridional average, which we denote via an overbar:

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{‾}}{•}={\displaystyle \frac{\mathrm{1}}{L_y}}\underset{\mathrm{0}}{\overset{L_y}{\int }}\cdot \mathrm{d}y.\end{array}$$

Here, L_y is the meridional length of the region of interest and y is the meridional coordinate. Though we refer to this average as “meridional” throughout the text, for the purpose of comparison with EBUSs in nature, this average might be thought of instead as an along-coast average or as an average following isobaths, under the assumption that the additional metric terms introduced by such coordinate transformations are negligible. We next perform a Reynolds decomposition of the velocity and tracer fields:

$$\begin{array}{}\text{(5a)}& {\displaystyle }\mathit{u}& {\displaystyle }=\stackrel{\mathrm{‾}}{\mathit{u}}+{\mathit{u}}^{\prime },\text{(5b)}& {\displaystyle }c& {\displaystyle }=\stackrel{\mathrm{‾}}{c}+c^{\prime },\end{array}$$

where primes ^′ denote perturbations from the meridional average. Taking a meridional average of Eq. (3) then yields

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}\right|}_{\text{phys}}& =\underset{\text{mean advection}}{\underbrace{-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{\mathit{u}}\stackrel{\mathrm{‾}}{c}\right)}}\\ & \underset{\text{eddy flux}}{\underbrace{-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\right)}}\underset{\text{meridional advection}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{L_y}}{\left[vc\right]}_{\mathrm{0}}^{L_y}}}+\underset{\text{mixing}}{\underbrace{{\displaystyle \frac{\partial }{\partial z}}\left(\stackrel{\mathrm{‾}}{{\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial c}{\partial z}}}\right)}}.\end{array}\end{array}$$

Here, we have used Eq. (5a)–(5b) and the property that perturbations vanish under the average, i.e., $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }}=\stackrel{\mathrm{‾}}{c^{\prime }}=\mathrm{0}$. We further define $\mathrm{\nabla }\equiv {\partial }_x\widehat{\mathit{x}}+{\partial }_z\widehat{\mathit{z}}$ as the zonal–vertical gradient operator, and $\mathit{u}=u\widehat{\mathit{x}}+w\widehat{\mathit{z}}$ as the zonal–vertical velocity vector. The square bracket indicates the difference between vc at the northern and southern boundaries of the domain of integration, i.e.,

$$\begin{array}{}\text{(7)}& {\left[vc\right]}_{\mathrm{0}}^{L_y}={\left.vc\right|}_{y=L_y}-{\left.vc\right|}_{y=\mathrm{0}}.\end{array}$$

In its current form, Eq. (6) cannot be solved prognostically for $\stackrel{\mathrm{‾}}{c}$ because it includes correlations between perturbation quantities, i.e., the eddy tracer flux $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}$. Assuming that these perturbations are associated with mesoscale eddies, we parameterize the eddy tracer flux following Gent and McWilliams (1990) and Redi (1982). Specifically, we decompose the eddy tracer flux into advection of the mean tracer $\stackrel{\mathrm{‾}}{c}$ by “eddy-induced velocity” u^⋆ and diffusion of $\stackrel{\mathrm{‾}}{c}$ along the mean buoyancy surfaces (see Burke et al., 2015):

$$\begin{array}{}\text{(8)}& \mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\right)=\mathrm{\nabla }\cdot \left({\mathit{u}}^{\star }\stackrel{\mathrm{‾}}{c}\right)-\mathrm{\nabla }\cdot \left({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}\right).\end{array}$$

Here, ∇_∥ denotes the gradient along mean buoyancy surfaces (see Sect. 2.2. A more detailed derivation of Eq. (8) is given in Appendix A. We additionally simplify the meridional tracer advection term by assuming that $\stackrel{\mathrm{‾}}{\partial v/\partial y}\approx \mathrm{0}$, i.e., that the meridional tracer flux convergence is dominated by meridional tracer gradients, and that correlations between κ_dia and c are negligible, i.e., that the meridionally averaged vertical diffusive tracer flux serves to diffuse $\stackrel{\mathrm{‾}}{c}$ downgradient. With these simplifications, the full equation for the physical evolution of tracers is given by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}\right|}_{\text{phys}}& =\underset{\text{mean advection}}{\underbrace{-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{\mathit{u}}\stackrel{\mathrm{‾}}{c}\right)}}\underset{\text{meridional advection}}{\underbrace{-{\displaystyle \frac{\stackrel{\mathrm{‾}}{v}}{L_y}}{\left[c\right]}_{\mathrm{0}}^{L_y}}}\\ & \underset{\text{eddy advection}}{\underbrace{-\mathrm{\nabla }\cdot \left({\mathit{u}}^{\star }\stackrel{\mathrm{‾}}{c}\right)}}\underset{\text{eddy stirring}}{\underbrace{-\mathrm{\nabla }\cdot \left({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}\right)}}+\underset{\text{mixing}}{\underbrace{{\partial }_z\left({\mathit{\kappa }}_{\text{dia}}{\partial }_z\stackrel{\mathrm{‾}}{c}\right)}}.\end{array}\end{array}$$

The terms on the right-hand side of Eq. (9) are discussed further in the following sections: in Sect. 2.2, we discuss the evolution of the mean velocity $\stackrel{\mathrm{‾}}{\mathit{u}}$ via the momentum equations, and in Sect. 2.3 we discuss the subgrid-scale parameterizations, i.e., eddy advection, eddy stirring, and mixing.

2.2 Momentum evolution equations

To evolve a meridionally averaged tracer $\stackrel{\mathrm{‾}}{c}$ using Eq. (9), the meridionally averaged velocity field ${\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{3}}$ is required. This velocity field is evolved in MAMEBUS by solving a simplified form of the hydrostatic Boussinesq momentum and continuity equations with a linear equation of state (Vallis, 2017):

$$\begin{array}{}\text{(10a)}& {\displaystyle \frac{\partial u}{\partial t}}& {\displaystyle }=-{\mathit{u}}_{\mathrm{3}}\cdot {\mathrm{\nabla }}_{\mathrm{3}}u+fv-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial u}{\partial z}}\right),\text{(10b)}& {\displaystyle \frac{\partial v}{\partial t}}& {\displaystyle }=-{\mathit{u}}_{\mathrm{3}}\cdot {\mathrm{\nabla }}_{\mathrm{3}}v-fu-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial y}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial v}{\partial z}}\right),\text{(10c)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}& {\displaystyle }=b,\text{(10d)}& {\displaystyle }{\mathrm{\nabla }}_{\mathrm{3}}\cdot {\mathit{u}}_{\mathrm{3}}& {\displaystyle }=\mathrm{0},\text{(10e)}& {\displaystyle }b& {\displaystyle }=g\mathit{\alpha }\mathit{\theta }.\end{array}$$

Here, $\mathit{\varphi }=p/{\mathit{\rho }}_{\mathrm{0}}$ is the dynamic pressure, where ρ₀ is an arbitrary reference density, b is the buoyancy, θ is the potential temperature, α is the thermal expansion coefficient (assumed constant), g is the gravitational constant, and f is the Coriolis parameter. Note that we have assumed that momentum is mixed by microscale turbulence following the same diffusivity κ_dia as tracers (see Sect. 2.1), i.e., that the turbulent Prandtl number (e.g., Kays, 1994) is exactly equal to 1.

As in Sect. 2.1, we now meridionally average Eq. (10a)–(10e) to obtain evolution equations for $\stackrel{\mathrm{‾}}{u}$ and $\stackrel{\mathrm{‾}}{v}$, and thus implicitly also for $\stackrel{\mathrm{‾}}{w}$. This yields the following set of averaged equations:

$$\begin{array}{}\text{(11a)}& {\displaystyle \frac{\partial \overline{u}}{\partial t}}& {\displaystyle }=f\stackrel{\mathrm{‾}}{v}-{\displaystyle \frac{\partial \overline{\mathit{\varphi }}}{\partial x}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial \overline{u}}{\partial z}}\right),\text{(11b)}& {\displaystyle \frac{\partial \overline{v}}{\partial t}}& {\displaystyle }=-f\stackrel{\mathrm{‾}}{u}-{\displaystyle \frac{\mathrm{1}}{L_y}}{\left[\mathit{\varphi }\right]}_{\mathrm{0}}^{L_y}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial \overline{v}}{\partial z}}\right),\text{(11c)}& {\displaystyle }\mathrm{\nabla }\cdot \overline{\mathit{u}}& {\displaystyle }=\mathrm{0},\text{(11d)}& {\displaystyle \frac{\partial \overline{\mathit{\varphi }}}{\partial z}}& {\displaystyle }=\mathit{u},\text{(11e)}& {\displaystyle }\stackrel{\mathrm{‾}}{b}& {\displaystyle }=g\mathit{\alpha }\stackrel{\mathrm{‾}}{\mathit{\theta }}.\end{array}$$

Here, we have made the frictional–geostrophic approximation (e.g., Edwards et al., 1998), assuming that the Rossby number of the flow is small (e.g., Vallis, 2017), and thus that momentum advection (second terms from the left in Eq. 10a–10b) is negligible compared to other terms in the momentum equation. This assumption may indeed have some limitations in upwelling regions with steep topography and strong stratification. Lentz and Chapman (2004) show that in the cross-shelf momentum flux divergence balances the wind stress and supports an on shore return flow, which can impact nitrate concentrations on the shelf (Jacox and Edwards, 2011).

On the other hand, we have retained the time-evolution terms (leftmost terms in Eq. 10a–10b) to allow forward evolution of the horizontal velocity fields; if these terms were neglected, then these terms would need to be computed diagnostically at each time step. The resulting system is almost identical to the T3W equations (Dauhajre and McWilliams, 2018), a time-varying extension of the turbulent thermal wind balance (Gula et al., 2014), which was developed to explain the circulation of submesoscale fronts. The meridional pressure gradient in Eq. (11a) is imposed, rather than solved for prognostically, and is assumed to be set by the larger-scale subtropical gyre circulation encompassing the EBUS, which explicitly differs from the work done in Dauhajre and McWilliams (2018) which focuses on more rapid time-varying evolution on smaller scales. Together with the tracer advection equation for potential temperature (i.e., Eq. 9 with c=θ), Eq. (11a)–(11e) comprise a closed set of equations for the physical evolution of MAMEBUS.

In Eq. (11c), we have invoked the earlier assumption that $\stackrel{\mathrm{‾}}{\partial v/\partial y}\approx \mathrm{0}$ (see Sect. 2.1), such that the averaged velocity field is non-divergent in the x–z plane. This implies that the zonal–vertical velocity field can be related to a mean stream function $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ via

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{u}=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial z}},\stackrel{\mathrm{‾}}{w}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial x}}.\end{array}$$

These relationships allow us to calculate $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ and thus $\stackrel{\mathrm{‾}}{w}$, from $\stackrel{\mathrm{‾}}{u}$, subject to the boundary conditions

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{‾}}{\mathit{\psi }}=\mathrm{0}\text{ at }z=\mathrm{0},z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right).\end{array}$$

Here, $z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)$ is the mean sea floor elevation.

Additional boundary conditions are required to solve Eq. (11a)–(11e) prognostically. Specifically, we require that the vertical turbulent stress in Eq. (11a)–(11b) matches the wind stress applied at the sea surface and the drag stress at the sea floor, with the latter formulated via a linear drag law. Formally, these boundary conditions are

$$\begin{array}{}\text{(14a)}& {\displaystyle }{\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u}}{\partial z}}=\mathrm{0},{\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v}}{\partial z}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\tau }}}^y}{{\mathit{\rho }}_{\mathrm{0}}}}\text{ at }z=\mathrm{0},\text{(14b)}& {\displaystyle }{\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u}}{\partial z}}=r\stackrel{\mathrm{‾}}{u},{\mathit{\kappa }}_{\text{dia}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v}}{\partial z}}=r\stackrel{\mathrm{‾}}{v}\text{ at }z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right).\end{array}$$

Here, r is a linear drag coefficient and ${\overline{\mathit{\tau }}}^y$ is the meridional wind stress.

2.3 Physical parameterizations

In this section, we describe the parameterization of unresolved microscale mixing in the tracer evolution (Eq. 9) and the horizontal momentum (Eq. 11a–11b), and of mesoscale eddy advection and stirring in Eq. (9). This amounts to parameterizing the diapycnal diffusivity κ_dia, the isopycnal diffusivity κ_iso, and the eddy velocity u^⋆.

2.3.1 Diapycnal mixing

We formulate the diapycnal mixing coefficient κ_dia as a sum of four distinct contributing processes: surface mixed layer turbulence (κ_sml), bottom boundary layer turbulence (κ_bbl), turbulence due to convective overturns within the water column (κ_conv), and background mixing due to internal wave breaking (κ_bg). Formally, we write

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathit{\kappa }}_{\text{dia}}(x,z,t)& ={\mathit{\kappa }}_{\text{sml}}(x,z)+{\mathit{\kappa }}_{\text{bbl}}(x,z)\\ & +{\mathit{\kappa }}_{\text{conv}}(x,z,t)+{\mathit{\kappa }}_{\text{bg}}(x,z).\end{array}\end{array}$$

The terms on the right-hand side of Eq. (15) are discussed in turn in the following paragraphs.

The diapycnal diffusivity in the surface mixed layer, κ_sml, is prescribed to have the same structure as that used in the κ-profile parameterization (KPP) of Large et al. (1994). However, for simplicity, the mixed layer depth H_sml(x) and maximum magnitude κ_sml(x) are prescribed functions, rather than depending on the local surface forcing. The vertical profile of κ_dia in the surface mixed layer, i.e., $-H_{\text{sml}}<z<\mathrm{0}$, is given by

$$\begin{array}{}\text{(16)}& {\mathit{\kappa }}_{\text{sml}}(x,z)={\mathit{\kappa }}_{\text{sml}}^{\mathrm{0}}{G}_{\text{KPP}}\left({\mathit{\sigma }}_{\text{sml}}\right),\end{array}$$

where the dimensionless surface mixed layer vertical coordinate ${\mathit{\sigma }}_{\text{sml}}=-z/H_{\text{sml}}$ is defined such that $\mathrm{0}\le {\mathit{\sigma }}_{\text{sml}}\le \mathrm{1}$ within the mixed layer. The structure function G_KPP(σ_sml) is given by

$$\begin{array}{}\text{(17)}& G_{\text{KPP}}\left(\mathit{\sigma }\right)=\left\{\begin{array}{rlr}{\displaystyle }& {\displaystyle \frac{\mathrm{27}}{\mathrm{4}}}{\mathit{\sigma }}_{\text{sml}}(\mathrm{1}-{\mathit{\sigma }}_{\text{sml}}{)}^{\mathrm{2}},& {\displaystyle }\mathrm{0}\le {\mathit{\sigma }}_{\text{sml}}\le \mathrm{1},\\ {\displaystyle }& {\displaystyle }\mathrm{0},& {\displaystyle }{\mathit{\sigma }}_{\text{sml}}\ge \mathrm{1},\end{array}\right.\end{array}$$

following Large et al. (1994) and Troen and Mahrt (1986). The scaling factor of 27∕4 ensures that G_KPP(σ_sml) has a maximum of 1 for $\mathrm{0}<{\mathit{\sigma }}_{\text{sml}}<\mathrm{1}$.

The diapycnal diffusivity in the bottom boundary layer, κ_bbl, is prescribed in the same way as κ_sml but over the depth range ${\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}<z<{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}+H_{\text{bbl}}\left(x\right)$. Thus, analogous to Eq. (16), we prescribe

$$\begin{array}{}\text{(18)}& {\mathit{\kappa }}_{\text{bbl}}(x,z)={\mathit{\kappa }}_{\text{bbl}}^{\mathrm{0}}{G}_{\text{KPP}}\left({\mathit{\sigma }}_{\text{bbl}}\right),\end{array}$$

where the dimensionless bottom boundary layer vertical coordinate is defined as ${\mathit{\sigma }}_{\text{bbl}}=(z-{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}})/H_{\text{bbl}}$.

At any point in space and time at which the water column is statically unstable, i.e., when N²<0, we increase the value of κ_dia is increased locally to parameterize the effect of density-driven convection. That is, we prescribe κ_conv following

$$\begin{array}{}\text{(19)}& {\mathit{\kappa }}_{\text{conv}}=\left\{\begin{array}{rlr}{\displaystyle }& {\displaystyle }{\mathit{\kappa }}_{\text{conv}}^{\mathrm{0}},& {\displaystyle }{N}^{\mathrm{2}}<\mathrm{0},\\ {\displaystyle }& {\displaystyle }\mathrm{0},& {\displaystyle }{N}^{\mathrm{2}}\ge \mathrm{0}.\end{array}\right.\end{array}$$

Finally, the background diapycnal mixing, κ_bg(x,z), is simply prescribed as a constant background diffusivity. There are others that can be used (e.g., St. Laurent et al., 2002), but we opt for simplicity in the first version of this model.

2.3.2 Eddy advection and isopycnal mixing

We now discuss the formulation of the eddy advection and isopycnal mixing terms in Eq. (9). As discussed in Sect. 2.1, we follow the assumptions and formalism of the Gent and McWilliams (1990) and Redi (1982) parameterizations, which are commonly used in ocean models that do not explicitly resolve mesoscale eddies (e.g., Gent, 2011). These parameterizations assume that eddy-induced fluxes of buoyancy and tracer diffusion are directed along isopycnal slopes and so must be augmented in the ocean's surface mixed layer (SML) and bottom boundary layer (BBL). Here, the isopycnal slopes become very steep and isopycnals intersect the sea surface and floor (Tréguier et al., 1997). MAMEBUS therefore uses a modified form of the Ferrari et al. (2008) boundary layer parameterization, in which eddy buoyancy and tracer fluxes are rotated through the SML and BBL in order to enforce vanishing eddy-induced mass and tracer fluxes through the boundaries. Here, we summarize salient properties of this scheme, and in Appendix C we highlight differences between our scheme and that of Ferrari et al. (2008).

The eddy-induced velocity ${\mathit{u}}^{\star }=(u^{\star },w^{\star })$, introduced in Eq. (8), is non-divergent by construction (see Appendix A) and so we write it as

$$\begin{array}{}\text{(20)}& u^{\star }=-{\displaystyle \frac{\partial {\mathit{\psi }}^{\star }}{\partial z}},w^{\star }={\displaystyle \frac{\partial {\mathit{\psi }}^{\star }}{\partial x}},\end{array}$$

where ψ^⋆ is the “eddy stream function”. This advecting stream function is assumed to be the same for all tracers, which is accurate in the limit of small-amplitude fluctuations of the velocity and tracer fields (Plumb, 1979), and takes the form

$$\begin{array}{}\text{(21)}& {\mathit{\psi }}^{\star }={\mathit{\kappa }}_{\text{gm}}{S}_{\text{gm}}.\end{array}$$

Here, κ_gm is the Gent–McWilliams diffusivity and the S_gm is the is the Gent–McWilliams slope. The latter is conventionally set equal to the mean isopycnal slope (Gent and McWilliams, 1990):

$$\begin{array}{}\text{(22)}& S_{\text{int}}=-{\partial }_x\stackrel{\mathrm{‾}}{b}/{\partial }_z\stackrel{\mathrm{‾}}{b}.\end{array}$$

However, we allow S_gm to diverge from S_int in the SML and BBL, in part to ensure that the no-flux surface and bottom boundary conditions are satisfied (Ferrari et al., 2008):

$$\begin{array}{}\text{(23)}& {\mathit{\psi }}^{\star }=\mathrm{0}\text{ at }z=\mathrm{0},z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right).\end{array}$$

Specifically, we prescribe

$$\begin{array}{}\text{(24)}& S_{\text{gm}}=\left\{\begin{array}{ll}{S}_{\text{sml}},& -H_{\text{sml}}\left(x\right)<z<\mathrm{0},\\ S_{\text{int}},& {\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)+H_{\text{bbl}}\left(x\right)<z<-H_{\text{sml}}\left(x\right),\\ S_{\text{bbl}},& {\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)<z<{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)+H_{\text{bbl}}\left(x\right).\end{array}\right.\end{array}$$

The formulations of the modified slopes S_sml and S_bbl are discussed below in the following sections (“surface mixed layer” and “bottom boundary layer” subsections).

The isopycnal mixing operator serves to mix tracers down their mean gradients, in a direction that is parallel to mean isopycnal surfaces in the ocean interior, following (Redi, 1982). This may be written componentwise as

$$\begin{array}{}\text{(25)}& \begin{array}{rl}\mathrm{\nabla }\cdot \left({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}\right)& ={\displaystyle \frac{\partial }{\partial x}}\left({\mathit{\kappa }}_{\text{iso}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}+{\mathit{\kappa }}_{\text{iso}}{S}_{\text{iso}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\text{iso}}{S}_{\text{iso}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}+{\mathit{\kappa }}_{\text{iso}}{S}_{\text{iso}}^{\mathrm{2}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}\right),\end{array}\end{array}$$

where S_iso denotes the slope of the surface along which the tracer is to be mixed and is assumed to be small (S_iso≪1). Similar to S_gm, this slope is conventionally set equal to the mean isopycnal slope S_int, but we apply modifications to the formulation of S_iso in the SML and BBL to ensure that there is zero eddy-induced tracer flux through the domain boundaries, i.e.,

$$\begin{array}{}\text{(26)}& {\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}\cdot \widehat{\mathit{n}}=\mathrm{0}\text{ at }z=\mathrm{0},z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right),\end{array}$$

where $\widehat{\mathit{n}}$ is a unit vector oriented perpendicular to the sea surface or sea floor. Specifically, we prescribe

$$\begin{array}{}\text{(27)}& S_{\text{iso}}=\left\{\begin{array}{ll}{S}_{\text{sml}},& -H_{\text{sml}}\left(x\right)<z<\mathrm{0},\\ S_{\text{int}},& {\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)+H_{\text{bbl}}\left(x\right)<z<-H_{\text{sml}},\\ {\stackrel{\mathrm{̃}}{S}}_{\text{bbl}},& {\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)<z<{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)+H_{\text{bbl}}\left(x\right).\end{array}\right.\end{array}$$

Thus, S_gm and S_iso are identical everywhere above the BBL. The need for a distinction within the BBL is explained below in the “surface mixed layer” and “bottom boundary layer” subsections.

Surface mixed layer

We now discuss the formulation of S_sml, the effective isopycnal slope in the surface mixed layer. Following Ferrari et al. (2008), we construct S_sml in a way that avoids singularities due to the vanishingly small vertical buoyancy gradients, and thus near-infinite isopycnal slopes, that occur in the mixed layer. This is achieved by using the vertical buoyancy gradient at the base of the mixed layer to define the effective slope as

$$\begin{array}{}\text{(28)}& S_{\text{sml}}=-G_{\text{sml}}\left({\mathit{\sigma }}_{\text{sml}}\right){\displaystyle \frac{{\displaystyle }{\partial }_x\stackrel{\mathrm{‾}}{b}}{{\displaystyle }{\partial }_z\stackrel{\mathrm{‾}}{b}{|}_{z=-H_{\text{sml}}}}},\end{array}$$

where ${\mathit{\sigma }}_{\text{sml}}=-z/H_{\text{sml}}$ is a dimensionless vertical coordinate for the SML, as in Sect. 2.3.1. The corresponding eddy stream function (Eq. 21) is identical to that of Ferrari et al. (2008):

$$\begin{array}{}\text{(29)}& {\mathit{\psi }}^{\star }=-{\mathit{\kappa }}_{\text{gm}}{G}_{\text{sml}}\left({\mathit{\sigma }}_{\text{sml}}\right){\displaystyle \frac{{\displaystyle }{\partial }_x\stackrel{\mathrm{‾}}{b}}{{\displaystyle }{\partial }_z\stackrel{\mathrm{‾}}{b}{|}_{z=-H_{\text{sml}}}}},z\ge -H_{\text{sml}}.\end{array}$$

The structure function G_sml(z) is required to enforce continuity of the vertical tracer fluxes and flux divergences at the surface and at the base of the mixed layer. For example, Eq. (23) requires that G_sml vanishes at the surface:

$$\begin{array}{}\text{(30)}& G_{\text{sml}}\left(\mathrm{0}\right)=\mathrm{0}.\end{array}$$

We further require that the eddy stream function and eddy residual tracer fluxes be continuous at the base of the SML, i.e., that S_sml=S_int, which requires that

$$\begin{array}{}\text{(31)}& G_{\text{sml}}\left(\mathrm{1}\right)=\mathrm{1}.\end{array}$$

Finally, we require continuity of the divergence of the eddy tracer flux in order to avoid producing singularities at the SML base. The zonal and vertical components of the eddy tracer flux are

$$\begin{array}{}\text{(32a)}& {\displaystyle }\stackrel{\mathrm{‾}}{u^{\prime }{c}^{\prime }}& {\displaystyle }={\mathit{\kappa }}_{\text{gm}}{S}_{\text{gm}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}-{\mathit{\kappa }}_{\text{iso}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}+S_{\text{iso}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}\right),\text{(32b)}& {\displaystyle }\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}& {\displaystyle }=-{\mathit{\kappa }}_{\text{gm}}{S}_{\text{gm}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}-{\mathit{\kappa }}_{\text{iso}}{S}_{\text{iso}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}+S_{\text{iso}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}\right).\end{array}$$

It may be shown that continuity of $\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}$ across $z=-H_{\text{sml}}$ is guaranteed if

$$\begin{array}{}\text{(33)}& {\left.{\displaystyle \frac{\partial S_{\text{sml}}}{\partial z}}\right|}_{z=-H_{\text{sml}}^{+}}={\left.{\displaystyle \frac{\partial S_{\text{int}}}{\partial z}}\right|}_{z=-H_{\text{sml}}^{-}}\Rightarrow G_{\text{sml}}^{\prime }\left(\mathrm{1}\right)={\displaystyle \frac{H_{\text{sml}}}{{\mathit{\lambda }}_{\text{sml}}}},\end{array}$$

where ${\mathit{\lambda }}_{\text{sml}}={\partial }_{zz}\stackrel{\mathrm{‾}}{b}/{\partial }_z\stackrel{\mathrm{‾}}{b}{|}_{z=-H_{\text{sml}}}$ is a vertical length scale for eddy motions at the base of the mixed layer.

The simplest form for G_sml(z) that satisfies conditions for Eqs. (30), (31), and (33) is a quadratic function of depth:

$$\begin{array}{}\text{(34)}& G_{\text{sml}}\left({\mathit{\sigma }}_{\text{sml}}\right)=-\left(\mathrm{1}-{\displaystyle \frac{H_{\text{sml}}}{{\mathit{\lambda }}_{\text{sml}}}}\right){\mathit{\sigma }}_{\text{sml}}^{\mathrm{2}}+\left(\mathrm{2}-{\displaystyle \frac{H_{\text{sml}}}{{\mathit{\lambda }}_{\text{sml}}}}\right){\mathit{\sigma }}_{\text{sml}}.\end{array}$$

Equation (34) is currently implemented in MAMEBUS. A more sophisticated form of G_sml that arguably has stronger physical motivation is given by Ferrari et al. (2008). They split the SML into a true mixed layer, in which G_sml varies linearly (and so the eddy velocity is approximately uniform), overlying a transition layer, in which G_sml(σ_sml) varies quadratically.

Bottom boundary layer

The scheme described above for the SML relies on the fact that the ocean surface is approximately flat, which allows the same effective slopes S_sml to be used for S_gm and S_iso. The sloping sea floor requires separate BBL slopes, S_bbl and ${\stackrel{\mathrm{̃}}{S}}_{\text{bbl}}$, and structure functions, G_bbl and ${\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}$, to satisfy the required conditions of no volume nor tracer flux through the boundary, i.e., Eqs. (23) and (26).

Analogous to the SML, we define the effective slope S_bbl as

$$\begin{array}{}\text{(35)}& S_{\text{bbl}}=-G_{\text{bbl}}\left({\mathit{\sigma }}_{\text{bbl}}\right){\displaystyle \frac{{\displaystyle }{\partial }_x\stackrel{\mathrm{‾}}{b}}{{\displaystyle }{\partial }_z\stackrel{\mathrm{‾}}{b}{|}_{z={\mathit{\eta }}_{\mathrm{b}}+H_{\text{bbl}}}}},\end{array}$$

where ${\mathit{\sigma }}_{\text{bbl}}=(z-{\mathit{\eta }}_{\mathrm{b}}(x\left)\right)/H_{\text{bbl}}\left(x\right)$ is the BBL vertical coordinate, as in Sect. 2.3.1. The eddy stream function in the BBL is therefore

$$\begin{array}{}\text{(36)}& {\mathit{\psi }}^{\star }=-{\mathit{\kappa }}_{\text{gm}}{G}_{\text{bbl}}\left({\mathit{\sigma }}_{\text{bbl}}\right){\displaystyle \frac{{\displaystyle }{\partial }_x\stackrel{\mathrm{‾}}{b}}{{\displaystyle }{\partial }_z\stackrel{\mathrm{‾}}{b}{|}_{z={\mathit{\eta }}_{\mathrm{b}}+H_{\text{bbl}}}}},z\le {\mathit{\eta }}_{\mathrm{b}}+H_{\text{bbl}}.\end{array}$$

To satisfy the condition of zero volume flux through the sea floor, Eq. (23), the effective slope must vanish at z=η_b(x), which requires

$$\begin{array}{}\text{(37)}& G_{\text{bbl}}\left(\mathrm{0}\right)=\mathrm{0}.\end{array}$$

To ensure continuity of the eddy stream function at the top of the BBL, we require that S_bbl approach S_int, i.e.,

$$\begin{array}{}\text{(38)}& G_{\text{bbl}}\left(\mathrm{1}\right)=\mathrm{1}.\end{array}$$

Finally, to ensure continuity of the eddy bolus velocity, we require that the gradient of S_gm be continuous at $z={\mathit{\eta }}_{\mathrm{b}}+H_{\text{bbl}}$. This imposes a constraint analogous to Eq. (33) on G_bbl:

$$\begin{array}{}\text{(39)}& G_{\text{bbl}}^{\prime }\left(\mathrm{1}\right)=-{\displaystyle \frac{H_{\text{bbl}}}{{\mathit{\lambda }}_{\text{bbl}}}},\end{array}$$

where ${\left.{\mathit{\lambda }}_{\text{bbl}}={\stackrel{\mathrm{‾}}{b}}_{zz}/{\stackrel{\mathrm{‾}}{b}}_z\right|}_{z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}+H_{\text{bbl}}}$ is a vertical length scale for eddies at the top of the BBL. To satisfy Eqs. (37)–(39), we select a quadratic form for the structure function G_bbl(σ_bbl):

$$\begin{array}{}\text{(40)}& G_{\text{bbl}}\left({\mathit{\sigma }}_{\text{bbl}}\right)=-\left(\mathrm{1}+{\displaystyle \frac{H_{\text{bbl}}}{{\mathit{\lambda }}_{\text{bbl}}}}\right){\mathit{\sigma }}_{\text{bbl}}^{\mathrm{2}}+\left(\mathrm{2}+{\displaystyle \frac{H_{\text{bbl}}}{{\mathit{\lambda }}_{\text{bbl}}}}\right){\mathit{\sigma }}_{\text{bbl}}.\end{array}$$

However, the effective slope S_bbl can no longer be used to define S_iso in the BBL: Eq. (26) requires that the effective slope be aligned with the bottom slope at the sea floor, $S_{\mathrm{b}}={\partial }_x{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}$ at $z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}$. We must therefore employ a modified effective slope ${\stackrel{\mathrm{̃}}{S}}_{\text{bbl}}$ in the isopycnal mixing operator, as expressed in Eq. (27). We define ${\stackrel{\mathrm{̃}}{S}}_{\text{bbl}}$ as

$$\begin{array}{}\text{(41)}& {\stackrel{\mathrm{̃}}{S}}_{\text{bbl}}=S_{\text{bbl}}+\left(\mathrm{1}-{\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}\left(z\right)\right)S_{\mathrm{b}},\end{array}$$

where ${\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}\left({\mathit{\sigma }}_{\text{bbl}}\right)$ is a modified structure function that also vanishes at the ocean bed:

$$\begin{array}{}\text{(42)}& {\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}\left(\mathrm{0}\right)=\mathrm{0}.\end{array}$$

Continuity of the eddy tracer fluxes at the top of the BBL requires that

$$\begin{array}{}\text{(43)}& {\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}\left(\mathrm{1}\right)=\mathrm{1}.\end{array}$$

Finally, continuity of the eddy flux divergence is enforced by

$$\begin{array}{}\text{(44)}& {\left.{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}}{\partial z}}\right|}_{z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}+H_{\text{bbl}}}=\mathrm{0}.\end{array}$$

To satisfy Eqs. (42)–(44), we select a quadratic form for the structure function ${\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}\left({\mathit{\sigma }}_{\text{bbl}}\right)$:

$$\begin{array}{}\text{(45)}& {\stackrel{\mathrm{̃}}{G}}_{\text{bbl}}\left({\mathit{\sigma }}_{\text{bbl}}\right)={\mathit{\sigma }}_{\text{bbl}}(\mathrm{2}-{\mathit{\sigma }}_{\text{bbl}}).\end{array}$$

2.4 Biogeochemical model formulation

The current biogeochemical model implemented in MAMEBUS is an NPZD (nutrient–phytoplankton–zooplankton–detritus) model. This NPZD model is modeled after the size-structured AstroCAT (Banas, 2011) and Darwin models (Ward et al., 2012). For the purpose of this paper, we reduced the size-structured ecosystem model to single phytoplankton and zooplankton size classes, while preserving the option to run multiple size classes in future versions of the model. We also include a detritus variable, which allows for sinking and export of organic matter away from the euphotic zone and redistribution of nutrients in the water column.

The biogeochemical equations in MAMEBUS are formulated similarly to previous NPZD models but cast in terms of the meridionally averaged nutrient, phytoplankton, zooplankton, and detritus concentrations. We neglect additional terms that would be introduced by first formulating the equations and then taking the meridional average, e.g., covariances of the type $\stackrel{\mathrm{‾}}{P^{\prime }{Z}^{\prime }}$. This assumption is partially predicated on the idea that zonal gradients in biogeochemical tracers (e.g., nutrients and chlorophyll) are much stronger than meridional gradients, as supported by observations and models (Fiechter et al., 2018). For example, Venegas et al. (2008) show that average chlorophyll concentrations during the upwelling season vary approximately 2-fold in the northern California Current System, whereas observations from the California Cooperative Oceanic Fisheries Investigations (CalCOFI) (Fig. 7) show variations by an order of magnitude between nearshore and offshore stations. Along-shore gradients in chlorophyll are observed along the coast, where they are driven by wind and topographic variations; however, they are generally much smaller than the gradient between the coast and the offshore region (Fiechter et al., 2018). We recognize that this is a simplification of the true variability in EBUSs, but we consider it appropriate on average over the entire upwelling system, in particular within the idealized MAMEBUS framework, and plan to reassess it in future work.

We drop the bar notation indicating a meridional average for this section, with the understanding that all variables denote meridionally averaged quantities. In the following, we include size-dependent uptake and grazing, along with variable sinking speeds for detritus, to retain essential size-dependent biogeochemical interactions and export fluxes. This will facilitate a future introduction of multiple size classes in the model. All variables and coefficients are given in Table 1. We note that all of the parameter values and equations described below measure time in days, whereas more generally MAMEBUS measures time in seconds; appropriate conversions are made in the model code to ensure dimensional consistency. The main conservation equations for biogeochemical tracers are

$$\begin{array}{}\text{(46a)}& {\displaystyle \frac{\partial N}{\partial t}}{|}_{\text{bio}}=-\mathcal{U}(N,I,T,P)+\mathcal{R}\left(D\right),\text{(46b)}& {\displaystyle \frac{\partial P}{\partial t}}{|}_{\text{bio}}=\mathcal{U}(N,I,T,P)-\mathcal{G}(P,Z)-\mathcal{M}\left(P\right),\text{(46c)}& {\displaystyle \frac{\partial Z}{\partial t}}{|}_{\text{bio}}=\mathit{\lambda }\mathcal{G}(P,Z)-\mathcal{M}\left(Z\right),\text{(46d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial D}{\partial t}}{|}_{\text{bio}}& =\mathcal{M}\left(P\right)+\mathcal{M}\left(Z\right)+(\mathrm{1}-\mathit{\lambda })\mathcal{G}(P,Z)\\ & -{\displaystyle \frac{\partial }{\partial z}}{w}_{\text{sink}}D-\mathcal{R}\left(D\right),\end{array}\end{array}$$

where T (^∘C) is the model temperature, I (W m⁻²) is the local irradiance profile, N (mmol N m⁻³) is nitrate concentration, P (mmol N m⁻³) is phytoplankton concentration, Z (mmol N m⁻³) is zooplankton concentration, and D (mmol N m⁻³) is the detritus concentration. The terms on the right-hand sides of Eq. (46a)–(46d) are explained in the following subsections.

Moore et al. (2001)Moore et al. (2001)Banas (2011)Moore et al. (2001)Banas (2011)Moore et al. (2001)Ward et al. (2012)Ward et al. (2012)Banas (2011)Edwards and Bees (2001)

Table 1 Parameters and values used in the ecosystem model implemented in MAMEBUS. Coefficients without explicit references are chosen by the user.

2.4.1 Nutrient uptake

Common controls on phytoplankton population are bottom-up limitation (i.e., nutrient control), and top-down grazing by zooplankton (Sarmiento and Gruber, 2006). We formulate bottom-up controls using typical choices for light- and temperature-dependent terms, and Michaelis–Menten uptake (Sarmiento and Gruber, 2006). The functional form of uptake is given by

$$\begin{array}{}\text{(47)}& \mathcal{U}(N,I,T,P)=\mathit{\phi }\left(I\right)\mathit{\phi }\left(T\right)U^{\text{max}}{\displaystyle \frac{N}{N+k_N}}P,\end{array}$$

where φ(I) and φ(T) are light- and temperature-limiting functions, respectively. The light attenuation is modeled by integrating the Beer–Lambert law, following Moore et al. (2001):

$$\begin{array}{}\text{(48a)}& {\displaystyle \frac{\partial I\left(z\right)}{\partial z}}& {\displaystyle }=-k_{\text{par}}I\left(z\right),\text{ where }{I}_{\mathrm{0}}=I(z=\mathrm{0})=Q_{\mathrm{sw}}{I}_{\mathrm{p}},\text{(48b)}& {\displaystyle }{k}_{\text{par}}& {\displaystyle }=k_{\mathrm{w}}+P\cdot k_{\mathrm{c}},\end{array}$$

and the light-dependent uptake function is modeled following Sarmiento and Gruber (2006),

$$\begin{array}{}\text{(49)}& \mathit{\phi }\left(I\right)={\displaystyle \frac{I\left(z\right)}{\sqrt{I_{\mathrm{0}}^{\mathrm{2}}+I(z{)}^{\mathrm{2}}}}}.\end{array}$$

The temperature component of the uptake function is

$$\begin{array}{}\text{(50)}& \mathit{\phi }\left(T\right)=e^{-r_T\left(T-T_{\mathrm{0}}\right)}.\end{array}$$

The maximum uptake rate is an allometric relationship defined as

$$\begin{array}{}\text{(51)}& U^{\text{max}}=a_{\mathrm{u}}{\left({\displaystyle \frac{{\mathrm{\ell }}_{\mathrm{p}}}{{\mathrm{\ell }}_{\mathrm{0}}}}\right)}^{b_{\mathrm{u}}},\end{array}$$

where ℓ_p is the user-determined phytoplankton size expressed as equivalent spherical diameter (ESD), and ℓ₀=1 µm is a normalized length scale, with all allometrically defined variables listed in Table 2. While there are other options for the bases of these allometric relationships outlined in this section (e.g., cell volume), we make the decision to use ESD as a measure of cell size. Finally, the half-saturation coefficient is k_N=0.1 mmol N m⁻³.

Tang (1995)Tang (1995)Hansen et al. (1994)Hansen et al. (1994)Hansen et al. (1994)Hansen et al. (1994)

Table 2 Parameters and values used in the ecosystem model implemented in MAMEBUS. Coefficients without explicit references are chosen by the user.

2.4.2 Grazing

Top-down processes are represented by zooplankton grazing on phytoplankton. Andersen et al. (2016) noted that there is an optimal length scale for active predation and grazing, as a strategic trade-off for optimal biomass assimilation. We make the assumption that the biomass assimilation of phytoplankton by zooplankton also follows Michaelis–Menten dynamics, then the functional form of grazing is given by

$$\begin{array}{}\text{(52)}& \mathcal{G}(P,Z)=G^{\text{max}}{\displaystyle \frac{\mathit{\vartheta }P}{k_P+\mathit{\vartheta }P}}Z,\end{array}$$

where the maximum grazing rate is defined by an allometric relationship defined as

$$\begin{array}{}\text{(53)}& G^{\text{max}}=a_{\mathrm{g}}{\left({\displaystyle \frac{{\mathrm{\ell }}_z}{{\mathrm{\ell }}_{\mathrm{0}}}}\right)}^{b_{\mathrm{g}}},\end{array}$$

where d⁻¹ represents a “per day” quantity. We define a Gaussian distribution about an optimal grazing length-scale following Banas (2011):

$$\begin{array}{}\text{(54)}& \mathit{\vartheta }=\exp \left(-{\displaystyle \frac{{\log }_{\mathrm{10}}\left({\mathrm{\ell }}_{\mathrm{p}}\right)-{\log }_{\mathrm{10}}\left({\mathrm{\ell }}_{\text{opt}}\right)}{\mathrm{\Delta }\mathrm{\ell }}}\right),\end{array}$$

where Δℓ sets the width of the optimal grazing profile and defines a band of grazing about the optimal prey size, ℓ_opt. By allowing for a variable band of grazing, we are able to control the assimilation efficiency of phytoplankton by zooplankton through direct preferential grazing. Accordingly, we model the optimal prey size based on a preferential grazing profile centered about an optimal predator–prey length scale:

$$\begin{array}{}\text{(55)}& {\mathrm{\ell }}_{\text{opt}}=a_{\mathrm{o}}{\left({\displaystyle \frac{{\mathrm{\ell }}_z}{{\mathrm{\ell }}_{\mathrm{0}}}}\right)}^{b_{\mathrm{o}}}.\end{array}$$

2.4.3 Mortality

Mortality closure terms often set important internal dynamics in ecosystem models (Poulin and Franks, 2010). While linear mortality terms are generally used for phytoplankton, zooplankton mortality is often modeled via a quadratic term to avoid unrealistic oscillations and stabilize the solution (Poulin and Franks, 2010). The quadratic mortality term may be rationalized as a representation of mixotrophic grazing, zooplankton self-grazing, and higher-order grazing in NPZD models (Raick et al., 2006). Therefore, we model phytoplankton mortality as

$$\begin{array}{}\text{(56)}& \mathcal{M}\left(P\right)={\mathit{\mu }}_p{U}^{\text{max}}P,\end{array}$$

and zooplankton mortality as

$$\begin{array}{}\text{(57)}& \mathcal{M}\left(Z\right)={\mathit{\mu }}_z{Z}^{\mathrm{2}}.\end{array}$$

2.4.4 Remineralization and particle sinking

Sinking particles are an essential component of the vertical transport of nutrients from the surface to the deep ocean (Sarmiento and Gruber, 2006). Once particles sink past the euphotic zone, they are remineralized and returned to the subsurface nutrient pool. In this model, we represent remineralization processes via a linear rate, i.e.,

$$\begin{array}{}\text{(58)}& \mathcal{R}\left(D\right)=r_{\text{remin}}D,\end{array}$$

where r_remin is the specific remineralization rate.

Particles sink at a constant average speed in the water column, following Eq. (46d). At the bottom boundary, we impose zero sinking flux, i.e., w_sink=0 at $z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)$. Thus, any nutrients that sink to the sea floor as detritus must remineralize there. This allows for redistribution of nutrients by mixing within the bottom boundary layer, diffusion into the interior, and transport via upwelling onto the shelf.

2.5 Non-conservative terms

In this section, we describe the treatment of all non-conservative terms in the tracer evolution equation. MAMEBUS allows arbitrary restoring of all tracers, which may be used, for example, to impose offshore boundary conditions or to impose restoring at the sea surface. Fixed fluxes of all tracers may also be imposed through the surface. More precisely, we formulate the non-conservative tracer tendency as

$$\begin{array}{}\text{(59)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{nct}}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{restore}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{flux}}.\end{array}$$

The restoring and surface flux components of this tendency are discussed separately below.

2.5.1 Restoring

The restoring of a tracer is represented as an exponential decay to a prescribed, spatially varying tracer field, ${\stackrel{\mathrm{‾}}{c}}_{\mathrm{r}}(x,z)$, with timescale t_r(x,z). The tracer restoring is then formulated as

$$\begin{array}{}\text{(60)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{restore}}=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{c}-{\stackrel{\mathrm{‾}}{c}}_r}{t_r}}.\end{array}$$

2.5.2 Tracer fluxes

Surface fluxes are represented as a tendency in the tracer concentration in the surface grid boxes. For an arbitrary tracer c, we formulate the surface flux term as follows:

$$\begin{array}{}\text{(61)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}{|}_{\text{flux}}={\displaystyle \frac{\partial F_{\text{flux}}^c}{\partial z}},F_{\text{flux}}=\left\{\begin{array}{ll}{F}_{\text{flux},\mathrm{0}}^c,& z=\mathrm{0},\\ \mathrm{0},& z<\mathrm{0}.\end{array}\right.\end{array}$$

Here, $F_{\text{flux},\mathrm{0}}^c$ is the downward flux of c (units of [c]m s⁻¹) at the surface. For the case of buoyancy, the surface flux is imposed as a surface energy flux, Q_s (W m⁻²), with

$$\begin{array}{}\text{(62)}& F_{\text{flux},\mathrm{0}}^b={\displaystyle \frac{g\mathit{\alpha }{Q}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{0}}{C}_{\mathrm{p}}}},\end{array}$$

where C_p=4000 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$ is the specific heat capacity.

3 MAMEBUSv1.0 algorithm

In this section, we discuss the numerical solution of the model equations presented in Sect. 2. This entails a recasting of the equations in terrain-following, or “sigma”, coordinates (e.g., Song and Haidvogel, 1994; Shchepetkin and McWilliams, 2003), followed by the spatial discretization of the equations and algorithms for numerical time stepping.

3.1 Formulation in terrain-following coordinates

We solve the model equations presented in Sect. 2 in a coordinate system that “stretches” in the vertical to follow the shape of the sea floor. Such a coordinate system avoids “steps” in the sea floor that arise, for example, when using geopotential vertical coordinates and allows fine vertical resolution of the bottom boundary layer (e.g., Song and Haidvogel, 1994; Shchepetkin and McWilliams, 2003). Formally, we make a coordinate transformation $(x,z)\to (x,\mathit{\sigma })$, where σ is a dimensionless vertical coordinate and is defined such that σ=0 at z=0 and $\mathit{\sigma }=-\mathrm{1}$ at $z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)$. This transformation requires a relationship between z and σ via a transformation function

$$\begin{array}{}\text{(63)}& z=\mathit{\zeta }(x,\mathit{\sigma }).\end{array}$$

For example, a “pure” σ coordinate corresponds to the choice

$$\begin{array}{}\text{(64)}& \mathit{\zeta }(x,\mathit{\sigma })=-\mathit{\sigma }{\stackrel{\mathrm{‾}}{h}}_{\mathrm{b}}\left(x\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{b}}\left(x\right)=-{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}\left(x\right)$ is the meridionally averaged water column thickness. However, this is not necessarily the most practical choice for numerical applications, in which it is useful to focus the vertical resolution over certain depth ranges (especially those close to the top and bottom boundaries of the ocean). MAMEBUS currently implements the UCLA-ROMS (Shchepetkin and McWilliams, 2005) transformation function:

$$\begin{array}{}\text{(65)}& \mathit{\zeta }(x,\mathit{\sigma })=h_{\mathrm{b}}\left(x\right)\left[{\displaystyle \frac{h_{\mathrm{c}}\mathit{\sigma }+h_{\mathrm{b}}\left(x\right)C\left(\mathit{\sigma }\right)}{h_{\mathrm{c}}+h_{\mathrm{b}}\left(x\right)}}\right].\end{array}$$

Here, C(σ) is the stretching function, defined as

$$\begin{array}{}\text{(66)}& C\left(\mathit{\sigma }\right)=\left\{\begin{array}{ll}{\displaystyle \frac{\exp \left({\mathit{\theta }}_{\mathrm{b}}\stackrel{\mathrm{̃}}{C}\left(\mathit{\sigma }\right)\right)-\mathrm{1}}{\mathrm{1}-\exp (-{\mathit{\theta }}_{\mathrm{b}})}},& {\mathit{\theta }}_{\mathrm{b}}>\mathrm{0},\\ \stackrel{\mathrm{̃}}{C}\left(\mathit{\sigma }\right),& {\mathit{\theta }}_{\mathrm{b}}\le \mathrm{0},\end{array}\right.\end{array}$$

where

$$\begin{array}{}\text{(67)}& \stackrel{\mathrm{̃}}{C}\left(\mathit{\sigma }\right)=\left\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}-\cosh \left({\mathit{\theta }}_{\mathrm{s}}\mathit{\sigma }\right)}{\cosh \left({\mathit{\theta }}_{\mathrm{s}}\right)-\mathrm{1}}},& {\mathit{\theta }}_{\mathrm{s}}>\mathrm{0},\\ -{\mathit{\sigma }}^{\mathrm{2}},& {\mathit{\theta }}_{\mathrm{s}}\le \mathrm{0}.\end{array}\right.\end{array}$$

Here, C and $\stackrel{\mathrm{̃}}{C}$ are the bottom and surface components of the stretching function, respectively. The parameters ${\mathit{\theta }}_{\mathrm{s}}\in [\mathrm{0},\mathrm{10}]$ and ${\mathit{\theta }}_{\mathrm{b}}\in [\mathrm{0},\mathrm{4}]$ are surface and bottom stretching parameters; larger values cause the near-surface and near-bottom portions of the domain to occupy larger fraction of σ space. The parameter h_c defines a surface layer thickness, in which the coordinate system is approximately aligned with geopotentials, provided that h_b≫h_c.

We now write the physical tracer evolution (Eq. 9) in σ coordinates. For a given function $f=f(x,z(x,\mathit{\sigma }\left)\right)$, we can write derivatives with respect to x and σ as

$$\begin{array}{}\text{(68a)}& {\displaystyle }{\left.{\displaystyle \frac{\partial f}{\partial x}}\right|}_{\mathit{\sigma }}& {\displaystyle }={\left.{\displaystyle \frac{\partial f}{\partial x}}\right|}_z+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}{\left.{\displaystyle \frac{\partial f}{\partial z}}\right|}_x,\text{(68b)}& {\displaystyle }{\left.{\displaystyle \frac{\partial f}{\partial \mathit{\sigma }}}\right|}_x& {\displaystyle }={\displaystyle \frac{\partial \mathit{\zeta }}{\partial \mathit{\sigma }}}{\left.{\displaystyle \frac{\partial f}{\partial z}}\right|}_x.\end{array}$$

Using these identities, we may write the divergence of an arbitrary vector F, with components F^(x) and F^(z) in the $\widehat{\mathit{x}}$ and $\widehat{\mathit{z}}$ directions, respectively, as

$$\begin{array}{}\text{(69)}& \begin{array}{rl}\mathrm{\nabla }\cdot \mathit{F}& ={\left.{\displaystyle \frac{\partial }{\partial x}}\right|}_z{F}^{\left(x\right)}+{\left.{\displaystyle \frac{\partial }{\partial z}}\right|}_x{F}^{\left(z\right)}\\ & ={\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial x}}\right|}_{\mathit{\sigma }}\left({\mathit{\zeta }}_{\mathit{\sigma }}{F}^{\left(x\right)}\right)\\ & +{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left(F^{\left(z\right)}-{\mathit{\zeta }}_x{F}^{\left(x\right)}\right).\end{array}\end{array}$$

Equation (69), combined with the definition of the mean stream function (Eq. 12), allows us to write the mean advection term in Eq. (9) as

$$\begin{array}{}\text{(70)}& \begin{array}{rl}\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{\mathit{u}}\stackrel{\mathrm{‾}}{c}\right)& =\mathrm{\nabla }\cdot \left(-\stackrel{\mathrm{‾}}{c}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial z}}\right|}_x,\stackrel{\mathrm{‾}}{c}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial x}}\right|}_z\right)\\ & ={\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial x}}\right|}_{\mathit{\sigma }}\left(-\stackrel{\mathrm{‾}}{c}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial \mathit{\sigma }}}\right|}_x\right)\\ & +{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left(\stackrel{\mathrm{‾}}{c}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial x}}\right|}_{\mathit{\sigma }}\right).\end{array}\end{array}$$

An analogous expression may be obtained for the eddy advection term, $\mathrm{\nabla }\cdot \left({\mathit{u}}^{\star }\stackrel{\mathrm{‾}}{c}\right)$, in Eq. (9), using the definition (Eq. 21) of the eddy stream function. Next, we apply Eq. (69) to the isopycnal mixing operator, defined by Eq. (25), in Eq. (9) to obtain

$$\begin{array}{}\text{(71)}& \begin{array}{rl}\mathrm{\nabla }\cdot \left({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}\right)& ={\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial x}}\right|}_{\mathit{\sigma }}\left[{\mathit{\zeta }}_{\mathit{\sigma }}{\mathit{\kappa }}_{\text{iso}}\left({\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}\right|}_{\mathit{\sigma }}\right.\right.\\ & \left.\left.+(S_{\text{iso}}-S_{\mathit{\sigma }}){\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial \mathit{\sigma }}}\right|}_x\right)\right]\\ & +{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left[{\mathit{\kappa }}_{\text{iso}}(S_{\text{iso}}-S_{\mathit{\sigma }})\left({\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}\right|}_{\mathit{\sigma }}\right.\right.\\ & \left.\left.+(S_{\text{iso}}-S_{\mathit{\sigma }}){\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial \mathit{\sigma }}}\right|}_x\right)\right],\end{array}\end{array}$$

where S_σ=ζ_x is the slope of surfaces of constant σ in x–z space, i.e., the slope of the σ coordinate grid lines. Thus, the isopycnal mixing operator is essentially just modified by subtracting S_σ from S_iso to obtain the mixing slope relative to the slope of the σ coordinate grid. Over most of the water column, S_iso is equal to the isopycnal slope S_int, given by

$$\begin{array}{}\text{(72)}& \begin{array}{rl}{S}_{\text{int}}& =-{\displaystyle \frac{{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial x}}\right|}_z}{{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial z}}\right|}_x}}=-{\displaystyle \frac{{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial x}}\right|}_{\mathit{\sigma }}-S_{\mathit{\sigma }}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial z}}\right|}_x}{{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial z}}\right|}_x}}\\ & =-{\displaystyle \frac{{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial x}}\right|}_{\mathit{\sigma }}}{{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial \mathit{\sigma }}}\right|}_x}}+S_{\mathit{\sigma }}.\end{array}\end{array}$$

Thus, the quantity S_int−S_σ can actually be computed more directly than the true isopycnal slope, as

$$\begin{array}{}\text{(73)}& S_{\text{int}}-S_{\mathit{\sigma }}=-{\displaystyle \frac{{\displaystyle }{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial x}}\right|}_{\mathit{\sigma }}}{{\displaystyle }{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{b}}{\partial \mathit{\sigma }}}\right|}_x}}.\end{array}$$

Finally, the σ coordinate transformation of the vertical (quasi-diapycnal) mixing operator is

$$\begin{array}{}\text{(74)}& {\left.{\displaystyle \frac{\partial }{\partial z}}\right|}_x\left({\mathit{\kappa }}_{\text{dia}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial z}}\right|}_x\right)={\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left({\mathit{\kappa }}_{\text{dia}}{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial \mathit{\sigma }}}\right|}_x\right).\end{array}$$

To summarize, we write Eq. (9) in σ coordinates as

$$\begin{array}{}\text{(75)}& {\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}\right|}_{\text{phys}}=G_{\text{adv}}+G_{\text{iso}}+G_{\text{dia}}+G_{\text{lat}},\end{array}$$

where the tendency terms are

$$\begin{array}{}\text{(76a)}& {\displaystyle }\begin{array}{rl}{G}_{\text{adv}}& =-{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial x}}\right|}_{\mathit{\sigma }}\left[{\mathit{\zeta }}_{\mathit{\sigma }}\left(-\stackrel{\mathrm{‾}}{c}{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial {\mathit{\psi }}^{†}}{\partial \mathit{\sigma }}}\right|}_x\right)\right]\\ & -{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left(\stackrel{\mathrm{‾}}{c}{\left.{\displaystyle \frac{\partial {\mathit{\psi }}^{†}}{\partial x}}\right|}_{\mathit{\sigma }}\right),\end{array}\text{(76b)}& {\displaystyle }\begin{array}{rl}{G}_{\text{iso}}& ={\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial x}}\right|}_{\mathit{\sigma }}\left[{\mathit{\zeta }}_{\mathit{\sigma }}{\mathit{\kappa }}_{\text{iso}}\left({\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}\right|}_{\mathit{\sigma }}\right.\right.\\ & \left.\left.+(S_{\text{iso}}-S_{\mathit{\sigma }}){\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial \mathit{\sigma }}}\right|}_x\right)\right],\\ & +{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left[{\mathit{\kappa }}_{\text{iso}}(S_{\text{iso}}-S_{\mathit{\sigma }})\left({\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial x}}\right|}_{\mathit{\sigma }}\right.\right.\\ & \left.\left.+(S_{\text{iso}}-S_{\mathit{\sigma }}){\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial \mathit{\sigma }}}\right|}_x\right)\right],\end{array}\text{(76c)}& {\displaystyle }{G}_{\text{dia}}={\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}\right|}_x\left({\mathit{\kappa }}_{\text{dia}}{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial \mathit{\sigma }}}\right|}_x\right)\text{(76d)}& {\displaystyle }{G}_{\text{lat}}=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{v}}{L_y}}{\left[c\right]}_{\mathrm{0}}^{L_y}.\end{array}$$

Here, we define

$$\begin{array}{}\text{(77)}& {\mathit{\psi }}^{†}=\stackrel{\mathrm{‾}}{\mathit{\psi }}+{\mathit{\psi }}^{\star }\end{array}$$

as the total advective or “residual” stream function (Plumb and Ferrari, 2005b), and we have added factors of ζ_σ in Eq. (76a) so that the fluxes can be directly identified with the zonal velocity, $u^{†}=-{\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}\partial {\mathit{\psi }}^{†}/\partial \mathit{\sigma }{|}_x$, and the dia-σ velocity, ${\mathit{\varpi }}^{†}=\partial {\mathit{\psi }}^{†}/\partial x{|}_{\mathit{\sigma }}$. Note also that every derivative with respect to σ is multiplied by ${\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}$, and that their product ${\mathit{\zeta }}_{\mathit{\sigma }}^{-\mathrm{1}}{\partial }_{\mathit{\sigma }}$ is equivalent to a derivative with respect to z. This allows us to simplify the numerical discretization by avoiding explicit references to σ coordinates and computing these derivatives via finite differencing in z coordinates.

Figure 2 Illustration of the numerical grid used to compute solutions to the model equations.

3.2 Spatial discretization of the tracer evolution equation

We solve Eq. (9) using the slope-limited finite-volume scheme of Kurganov and Tadmor (2000) for systems of conservation laws. We divide the domain into a grid of N_x by N_ζ cells, with uniform side lengths Δ_x and Δ_ζ in x/ζ space, as shown in Fig. 2. We store the cell-averaged value of $\stackrel{\mathrm{‾}}{c}$ at the center of the (j,k)^th grid cell, which we denote as ${\stackrel{\mathrm{‾}}{c}}_{j,k}\left(t\right)$. The mean, eddy, and residual stream functions are most naturally defined at the cell corners, as this allows a straightforward calculation of the residual velocities at the cell edges:

$$\begin{array}{}\text{(78a)}& {\displaystyle }& {\displaystyle }{u}_{j+\mathrm{1}/\mathrm{2},k}^{†}=-{\displaystyle \frac{{\mathit{\psi }}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}^{†}-{\mathit{\psi }}_{j+\mathrm{1}/\mathrm{2},k-\mathrm{1}/\mathrm{2}}^{†}}{({\mathrm{\Delta }}_z{)}_{j+\mathrm{1}/\mathrm{2},k}}},\text{(78b)}& {\displaystyle }& {\displaystyle }{\mathit{\varpi }}_{j,k+\mathrm{1}/\mathrm{2}}^{†}={\displaystyle \frac{{\mathit{\psi }}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}^{†}-{\mathit{\psi }}_{j-\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}^{†}}{{\mathrm{\Delta }}_x}}.\end{array}$$

Here, we use Δ_z as a shorthand for the spatially varying vertical grid spacing, defined as a centered difference between adjacent grid points. The vertical grid spacing is defined for all cell centers, corners, and faces:

$$\begin{array}{}\text{(78c)}& {\displaystyle }({\mathrm{\Delta }}_z{)}_{j,k}=z_{j,k+\mathrm{1}/\mathrm{2}}-z_{j,k-\mathrm{1}/\mathrm{2}},\text{(78d)}& {\displaystyle }({\mathrm{\Delta }}_z{)}_{j+\mathrm{1}/\mathrm{2},k}=z_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}-z_{j+\mathrm{1}/\mathrm{2},k-\mathrm{1}/\mathrm{2}},\text{(78e)}& {\displaystyle }({\mathrm{\Delta }}_z{)}_{j,k+\mathrm{1}/\mathrm{2}}=z_{j,k+\mathrm{1}}-z_{j,k},\text{(78f)}& {\displaystyle }({\mathrm{\Delta }}_z{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}=z_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}}-z_{j+\mathrm{1}/\mathrm{2},k},\end{array}$$

where $z_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}$ denotes the physical elevation of each grid point. Note again that ϖ is the velocity normal to the upper and lower faces of the grid cell and so differs slightly from the true vertical velocity w.

To compute the advective tendency, G_adv, the Kurganov and Tadmor (2000) scheme requires a linear interpolation of $\stackrel{\mathrm{‾}}{c}$ over each grid cell. The linear slopes in the x and ζ directions around ${\stackrel{\mathrm{‾}}{c}}_{j,k}\left(t\right)$ are calculated via slope-limited finite differences between ${\stackrel{\mathrm{‾}}{c}}_{j,k}\left(t\right)$ and its adjacent grid points:

$$\begin{array}{}\text{(79)}& {\displaystyle }\begin{array}{rl}({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j,k}& =\mathrm{minmod}\left(\mathit{\theta }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1},k}-{\stackrel{\mathrm{‾}}{c}}_{j,k}}{{\mathrm{\Delta }}_x}},{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1},k}-{\stackrel{\mathrm{‾}}{c}}_{j-\mathrm{1},k}}{\mathrm{2}{\mathrm{\Delta }}_x}},\right.\\ & \left.\mathit{\theta }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k}-{\stackrel{\mathrm{‾}}{c}}_{j-\mathrm{1},k}}{{\mathrm{\Delta }}_x}}\right),\end{array}\text{(80)}& {\displaystyle }\begin{array}{rl}({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j,k}& =\mathrm{minmod}\left(\mathit{\theta }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}}-{\stackrel{\mathrm{‾}}{c}}_{j,k}}{({\mathrm{\Delta }}_z{)}_{j,k+\mathrm{1}/\mathrm{2}}}},{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}}-{\stackrel{\mathrm{‾}}{c}}_{j,k-\mathrm{1}}}{z_{j,k+\mathrm{1}}-z_{j,k-\mathrm{1}}}},\right.\\ & \left.\mathit{\theta }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k}-{\stackrel{\mathrm{‾}}{c}}_{j,k-\mathrm{1}}}{({\mathrm{\Delta }}_z{)}_{j,k-\mathrm{1}/\mathrm{2}}}}\right),\end{array}\end{array}$$

with parameter $\mathrm{1}<\mathit{\sigma }<\mathrm{2}$. The “minmod” function evaluates to zero if its arguments have differing signs and otherwise evaluates to its argument with smallest modulus. The cell center estimates of the derivatives are then used to construct two different estimates of $\stackrel{\mathrm{‾}}{c}$ at each cell face via

$$\begin{array}{}\text{(81a)}& {\displaystyle }{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1}/\mathrm{2},k}^{(-)}& {\displaystyle }={\stackrel{\mathrm{‾}}{c}}_{j,k}+\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Delta }}_x({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j,k},\text{(81b)}& {\displaystyle }{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1}/\mathrm{2},k}^{(+)}& {\displaystyle }={\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1},k}-\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Delta }}_x({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j+\mathrm{1},k},\text{(81c)}& {\displaystyle }{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}/\mathrm{2}}^{(-)}& {\displaystyle }={\stackrel{\mathrm{‾}}{c}}_{j,k}+(z_{j,k+\mathrm{1}/\mathrm{2}}-z_{j,k})({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j,k},\text{(81d)}& {\displaystyle }{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}/\mathrm{2}}^{(+)}& {\displaystyle }={\stackrel{\mathrm{‾}}{c}}_{j,k}-(z_{j,k+\mathrm{1}}-z_{j,k+\mathrm{1}/\mathrm{2}})({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j,k+\mathrm{1}/\mathrm{2}}.\end{array}$$

Finally, advective fluxes are determined at the cell faces using an estimate of the maximum propagation speed in the system, which in our case is simply the residual velocity, and thus the fluxes reduce to an upwind approximation:

$$\begin{array}{}\text{(82a)}& {\displaystyle }\begin{array}{rl}{F}_{j+\mathrm{1}/\mathrm{2},k}^{\left(u\right)}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathrm{\Delta }}_z{)}_{j+\mathrm{1}/\mathrm{2},k}\left[u_{j+\mathrm{1}/\mathrm{2},k}^{†}\left({\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1}/\mathrm{2},k}^{(+)}+{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1}/\mathrm{2},k}^{(-)}\right)\right.\\ & \left.-\left|u_{j+\mathrm{1}/\mathrm{2},k}^{†}\right|\left({\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1}/\mathrm{2},k}^{(+)}-{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1}/\mathrm{2},k}^{(-)}\right)\right],\end{array}\text{(82b)}& {\displaystyle }\begin{array}{rl}{F}_{j,k+\mathrm{1}/\mathrm{2}}^{\left(\mathit{\varpi }\right)}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\mathit{\varpi }}_{j,k+\mathrm{1}/\mathrm{2}}^{†}\left({\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}/\mathrm{2}}^{(+)}+{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}/\mathrm{2}}^{(-)}\right)\right.\\ & \left.-\left|{\mathit{\varpi }}_{j,k+\mathrm{1}/\mathrm{2}}^{†}\right|\left({\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}/\mathrm{2}}^{(+)}-{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}/\mathrm{2}}^{(-)}\right)\right].\end{array}\end{array}$$

For this version of the model, the formulation of the Kurganov and Tadmor (2000) scheme considers only the maximum propagation speed of the momentum, u, and excludes the internal gravity wave speed which is supported with the momentum calculation in Sect. 2.2. As a result, this would alter the overall advective fluxes; however, we omit this in the current version of the model and note that the full formulation can be implemented here, but we choose to leave this calculation to be updated in a future version of the model. The advective tendency in (x,z) space is then computed via straightforward finite differencing of these fluxes:

$$\begin{array}{}\text{(83)}& \begin{array}{rl}(G_{\text{adv}}{)}_{j,k}& =-{\displaystyle \frac{F_{j+\mathrm{1}/\mathrm{2},k}^{\left(u\right)}-F_{j-\mathrm{1}/\mathrm{2},k}^{\left(u\right)}}{{\mathrm{\Delta }}_x({\mathrm{\Delta }}_z{)}_{j,k}}}\\ & -{\displaystyle \frac{F_{j,k+\mathrm{1}/\mathrm{2}}^{\left(\mathit{\varpi }\right)}-F_{j,k-\mathrm{1}/\mathrm{2}}^{\left(\mathit{\varpi }\right)}}{({\mathrm{\Delta }}_z{)}_{j,k}}}-{\displaystyle \frac{F_{j,k}^v}{L_y}}.\end{array}\end{array}$$

The advective discretization (83) requires the residual stream function to be known on all grid cell corners, which allows the numerical fluxes to be computed at the cell edges. The mean stream function $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ is computed from the mean velocity field via Eq. (12):

$$\begin{array}{}\text{(84)}& {\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}=-\sum _{m=\mathrm{0}}^k{\stackrel{\mathrm{‾}}{u}}_{j+\mathrm{1}/\mathrm{2},k}({\mathrm{\Delta }}_z{)}_{j+\mathrm{1}/\mathrm{2},k}.\end{array}$$

The eddy stream function (Eq. 21) depends on the “true” slope (Eq. 72) of the local $\stackrel{\mathrm{‾}}{b}$ contours, which we discretize as

$$\begin{array}{}\text{(85)}& (S_{\text{int}}{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}=-{\displaystyle \frac{({\partial }_x\stackrel{\mathrm{‾}}{b}{|}_z{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}}{({\partial }_z\stackrel{\mathrm{‾}}{b}{|}_x{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}}}.\end{array}$$

The calculation of the derivatives with respect to x and z is described in Sect. 3.5. We then construct ψ^⋆ on cell corners as

$$\begin{array}{}\text{(86)}& {\mathit{\psi }}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}^{\star }=\left({\mathit{\kappa }}_{\text{gm}}{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}\right(S_{\text{int}}{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}.\end{array}$$

The tracer tendency due to isopycnal diffusion, G_iso, is discretized following the formulation of Kurganov and Tadmor (2000) for parabolic operators.

$$\begin{array}{}\text{(87a)}& {\displaystyle }\begin{array}{rl}{H}_{j+\mathrm{1}/\mathrm{2},k}^{\left(x\right)}& =\left({\mathrm{\Delta }}_z{)}_{j+\mathrm{1}/\mathrm{2},k}\right({\mathit{\kappa }}_{\text{iso}}{)}_{j+\mathrm{1}/\mathrm{2},k}\\ & \cdot \left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j+\mathrm{1},k}-{\stackrel{\mathrm{‾}}{c}}_{j,k}}{{\mathrm{\Delta }}_x}}+(S_{\text{iso}}-S_{\mathit{\sigma }}{)}_{j+\mathrm{1}/\mathrm{2},k}\right.\\ & \left.\cdot {\displaystyle \frac{({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j,k}+({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j+\mathrm{1},k}}{\mathrm{2}}}\right],\end{array}\text{(87b)}& {\displaystyle }\begin{array}{rl}{H}_{j,k+\mathrm{1}/\mathrm{2}}^{\left(\mathit{\sigma }\right)}& =\left({\mathit{\kappa }}_{\text{iso}}{)}_{j,k+\mathrm{1}/\mathrm{2}}\right(S_{\text{iso}}-S_{\mathit{\sigma }}{)}_{j,k+\mathrm{1}/\mathrm{2}}\\ & \cdot \left[{\displaystyle \frac{({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j,k}+({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j,k+\mathrm{1}}}{\mathrm{2}}}\right.\\ & \left.+(S_{\text{iso}}-S_{\mathit{\sigma }}{)}_{j,k+\mathrm{1}/\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}}-{\stackrel{\mathrm{‾}}{c}}_{j,k}}{(\mathrm{\Delta }z{)}_{j,k+\mathrm{1}/\mathrm{2}}}}\right],\end{array}\end{array}$$

where $({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j,k}$, $({\partial }_z\stackrel{\mathrm{‾}}{c}{)}_{j+\mathrm{1},k}$, $({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j,k}$, and $({\partial }_x\stackrel{\mathrm{‾}}{c}{)}_{j,k+\mathrm{1}}$ are computed via Eq. (81a)–(81d). In the interior, the isopycnal diffusion slope S_iso=S_int and is calculated on cell corners via Eq. (85) and interpolated to cell faces via

$$\begin{array}{}\text{(88a)}& {\displaystyle }\begin{array}{rl}(S_{\text{iso}}{)}_{j+\mathrm{1}/\mathrm{2},k}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left((S_{\text{iso}}{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}\right.\\ & \left.+(S_{\text{iso}}{)}_{j+\mathrm{1}/\mathrm{2},k-\mathrm{1}/\mathrm{2}}\right),\end{array}\text{(88b)}& {\displaystyle }\begin{array}{rl}(S_{\text{iso}}{)}_{j,k+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left((S_{\text{iso}}{)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}\right.\\ & \left.+(S_{\text{iso}}{)}_{j-\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}\right).\end{array}\end{array}$$

The diffusive tendency is then computed via straightforward finite differencing of the H fluxes:

$$\begin{array}{}\text{(89)}& \begin{array}{rl}(G_{\text{iso}}{)}_{j,k}& =-{\displaystyle \frac{H_{j+\mathrm{1}/\mathrm{2},k}^{\left(x\right)}-H_{j-\mathrm{1}/\mathrm{2},k}^{\left(x\right)}}{{\mathrm{\Delta }}_x({\mathrm{\Delta }}_z{)}_{j,k}}}\\ & -{\displaystyle \frac{H_{j,k+\mathrm{1}/\mathrm{2}}^{\left(\mathit{\sigma }\right)}-H_{j,k-\mathrm{1}/\mathrm{2}}^{\left(\mathit{\sigma }\right)}}{({\mathrm{\Delta }}_z{)}_{j,k}}}.\end{array}\end{array}$$

The tracer tendency due to diapycnal mixing, G_dia, is computed implicitly. During each time step, all other physical and biogeochemical tendencies are computed and used to advance ${\stackrel{\mathrm{‾}}{c}}_{j,k}$ forward one time step Δ_t, i.e.,

$$\begin{array}{}\text{(90)}& {\stackrel{\mathrm{‾}}{c}}_{j,k}^{\star }={\stackrel{\mathrm{‾}}{c}}_{j,k}^n+\mathcal{F}\left[{\stackrel{\mathrm{‾}}{c}}^n\right].\end{array}$$

Here, n denotes the time step number, and ${\stackrel{\mathrm{‾}}{c}}^{\star }$ denotes an estimate of $\stackrel{\mathrm{‾}}{c}$ at t+Δ_t (see Sect. 3.3 for details of the time-stepping schemes). The updated tracer concentration is then further modified via the addition of a “correction” due to diapycnal diffusion. At each longitude, or for each j, we solve

$$\begin{array}{}\text{(91)}& \begin{array}{rl}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k}^{n+\mathrm{1}}-{\stackrel{\mathrm{‾}}{c}}_{j,k}^{\star }}{{\mathrm{\Delta }}_{\mathrm{t}}}}& ={\displaystyle \frac{\mathrm{1}}{z_{j,k+\mathrm{1}/\mathrm{2}}-z_{j,k-\mathrm{1}/\mathrm{2}}}}\\ & \cdot \left[({\mathit{\kappa }}_{\text{dia}}{)}_{j,k+\mathrm{1}/\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k+\mathrm{1}}^{n+\mathrm{1}}-{\stackrel{\mathrm{‾}}{c}}_{j,k}^{n+\mathrm{1}}}{z_{j,k+\mathrm{1}}-z_{j,k}}}\right.\\ & \left.-({\mathit{\kappa }}_{\text{dia}}{)}_{j,k-\mathrm{1}/\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{c}}_{j,k}^{n+\mathrm{1}}-{\stackrel{\mathrm{‾}}{c}}_{j,k-\mathrm{1}}^{n+\mathrm{1}}}{z_{j,k}-z_{j,k-\mathrm{1}}}}\right].\end{array}\end{array}$$

Equation (91) defines a tridiagonal matrix system of algebraic equations for the unknowns $\mathit{\left\{}{\stackrel{\mathrm{‾}}{c}}_{j,k}^{n+\mathrm{1}}\right|k=\mathrm{1}\mathrm{\dots }{N}_z\mathit{\}}$, which is inverted using the Thomas algorithm.

Finally, the meridional advection is discretized via a straightforward upwind advection scheme:

$$\begin{array}{}\text{(92)}& \begin{array}{rl}& (G_{\text{lat}}{)}_{j,k}={\displaystyle \frac{v_{j,k}^{\left(c\right)}}{L_y}}\left(c_{j,k}-c_{j,k}^u\right),\\ & v_{j,k}^{\left(c\right)}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(v_{j+\mathrm{1}/\mathrm{2},k}^{†}+v_{j-\mathrm{1}/\mathrm{2},k}^{†}\right),\end{array}\end{array}$$

where v^(c) denotes the meridional velocity on tracer points and c^u denotes the upstream tracer concentration, defined as

$$\begin{array}{}\text{(93)}& c_{j.k}^{\mathrm{u}}=\left\{\begin{array}{ll}{c}_{j,k}^{\mathrm{N}},& v<\mathrm{0},\\ c_{j,k}^{\mathrm{S}},& v>\mathrm{0}.\end{array}\right.\end{array}$$

Here, c^N and c^S are the tracer concentrations at the northern and southern ends of the domain, respectively. In all of the steps listed above, conditions of zero residual stream function and zero normal tracer flux are applied at the domain boundaries. These conditions are imposed by simply setting ψ^† to zero on all boundary points and by setting the numerical fluxes (F, H, etc.) to zero on the boundary cell faces.

3.3 Temporal discretization

MAMEBUS evolves the model equations forward in time using Adams–Bashforth (AB) methods Durran (e.g., 1991) modified to allow for adaptive-time-step sizes. In this section, we outline the derivation of these methods and formally show the derivation in Appendix B. We implement the adaptive-time-step AB methods because this family of methods is numerically stable with our scheme for the momentum equations (see Sect. 2.2). We then describe the constraints on the time step imposed by the Courant–Friedrichs–Lewy (CFL) condition.

3.3.1 Adaptive-time-step Adams–Bashforth methods

Our time-integration scheme uses a family of time-step-variable Adams–Bashforth integrative methods. This specific formulation of the AB methods allows for the model time step to be adjusted dynamically following the CFL conditions described in Sect. 3.3.2. Consider a tracer quantity c that evolves according to

$$\begin{array}{}\text{(94)}& {\displaystyle \frac{\partial c}{\partial t}}=f(t,c(t\left)\right).\end{array}$$

Here, the function f conceptually represents the entire model state, including the physical, biogeochemical, and non-conservative tendencies. We make a note here that the diffusive component of the time-integration step is calculated implicitly and not included in the ABIII integration step (see Eq. 91). We implement the third-order Adams–Bashforth, or ABIII method, in this version of the model as the default option for time integration.

$$\begin{array}{}\text{(95)}& \begin{array}{rl}\text{ABIII:}c\left(t_{n+\mathrm{1}}\right)& =c\left(t_n\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\left(f^{n-\mathrm{2}}{\displaystyle \frac{\mathrm{\Delta }{t}_{n+\mathrm{1}}^{\mathrm{2}}(\mathrm{2}\mathrm{\Delta }{t}_{n+\mathrm{1}}+\mathrm{3}\mathrm{\Delta }{t}_n)}{\mathrm{\Delta }{t}_{n-\mathrm{1}}(\mathrm{\Delta }{t}_n+\mathrm{\Delta }{t}_{n-\mathrm{1}})}}\right.\\ & -f^{n-\mathrm{1}}{\displaystyle \frac{\mathrm{\Delta }{t}_{n+\mathrm{1}}^{\mathrm{2}}(\mathrm{2}\mathrm{\Delta }{t}_{n+\mathrm{1}}+\mathrm{3}\mathrm{\Delta }{t}_n+\mathrm{3}\mathrm{\Delta }{t}_{n-\mathrm{1}})}{\mathrm{\Delta }{t}_{n-\mathrm{1}}\mathrm{\Delta }{t}_n}}\\ & +\left.f^n{\displaystyle \frac{\begin{array}{c}\mathrm{\Delta }{t}_{n+\mathrm{1}}(\mathrm{2}\mathrm{\Delta }{t}_{n+\mathrm{1}}^{\mathrm{2}}+\mathrm{6}\mathrm{\Delta }{t}_{n+\mathrm{1}}\mathrm{\Delta }{t}_n+\mathrm{3}\mathrm{\Delta }{t}_{n+\mathrm{1}}\\ \mathrm{\Delta }{t}_{n-\mathrm{1}}+\mathrm{6}\mathrm{\Delta }{t}_n^{\mathrm{2}}+\mathrm{6}\mathrm{\Delta }{t}_n\mathrm{\Delta }{t}_{n-\mathrm{1}})\end{array}}{(\mathrm{\Delta }{t}_n+\mathrm{\Delta }{t}_{n-\mathrm{1}})\mathrm{\Delta }{t}_n}}\right)\end{array}\end{array}$$

The first two time steps require the lower-order methods. We implement a forward Euler for the first time step and a second-order AB scheme (defined below) for the second time step:

$$\begin{array}{}\text{(96)}& \begin{array}{rl}\text{ABII:}c\left(t_{n+\mathrm{1}}\right)& =c\left(t_n\right)+{\displaystyle \frac{\mathrm{\Delta }{t}_{n+\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }{t}_n}}\left(\mathrm{2}f(t_n,c(t_n\left)\right)\mathrm{\Delta }{t}_n\right.\\ & +f(t_n,c(t_n\left)\right)\mathrm{\Delta }{t}_{n+\mathrm{1}}\\ & \left.-f(t_{n-\mathrm{1}},c(t_{n-\mathrm{1}}\left)\right)\mathrm{\Delta }{t}_{n+\mathrm{1}}\right).\end{array}\end{array}$$

Here, the notation Δt_n indicates the nth time step. Derivations for the adaptive-time-stepping ABII and ABIII methods are given in Appendix B.

3.3.2 CFL conditions

MAMEBUS selects each model time step adaptively to ensure that time stepping is numerically stable. The time step is chosen to ensure that the CFL conditions for each of MAMEBUS's various advective and diffusive operators, described in preceding subsections, are satisfied.

The time step for advection of tracers is limited by the timescale associated with advective propagation across the width of a grid box (Δ_x or Δ_z). These constraints can approximately be expressed as

$$\begin{array}{}\text{(97a)}& {\displaystyle }\mathrm{\Delta }t<& {\displaystyle \frac{{\mathrm{\Delta }}_x}{\left|\stackrel{\mathrm{‾}}{u}\right|+u_{\text{igw}}}},\text{(97b)}& {\displaystyle }\mathrm{\Delta }t<& {\displaystyle \frac{{\mathrm{\Delta }}_z}{\left|\stackrel{\mathrm{‾}}{w}\right|}}\end{array}$$

(Durran, 2010). Here, u_igw is the maximum horizontal propagation speed of internal gravity waves (Chelton et al., 1998):

$$\begin{array}{}\text{(98)}& c_{\text{igw}}={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\int N\mathrm{d}z.\end{array}$$

Particulate sinking in the NPZD model is also calculated explicitly and constrains the time step via a similar CFL criterion:

$$\begin{array}{}\text{(99)}& \mathrm{\Delta }t<{\displaystyle \frac{\left({\mathrm{\Delta }}_z\right)}{\left|w_{\text{sink}}\right|}},\end{array}$$

where w_sink is the sinking speed of the particles.

We apply additional constraints on the time step to ensure that diffusive operators are stable. The standard numerical stability criterion for a Laplacian diffusion operator is (Griffies, 2018)

$$\begin{array}{}\text{(100)}& \mathrm{\Delta }t<{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{s}}^{\mathrm{2}}}{\mathit{\kappa }}},\end{array}$$

where κ is a diffusion coefficient and Δ_s is the spatial grid spacing. In the horizontal (Δ_s=Δ_x, the diffusion coefficients that determine the diffusive time step are the eddy diffusion and isopycnal diffusion coefficients, when κ=κ_gm and κ=κ_iso, respectively. ¹ In the vertical (Δ_s=Δ_z), the diffusive time step is determined by the diapycnal diffusivity, κ=κ_dia, and by the vertical component of the eddy and isopycnal diffusion operators, $\mathit{\kappa }={\mathit{\kappa }}_{\text{gm}}{S}_{\text{int}}^{\mathrm{2}}$ and $\mathit{\kappa }={\mathit{\kappa }}_{\text{iso}}{S}_{\text{int}}^{\mathrm{2}}$, respectively (see, e.g., Ferrari et al., 2008).

3.4 Discrete momentum equations and barotropic pressure correction

In this section, we describe the discretization of the momentum equations presented in Sect. 2.2, specifically in Eq. (11a), (11b), and (11c). To facilitate our discretization, we split the pressure, ϕ, into barotropic and baroclinic components:

$$\begin{array}{}\text{(101)}& \mathit{\varphi }=\underset{\text{barotropic}}{\underbrace{\mathrm{\Pi }}}+\underset{\text{baroclinic}}{\underbrace{{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{z}{\overset{\mathrm{0}}{\int }}\mathit{\rho }\mathrm{d}z}}.\end{array}$$

The barotropic and baroclinic components correspond to the pressure at the surface and the hydrostatic pressure variation with depth, respectively.

The numerical time integration is calculated in a series of steps which include an explicit calculation of the non-diffusive time step, an implicit calculation of the vertical diffusion, and a barotropic corrector step in order to ensure that the flow is non-divergent. The calculation of the explicit time step is outlined in Sect. 3.3 and Appendix B. In order to be numerically consistent with the calculation of the stream function, the mean horizontal velocities $\stackrel{\mathrm{‾}}{u}$ and $\stackrel{\mathrm{‾}}{v}$ are stored on the western face of each grid cell. In Fig. 2, these are labeled as the u points. The time step calculation is shown below, noting that the explicit components of the time step, ℰ{⋅} are calculated following the ABIII methods outlined in Sect. 3.3, and the implicit diffusion ℐ{⋅} is calculated following Eq. (91).

Given the mean momentum at time step n, ${\stackrel{\mathrm{‾}}{\mathit{u}}}^n$, we first perform the explicit component of the time step to construct an estimate of ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{n+\mathrm{1}}$, denoted as ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star }$:

$$\begin{array}{}\text{(102)}& {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star }={\stackrel{\mathrm{‾}}{\mathit{u}}}^n+\mathcal{E}\left\{-f\widehat{z}\times {\mathit{u}}^n-{\displaystyle \frac{\partial \mathrm{\Pi }}{\partial y}}\widehat{\mathit{y}}-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }\underset{z}{\int }{\mathit{\rho }}^n\mathrm{d}z\right\}.\end{array}$$

Note that the zonal barotropic pressure gradient, ${\partial }_x\mathrm{\Pi }\widehat{\mathit{y}}$, is excluded from this equation; this will be revisited in the final component of the time step. The discretization of the horizontal pressure-gradient terms in Eq. (102) is described in Sect. 3.5.

We next compute the tendency due to vertical viscosity following Eq. (91), which we denote via the operator ℐ. We thereby construct a second estimate of the velocity at time step n+1, denoted as ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star \star }$:

$$\begin{array}{}\text{(103)}& {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star \star }=\mathcal{I}\mathit{\left\{}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star }\mathit{\right\}}.\end{array}$$

Finally, we apply a tendency due to the zonal barotropic pressure gradient, ensuring that mass is conserved in each vertical fluid column (Dauhajre and McWilliams, 2018):

$$\begin{array}{}\text{(104)}& \underset{z}{\int }{\stackrel{\mathrm{‾}}{u}}^{n+\mathrm{1}}\mathrm{d}z=\underset{z}{\int }{\stackrel{\mathrm{‾}}{u}}^n\mathrm{d}z=\mathrm{0},\end{array}$$

as required by Eq. (11c). We formulate the barotropic pressure correction as

$$\begin{array}{}\text{(105)}& {\stackrel{\mathrm{‾}}{u}}^{n+\mathrm{1}}={\stackrel{\mathrm{‾}}{u}}^{\star \star }-\mathrm{\Delta }t{\displaystyle \frac{\partial \mathrm{\Pi }}{\partial x}}.\end{array}$$

Substituting Eq. (105) into (104), we obtain

$$\begin{array}{}\text{(106)}& \mathrm{\Delta }t{\displaystyle \frac{\partial \mathrm{\Pi }}{\partial x}}\widehat{\mathit{x}}={\displaystyle \frac{\mathrm{1}}{\left|{\mathit{\eta }}_{\mathrm{b}}\right(x\left)\right|}}\underset{z}{\int }{\stackrel{\mathrm{‾}}{u}}^{\star \star }\mathrm{d}z.\end{array}$$

This implies that the tendency in the mean zonal velocity due to the zonal barotropic pressure gradient must serve to bring the depth-integrated zonal velocity to zero, i.e.,

$$\begin{array}{}\text{(107)}& {\stackrel{\mathrm{‾}}{u}}^{n+\mathrm{1}}={\stackrel{\mathrm{‾}}{u}}^{\star \star }-{\displaystyle \frac{\mathrm{1}}{\left|{\mathit{\eta }}_{\mathrm{b}}\right(x\left)\right|}}\underset{z}{\int }{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star \star }\mathrm{d}z.\end{array}$$

The calculation of the vertical integral of ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{\star \star }$ is computed in the model using a Kahan sum (Kahan, 1965).

3.5 Horizontal pressure- and buoyancy-gradient calculations

Pressure-gradient calculations in σ coordinates have been long known to produce discretization errors from the misalignment of geopotential and σ coordinate surfaces and rely on large cancelations in the vertical gradient near steep slopes (Arakawa and Suarez, 1983; Haney, 1991; Mellor et al., 1994, 1998). We follow Shchepetkin and McWilliams (2003) to calculate the horizontal-pressure-gradient force and reduce the errors in horizontal gradient calculations, which otherwise produce large spurious along-slope currents in MAMEBUS (not shown). This algorithm has been extensively tested via its implementation in ROMS (Shchepetkin and McWilliams, 2003, 2005), so we omit our own tests of the pressure-gradient calculation scheme in this study. For numerical consistency, we also calculate horizontal buoyancy gradients, required to evaluate the isopycnal slope (see Sect. 3.2), using the same algorithm.

Figure 3 Stencil for the isopycnal slope and pressure-gradient scheme given by Shchepetkin and McWilliams (2003). The points indicate the buoyancy (density) points. The solid lines are the reconstructed coordinate lines used in the horizontal calculation, and the shaded area shows the area integral of the horizontal buoyancy gradient.

3.5.1 Zonal pressure gradients

In this section, we outline the implementation of the zonal baroclinic pressure-gradient calculation used in MAMEBUS following Shchepetkin and McWilliams (2003). The ultimate goal of the algorithm is to calculate the following baroclinic pressure gradient at a cell center (see Eq. 102):

$$\begin{array}{}\text{(108)}& -{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial x}}{|}_z^{\text{baroclinic}}=-{\displaystyle \frac{\mathrm{1}}{A}}{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\left(\underset{A}{\iint }{\displaystyle \frac{\partial \mathit{\rho }}{\partial x}}{|}_z\mathrm{d}x\mathrm{d}z\right),\end{array}$$

where ρ, like the linear case, is the density anomaly, ρ₀ is a reference density, and A is the area between four adjacent buoyancy grid points (shaded area in Fig. 3). Shchepetkin and McWilliams (2003) calculate the second term by implementing a Lagrangian polynomial reconstruction of the z and ρ fields. By Green's theorem, we write the integrated horizontal density gradient as

$$\begin{array}{}\text{(109)}& \underset{A}{\iint }{\displaystyle \frac{\partial \mathit{\rho }}{\partial x}}{|}_z\mathrm{d}x\mathrm{d}z=\oint \mathit{\rho }\mathbf{d}\mathit{z}.\end{array}$$

The two algorithms differ on the treatment of the vertical integration of pressure. The density Jacobian algorithm first interpolates the density field onto the σ grid and calculates the values of ρ and z along the solid lines in Fig. 3, then integrates to obtain the pressure. The σ coordinate primitive form algorithm first integrates the pressure and calculates the gradients using a vertical correction term. We opt for the second algorithm which results in an equation for the horizontal gradient which closely resembles Eq. (68a). Furthermore, this algorithm tends to be more stable in our model. The algorithm is calculated as follows:

First, we calculate all elementary differences in ρ and z:

$$\begin{array}{}\text{(110a)}& {\displaystyle }\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{j+\mathrm{1}/\mathrm{2},k}& {\displaystyle }={\mathit{\rho }}_{j+\mathrm{1},k}-{\mathit{\rho }}_{j,k},\text{(110b)}& {\displaystyle }\left({\mathrm{\Delta }}_z\right){\mathit{\rho }}_{j,k+\mathrm{1}/\mathrm{2}}& {\displaystyle }={\mathit{\rho }}_{j,k+\mathrm{1}}-{\mathit{\rho }}_{j,k},\text{(110c)}& {\displaystyle }\left({\mathrm{\Delta }}_x\right)z_{j+\mathrm{1}/\mathrm{2},k}& {\displaystyle }=z_{j+\mathrm{1},k}-z_{j,k},\text{(110d)}& {\displaystyle }\left({\mathrm{\Delta }}_z\right)z_{j,k+\mathrm{1}/\mathrm{2}}& {\displaystyle }=z_{j,k+\mathrm{1}}-z_{j,k},\end{array}$$

where z_j,k is the depth value at the cell centers, where the density tracer is located. Note that the edges of Fig. 3 correspond to the cell centers in MAMEBUS; this requires some extrapolation at the boundaries so that the elementary differences are fully calculated throughout the domain. For all variables, we assume that the elementary differences at the boundary are zero.

We then calculate the hyperbolic differences of all variables. This step calculates an estimate of the derivatives following a cubic spline formalism outlined in more detail in Shchepetkin and McWilliams (2003). The derivatives are then given by

$$\begin{array}{}\text{(111a)}& {\displaystyle }{h}_x{\mathit{\rho }}_{j,k}& {\displaystyle }={\displaystyle \frac{\mathrm{2}\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{j+\mathrm{1}/\mathrm{2},k}\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{j-\mathrm{1}/\mathrm{2},k}}{\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{j+\mathrm{1}/\mathrm{2},k}+\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{j-\mathrm{1}/\mathrm{2},k}}},\text{(111b)}& {\displaystyle }{h}_x{z}_{j,k}& {\displaystyle }={\displaystyle \frac{\mathrm{2}\left({\mathrm{\Delta }}_x\right)z_{j+\mathrm{1}/\mathrm{2},k}\left({\mathrm{\Delta }}_x\right)z_{j-\mathrm{1}/\mathrm{2},k}}{\left({\mathrm{\Delta }}_x\right)z_{j+\mathrm{1}/\mathrm{2},k}+\left({\mathrm{\Delta }}_x\right)z_{j-\mathrm{1}/\mathrm{2},k}}}.\end{array}$$

The vertical hyperbolic differences for ρ and z are calculated similarly using Eq. (110b) and (110d). Again, the hyperbolic differences on the boundaries are not defined using Eq. (111a) and (111b), so we extrapolate the hyperbolic averages of density. For example, at the western edge of the domain, we define

$$\begin{array}{}\text{(112)}& \begin{array}{rl}{h}_x{\mathit{\rho }}_{N,k}& ={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left(\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{N-\mathrm{1}/\mathrm{2},k}-\left({\mathrm{\Delta }}_x\right){\mathit{\rho }}_{N-\mathrm{3}/\mathrm{2},k}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(h_x{\mathit{\rho }}_{N-\mathrm{1},k}\right).\end{array}\end{array}$$

Analogous extrapolation schemes are applied at all domain boundaries.

We the calculate the pressure field using the hydrostatic relationship. This is done via a vertical integration of the density reconstructed along the vertical lines in Fig. 3. The pressure field is calculated from the surface down. The pressure is calculated in the surface grid cells as

$$\begin{array}{}\text{(113)}& {\mathit{\varphi }}_{j,N}={\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\mathit{\rho }}_{j,N}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\zeta }}_j-z_{j,N}\right){\displaystyle \frac{{\mathit{\rho }}_{j,N}-{\mathit{\rho }}_{j,N-\mathrm{1}}}{z_{j,N}-z_{j,N-\mathrm{1}}}}\right)\left({\mathit{\zeta }}_j-z_{j,N}\right),\end{array}$$

where ζ_j=0 from the rigid lid assumption in MAMEBUS. Then the pressure is calculated at successively deeper grid levels as

$$\begin{array}{}\text{(114)}& \begin{array}{rl}{\mathit{\varphi }}_{j,k}={\mathit{\varphi }}_{j,k+\mathrm{1}}& +{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{j,k+\mathrm{1}}+{\mathit{\rho }}_{j,k}}{\mathrm{2}}}\left(z_{j,k+\mathrm{1}}-z_{j,k}\right)\right.\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{10}}}\left\{(h_z{\mathit{\rho }}_{j,k+\mathrm{1}}-h_z{\mathit{\rho }}_{j,k})\left[z_{j,k+\mathrm{1}}\right.\right.\\ & \left.-z_{j,k}-{\displaystyle \frac{h_z{z}_{j,k+\mathrm{1}}+h_z{z}_{j,k}}{\mathrm{12}}}\right]\\ & -(h_z{z}_{j,k+\mathrm{1}}-h_z{z}_{j,k})\left[{\mathit{\rho }}_{j,k+\mathrm{1}}\right.\\ & \left.\left.\left.-{\mathit{\rho }}_{j,k}-{\displaystyle \frac{h_z{\mathit{\rho }}_{j,k+\mathrm{1}}+h_z{\mathit{\rho }}_{j,k}}{\mathrm{12}}}\right]\right\}\right).\end{array}\end{array}$$

We then correct for the iso-σ pressure gradient introduced by the slope of the σ coordinate grid, analogous to the continuous expression in Eq. (68a). This step calculates the product of ρ and the local slope of the σ coordinate and corrects for the interpolation errors from the coordinate transformation. Following the notation used in Shchepetkin and McWilliams (2003),

$$\begin{array}{}\text{(115)}& \begin{array}{rl}F{C}_{j+\mathrm{1}/\mathrm{2},k}& ={\displaystyle \frac{{\mathit{\rho }}_{j+\mathrm{1},k}+{\mathit{\rho }}_{j,k}}{\mathrm{2}}}(z_{j+\mathrm{1},k}-z_{j,k})\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{10}}}\left\{(h_x{\mathit{\rho }}_{j+\mathrm{1},k}-h_x{\mathit{\rho }}_{j,k})\right.\\ & \cdot \left[z_{j+\mathrm{1},k}-z_{j,k}-{\displaystyle \frac{h_x{z}_{j+\mathrm{1},}+h_x{z}_{j,k}}{\mathrm{12}}}\right]\\ & -(h_x{z}_{j+\mathrm{1},k}-h_x{z}_{j,k})\\ & \left.\left[{\mathit{\rho }}_{j+\mathrm{1},k}-{\mathit{\rho }}_{j,k}-{\displaystyle \frac{h_x{\mathit{\rho }}_{j+\mathrm{1},k}+h_x{\mathit{\rho }}_{j,k}}{\mathrm{12}}}\right]\right\}.\end{array}\end{array}$$

Finally, we use Eqs. (114) and (115) to calculate the pressure gradients:

$$\begin{array}{}\text{(116)}& \begin{array}{rl}& {\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}{|}_z\right)}_{j+\mathrm{1}/\mathrm{2},k}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}\left({\mathit{\varphi }}_{j+\mathrm{1},k}-{\mathit{\varphi }}_{j,k}+{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\cdot F{C}_{j+\mathrm{1}/\mathrm{2},k}\right).\end{array}\end{array}$$

Buoyancy gradients

The buoyancy gradient is calculated similarly to the pressure gradient. However, because we do not vertically integrate the buoyancy term, we opt to use the density Jacobian algorithm described in Shchepetkin and McWilliams (2003). The pressure-gradient algorithm described above integrates the pressure and then corrects for the pressure gradient in σ coordinates. The density gradient algorithm described below calculates the line integral about the area enclosed by the ϕ points where the buoyancy gradient is located (see Fig. 3). Therefore, we use the following form to calculate the buoyancy gradient:

$$\begin{array}{}\text{(117)}& \begin{array}{rl}\underset{A}{\iint }{\displaystyle \frac{\partial b}{\partial x}}\mathrm{d}x\mathrm{d}z& =\oint b\mathrm{d}\mathit{z}\\ & =F{X}_{j+\mathrm{1},k+\mathrm{1}/\mathrm{2}}+F{C}_{j+\mathrm{1}/\mathrm{2},k}-F{X}_{j,k+\mathrm{1}/\mathrm{2}}\\ & -F{C}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}},\end{array}\end{array}$$

where $F{X}_{j,k+\mathrm{1}/\mathrm{2}}$ is the value of the integral (Eq. 117) along the vertical sides, and $F{C}_{j+\mathrm{1}/\mathrm{2},k}$ is the value of the integral along the horizontal sides. This calculation follows a similar procedure to the pressure gradient.

First, we calculate the elementary differences, and the hyperbolic averages in b and z, given by Eqs. (110a)–(111b). Then we calculate the value of the integral along the upper and lower sides of the domain following

$$\begin{array}{}\text{(118)}& \begin{array}{rl}F{C}_{j+\mathrm{1}/\mathrm{2},k}& ={\displaystyle \frac{b_{j+\mathrm{1},k}+b_{j,k}}{\mathrm{2}}}(z_{j+\mathrm{1},k}-z_{j,k})\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{10}}}\left\{(h_x{b}_{j+\mathrm{1},k}-h_x{b}_{j,k})\right.\\ & \left[z_{j+\mathrm{1},k}-z_{j,k}-{\displaystyle \frac{h_x{z}_{j+\mathrm{1},}+h_x{z}_{j,k}}{\mathrm{12}}}\right]\\ & -(h_x{z}_{j+\mathrm{1},k}-h_x{z}_{j,k})\\ & \cdot \left.\left[b_{j+\mathrm{1},k}-b_{j,k}-{\displaystyle \frac{h_x{b}_{j+\mathrm{1},k}+h_x{b}_{j,k}}{\mathrm{12}}}\right]\right\}.\end{array}\end{array}$$

Note that this formulation is the same as Eq. (115) but with buoyancy instead of pressure. Then we calculate the value of the line integral along the vertical components of the cell:

$$\begin{array}{}\text{(119)}& \begin{array}{rl}F{X}_{j,k+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{b_{j,k+\mathrm{1}}+b_{j,k}}{\mathrm{2}}}\left(z_{j,k+\mathrm{1}}-z_{j,k}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{10}}}\left\{(h_z{b}_{j,k+\mathrm{1}}-h_z{b}_{j,k})\right.\\ & \cdot \left[z_{j,k+\mathrm{1}}-z_{j,k}-{\displaystyle \frac{h_z{z}_{j,k+\mathrm{1}}+h_z{z}_{j,k}}{\mathrm{12}}}\right]\\ & -(h_z{z}_{j,k+\mathrm{1}}-h_z{z}_{j,k})\\ & \cdot \left.\left[b_{j,k+\mathrm{1}}-b_{j,k}-{\displaystyle \frac{h_z{b}_{j,k+\mathrm{1}}+h_z{b}_{j,k}}{\mathrm{12}}}\right]\right\}.\end{array}\end{array}$$

Shchepetkin and McWilliams (2003) write Eq. (68a) as

$$\begin{array}{}\text{(120)}& {\displaystyle \frac{\partial b}{\partial x}}{|}_z=\mathcal{J}(b,z)={\displaystyle \frac{\partial b}{\partial x}}{|}_{\mathit{\sigma }}{\displaystyle \frac{\partial z}{\partial \mathit{\sigma }}}-{\displaystyle \frac{\partial b}{\partial \mathit{\sigma }}}{\displaystyle \frac{\partial z}{\partial x}}{|}_{\mathit{\sigma }}.\end{array}$$

This allows us to numerically integrate the buoyancy gradient in the cell as

$$\begin{array}{}\text{(121)}& \begin{array}{rl}\oint b\mathrm{d}\mathit{z}& ={\left(A{\displaystyle \frac{\partial b}{\partial x}}\right)}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}\\ & =F{X}_{j+\mathrm{1},k+\mathrm{1}/\mathrm{2}}+F{C}_{j+\mathrm{1}/\mathrm{2},k}\\ & -F{X}_{j,k+\mathrm{1}/\mathrm{2}}-F{C}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}},\end{array}\end{array}$$

where A, again, is the area of the cell. At the surface, the boundary condition is given that $F{C}_{j+\mathrm{1}/\mathrm{2},N+\mathrm{1}}\equiv \mathrm{0}.$

Finally, in order to calculate the horizontal buoyancy gradient, we divide by the area. Since the area of each cell is defined by the cell-centered ϕ points, we implement Gauss' area formula:

$$\begin{array}{}\text{(122)}& \begin{array}{rl}{A}_{j+\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|x_{j,k+\mathrm{1}}{z}_{j,k}\right.+x_{j,k}{z}_{j+\mathrm{1},k}\\ & +x_{j+\mathrm{1},k}{z}_{j+\mathrm{1},k+\mathrm{1}}+x_{j+\mathrm{1},k+\mathrm{1}}{z}_{j,k+\mathrm{1}}\\ & -x_{j,k}{z}_{j,k+\mathrm{1}}-x_{j+\mathrm{1},k}{z}_{j,k}\\ & \left.-x_{j+\mathrm{1},k+\mathrm{1}}{z}_{j+\mathrm{1},k}-x_{j,k+\mathrm{1}}{z}_{j+\mathrm{1},k+\mathrm{1}}.\right|\end{array}\end{array}$$

3.5.2 Meridional pressure gradients

The along-shore pressure gradient in Eq. (11b), denoted by ${\left[\mathit{\varphi }\right]}_{\mathrm{0}}^{L_y}/L_y$, is determined by along-shore gradients in the surface pressure and buoyancy/density that are imposed as model input parameters. We integrate the profiles of pressure following the hydrostatic relationship. We define ρ^N and ρ^S as the densities at the northern and southern ends of the domain, and ${\mathrm{\Pi }}^y={\mathrm{\Pi }}^{N,S}$ as the surface pressures at northern and southern ends of the domain. Then the pressure is given by

$$\begin{array}{}\text{(123)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{L_y}}{\left[\mathit{\varphi }\right]}_{\mathrm{0}}^{L_y}& ={\displaystyle \frac{\partial \mathrm{\Pi }}{\partial y}}-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{z}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{\partial \mathit{\rho }}{\partial y}}\mathrm{d}z,\\ & ={\displaystyle \frac{{\mathrm{\Pi }}^N-{\mathrm{\Pi }}^S}{{\mathit{\rho }}_{\mathrm{0}}{L}_y}}+{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\underset{z}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{{\mathit{\rho }}^N-{\mathit{\rho }}^S}{L_y}}\mathrm{d}z.\end{array}\end{array}$$

Here, the along-shore variations in sea surface pressure and density are both model inputs. We discretize the meridional pressure gradient as

$$\begin{array}{}\text{(124)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{L_y}}{\left[\mathit{\varphi }\right]}_{\mathrm{0}}^{L_y}{|}_{j,k}& ={\displaystyle \frac{\mathrm{1}}{L_y}}\left[\left({\mathrm{\Pi }}_{j,N}^N-{\mathrm{\Pi }}_{j,N}^S\right)\right.+{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\sum _{k=k^{\prime }}^{N-\mathrm{1}}\\ & \cdot \left({\displaystyle \frac{{\mathit{\rho }}_{j,k^{\prime }+\mathrm{1}}^N+{\mathit{\rho }}_{j,k^{\prime }}^N}{\mathrm{2}}}-{\displaystyle \frac{{\mathit{\rho }}_{j,k^{\prime }+\mathrm{1}}^S+{\mathit{\rho }}_{j,k^{\prime }}^S}{\mathrm{2}}}\right)\\ & \left.\cdot \left(z_{j,k^{\prime }+\mathrm{1}}-z_{j,k^{\prime }}\right)\right].\end{array}\end{array}$$

Though MAMEBUS allows meridional pressure gradients to be imposed, we have excluded them from our reference solutions in the interest of simplicity. However, previous studies have highlighted the importance of meridional pressure gradients in supporting interior cross-slope transport and in driving poleward undercurrents Connolly et al. (2014). We plan to address the effects of meridional pressure gradients on EBUS ecosystem dynamics in future scientific studies using MAMEBUS.

4 Implementation details

In this section, we outline the details for implementation in MAMEBUS. The model code is written in the C programming language. The model expects various user inputs that include initial conditions, along with user-defined model calculation details in Table 4 that include, but are not limited to, the momentum calculation scheme and the time-stepping scheme. The MAMEBUS distribution includes sample MATLAB codes that package these user inputs.

The software needed to run this model includes

1. MATLAB (2016) or later and

2. a C compiler (e.g., GCC).

Our provided setup also includes example scripts for running the model on a cluster; however, this model can be easily run locally on a laptop or desktop on any operating system so long as the necessary software is installed. Table 6 shows run times for the model on both the cluster and a 2015 Mac laptop.

MAMEBUS has three active physical variables: the zonal and meridional momenta, and the temperature (buoyancy). The current implementation of the biogeochemical model has four active variables: nitrate (N), phytoplankton (P), zooplankton (Z), and detritus (D). A variable number of additional passive tracers may also be included.

Table 3 Input parameters expected by the MAMEBUS model code. All parameters listed in this table are chosen by the user. The sample values listed in this table are those used in the reference experiments described in Sect. 5.

4.1 Expected user inputs and options available

MAMEBUS expects a list of parameters given in Table 3 that control the physical components of the model, the model run details, and the grid setup. Other identifiers included in this model, given in Table 4, determine which internal schemes the model uses for each specific run. Furthermore, MAMEBUS expects a set of input parameters from physical tracers, forcing, diffusivity, and restoring, along with initial profiles of biogeochemical tracers that are listed in Table 5.

For the solutions shown in Sect. 5, the following initial conditions are detailed in Sect. 5.1.

Table 4 MAMEBUS numerical scheme options and descriptors.

Table 5 A table outlining the initial profiles that MAMEBUS expects during initialization. To visualize the grid locations, see Fig. 2. Each initial profile is included in “all modelTypes” unless otherwise stated. Note that N_x is the number of zonal domain points and N_z is the number of vertical domain points given in Table 3.

4.2 Model run details

The main function of the mamebus.c file has five major components and steps:

1. Calculate the time tendency of each tracer. The time step is calculated using the “tderiv” function detailed in Fig. 4. The explicit tendencies are calculated following Sect. 2.

2. Add implicit vertical diffusion and remineralization (Eq. 91).

3. Apply zonal barotropic pressure-gradient correction if the “momentumScheme” is MOMENTUM_TTW (Sect. 3.4).

4. Enforce zero tendency where relaxation time is zero (Sect. 2.5.1).

5. Write model state (Sect. 4.3).

4.3 Model data

All of the model input and output are saved in binary files. Depending on the “monitorFreq” or the frequency of output, the model will interpolate the between time steps, calculate the correct model state if necessary, and write the data to file. The following list contains all files that are written to file during the time-integration step. For each model, there is an option to include an arbitrary number of passive tracers; however, this is the standard list of tracers that are included in the indicated modelTypes.

• Residual stream function, ψ^† (all modelTypes)

• Mean stream function, $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ (all modelTypes)

• Eddy stream function, ψ^⋆ (all modelTypes)

• Temperature field (all modelTypes)

• Nitrate (NPZD model)

• Phytoplankton (NPZD model)

• Zooplankton (NPZD model)

• Detritus (NPZD model)

5 Reference solution and model validation

In this section, we present reference solutions for MAMEBUS. Below, we discuss the choice of parameters, the non-conservative forcing, and profiles of restoring. We focus predominantly on the output of a single run and plan in the future to run parameter sweeps to better understand the response of the ecosystem dynamics to the physical forcing.

5.1 Model geometry, initial conditions, and forcing

The model is configured to represent an idealized California Current System (CCS). While the model can be formulated to represent a general EBUS, we use the California Current System as a test case because this allows comparison of our results with measurements from McClatchie (2016). Note that we exclude salinity as a physical tracer; while it may be important in determining the structure of the California undercurrent (Connolly et al., 2014), we find that the main features of stratification can be well described by temperature.

A list of input fields that MAMEBUS expects is given in Table 5, with a subset illustrated in Fig. 5. The solutions shown in Sect. 5 use the following choices for these input fields. The wind-stress profile is given by

$$\begin{array}{}\text{(125)}& \mathit{\tau }\left(x\right)={\mathit{\tau }}_{\mathrm{0}}\tanh \left({\mathit{\lambda }}_{\mathit{\tau }}{\displaystyle \frac{L_x-x}{L_x}}\right),\end{array}$$

where L_x is the width of the computational profile given in Table 3, and λ_τ=4 is a tuning parameter that controls the horizontal width of the wind-stress drop-off, or wind-stress curl. We tune the offshore maximum to approximate values reported by Castelao and Luo (2018). While this is the example of wind-stress forcing we choose to use to validate our model, any form of wind-stress forcing can be defined by the user.

The topography for the reference solutions is

$$\begin{array}{}\text{(126)}& {\mathit{\eta }}_{\mathrm{b}}\left(x\right)=H-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(H-H_{\mathrm{s}}\right)\tanh \left({\displaystyle \frac{x-x_{\mathrm{t}}}{L_{\mathrm{t}}}}\right),\end{array}$$

where H is depth of the computational domain, H_s is the slope depth, x_t is the location of the continental slope in the computational domain, and L_t is the width of the continental slope given from the topographic slope parameter. All parameters are given in Table 3. The topography is tuned to represent an idealized profile of bathymetry (ETOPO5; Eto, 1988) taken from the geographic coordinates given from line 80 in the CalCOFI data (McClatchie, 2016).

Figure 4 The call tree from the main function of mamebus.c.

Figure 5 Initial temperature profile with a profile of offshore restoring which is modeled as a sponge layer on the western side of the boundary, and at the surface, there is a surface restoring to an atmospheric profile, idealized to a profile of temperature from CalCOFI. The northward wind stress is shown at the top of the figure. The white lines in the temperature field are a few lines of constant initial temperature.

The initial conditions for the tracers in the model are the initial temperature profile, including timescales and inputs for restoring, and initial conditions for the NPZD model, which are tuned to give an approximate concentration of 30 mmol m⁻³ in the deep ocean. The biogeochemical tracers are not restored in this set of reference solutions. The initial profile of temperature is shown in Fig. 5 and given by

$$\begin{array}{}\text{(127)}& \begin{array}{rl}{T}_{\text{init}}(x,z)& =T_{\text{min}}+\left(T_{\text{max}}-T_{\text{min}}\right)\\ & \cdot {\displaystyle \frac{\exp \left(\frac{z}{H^{\star }}+\mathrm{1}\right)-\exp \left(-\frac{H}{H^{\star }}+\mathrm{1}\right)}{\exp \left(\mathrm{1}\right)-\exp \left(-\frac{H}{H^{\star }}+\mathrm{1}\right)}},\end{array}\end{array}$$

where the minimum and maximum temperatures in the domain are T_min=4 ^∘C, $T_{\text{max}}=T_{\text{max}}^{\mathrm{s}}-(T_{\text{max}}^{\mathrm{s}}-T_{\text{min}}^{\mathrm{s}})x/L_x$. The maximum and minimum surface temperatures are $T_{\text{max}}^{\mathrm{s}}=\mathrm{22}$ ^∘C and $T_{\text{min}}^{\mathrm{s}}=\mathrm{18}$ ^∘C, respectively. H^⋆ is a decay scale for the temperature from the surface. This profile is tuned so that the temperature profile on the western side of the domain approximately matches the profile of temperature from CalCOFI (McClatchie, 2016) in Fig. 7. We initialize the temperature field with a small tilt in the isosurfaces to speed up the spin-up process. This same initial condition is used as the reference for temperature restoring. The timescale for restoring is given by

$$\begin{array}{}\text{(128)}& R_T^{\text{west}}(x,z)={\left({\displaystyle \frac{\mathrm{1}}{R_T^{\text{max}}}}{\displaystyle \frac{L_{\mathrm{r}}-x}{L_{\mathrm{r}}}}\right)}^{-\mathrm{1}},x<L_{\mathrm{r}},\end{array}$$

where L_r=50 km is the width of the sponge layer on the western side of the domain, and $R_T^{\text{max}}=\mathrm{30}$ d is the fastest relaxation timescale for temperature. In the surface grid boxes, the restoring timescale is given by

$$\begin{array}{}\text{(129)}& R_T^{\text{surf}}\left(x\right)=\mathrm{1}\mathrm{d},\end{array}$$

which is consistent with the formulation of Haney (1971) for surface grid box thicknesses of approximately 1 m. The restoring at the surface grid box is set to the initial profile of temperature given in Eq. (127).

The initial conditions for NPZD tracers include a constant concentration of nitrate, N_max=30 mmol m⁻³, phytoplankton, P_max=0.02 mmol m⁻³, zooplankton, Z_max=0.01 mmol m⁻³, and an initial profile of detritus of zero. This choice allows for the internal ecosystem dynamics to control the biogeochemical solutions. Finally, the cell size we choose for the phytoplankton cell is ℓ_p=1 µm. The zooplankton cell is optimized to give the optimal predator–prey length scale between the phytoplankton and zooplankton interactions, i.e.,

$$\begin{array}{}\text{(130)}& {\mathrm{\ell }}_z=\exp \left({\displaystyle \frac{\mathrm{1}}{\mathrm{0.56}}}\log \left({\displaystyle \frac{{\mathrm{\ell }}_p}{\mathrm{0.65}}}\right)\right).\end{array}$$

5.2 Isopycnal, buoyancy, and diapycnal mixing

The unresolved mesoscale and microscale mixing in the tracer evolution Eq. (9) are detailed in Sect. 2.3.1 and 2.3.2, respectively. The diapycnal diffusivities are independent of wind stress and are determined by user-input mixed layer depth and maximum magnitudes. The isopycnal and buoyancy diffusivities are time-invariant fields whose spatial structure is prescribed by the user.

In our model reference configuration, the eddy and buoyancy diffusivities are functions of the baroclinic radius of deformation – the preferential length scale at which baroclinic instability occurs and closest to the fastest growing mode in the Eady model (Eady, 1949). In MAMEBUS, these diffusivities also exponentially decreases with depth. There are choices for more sophisticated parameterizations of eddy transfer across continental slopes (Wang and Stewart, 2018, 2020), but in this current version of the model, we opt for a simpler description. For example, the buoyancy diffusivity coefficient is defined as

$$\begin{array}{}\text{(131)}& {\mathit{\kappa }}_{\text{gm}}={\displaystyle \frac{{\mathit{\kappa }}_{\text{gm}}^{\mathrm{0}}{R}_d}{H_{\text{max}}}}\exp \left(\mathit{\lambda }{\displaystyle \frac{z}{{\mathit{\eta }}_{\mathrm{b}}}}\right),\end{array}$$

where λ<1 is a tuning coefficient that allows for adjustment of the depth of the exponential profile of diffusivity, H_max is the maximum depth of the topography offshore, and η_b is the depth of the topography. For all solutions shown in this section, λ=0.25. Note that for this formulation, we assume that z<0. The maximum buoyancy diffusivity is ${\mathit{\kappa }}_{\text{gm}}^{\mathrm{0}}=\mathrm{1200}$ m² s⁻¹. Furthermore, κ_iso=2κ_gm, following Smith and Marshall (2009) and Abernathey and Marshall (2013). The isopycnal and buoyancy diffusivity profiles are shown in the left and center panels of Fig. 6, respectively.

The diapycnal diffusivities shown in the right panel of Fig. 6, with structure function described in Eq. (17), are set so that the maximum diffusivity in the mixed layers are given by ${\mathit{\kappa }}_{\text{sml}}^{\mathrm{0}}={\mathit{\kappa }}_{\text{bbl}}^{\mathrm{0}}=\mathrm{0.1}$ m² s⁻¹; otherwise, the ambient diffusivity in the interior is given by ${\mathit{\kappa }}_{\text{bg}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹. In the case where the mixed layers join at the eastern edge of the domain, the profiles of diffusivity are simply added.

5.3 Model validation

We run the reference solutions of MAMEBUS for 25 model years, with initial conditions and physical forcing described in Sect. 5.1. We validate the model against observations of temperature, nitrate, and chlorophyll a concentration in the euphotic zone, based on observations from the CalCOFI program (McClatchie, 2016). For this comparison, we interpolate a typical CalCOFI section (line 80) to a σ coordinate grid with realistic topography from the ETOPO database (Eto, 1988). We chose to validate our model with a single transect of from CalCOFI instead of several transects along the same line because averaging over time smooths over the deep chlorophyll maximum.

Furthermore, we prescribed a continental shelf that is deeper than in nature in order to reduce the model's computation time. Further shallowing the continental shelf is possible, but the CFL constraint imposed by the finer vertical resolution on the shelf extends the computation time.

While the continental slope is tuned to have a similar slope as observations in central California near the shelf break, the mixed layers in this model run are set to a constant depth zonally and overlap on the shelf. This choice has been made for simplicity and could be refined via zonally varying mixed layer depths to improve agreement with specific EBUSs. The well-mixed area on the shelf is an analogue to the inner shelf, albeit somewhat deeper than those found in nature (Lentz and Fewings, 2012). In our model comparison, we neglect the inner shelf region in the model and compare the solutions and starting approximately 50 km from the coast.

Figure 6 Inputs of buoyancy diffusivity (a), isopycnal diffusivity (b), and diapycnal diffusivity (c) used in the reference solution to MAMEBUS shown in Sect. 5. Note that the isopycnal and diapycnal diffusivities are shown over the entire domain, and the diapycnal diffusivities are shown over the upper 75 m of the domain to highlight the boundary layer mixing and the mixing on the eastern side of the domain on the shelf where the boundary layers merge.

The model temperature is generally in good agreement with observations for the upper ocean, reproducing sloping isotherms towards the coast, and realistic surface values. We observe a cold bias near the coast, which could be a result of the constant wind-stress curl forcing over the domain, inducing upwelling that is too strong in the model. A cold bias observed in the surface just outside the shelf, and a warm bias offshore, are likely caused by the prescription of a constant mixed layer depth, which may be too deep in the model for this particular section and time of the year.

As shown by the middle row of Fig. 7, model nitrate agrees reasonably well with observations in the upper layers, although biases remain, in particular in deeper layers. This may be caused by several factors, including biases in the cross-shore and vertical circulation, and in the cycling of inorganic nutrients and organic matter. For example, remineralization processes are simplified in the model, which does not include dissolved organic matter, and represents export by a single particle size class with a constant sinking speed that was not explicitly tuned to match nutrients.

The bottom row of Fig. 7 shows that the model captures the main features of the observed chlorophyll distribution (here calculated based on a fixed chlorophyll to phytoplankton nitrogen ratio of 4:5 mg, chl m⁻³ : mmol N m⁻³ following Furuya, 1990). High surface concentrations are reproduced near the shelf, with values decreasing further offshore. A deep chlorophyll maximum develops in the lower euphotic zone, at depths between 40 and 80 m, progressively deepening from the coastal to the oligotrophic region offshore. While these patterns are fairly realistic, we note that the very high chlorophyll concentrations observed near the shelf are missing from the model. This underestimation may be caused by the oversimplification of the ecosystem structure in the NPZD model, which only includes a single phytoplankton group, while multiple groups are likely required for a more correct representation of enhanced coastal phytoplankton biomass (Van Oostende et al., 2018). Furthermore, aspects of these differences could be caused by the idealized nature of the 2-D circulation simulated by the physical model.

Figure 7 Model validation against in situ CalCOFI data taken along line 80 (point conception) during July 2015. The column on the left shows output from the model under constant wind forcing and is averaged over the last 5 model years. The column on the right shows values taken from CalCOFI and interpolated onto a σ coordinate grid to allow for direct comparison. The dots in the panels are locations where the data are sampled. This figure shows the comparison between potential temperature, θ (a, b), nitrate (c, d), and chlorophyll concentration (e, f).

In order to compare physical solutions, we also include solutions which show the residual stream function, including the mean and eddy components in Fig. 8. The mean stream function is calculated via the momentum equations given in Sect. 2.2, whereas the eddy stream function is described in Sect. 2.3. The positive values indicate clockwise circulation, which, in this case, is indicative of eddy re-stratification opposing the mean upwelling branch (Colas et al., 2013). The negative values indicate counterclockwise circulation. Figure 8 shows that residual upwelling of waters onto the continental shelf via the bottom boundary layer, as interior transport onto the shelf is compensated by eddies. In the deep ocean (>500 m), there is a relatively strong residual overturning circulation that is likely associated with bottom intensification of the diapycnal mixing coefficient (e.g., McDougall and Ferrari, 2017).

5.4 Resolution parameter sweep

In this section, we describe the changes in solutions due to model resolution. We chose four different resolutions and explored the results. Figure 9 shows the solutions of MAMEBUS after 30 model years. Each panel in Fig. 9 shows the model state in the euphotic zone, averaged over the final 10 years of integration. All resolutions have the same setup and forcing as described in Sects. 4 and 5. The top row shows the potential temperature (θ), the middle row shows the nitrate concentration, and the bottom row shows the phytoplankton concentration. The model grid resolution increases from left to right, with the coarsest simulation run on a grid of 32 points horizontally and vertically, and the highest-resolution simulation run on a grid of 128 points horizontally and vertically.

Increasing the resolution leads to an overall shoaling of nutrients toward the surface. The largest overall change in near-slope nutrient concentration occurs when the resolution doubles from 32 to 64 horizontal points and vertical levels. Increasing the resolution beyond a 64×64 grid does not substantially change the horizontal distribution of phytoplankton. As referenced in Table 6, doubling the resolution increases the model run time by a multiple of approximately 20. Thus, while the model can practically be run at higher resolution, our tests show that intermediate resolution (64 horizontal and vertical levels) is sufficient to produce a favorable comparison with in situ data, without substantially increasing the computation time.

Figure 8 Stream functions calculated by MAMEBUS. This figure shows the residual stream function (a) ${\mathit{\psi }}^{†}=\stackrel{\mathrm{‾}}{\mathit{\psi }}+{\mathit{\psi }}^{\star }$, the mean stream function (b) as calculated in Sect. 2.2, and the eddy stream function (c) as described in Sect. 2.3. Note that positive values indicate clockwise circulation, whereas negative values indicate counter-clockwise circulation.

Table 6 A table outlining model run times of varying resolution between a computational cluster comprised of Intel Xeon E5-2650 v3 CPUs and a 2015 Mac laptop running macOS Catalina (version 10.15.7) for 20 model years, for both computing systems, the model is run on a single core. The highest resolution simulation (128×128 horizontal and vertical levels) was conducted on the cluster only due to computational constraints on a laptop.

Figure 9 This figure shows the model output of temperature, nitrate, and chl with varying resolution. The model was run for 30 years, and solutions shown are averaged over the final 30 years of the model run.

6 Discussion and future work

In this paper, we described the formulation, implementation, and main features of MAMEBUS, an idealized, meridionally averaged model of eastern boundary upwelling systems. The solutions are determined by a general evolution equations for materially conserved tracers (Sect. 2) and the fluid momentum equations under the time-dependent turbulent thermal wind (T3W) approximation (Dauhajre and McWilliams, 2018). It includes parameterizations of mesoscale eddy transfer and surface and bottom boundary layer mixing (Sect. 2.3), and a simple ecosystem formulation (Sect. 2.4). We further detailed the algorithms and discretizations implemented in the model (Sect. 3) and discussed reference model inputs and solutions (Sect. 4). Finally, we performed a preliminary validation based on observations from the California Current System, and we discussed the sensitivity of the model to horizontal and vertical resolution (Sect. 5).

MAMEBUS represents a simple, physically consistent tool in which to test and tune a variety of physical parameterizations and ecosystem model formulations. The ultimate goals of this research include exploration of physical–biogeochemical interactions in EBUS, mechanistic understanding of the factors that control cross-shore gradients in biogeochemical and ecological properties, and investigation of the processes that drive differences between distinct EBUSs.

Because of the 2-D framework, we acknowledge shortcomings to the model formulation, including physical aspects like intensification of upwelling around topographic features, for example, resulting from variations in the wind-stress curl (Castelao and Luo, 2018) or fine-scale ocean dynamics. Furthermore, while we parameterize the effect of mesoscale eddies on circulation, we do not account for submesoscale eddies on the shelf, which could play an important role in tracer transport (Dauhajre and McWilliams, 2018). We also do not explicitly represent breaking internal waves and tides on the shelf, which may play an important role in dissipating energy and mixing tracers when the water column is shallow (Lamb, 2014).

In future studies, we plan to use MAMEBUS to explore the effect of physical drivers such as wind stress, bathymetry, stratification, and eddies in controlling the zonal distribution of phytoplankton and food web processes, as informed by a size-structured ecosystem model. Furthermore, we plan to expand upon the physical framework in this paper by expanding eddy parameterizations to include the effect of submesoscale eddies on the shelf, where the mesoscale eddy activity is inhibited. An aspect of MAMEBUS that requires further investigation is the effect of meridional pressure gradients, which we neglected in our reference solutions in Sect. 5. In reality, the presence of along-shore pressure gradients may support interior across-shore transport away from the surface and bottom boundary layers, with the potential to reshape the coastal ecosystem.

With its limited computational cost, MAMEBUS can be used to investigate a wide parameter space in EBUSs and determine their sensitivity to a range of perturbations in major physical forcings, from changes in wind stress to increasing buoyancy forcing associated with climate change (Rykaczewski and Dunne, 2010; Sarmiento et al., 1998). Furthermore, by allowing coupling to a variety of biogeochemical and ecosystem models, MAMEBUS can be used to inform comprehensive regional models (Shchepetkin and McWilliams, 2005), for which computational costs preclude exhaustive sensitivity studies.

Appendix A: Decomposing mesoscale eddy advective/diffusive fluxes

In this Appendix, we discuss the partitioning of the mesoscale eddy tracer flux into components due to advection and isopycnal diffusion, used in Sect. 2.1 to derive MAMEBUS's central tracer evolution equation, Eq. (9). We show that the eddy tracer flux, $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}$, can be arbitrarily decomposed into components directed along mean buoyancy surfaces and along mean tracer surfaces. These components will later be associated with eddy advection and isopycnal stirring, respectively.

The effect of mesoscale eddies on the averaged tracer concentrations is given by the convergence of the eddy tracer flux (Eq. 8):

$$\begin{array}{}\text{(A1)}& {\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}\right|}_{\text{eddies}}=-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\right),\end{array}$$

and appears on the right-hand side of Eq. (6). Being quasi-adiabatic flows, mesoscale eddies serve to stir material tracers along isopycnal surfaces; this corresponds to an eddy tracer flux directed along buoyancy surfaces (Redi, 1982). Eddies also induce a “bolus” advective transfer of tracers, a generalized “Stokes drift” that corresponds to an eddy tracer flux directed along mean tracer isosurfaces (Gent and McWilliams, 1990). Both of these effects are routinely parameterized in general circulation models (Griffies, 1998, 2018). To partition the eddy tracer flux between isopycnal stirring and bolus advection, we therefore pose a decomposition of $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}$ into components directed along mean isopycnals and along mean tracer surfaces, respectively:

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}={\mathit{\alpha }}_c{\widehat{\mathit{\tau }}}_c+{\mathit{\alpha }}_{\mathrm{b}}{\widehat{\mathit{\tau }}}_{\mathrm{b}}.\end{array}$$

Here, ${\widehat{\mathit{\tau }}}_{\mathrm{c}}$ and ${\widehat{\mathit{\tau }}}_{\mathrm{b}}$ are unit vectors that point along mean $\stackrel{\mathrm{‾}}{c}$ surfaces and along mean $\stackrel{\mathrm{‾}}{b}$ surfaces, respectively:

$$\begin{array}{}\text{(A3)}& {\widehat{\mathit{\tau }}}_{\mathrm{c}}=\widehat{\mathit{y}}\times {\displaystyle \frac{\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\mathrm{|}\mathrm{|}}},{\widehat{\mathit{\tau }}}_{\mathrm{b}}=\widehat{\mathit{y}}\times {\displaystyle \frac{\mathrm{\nabla }\stackrel{\mathrm{‾}}{b}}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{b}\mathrm{|}\mathrm{|}}}.\end{array}$$

Note that the x components of ${\widehat{\mathit{\tau }}}_{\mathrm{b}}$ and ${\widehat{\mathit{\tau }}}_c$ are positive provided that $\stackrel{\mathrm{‾}}{b}$ and $\stackrel{\mathrm{‾}}{c}$ increase monotonically upward. By taking the vector cross products ${\widehat{\mathit{\tau }}}_{\mathrm{c}}\times \text{Eq. (A2)}$ and ${\widehat{\mathit{\tau }}}_{\mathrm{b}}\times \text{Eq. (A2)}$, we can solve for the vector lengths α_c and α_b:

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{c}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\times {\widehat{\mathit{\tau }}}_{\mathrm{b}}}{{\widehat{\mathit{\tau }}}_{\mathrm{c}}\times {\widehat{\mathit{\tau }}}_{\mathrm{b}}}},{\mathit{\alpha }}_{\mathrm{b}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\times {\widehat{\mathit{\tau }}}_{\mathrm{c}}}{{\widehat{\mathit{\tau }}}_{\mathrm{b}}\times {\widehat{\mathit{\tau }}}_{\mathrm{c}}}}.\end{array}$$

Then, using Eqs. (A2)–(A4), we write Eq. (A1) as

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{c}}{\partial t}}\right|}_{\text{eddies}}& =-\left(\mathrm{\nabla }\times {\displaystyle \frac{{\mathit{\alpha }}_c}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\mathrm{|}\mathrm{|}}}\widehat{\mathit{y}}\right)\cdot \mathrm{\nabla }\stackrel{\mathrm{‾}}{c}-\mathrm{\nabla }\\ & \cdot \left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{b}}}{(\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\cdot {\widehat{\mathit{\tau }}}_{\mathrm{b}})}}(\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\cdot {\widehat{\mathit{\tau }}}_{\mathrm{b}}){\widehat{\mathit{\tau }}}_{\mathrm{b}}\right).\end{array}\end{array}$$

The first term on the right-hand side of Eq. (A5) takes the form of an advection operator, in which we can identify the eddy stream function

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\mathit{\psi }}^{\star }={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{c}}}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\mathrm{|}\mathrm{|}}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\cdot \mathrm{\nabla }b}{\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\times \mathrm{\nabla }b}}.\end{array}$$

Note that this definition is ill-defined in the limit $\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\times \mathrm{\nabla }b\to \mathrm{0}$; in this limit, ${\widehat{\mathit{\tau }}}_{\mathrm{b}}$ and ${\widehat{\mathit{\tau }}}_{\mathrm{c}}$ are parallel, the eddy tracer flux is purely advective, and the stream function becomes

$$\begin{array}{}\text{(A7)}& {\mathit{\psi }}^{\star }={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\times \mathrm{\nabla }\stackrel{\mathrm{‾}}{c}}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}}}.\end{array}$$

The second term on the right-hand side of Eq. (A5) has been written in the form of the divergence of a flux along mean buoyancy surfaces, with the isopycnal gradient operator (see Sect. 2.1) appearing explicitly as

$$\begin{array}{}\text{(A8)}& {\mathrm{\nabla }}_{\parallel }\equiv (\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\cdot {\widehat{\mathit{\tau }}}_{\mathrm{b}}){\widehat{\mathit{\tau }}}_{\mathrm{b}}.\end{array}$$

We can then identify the isopycnal diffusivity κ_iso as

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{\mathit{\kappa }}_{\text{iso}}& =-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{b}}}{(\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\cdot {\widehat{\mathit{\tau }}}_{\mathrm{b}})}}=-(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\cdot \mathrm{\nabla }\stackrel{\mathrm{‾}}{c}){\displaystyle \frac{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{b}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{b}\times \mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}}}\\ & =-{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{c}^{\prime }}\cdot \mathrm{\nabla }\stackrel{\mathrm{‾}}{c}}{\mathrm{|}\mathrm{|}\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }}},\end{array}\end{array}$$

where θ is the angle between the vectors $\mathrm{\nabla }\stackrel{\mathrm{‾}}{b}$ and $\mathrm{\nabla }\stackrel{\mathrm{‾}}{c}$.

While the above derivation is general, for application in MAMEBUS, we must make assumptions about the eddy tracer fluxes. Specifically, we assume (i) that approximately identical eddy stream functions ψ^⋆ advect each different model tracer, (ii) that the isopycnal diffusivity is positive (i.e., that eddy tracer fluxes are always directed down the mean tracer gradients), and (iii) that the isopycnal diffusivity is approximately equal for different model tracers. These assumptions are satisfied in the limit of small-amplitude fluctuations (u^′ and c^′) (Plumb, 1979; Plumb and Ferrari, 2005b).

Appendix B: Derivation of time-variable Adams–Bashforth methods

For a given tracer defined with an associated time tendency equation of the form

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\partial c}{\partial t}}=f(t,c(t\left)\right).\end{array}$$

We integrate Eq. (B1) in time from $[t_{n+\mathrm{2}},t_{n+\mathrm{1}}]$:

$$\begin{array}{}\text{(B2)}& \underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}{\displaystyle \frac{\partial c}{\partial t}}\mathrm{d}\mathit{\tau }=\underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}f(\mathit{\tau },c(\mathit{\tau }\left)\right)\mathrm{d}\mathit{\tau }.\end{array}$$

By the fundamental theorem of calculus,

$$\begin{array}{}\text{(B3)}& c\left(t_{n+\mathrm{2}}\right)-c\left(t_{n+\mathrm{1}}\right)=\underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}f(\mathit{\tau },c(\mathit{\tau }\left)\right)\mathrm{d}\mathit{\tau }.\end{array}$$

We interpolate the right-hand side using a Lagrangian polynomial of the form

$$\begin{array}{}\text{(B4)}& p\left(\mathit{\tau }\right)={\displaystyle \frac{\mathit{\tau }-t_{n+\mathrm{1}}}{t_n-t_{n+\mathrm{1}}}}f(t_n,c(t_n\left)\right)+{\displaystyle \frac{\mathit{\tau }-t_n}{t_{n+\mathrm{1}}-t_n}}f(t_{n+\mathrm{1}},c(t_{n+\mathrm{1}}\left)\right).\end{array}$$

Then using Eq. (B4) and substituting it into Eq. (B3), we have

$$\begin{array}{rl}\underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}f(\mathit{\tau },c(\mathit{\tau }\left)\right)\mathrm{d}\mathit{\tau }& =\underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}p\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\\ & =\underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}\left({\displaystyle \frac{\mathit{\tau }-t_{n+\mathrm{1}}}{t_n-t_{n+\mathrm{1}}}}f(t_n,c(t_n\left)\right)\right.\\ & \left.+{\displaystyle \frac{\mathit{\tau }-t_n}{t_{n+\mathrm{1}}-t_n}}f(t_{n+\mathrm{1}},c(t_{n+\mathrm{1}}\left)\right)\right)\mathrm{d}\mathit{\tau }\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{(\mathit{\tau }-t_{n+\mathrm{1}}{)}^{\mathrm{2}}}{t_n-t_{n+\mathrm{1}}}}f(t_n,c(t_n\left)\right)\right.\\ & {\left.+{\displaystyle \frac{(\mathit{\tau }-t_n{)}^{\mathrm{2}}}{t_{n+\mathrm{1}}-t_n}}f(t_{n+\mathrm{1}},c(t_{n+\mathrm{1}}\left)\right)\right]}_{t_{n+\mathrm{1}}}^{t_{n+\mathrm{2}}}.\end{array}$$

Defining $\mathrm{\Delta }{t}_{n+\mathrm{1}}=t_{n+\mathrm{1}}-t_n$, we obtain

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}\underset{t_{n+\mathrm{1}}}{\overset{t_{n+\mathrm{2}}}{\int }}f(\mathit{\tau },c(\mathit{\tau }\left)\right)\mathrm{d}\mathit{\tau }& ={\displaystyle \frac{\mathrm{\Delta }{t}_{n+\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }{t}_{n+\mathrm{1}}}}\left(\mathrm{2}f(t_{n+\mathrm{1}},c(t_{n+\mathrm{1}}\left)\right)\mathrm{\Delta }{t}_{n+\mathrm{1}}\right.\\ & +f(t_{n+\mathrm{1}},c(t_{n+\mathrm{1}}\left)\right)\mathrm{\Delta }{t}_{n+\mathrm{2}}\\ & \left.-f(t_n,c(t_n\left)\right)\mathrm{\Delta }{t}_{n+\mathrm{2}}\right).\end{array}\end{array}$$

Substituting Eq. (B3) into Eq. (B5) yields the full ABII time-stepping scheme, given by Eq. (96).

For higher-order AB methods, we consider a sth-order Lagrangian polynomial of the form

$$\begin{array}{}\text{(B6)}& {\displaystyle }p\left(\mathit{\tau }\right)& {\displaystyle }=\sum _{m=\mathrm{0}}^{s-\mathrm{1}}{p}_m\left(\mathit{\tau }\right)f(t_{n+m},c(t_{n+m}\left)\right)\text{(B7)}& {\displaystyle }{p}_m\left(\mathit{\tau }\right)& {\displaystyle }=\prod _{\begin{array}{c}l=\mathrm{0}\\ l\ne m\end{array}}^{s-\mathrm{1}}{\displaystyle \frac{\mathit{\tau }-t_{n+l}}{t_{n+m}-t_{n+l}}},\end{array}$$

where setting s=3 as the number of known points in the interpolating polynomial results in the ABIII method. Then, s−1 is the degree of the polynomial. The general form of higher-order AB methods is

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}& c\left(t_{n+s}\right)-c\left(t_{n+s-\mathrm{1}}\right)\\ & =\underset{t+s-\mathrm{1}}{\overset{t+s}{\int }}\sum _{m=\mathrm{0}}^{s-\mathrm{1}}{p}_m\left(\mathit{\tau }\right)f(t_{n+m},c(t_{n+m}\left)\right)\mathrm{d}\mathit{\tau }.\end{array}\end{array}$$

The algebra to solve for the full discrete form of the ABIII method follows the derivation of the ABII method above. The solution to the integration in Eq. (B8) is given by Eq. (95).

Appendix C: Comparison of boundary layer parameterizations with Ferrari et al. (2008)

Our representation of eddy advection and isopycnal stirring in the surface mixed layer (SML) and bottom boundary layer (BBL) is adapted from Ferrari et al. (2008) and is described in Sect. 2.3.2. We now directly compare our SML/BBL scheme against that of Ferrari et al. (2008) to highlight the key differences.

As discussed in the “surface mixed layer” subsection, our SML scheme leads to the same eddy stream function as that of Ferrari et al. (2008), given by Eq. (29). In contrast, the residual eddy tracer flux in the SML differs as follows:

$$\begin{array}{}\text{(C1a)}& {\displaystyle }\begin{array}{rl}({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}{)}_{\text{FMCD}}& ={\mathit{\kappa }}_{\text{iso}}\left({\stackrel{\mathrm{‾}}{c}}_x+S_{\text{sml}}{\stackrel{\mathrm{‾}}{c}}_z\right)\widehat{\mathit{x}}\\ & +{\mathit{\kappa }}_{\text{iso}}\left(S_{\text{sml}}{\stackrel{\mathrm{‾}}{c}}_x+S_{\text{sml}}{S}_{\text{int}}{\stackrel{\mathrm{‾}}{c}}_z\right)\widehat{\mathit{z}},\end{array}\text{(C1b)}& {\displaystyle }\begin{array}{rl}{\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}& ={\mathit{\kappa }}_{\text{iso}}\left({\stackrel{\mathrm{‾}}{c}}_x+S_{\text{sml}}{\stackrel{\mathrm{‾}}{c}}_z\right)\widehat{\mathit{x}}\\ & +{\mathit{\kappa }}_{\text{iso}}\left(S_{\text{sml}}{\stackrel{\mathrm{‾}}{c}}_x+S_{\text{sml}}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{c}}_z\right)\widehat{\mathit{z}},\end{array}\end{array}$$

where “FMCD” denotes the formulation of Ferrari et al. (2008). Thus, our Eq. (C1b) differs from Eq. (C1a) only by the replacement S_int by S_sml in the vertical eddy residual tracer flux. A drawback of using S_int is that typically the vertical buoyancy gradient is very small in the SML, so the form (Eq. C1a) may not be numerically stable. Ferrari et al. (2008) propose a modification of the vertical component of the tracer residual eddy flux to avoid dividing by small vertical buoyancy gradients in the mixed layer:

$$\begin{array}{}\text{(C2)}& \begin{array}{rl}& ({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{c}{)}_{\text{FMCD}}\cdot \widehat{\mathit{z}}\\ & =-{\mathit{\kappa }}_{\text{iso}}\widehat{G}\left(z\right)\left(-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_x{\stackrel{\mathrm{‾}}{b}}_z}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{c}}_x+{\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_x^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{c}}_z\right)\widehat{\mathit{z}}.\end{array}\end{array}$$

However, this alternative breaks the symmetry of the diffusion tensor, and requires the introduction of an additional vertical structure function, $\widehat{G}\left(z\right)$. Our formulation, Eq. (C3b), retains the symmetry of the stress tensor and preserves continuity of the vertical flux and its derivative at $z=-H_{\text{sml}}$ with the same structure function G_sml(z) (see Sect. 2.3.2). It is also simpler to implement, as both the stream function (Eq. 29) and the residual eddy flux (Eq. C1b) can be written succinctly in terms of the effective slope (Eq. 28).

Another difference between our formulation and that of Ferrari et al. (2008) arises in the eddy stirring of buoyancy in the SML:

$$\begin{array}{}\text{(C3a)}& {\displaystyle }({\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{b}{)}_{\text{FMCD}}={\mathit{\kappa }}_{\text{iso}}\left({\stackrel{\mathrm{‾}}{b}}_x+S_{\text{sml}}{\stackrel{\mathrm{‾}}{b}}_z\right)\widehat{\mathit{x}},\text{(C3b)}& {\displaystyle }\begin{array}{rl}{\mathit{\kappa }}_{\text{iso}}{\mathrm{\nabla }}_{\parallel }\stackrel{\mathrm{‾}}{b}& ={\mathit{\kappa }}_{\text{iso}}\left({\stackrel{\mathrm{‾}}{b}}_x+S_{\text{sml}}{\stackrel{\mathrm{‾}}{b}}_z\right)\widehat{\mathit{x}}\\ & +{\mathit{\kappa }}_{\text{iso}}\left(S_{\text{sml}}{\stackrel{\mathrm{‾}}{b}}_x+S_{\text{sml}}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{b}}_z\right)\widehat{\mathit{z}}.\end{array}\end{array}$$

The Ferrari et al. (2008) residual eddy buoyancy flux has no vertical component, whereas ours does. This impacts the rate of available potential energy release in the SML by modifying the total vertical eddy buoyancy flux:

$$\begin{array}{}\text{(C4a)}& {\displaystyle }{\stackrel{\mathrm{‾}}{w^{\prime }{b}^{\prime }}}_{\text{FMCD}}={\mathit{\kappa }}_{\text{gm}}{G}_{\text{sml}}\left(z\right){\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_x^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}}},\text{(C4b)}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{‾}}{w^{\prime }{b}^{\prime }}& ={\mathit{\kappa }}_{\text{gm}}{G}_{\text{sml}}\left(z\right){\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_x^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}}}\\ & +{\mathit{\kappa }}_{\text{iso}}{G}_{\text{sml}}\left(z\right){\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_x^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}}}\\ & \cdot \left(\mathrm{1}-G_{\text{sml}}\left(z\right){\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_z}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}}}\right).\end{array}\end{array}$$

The key difference here is that our version (Eq. C4b) typically releases more potential energy and is not strictly positive definite; if ${\stackrel{\mathrm{‾}}{b}}_z\gg {\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}$, then in principle $\stackrel{\mathrm{‾}}{w^{\prime }{b}^{\prime }}$ may be negative. This corresponds to creation of potential energy, whereas previous studies suggest that potential energy should be consistently released in the SML (Colas et al., 2013). However, by construction, the vertical derivative of this term is zero at $Z=-H_{\text{sml}}$, and in any practical case ${\stackrel{\mathrm{‾}}{b}}_z$ will be smaller than ${\stackrel{\mathrm{‾}}{b}}_z{|}_{z=-H_{\text{sml}}}$ throughout the boundary layer. This suggests that if the vertical eddy length scale λ is positive ${\stackrel{\mathrm{‾}}{b}}_{z}z<\mathrm{0}$, then our scheme releases potential energy everywhere. Note also that the GM component of the vertical eddy buoyancy flux always releases potential energy.

Finally, we compare the horizontal component of the eddy buoyancy flux in the SML:

$$\begin{array}{}\text{(C5a)}& {\displaystyle }{\stackrel{\mathrm{‾}}{u^{\prime }{b}^{\prime }}}_{\text{FMCD}}=-{\mathit{\kappa }}_{\text{gm}}{\stackrel{\mathrm{‾}}{b}}_x,\text{(C5b)}& {\displaystyle }\stackrel{\mathrm{‾}}{u^{\prime }{b}^{\prime }}=-{\mathit{\kappa }}_{\text{iso}}{\stackrel{\mathrm{‾}}{b}}_x+({\mathit{\kappa }}_{\text{gm}}-{\mathit{\kappa }}_{\text{iso}})S_{\text{sml}}{\stackrel{\mathrm{‾}}{b}}_z.\end{array}$$

Whereas the Ferrari et al. (2008) scheme preserves strict lateral downgradient diffusion, this is only true in our scheme if κ_gm=κ_iso.

Further to this comparison with the formulation of Ferrari et al. (2008), we note that the fluxes discussed above differ substantially in the BBL over sloping topography. For example, the vertical buoyancy flux becomes

$$\begin{array}{}\text{(C6)}& \stackrel{\mathrm{‾}}{w^{\prime }{b}^{\prime }}={\mathit{\kappa }}_{\text{gm}}{G}_{\text{bbl}}\left(z\right){\displaystyle \frac{{\stackrel{\mathrm{‾}}{b}}_x^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{b}}_z{|}_{z={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{b}}+H_{\text{bbl}}}}}-{\mathit{\kappa }}_{\text{iso}}{\stackrel{\mathrm{̃}}{S}}_{\mathrm{e}}{\stackrel{\mathrm{‾}}{b}}_z({\stackrel{\mathrm{̃}}{S}}_{\mathrm{e}}-S_{\text{int}}).\end{array}$$

Thus, in general, the eddy buoyancy flux will act to create potential energy ($\stackrel{\mathrm{‾}}{w^{\prime }{b}^{\prime }}<\mathrm{0}$) unless the isopycnal slope S_int is of the same sign as the bottom slope and larger in magnitude. In order to avoid this, it would be necessary to set κ_iso=0 throughout the BBL. This is a separate consideration from the orientation of the residual flux vector, which must certainly lie parallel to the topography if the diffusivity is nonzero.