Biogeochemical versus ecological consequences of modeled ocean physics

This study compares the simulated biological differences in a “diverse ecosystem” model coupled to two different global physical model configurations, one that has a relatively high resolution and is eddy permitting and one that is much coarser and does not resolve eddies. This is an interesting starting point and the abstract of this manuscript is quite


Introduction
Ocean general circulation models have proved an invaluable tool for studying the role of phytoplankton in the global biogeochemical cycles of climatically important elements.Recent advances have resulted in ever-higher-resolution physical models of the ocean circulation (Menemenlis et al., 2008) and more complex ecological models incorporating larger numbers of phytoplankton functional groups and even individual phytoplankton phenotypes (Follows and Dutkiewicz, 2011).This trend for increasing resolution and complexity is aimed at creating model systems which incorporate some of the complexity seen in reality, with the hope of better resolving biogeochemical processes.In this study we explore how differences in physical resolution affect, differently, the biogeochemistry and ecosystem of a coupled model.
The physical framework of an ocean model is fundamental to accurately model biogeochemical cycles and phytoplankton ecology (Doney, 1999;Anderson, 2005).Observations have shown that phytoplankton biomass and community structure have characteristic temporal and spatial scales Published by Copernicus Publications on behalf of the European Geosciences Union.
corresponding most closely with the oceanic mesoscale and submesoscale (Platt, 1972;Strass, 1992;Abbott and Letelier, 1998;Doney et al., 2003;Cotti-Rausch et al., 2016).However, most global ocean models incorporating biogeochemical and ecological processes (especially those coupled in climate studies) do not resolve scales less than ∼ 1 • .In the ocean, the characteristic temporal and spatial scales of biology coincide with those of mesoscale and submesoscale physical dynamics.Spall and Richards (2000) showed in a high-resolution model of an unstable frontal jet that spatial heterogeneity in primary production occurred on scales of a few to tens of kilometers and that primary production could increase locally by up to 100 %.Coarse-resolution global biogeochemical models do not resolve these dynamics.Lévy (2008) and Mahadevan (2016) review the consequences of this for biogeochemical models.Previous studies have found that neglecting to resolve the mesoscale could result in errors of up to 30 % in the estimates of primary production (Lévy et al., 1998;Oschlies and Garçon, 1998;Mahadevan and Archer, 2000;McGillicuddy et al., 2003).Furthermore, Lévy et al. (2001) found discrepancies of up to 50 % in integrated primary production comparing a coarse-resolution model with one that resolved submesoscale dynamics.These studies found that, in the oligotrophic subtropical gyres, mesoscale (and submesoscale) dynamics drove an increased nutrient supply to the surface mixed layer, which enhanced rates of primary production.In the subpolar gyres, eddies appear to have a different effect on ocean biology.McGillicuddy et al. (2003) found in a 0.1 • -resolution model of the North Atlantic that mesoscale processes drove a geostrophic adjustment to deep winter convection, which reduced nutrient supply.However, nutrients are less likely to be limiting than light in the subpolar gyres, so this may also have a positive effect on rates of primary production.Although these previous studies predict the role of the mesoscale in modulating biological and biogeochemical responses in different regional settings, it is unclear what the integrated effect will be globally and what downstream effects might result.
The above studies, although they resolved higherresolution physics, typically employed rather simple biogeochemical models incorporating only one or two phytoplankton functional types (PFTs).However, marine microbial communities are known to be incredibly diverse, and this diversity plays an important role in mediating global biogeochemical cycles.Thanks to the continuing expansion of computing resources, diversity has been included in global biogeochemical models which resolve several phytoplankton functional groups (Chai et al., 2002;Gregg et al., 2003;Le Quéré et al., 2005), and even several tens of phytoplankton phenotypes within multiple functional groups have been developed (Follows et al., 2007;Ward et al., 2012).It is unknown whether a change in physical resolution will result in any changes in the emergent community structure of one of these diverse models.Sinha et al. (2010) found that an intermediate-complexity ecosystem model which re-solved five PTFs, run with two different physical models at a similar (coarse) resolution, could result in regional changes in the modeled phytoplankton communities.However, it is unclear whether the emergent community resulting from a more complex ecosystem model will be more or less robust to changes in the physical forcing.
Here we present a process study that explores the effect of refining the physical resolution on a global, diverse ecosystem model which incorporates 78 distinct phytoplankton phenotypes.We explore both the effect on the bulk biogeochemical properties and the community structure of the model solutions.The objective of this study is not to assess which model performs best with respect to reality or to suggest an optimal resolution.Rather, we aim to examine how changes in the resolution and parameterization of subgrid scale processes of the model domain alter the emergent biogeochemical and ecological properties of this diverse ecosystem model and specifically to identify tools to help understand these differences and similarities.

Method
This study is based upon numerical simulations of global ocean circulation, biogeochemical cycles and diverse phytoplankton populations.We have employed the Massachusetts Institute of Technology General Circulation Model (MITgcm; Marshall et al., 1997) and biogeochemical and ecological model components as detailed in Dutkiewicz et al. (2009).Below, we describe these physical circulation and biogeochemical-ecological models in more detail.

Physical model configurations
We used two physical model configurations developed by the Estimating the Circulation and Climate of the Oceans (ECCO) project (Wunsch et al., 2009).Both span the period 1992-1999 and were constrained to be consistent with observed altimetry and hydrography.The high-resolution ECCO2 configuration, with an effective resolution of 1/6 • and horizontal grid dimensions of ∼ 18 km (ECCO2; Menemenlis et al., 2008), is referred to as eddy permitting, as it resolves large eddies at low latitudes, but remains below the first baroclinic Rossby radius at high latitudes (Chelton et al., 1998).We compare this high-resolution configuration with the ECCO-Global Ocean Data Assimilation Experiment state estimate (ECCO-GODAE, also referred to as ECCO version 3 in the lineage of ECCO production solutions; Wunsch andHeimbach, 2007, 2013), which has a coarser grid resolution of 1 • .Both physical models are data assimilation products.The ECCO-GODAE product is based on the Lagrange multiplier method, and the ECCO2 product employs a simplified Green's function method.It should be noted that unlike most so-called ocean "reanalysis" products, both of the ECCO products are obtained using "smoother" methods, which avoid the artificial adjustment motions that can be triggered during what is called "analysis increments" in data assimilative models which use "filtering" methods (Wunsch and Heimbach, 2013;Stammer et al., 2016).Because the ECCO products are free of such artificial adjustments (expressed in particular in strong vertical adjustment motions that may confound vertical transport of biogeochemical tracers), they also conserve tracer and momentum budgets exactly.
A core focus of this study is comparing and contrasting the modeled ecological and biogeochemical behavior in a mesoscale eddy-permitting model (ECCO2) with a noneddying model (ECCO-GODAE).Accordingly we limit our analysis latitudinally to the region between 60 • S and 60 • N. Within this latitudinal band, the first baroclinic radius of deformation (Chelton et al., 1998) is larger than the ECCO2 model horizontal resolution, so the ECCO2 model admits corresponding mesoscale eddy dynamics.North and south of this latitude line, mesoscale eddies are not well resolved in either model, so neither model is in a so-called large-eddy regime (Smagorinsky, 1963).In the excluded high-latitude regions the two physical model mixed layers differ systematically.These differences result from the absence of parameterized mesoscale eddy dynamics in the ECCO2 model (Danabasoglu et al., 1994).This contrast in mixed layer depth, rather than behaviors due to mesoscale eddies, sets the differences in the biogeochemical response seen in ice-free high latitudes.For simplicity, we will refer to the ECCO-GODAE simulation as CR and the ECCO2 simulation as HR.

Ecological and biogeochemical model
The ecological model used in this study has previously been discussed in Follows et al. (2007) and Dutkiewicz et al. (2009).Briefly, we transport inorganic and organic forms of nitrogen, phosphorous, iron and silica, and resolve 78 phytoplankton phenotypes and two simple grazers.The biogeochemical and biological tracers interact through the formation, transformation and remineralization of organic matter.Excretion and mortality transfer living organic material into sinking particulate and dissolved organic detritus which are transpired back to inorganic form.The time rate of change in the biomass of each of the modeled phytoplankton types, P j , is described in terms of a light-, temperature-and nutrientdependent growth, sinking, grazing, mortality and transport by the fluid flow.Many realizations of this ecological model, coupled to the ECCO-GODAE physical circulation, have been used to study a range of ecological questions, e.g., the role of top-down controls in setting patterns of phytoplankton diversity (Prowe et al., 2012;Vallina et al., 2014), the biogeography of nitrogen-fixing phytoplankton (Monteiro et al., 2011) and role of transport in setting patterns in phytoplankton diversity (Clayton et al., 2013).
In this study, the ecological model was initialized with 78 phytoplankton phenotypes with a broad range of phys-iological attributes.Phytoplankton were assigned to one of two broad size classes by random draw at the initialization of the model, and a set of physiological tradeoffs that reflect empirical observations were imposed accordingly.We stochastically assigned plausible values for the nutrient half-saturation constants (κ N ), light and temperature sensitivities within each phytoplankton functional group (Prochlorococcus-like, picophytoplankton, diatoms and large phytoplankton).Both the HR and CR simulations were initialized with an identical set of phytoplankton phenotypes and initial conditions for all phytoplankton phenotypes, and in both cases the ecosystem model was forced with identical initial conditions for all variables, light forcing and dust inputs.Initial conditions for nutrients were taken from the January climatological values given in the World Ocean Atlas 2005 (Garcia et al., 2006).Interactions with the environment, competition with other phytoplankton, and grazing determine the composition of the phytoplankton communities that persist in the model solutions.
We run both models for a total of 8 years, for the period from 1992 to 1999 from the same initial conditions.Although this is a relatively short run for a biogeochemical model, it currently would not be feasible to spin up the HR simulation for much longer.In several previous studies using the CR model (Dutkiewicz et al., 2009(Dutkiewicz et al., , 2012) ) and in this study for both CR and HR, we find that it takes only 3 years for the ecosystem to reach a stable annual cycle, so a run of 7-8 years is sufficient to produce a repeating seasonal cycle in the modeled phytoplankton community.There was no discernible trend in the globally integrated annual average phytoplankton biomass in the last few years of either of the simulations.The regional nutrient drifts that are linked to slow changes in the deep nutrient concentrations are small and significantly less than the seasonal and interannual variability.
We compare the model results for the last year (1999) from both model configurations to test the sensitivity of the ecosystem to the modeled ocean physics.The results of that comparison are described in the following section.

Results
We describe differences in some of the physical properties most directly relevant to biogeochemical processes: sea surface temperature (SST) and mixed layer depth (MLD).Differences in other physical fields between the two model configurations have previously been discussed in Clayton et al. ( 2013); here we examine only those physical processes most directly relevant to biogeochemical processes not considered in that previous work.The SST patterns in both models are broadly similar, but we did find local differences in SST of up to 3 • C in some regions (Fig. 1).The HR simulation appeared to have a slight cool bias relative to the CR simulation, which is an indicator of enhanced upwelling and/or vertical mixing in the higher-resolution configuration.There  were marked regional patterns in the MLD differences between model configurations.Annual average MLDs were consistently deeper in the low latitudes in the HR simulation, whereas they tended to be slightly shallower in the midlatitudes.As mentioned above, we exclude the high latitudes (> 60 • ) from our analysis.

Primary production and biomass
Both model configurations result in largely similar patterns in phytoplankton biomass and primary production, with low biomass and productivity associated with the subtropical gyres, and with higher biomass and productivity found in the mid-latitudes and upwelling zones.Although globally integrated primary production is very similar between both models, with only a 0.6 % overall increase in primary production in the HR simulation, there are clear regional differences (Fig. 2), with higher rates of production in the low latitudes and lower rates of production in the mid-latitudes, in the HR simulation.We find a ∼ 20 % global increase in the standing stock of phytoplankton biomass in the HR simulation, mostly accounted for by higher phytoplankton biomass between 40 • S and 40 • N. The largest differences in both phytoplankton biomass and primary production are associated with the western boundary currents and the equatorial upwelling zone, with higher values in the HR simulation.We also see decreases in phytoplankton biomass and productivity in the HR simulation associated with the boundaries between different biogeographical provinces.

Phytoplankton community structure
The ecological model represents a diverse community of phytoplankton, made up of individual phenotypes (analogous to species or ecotypes), which are defined by their nominal size, their nutrient requirements, their optimal light and temperature ranges, and their palatability to grazers.Each phytoplankton phenotype belongs to one of four possible modeled phytoplankton functional groups.Phenotypes within each functional group have similar cell sizes and nutrient requirements but differ in their optimal light and temperature ranges.To understand the impact of model resolution on the modeled phytoplankton community structure, we examine the regional patterns in the dominant modeled phytoplankton functional groups (Fig. 3a and b) and dominant modeled phytoplankton phenotypes (Fig. 3c and d).We find that the emergent phytoplankton community structure is remarkably similar between the two simulations.The low latitudes are dominated by Prochlorococcus analogues, whereas the mid-latitudes are dominated by diatoms and large phytoplankton (Fig. 3a and b), with very little difference between models.What is even more striking, is that the regional patterns in the dominant phytoplankton phenotypes are largely unchanged between simulations.The only region that shows Other regions where the models do not agree occur mainly at the borders between different biogeochemical provinces.This suggests that changes in model resolution do not drive changes in the geographical range of dominant phenotypes, but this also suggests that they do not result in ecological shifts where different modeled phenotypes might become dominant.Given the big differences in modeled primary production and biomass described above, this is an unexpected result.
Differences in the patterns of phytoplankton biodiversity between the model simulations have previously been described in Clayton et al. (2013).Although we find an overall increase in phytoplankton biodiversity in the HR simulation (not shown here), those differences are mainly driven by the persistence of more rare species in the HR simulation with respect to CR and are thus unlikely to have a significant effect on integrated bulk biogeochemical processes.

Nutrients and nutrient limitation
Differences in model resolution are expected to drive changes in the supply of nutrients to the surface photic zone.We do find differences in the concentration of surface macroand micronutrients (Fig. 4).Nitrate and dissolved iron both show marked patterns in the difference of their surface distributions between simulations.The concentration of nitrate remains unchanged in the Atlantic, the Indian Ocean, and in some parts of the North and South Pacific (Fig. 4a).However, there is a decrease in surface nitrate concentrations in the HR simulation in the mid-latitudes of the North Pacific and the Southern Ocean as well as in the subtropical and equatorial Pacific.Conversely, we see an increase in surface nitrate in the HR simulation in the Brazil-Malvinas Confluence Zone and the region of the Subantarctic Front in the Southern Ocean.Differences in surface dissolved iron concentrations exhibit markedly different patterns (Fig. 4b).The surface concentration of dissolved iron remains unchanged between model simulations in the subtropical and equatorial Pacific and parts of the subtropical North Atlantic.We see a general decrease in surface iron in the HR simulation in the Subantarctic Front in the Southern Ocean and the southern part of the Indian Ocean and an increase in dissolved iron in the mid-latitudes of the North Pacific, the equatorial Atlantic, and the subpolar North Atlantic.
In addition to assessing the differences in surface nutrient concentrations, we determined the locally limiting nutrient for both models (Fig. 5).Nitrate and dissolved iron are the spatially dominant limiting nutrients in both model simulations, with almost identical regional patterns of limitation.However, we do find small but significant regions where the dominant limiting nutrients in both simulations do not correspond (Fig. 5b), primarily located at the boundaries between biogeographical provinces.

Discussion
In this study, we have assessed the effect of explicitly representing mesoscale dynamics in a global ocean ecologicalbiogeochemical model.Strikingly, we find that the realized phytoplankton communities in both simulations are remarkably similar but that there are marked regional variations in bulk ecosystem properties such as primary production and phytoplankton biomass.We also find that although the general regional patterns of nutrient limitation remain unchanged, surface concentrations of nitrate and dissolved iron can vary markedly between simulations in some regions but are almost completely unchanged in others.Why do we see such marked changes in the distribution of biogeochemical properties of the models but not in the modeled phytoplankton community?Here we explore what drives these similarities and differences, frequently using simple theoretical constructs.The goal here is to provide a framework that other modellers can use to help understand some of the implications of the physical-biogeochemical and ecological consequences of different forcings and resolutions.
Ultimately, any differences in the biogeochemical and ecological properties between the two model solutions occur either because there are differences in the local physical fields (e.g., MLDs and transport) or because physical features of the ocean circulation such as the western boundary currents and gyre boundaries are realized in different locations.We identify three main regions that are affected in different ways by the addition of mesoscale dynamics: the tropical and subtropical regions, the subpolar gyres and the boundaries between biogeographical provinces.We explore the linkages between the physical and biological components of the system to put these differences into context.

Tropical and subtropical regions
In the low latitudes, we found that the addition of mesoscale dynamics resulted in an overall deepening of the annual mean mixed layer depth (Fig. 1).We found an increase in phyto- plankton biomass and primary production, but the dominant phytoplankton functional group and phenotypes remained unchanged.The increase in primary production and phytoplankton biomass in the HR simulation was not entirely surprising, as it has been predicted by previous regional and idealized studies on the effect of increasing model resolution on primary production (Lévy et al., 2001;Oschlies and Garçon, 1998;Mahadevan and Archer, 2000;Spall and Richards, 2000;McGillicuddy et al., 2003).However, we may have expected a large increase in the biomass of large phytoplankton thanks to more episodic eddy-driven nutrient injections to the surface mixed layer, as observed by Benitez-Nelson et al. ( 2007) and Brown et al. (2008).In fact, the increase in phytoplankton biomass was driven by an increase in the abundance of the dominant phytoplankton functional groups, Prochlorococcus analogues and picophytoplankton, which are both small, gleaner types.We also found that where a particular nutrient was limiting in both model simulations, its surface concentration remained unchanged.
The subtropical gyres are steady, stable regions where seasonality is low and phytoplankton growth is nutrient limited.In this context, we can apply the simple resource competition framework introduced by Tilman et al. (1982) and applied to this system by Dutkiewicz et al. (2009) to understand the simultaneous differences and similarities in the model results.We set up a simple model system, where phytoplankton growth is nutrient limited and balanced by a simple linear mortality term: where R is the limiting resource, P is the phytoplankton biomass, µ is the maximum phytoplankton growth rate, k is the resource half-saturation constant, m is the phytoplankton mortality rate and S R is the rate of resource supply into the system.The steady state solution for this system is In this framework, the steady state concentration of the limiting resource, R * , is controlled not by the magnitude of the nutrient supply but by the physiological attributes of the dominant phytoplankton phenotype.In this type of stable system, the emergent dominant phytoplankton phenotype will be the one that can draw the resource down to the lowest R * value.As we have set the phytoplankton communities to be composed of the same phenotypes in both simulations, and R * is set by the physiological traits of the phytoplankton (Dutkiewicz et al., 2009), we select for the same dominant phenotype in both simulations.Conversely, the steady state concentration of the phytoplankton, P * , is a function of the resource supply, S R , as well as the mortality rate.Thus, an increase in S R results in increased phytoplankton biomass, but not in a shift in the dominant phytoplankton type, as R * is only a function of the physiology of the fittest phytoplankton type.Enhanced mixing in the HR simulation drives an increase in S R but does not change the dominant phytoplankton phenotype.Rather, it results in an increase in the biomass of the dominant type and an increase in productivity (µP * ).At the same time, the concentration of the limiting nutrient, R * , remains unchanged (as seen in nitrate and dissolved iron, Fig. 4), but the non-limiting nutrients decrease due to increased removal by the increased productivity (Fig. 4).This may also result in a decrease in the supply of resources downstream from the subtropical gyres.
For instance, we see a dramatic increase in phytoplankton biomass and primary production in the equatorial region of the HR model simulation.We evaluated the annual average vertical advective flux of nitrate (wNO 3 ) at 100 m between 5 • N and 5 • S for both models.The regionally integrated mean annual vertical nitrate fluxes evaluated at 100 m for model year 1999 were 453.1 and 383.7 mmol NO 3 m −2 yr −1 for the HR and CR simulations, respectively.There is a clear increase in vertical nutrient supply in the equatorial upwelling zone in HR that can account for the dramatic increase in phytoplankton stock and primary productivity in the equatorial Pacific.

Subpolar regions
We see a somewhat different picture in the mid-to high latitudes, where there is an overall decrease in annual mean primary production, and to a lesser degree in phytoplankton biomass in HR compared to CR (Fig. 2).Again, the dominant phytoplankton functional groups and phenotypes remain largely unchanged between models in these regions (Fig. 3).What is most striking is the much higher abundance of zooplankton found in the northern subpolar gyres of the CR model solution (Fig. 6).In order to explain these differences, we must look more closely at the evolution of the seasonal cycle in both models (Fig. 6).
First we consider the different behaviors of the models during the winter months.HR has deeper mixed layer depths from November through March or April, depending on latitude.During this period, higher light limitation explains the lower primary production levels.In the context of this model, the limiting factors on the phytoplankton growth rate are multiplicative -e.g., for each modeled phytoplankton phenotype j , the growth term is given by where γ T and γ I are the temperature and light limitation terms, respectively.During the winter, when nutrients are replete (and the nutrient limitation term approaches 1), if either the light or temperature fields experienced by the modeled phytoplankton are consistently less favorable, then µ will be decreased.Consistently deeper winter MLDs in the HR simulation during the winter months result in lower winter primary production in HR than CR.
The onset of the spring bloom occurs roughly 1 month earlier (March-April) in the HR simulation than in the CR simulation (April-May).This is reflected in the higher primary production in the HR simulation in March and April, followed by a reversal with higher primary production in the CR simulation in May.This is driven by differences in the timing of the shoaling of the MLDs, which occurs earlier in the HR model due to the action of mesoscale eddies.
During the summer and into the early autumn, the shallower mixed layer depths restrict nutrient resupply and nutrients become limiting.As the physical and biogeochemical properties of both models are relatively stable over several months during the summer period, we represent them as an idealized steady state system (see e.g., Dutkiewicz et al., 2009, Fig. 4) to better understand the processes driving the differences between the two models.In these higher latitudes, grazer dynamics become uncoupled from their prey.Thus we construct an idealized model system similar to Eqs. ( 1) and ( 2) but that now includes explicit grazing by zooplankton.
where R is the resource, P is the phytoplankton biomass, µ MAX is the maximum phytoplankton growth rate, k is the nutrient half-saturation constant, g is the grazing rate of phytoplankton by zooplankton, m z is the zooplankton mortality S. Clayton et al.: Biogeochemical consequences of ocean model resolution rate and S R is the rate of resource supply into the system.The steady state solution for this system is For this system, P * , the phytoplankton standing stock, is controlled by the physiological traits of the zooplankton which are constant, so regardless of changes in S R phytoplankton biomass will remain the same.However, Z * , the zooplankton standing stock, is directly proportional to S R , so increased nutrient supply results in higher zooplankton abundance.Applying this framework to our model results, we are able to explain the increase in zooplankton biomass, the roughly unchanged phytoplankton biomass and the increased primary production in CR compared to HR.In contrast to the situation during the winter, in this case it is driven primarily by differences in nutrient supply.We evaluated the annual mean vertical nitrate supply by advection, which was 142.1 and 108.1 mmol NO 3 m −2 yr −1 in CR and HR, respectively.That, coupled with the deeper summer MLDs in CR, results in an increased nutrient supply compared to HR. Zhang et al. (2013) find a similar response to nutrient injections at a shelf break front, where increases in primary production are funneled up the trophic levels and result in increased zooplankton biomass but are not reflected in increased phytoplankton biomass.On a larger scale, Ward et al. (2012) found in a global, size-structured ecosystem model that, in regions of higher nutrient supply, top-down control drives an increase in the biomass of larger plankton size classes relative to oligotrophic regions rather than an unchecked increase in smaller size classes.Similar patterns can be seen in the Southern Hemisphere, where the MLD shoals earlier in the spring in the HR simulation, and primary production and phytoplankton biomass are both lower in the HR simulation during the summer.

Boundaries between biogeographical provinces
We found marked differences in biological fields which appear to be due to differences in the geographical extent of biogeographical provinces between simulations.Differences in modeled phytoplankton biomass and primary production in the North and South Pacific subtropical gyres coincide exactly with geographical differences in regions defined by different limiting nutrients in the simulations (Fig. 2).We attribute this decrease in production and biomass to a decrease in the lateral supply of nutrients from the equatorial upwelling region.This decreased lateral supply out of the equatorial region in the HR simulation is due to more local primary production in the equatorial region (Fig. 2).Although dissolved iron is drawn down to the same concentra-tion in both simulations, primary production and phytoplankton biomass are higher, and nitrate is much lower in the HR simulation.As discussed above, resource competition theory predicts that the concentration of the limiting nutrient, R * (in this case dissolved iron; Fig. 5), is unchanged when the phytoplankton community remains unchanged.Increased dissolved iron supply associated with deeper MLDs would drive the increase in primary production and phytoplankton biomass observed in the in the HR simulation.However, this increase in local primary production, consuming higher levels of non-limiting macronutrients locally, reduces the Ekman transfer of non-limiting nutrients to the neighboring subtropical gyres (Dutkiewicz et al., 2005).A similar shift in biological transition zones, associated with model resolution, was found in a regional study of the California Current system.Fiechter et al. (2014) found that increasing the resolution of their model from 1/3 to 1/30 • resulted in a marked shift in the location of the transition between nearshore outgassing and offshore absorption of CO 2 .These differences in gas fluxes were driven, in their study, by different patterns in nutrient upwelling and transport offshore.

Stability of the phytoplankton community structure
As we have discussed in the previous sections, one of the most striking results of our model comparison is that despite the difference in the modeled physics, the emergent phytoplankton communities in both simulations are almost identical.We only know of one comparable study where the response of a complex ecosystem model to different model physics has been evaluated (Sinha et al., 2010).In that study, although there was a difference in resolution between the two physical models used to force the ecosystem model (1 • vs. 2 • ), unlike our study, both of these physical models were too coarse to resolve eddies.Although the physical models in the (Sinha et al., 2010) study were initialized and forced in essentially the same way, they were run out without being constrained to observations.This resulted in large differences in their model solutions for physical fields (e.g., seasonal SST, see their Fig.3), which were used to force the ecosystem model.One very key difference in our study with respect to the earlier work of Sinha et al. (2010) is that both of the physical models used to force our ecosystem model, ECCO-GODAE and ECCO2, were state estimates constrained to converge as closely as possible to observational data.Although the state estimates were constrained in different ways, and differences between them remain, compared to freerunning models without any assimilation, these differences are comparatively small.In the end, the differences in the ECCO state estimates used in this study are largely a consequence of whether or not they represent mesoscale dynamics.

Conclusions
In this study we have compared and contrasted the behavior of a complex ecological-biogeochemical model when coupled to either a mesoscale eddy-permitting high-resolution physical model (HR, ECCO2) or a non-eddying coarseresolution physical model (CR, ECCO-GODAE).We found that increasing the model resolution to include mesoscale dynamics does not greatly affect the structure of our modeled phytoplankton ecosystem.However it does have a significant effect on the regional distribution of the bulk properties of the ecosystem: primary production and phytoplankton and zooplankton biomass.We have used theoretical constructs to help us pull apart the reasons for these differences and similarities.
One of the most striking results of this study is the robustness of the emergent modeled phytoplankton community.We found that the dominant phytoplankton functional groups and phenotypes remained unchanged between simulations, despite differences in SST and MLDs.In contrast, we found marked regional differences in phytoplankton and zooplankton biomass and in primary production.By applying concepts from resource competition theory (Tilman et al., 1982;Dutkiewicz et al., 2009), we can explain why in the subtropical gyres, despite increased primary production, the dominant phytoplankton phenotype and the surface concentration of the limiting nutrient remain the same in both simulations.The combination of low seasonality and low grazing pressure means that despite an increased nutrient supply, phytoplankton with the lowest R * will always be selected for in this region.Conversely, deeper winter and shallower summer MLDs in the HR model, along with differences in the timing of the spring shoaling of the MLD, result in lower primary production and lower zooplankton biomass in HR than in the CR model.
Given the complexity of our ecosystem model, which incorporates 78 individual phytoplankton types, it may seem surprising that our modeled phytoplankton community structure is so similar in both cases.This presents interesting implications for marine biogeochemical and ecological modeling.It is clear that accuracy in the representation of the physical dynamics of the environment is necessary for effectively modeling ocean biogeochemistry.The higher taxonomic resolution of our ecological model may in fact allow for subtler gradations of change in the phytoplankton community when the environment is changed, whereas in coarser ecological models regime shifts could easily result due to larger differences between modeled phenotypes.We do not advocate simply tuning parameters to get the "right" result, but rather increasing the physiological parameter space constrained by laboratory and observational work in order to create a more robust and representative model of the phytoplankton community.
We have shown that the bulk biogeochemical properties of this ecological model are more sensitive to differences in modeled ocean physics than the structure of the ecosystem itself.Given that this model has many similarities to other widely used biogeochemical models, which also resolve multiple phytoplankton PFTs, this study provides important insights into how these models might behave under different physical conditions.Our results show that the physics resolved by the models do matter, regardless of the scope of the question being addressed.The modeled MLD plays a central role in mediating bulk biogeochemical processes, specifically through vertical nutrient supply and the modulation of the light environment, which ultimately control the magnitude of bulk biological rate processes.Notably, although global values of primary production, and to a lesser extent phytoplankton biomass, are similar between models, the bulk biogeochemical properties of the model solutions differ regionally.This shows that regardless of whether one is interested in global or regional questions, model resolution is crucially important.
It may be tempting to conclude from the results presented here that, in fact, model resolution is less important when considering ecological problems.However, we would strongly caution against this as although we have shown that the dominant constituents of the phytoplankton community remain largely unchanged regardless of model resolution, significant differences in phytoplankton productivity, diversity (Clayton et al., 2013) and overall community composition do result from differences in model resolution.
We have shown in this study that, even for complex ecological models, it is possible to explain differences in model solutions, and a similar approach could be taken to evaluate the effect of different model physics on different ecological model formulations.
Data availability.Due to the large volume of the model output data used in this study, we are unable to make it publicly available.However, the corresponding author is happy to provide selected model fields upon request.The modeling framework used in this study, MITgcm, can be downloaded from http://mitgcm.org.

Figure 1 .
Figure 1.(a) Annual average sea surface temperature (SST) and (c) annual average mixed layer depths (MLD) in the HR simulation for the 1999 model year; (b) the difference in SST and (d) in MLD between the two simulations.Positive values indicate deeper MLD in the HR simulation, and negative values indicate deeper MLD in the CR simulation

Figure 2 .
Figure 2. (a) Annual average total phytoplankton biomass in g C m −2 and (c) annual primary production in g C m −2 yr −1 in the HR model solution for 1999.(b) The difference in phytoplankton biomass and (d) the difference in annual primary production.Positive values indicate higher values in the HR simulation, and negative values indicate higher values in the CR simulation.The solid black contour lines in (b, d) indicate the region where the limiting nutrient differs between the two model simulations

Figure 3 .
Figure 3. Dominant phytoplankton functional group (a) in the HR simulation and (b) in the CR simulation.Diatoms are shown in red, large phytoplankton in yellow, picophytoplankton in green and Prochlorococcus-like phenotypes in blue.Dominant phytoplankton phenotype (c) in the HR simulation and (d) in the CR simulation.Each color represents a different phenotype.

Figure 4 .
Figure 4. Annual average surface concentrations of (a) nitrate in mmol N m −3 and (c) dissolved iron in mmol Fe m −3 in the HR model solution for model year 1999; (b) the difference in nitrate and (d) the difference in dissolved iron.Positive values indicate higher values in the HR simulation, and negative values indicate higher values in the CR simulation.

Figure 5 .
Figure 5. Limiting nutrients were assessed based on a biomass weighted average of the most limiting nutrient for each of the 78 total phytoplankton types (a) in the HR model and (c) in the CR model.Iron-limited regions are shown in orange, nitrate-limited in light blue, phosphate-limited in dark blue and silica-limited in red.(b) Regions where the model simulations predict different limiting nutrients are shown in dark blue, whereas pink regions have the same limiting nutrient in both simulations.

Figure 6 .
Figure 6.Hovmöller diagram showing the seasonal evolution of zonally averaged model properties for the HR simulation (left column), the CR simulation (middle column) and their difference (right column).First row: mixed layer depth; second row: phytoplankton biomass; third row: primary production; fourth row: zooplankton biomass; fifth row: surface nitrate concentration.The difference is calculated as HR − CR, so positive values indicate larger values in the HR simulation.