The assessment of a global marine ecosystem model on the basis of emergent properties and ecosystem function : a case study with ERSEM

Ecosystem models are often assessed using quantitative metrics of absolute ecosystem state, but these model– data comparisons are disproportionately vulnerable to discrepancies in the location of important circulation features. An alternative method is to demonstrate the models capacity to represent ecosystem function; the emergence of a coherent natural relationship in a simulation indicates that the model may have an appropriate representation of the ecosystem functions that lead to the emergent relationship. Furthermore, as emergent properties are large-scale properties of the system, model validation with emergent properties is possible even when there is very little or no appropriate data for the region under study, or when the hydrodynamic component of the model differs significantly from that observed in nature at the same location and time. A selection of published meta-analyses are used to establish the validity of a complex marine ecosystem model and to demonstrate the power of validation with emergent properties. These relationships include the phytoplankton community structure, the ratio of carbon to chlorophyll in phytoplankton and particulate organic matter, the ratio of particulate organic carbon to particulate organic nitrogen and the stoichiometric balance of the ecosystem. These metrics can also inform aspects of the marine ecosystem model not available from traditional quantitative and qualitative methods. For instance, these emergent properties can be used to validate the design decisions of the model, such as the range of phytoplankton functional types and their behaviour, the stoichiometric flexibility with regards to each nutrient, and the choice of fixed or variable carbon to nitrogen ratios.


Introduction
Numerical models of the environment are used frequently for informing policy decisions, for forecasting the impact of climate change, and to obtain a deeper understanding of nature.
In order for a model to be used for any of these purposes, the model must first be shown to be a valid representation of the system under study.There are two objective strategies available to demonstrate that the model is a valid representation of the system under study.The first strategy is to reproduce the spatial and temporal distributions of historic observations, and the second is to reproduce the functional relationships.
There is a long history of works that demonstrate model validation using static fields, spatial distributions and dynamic variability, including Droop (1973), Fasham et al. (1990), Taylor (2001), Blackford (2004), Allen et al. (2007), Jolliff et al. (2009), Shutler et al. (2011), Saux Picart et al. (2012), de Mora et al. (2013), and Kwiatkowski et al. (2014).However, validating a modern biogeochemical model using static fields and spatial distributions may give an appropriate assessment of the coupled biogeochemical and hydrodynamic modelled system, but the performance of the biogeochemical model may be obscured by deficiencies in the modelled circulation.For instance, the point-to-point analysis described in de Mora et al. (2013) is vulnerable to discrepancies between the model and the observations in the location of important circulation features such as fronts, coastlines or upwelling regions.These problems in the physical model may needlessly penalise the performance of the biogeochemical model when validating using point-to-point matching.Validation methods that use historic data are also sensitive to Published by Copernicus Publications on behalf of the European Geosciences Union.

L. de Mora et al.:
The assessment of a global marine ecosystem model using emergent properties initial conditions and the boundary conditions of the model.These problems are amplified for models with coarse spatial and temporal resolution, such as the monthly mean of a 1 • global model.The disentanglement of the performance of the biogeochemical model from that of the physics is a major challenge in marine ecosystem modelling (Holt et al., 2014).
Direct comparisons of model to data inform only about how similar the model is to the observations.Such methods also risk compartmentalising the validation of ecosystems and may not cover the interaction of their parts.The ability of a model to represent present-day measurements is important, but it does not inform about the models representation of the behaviour of the ecosystem as a whole.Furthermore, historical static fields may not necessarily validate a model, which is subjected to a changing climate due to the scarce availability of long-term time-series data sets.
As a solution to the problem of the absence of data and presence of poorly constrained physics, Holt et al. (2014) wrote "there is a need for metrics that assess the fidelity of the biogeochemical processes independently of the physics".In that work, they identify one such functional relationship: the link between diatom chlorophyll and total community chlorophyll.They demonstrated that the fraction of the community chlorophyll that originates in diatoms increases with total chlorophyll in multiple models.A relationship between diatom concentration and total chlorophyll was also observed in nature in Hirata et al. (2011).In effect, a relationship seen in in situ observations also appeared in multiple biogeochemical models.In addition, the relationship observed in Holt et al. (2014) was a widespread general response in all the plankton models that was independent of local hydrodynamic conditions.Furthermore, this relationship is a largescale property of the marine ecosystem, and is expected to hold true even in regions with few historical measurements.This relationship is important because it occurred independently of the hydrodynamic model, and because it reflected the functioning of the modelled ecosystem in a way that would not be visible in a simple point-to-point comparison of ecosystem state.
Beyond Hirata et al. (2011), there are many works that link a large phytoplankton size class with the community structure: Devred et al. (2011), Brewin et al. (2012Brewin et al. ( , 2014Brewin et al. ( , 2015)).Many features of an ecosystem can affect the balance of large phytoplankton chlorophyll against the total community chlorophyll, such as the large phytoplankton response to nutrients, light and temperature, competition for light and nutrients from other phytoplankton, and predation on large phytoplankton relative to other phytoplankton classes.Each of these interactions between two or more components of the ecosystem are examples of ecosystem functions.In the context of marine ecosystem modelling, ecosystem functions are physical, biological or geochemical interactions, processes or relationships that take place within a model's framework.Ecosystem function in models can be both explicitly enforced during the model development and tuning, or they can occur without being explicitly parameterised.The interplay of multiple ecosystem functions can result in the emergence of observable relationships.The link between diatom chlorophyll and total community chlorophyll, as shown by Holt et al. (2014), is an example of such an emergent relationship.These emergent relationships can be used to characterise and validate the ecosystem and its functioning.As in the example from Holt et al. (2014), emergent relationships are important because they occur independently of the hydrodynamic model, and because they reflect the functioning of the modelled ecosystem in a way that would not be visible in a simple point-to-point comparison of ecosystem state.
A selection of historically published large-scale emergent relationships are proposed to illustrate the validation of a complex ecosystem model and demonstrate the power of emergent property validation.The example ecosystem model used here is European Regional Seas Ecosystem Model (ERSEM), and it is run coupled with the Nucleus for European Modelling of the Ocean (NEMO) in a global hindcast scenario.The emergent relationships shown here are the community structure, the carbon to chlorophyll ratio, the ratio of particulate organic carbon against particulate organic nitrogen and stoichiometric balance.
After this introductory section, Sect. 2 contains a brief description of the circulation model, NEMO, the ecosystem model, ERSEM, and the sea ice model, CICE, used in this study.Section 3 is a non-exhaustive list of examples of ecosystem relationships that have been investigated in the ERSEM ecosystem model.An expanded version of the community structure relationship described by Holt et al. (2014) is included in Sect.3.1.Section 3.2 shows the ratio of carbon to chlorophyll in phytoplankton and particulate organic matter.Section 3.3 demonstrates how the model reproduces the ratio of particulate organic carbon and nitrogen as described by Redfield (1934), Martiny et al. (2013).Section 3.4 illustrates the internal stoichiometric relationships for ERSEM for each of the nutrient currencies modelled.Finally, Sect. 4 discusses the successes, potential and limitations of these methods.

The ERSEM and NEMO models
ERSEM is a lower trophic level biogeochemical cycling, carbon-based process model that uses the functional-group approach (Baretta et al., 1995;Blackford, 2004;Butenschön et al., 2015).The carbon, nitrogen, phosphorus, silicon and iron cycles are explicitly resolved.The pelagic ERSEM model simulates four phytoplankton functional types, three zooplankton functional types, one bacterial functional type and six detritus size classes.It contains variable stoichiometric ratios for each of the plankton functional types.A diagram showing the major organisms, nutrients, chemical systems, organic matter and fluxes modelled in ERSEM is shown in Fig. 1.ERSEM can be run in parallel with any one of a range of benthic models of varying complexity, the parameterisation used for this publication uses a simple parameterisation of remineralisation where sedimented organic matter is recycled to the water column in inorganic form.
The initial nutrient conditions for nitrate, phosphorus and silicate were taken from the World Ocean Atlas database (Garcia et al., 2010).The initial iron concentrations were zonally averaged interpolations of the iron data from Tagliabue et al. (2012).The iron dust surface deposition climatology was based on Mahowald et al. (2005).The remaining biogeochemical fields were initialised to low concentrations.
In order to be a better representation of nature, marine ecosystem models like ERSEM are typically run in conjunction with an ocean circulation model, such as the NEMO (Madec, 2008).NEMO is a framework of ocean related engines, ocean dynamics, thermodynamics, sinks and sources, and transport.It was designed to be a flexible tool for studying the ocean and its interactions with the other components of the earth climate system over a wide range of space and timescales.The version of NEMO used in this study was version 3.2 and the ocean circulation model was interfaced with the CICE (Hunke et al., 2013).CICE has several interacting components: a thermodynamic sub-model, an ice dynamics sub-model, a vertical and a horizontal transport sub-model.
The coupled NEMO-ERSEM-CICE models were run simultaneously in the ORCA1 1 • global tripolar grid, with 75 fixed depth layers on the UK Met Office super-computing system (MONSooN).The atmospheric boundary conditions were taken from the CORE2 global air-sea flux data set (Large and Yeager, 2009).The coupled model was run for 117 years, from 1890 to 2007.The initial conditions for the circulation model at 1890 were taken from a 60 year physicsonly climatological spin-up from a still global ocean.The first 60 years  of coupled NEMO-ERSEM running were spun-up using climatological CORE2 forcing.After 1950, the remaining 57 years were run with the interannually variable version of the CORE2 forcing.
The simulation was run as part of the iMarNet project: an inter-comparison project of six UK biogeochemical models (Kwiatkowski et al., 2014).The ERSEM model run shown here is an updated parameterisation relative to the ERSEM model data used in that study.A unique requirement of this project was that the six biogeochemical models were coupled to identical physical models.All six biogeochemical models were required to use identical model parameters and settings to run the physical component of the simulation, NEMO.Furthermore, all six models used the same computing resource, the same coding framework and the same initial conditions for nutrients.In other words, the physical parameterisations were prescribed according to specific pre-defined settings and no further changes to the physical settings were permitted.
www.geosci-model-dev.net/9/59/2016/Geosci.Model Dev., 9, 59-76, 2016 biogeochemical models The use of emergent relationships as a tool to assess the quality of a marine biogeochemistry model is the central thesis presented in this work.However, there are some important caveats that should be stated.First, this method is not intended to replace current validation methods, but to complement them.The standard methods such as a point-to-point comparison should remain valid tools for objectively determining model quality.Second, the emergence of a coherent natural relationship in a simulation is an indication that the model has an appropriate representation of the ecosystem functions that create the emergent relationship.However, the ability to reproduce an emergent property does not guarantee the accuracy of the model choices and parameterisation.
Only through a thorough investigation of the model structure and behaviour can we know whether the model reproduces the emergent property for the right reasons.Nevertheless, we aim to demonstrate that the reproduction of emergent relationships by the model are a diagnostic tools, which may be used to identify the origins of model strengths and weaknesses.
In an ideal world, this work would present a comprehensive set of emergent relationships to assist the validation of any marine biogeochemical model.However, models can differ enormously in their complexity, design choices, selection of ecosystem functions implemented, parameterisations, location and scope, such that each unique model would need to be validated using a different set of emergent properties.The emergent relationships shown here are tailored to ERSEM, but many of these relationship may also appear in other models.
The use of emergent relationships to validate models works best under certain conditions.The desired scenario is when the emergent relationship has been observed multiple times in nature with several independent data sets, is valid over large spatial and temporal ranges, and has been reported with a high statistical probability, if fitted to a mathematical function.The ideal scenario regarding the model version of the emergent relationship is that it is not purposefully and obviously imposed a priori by the choice of model parameterisation, it is not a direct extrapolation of the choices made in model design, the emergent property should not be reproduced in the simulation by tuning a small number of parameters, and the model should not be explicitly tuned to match the data.
As emergent relationships are large-scale properties of the system, validation with emergence is possible even when there is very little or no appropriate data for the region under study, or when the physical circulation component of the model differs significantly from that observed in nature at the same location and time.This is one of the strengths of the validation with emergent properties method: it can be used to demonstrate ecosystem model quality in the absence of perfectly simulated physical marine environment or extensive local measurements.Nevertheless, it is important that the model and the data originate from similar marine environments.For instance, emergent property validation cannot compensate for a catastrophic failure of the hydrodynamic model, nor can it be used to validate the model in regions with unusual and understudied behaviour.However, it is not crucial to match up the exact locations in model and data, as used in the point-to-point study in de Mora et al. (2013).
A full understanding of the causal nature of the relationship is not a strict requirement in order to be informative.There are many naturally occurring phenomena for which an explicit explanation is not immediately available, but for which the relationship is nevertheless stable.A well-known example is the Redfield ratio (Redfield, 1934;Arrigo, 2005).Unexplained relationships can inform about the validity of some aspect of the model behaviour.Furthermore, if the model could reproduce a natural emergent property in the absence of a purposefully prescribed functional relationship, then it may be possible to use the model to discover a causal relationship.
The emergence of a coherent natural relationship in a simulation is an indication that the model has a appropriate representation of the ecosystem functions that create the emergent relationship.If the emergent relationship is not seen in the model, this implies that the ecosystem functions that bring about the emergence are not correctly implemented in the model.If the emergent relationship is present in the model, but breaks down under certain conditions, this means that those conditions are not correctly modelled.Such break down points can be a combination of extreme physical, biogeochemical, or ecosystem conditions, and can be used to pinpoint the limitations of the simulation, allowing for powerful and precise criticisms of the model.However, some caution is required as there may also be break down conditions for the emergent property in nature.Nevertheless, every emergent property that the model fails to reproduce is a new direction for future model development.
The remainder of this section is a non-exhaustive list of the emergent properties that can be used to assess the quality of the NEMO-ERSEM iMarNet hindcast.The first emergent property shown here is an extension of the diatom chlorophyll relationship previously described in Holt et al. (2014).

Phytoplankton community structure
In modelling and observational marine biology, phytoplankton are often grouped into size-or function-based classifications (Woodward et al., 2005).The size of a phytoplankter can directly influence a range of physiological processes, which combined together with other individuals can have an ecosystem-wide effect.Physiologically, a cells size may affect its nutrient uptake, metabolic rates, physiology and light absorption (Nelson et al., 1993;Stolte and Riegman, 1995;Huete-Ortega et al., 2012;Zhang et al., 2012).At the  (c) are diatoms and large phytoplankton, nanophytoplankton and picophytoplankton, respectively.The least-squared fit to the three-population absorption model to ERSEM is shown as a green line.The other coloured lines are the five fits from Hirata et al. (2011), Devred et al. (2011), and Brewin et al. (2012, 2014, 2015).The dashed vertical line indicates a typical detection limit of HPLC and SFF methods.ecosystem scale, these individual-level effects combine to influence large-scale observable properties such as the community primary production, export, the food web, and the light penetration depth into seawater (Riegman et al., 1993;Finkel et al., 2010;Huete-Ortega et al., 2012).In addition to size classes, phytoplankton function-based classifications are also used in modelling and observational marine biology.Phytoplankton functional types allow for a grouping of phytoplankton species by their role in the ecosystem, their preferred nutrient sources, and their production, export and sinking rates.
Many marine ecosystem models use a size-or functionbased classification.The phytoplankton classification in ERSEM follows both size-and function-based classifications.There are four plankton functional types (PFTs) in ERSEM.Three of the groups are sized based: nanophytoplankton, picophytoplankton and large phytoplankton.The fourth group, diatoms, are between nanophytoplankton and large phytoplankton in size but include a silicon component.Each PFT has different nutritive affinities and requirements, metabolic rates, and different palatability as a food source for predators.
In both the ecosystem and in models, the relative abundance of each class is referred to as the community structure.The relative abundance is often measured in terms of chlorophyll, but it also can be gauged in units of cell count, total cell volume, or carbon or nitrogen biomass.The community structure in the model is dependent on a large number of factors, including the nutritive affinity and nutritional storage capacity of each PFT relative to each nutrient, the palatability as a food source of each PFT for each zooplankton functional type, and local environmental conditions such as light, temperature and nutritive up-welling.In ERSEM, there is no explicit parameterisation of their absolute or relative abundances; it is a property that arises out of a combination of many ecosystem functions.
The relationship between the abundance of each plankton functional type and the total community chlorophyll has been presented repeatedly with data from both high performance liquid chromatography (HPLC) and size fractionated filtration (SFF) measurements.Five examples of this relationship are Hirata et al. (2011), Devred et al. (2011), and Brewin et al. (2012, 2014, 2015).Each of these fits is prepared using a different data set.Hirata et al. (2011) used multiple HPLC databases from around the world.Devred et al. (2011) used chlorophyll and absorption data from the northwest Atlantic region that was collected between 1996 and 2003.Brewin et al. (2012) used HPLC data in the Indian Ocean taken between 1995 and 2007.Brewin et al. (2014) used SFF from the Atlantic Meridianal Transect (AMT) cruises in the Atlantic Ocean between 1996 and 2012.Brewin et al. (2015) used a aggregation of 16 unique globally distributed databases.
Figure 2 shows the five fits of in situ community structure (Hirata et al., 2011;Devred et al., 2011;Brewin et al., 2012Brewin et al., , 2014Brewin et al., , 2015)).The fit of the ERSEM simulation to the threepopulation absorption model of Brewin et al. (2010) using a least-squared fit is also shown.This fit was performed using model data after the simulation had been completed and did not influence the relationship of the plankton functional types during the simulation.For all three panels, the x axes are the total community chlorophyll, and the y axes are the percentage of the total chlorophyll that came from each PFT.The diatom and large phytoplankton functional types are summed together in the left panel.toplankton and the right panel shows picophytoplankton.The dashed vertical line indicates a typical detection limit of HPLC and SFF measurements of 0.1 mg Chl m −3 .The Devred et al. (2011), and Brewin et al. (2012, 2014, 2015) are also fitted to the three-population absorption model, but the Hirata et al. ( 2011) is fitted with its own community structure model.The difference between the Hirata model and the three-population model is that the Hirata model is a best fit to the data, whereas the three-population model is derived based on empirical principles.Despite their differences, these fits have little variability in the overall shape of the fit between data sets and methodologies.Note that only the fits are shown in Fig. 2, the in situ data itself are not shown, nor is the model data.Showing the fit as a smooth line hides the substantial spread of both the in situ and model data.For instance, the data shown in the widely published Fig. 2 of Hirata et al. (2011) has already been smoothed with a 5-point-running mean, and that running mean is further smoothed to the fit line shown in Fig. 2.
Figure 3 is a three panel plot of phytoplankton community structure showing the model data as a logarithmically scaled two-dimensional data histogram distribution in blue scale.This figure also shows a fit to the model data, and the in situ data fits from Brewin et al. (2015), and Hirata et al. (2011).The fit to ERSEM, the Brewin 2015 and the Table 1.Parameters of the fits to the three-population absorption model for the Brewin et al. (2015), and for ERSEM.Two fits to ERSEM are shown: the first is the fit to the ERSEM data set excluding the polar, shallow and inland seas, and the second includes all these regions to a depth of 40 m.The parameters are C m p, n : maximum piconano chlorophyll, S p, n is the slope for piconano chlorophyll, C m p is the maximum picophytoplankton chlorophyll and S p is the slope for picophytoplankton.The model data in Fig. 3 were taken from the top 40 m of the global ocean, excluding the shallow seas, inland seas, the Southern Ocean from 45 • south and the Arctic Ocean from 55 • north.The Arctic and Southern oceans were excluded from this analysis for two reasons: first, there were no Arctic community structure data in four of the studies mentioned above, and Brewin et al. (2015) only includes a small number of Southern Ocean data.Second, they were excluded because this is a region where the physics in this model run was particularly troublesome.As mentioned in Sect.2, the physics model was fully prescribed as part of the iMarNet inter-comparison project, and no tuning of the physics was permitted.The shallow and inland seas were excluded because a coarse 1 • global model cannot resolve the complex dynamics needed to reproduce these regions.The model data and in situ data were both fit to the same three-population absorption model (Brewin et al., 2010), and both sets of fit parameters are shown in Table 1.
The model reproduces the overall shape of the picophytoplankton community structure in Fig. 3, in that the picophytoplankton contribution to community chlorophyll is higher at low chlorophyll concentrations.However, the picophytoplankton community structure is concave in the model and convex in the three-population model.At low chlorophyll concentrations, the model picophytoplankton commu-65 nity structure is more similar in shape to the Hirata et al. (2011) fit, which is the only in situ data fit from Fig. 2 not based on the three-population model.
The three-population absorption model of Brewin et al. (2010) can be summarised in the following equations.The picophytoplankton component of the total community chlorophyll, chl p , is calculated as where chl is the total community chlorophyll, C m p is the maximum picophytoplankton chlorophyll, and S p is the initial slope of the exponential function for picophytoplankton.Note that Figs. 2 and 3 show the phytoplankton functionaltype contribution as a percentage of total community chlorophyll instead of as the chlorophyll concentration used in Eq. ( 1).
where C m p, n is the maximum piconano chlorophyll and S p, n is the initial slope of the exponential function for piconano.The chl p, n function is not shown explicitly in either Fig. 2 or Fig. 3, but is used in the calculation of the nanophytoplankton, diatoms and large phytoplankton functional-type community structure.
Unlike chl p and chl p, n , the nanophytoplankton chlorophyll, chl n , and microphytoplankton chlorophyll, chl m , are not explicitly fitted.Instead, their contribution to total community chlorophyll is determined from a combination of Eqs. ( 1) and (2).The nanophytoplankton chlorophyll, chl n , shown in Eq. ( 3), is the difference between the piconano group chlorophyll and the picophytoplankton chlorophyll.This is shown in the middle panel of Figs. 2 and 3.
The diatoms and large phytoplankton, (together also known as microphytoplankton) chlorophyll, chl m , are the remainder of the total chlorophyll that is not accounted for by the piconano component.It is calculated as the difference between the piconano chlorophyll and the total chlorophyll.This is shown in the top pane of Fig. 3.
Just to explicitly state the overarching assumption used here, the total chlorophyll, chl, is equal to the sum of the threecomponent functional-type chlorophyll.
Note that the three fits of Fig. 3 are not free to vary independently.For any given value of total chlorophyll on the x axis, the sum of the three populations much be equal to 100 %.Furthermore, the fits are not free to vary to any shape, they are limited by the structure available to Eqs. ( 1), ( 3) and (4).Equations ( 1) and ( 2) are fitted to the simulated ERSEM data using in the least-squares fit.In addition, the fit to ERSEM was performed such that each model point had equal weight, while the underlying histogram distribution is shown with a log scale.This means that the fits are logarithmically more influenced by the high data density regions in this figure.For these reasons, the fits may not appear to match the overall shape of the model distribution, while still being an acceptable fit.

The carbon to chlorophyll ratio in particulate organic matter and phytoplankton
The ratio of carbon to chlorophyll is an important indicator of ecosystem state.This ratio has ecosystem relevance on any scale, from an individual cellular level up to the entire community.Within a given cell, the carbon to chlorophyll ratio is a sensitive indicator of a cells physiological state (Geider, 1987) and is affected by the cells response to irradiation, temperature and nutrient supply.In remote sensing and modelling usage, this ratio also plays a significant role in the calculation of phytoplankton biomass from ocean colour and in the modelling of primary production and in Sathyendranath et al. (2009) and Geider et al. (1997).For these reasons, much effort has been dedicated to the study and modelling of the carbon to chlorophyll ratio in modelling and in situ measurements.
It has long been known that a single value for the carbon to chlorophyll ratio under all environmental conditions is inappropriate for ecological studies (Geider, 1987).In a low-light environment, cells can produce excess chlorophyll relative to their carbon content in order to maximise light harvesting.In contrast, the cells invest in different compounds to protect the fragile parts of the cell from irradiation damage in a high light environment (MacIntyre et al., 2002;Polimene et al., 2014).This behaviour is known as photo-acclimation.Similarly, it has long been known that the carbon to chlorophyll ratio increases with decreasing temperature (Eppley, 1972).While the exact pathway is not fully understood, it has been suggested that this shift is limited by physical properties of the components of photosynthetic apparatus.For instance, an increase in carbon relative to chlorophyll may be required to maintain fluidity at cold temperature, or low temperature might impose constraints on the enzyme catalytic reactions (Geider, 1987).Similarly, each plankton functional type responds in a different way to changes in temperature, affecting the community structure, which may in turn influence the carbon to chlorophyll ratio.
While there are many techniques available for measuring total chlorophyll concentration, the direct measurement of phytoplankton carbon biomass is difficult because of the challenge of separating the phytoplankton from the www.geosci-model-dev.net/9/59/2016/Geosci.Model Dev., 9, 59-76, 2016 other particulate organic matter (Graff et al., 2015).Instead, chlorophyll concentration is typically used as a proxy for phytoplankton biomass.For obvious reasons, comparing biomass derived from total chlorophyll to total chlorophyll is not an independent method for obtaining the carbon to chlorophyll ratio.However, it is possible to measure total community chlorophyll and particulate organic carbon (POC) simultaneously and independently.
A meta-study of simultaneous particulate organic carbon and chlorophyll measurements was published in Sathyendranath et al. (2009).The chlorophyll measurements were divided according to two measurement technologies: high performance liquid chromatography (HPLC) chlorophyll, and Turner fluorometer data.The data used were taken from 16 cruises in the Labrador Sea, the Scotian shelf and the Arabian Sea, between 1993 and 2001.Then they used the accumulated data to produce a fit to the following relationship: where POC is the particulate organic carbon, chl is the total community chlorophyll and m and p are the fitted parameters.
In order to study the relationship between phytoplankton carbon and total chlorophyll, Sathyendranath et al. (2009) applied the assumption that "at any given chlorophyll concentration, the lowest particulate carbon observed represents the phytoplankton carbon associated with that chlorophyll concentration."Then, they used a 1 % quantile regression to determine a relationship between the lowest particulate organic carbon observed for each total chlorophyll concentration.Although the methodology was more complex, they produced a fit that was mathematically similar to Eq. ( 6): where C is the phytoplankton carbon, chl is the total community chlorophyll and n and q are the fitted parameters.ERSEM employs a variant of the Geider model (Geider et al., 1997(Geider et al., , 1998) ) for describing the carbon to chlorophyll ratio in phytoplankton.Each PFT in ERSEM has a different parameterisation of the Geider model, allowing for phytoplankton carbon to chlorophyll behaviour to differ for each plankton functional type and to vary according to local conditions.In practice, the range of possible values of the phytoplankton carbon to chlorophyll ratio has a minimum and maximum value for each functional type, but each phytoplankton carbon to chlorophyll ratio is free to vary inside that envelope.There are no limits imposed on the community phytoplankton carbon to chlorophyll ratio.
The total community chlorophyll in ERSEM is the sum of the chlorophyll concentration of each phytoplankton functional type.The modelled POC is the sum of the carbon components of all four phytoplankton functional types, all three zooplankton functional types, and the particulate detrital carbon fields.Each one of these carbon compartments is free to vary independently of total chlorophyll.For these reasons, the ratio of POC to community chlorophyll is not an obvious and purposefully prescribed consequence of model choices, and cannot be tuned with a small number of parameters.The ratio of the particulate organic carbon against the total chlorophyll is not explicitly restrained or parameterised in any way.
Figure 4 shows the relationship between particulate organic carbon and total chlorophyll.The model data are shown as a logarithmically scaled two-dimensional histogram distribution with in blue scale, and the model distribution is taken as all model points in the top 40 m of the monthly climatology of the final 10 years of the simulation, excluding shallow seas, inland seas, the Southern Ocean from 45 • south, and the Arctic Ocean from 55 • north.A fit of the models relationship between particulate organic carbon and total chlorophyll to Eq. ( 6) is also shown as a full green line.The two Sathyendranath et al. (2009) fits of POC to chlorophyll are also shown as full lines and the two 1 % quantile regression fits representing phytoplankton carbon against chlorophyll are shown as dashed lines.The parameters of these fits are shown in Table 2.
The two 1 % quantile regressions fits are included in Fig. 4 because they indicate a theoretical lower bound for the modelled POC : Chl field.Approximately 3 % of the model data fall below the theoretical lower limit indicated by the 1 % Sathyendranath et al. (2009) Turner line.This is an acceptable fraction, as the 1 % quantile regression was an arbitrary cut-off point for the minimum POC concentration for each chlorophyll range in the data.However, the data below this line occur in the model only on the edges of the Arctic and Southern oceans in the winter.When the entire model domain is included down to 40 m deep, including the Arctic and Table 2.The parameters of the particulate organic carbon to total chlorophyll and phytoplankton carbon to chlorophyll fits to the Eqs.( 6) and (7).Southern oceans, and inland and shallow seas, the fraction of data becomes as high as 11 %.
Figure 5 shows the total phytoplankton carbon against total chlorophyll.This figure also shows the model data as a two-dimensional histogram distribution logarithmically scaled in blue scale and a least-squares fit of the data to Eq. ( 7).In a modelling analysis, it is unnecessary to apply the 1 % quantile regression, as it is straightforward to extract phytoplankton carbon.Figure 5 also shows the two 1 % quantile regression fits from Sathyendranath et al. (2009) as dashed lines as in Fig. 4.However, the model data shown in this figure is phytoplankton carbon against total chlorophyll; accordingly, the model data and fit are expected to match these dashed lines.As before, the model distribution is the top 40 m of the monthly climatology of the final 10 years of the simulation, excluding shallow seas, inland seas and the Arctic and Southern oceans.Both the fits and the bulk  2009) phytoplankton carbon to total chlorophyll relationship.This indicates that the combination of temperature, light and nutrients in the modelled global ocean result in a reasonable and natural response of the modelled phytoplankton.This natural response in the phytoplankton indicates that the combination of model parameters describing the phytoplankton community in this simulation are a valid parameter set for the modelled physical and biogeochemical environment they inhabit.

Particulate organic carbon and particulate organic nitrogen
The ratio of carbon to nitrogen in the ocean has long been a historically important measurement and was first identified by Redfield (1934).The interplay of carbon and nitrogen has since been a major component of our understanding of ocean biogeochemistry, underlying modern theories of nutrient cycles, nutrient limitation in phytoplankton and stoichiometric variability.The balance of carbon to nitrogen in the ocean has been historically measured in the ratio of carbon to nitrogen in particulate organic matter (POM).
A meta-study of the ratio of particulate organic carbon (POC) against particulate organic nitrogen (PON) was published in Martiny et al. (2013).Their paper used over 40 000 globally distributed POC-PON paired samples from 5383 unique stations in the upper 200 m of the ocean water column.They produced a frequency distribution of POC against PON in the global ocean.The power of this figure was that it was a clear demonstration of the discrepancy between a modern understanding of the stoichiometric balance of carbon to nitrogen and the canonical Redfield ratio.Instead of a static Redfield ratio, the observed POC : PON ratio varied between 2 and 20, and the canonical Redfield ratio was closer to the median of the data set than to the mean or the mode.
The POC : PON frequency distribution from Martiny et al. ( 2013) was reproduced using simulated data in Fig. 6.This figure shows the POC : PON ratio calculated in the model, not only the POC : PON ratio of the in situ data from Martiny et al. (2013) but also the canonical Redfield ratio.Both the model and the in situ histograms in Fig. 6 were normalised to unity area.A summary of the statistical analysis of Fig. 6 is shown in Table 3.The mode shown in Table 3 was calculated by finding the most populous group after binning the model ratio in bins of width 0.1.The bin widths used to calculate the mode are finer than those shown in Fig. 6.
The POC of the model in Fig. 6 was calculated as the sum of the carbon components of all phytoplankton functional types, all zooplankton functional types and the particulate detritus groups.These are the same groups that were used to calculate POC in Sect.3.2.The PON was calculated as the sum of the nitrogen components of the same groups.An artificial detection limit of 0.1 µmol m −3 was applied to the modelled PON component of the ratio.The model data were a monthly climatology of the final 10 years of the simulation (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007), excluding the Arctic Ocean and all model data below 200 m.Unlike the data selections of Sect.3.1 and 3.2, the Antarctic domain, inland seas and shallow seas are included in this data.This is because the Martiny et al. (2013) data set does not contain any measurements from the Arctic Ocean, but it does include many data from the Southern Ocean, inland seas and shallow seas.
In nature as in models, particulate organic matter is usually composed of a combination of phytoplankton, zooplankton and detritus.Experimentally, POM is typically measured as all organic matter over a certain size that can accumulate in a filter.In ERSEM, POM is the sum of the all four phytoplankton functional types, all three zooplankton functional types, and the particulate detrital fields.Each phytoplankton functional type has an internal variable stoichiometry and the ability to accumulate a luxury buffer of nitrogen.The micro-zooplankton and heterotrophic nanoflagellates also have variable stoichiometry that are influenced by their food source.The mesozooplankton has fixed stoichiometry and exudes the excess nitrogen back into detritus.There are six detritus size classes in ERSEM: three dissolved and three particulate classes, but the dissolved organic matter fields did not contribute to the POM shown here.

Intracellular elemental stoichiometry
Stoichiometry is the balance of each element in organisms and in the ecosystem (Sterner et al., 2002).As mentioned previously, Redfield observed co-variability in the concentrations of dissolved nitrate and phosphate in seawater and in the composition of plankton (Redfield, 1934).While further co-variation has since been observed, considerable variability in the balance of elements in particulate organic matter and intracellular material has also been observed (Martiny et al., 2013).The ratio of each element against carbon has been shown to vary significantly between region, taxa, ecosystem role and physiological status (Moore et al., 2013).The co-variability of nutrients and carbon in the ocean is closely linked with many important metrics of ecosystem behaviour, such as primary productivity, community structure, export and growth rates, and nutrient limitation.
Despite the significance of marine nutrient cycles, colimitation by two or more nutrients is still poorly understood, and appears infrequently in models (Moore et al., 2013).Furthermore, variable stoichiometry and co-limitation are required features in order to represent the spatial distribution of nutrient limitation.ERSEM is one of the few models that does implement variable stoichiometry.In addition to its carbon cycle, ERSEM resolves four nutrient cycles: nitrogen, phosphorus, iron, and silicon.All four nutrients can become limiting or co-limiting for any given phytoplankton functional type, with the exception that only the diatoms interact with silicon.In addition, nitrogen, phosphorus and iron uptake and limitation in ERSEM are based on the Droop model (Droop, 1973), which uses the internal nutrient concentration to carbon ratio rather than external inorganic nutrient concentrations to determine phytoplankton behaviour.The silicate limitation and uptake for diatoms is computed from the external availability of dissolved silicate, based on Michaelis-Menten kinetics (Johnson and Goody, 2011).
In their meta-study, Moore et al. (2013) compiled multiple data sets of simultaneous measurements of particulate organic carbon and nutrients.To visualise this data, they presented a comparison of the ratio of each element against carbon.For each nutrient, they plotted the typical organic nutrient : carbon ratio on the x axis against the typical inorganic nutrient : carbon ratio on the y axis.Figure 7 shows a ver-Figure 7. A comparison of the ratio of each modelled nutrient to carbon ratio in organic matter against the dissolved inorganic nutrient to carbon ratios.This figure shows a colour-coded distribution of the modelled nitrogen, phosphorus, iron and silicon to carbon ratios.The model distribution means are indicated by circular markers, and the typical in situ value and observed range from Moore et al. (2013) are shown as square markers with horizontal bars.The vertical bars associated with the square markers were calculated from World Ocean Atlas data (Garcia et al., 2010) and GEOTRACES (Henderson et al., 2007).sion of this figure produced using the four ERSEM nutrients.This figure shows the model distribution with the mean of the model data in circular markers, and the typical in situ value and observed range from Moore et al. (2013) in square markers with black bars.
The organic component of Fig. 7 was calculated as the ratio of particulate organic nutrient to particulate organic carbon for each nutrient.The POC was calculated in the same way as in Sect.3.2 and 3.3: where POC is the sum of the carbon components of all phytoplankton functional types, all zooplankton functional types and the particulate detritus groups.The particulate organic nitrogen, phosphorus and iron were calculated as the sum of the nitrogen, phosphorus and iron components of the same groups.Although ERSEM includes four pelagic silicon fields (diatom silicon, inorganic silicate, and medium and large particulate detrital silicon), the calculation of particulate organic silicon follows the methods of Moore et al. (2013), and only uses the silicon component of diatoms.
As previously described, the modelled ratio of particulate organic nutrient to carbon can vary according to local conditions such as light, temperature and predation.For all of these interlocking and competing components to reproduce the stoichiometric variability of the global ocean, all phyto- All observation data are taken from Moore et al. (2013) except: a indicates values from Martiny et al. (2013), which were not taken into account in Moore et al. (2013); b indicates data calculated using World Ocean Atlas data (Garcia et al., 2010); and c indicates data calculated using GEOTRACES data (Henderson et al., 2007).plankton functional types, zooplankton functional types and detrital fields and their interactions need to be parameterised sensibly.The overall particulate organic matter stoichiometry in ERSEM is not susceptible to tuning via a small number of parameters.
The inorganic component of the model data were taken directly from the output of the simulation.The dissolved inorganic carbon (DIC) cycle in the model is described in Artioli et al. (2012).There are no explicit limitations of the upper or lower limits of inorganic nitrogen, phosphorus or silicon in the model.However, there is a soft cut-off for iron at 0.6 nM to take into account for hydroxide precipitation (Aumont et al., 2003) and a firmer upper limit of 2 nM to take into account for the saturation concentration.
The minimum, maximum and mean values of each ratio against carbon in Fig. 7 are shown in Table 4.This table includes the typical Moore et al. (2013) data mean and range, and the data from the model run.Moore et al. (2013) did not include any measure of variability in the inorganic data; the inorganic range shown in Fig. 7 and Table 4 was included here as another test of the model.The variability in the inorganic axis is due to the co-variability in the dissolved inorganic carbon and the dissolved inorganic nutrient.However, the range of variability in DIC is usually of the order of 15 %, but dissolved inorganic nitrogen, phosphorus, iron and silicon can vary by several orders of magnitude.This means that most of the variability in the y axis is mostly due to the inorganic nutrient, not the DIC.The ranges of the ratio were estimated using the Moore et al. (2013)  constant value for the DIC, and the minimum and maximum values of the nutrients contribution to the inorganic nutrient to DIC ratio were taken from data from the World Ocean Atlas (Garcia et al., 2010) for nitrogen, phosphorus and silicon, and GEOTRACES for iron (Henderson et al., 2007).This means that the ranges are not the most extreme values ever recorded, but reflect the most extreme values of inorganic nutrient concentration seen on a gridded climatological data product.Furthermore, the range of variability in the organic nutrient : carbon ratios might not reflect the most current knowledge.As stated in Moore et al. (2013), "maximum and minimum values will typically correspond to nutrient replete or limited cultures respectively and ranges could potentially be extended through observations of other taxa and growth conditions."For instance, Moore et al. (2013) cited an observed maximum phytoplankton N : C quota of 0.169, but the Martiny et al. (2013) data set has observational data all the way down to their cut-off point, which was a C : N ratio of 2 : 1.Note that Moore et al. (2013) used the N : C ratio, but Martiny et al. (2013) used a C : N ratio.Also note that the model data do not show any histogram distribution, only presence/absence.The presence/absence model data in Fig. 7 is extracted from the top 40 m of the global model between 1997 and 2007 and excludes the shallow seas, inland seas and the Arctic and Southern oceans.The model data were not further sub-sampled to match the sampling locations of the in situ data.This data have a monthly time resolution and are not climatologically averaged.
In addition, it has been shown by Karl et al. (1998) that measurements of particulate organic matter may be exaggerated by contributions from dissolved organic matter.Whereas the ratio calculated in the model does not include any contribution from dissolved organic matter.

Discussions
It has been shown that the ERSEM global hindcast successfully reproduced many natural behaviours of the ecosystem.Each of these behaviours has covered a different aspect of ecosystem function, and when combined together they illustrate the potential of model validation with emergent properties and ecosystem function.

Phytoplankton community structure
In the top panel of Fig. 3, there is a cluster of points where the diatom and large phytoplankton functional types unexpectedly dominate the community structure at low total chlorophyll.These points account for less than 0.2 % of the data, after the cuts described above were applied.Furthermore, they only appear adjacent to the excluded shallow and polar seas.While there is also some evidence of the proportion of large cell increase in the polar regions (Sosik and Olson, 2002), we postulate that it is a combination of multiple factors that creates this excess microphytoplankton.First, in the polar regions there is an abundance of nutrients, and especially silicon, caused by excessive mixing in the physical model.Second, the model is parameterised to favour diatoms in lowlight regions.These factors collude to create an abundance of diatoms and large phytoplankton at low total chlorophyll.When the polar, shallow and inland sea regions are included in the model, the number of points included in these regions of the figure increases.As an example, the fit to the threepopulation model was performed to all the model data from the top 40 m of the surface.The results of this fit are shown in the ERSEM (top 40 m) column of Table 1.Relative to the Brewin et al. (2015) fit, the fit of the ERSEM model data to the three-component model shows an overabundance of diatoms and large phytoplankton and an underestimate of picophytoplankton at almost all chlorophyll concentrations.Nevertheless, it is important to stress that the three-population community model is an appropriate emergent property for open ocean outside the polar regions.
Similarly, both Figs. 2 and 3 show that the model also has a modest surplus of diatoms and large phytoplankton at low chlorophyll concentrations in this simulation, which coincides with a low proportion of picophytoplankton.It is likely that this fault is caused by the same factors that cause the microphytoplankton PFT to dominate the community in a small number of cases.However, it is likely to be mitigated in most of the ocean by a lower silicate concentrations leading to slightly stronger silicate limitation for diatoms.
The model data were limited to the top 40 m of the surface ocean, and the relationship breaks down in Arctic waters in the model.Hints of the breakdown of the community structure appear in the in situ data (Brewin et al., 2015), but are seen clearly in the model.It became clear that the large phytoplankton and diatom functional types are in excess at low chlorophyll concentrations in the modelled community structure.In addition, the picophytoplankton chlorophyll as a fraction of the community behaves much more like the concave empirical fit to data seen in Hirata et al. (2011) than like the convex three-population model of Brewin et al. (2010) at low chlorophyll concentrations.
Despite these limitations, ERSEM was very successful at reproducing the overall shape of the community structure and natural balance of phytoplankton abundance between the four PFTs.This means that the combination of the nutrient affinity, growth rates, photosynthetic behaviour and predation rates ecosystem functions were modelled appropriately enough to bring out a natural emergent community structure.This is a robust and well-known emergent property that can be reproduced successfully by ERSEM, despite the problems in the prescribed physical simulation.

The carbon to chlorophyll ratio
Together, the particulate organic carbon to chlorophyll and phytoplankton carbon to chlorophyll ratios from Sect.3.2 demonstrate that the phytoplankton biomass forms an appropriate fraction of the particulate organic matter.This means that the balance of producers to the rest of the organic matter, including zooplankton and detritus, is similar to nature over the range of observed total community chlorophyll.
In both Figs. 4 and 5, there is a region where the phytoplankton carbon is much less than the Sathyendranath et al. (2009) fit, and this group of data appear to have a different slope than the bulk of the data.The phytoplankton carbon in the data points are almost entirely composed of diatoms and large phytoplankton, near the polar regions, and account for approximately 2.5 % of the data set.Similarly to Sect.3.1, this region highlights that the issues of the abundance of silicate caused by excessive mixing, the favouring of diatoms in low-light regions and the relatively low grazing pressure on microphytoplankton from zooplankton at low phytoplankton biomass concentrations.When the polar, shallow and inland sea regions to a depth of 40 m are included in the analysis, the number of points included in these regions increases up to approximately 11 % of the data set.This is another indication that either the model has not captured the behaviour of the high latitude regions, or the emergent property breaks down in these regions.Unfortunately, Sathyendranath et al. (2009) did not include any data from Polar regions in the winter that could be used to test this hypothesis.
The carbon to chlorophyll relationship has many knockon effects in the model: it influences the entire carbon cycle and has a huge impact on the calculation of total global primary production.When the model successfully reproduces the carbon to chlorophyll ratio in a global ocean simulation, this is an indicator that it has a good enough approximation of the roles of light, temperature and nutrient limitation in each of the plankton functional types.The fact that the model reproduces the natural range of behaviours of the carbon to chlorophyll ratio highlights that the ecosystem model functions appropriately over a range of environments.However, there may be other ecosystem functions that are currently absent in the model that could affect this ratio, such as the presence of higher trophic levels, anthropogenic detritus sources, or an improved model of photo-inhibition.

Particulate organic carbon and particulate organic nitrogen
Figure 6 shows a comparison of the ratio of particulate organic carbon to particulate organic nitrogen in the model and in an distribution of in situ measurements.In ERSEM, none of the detritus fields have any limitations on their stoichiometric variability.This means that the models POC : PON ratio can vary according to local conditions and predation, and the overall particulate organic matter stoichiometry in ERSEM is not susceptible to tuning via a small number of parameters.In order for all of these interlocking and competing components to reproduce the POC : PON ratio variability in the global ocean, it requires all the phytoplankton functional types, zooplankton functional types and detrital fields to be balanced and healthy.While the model captures the central tendency of the in situ data, it does not capture the range of observed POC : PON ratios or the shape of the distribution seen in the data.The model underestimates the frequency of POC : PON ratios below 5 and above 6.5.It is possible that some of this difference is due to spatial bias and uneven sampling of the in situ data.In that case, it may be possible to capture more of the shape of the Martiny et al. (2013) data by sub-sampling the model data to match the distribution used to produce their data.However, the model data shown here is a monthly mean of a 1 • by 1 • square of ocean.The variability seen when taking an instantaneous measurement of the concentration in a bottle of seawater will always be more extreme that the mean value of a 1 • by 1 • square of ocean.Second, in this work, we attempt to validate the models ecosystem function over a large scale without the use of point-to-point matching.
The model's narrow POC : PON distribution is reflected in its the standard deviation (0.61), which is much lower than that seen in data (2.46).However, the range of stoichiometric variability seen in measured POC : PON data are underestimated in the model.Furthermore, there are precisely zero model data with a POC : PON ratio below 4.3 or above 16.5, whereas the Martiny et al. (2013) data range from 2.0 to 20.This behaviour in the model is linked to the fixed maximum luxury buffer of nitrogen relative to carbon in all the phytoplankton functional types.This maximum nitrogen buffer translates to a fixed minimum value of the POC : PON ratio, which is maintained as it cascades through the trophic levels.On the other end of the scale, there is a minimum requirement of nitrogen to carbon, below which the phytoplankton are nitrogen limited and do not grow.
One downside of the simultaneous analysis of all PFTs is the introduction of the risk of compensating errors.This occurs when an error in the C : N ratio from one group is balanced by a similar but opposite error in another group.For this reason, we recommend that model validations also check each the C : N ratio for each PFT individually.
While the model does not reproduce the standard deviation or the tails of the distribution seen in data, the ERSEM simulation was particularly successful at reproducing the mode of the Martiny et al. (2013)  of particulate organic detritus, all of which have variable stoichiometry (except mesozooplankton), is a sign of the validity of models parameterisation of the balance of carbon to nitrogen.

Intracellular elemental stoichiometry
The stoichiometric variability of particulate organic matter against the ratio of inorganic nutrients to DIC was shown in Fig. 7.This figure showed the typical organic nutrient : carbon ratio on the x axis against the typical inorganic nutrient : carbon ratio on the y axis.In addition to the typical values and the models distributions, this figure also showed a range of values that have been observed for each element.This figure demonstrated that the range of stoichiometric behaviours present in the model match those measured in situ.
The ratios of the inorganic nutrients to carbon in the model do not typically get as low as the estimates of this ratio, except iron.This is a problem with the model parameterisations that has also been seen in Sect.4.3, which will need to be addressed in future parameterisations.This is an example of the use to emergent property validation as a tool to direct future model development efforts.
The ratio of nitrogen to carbon is shown in blue in Fig. 7.The model captured the mean organic N : C ratio, but had a wider range of values than that quoted by Moore et al. (2013).However, the maximum value of the N : C ratio has been extended from 0.169 to 0.5 after the Martiny et al. (2013) results were included.The model underestimated both the mean inorganic ratio and the range of variability in the inorganic N : C ratio.The model appears to have a fixed lower limit of the dissolved inorganic nitrogen : DIC ratio that is higher than the minimum nitrogen : carbon ratio estimated from the World Ocean Atlas (WOA) data set (Garcia et al., 2010).
The mean organic phosphorus : carbon ratio is overestimated by the model, but the mean inorganic P : C ratio is underestimated.The range of the inorganic and organic P : C ratios were underestimated by the model relative to the Moore et al. (2013) data.However, both the organic and inorganic phosphorus in the model show a wide range of behaviours, reflecting those seen in nature.Similarly to the nitrogen case, the lowest values of the models dissolved inorganic phosphorus to DIC ratio is higher than the lowest values seen in the estimated DIP : DIC range.This means that the inorganic phosphorus in the model never gets as depleted as the minimum observed ratio estimate.However, the range of the inorganic P : C ratio from the WOA estimates has no indication of the frequency distributions of the P : C ratio range, so it is entirely possible that this extremely low inorganic P : C ratio is also relatively rare.In addition, the model data are the mean of a 1 • by 1 • patch of ocean, whereas the observational data originates from the mean of a 1 litre bottle, which would imply that we can expect fewer extreme values in the model.The modelled P : C ratio also illustrates the effect of external resource limitation: as the inorganic P : C ratio decreases, the organic P : C ratio simultaneously decreases.This figure cannot indicate whether the drop in the organic P : C ratio occurs in the phytoplankton, zooplankton, or detritus, or in some combination of all three groups.Due to the trophic cascade of the POC : PON ratio described in Sect.4.3, we postulate that this effect is caused by the modelled phytoplankton becoming nutrient stressed in low phosphorus environments.
The ERSEM model has four pelagic silicon fields: diatom silicon, inorganic silicate, and medium and large detritus silicon.The Si : C ratio in Moore et al. (2013) and in Fig. 7 are strictly limited to diatoms; there are no quotas associated with particulate organic silicon in other components of the ecosystem shown here.The modelled ratio of silicon to carbon, shown in purple in Fig. 7, captured the range of variability in the inorganic version of the ratio.The range of the dissolved inorganic silicon to DIC ratio was estimated from World Ocean Atlas, Garcia et al. (2010), and is shown as a vertical purple dashed line in Fig. 7.As only diatom silicon are included in this figure and ERSEM has very little variability in the silicate stoichiometry for diatoms, there is no variability in the organic component for silicate.There is only a very slim range of Si : C ratios allowed in the modelled diatoms because the diatoms Si : C quota is set close to the optimal Si : C quota.External silicate limitation directly reduces carbon assimilation and respiration losses may cause small fluctuations in the diatoms Si : C quota (Butenschön et al., 2015).This means that the nutrient stress effect seen in the P : C ratio is not seen in the Si : C ratio.Instead of regulating the internal quota at low inorganic Si : C ratios, diatom growth is reduced and the community structure changes to disfavour diatoms.
The mean organic iron : carbon ratio in the model is lower than the same ratio in Moore et al. (2013): the model underestimated the mean organic ratio by an order of magnitude.However, the Fe : C ratio is the only inorganic nutrient : carbon ratio where the model captured the estimated range.While there is an atmospheric iron source from dust, the model does not include any atmospheric, riverine or hydrothermal sources of nitrogen, silicon or phosphorus.The nitrogen, silicon and phosphorus shown in this paper have been circulated, consumed and recycled for more than 100 simulated years and the relationship between organic and inorganic nutrients, and nutrients against carbon are still all representative of nature.On the other hand, the iron cycle is nudged towards what is observed in nature by an climatological surface deposition, and through hydroxide precipitation and saturation removal of excess iron.This means that the distribution of inorganic iron is not an emergent property of the model, but rather a tuned outcome.These nudges are needed because the iron cycle of ERSEM is much less complex than that seen in a nature.An example of a more natural iron model is Tagliabue et al. (2009), which has three bioavailable forms of iron and two complexed forms of iron.Despite this, as the inorganic Fe : C ratio in our model de-creases, the organic ratio also decreases, indicating that the organic community become increasingly nutrient stressed in low iron environments.
Anthropogenic nutrient loading is expected to increasingly influence nutrient cycles in the ocean and this may lead to shifts in the nutrient balance (Paerl, 1997).Unfortunately, this model does not include any anthropogenic nutrient loading, or indeed any flux of riverine nutrients, and the iron dust deposition is forced with an annual climatology.
Overall, Fig. 7 informs one about the relationship between the inorganic and organic component of the stoichiometric balance.Effectively, it illustrates whether nutrient limitation and nutrient stress are parameterised in a way that reflects nature.Much of the modelled organic matter appears to be iron poor and phosphorus rich relative to nature.The model never captures the lowest dissolved inorganic nitrate, phosphate, or silicate concentrations.It might be expected that the model will produce a wider range of quotas than the historic data sets as the ocean is vastly under-sampled relative to the model.On the other hand, the model is the mean of a 1 • by 1 • patch of ocean, and the data are typically the mean of a 1 litre bottle, which would imply less variability in the model.Furthermore, some of the in situ data may originate in coastal data sets, which have a higher spatial variability than would be seen in a coarse global model.

The role of emergent properties
Combined together, these relationships have provided information about community structure and balance of C : Chl in phytoplankton, the ratio of POC : PON in particulate organic matter, the stoichiometric flexibility of POM and dissolved inorganic nutrients.While some selection cuts have been necessary to reduce the impact of non-physical behaviour, the combination of the relationships can be used to validate the ecosystem model without relying on the model to reproduce an historic measurement at exactly the right place and time.These methods should allow for a validation of the behaviour of the BGC model irrespective of the quality of localised features of the physical model.In this way, these tools compliment current validation methods, such as the point to point, that may not function ideally in an inappropriately parameterised physical ocean model.Furthermore, many of the features seen here, such as the community structure, would not be explicitly apparent in a point-to-point comparison of model to observations.
While it has since expanded beyond its original remit, the European Regional Seas Ecosystem Model was originally built as a model for simulating temperate shelf seas.This work confirms the work of Vichi et al. (2007) by demonstrating using emergent properties that many of ERSEM's design choices and parameterisation are still appropriate in a global context.Furthermore, these emergent relationships were not explicitly parameterised in model development; all of them arise naturally out of a combination of many other ecosystem functions.
The combination of these well-known phenomena have allowed for a test of the majority of the modelled fields throughout the surface ocean.However, these relationships do not cover all aspects of the model.They do not give information about the food web such as the balance of zooplankton and detritus functional types to each other and to the phytoplankton functional types, or about the bacterial community or benthic environment in the model.Also these relationship do not cover important fluxes such as primary production, air-sea gas exchange, or export from the surface layers.
The nutrient cycles of carbon, nitrogen, phosphorus and iron are all influenced by the bacterial loop.The bacterial biomass does not contribute their biomass to the calculation of particulate organic matter used in Sect.3.3 or 3.4, but the bacterial functional type competes with the phytoplankton for the inorganic nitrogen and phosphorus.The bacteria is an additional food source for zooplankton whilst also competing with the zooplankton by scavenging non-silicon particulate detritus.In addition, the bacterial group only excretes to the dissolved organic matter detrital fields.This means that the bacterial functional type wields influence over a significant portion of the pelagic model.Unfortunately, this work could not locate an emergent property to investigate the bacterial behaviour in the model.Some potential avenues where emergent properties may be found in the future include the ratio of primary production to bacterial production or the bacterial growth efficiency.
The models grazers biomass were implicitly included in the POC : PON, POC : Chl and the stoichiometric relationships.However, there are no metrics included here to study the zooplankton by themselves.The authors are not aware of any stable global emergent property for describing the grazer community.
These emergent relationships were selected to reduce the impact of spatial biases, but these relationships may still be influenced by uneven in situ data spatial and temporal coverage.This bias could potentially be resolved by using a pointto-point analysis for the emergent properties; however, this may limit the scope and the power of the emergent property validation.These emergent properties also require the assumption that the property can be extended to cover the entire ocean.Some emergent properties have not been tested in shallow or high latitude seas, and may not hold across all marine environments.

Conclusions
Ecosystem relationships are coherent structures, patterns and properties that are observed to be robust in nature and can be reproduced by a sufficiently complex model.As ecosystem functions arise independently of local physical conditions of the ocean, they can be used to demonstrate model quality in the case of when the physical features of a sea are not co-located in the model and in nature.While these methods cannot compensate for a catastrophic failure of the circulation model, these methods may overcome limited weaknesses in the ocean physics such as a displacement in ocean fronts or the mixed layer depth.These limited uncertainties have the potential of spoiling a point-to-point metric even though they may only be a minor error in the overall picture.
Most importantly, an explicit analysis of the reproduction of ecosystem functions by the model is the only way to demonstrate the models capacity to represent ecosystem function, as opposed to a demonstration of quantitative metrics of absolute ecosystem state.None of the features shown here would be visible in a model to data comparison of static historical concentration measurements.For these reasons, ecosystem functions are a critical tool for the validation of marine ecosystem models.
In future works, it should be possible to use emergent properties to test hypotheses during biogeochemical model development.As an example, if some process X is expected to influence an emergent property Y, a comparison of the emergent property Y in the presence and absence of X may yield insight as to the value of that process in models.Similarly, a comparison of the emergent property Y in two models with alternative parameterisations for process X may facilitate the selection of a parameterisation.A second potential avenue for future research would be to test the independence of the biogeochemical model from the physical model by comparing the emergent properties from the same BGC model coupled against multiple physical models.Finally, it is important to remember that each emergent property that the model fails to reproduce is a new direction for future model development.

Figure 1 .
Figure 1.Schematic representing the major organisms, nutrients, chemical systems, organic matter and fluxes modelled in ERSEM.Blue connectors represent inorganic carbon, red represents nutrient fluxes, black represents predator-prey interactions and green represents the non-living organics.

Figure 2 .
Figure 2. Phytoplankton community structure for all fits including the fit to ERSEM.The panes labelled (a), (b) and (c) are diatoms and large phytoplankton, nanophytoplankton and picophytoplankton, respectively.The least-squared fit to the three-population absorption model to ERSEM is shown as a green line.The other coloured lines are the five fits from Hirata et al. (2011),Devred et al. (2011), and Brewin et al.  (2012, 2014, 2015).The dashed vertical line indicates a typical detection limit of HPLC and SFF methods.

Figure 3 .
Figure 3. Phytoplankton community structure.The model data are shown as the logarithmically scaled two-dimensional data histogram distribution in blue scale.A least-squares fit of the model data to the three-population absorption model of Brewin et al. (2010) is shown as a full green line, and the fit of historic in situ data to the three-population absorption model from Brewin et al. (2015) is shown in a purple line.A fit to data from Hirata et al. (2011) is shown in a black line.The dashed vertical line indicates a typical detection limited of HPLC and SFF methods.
Brewin et al. (2015) ERSEM ERSEM (top 40 m) are identical in Figs. 2 and 3.For clarity, the other fits are not shown in Fig. 3.For all three panels, the shared x axis is the total community chlorophyll, and the y axes are the percentage of the total chlorophyll that came from each PFT.The model data are shown as a logarithmically scaled two-dimensional data histogram distribution in blue scale.The diatom and large phytoplankton functional types are summed together in the top panel.The middle panel shows nanophytoplankton and the bottom panel shows picophytoplankton.The dashed vertical line indicates a typical detection limited of HPLC and SFF measurements.Note that only the fits to the in situ data are shown, the in situ data itself are not shown.

Figure 4 .
Figure 4.The ratio of particulate organic carbon to total chlorophyll.The model data are shown as the logarithmically scaled twodimensional data histogram distribution in blue scale.The full lines indicate the two Sathyendranath et al. (2009) fits to data, and a fit of the model to Eq. (6).The dashed lines show the two 1 % quantile regression fit from the data and they indicate a theoretical lower bound for the modelled POC : Chl field.
* These values were calculated using a 1 % quantile regression.

Figure 5 .
Figure 5.The ratio of phytoplankton carbon to total chlorophyll.The model data are shown as the logarithmically scaled twodimensional data histogram distribution in blue scale.The dashed lines show fit to data, a fit of the model to Eq. (7), and the results of the two 1 % quantile regression fromSathyendranath et al. (2009).

Figure 6 .
Figure 6.The ratio of particulate organic carbon to particulate organic nitrogen in the Martiny et al. (2013) in situ data set and in the model.The Redfield ratio is also shown as a red vertical line.The model data were taken from a monthly climatology of the top 200 m of the final 10 years of running, excluding the Arctic Ocean.Both the ERSEM and the Martiny et al. (2013) histogram were normalised to unity area.
POC : PON ratio.This means that the most common values in the modelled POC : PON ratio are the same as the most common values in the in situ measurement of POC : PON.The reproduction of the mode of the data set by the ERSEM model is a strong indication that the most common behaviour of the POC : PON relationship is appropriately simulated.The ability to reproduce this ratio from the combination of four phytoplankton functional types, three zooplankton functional types and three classes www.geosci-model-dev.net/9

Table 3 .
Martiny et al. (2013)e the POC : PON distribution as reported byMartiny et al. (2013), compared to the result of this study and to the canonical Redfield Ratio.
of the model data distribution in this figure coincide with the two Sathyendranath et al. (

Table 4 .
Table showing the typical, minimum and maximum organic and inorganic ratios against carbon for Moore et al. (2013) and for the ERSEM simulation.
They allow us to demonstrate how natural behaviour emerges from the As they are well-established functional relationships that hold true over large regions of the global ocean, they are a valuable tool for validating ecosystem model in data sparse regions.