BFM17 v1.0: Reduced-Order Biogeochemical Flux Model for Upper Ocean Biophysical Simulations

We present a newly developed reduced-order biogeochemical flux model that is complex and flexible enough to capture open-ocean ecosystem dynamics, but reduced enough to incorporate into highly resolved numerical simulations with limited additional computational cost. The reduced-order model, which is derived from the full 56 state variable Biogeochemical Flux Model (BFM56; Vichi et al. (2007)), follows a biological and chemical functional group approach and allows for the development of critical non-Redfield nutrient ratios. Matter is expressed in units of carbon, nitrogen, and phosphate, following 5 techniques used in more complex models. To reduce the overall computational cost and to focus on open-ocean conditions, the reduced model eliminates certain processes, such as benthic, silicate, and iron influences, and parameterizes others, such as the bacterial loop. The model explicitly tracks 17 state variables, divided into phytoplankton, zooplankton, dissolved organic matter, particulate organic matter, and nutrient groups. It is correspondingly called the Biogeochemical Flux Model 17 (BFM17). After providing a detailed description of BFM17, we couple it with the one-dimensional Princeton Ocean Model (POM) for 10 validation using observational data from the Sargasso Sea. Results show good agreement with the observational data and with corresponding results from BFM56, including the ability to capture the subsurface chlorophyll maximum and bloom intensity. In comparison to previous reduced-order models of similar size, BFM17 provides improved correlations between model output and field data, indicating that significant improvements in the reproduction of in situ data can be achieved with a low number of variables, while maintaining the functional group approach. 15


Introduction
Biogeochemical (BGC) tracers and their interactions with upper-ocean physical processes, from basin-scale circulations to millimeter-scale turbulent dissipation, are critical for understanding the role of the ocean in the global carbon cycle. These interactions cause multi-scale spatial and temporal heterogeneity in tracer distributions (Strass, 1992;Yoder et al., 1992;Jr. et al., 2001;Gower et al., 1980;Denman and Abbott, 1994;Strutton et al., 2012;Clayton, 2013;Abraham, 1998;Bees, 1998; importance of which has historically been under-appreciated even though recent observational data has indicated its potential importance as a limiting nutrient, particularly in the Atlantic Ocean (Ammerman et al., 2003). To reduce model complexity, we parameterize certain processes for which field data are lacking, such as bacterial remineralization.
In the present study, we outline, in detail, the formulation of BFM17 and its development from BFM56. We couple BFM17 to the one-dimensional Princeton Ocean Model (POM) and validate the model for open-ocean conditions using observational 75 data from the Sargasso Sea. We also compare results from BFM17 and the larger BFM56 for the same open-ocean conditions. As a result of the focus on open-ocean conditions, further assumptions have been made in obtaining BFM17 from BFM56, such as the exclusion of any representation for the benthic system and the absence of limiting nutrients such as iron and silicate.
It should be noted that the primary focus in the present study is to demonstrate the accuracy of BFM17 when compared to results from observations and BFM56; as such, here we only consider one open-ocean location (i.e., the Sargasso Sea). 80 However, the correspondence between BFM17 and the more general BFM56 provides confidence that the reduced model will also prove effective at modeling other ocean locations and conditions, and exploring the range of applicability of BFM17 remains a subject of future research. We also emphasize that relatively limited calibration of BFM17 parameters has been performed in the present study. Most parameters are set to their values used in the larger BFM56 (Vichi et al., 2007(Vichi et al., , 2013, and optimization of these parameters over a range of ocean conditions is another important direction of future research. 85 Finally, we note that other reduced-order BGC models have been calibrated using data from the Sargasso Sea, such as those developed in Levy et al. (2005), Ayata et al. (2013), Spitz et al. (2001), Doney et al. (1996), Fasham et al. (1990), Fennel et al. (2001, Hurtt and Armstrong (1996), Hurtt and Armstrong (1999), and Lawson et al. (1996). However, each of these models employs less than 10 species and none uses a CFF approach. Although some of these models employ data assimilation techniques (e.g., Spitz et al. (2001)) and produce relatively accurate results, most leave room for improvement. With a minimal 90 increase in the number and complexity of the model equations, such as those associated with tracking phosphate in addition to carbon and nitrate, and including both particulate and dissolved organic nutrient budgets, we postulate that a significant increase in model accuracy can be achieved over previous reduced-order models. Additionally, with this increase in model complexity, the disparate gap between the complexity of BGC models used in small-and global-scale studies is reduced, thereby simplifying up-and down-scaling efforts. This last point is emphasized here by the good agreement between results 95 from BFM17 and BFM56.
In the following, the zero-dimensional (0D) BFM17 model is introduced in Section 2. In Section 3, BFM17 is coupled to the one-dimensional (1D) POM physical model. A discussion of the methods used to calibrate and validate the model with field data collected within the Sargasso Sea is presented in Section 4. Model results, a skill assessment, a comparison to results from BFM56, and a brief comparison to other similar BGC models are discussed in Section 5. The 17 state equation BFM17 is a reduced-order BGC model derived from the original 56 state equation BFM56 (Vichi et al., 2007(Vichi et al., , 2013, which is based on the CFF approach. In this approach, functional groups are partitioned into living organic, nonliving organic, and non-living inorganic CFFs, and exchange of matter occurs through constituent units of carbon, nitrogen, and phosphate. To date, there are no other BGC models with this order of reduced complexity using the CFF approach, making 105 BFM17 unique and able to accurately reproduce complex ecosystem dynamics.
The reduced-order BFM17 is a pelagic model intended for upper thermocline, open-ocean, oligotrophic regions and is obtained from the more-complete BFM56 by omitting quantities and processes of lesser significance in these regions, subject to the constraint that variable internal nutrient dynamics are of continued importance. In BFM17, the living organic CFF is comprised of single phytoplankton and zooplankton living functional groups (LFGs); these two groups are the bare minimum 110 needed within a BGC model and already account for six state equations (corresponding to carbon, nitrogen, and phosphate constituents of both groups). The baseline parameters used to model phytoplankton loosely correspond to the flagellate LFG in BFM56, while the zooplankton parameters correspond to the micro-zooplankton LFG (Vichi et al., 2007(Vichi et al., , 2013. Compared to BFM56, some of the parameter values in BFM17 were altered to represent general phytoplankton and zooplankton LFGs and to improve agreement with observational data, although most parameters retain the same values as in BFM56. Dissolved 115 and particulate organic matter, each with their own partitions of carbon, nitrogen, and phosphate, are also included to account for nutrient recycling and carbon export due to particle sinking, both of which are important in upper thermocline, open-ocean, oligotrophic regions. Remineralization of nutrients is provided by parameterized bacterial closure terms, thereby reducing complexity while still maintaining critical nutrient recycling. Lastly, we track chlorophyll, dissolved oxygen, phosphate, nitrate, and ammonium, since their distributions and availability can greatly enhance or hinder important biological and chemical 120 processes. Iron is omitted from BFM17, limiting the applicability of the model in regions where iron components are important, such as the Southern Ocean. Top-down control of the ecosystem is also not included. Instead, a simple constant zooplankton mortality is used, as this is a complicated process and understanding where to add this closure and where to feed the particulate and dissolved nutrients from this process in a reduced-order model is not well understood. However, the addition of a top-down 125 closure term was tested in various ways and no major differences were observed in the model results. Consequently, it was 4 https://doi.org/10.5194/gmd-2020-134 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License. assumed that the constant mortality term was sufficient for this model, similar to other reduced-order models (Fasham et al., 1990;Lawson et al., 1996;Clainche et al., 2004). Additionally, the benthic system within BFM56 (Mussap et al., 2016) has been removed. It is assumed that within the upper thermocline of the open ocean, the ecosystem is not substantially influenced by a benthic system and any water-column influences from depth can be taken into account using boundary conditions (such 130 as those discussed in Section 4).
In summary, notable novel attributes of BFM17, in comparison with other reduced-order models of comparable complexity, are the use of (i) CFFs for living organisms, including two LFGs for phytoplankton and zooplankton, (ii) CFFs for both particulate and dissolved organic matter, and (iii) a full nutrient profile (i.e., phosphate, nitrate, and ammonium).

135
In the following, the detailed equations for each of the 17 state variables that comprise BFM17 are outlined. A summary of the 17 state variables is provided in Table 1 and a schematic of the CFFs and LFGs used in BFM17, along with their interactions, is shown in Figure 1.

Environmental parameters
The BFM17 interacts with the environment through temperature and irradiance inputs. Temperature directly affects all physi-140 ological processes and is represented in the model by introducing the non-dimensional parameter f where T * is a base temperature and Q 10,j is a coefficient that may differ for the phytoplankton and zooplankton LFGs, denoted P i and Z i , respectively. Here, the subscript i is used to denote different chemical constituents (i.e., C, N, and P) and j is used to denote different LFGs. Base values used for T * and Q 10,j are shown in Table 2. The model additionally employs a 145 temperature-dependent nitrification parameter f (T ) N , which is defined similarly to Eq.
(1) as where Q 10,N is given in Table 2 In contrast to temperature, irradiance only directly affects phytoplankton, serving as the primary energy source for phytoplankton growth and maintenance. Irradiance is a function of the incident solar radiation at the sea surface. Within BFM17, the 150 amount of photosynthetically active radiation (PAR) at any given location z is parameterized according to the Lambert-Beer model as where Q S is the short-wave surface irradiance flux, which is typically obtained from real-world measurements of the atmospheric radiative transfer, ε PAR is the fraction of PAR within Q S , λ w is the background light extinction due to water, and λ bio is 155   The living organic CFF is further subdivided into phytoplankton and zooplankton living functional groups (LFGs). the light extinction due to suspended biological particles. Values for ε PAR and λ w are given in Table 2. The extinction coefficient due to particulate matter, λ bio , is dependent on phytoplankton chlorophyll, P chl , and particulate detritus, R C , and is written as where c P and c R (2) are the specific absorption coefficients of phytoplankton chlorophyll and particulate detritus, respectively, with values given in Table 2.

Phytoplankton equations
The phytoplankton LFG in BFM17 is part of the living organic CFF and is composed of separate state variables for the constituents carbon, nitrogen, phosphorous, and chlorophyll, denoted P C , P N , P P , and P chl respectively (see also 2. Phytoplankton functional group in the living organic CFF, nitrogen constituent (state variable P N ): 3. Phytoplankton functional group in the living organic CFF, phosphorus constituent (state variable P P ): 170 4. Phytoplankton functional group in the living organic CFF, chlorophyll constituent (state variable P chl ): where the descriptions of each of the source and sink terms are provided in Table 3. The subscript "bio" on the left-hand side terms indicates that these are the total rate expressions associated with all biological processes.
For the evolution of the phytoplankton carbon constituent given by Eq. (5), gross primary production depends on the non-175 dimensional regulation factors for temperature and light as well as on the maximum photosynthetic growth rate and the phytoplankton carbon instantaneous concentration. This then gives 8 https://doi.org/10.5194/gmd-2020-134 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
where r (0) P is the maximum photosynthetic rate for phytoplankton (reported in Table 4) and f (T ) P is the temperature regulation factor for phytoplankton given by Eq.
(1). The term f (E) P is the light regulation factor for phytoplankton, which is defined 180 following (Jassby and Platt, 1976) as where E PAR is defined in Eq.
(3) and E K (the "optimal" irradiance) is given by The parameter α (0) chl is the maximum light utilization coefficient and is defined in Table 4.

185
Phytoplankton respiration is parameterized in Eq. (5) as the sum of the basal respiration and activity respiration rates, namely where b P is the basal specific respiration rate, γ P is the respired fraction of the gross primary production, the gross primary production term is given by Eq. (9), and the exudation term is defined below in Eq. (18). Values and descriptions for b P and 190 γ P are given in Table 4.
Both phytoplankton exudation and lysis, defined below, depend on a multiple nutrient limitation term f 195 f (P) P = min 1, max 0, The parameters φ  Table 4.
Phytoplankton lysis includes all mortality due to mechanical, viral, and yeast cell disruption processes, and is partitioned 200 between particulate and dissolved detritus. The internal cytoplasm of the cell is released to dissolved detritus, denoted by R (1) i , while structural parts of the cell are released to particulate detritus, denoted by the state variable R (2) i , where i = C, N, P (see also Table 1). The resulting lysis terms in Eqs. (5)-(7) are then given by 205 is the nutrient stress threshold and d P is the maximum specific nutrient-stress lysis rate, both of which are given in Table 4. The term ε (N,P) P is a fraction that ensures nutrients within the structural parts of the cell, which are less degradable, are always released as particulate detritus. This fraction is determined by the expression where φ (min) N and φ (min) P are given in Table 4.

210
If phytoplankton cannot equilibrate their fixed carbon with sufficient nutrients, this carbon is not assimilated and is instead released in the form of dissolved carbon, denoted by state variable R (1) C , in a process known as exudation. The exudation term in Eq. (5) is parameterized as where β P is the excreted fraction of gross primary production, defined in Table 4, and the gross primary production term is 215 again given by Eq. (9).
The nutrient uptake of Eqs. (6) and (7) combines both the intracellular quota (i.e., Droop) and external concentration (i.e., Monod) approaches Baretta-Bekker et al. (1997). The total phytoplankton uptake of nitrogen, represented by the combination of the two uptake terms in Eq. (6), is the minimum of a diffusion-dependent uptake rate when internal nutrient quotas are low and a rate that is based upon balanced growth needs and any excess uptake, namely where a (N) P is the specific affinity for nitrogen, h P is the half saturation constant for ammonium uptake, and φ (max) N is the maximum nitrogen quota; base values for these three parameters are given in Table 4. The net primary productivity G P in Eq. (19) is given as The specific uptake rate ν P appearing in Eq. (19) is given by It should be noted that only the sum of the two uptake terms, represented by Eq. (19), is required in the governing equation for P N given by Eq. (6). However, in the governing equations for nitrate and ammonium, denoted N (2) and N (3) (see Table   1) that will be presented later, expressions are required for the individual uptake portions from nitrate and ammonium. When 230 the total phytoplankton nitrogen uptake rate from Eq. (19) is positive, the individual portions from nitrate and ammonium are determined by where the rates on the right-hand sides are obtained from Eq. (19), and ε P is given as The preference for ammonium is defined by the saturation function s N and is given by When the phytoplankton nitrogen uptake rate from Eq. (19) is negative, however, the entire nitrogen uptake goes to the dissolved organic nitrogen pool, R (1) As with the uptake of nitrogen, phytoplankton uptake of phosphorus in Eq. (7) is the minimum of a diffusion-dependent rate and a balanced growth/excess uptake rate. This uptake comes entirely from one pool and the uptake term in Eq. (7) is correspondingly given by 11 https://doi.org/10.5194/gmd-2020-134 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
where a (P) P is the specific affinity constant for phosphorous and φ (max) P is the maximum phosphorous quota. Values for both 245 parameters are given in Table 4. If the uptake rate is negative, the entire phosphorus uptake goes to the dissolved organic phosphorus pool, R P . Predation of phytoplankton within BFM17 is solely performed by zooplankton, and each of the predation terms appearing in Eqs. (5)-(7) are equal and opposite to the zooplankton predation terms, namely 250 Equations for the zooplankton predation terms are given in the next section.
Finally, phytoplankton chlorophyll, denoted P chl with the rate equation given by Eq. (8), contributes to the definition of the optimal irradiance value in Eq. (11) and of the phytoplankton contribution to the extinction coefficient in Eq. (4). The phytoplankton chlorophyll source term in Eq. (8) is made up of only two terms: chlorophyll synthesis and loss. Net chlorophyll synthesis is a function of acclimation to light n conditions, availability of nutrients, and turnover rate, and is given by where ρ chl regulates the amount of chlorophyll in the phytoplankton cell and all other terms in the above expression have been defined previously. The term ρ chl is computed according to a ratio between the realized photosynthetic rate (i.e., gross primary production) and the maximum potential photosynthesis Geider et al. (1997), and is correspondingly given as chl is the maximum chlorophyll to carbon quota and α (0) chl is the maximum light utilization coefficient, both of which can be found in Table 4. Chlorophyll loss in Eq. (8) is simpler and is just a function of predation, where the amount of chlorophyll transferred back to the infinite sink is proportional to the carbon predated by zooplankton, giving

Zooplankton equations 265
The zooplankton LFG group in BFM17 is part of the living organic CFF and is composed of separate state variables for carbon, nitrogen, and phosphorous, denoted Z C , Z N , and Z P , respectively (see also Table 1). The governing equations for the constituent state variables are given by: 5. Zooplankton functional group in the living organic CFF, carbon constituent (state variable Z C ): 6. Zooplankton functional group in the living organic CFF, nitrogen constituent (state variable Z N ): 7. Zooplankton functional group in the living organic CFF, phosphorus constituent (state variable Z P ): where, once more, descriptions of each of the source and sink terms are provided in Table 3.

275
Zooplankton predation of phytoplankton, which appears as the first term in each of Eqs. (31)-(33), primarily depends on the availability of phytoplankton and their capture efficiency, and is expressed as where r Z is the potential specific growth rate and h (F ) Z is the Michaelis constant for total food ingestion. These parameters and their base values are included in Table 5. Here, f (T ) Z is the temperature regulating factor for zooplankton growth given by 280 Eq. (1). The total food availability can be expressed as δ Z,P P C , where δ Z,P is the prey availability of phytoplankton and is included in Table 5.
Zooplankton respiration is the sum of active and basal metabolism rates, where active respiration is the cost of nutrient ingestion, or predation. The resulting respiration rate is given by

285
where η Z is the assimilation efficiency, β Z is the excreted fraction uptake, and b Z is the basal specific respiration rate. All three parameters are included in Table 5.
The biological release terms in Eqs. (31)-(33) are the sum of zooplankton excretion, egestion, and mortality. Excretion and egestion are the portions of ingested nutrients, resulting from predation, that have not been assimilated or used for respiration.
Zooplankton mortality is parameterized as the sum of a constant mortality rate and an oxygen-dependent regulation factor 290 given by where O represents the oxygen constituent of the dissolved gas in the inorganic CFF and h O is the half saturation coefficient for chemical processes given in Table 6. The total biological release is then partitioned into particulate and dissolved organic matter, giving 13 https://doi.org/10.5194/gmd-2020-134 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
Z is the oxygen dependent specific morality rate. Base values for each parameter are given in Table 5.
The zooplankton also excrete into the nutrient pools of phosphate, N (1) , and ammonium, N (3) . These effects are represented 300 by the final terms of Eqs. (32) and (33), which are parameterized by where ν (N) Z and ν Z are specific rate constants and ϕ (opt) N and ϕ (opt) P are the optimal zooplankton quotas for nitrogen and phosphorous, respectively. All four parameters are included in Table 5.

Dissolved organic matter equations
The governing equations for the three constituents of dissolved organic matter are given by: 8. Dissolved matter in non-living organic CFF, carbon constituent [state variable R (1) C ]: 10. Dissolved matter in non-living organic CFF, phosphorus constituent [state variable R (1) P ]: All terms except for the last terms in each of these equations, representing remineralization, have been defined in previous sections. Remineralization of dissolved organic matter by bacteria is parameterized within BFM17 as a rate that is proportional 315 to the local concentration of that dissolved constituent. In Eq. (41), remineralization is parameterized as α (sinkC) is a constant that controls the rate at which dissolved carbon is remineralized and returned to the pool of carbon; this constant is given in Table 6. In Eqs. (42) and (43), remineralization is represented by the parameters ζ N (3) and ζ N (1) , which are the specific remineralization rates of dissolved ammonium and phosphate, respectively. These rates are also included in Table 6 2.1.5 Particulate organic matter equations 320 The governing equations for the three constituents of particulate organic matter are given by: 11. Particulate matter in non-living organic CFF, carbon constituent [state variable R (2) C ]: Once again, all terms except for the final remineralization terms in each equation have been defined in previous sections.
Remineralization of particular organic matter by bacteria is parameterized within BFM17 as a rate that is proportional to the local concentration of that particulate constituent. In Eq. (44), remineralization is parameterized by α (sinkC) is a constant that controls the rate at which the particulate carbon is remineralized. The base value for this constant is provided in Table 6. The parameters ξ N (3) and ξ N (1) are the specific remineralization rates of particulate ammonium and phosphate, respectively. The specific remineralization rates for particulate organic matter are also presented in Table 6.

Dissolved gas and nutrient equations
The only dissolved gas resolved by BFM17 is oxygen, O, (carbon dioxide is treated as an infinite source/sink) and the dissolved 335 nutrients in the model are phosphate, N (1) , nitrate, N (2) , and ammonium, N (3) (see also 340 15. Dissolved nutrient in the inorganic CFF, phosphate constituent (state variable N (1) ): 16. Dissolved nutrient in the inorganic CFF, nitrate constituent (state variable N (2) ): 17. Dissolved nutrient in the inorganic CFF, ammonium constituent (state variable N (3) ): Aeration of the surface layer by wind, ∂O/∂t| wind , is parameterized as described in Refs. (Wanninkhof, 1992(Wanninkhof, , 2014. In a zero-dimensional model it is a source term for dissolved oxygen and so belongs in Eq. (47). However, in any model of one dimension or more it should be treated as a surface boundary condition for dissolved oxygen and so belongs in Eq. (69)  to units of oxygen and nitrogen, respectively. All terms in the above equations have been defined in previous sections, except for nitrification. Nitrification is a source term for nitrate and is parameterized as a sink of ammonium and oxygen as where Λ (nit) N (3) is the specific nitrification rate, given in Table 6. The terms f

Zero-Dimensional Test of the BFM17
As an initial test of BFM17, the model was integrated in a 0D (i.e., time only) test for 10 years using sinusoidal forcing for the temperature (in units of • C), salinity (psu), 10 m wind-speed (m s −1 ), and PAR (W m −2 ) cycles. This forcing is implemented as

365
Note that, in the 0D framework, the wind forcing does not constrain the biogeochemical dynamics, but does play a role in oxygen exchange with the atmosphere, defined according to Wanninkhof (Wanninkhof, 1992(Wanninkhof, , 2014. Initial values for chlorophyll, oxygen, phosphate, and nitrate were taken to be similar to values from the observational data, with P chl = 0.2 mg Chl-a m −3 , O = 230 mmol O 2 m −3 , N (1) = 0.06 mmol P m −3 , and N (2) = 1.0 mmol N m −3 ). Phytoplankton carbon was calculated using the maximum chlorophyll to carbon ratio, θ chl in Table 4. Initial values for zooplankton carbon, 370 dissolved carbon, and particulate organic carbon were assumed to be the same as the phytoplankton carbon. Ammonium was assumed to have the same initial concentration as phosphate. All other constituents were calculated using their respective optimal ratios in Tables 4 and 5.

Coupled Physical-Biogeochemical Flux Model
In the following sections, we describe the coupled 1D physical and biogeochemical model used for the comparisons with the 380 observational data. The coupled physical and biogeochemical model is a time-depth model that integrates in time the generic equation for all biological state variables, denoted A j , given by where the first term on the right-hand side accounts for sources and sinks within each species due to biological and chemical reactions. Although the BFM17 formulation and model results are the primary focus of the present study, we also perform 385 coupled physical-biogeochemical simulations using BFM56 for comparison. Equation (53) C,N,P , and its value is given in Table 7. We assume v (set) = 0 for zooplankton, since zooplankton To obtain the complete 1D biophysical model, BFM17 has been coupled with a modification of the three-dimensional Princeton Ocean Model (POM) (Blumberg and Mellor, 1987) that considers only the vertical and time dimensions; that is, the 400 evolution of the system in the (z, t) space. It is well known that the primary calibration dimension in marine ocean biogeochemistry is along the vertical direction, as shown in several previous calibration and validation exercises (Vichi et al., 2003;Triantafyllou et al., 2003;Mussap et al., 2016).
The 1D POM solver (POM-1D) is used to calculate the vertical structure of the two horizontal velocity components, denoted U and V , the potential temperature, T , salinity, S, density, ρ, turbulent kinetic energy, q 2 /2, and mixing length scale, .

405
In this model adaptation, vertical temperature and salinity profiles are imposed from given climatological monthly profiles, as previously done in Mussap et al. (2016) and Bianchi et al. (2005). The model computes only the time evolution of the horizontal velocity components, the turbulent kinetic energy and the mixing length scale, all of which are used to compute the turbulent diffusivity term, K H , required in Eq. (53). In this configuration, POM-1D is called "diagnostic" since temperature and salinity are prescribed. Furthermore, pressure effects are neglected in the density equation and the buoyancy gradients and temperature 410 are used in place of potential temperature since we consider only the upper water column.
In diagnostic mode, POM-1D solves the momentum equations for U and V given by ∂V ∂t where f = 2Ω sin φ is the Coriolis force, Ω is the angular velocity of the Earth, and φ is the latitude. The vertical viscosity K M

415
and diffusivity K H are calculated using the closure hypothesis of Mellor and Yamada (1982) as where q is the turbulent velocity and S H and S M are stability functions written as The coefficients in the above expressions are (A 1 , B 1 , A 2 , B 2 , C 1 ) = (0.92, 16.6, 0.74, 10.1, 0.08), with where ρ 0 = 1025 kg m −3 , g = 9.81 m s −2 . Following Mellor (2001), G H is limited to have a maximum value of 0.028. The equation of state relating ρ to T and S is nonlinear (Mellor, 1991) and given by 425 ρ = 999.8 + (6.8 × 10 −2 − 9.1 × 10 −3 T + 1.0 × 10 −4 T 2 − 1.1 × 10 −6 T 3 + 6.5 × 10 −9 T 4 ) T where the polynomial constants have been written only up to the first digit. For a more precise reproduction of these constants, the reader is referred to Mellor (1991). Finally, the governing equations solved to obtain the turbulence variables q 2 /2 and 430 are ∂ ∂t where K q = κ K H is the vertical diffusivity for turbulence variables, κ = 0.4 is the von Karman constant, and W = 1 + E 2 2 /κ 2 (1/|z| + 1/|z − H|) 2 with (E 1 , E 2 ) = (1.8, 1.33). In Eqs. (62) and (63), the time rate of change of the 435 turbulence quantities is equal to the diffusion of turbulence (the first term on the right hand side of both equations), the shear and buoyancy turbulence production (second and third terms), and the dissipation (the fourth term). This is a second-order turbulence closure model that was formulated by Mellor (2001) as a particular case of the Mellor and Yamada (1982) model for upper ocean mixing.
Boundary conditions for the horizontal velocities U = (U, V ) and the turbulence quantities are where τ w = C d |u w |u w is the surface wind stress, u w is the surface wind vector, C d is a constant drag coefficient chosen to be 445 2.5 × 10 −3 , and z = 0 and z = z end denote the locations of the surface and the greatest depth modeled, respectively.
For all variables except oxygen, surface boundary conditions for the coupled model variable A j are By contrast, the surface boundary condition for oxygen has the form

450
where Φ O is the air-sea interface flux of oxygen computed according to Wanninkhof (1992Wanninkhof ( , 2014. The bottom (i.e., greatest depth) boundary conditions for phytoplankton, zooplankton, dissolved organic matter, and particulate organic matter are This boundary condition was chosen since it allows removal of the scalar quantity A j through the bottom boundary of the domain. This can be seen by integrating Eq. (53) over the boundary layer depth using the boundary condition above, giving where the biological part of Eq. (53) has been neglected and the resulting temporal change in the integrated scalar A j is negative since |(W + W E )| < |v (set) |, as shown in Table 7. For oxygen, phosphate, and nitrate, the bottom boundary conditions are where λ j and A * j are the corresponding relaxation velocity and observed at-bottom boundary climatological field data value, 460 respectively, of that species. Base values for the relaxation velocities are included in Table 7. Lastly, the bottom boundary condition for ammonium is dependent on the gradient of particulate organic nitrogen as where κ N (3) is a relaxation diffusivity. In general, ammonium is not often included in observational measurements, so the gradient in particulate organic nitrogen is used to approximate the bottom boundary condition for ammonium. The relaxation 465 diffusivity for ammonium at the bottom, κ N (3) , is included in Table 7.

Study Site Description
Field data for calibration and validation of BFM17 are taken from the Bermuda Atlantic Time-series Study (BATS) (Steinberg et al., 2001) and the Bermuda Testbed Mooring (BTM) (Dickey et al., 2001) sites, which are located in the Sargasso Sea and, in contrast to many oligotrophic regimes, phosphate is the dominant limiting nutrient (Fanning, 1992;Michaels et al., 1993;Cavender-Bares et al., 2001;Steinberg et al., 2001;Ammerman et al., 2003;Martiny et al., 2013;Singh et al., 2015).

Data Processing
The region encompassing the BATS and BTM sites is characterized as an open ocean, oligotrophic region that is phosphate 480 limited. This region has thus been chosen for initial calibration and validation of BFM17 due to the prevalence of oligotrophic regimes in the open ocean and to demonstrate the ability of BFM17 to capture difficult non-Redfield ratio regimes (which occur in phosphate-limited regions). The BATS/BTM data have also been collected over many years, providing long time series for model calibration and validation.
Data from the BATS/BTM area is used in the present study for two purposes: (i) as initial, boundary, and forcing conditions 485 for the POM-1D biophysical simulations with BFM17 and BFM56, and (ii) as target fields for validation of the simulations.
In addition to the subsurface BATS data, we also use BTM surface data, such as the 10 m wind speed and PAR. For each observational quantity, we compute monthly averages over 27 years for the BATS data and 23 years (not continuous) for the BTM data. Additionally, we interpolate the BATS data to a vertical grid with 1m resolution. We subsequently smooth the interpolated data to maintain a positive buoyancy gradient, thereby eliminating any spurious buoyancy-driven mixing due to 490 interpolation and averaging. Figure 3 shows the monthly climatological profiles of temperature and salinity from the BATS data, as well as the PAR and 10 m wind speed from the BTM data. Similar processing is also performed on biological variables, which largely serve as target fields for the validation of BFM17. Figure 4 shows monthly-averaged vertical profiles of chlorophyll, oxygen, nitrate, phosphate, particulate organic nitrogen, and net primary production from the BATS data.

Inputs to the Physical Model
The physical model computes density from the prescribed temperature and salinity, and surface wind stress from the 10 m wind speed; temperature, salinity, and wind speed are all provided by the BATS/BTM data. The model also uses this data in the turbulence closure to compute the turbulent viscosity and diffusivity. This diagnostic approach eliminates any drifts in temperature and salinity that might occur due to improper parameterizations of lateral mixing in a 1D model, therefore 500 providing greater reliability. In addition to the 10 m wind speed, temperature, and salinity, BFM requires monthly varying PAR at the surface. For all the monthly mean input data sets, a correction (Killworth, 1995) is applied to the monthly averages to account for monthly mean errors due to linear interpolation to the model time step. We imposed both general circulation, W , and mesoscale eddy, W E , vertical velocities in the simulations. The imposed vertical profiles of these velocities have been adapted from Bianchi et al. (2005), where the the velocities are assumed to be 505 zero at the surface and reach their maxima near the base of the Ekman layer, which is assumed to be at or below the bottom boundary of the simulations. The general large-scale upwelling/downwelling circulation, W , is due to Ekman pumping and is correspondingly given as wherek denotes the unit vector in the vertical direction. The monthly average value and sign of the wind stress curl, ∇×τ w , for 510 the general BATS/BTM region was taken from the Scatterometer Climatology of Ocean Winds database Chelton, 2008, 2011). The monthly value of W from Eq. (74) is then assumed to be the maximum, occurring at the base of the Ekman layer, for that particular month. Given the sign of the wind stress curl for the BATS/BTM region, a negative W was calculated, indicating general downwelling processes in this region. Seasonal profiles of W are shown in Figure 3(e). Due to the prevalence of mesoscale eddies within the BATS/BTM region (Hua et al., 1985), which can provide episodic 515 upwelling of nutrients to the upper water column, we also include an additional positive upwelling vertical velocity, W E , which has a timescale of 15 days. The general profile of W E is assumed to be the same as for W , with a value of zero at the surface and a maximum value at depth. However, there is no linear interpolation between each 15-day period and the maximum magnitude of W E is randomized between 0 and 0.1 m d −1 , as shown in Figure 3(f) for each 15-day period.

520
Although the BATS/BTM data includes information on many biological variables, initial conditions for only 5 of the 17 species within BFM17 could be extracted from the data. Similarly to the temperature and salinity, the initial chlorophyll, particulate organic nitrogen, oxygen, nitrate, and phosphate were interpolated to a mesh with 1 m vertical grid spacing, averaged over the initial month of January, and smoothed vertically in space to give the initial profiles seen in Figure 5(a). The remaining 12 state variable initial conditions were determined either through the adoption of the Redfield ratio C:N:P ≡ 106:16:1 (Redfield et al.,525 2005), or assuming a reasonably low initial value. Since the 1D simulations were run to steady state over 10 years, memory of these initial states was assumed to be lost, with little effect on the results.

Validation Results
The coupled BFM17-POM1D model was run using the parameter values from Tables 4-7. The empirical values were decided on the basis of standard literature values Vichi et al. (2007Vichi et al. ( , 2003Vichi et al. ( , 2013 with some adjustments to improve agreement with 540 observational data. The simulations were allowed to run out to steady-state and multi-year monthly means were calculated as functions of depth for chlorophyll, oxygen, nitrate, phosphate, particulate organic nitrogen, and net primary production, each of which were measured at the BATS/BTM site. Figure 6 qualitatively compares the BATS data (top row) with the results of from BFM17 (middle row). The model is able to capture the initial spring bloom between January and March brought on by physical entrainment of nutrients, the corresponding 545 peak in net primary production and PON around the same time, and the subsequent subsurface chlorophyll maxima during the summer. The predicted oxygen levels are lower than observed values, while nitrate and phosphate levels are generally higher than observed values. The overall structures of oxygen, nitrate, and phosphate predicted by BFM17 are, however, quite similar to that of the BATS target fields. These trends are consistent with the trends seen in the results from BFM56 (bottom row of Figure 6), suggesting that the two models are in generally close agreement.

550
To quantitatively evaluate BFM17, a model skill assessment was performed for each target field. The same skill assessment was performed for BFM56 to compare the two models. The results are summarized by the Taylor diagram in Figure 7. This diagram can be used to assess the extent of misfit between the models and observations by showing the normalized root mean square (RMS) errors, normalized standard deviation, and the correlation coefficient between each of the model outputs and the BATS target fields.

555
The normalized RMS errors were calculated as ε rms /σ obs , where ε rms is the RMS error between the model and the observation fields and σ obs is the standard deviation of the observation field. The normalized standard deviation was calculated as σ mod /σ obs where σ mod is the standard deviation of the model fields. The normalized RMS errors, normalized standard deviation, and the correlation coefficients each give an indication of the relative similarities in amplitude, variations in amplitude, and structure of each modeled field compared to the BATS target fields, respectively. For each variable, these statistics were 560 calculated over all months and all depths shown in Figure 6.
The Taylor diagram in Figure 7 shows that BFM17 and BFM56 produce similar results. For all variables except oxygen, errors in the amplitudes are within roughly one standard deviation of the observations. Additionally, the structure of the model    Fasham et al. (1990) studies used only optimization to determine a select set of parameters. in the oxygen values. The overall slightly better agreement of BFM17 with the BATS data results from the adjustment of some model parameters in BFM17 due to the removal of specific phytoplankton and zooplankton species in favor of general LFGs, and to the parameterization of remineralization using new closure terms that were calibrated to give reasonable agreement with the observational data. By contrast, BFM56 was run using the baseline parameter values that were not adjusted to improve agreement with the observations.

580
The correlation coefficients and RMS errors for both BFM17 and BFM56 are also comparable with the Ayata et al. (2013) and Fasham et al. (1990)  These results show that, with a relatively small increase in the number of biological tracers as compared to similar models, BFM17 is generally able to increase correlation coefficient values and decrease RMS error values for each target field in comparison to similar models. Moreover, BFM17 approaches the accuracy of models that use data assimilation to improve agreement with the observations, such as the Spitz et al. (2001)   to phosphorous (bottom row) for phytoplankton (first column), dissolved organic detritus (second column), and particulate organic detritus (third column). Each field is normalized by the respective Redfield ratio.
Finally, a key benefit of the chemical functional family approach used by BFM17 is the ability of the model to predict non-Redfield nutrient ratios. Figure 8 shows the constituent component ratios normalized by the respective Redfield ratios for BFM17. The figure includes the component ratios of carbon to nitrogen, carbon to phosphorous, and nitrogen to phospho-595 rous for phytoplankton, DOM, and POM. Zooplankton nutrient ratios were not included because the parameterization of the zooplankton relaxes the nutrient ratio back to a constant value. The normalized ratio values are uniform non-unity valued fields.
Ultimately, Figure 8 shows that BFM17 is able to capture the phosphate-limited dynamics that characterize the BATS/BTM region (Fanning, 1992;Michaels et al., 1993;Cavender-Bares et al., 2001;Steinberg et al., 2001;Ammerman et al., 2003;Martiny et al., 2013;Singh et al., 2015). In particular, Figure 8 shows that all results comparing carbon or nitrogen to phos-600 phorous for BFM17 produce normalized vales greater than 1, where the normalization is carried out using the Redfield ratio (i.e., a normalized value greater than 1 indicates that the field is denominator limited). Figure 8 also shows that the ratios are To calibrate and test the model, it was coupled to the one-dimensional Princeton Ocean Model (POM-1D) and forced using field data from the Bermuda Atlantic test site area. The full 56 state variable Biogeochemical Flux Model (BFM56) was also run using the same forcing. Results were compared between the two models and all six of the BATS target fields-chlorophyll, oxygen, nitrate, phosphate, PON, and NPP-and a model skill assessment was performed, concluding that the BFM17 does well at reproducing observations and produces comparable results to BFM56. In comparison with similar studies using slightly 615 less complex models, BFM17-POM1D performs on par with, or better than, those studies.
In the future, a sensitivity study is necessary to assess the most sensitive model parameters, both in BFM17 as well as in the 1D physical model. After identification of these most sensitive model parameters, an optimization can be performed to reduce discrepancies between the BATS observation biology fields and the corresponding model output fields. Finally, BFM17 is now of a size that it can be efficiently integrated in a large-scale LES, and future work will examine model results in 3D simulations 620 of the upper ocean.
Code and data availability. Current versions of BFM17 and BFM56 are at https://github.com/marco-zavatarelli/BFM17-56/tree/BFM17-56 under the GNU General Public License version 3. The exact versions of the models used to produce the results used in this paper are archived on Zenodo at http://doi.org/10.5281/zenodo.3839984. Data and scripts used to produce all figures in this paper are archived on Zenodo at http://doi.org/10.5281/zenodo.3840562.

625
Author contributions. KMS, SK, PEH, MZ, and NP formulated the reduced order model, EFK and KEN verified and debugged the model, KMS and SK ran the model simulations and produced the results presented in the paper, KMS and PEH prepared the initial draft of the paper, and all authors edited the paper to produce the final version.
Competing interests. The authors declare that they have no competing interests or other conflicts of interest.
was supported by an ANSEP Alaska Grown Fellowship and by an NSF Graduate Research Fellowship. EFK and KEN were supported by NSF OAC-1535065. The data analyzed in this paper are available from Mendeley Data (https://data.mendeley.com). This work utilized the RMACC Summit supercomputer supported by NSF (ACI-1532235, ACI-1532236), CU-Boulder, and CSU, as well as the Yellowstone (ark:/85065/d7wd3xhc) and Cheyenne (doi:10.5065/D6RX99HX) supercomputers provided by NCAR CISL, sponsored by NSF.