Using Global-Scale Earth System Models for Regional Fisheries Applications

Climate change may impact ocean ecosystems through a number of mechanisms, including shifts in primary productivity or plankton community structure, ocean acidification, and deoxygenation. These processes can be simulated with global Earth system models (ESMs), which are increasingly being used in the context of fisheries management and other living marine resource (LMR) applications. However, projections of LMR-relevant metrics such as net primary production can vary widely between ESMs, even under identical climate scenarios. Therefore, the use of ESM should be accompanied by an understanding of the structural differences in the biogeochemical sub-models within ESMs that may give rise to these differences. This review article provides a brief overview of some of the most prominent differences among the most recent generation of ESM and how they are relevant to LMR application.


INTRODUCTION
Environmental conditions can affect living marine resources (LMRs) through a number of mechanisms. Physical and chemical properties such as temperature, salinity, oxygen concentration, and pH can directly influence the vital rates of many organisms. These direct effects may manifest themselves within specific populations through increasing or decreasing growth, reproduction, and mortality (Pörtner, 2012); changes in geographic distribution as populations shift toward more favorable habitat; or shifts in the phenological timing of environmentally influenced events such as phytoplankton blooms (Karp et al., 2019). Food web interactions further complicate the potential influence of even small changes within one population Fay et al., 2017;Marshall et al., 2017;Masi et al., 2018).
The integration of these environmental and ecological processes has become an explicit aim of regional ecosystembased fisheries management (EBFM) frameworks (Marshall et al., 2017;Holsman et al., 2020), and global Earth system models (ESMs) offer an enticing tool to the regional fisheries scientist to potentially quantify these many impacts of environmental change on LMRs. The representation of the biosphere within these coupled models has expanded rapidly over the past several iterations of the Intergovernmental Panel of Climate Change (IPCC) assessment reports. Most general circulation models (GCMs) that participated in the Coupled Model Intercomparison Project Phase 3 (CMIP3) (Meehl et al., 2007) coupled only physical ocean and atmosphere components and did not include any representation of ocean or land biology. By CMIP5 (Taylor et al., 2012), a dozen modeling centers included an ESM variant (i.e., a GCM that also simulates chemical and biological components of the Earth system) with some form of ocean biogeochemistry (Bopp et al., 2013). The recent simulations released under CMIP6 (Eyring et al., 2016) continue this trend, with most modeling centers further developing their biogeochemical models to quantify the impact of biogeochemical and lower-trophic-level processes on both the carbon cycle and on LMRs (Kwiatkowski et al., 2020).
However, there are several complexities involved in connecting ESM output to the regional marine ecosystem models that support LMR management. Similarly to GCMs, ESMs are designed to emphasize global scale dynamics over regional processes (Stock et al., 2011). Bias correction and/or downscaling may, therefore, be required before ESMs can sufficiently capture regional patterns of greatest relevance to LMRs, particularly in productive coastal waters (Holt et al., 2009;Brown et al., 2016;Muhling et al., 2018;Xiu et al., 2018;Echevin et al., 2020;Drenkard et al., 2021). The taxonomic diversity of phytoplankton and zooplankton in ESMs may also be insufficient to represent characteristics important to higher-order consumers, such as energy density or lipid content of various zooplankton assemblages or the relative effects of ocean acidification on different plankton communities (e.g., Rose et al., 2010;Miller et al., 2017;Gao et al., 2019). In addition, even in downscaled ESMs the spatio temporal resolution may be insufficient for linking life stage-specific, highly localized responses of higher-trophic-level LMRs to their biophysical environment (Petitgas et al., 2013;Hollowed et al., 2020). As marine ecosystem models are built to capture different aspects of these responses, they can also have different theoretical frameworks, objectives, and modeling structures, leading to difficulties in the direct comparison of projections among different models (Tittensor et al., 2018;Lotze et al., 2019). Projections from marine ecosystem models coupled to ESMs can, therefore, be widely divergent, incorporating considerable uncertainty from both model types (Lotze et al., 2019). Lastly, projections of managed resources such as fish stocks are complicated by the effects of fishing removals and management structures that interact with bottom-up drivers of stock size (Woodworth-Jefcoats et al., 2015;Barange, 2019;Lotze et al., 2019). Successfully working through these issues typically requires multidisciplinary collaborations with expertise from physical and biological oceanography up through resource assessment and management (Stock et al., 2011;Hollowed et al., 2020). Despite these challenges, ESMs offer an opportunity to provide LMR managers with projections of future ecosystem states, facilitating the development of climate-resilient management strategies.
While a number of existing studies discuss, in comprehensive detail, the predicted changes in lower trophic level variables across the suite of CMIP5 and CMIP6 ESM model simulations and the potential mechanisms underlying both the mean trends and variations between models (Bopp et al., 2013;Laufkotter et al., 2015;Kwiatkowski et al., 2020), these comparisons typically focus on large-scale patterns in nutrient cycling and are primarily targeted toward a biogeochemistry audience. For the reasons detailed above, the LMR end user typically approaches these models from a different perspective (Rose et al., 2010), and the output variables in which they are most interested may show different trends, variabilities, and skills across the spatiotemporal scales of interest (Table 1). In this study, we follow in the footsteps of these existing reviews, particularly Laufkotter et al. (2015) and Séférian et al. (2020), but from the point of view of a regional fisheries end user. The performance of each model may vary widely depending on the region of focus and the scientific question under consideration, and attempting to perform a full intercomparison of global ESM simulation output on a regional scale is well beyond the scope of this study. Instead, we highlight key structural differences within the CMIP6 suite of ESMs and discuss the ways in which these differences may lead to variability among the predicted ecological indicators most relevant for higher trophic levels. We hope that by highlighting these structural differences, we provide an entry point for LMR end users to better diagnose the potential drivers of biogeochemical and ecosystem uncertainties within their regions of interest.

ACCESS TO, USE OF, AND INTERPRETATION OF GLOBAL-SCALE ESMS
While an ESM is often referred to as a single entity, each ESM is made up of many coupled sub-models representing the different components of the Earth system. Within modern ESMs, these components typically include the atmosphere, ocean, land, ocean biogeochemistry, and sea ice. While the primitive equations governing the atmospheric and ocean components of ESMs are well established, ocean biogeochemical models are less mature, often based on empirical relationships and with far fewer observations constraining the underlying equations and parameters than in their physical counterparts. As a result, the state variables and process equations included in both global and regional scale biogeochemical models vary widely, with model complexity ranging from three-to four-box nutrient/phytoplankton/zooplankton/detritus models to complex multi resource, multiple-plankton functional type models with dozens of state variables (Friedrichs et al., 2007). Even within model intercomparison experiments, there is little standardization of prognostic biogeochemical models; for example, only gas exchange and carbonate chemistry are standardized within the CMIP6 protocol, with all other processes (e.g., primary production, nutrient cycling, and plankton community composition) left to the discretion of individual modeling centers (Orr et al., 2017). Substantial differences also exist in the parameterization of unresolved physics among models (due to their computationally limited spatial/temporal resolution) with potentially strong impacts on biogeochemistry.  (Orr et al., 2017), without any explicit simulation of biotic processes, to the full ecosystem models. Given time and resource constraints, using this entire suite of ESMs is impractical for a regional application. Therefore, the first step in using global ESMs in a regional context typically involves choosing the subset of ESMs best suited to answering the research question of interest. However, making an informed choice on this front can be a daunting task due to the high complexity of the numerous ESMs.
A variety of best practices have been suggested regarding choosing global ESMs to address regional questions, either through direct use or via statistical or dynamical downscaling (Drenkard et al., 2021). One suggestion is to choose models whose outputs span the envelope of uncertainty for a key variable Pozo Buil et al., 2021). Alternatively, one may limit use to models that faithfully capture the historical observed state of a particular quantity that is known to affect the ecosystem of interest, e.g., seasonal sea ice in the Bering Sea (Overland et al., 2011), upwelling and horizontal transport off the California coast (Combes et al., 2013;Di Lorenzo et al., 2013), or hypoxia due to riverine influence in the Gulf of Mexico (de Mutsert et al., 2016;Fennel et al., 2016). Between these two alternatives is a model choice based on emerging constraints (Eyring et al., 2019;Hall et al., 2019), which is a process-based approach that acknowledges that some models are likely better suited than others for a given climate change application but also that historical model fidelity does not necessarily translate to accurate climate sensitivity under future forcing (Pierce et al., 2009).
Consideration of the broad array of biogeochemical model specifications can be beneficial when using ESMs within an LMR context, both for informing the choice of which models can best address the relevant questions and interpreting the potential underlying mechanisms leading to inter-ESM differences in output. This is particularly important given the large influence of biogeochemical model structure and parameterization on an long-term predictions of LMR-relevant variables of ESM such as pH, oxygen, and net primary productivity. Within longterm climate change simulations, Frölicher et al. (2016) found that while uncertainties related to ESM internal variability and emissions scenario choice dominate the predicted end-of-century uncertainty range for physical ocean variables like temperature and carbon chemistry variables like pH, biological output variables like net primary production were far more sensitive to model uncertainty (i.e., uncertainty resulting from the use of different numerical formulations and parameterizations for physical and biological processes).
However, from a practical standpoint, making an informed choice regarding which ESMs to use for a regional application can be difficult. While global model documentation has become more common and more cohesive in recent years (Journal of Advances in Modeling Earth Systems, 2018-2020), and model cross-comparisons are available in the literature (e.g., Kwiatkowski et al., 2020;Séférian et al., 2020), this model documentation is often disconnected from common data access points and, therefore, difficult for end users unfamiliar with the biogeochemical literature to locate. Major modeling experiments (e.g., CMIP3, CMIP5, and CMIP6) do not include "release notes" to highlight how models have changed from one iteration to the next, and often the specific details of a model are spread over numerous publications, reflecting the iterative development of each model. The lack of centralized documentation or crosscomparison of the many ESM options available to LMR end users often presents a barrier to making an informed choice regarding model selection. Additionally, LMR researchers often encounter challenges in accessing model output. The storing and sharing of very large datasets is a challenge for the major modeling centers; their choices of which variables to archive and at what resolution are often driven by limited space and do not always meet LMR needs for high spatiotemporal resolution. Public datasets often include only a small subset of the depth-resolved ocean and biogeochemical variables at    In the past, LMR scientists have needed to establish personal relationships within major modeling centers in order to get access to non-publicly archived variables. In cases where the variables of interest were not archived, either at all or at sufficient resolution, new ESM simulations may be needed, increasing the need for direct collaboration even further. This level of collaboration is difficult to maintain with even one modeling center, much less multiple; this high barrier to access further complicates the process of making an informed selection of ESMs. In the following sections, we compare and contrast how key biogeochemical processes are represented within the CMIP6 models with explicit biogeochemical models. We choose this subset from the full list in Table 2 as a practical consideration, given that these models are most likely to be of interest to LMR end users in the coming years. We do not attempt to analyze differences in output values or how skillfully each model performs when compared to observations. Particularly when extracting results within a small region, as is typical in LMR applications, which models are most appropriate can be dependent on the exact research question being asked. For an overview of model performance across the CMIP5 and CMIP6 generations of models, we instead refer readers to Séférian et al. (2020). This study instead is intended to demystify some of the key processes being simulated within complex ESMs, allowing LMR end users to make more informed choices with regard to the models that may best suit their applications and to diagnose potential sources of model spread that influence predictions of LMR-relevant metrics.

A COMPARISON OF CMIP6 ESMS
We focus on a few common features of the CMIP6 ESMs for this intercomparison: the structure of each model (i.e., the state variables and biological processes included), the formulations used to calculate primary production, the role of temperature in mediating various biological rates, zooplankton predator-prey function responses, detrital remineralization formulations, and river runoff implementations. These particular features were chosen as key components shared across the models that are most likely to influence the metrics, such as primary and secondary production, zooplankton biomass and community composition, and trophic transfer efficiency, that LMR end users may be most interested in. These metrics are typically those that can be used to predict species biomass, recruitment, and survival; shifts in spatial distribution; and trophic interactions, particularly with respect to target species of commercial, recreational, and subsistence harvest (Tommasi et al., 2021). The highlighted dissimilarities between models may play varying roles in contributing to inter-model spread, depending on the particular region, ecosystem, and metric of interest, but we hope that highlighting these model differences will provide a starting point for the often complex task of connecting model uncertainty to specific physical and biogeochemical drivers.
We note that there are also a few processes that may be of potential interest to LMR researchers that we chose not to highlight. We do not address processes related to gas exchange or carbonate chemistry; these processes are explicitly specified as part of the CMIP6 Ocean Model Intercomparison Project FIGURE 1 | Schematic for the BLING model. Figures 1-12. Schematics for the CMIP6 ESMs. Boxes indicate state variables, and solid arrows represent fluxes between state variables; dashed lines indicate that the rate of change of one state variable is calculated proportionately to another. Colors indicate the base element for each state variable, with blue, green, pink, orange, and gray representing nitrogen, phosphorous, silicon, iron, and carbon, respectively. State variable boxes are positioned vertically based on functional role, indicating whether the functional group includes producers (and within that, subcategories of nitrogen fixers, Si-users, and CaCO3-users), consumers, or detritus. Horizontal position indicates the approximate mean size of the cells/bodies/particles represented by each state variable.
biogeochemical (OMIP-BGC) simulation protocols (Orr et al., 2017) and as such are implemented in the same way across all models within the CMIP6 suite. We also do not address certain processes that involve the coupling of the biogeochemical model components to other parts of the earth system, such as atmospheric deposition of iron. We avoid these processes primarily as a practical consideration, since they are not part of the biogeochemical models themselves, but rather a function of the other model components (such as the atmospheric and land models) that each biogeochemical model may be coupled to. We make an exception for river runoff, given its often large impact on the near-shore ecosystems that are of interest to LMR scientists.

Functional Groups and Processes
The CMIP6 models vary widely in structural complexity. This variety reflects the numerous applications for which different models were designed. Many ESMs were originally designed primarily to "close" the carbon cycle, i.e., to allow the full inventory of the oceanic processes influencing oceanic sequestration of carbon emissions. As such, their focus was less on resolving the ecosystem components themselves but rather on constraining the biologically mediated fluxes of carbon between the surface waters and deep ocean. More recently, ESMs have developed a dual role, serving both the original carbon system closure purpose and the ability to resolve energy transfer within more realistic and complex planktonic food webs. In terms of structural complexity, the CMIP6 biogeochemical models can be loosely grouped into a few categories.
At the simple end of the spectrum lies BLING, the biogeochemical model that runs within the CM4 model of GFDL (Figure 1). With only three explicitly tracked state variables and implicit treatment of the living biota, the low complexity of BLING aims to capture the impact of biology on global nutrient cycles while adding minimal computational overhead (Galbraith et al., 2010;Dunne et al., 2020). Its lack of explicitly resolved primary or secondary producers likely limits its use in the living marine resource context. The tradeoff, though, is that it can be run within higher-resolution simulations than its more complex and computationally heavy counterparts. This, along with its diagnostic estimates of biologically relevant outputs like chlorophyll and carbon-system quantities, may be useful to coastal regional applications where increased horizontal resolution outweighs the need for greater state variable resolution.
The second category of models includes traditional NPZD models that include single phytoplankton, zooplankton, and detrital compartments. This category includes CMOC of CCCMA (Figure 2), WOMBAT of CSIRO (Figure 3), and MRI-NPZD model of the Japanese Meteorological Research Institute (Figure 4). The HAMOCC model (Figure 5), variants of which are run within the models of Norwegian Climate Center and Max Planck Institute, also uses a single state variable each for phytoplankton and zooplankton but separates detritus into dissolved and particulate detritus. The OECO of MIROC model (Figure 6) uses two phytoplankton functional groups to distinguish nitrogen-fixing diazotrophs from non-nitrogenfixing producers but otherwise falls into this simpler category of model. While several models in this category consider multiple limiting macronutrients (nitrogen, phosphorous, and silicon) and micronutrients (iron), they assume constant stoichiometric ratios of these nutrients within the food web.
On the more complex end of the spectrum, the remaining models partition phytoplankton, zooplankton, and detritus into multiple state variables. This group of models includes PISCES (Figure 7), which is used in both the CNRM and IPSL models; COBALT of GFDL (Figure 8); CanOE of CCCMA (Figure 9); NOBM of NASA Goddard (Figure 10); MEDUSA of the UK (Figure 11); and MARBL of the CESM model (Figure 12). Phytoplankton within these models are separated into 2-4 functional groups; functional groups are typically defined based on their interactions with different nutrient pools, their size, and their accompanying influence on metabolic rates, and in the case of the NOBM model, by optical properties. Size is also one of the primary features used to separate zooplankton functional groups due to the role of body size in determining prey preferences, metabolic rates, and the partitioning of mortality and metabolic byproducts into the different detrital pools. Nonliving organic material is categorized by either size or lability. This more complex category of models typically allows for variable stoichiometry within the primary producer and detrital state variables; most assume a fixed stoichiometry for consumers, though the MARBL model is an exception to this. Allowing for variable vs. fixed stoichiometry has been demonstrated to significantly impact simulated trophic transfer efficiency under climate change scenarios (Kwiatkowski et al., 2019), so this additional complexity may be particularly important for LMR applications.
With an increase in the number of state variables in a model comes an increase in the complexity of connections between those state variables, and the number of parameters (each with its own uncertainty) necessary to constrain them. These connections include both the pathways that pass biomass and energy to higher trophic levels and the pathways that lead to export from or recycling of nutrients in the surface waters. A more complex model may be better able to capture the shifts in community composition and prey quality that are key to simulating upper trophic level recruitment and survival. For example, environmental conditions that favor larger phytoplankton species often lead to higher transfer efficiencies to upper trophic level species, while small phytoplankton production is often directed toward smaller grazers and the microbial loop (Azam et al., 1983;Armengol et al., 2019). These potential shifts in energy pathways are of particular interest in LMR applications that focus on plankton groups in their role as "fish food" rather than as regulators of biogeochemical pathways (Rose et al., 2010). Van Oostende et al. (2018) also found that adding an additional phytoplankton group to their model of the California Current system led to better resolution of both chlorophyll maxima and coastal hypoxia. Likewise, explicitly resolving different zooplankton size classes can increase the types of ecosystem analyses that are possible based on ESM output. For example, the GFDL COBALT model is essentially a more ecosystem-heavy version of its predecessor, TOPAZ . By adding three explicit zooplankton groups (and bacteria) in place of the implicit zooplankton parameterization found in TOPAZ, COBALT has been able to be used for a diverse set of LMR applications: to model and predict fisheries catch in LMEs Park et al., 2019), to couple to offline models to examine distributions of fish functional types (Petrik et al., 2019) and the impact of jellyfish on the global carbon cycle (Luo et al., 2020), and to examine the impacts of migration and light-dependent grazing on the vertical distribution of carbon, nutrients, and chlorophyll (Bianchi et al., 2013;Moeller et al., 2019). Greater resolution of detrital components may allow for better resolution of the different time scales over which material may be exported to deep water, transported into or out of a region of interest, or returned to the nutrient pools. These remineralization pathways can directly affect the local availability of nutrients that, in turn influence rates of primary production. This is particularly true in shallow shelf systems, where nutrients exported below the mixed layer may be reintroduced to surface waters by mixing events (Huang et al., 2015;Flynn et al., 2020). Particle export ratios and fluxes to the benthic ecosystems have also been shown to be better predictors of fisheries productivity (Friedland et al., 2012;Stock et al., 2017) and of the composition of fish communities (Van Denderen et al., 2018;Petrik et al., 2019) than net primary productivity alone. It is important to note, however, that while the state variables and processes included in more complex models often reflect our best estimation of the Earth system processes that influence both lower and upper trophic level dynamics, the added complexity does not necessarily impart greater skill (Séférian et al., 2020).

Phytoplankton Growth Limitation Terms
Primary productivity rates and phytoplankton biomass are some of the most commonly extracted ESM variables of interest to LMR end users. As the base of the food web, plankton production often determines the characteristics of the food web that can be supported in a given region. Changes to the magnitude, community composition, or timing of spring and fall blooms may impact the survival of upper trophic level species (Platt et al., 2003;Durant et al., 2007;Schweigert et al., 2013;Malick et al., 2015;Asch et al., 2019). Because the most important process controlling interior ocean pH is organic matter cycling, the processes limiting primary production also influence the pH (Lauvset et al., 2020). Acid-base reactions involved in nutrient uptake also have an impact on both pH and alkalinity (Wolf-Gladrow et al., 2007). Therefore, model structure related to the regulation of primary production may influence many of the key properties of interest to LMR end users.
Phytoplankton growth rates vary based on a number of factors, chief among them temperature, light, and availability of both the macronutrients and micronutrients necessary for growth. Within ESMs, maximum phytoplankton growth rate is typically modeled as a temperature-dependent maximum potential growth rate modulated by limitation terms (with values ranging from 0 to 1) that reduce that rate based on limited nutrient supplies and light. While this concept is consistent across all CMIP6 models, the exact implementation of each limitation term and the manner in which they interact vary from model to model. We will address variations in the temperaturedependent maximum growth rate in section 3.3 and focus here on the light-and nutrient-limitation terms. Nutrient-limitation terms typically fall into one of two broad categories. The first category calculates limitation as a function of ambient nutrient concentration. The majority of the CMIP6 ESMs use this approach for limitation terms related to the macronutrients (i.e., nitrogen and phosphate) (Tables 3, 4). Monod curves (Monod, 1942), also referred to as Michaelis-Menten curves due to their functional equivalence to the enzyme kinetics equation (Michaelis and Menten, 1913), define limitation in terms of a half-saturation constant (K X ). This is the most common nutrient-limitation term equation found in the CMIP6 ESMs. The MPI-NPZD model instead opts for an optimal uptake kinetics approach (Smith et al., 2009) that accounts for tradeoffs in nutrient encounter rates and assimilation rates; like the Monod curve, this approach assumes growth rate is a function of external nutrient concentration. The second category of limitation uses an internal cell quota model. A quota-model limitation term is based on the internally stored concentration of a particular element rather than the ambient concentration, and it allows for uptake and storage of one or more nutrients beyond immediate needs. This style of limitation term can be more appropriate when dealing with elements like silicon and iron; Fe:C and Si:C ratios can vary much more widely within cells relative to macronutrients like N and P. Within the CMIP6 suite, COBALT and PISCES use quota models for iron limitation. BLING also uses a variant on the quota model for iron limitation, adapted to infer cell quota from uptake rates. CanOE uses a quota model for both nitrogen and iron uptake.
Phytoplankton can acquire nitrogen from a variety of different organic and inorganic compounds, and the sources that contribute to the nitrogen limitation term of each model provide additional potential for model structural divergence (Table 3). Non-nitrogen-fixing phytoplankton typically acquire nitrogen from two primary forms: nitrate or ammonium. The relative use of new (nitrate-based) vs. regenerated (ammonium-based) nutrient sources for production reflects whether an ecosystem is dominated by larger phytoplankton species that pass energy up the food chain and contribute to higher export ratios of carbon from surface to deep water, or whether the ecosystem is dominated by smaller phytoplankton with production primary recycled within the microbial loop (Dunne et al., 2005). Due to its reduced state, ammonium is more readily utilized by phytoplankton for growth (Gruber, 2008). In addition, the presence of high concentrations of ammonium can directly inhibit the uptake of nitrate (O'Neill et al., 1989;Dortch, 1990;FIGURE 6 | Schematic for the OECO model. See Figure 1 caption for details. Frost and Franzen, 1992;Paulot et al., 2015). Therefore, the limitation terms of many of the CMIP6 ESMs are formulated to favor ammonium uptake over nitrate uptake and to adjust nitrate uptake in the presence of ammonium. PISCES, COBALT, CanOE, MARBL, and NOBM track the separate use of NH 4 vs. NO 3 , with total nitrogen limitation being a sum of the two individual limitation terms. COBALT and PISCES account for both ammonium preference and its inhibition on nitrate, as well as inhibition of ammonium uptake of nitrate in the presence of very high concentrations of nitrate (O'Neill et al., 1989). NOBM and CanOE impose a preference for NH 4 over NO 3 but do not apply any inhibition terms to the uptake of either nutrient. OECO, HAMOCC5.1, and HAMOCC6 track NO 3 only. CMOC and MEDUSA, meanwhile, track only a single pool of nitrogen. A number of the ESMs also account for nitrogen fixation. Nitrogen fixers are able to transform N 2 gas into fixed nitrogen (NH 4 ), and in models they are typically differentiated from other phytoplankton groups through the elimination of an N-limitation term in their growth equation. COBALT, OECO, HAMOCC6, and MARBL all simulate nitrogen-fixing diazotrophs in this manner. The NOBM model instead uses an additional factor to increase the N-limited growth rate of the cyanobacteria group when N is low, representing a shift within the functional group from non-nitrogen-fixing to nitrogen-fixing cyanobacteria. CMOC simulates only a single nutrient pool and single phytoplankton group but provides a flux of nitrogen into its surface nutrient pool based on implicit parameterizations of diazotrophic bacteria. Likewise, CanOE parameterizes nitrogen fixation as an input source to the ammonium field, dependent on light, temperature, and iron concentration and inhibited by high inorganic nitrogen.
Phytoplankton growth rates can also be very sensitive to even small variations in light and light sensitivity (e.g., Walsh et al., 2003). A number of different empirical curves, known as photosynthesis-irradiance (PI or PE) curves, have been proposed to quantify the relationship between solar irradiance and photosynthesis (e.g., Jassby and Platt, 1976;Platt and Jassby, 1976;Platt et al., 1980). Within the CMIP6 suite, NOBM uses the simplest form, identical to its Monod nutrient-limitation terms, with a functional-group-specific light half-saturation parameter. MEDUSA, OECO, and MPI-NPZD use a hyperbolic tangent curve (Smith, 1936;Jassby and Platt, 1976): (1) where µ max (T) is the light-and nutrient-replete growth rate at temperature T, α describes the initial slope of the PI curve, and I is the solar irradiance. The remaining models base their PI curves on that of Geider et al. (1997): where α chl is the chlorophyll-specific initial slope of the PI curve, θ chl is the chlorophyll-to-carbon (chl:C) ratio, and µ max (T, N) is temperature-and nutrient-mediated growth rate under light-replete conditions. In this form, the light-limitation curve is a function of nutrient and temperature limitation, such that nutrient concentration and temperature play both a direct and indirect role in modulating growth. The WOMBAT model simplifies this relationship by using µ max (T) in place of µ max (T, N). By introducing a chl:C parameter, this function captures the process of photoadaptation, where phytoplankton increase their chlorophyll levels to compensate for low light levels. COBALT and BLING account for the effects of both irradiance and nutrient limitation when calculating the θ chl value. MARBL and CanOE compute a prognostic chlorophyll with chlorophyll biosynthesis calculated as a function of light and nutrient availability (after Geider et al., 1998). Within CMOC, PISCES, NOBM, and MEDUSA2, θ chl is a function of irradiance only. WOMBAT, OECO, and HAMOCC do not account for variation in chl:C ratios in their light-limitation factors. A final structural difference arises in how each model applies colimitation from nutrients, light, and temperature. In general, primary productivity models use one of two models: a minimum limiting factor model (Liebig and Playfair, 1843) or a multiplicative model. All models in the CMIP6 suite use a minimum model across nutrients but vary regarding the colimitation with light. PISCES, COBALT, MARBL, and HAMOCC calculate a separate minimum limiting nutrient factor and light limitation factor, both applied multiplicatively to a temperature-dependent maximum uptake rate [note that for the models that use the Geider et al. (1997) nutrient-mediated light formulation, and this implementation in effect determines whether there is enough light to reach the nutrient-limited rate]. CanOE uses a similar implementation but with the nutrientlimitation factor calculated as the minimum of N-and Febased cell quota terms rather than using the same limitation term as for uptake. CMOC, NOBM, OECO, MRI-NPZD, and WOMBAT calculate a single limiting factor across both nutrient and light limitation and apply this to a temperature-dependent uptake rate.

Temperature Influence on Biological Rates
Variations in temperature dependence at various trophic levels, both within a single model and across models, can play a key role in determining how primary and secondary production rates in each ESM respond relative to each other both seasonally and under long-term climate change scenarios. This in turn can influence whether an increase in temperature leads to higher or lower trophic transfer efficiencies (i.e., fraction of production at one trophic level relative to one below it), which is of particular interest when quantifying how changes affect upper trophic level species. For example, Reum et al. (2020) found that temperature dependency assumptions were a key source of intermodel uncertainty in their examination of the impact of climate change on the Bering Sea food web.
The exponential scaling of maximum plankton growth rate with temperature was first reported by Eppley (1972), and the resultant "Eppley-curve" is an empirical function widely used to estimate primary productivity from satellites (Morel, 1991;Behrenfeld and Falkowski, 1997) and incorporated in models (e.g., Stock et al., 2014). The majority of the CMIP6 models formulate their temperature-dependent factors using an Eppleystyle exponential curve, with coefficients that reflect updated data compilations (e.g., Brush et al., 2002;Bissinger et al., 2008).
The Arrhenius-Van't Hoff equation (Arrhenius, 1915) is an alternative to the Eppley function that has traditionally been used to describe the rates of chemical reactions. However, it can also be used to describe biological rates, as metabolic rates are typically limited by a rate-limiting biochemical step, such as Rubisco carboxylation for autotrophic growth and adenosine triphosphate (ATP) synthesis for heterotrophic growth FIGURE 9 | Schematic for the CanOE model. See Figure 1 caption for details. (Gillooly et al., 2001;Brown et al., 2004). These ideas have been incorporated into the metabolic theory of ecology (MTE; Brown et al., 2004) that broadly describes how metabolic rates scale with biomass and temperature. MTE uses the activation energy (E a ) of the rate-limiting enzymatic reaction to describe the temperature dependence of metabolic rates. This mechanistic derivation of temperature limitation leads to a similarly-shaped temperature vs. rate curve as the more empirically derived Eppley formulation; though within the temperature range of most biological processes (0-40C), there can be up to a 10-15% difference in rates when using the Eppley formulation vs. the Arrhenius (or MTE) equation (Gillooly et al., 2002). Within the CMIP6 suite, CMOC and CanOE rely on the Arrhenius-Van't Hoff equation for their temperature sensitivities.
The sensitivity of the various temperature curves used within the CMIP6 models can be compared by looking at the Q10 value, i.e., the acceleration of a given reaction when the temperature is increased by 10 • C. While the exact equations used vary, all models within the CMIP6 suite apply a similar temperature dependence to phytoplankton growth rate, with similar Q10 values (Table 5).
However, the models differ when it comes to the effect of temperature on other processes. In PISCES, COBALT, NOBM, and MARBL, zooplankton grazing rates are also temperature dependent. Within PISCES, the zooplankton grazing rate is more strongly influenced by temperature than phytoplankton nutrient uptake, while in NOBM, the reverse is true; COBALT and MARBL apply the same Q10 factors to both processes. In the remaining models, the zooplankton grazing rate is unaffected by temperature. Detrital remineralization is also temperature dependent in all models except for MRI-NPZD and HAMOCC. Within CMOC, CanOE, and PISCES, the detrital remineralization Q10 values are higher than those for phytoplankton uptake rate, while in the remaining models the Q10 factors are the same. The temperature dependence of nutrient recycling and the microbial loop have been highlighted as a key source of uncertainty in projections of both the magnitude and direction of net primary productivity under climate change (Taucher and Oschlies, 2011). Finally, phytoplankton and zooplankton mortality rates are temperature dependent in a handful of models. WOMBAT, MARBL, BLING, and COBALT include temperature terms in their phytoplankton mortality rates, while PISCES, COBALT, and MARBL apply similar terms to zooplankton mortality.
From the perspective of modeling of higher trophic levels such as fish, the assumptions regarding temperature effects are critical. While warming may increase overall primary productivity, limits on the zooplankton grazing rate and increases in zooplankton mortality may serve as bottlenecks that prevent increased primary production from reaching fish and higher trophic levels. Likewise, the relative response of remineralization processes may alter the balance between benthic vs. pelagic energy pathways. Laufkötter et al. (2017) found that changes to the temperature dependence of remineralization rates within an ESM could strongly influence surface nutrient concentration and the extent of oxygen minimum zones. Therefore, assumptions of the ESM can strongly influence predictions when coupled to higher trophic level models.

Zooplankton Functional Response
Zooplankton functional response is one area where the CMIP6 models tend to converge structurally, with all models using either a type 2 or type 3 functional response ( Table 6). A type 2 functional response, which is used in the NOBM, CanOE, MARBL, and HAMOCC models, is defined as one where grazing rate increases linearly at low prey biomass and approaches a threshold rate at high prey biomass; the functional form is concave down everywhere (Gentleman et al., 2003). A type 3 functional response is similar, except for the presence of an inflection point that decreases predator grazing rates at low prey biomass relative to a type 2 function (Gentleman et al., 2003). Sigmoidal type 3 functional responses are used in WOMBAT, CMOC, MEDUSA, OECO, and MRI-NPZD; PISCES uses a type 3 threshold model instead. COBALT falls in between, using a type 2 response but with a small prey biomass adjustment to limit grazing at very low prey concentrations. This limited number of functional response forms contrasts with the many functional responses applied in ecosystem models for higher trophic levels, such as those summarized in Hunsicker et al. (2011). Where grazers can feed on multiple prey groups (including either phytoplankton, zooplankton, or detrital prey sources), many of the models calculate grazing rates using a single-resource model, where the grazing rate on a particular prey species is determined independently of whether other prey are present. PISCES, COBALT, NOBM, MEDUSA, and CanOE do consider total prey availability when calculating grazing rates. NOBM and CanOE calculate grazing as a function of total prey, with predation distributed across the prey species proportional to prey availability. MEDUSA and PISCES both use prey preferences, such that grazing rate is a function of both prey availability and the relative preference of predator for each prey species. COBALT uses a similar preference-based model but adds a switching function such that predator preference for prey increases when that prey is abundant.

Detrital Remineralization
The biological carbon pump is a collection of processes that determine the export and associated remineralization of particulate organic carbon (POC) produced by biological activity at the ocean surface (Passow and Carlson, 2012). These processes play a key role in the global carbon cycle; because carbon is assimilated by biological processes in the surface waters and then sequestered at depth, through both the soft-tissue pump (sinking of organic matter) and carbonate pump (sinking of carbonate mineral shells), the ocean takes up far more atmospheric CO 2 (Sarmiento and Gruber, 2006, and citations within) than through the solubility pump alone. The strength of the biological pump is estimated to be between 5 and 13 Pg C yr −1 (Laws et al., 2000;Dunne et al., 2007;Henson et al., 2011), though models from the CMIP5 era ranged in the lower end of those estimates (5-8 Pg C yr −1 ) (Laufkötter et al., 2016).
In addition to the role this process plays in the global carbon cycle, detrital remineralization can affect regional patterns in nutrient concentration. The speed at which material is transported from the surface waters to below the mixed layer may affect nutrient limitation and, therefore, primary production. Particularly in shallow coastal environments, how detrital remineralization is simulated and parameterized may influence whether organic material is recycled within the water column, is fed into a benthic food web, or is subject to burial.
Aside from the models that treat organic matter remineralization in a highly generalized fashion (e.g., BLING, CMOC, MRI-NPZD), CMIP6 models generally treat detrital remineralization in similar ways, though specific implementation may vary. In general, particulate organic matter in models is produced through a combination of phytoplankton loss and aggregation processes, and zooplankton grazing and mortality (Laufkötter et al., 2016). The specific type of detrital particle impacts the particle sinking speed or remineralization rates, which are correlated quantities. Particles that sink faster and/or are protected from bacterial remineralization have lower remineralization rates and, thus, longer remineralization length scales (defined as the depth at which 63% of the organic matter has been converted to inorganic forms). A key factor determining whether a detrital particle is protected or is faster-sinking is the presence of ballast materials (lithogenic particles, biogenic silica, and calcium carbonate), originating from organisms such as diatoms, coccolithophores, and foraminifera, and dust deposition (Armstrong et al., 2002;Francois et al., 2002;Klaas and Archer, 2002;Dunne et al., 2007). ESMs focus on characterizing ballast materials properly in order to achieve a representation of both ballasting and non-ballasting POM remineralization at depth, often relying on a variant of the Klaas and Archer (2002) mineral ballasting model.
The treatment of ballasting and non-ballasting particulate organic matter in CMIP6 models can then be divided into those that use a flux attenuation scheme (MARBL, HAMOCC, and OECO) and those that have sinking particles, whether in one or multiple size classes (COBALT, CanOE, PISCES, MEDUSA, and NOBM; also see the categorization of Séférian et al., 2020). In a flux attenuation scheme, organic material losses such as non-predatory mortality, egestion, and aggregation are collected into an implicit detrital pool that is then redistributed across depth to various inorganic nutrient pools based on depth profiles of remineralization rates. Treatment of ballast materials in this case results in the "hard POM" distributing more evenly throughout deeper depths, associated with their (implicitly) faster sinking speeds (e.g., Abell et al., 2000). In the remaining schemes, the organic material losses flow to one or more detrital state variables within the model grid cell in which the loss process occurs; the detrital state variables are then subject to vertical sinking. For the size-differentiated detrital models, loss fluxes are typically categorized based on the size and type of N, NO 3 , and NH 4 indicate concentrations of total nitrogen, nitrate, and ammonium, respectively. L X refers to the limitation term associated with tracer X. K X refers to a half-saturation constant for tracer X. α is a fitting constant associated with optimal uptake kinetics. source material, with larger and ballasted particles sinking more quickly than smaller and non-ballasted ones. In the majority of models, these size-and type-specific sinking speeds are constants though some models also take the effect of temperature on sinking speed and/or remineralization into account (e.g., NOBM and COBALT).

River Runoff
For regional LMR applications, river runoff can be important due to its impact on a coastal circulation, vertical stratification, salinity, nutrient and dissolved organic matter concentration, biological production, and carbon chemistry (e.g., Dagg and Breed, 2003;Lohrenz et al., 2008;Gomez et al., 2019). Physical impacts include buoyancy-driven coastal currents such as the Alaska Coastal Current and associated mesoscale eddies (Stabeno et al., 2016). Large riverine inputs of nutrients can lead to coastal eutrophication, and the subsequent remineralization of organic matter can promote hypoxia and enhanced acidity near bottom, especially in regions with poor ventilation and weak buffer capacity (Bianchi et al., 2010;Cai et al., 2011Cai et al., , 2017Fennel et al., 2016). In addition, the biogeochemical signature of land-influenced river water is distinct from that of ocean-derived water. The usually more acidic water from rivers determines decreased aragonite saturation levels (Duarte et al., 2013), increasing the vulnerability of coastal regions to ocean acidification. Riverine input of dissolved iron can have a substantial impact on coastal production (Coyle et al., 2019).
Within the CMIP6 suite of models, a wide variety of approaches are used to resolve riverine processes ( Table 7). Freshwater fluxes in many models are calculated as part of the coupled ESM. When this is the case, the river fluxes are a function of processes within the atmosphere, ocean, and land components of the model, and they respond to changes in the simulated climate. In other models, freshwater fluxes are derived from an external data set, such as the river runoff climatologies of Dai and Trenberth (2002) and PO 4 , Fe, Si, and DOP indicate concentrations of total phosphate, iron, silica, and dissolved organic phosphorus, respectively. L X refers to the limitation term associated with tracer X. K X refers to a half-saturation constant for tracer X. α is a fitting constant associated with optimal uptake kinetics. θ represents an Fe:X ratio, where X is either N, P, or C. Dai et al. (2009). Many of the global models add runoff as a surface flux to the ocean within a specified distance from the coast (Griffies et al., 2016). This technique only partially compensates for the severely limited resolution of the true estuarine circulation and mixing (and consequent boundary currents) by coarse ocean model grids. Some of the CMIP6 models include significant improvements to the surface flux approach used in CMIP5. These include the following: (1) the estuarine box model technique (Sun et al., 2017) now incorporated into CESM2 (Danabasoglu et al., 2020) and (2) enhanced vertical mixing in coastal grid cells receiving riverine input, now incorporated into MOM6 of GFDL (R. Hallberg, pers. comm.).
Overall, the CMIP6 suite of ESMs does not represent the spatiotemporal variability in regional-scale riverine fluxes of either freshwater or of nutrients, DIC, and alkalinity well; they are typically designed to capture global-scale trends that may not be appropriate when focusing on smaller regional scales. For example, the GFDL ESM4 input forcing (coupled to COBALT) assumes a river DIC of 320 mmol m −3 and a river total alkalinity (TA) of 420 meq m −3 ; these values are defined to compensate for global DIC and TA losses such that mass balance within the biogeochemical model domain is maintained. However, observations suggest that river DIC is usually greater than river alkalinity (see Figure 8 in Moore-Maley et al., 2018). With respect to global river discharge datasets, Kearney (2019) found that the algorithms used to fill spatiotemporal gaps in river gauge measurements within the Dai and Trenberth (2002) dataset were inappropriate for snowinfluenced rivers like those emptying into the Bering Sea. Fully coupling the biogeochemical processes connecting the ocean and land systems, including those associated with rivers, the benthic ecosystem, and the continental shelf, remains an outstanding task that ESMs have not yet been implemented. When applying global ESM results to highly river-influenced areas, LMR end users may need to carefully consider potential limitations like these. Parameters include g j = maximum grazing rate (per time), g j (T) = maximum grazing rate as a function of temperature, ǫ = capture rate (biomass −1 time −1 ), P i = biomass of prey, q ij = preference for prey i, K i = half-saturation constant (biomass), ij = switching function, = Ivlev constant (biomass −1 time −1 ), P thresh = prey threshold (biomass), and F(P k ) = total prey threshold function. Note that we list the models by biogeochemical module for consistency with other tables in this study, even though river input forcing is simulation-specific rather than being an intrinsic part of the biogeochemcial model.

CONCLUSIONS
Projections of changing climate, across both short and long timescales, are playing an increasingly important role in many LMR management strategies (Link et al., 2015). The ability of management frameworks to use ESM output may be dependent on the complexity of the ecosystem models used to support decision-making that can range from population dynamics models with no environmental covariates included, through to complex end-to-end ecosystem models incorporating movement, bioenergetics, and fishing fleet dynamics (Fulton et al., 2011;Marshall et al., 2019). However, as Ecosystem-Based Fishery Management becomes more widely implemented, assessment and management processes will likely become more able to ingest ecosystem indicators and ESM projections. Even if current frameworks are not able to easily include environmental covariates, they can still be considered using Management Strategy Evaluation or scenario planning (Haward et al., 2013;Punt et al., 2016;Surma et al., 2018).
To support this endeavor, LMR scientists and managers must not only have access to climate and biogeochemical model output but also understand the differences among models in terms of simulating key processes, resolution, etc. The complexity of the suite of models and the constraints on getting needed model output have thus far been a hindrance for using the model output for LMR science and management. This study provides a summary of some of the key points LMR end users may need to consider when using the CMIP6 suite of biogeochemical model output and also provides an entry point to the more in-depth descriptions found in the primary documenting literature for each model (see Table 2).

AUTHOR CONTRIBUTIONS
KK conceived of and organized the manuscript, drafted text, and created the figures. KK, JL, and FG researched and summarized the CMIP6 model equations. All authors assisted in organizing, writing, and editing text.

ACKNOWLEDGMENTS
Thank you to our two reviewers as well as our two NOAA internal reviewers, Charles Stock and Ryan Rykaczewski, for their input to make this paper as strong as possible.