A 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 ﬁsh catch. Previous studies have used models of varying complexity to study EBUS dynamics, ranging from minimal two-dimensional (2D) models to comprehensive regional and global models. An advantage of 2D models is that they are more computationally efﬁcient and easier to interpret than comprehensive regional models, but their key drawback is the lack of explicit representations of important three-dimensional processes 5 that control biology in upwelling systems. These processes include eddy quenching of nutrients and meridional transport of nutrients and heat. The authors present a Meridionally Averaged Model of Eastern Boundary Upwelling Systems (MAMEBUS) that aims at combining the beneﬁts of 2D and 3D approaches to modeling EBUSs by parameterizing the key 3D processes in a 2D 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 deﬁnes the equations that 10 describe the tracer, momentum, and biological evolution, along with physical parameterizations of eddy advection, isopycal 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 efﬁciently explore the wide space of physical and biogeochemical parameters that control the zonal variations in EBUSs.


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, 20 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 . 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. 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: where the bar indicates a meridional average. The key physical tracer that follows Equation (1) is temperature, θ, which serves 5 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, ∂ t c 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, ∂ t c 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 5 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, Here u 3 is the three-dimensional velocity vector, ∇ 3 is the three-dimensional gradient operator, and κ dia the microscale dif-10 fusivity. In (2) we have assumed that the velocity field is nondivergent, i.e. ∇ 3 · u 3 = 0. We further assume that u 3 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 simplifiy (2) by assuming that horizontal tracer gradients are small compared with vertical gradients, i.e. ∂ z c ∂ x c, ∂ 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., 15 ∂c ∂t phys ≈ −∇ 3 · (u 3 c) + ∂ ∂z κ dia ∂c ∂z . ( We now reduce the dimensionality of (3) by taking a meridional average, which we denote via an overbar, 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 20 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, 25 where primes denote perturbations from the meridional average. Taking a meridional average of (3) then yields Here we have used (5a)-(5b) and the property that perturbations vanish under the average, i.e. u = c = 0. We further define ∇ ≡ ∂ xx +∂ zẑ as the zonal/vertical gradient operator, and u = ux+wẑ 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.
In its current form, Equation (6) cannot be solved prognostically for c because it includes correlations between perturbation 5 quantities, i.e. the eddy tracer flux u c . 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 c by "eddy-induced velocity" u and diffusion of c along the mean buoyancy surfaces (see Burke et al., 2015), 10 Here ∇ denotes the gradient along mean buoyancy surfaces (see §2.2. A more detailed derivation of (8) is given in Appendix A. We additionally simplify the meridional tracer advection term by assuming that ∂v/∂y ≈ 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 c downgradient. With these simplifications, the full equation for the physical evolution of tracers is given by, The terms on the right-hand side of (9) are discussed further in the following sections: in Section 2.2 we discuss the evolution of the mean velocity u via the momentum equations, and in Section 2.3 we discuss the sub-gridscale parameterizations, i.e. eddy advection, eddy stirring and mixing.

20
To evolve a meridionally-averaged tracer c using (9), the meridionally-averaged velocity field u 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), 5 https://doi.org/10.5194/gmd-2020-173 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
Here, φ = p/ρ 0 is the dynamic pressure, where ρ 0 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 Section 2.1), i.e. that the turbulent Prandtl number (e.g. Kays, 1994) is exactly equal to one.
As in Section 2.1, we now meridionally average (10a)-(10e) to obtain evolution equations for u and v, and thus implicitly 5 also for w. This yields the following set of averaged equations: 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 (10a)-(10b)) is negligible compared to other terms in the momentum equation. However, we have retained the time-evolution terms (leftmost terms in (10a)-(10b)) to allow forward evolution of the horizontal velocity fields; if these terms were neglected then these terms would 15 need to be computed diagnostically at each time step. The resulting system is almost identical to the time-dependent turbulent thermal wind (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 (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 expicitly differs from the work done in Dauhajre and McWilliams (2018) which 20 focuses on more rapid time varying evolution on smaller scales. Together with the tracer advection equation for potential temperature (i.e. (9) with c = θ), (11a)-(11e) comprise a closed set of equations for the physical evolution of MAMEBUS.
In (11c) we have invoked the earlier assumption that ∂v/∂y ≈ 0 (see Section 2.1), such that the averaged velocity field is nondivergent in the x/z plane. This implies that the zonal/vertical velocity field can be related to a mean streamfunction ψ via 25 These relationships allow us to calculate ψ, and thus w, from u, subject to the boundary conditions Here z = η b (x) is the mean sea floor elevation.
Additional boundary conditions are required to solve (11a)-(11e) prognostically. Specifically, we require that the vertical turbulent stress in (11a)-(11b) matches the wind stress applied at the sea surface and the drag stress at the sea floor, with the 30 6 https://doi.org/10.5194/gmd-2020-173 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
latter formulated via a linear drag law. Formally, these boundary conditions are Here r is a linear drag coefficient andτ y is the meridional wind-stress.

Physical parameterizations 5
In this section we describe the parameterization of unresolved microscale mixing in the tracer evolution Equation (9) and the horizontal momentum Equations (11a)-(11b), and of mesoscale eddy advection and stirring in (9). This amounts to parameterizing the diapycnal diffusivity κ dia , the isopycnal diffusivity κ iso , and the eddy velocity u .

Diapycnal Mixing
We formulate the diapycnal mixing coefficient κ dia as a sum of four distinct contributing processes: surface mixed layer turbu-10 lence (κ 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 The terms on the right-hand side of (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 K- 15 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 sml < z < 0, is given by where the dimensionless surface mixed layer vertical coordinate σ sml = −z/H sml is defined such that 0 ≤ σ sml ≤ 1 within the 20 mixed layer. The structure function G KPP (σ sml ) is given by, following Large et al. (1994) and Troen and Mahrt (1986). The scaling factor 27/4 ensures that G KPP (σ sml ) has a maximum of 1 for 0 < σ sml < 1.
The diapycnal diffusivity in the bottom boundary layer, κ bbl , is prescibed in the same way as κ sml , but over the depth range 25 η b < z < η b + H bbl (x). Thus, analogous to (16), we prescribe where the dimensionless bottom boundary layer vertical coordinate is defined as σ bbl = (z − η b )/H bbl .
At any point in space and time at which the water column is statically unstable, i.e. when N 2 < 0, we increase the value of κ dia is increased locally to parameterize the effect of density-driven convection. That is, we prescribe κ conv following Finally, the background diapycnal mixing, κ bg (x, z), is simply prescribed as a constant background diffusivity. There are 5 other that can be used (e.g. St. Laurent et al. (2002)), but we opt for simplicity in the first version of this model.

Eddy advection and isopycnal mixing
We now discuss the formulation of the eddy advection and isopycnal mixing terms in (9). As discussed in Section 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 10 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 incrop at 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 15 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 u = (u , w ), introduced in (8), is nondivergent by construction (see Appendix A) and so we write it as where ψ is the "eddy streamfunction". This advecting streamfunction is assumed to be the same for all tracers, which is 20 accurate in the limit of small-amplitude fluctuations of the velocity and tracer fields (Plumb, 1979), and takes the form 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), 25 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) Specifically, we prescribe The formulation of the modified slopes S sml and S bbl are discussed below in Sections 2.3.2.1 and 2.3.2.2.
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 component-wise as 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.
10 wheren is a unit vector oriented perpendicular to the sea surface or sea floor. Specifically, we prescribe Thus S gm and S iso are identical everywhere above the BBL. The need for a distinction within the BBL is explained below in 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 nearinfinite 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 20 where σ sml = −z/H sml is a dimensionless vertical coordinate for the SML, as in Section 2.3.1. The corresponding eddy streamfunction (21) is identical to that of Ferrari et al. (2008), 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, (23) requires that G sml vanish at the surface: We further require that the eddy streamfunction and eddy residual tracer fluxes be continuous at the base of the SML, i.e. that S sml = S int , which requires that 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 It may be shown that continuity of ∇ · u c across z = −H sml is guaranteed if where λ sml = ∂ zz b/∂ z b| z=−Hsml is a vertical lengthscale for eddy motions at the base of the mixed layer.
The simplest form for G sml (z) that satisfies conditions (30), (31) and (33) is a quadratic function of depth,

15
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 20 effective slopes S sml to be used for S gm and S iso . The sloping sea floor requires separate BBL slopes, S bbl and S bbl , and structure functions, G bbl and G bbl to satisfy the required conditions of no volume nor tracer flux through the boundary, i.e. (23) and (26).
Analogous to the SML, we define the effective slope S bbl as where σ bbl = (z − η b (x))/H bbl (x) is the BBL vertical coordinate, as in Section 2.3.1. The eddy streamfunction in the BBL is https://doi.org/10.5194/gmd-2020-173 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
To satisfy the condition of zero volume flux through the sea floor, (23), the effective slope must vanish at z = η b (x), which requires G bbl (0) = 0.
To ensure continuity of the eddy streamfunction at the top of the BBL, we require that S bbl approach S int , i.e.
5 Finally, to ensure continuity of the eddy bolus velocity, we require that the gradient of S gm be continuous at z = η b +H bbl . This imposes a constraint analogous to (33) on G bbl , where λ bbl = b zz /b z z=η b +Hbbl is a vertical lengthscale for eddies at the top of the BBL. To satisfy (37)-(39), we select a quadratic form for the structure function G bbl (σ bbl ), However, the effective slope S bbl can no longer be used to define S iso in the BBL: (26) requires that the effective slope be aligned with the bottom slope at the sea floor, S b = ∂ x η b at z = η b . We must therefore employ a modified effective slope S bbl in the isopycnal mixing operator, as expressed in Equation (27). We define S bbl as 15 where G bbl (σ bbl ) is a modified structure function that also vanishes at the ocean bed, Continuity of the eddy tracer fluxes at the top of the BBL requires that Finally, continuity of the eddy flux divergence is enforced by To satisfy (42)-(44), we select a quadratic form for the structure function G bbl (σ bbl ),

Biogeochemical Model Formulation
The current biogeochemical model implemented in MAMEBUS is an NPZD (nutrient, phytoplankton, zooplankton, and 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 a single phytoplankton and zooplankton size classes, while preserving the option to run multiple size classes in future versions of the model. We also includes a 5 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 detrirus 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 10 P Z . 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 two-fold in the Northern California Current System, whereas observations from CalCOFI ( Figure 7) show variations by an order of magnitude between nearshore and offshore stations. Alongshore gradients in chlorophyll are observed along the coast, where 15 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 20 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 25 equations for biogeochemical tracers are: is zooplankton concentration, and D (mmol N/m 3 ) is the detritus concentration. The terms on the right-hand sides of (46a)-(46d) are explained in the following subsections.

Nutrient Uptake
Common controls on phytoplankton population are bottom-up limitation (i.e. nutrient control), and top-down grazing by zoo-5 plankton (Sarmiento and Gruber, 2006). We formulate bottom-up controls using typical choices for light-and temperaturedependent terms, and Michaelis-Menten uptake (Sarmiento and Gruber, 2006). The functional form of uptake is given by: where ϕ(I) and ϕ(T ) are light and temperature limiting functions, respectively. The light attenuation is modeled by integrating 10 the Beer-Lambert Law, following Moore et al. (2001), 13 https://doi.org/10.5194/gmd-2020-173 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.  Tang (1995) bu -0.45 Scaling parameter for uptake Tang (1995) ag 26 1/d Grazing rate Hansen et al. (1994) bg -0.4 Scaling parameter for grazing Hansen et al. (1994) ao 0.65 µm Optimal predator-prey length scale Hansen et al. (1994) bo 0.56 Scaling parameter for optimal predator-prey interaction Hansen et al. (1994) and the light-dependent uptake function is modeled following Sarmiento and Gruber (2006), The temperature component of the uptake function is, The maximum uptake rate is an allometric relationship defined as, where p is the user-determined phytoplankton size expressed as equivalent spherical diameter (ESD), and 0 = 1µm is a normalized length scale, with all allometrically defined variable listed in Table 2. While there are other options for the bases of these allometric relationships outlined in this section, (eg. 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 3 . 10

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 where the maximum grazing rate is defined by an allometric relationship defined as, where d −1 represents a "per day" quantity. We define a Gaussian distribution about an optimal grazing length-scale following Banas (2011), 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 5 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,

Mortality
Mortality closure terms often set important internal dynamics in ecosystem models (Poulin and Franks, 2010). While linear 10 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, 15 and zooplankton mortality as,

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 20 pool. In this model, we represent remineralization processes via a linear rate, i.e.: where r remin is the specific remineralization rate.
Particles sink at a constant average speed in the water column, following Equation (46d). At the bottom boundary we impose zero sinking flux, i.e. w sink = 0 at z = η b (x). Thus any nutrients that sink to the sea floor as detritus must remineralize there. 25 This allows for redistribution of nutrients by mixing within the bottom boundary layer, diffusion into the interior, and transport via upwelling onto the shelf.

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 nonconservative tracer tendency as The restoring and surface flux components of this tendency are discussed separately below.

Restoring
The restoring of a tracer is represented as an exponential decay to a prescribed, spatially-varying tracer field, c r (x, z), with time scale t r (x, z). The tracer restoring is then formulated as

Tracer fluxes
Surface fluxes are represented as a tendency in the tracer concentration in the surface gridboxes. For an arbitrary tracer c, we formulate the surface flux term as follows: 15 Here F c flux,0 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 2 ), with where C p = 4000J/ • C kg is the specific heat capacity.

MAMEBUSv1.0 Algorithm
In this section we discuss the numerical solution of the model equations presented in Section 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 5 We solve the model equations presented in Section 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) → (x, σ), where σ is a dimensionless vertical coordinate and is defined such that σ = 0 at z = 0 and σ = −1 at z = η b (x). This transformation requires a relationship between 10 z and σ via a transformation function For example, a "pure" sigma coordinate corresponds to the choice where is the merionally-averaged water column thickness. However, this is not necessarily the most practical 15 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, Here C(σ) is the stretching function, defined as wherẽ Here C andC are the bottom and surface components of the stretching function, respectively. The parameters θ s ∈ [0, 10] We now write the physical tracer evolution Equation (9) in σ-coordinates. For a given function f = f (x, z(x, σ)), we can write derivatives with respect to x and σ as Using these identities, we may write the divergence of an arbitrary vector F, with components F (x) and F (z) in thex andẑ directions, respectively, as Equation (69), combined with the definition of the mean streamfunction (12), allows us to write the mean advection term in (9) 10 as An analogous expression may be obtained for the eddy advection term, ∇ · (u c), in (9), using the definition (21) of the eddy streamfunction. Next, we apply (69) to the isopycnal mixing operator, defined by (25), in Equation (9) to obtain 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 20 Thus, the quantity S int − S σ can actually be computed more directly than the true isopycnal slope, as Finally, the σ-coordinate transformation of the vertical (quasi-diapycnal) mixing operator is To summarize, we write (9) in σ coordinates as where the tendency terms are Here we define 10 as the total advective or "residual" streamfunction (Plumb and Ferrari, 2005b), and we have added factors of ζ σ in (76a) so that the fluxes can be directly identified with the zonal velocity, u † = −ζ −1 σ ∂ψ † /∂σ| x , and the dia-σ velocity, † = ∂ψ † /∂x| σ . Note also that every derivative with respect to σ is multiplied by ζ −1 σ , and that their product ζ −1 σ ∂ σ 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. 15

Spatial discretization of the tracer evolution equation
We solve (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 c at the center of the (j, k) th grid cell, which we denote as c j,k (t). The mean, eddy and residual streamfunctions are most naturally defined at the cell corners, as this allows a straightforward calculation of the 20 residual velocities at the cell edges, 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, where z j+1/2,k+1/2 denotes the physical elevation of each gridpoint. 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 c over each grid cell. The linear slopes in the x and ζ directions around c j,k (t) are calculated via slope-limited finite differences 10 between c j,k (t) and its adjacent gridpoints, with parameter 1 < σ < 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 15 20 https://doi.org/10.5194/gmd-2020-173 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
estimates of c at each cell faces via 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 For this version of the model, the formulation of the Kurganov and Tadmor (2000) scheme considers only the maximum 10 propagation speed of the momentum, u, and excludes the internal gravity wave speed which is supported with the momentum calculation in Section 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, The advective discretization (83) requires the residual streamfunction to be known on all grid cell corners, which allows the numerical fluxes to be computed at the cell edges. The mean streamfunction ψ is computed from the mean velocity field via (12), 20 The eddy streamfunction (21) depends on the "true" slope (72) of the local b contours, which we discretize as The calculation of the derivatives with respect to x and z is described in Section 3.5. We then construct ψ on cell corners as ψ j+1/2,k+1/2 = (κ gm ) j+1/2,k+1/2 (S int ) j+1/2,k+1/2 .
The tracer tendency due to isopycnal diffusion, G iso , is discretized following the formulation of Kurganov and Tadmor (2000) for parabolic operators.
where (∂ z c) j,k , (∂ z c) j+1,k , (∂ x c) j,k , and (∂ x c) j,k+1 are computed via (81a)-(81d). In the interior, the isopycnal diffusion 5 slope S iso = S int and is calculated on cell corners via (85) and interpolated to cell faces via The diffusive tendency is then computed via straightforward finite-differencing of the H fluxes, 10 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 c j,k forward one time step ∆ t , i.e.
Here n denotes the time step number, and c denotes an estimate of c at t + ∆ t (see Section 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 15 diffusion. At each longitude, or for each j, we solve Equation (91) defines a tridiagonal matrix system of algebraic equations for the unknowns {c n+1 j,k |k = 1 . . . N z }, which is inverted using the Thomas algorithm.
Finally, the meridional advection is discretized via a straightforward upwind advection scheme, where v (c) denotes the meridional velocity on tracer points and c u denotes the upstream tracer concentration, defined as 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 streamfunction and zero normal tracer flux are applied at the domain boundaries. These

Temporal Discretization
MAMEBUS evolves the model equations forward in time using Adams-Bashforth methods (e.g. Durran, 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 are numerically stable with our scheme for the momentum equations (see Section 2.2). We then describe the constraints on the time-step imposed by the 5 Courant-Fredrichs-Lewy (CFL) condition.

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 Section 3.3.2. Consider a tracer quantity c that evolves according to 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 Equation (91)). We implement the third order Adams-Bashforth or ABIII method in this version of the model as the default option for time integration.
The first two time-steps require the lower order methods. We implement a forward Euler for the first time step and a secondorder AB scheme (defined below) for the second time step, 20 Here the notation ∆t n indicates the n th time step. Derivations for the adaptive time-stepping ABII and ABIII methods are given in Appendix B.

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 25 subsections, are satisfied.
The time step for advection of tracers is limited by the time scale associated with advective propagation across the width of a grid box (∆ x or ∆ z ). These constraints can approximately be expressed as (Durran, 2010). Here u igw is the maximum horizontal propagation speed of internal gravity waves (Chelton et al., 1998), Particulate sinking in the NPZD model is also calculated explicitly and constrains the time step via a similar CFL criterion 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 10 criterion for a Laplacian diffusion operator is (Griffies, 2018) ∆t < 1 2 where κ is a diffusion coefficient and ∆ s is the spatial grid spacing. In the horizontal (∆ s = ∆ x , the diffusion coefficients that determine the diffusive timestep are the eddy diffusion and isopycnal diffusion coefficients, when κ = κ gm and κ = κ iso respectively. 1 In the vertical (∆ s = ∆ z ), the diffusive time step is determined by the diapycnal diffusivity, κ = κ dia , and by the 15 vertical component of the eddy and isopycnal diffusion operators, κ = κ gm S 2 int and κ = κ iso S 2 int respectively (see e.g. Ferrari et al., 2008).

Discrete Momentum Equations and Barotropic Pressure Correction
In this section, we describe the discretization of the momentum equations presented in Section 2.2, specifically in Equations (11a), (11b), and (11c). To facilitate our discretization, we split the pressure φ, into barotropic and baroclinic components, 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 25 1 Note that although κ iso appears only in an advective operator in (83), this operator can be written as the divergence of a diffusion tensor (Griffies, 2018), and experience with MAMEBUS suggests that the more restrictive, diffusive formulation more accurately constrains the model time step.
non-divergent. The calculation of the explicit time-step is outlined in Section 3.3 and Appendix B. In order to be numerically consistent with the calculation of the streamfunction, the mean horizontal velocities u and v are stored on the western face of each grid cell. In Figure 2, these are labeled as the u points. The time-step calculation is shown below, noting that the explicit components of the time-step, E{·} are calculated following the ABIII methods outlined in Section 3.3, and the implicit diffusion I{·} is calculated following Equation (91).

5
Given the mean momentum at time step n, u n , we first perform the explicit component of the time-step to construct an estimate of u n+1 , denoted as u * , Note that the zonal barotropic pressure gradient, ∂ x Πŷ, 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 (102) described in Section 3.5.

10
We next compute the tendency due to vertical viscosity following equation (91), which we denote via the operator I. We thereby construct a second estimate of the velocity at time step n + 1, denoted as u * * , 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), as required by (11c). We formulate the barotropic pressure correction as Substituting Equation (105) into (104), we obtain 20 This implies that the tendency in the mean zonal velocity due to the barotropic zonal pressure gradient must serve to bring the depth-integrated zonal velocity to zero, i.e.
The calculation of the vertical integral of u * * is computed in the model using a Kahan sum (Kahan, 1965).

25
Pressure gradient calculations in sigma coordinates have been long known to produce discretization errors (Arakawa and Suarez, 1983;Haney, 1991;Mellor et al., 1994Mellor et al., , 1998. These errors arise from the misalignment of geopotential and sigma coordinate surfaces. We follow Shchepetkin and McWilliams (2003) to calculate the horizontal pressure gradient force. For numerical consistency, we also calculate horizontal buoyancy gradients, required to evaluate the isopycnal slope (see §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.

Zonal Pressure Gradients
In this section, we outline the implementation of the baroclinic zonal pressure gradient calculation used in MAMEBUS follow-5 ing Shchepetkin and McWilliams (2003). The ultimate goal of the algorithm is to calculate the following baroclinic pressure gradient at a cell center (c.f. (102)), where ρ, like the linear case, is the density anomaly, ρ 0 is a reference density, and A is the area between four adjacent buoyancy grid points (shaded area in Figure 3. Shchepetkin and McWilliams (2003) calculate the second term by implementing 10 a Lagrange polynomial reconstruction of the z and ρ fields. By Green's theorem, we write the integrated horizontal density gradient as The two algorithms differ on the treatment of the vertical integration of pressure. The Density Jacobian Algorithm first interpolates the density field onto the sigma grid and calculates the values of ρ and z along the solid lines in Figure 3, then 15 integrates to obtain the pressure. The Sigma 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 a Equation (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 where z j,k is the depth value at the cell centers, where the density tracer is located. Note that the edges of Figure 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, The vertical hyperbolic differences for ρ and z are calculated similarly using Equations (110b) and (110d). Again, the hyperbolic differences on the boundaries are not defined using Equations (111a) and (111b), so we extrapolate the hyperbolic averages of density. For example, at the western edge of the domain we define 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 Figure 3. The pressure field is calculated from the surface down. The pressure is calculated in the surface grid cells as where ζ j = 0 from the rigid lid assumption in MAMEBUS. Then the pressure is calculated at successively deeper grid levels as, 25 We then correct for the iso-σ pressure gradient introduced by the slope of the sigma coordinate grid, analogous to the continuous expression in Equation (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), Finally, we use Equations (114) and (115) to calculate the pressure gradients.

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 sigma coordinates.

10
The density gradient algorithm described below calculates the line integral about the area enclosed by the φ-points where the buoyancy gradient is located (see Figure 3). Therefore, we use the following form to calculate the buoyancy gradient, where F X j,k+1/2 are the value of the integral (Equation (117)) along the vertical sides, and F C j+1/2,k is the value of the integral along the horizontal sides. This calculation follows a similar procedure as the pressure gradient. 15 First, we calculate the elementary differences, and the hyperbolic averages in b and z, given by Equations (110a) through (111b). Then calculate the value of the integral along the upper and lower sides of the domain following, Note that this formulation is the same as Equation (115), but with buoyancy instead of pressure. Then we calculate the value 20 of the line integral along the vertical components of the cell, Shchepetkin and McWilliams (2003) write Equation 68a as, This allows us to numerically integrate the buoyancy gradient in the cell as, where, A, again is the area of the cell. At the surface, the boundary condition is given that F C j+1/2,N +1 ≡ 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,

Meridional Pressure Gradients
The alongshore pressure gradient in Equation (11b), denoted by φ Ly 0 /L y , and 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 10 following the hydrostatic relationship. We define ρ N and ρ S as the densities at the northern and southern ends of the domain, and Π y = Π N,S as the surface pressures at northern and southern ends of the domain. Then the pressure is given by, Here the along shore variations in sea surface pressure and density are both model inputs. We discretize the meridional pressure 15 gradient as, 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.

5
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. 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 are given in Table 4 which 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.

Expected User-inputs, and Options Available
For the solutions shown in Section 5, the following intial conditions are detailed in Section 5.1. 5

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 Figure 4. The explicit tendencies are calculated following Section 2.

5
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 timesteps, if necessary, calculate the correct model state, 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 these are the standard list of tracers that are included in the indicated modelTypes.

10
• Residual Streamfunction, ψ † , (all modelTypes) • Mean Streamfunction, ψ, (all modelTypes) • Eddy Streamfunciton, ψ , (all modelTypes) • Temperature field, (all modelTypes) • Nitrate, (NPZD model) 15 • Phytoplankton (NPZD model) • Zooplankton (NPZD model) • Detritus (NPZD model) 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.  (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 Figure 5. The solutions shown in Section 5 use he following choices for these input fields. The wind stress profile is given by

Model Geometry, Initial Conditions and Forcing
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. This profile is tuned to give an idealization of a wind stress curl profile 5 from maps shown in Castelao and Luo (2018).
The topography for the reference solutions is , where H is depth of the computational domain, H s is the slope depth, x t is the location of the continential 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).
The initial conditions for the tracers in the model are the intial temperature profile, including timescales and inputs for 5 restoring, and initial conditions for the NPZD model, which are tuned to give an approximate concentration of 30 mmol/m 3 in the deep ocean. The biogeochemical tracers are not restored in this set of reference solutions. The initial profile of temperature is shown in Figure 5 and given by,  Figure 7. We initialize the temperature field with a small tilt in the iso-surfaces 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 where L r = 50km is the width of the sponge layer on the western side of the domain, and R max T = 30 days is the fastest relaxation timescale for temperature. In the surface grid boxes, the restoring timescale is given by, which is consistent with the formulation of Haney (1971) for surface grid box thicknesses of approximately 1 m. 20 The initial conditions for NPZD tracers a a constant concentration of nitrate, N max = 30 mmol/m 3 , phytoplankton P max = 0.02 mmol/m 3 , zooplankton, Z max = 0.01 mmol/m 3 , and an initial profile of detritus of zero. This choice allows for the internal ecosystem dynamics to control the biogeochemical solutions. Finally, the cell sizes we choose for the phytoplankton cell is, The zooplankton cell is optimized to give the optimal predator-prey length scale between the phytoplankton and zooplankton interactions, ie, 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.

5
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 parameterzations of eddy transfer acros continential slopes 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: 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 κ 0 gm = 1200 m 2 /s. Furthermore, κ iso = 2κ gm , following Smith and Marshall (2009) and Abernathey and Marshall (2013). The isopycnal 15 and buoyancy diffusivity profiles are shown in the left and center panels of Figure 6, respectively.
The diapycnal diffusivities shown in the right panel of Figure 6, with structure function described in Equation (17), are set so that the maximum diffusivity in the mixed layers are given by, κ 0 sml = κ 0 bbl = 0.1 m 2 /s, otherwise the ambient diffusivity in the interior is given by, κ bg = 1e-5 m 2 /s. In the case where the mixed layers join at the eastern edge of the domain, the profiles of diffusivity are simply added.

Model Validation
We run the reference solutions of MAMEBUS for 25 model years, with initial conditions and physical forcing described in Section 5.1. We validate the model against observations of temperature, nitrate, and chlorophyll-a concentration in the 5 euphotic zone, based on observations from the CalCOFI program (McClatchie, 2016). For this comparison, we interpolate a typical CalCOFI section (Line 80) to a sigma 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. 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.

5
As shown by the middle row of Figure 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. In order to compare physical solutions, we also include solutions which show the residual streamfunction, including the mean and eddy components in Figure 8. The mean streamfunction is calculated via the momentum equations given in Section 20 2.2, whereas the eddy streamfunction is described in Section 2.3. The positive values indicate clockwise circulation, which, in this case, is indicative of eddy restratification 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 25 coefficient (McDougall and Ferrari, 2017, e.g.).

Model Verification (Resolution Tests)
In this section, we describe the changes in solutions due to model resolution. We chose four different resolutions, which include N x = N z = 32, 64, 96, 128, and explore the results. Figure 9 shows the solutions of MAMEBUS after 30 model years. Each panel is averaged over the final 10 model years, and has the same setup and forcing as described in Sections 4 and 5. Each  properties, and investigation of the processes that drive differences between distinct EBUS.
Because of the 2D 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 15 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 sizestructured ecosystem model. Furthermore, we plan to expand upon the physical framework in this paper by expanding eddy 5 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 Section 5. In reality, the presence of along-shore pressure gradients may supports interior acrossshore 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 10 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.
u c , 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 (8), and appears on the right-hand side of (6). Being quasi-adiabatic flows, mesoscale eddies serve to stir material tracers along 10 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 iso-surfaces (Gent and McWilliams, 1990). Both of these effects are routinely parameterized in general circulation models (Griffies, 1998(Griffies, , 2018. To partition the eddy tracer flux between isopycnal stirring and bolus advection, we therefore pose a decomposition of u c into components directed along mean isopycnals and along mean tracer surfaces, respectively, 15 u c = α cτ c + α bτ b .
Hereτ c andτ b are unit vectors that point along mean c surfaces and along mean b surfaces, respectively: Note that the x-components ofτ b andτ c are positive provided that b and c increase monotonically upward. By taking the vector cross productsτ c × (A2) andτ b × (A2), we can solve for the vector lengths α c and α b , Then, using Equations (A2)-(A4), we write Equation (A1) as The first term on the right-hand side of (A5) takes the form of an advection operator, in which we can identify the eddy streamfunction 25 ψ * = α c ∇c = u c · ∇b ∇c × ∇b .
Note that this definition is ill-defined in the limit ∇c × ∇b → 0; in this limitτ b andτ c are parallel, the eddy tracer flux is purely advective, and the streamfunction becomes The second term on the right-hand side of (A5) has been written in the form of the divergence of a flux along mean buoyancy surfaces, with the isopycnal gradient operator (see Section 2.1) appearing explicitly as We can then identify the isopycnal diffusivity κ iso as where θ is the angle between the vectors ∇b and ∇c.
While the above derivation is general, for application in MAMEBUS we must make assumptions about the eddy tracer 10 fluxes. Specifically, we assume: (i) that approximately identical eddy streamfunctions ψ advect each different model tracer, and (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).

15
For a given tracer defined with an associated time tendency equation of the form, ∂c ∂t = f (t, c(t)).
Then using B4 and substituting it into B3, we have, .
For higher order AB methods, we consider a s-th order Lagrange polynomial of the form, p m (τ )f (t n+m , c(t n+m ))dτ.
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 Equation (B8) is given by Equation (95). 15 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 Section 2.3.2. We now directly compare our SML/BBL scheme against that of Ferrari et al. (2008) to highlight the key differences.