
\documentclass[gmdd, online, hvmath]{copernicus}

\begin{document} %\hack{\sloppy}

% \title{The role of ecosystem function and emergent relationships in the assessment of global marine ecosystem models: a~case study with ERSEM}
\title{The assessment of a global marine ecosystem model on the basis of emergent properties and ecosystem function: a case study with ERSEM}

\Author{L.}{de Mora}
\Author{M.}{Butensch\"{o}n}
\Author{J.~I.}{Allen}

\affil{Plymouth Marine Laboratory, Prospect Place, The Hoe, Plymouth,  PL1 3DH, UK}



\correspondence{L.~de Mora (ledm@pml.ac.uk)}

\runningtitle{The assessment of a global marine ecosystem model using emergent properties}

\runningauthor{L.~de Mora et~al.}



\received{26~June~2015}
\accepted{20~July~2015}
\published{}


\firstpage{1}

\maketitle


\begin{abstract}
      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.
\end{abstract}


\hack{\newpage}


\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 which demonstrate model validation
      using static fields, spatial distributions and dynamic variability,
      including
      \citet{Droop1973,Fasham1990,Taylor2001,Blackford2004,Allen2007a,Jolliff2009,Shutler2011,SauxPicart2012,DeMora2013,Kwiatkowski2014}.
      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 \citet{DeMora2013} is vulnerable to discrepancies between
      the model and the observations in the location of important
      circulation features such as fronts, coastlines or up-welling 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
      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{\degree} global
      model. The disentanglement of the performance of the biogeochemical
      model from that of the physics is a~major challenge in marine
      ecosystem modelling \citep{Holt2014}.

      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
      datasets.

      As a~solution to the problem of the absence of data and presence of
      poorly constrained physics, \citet{Holt2014} 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 
      had also been observed in nature in \citet{Hirata2011}.
      In effect, a~relationship seen in in situ observations
      also appeared in multiple biogeochemical models. 
      In addition, the relationship observed in \citet{Holt2014} was a widespread general response
      in all the plankton models that was independent of local hydrodynamic conditions.
      Furthermore, this relationship is a large scale 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 \citet{Hirata2011}, there are many works that link a large phytoplankton
      size class with the community structure: \citet{ Devred2011,Brewin2012a,Brewin2014,Brewin2015a}.
      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~models framework.
      Ecosystem function in models can be both explicitly enforced during
      the model development and tuning, or they can occur without being
      explicitly parameterized.  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 \citet{Holt2014}, 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 \citet{Holt2014}, 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.

%       In addition, model validation with emergent properties can be more
%       valuable as a model test than a~direct comparison of model to laboratory-based experiments,
%       such as primary production bottle measurements. This is because emergent
%       property validation can be a~large scale test of many combinations of
%       factors, simultaneously testing the model in a~wide range of physical
%       and biogeochemical environments, where as laboratory experiments can
%       not simultaneously cover such a~diverse range of settings.

      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 the
      European Regional Seas Ecosystem Model (ERSEM), and it is
      run in a coupled NEMO-ERSEM 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.~\ref{sec:NEMOERSEM} contains
      a~brief description of the circulation model, NEMO, the ecosystem model,
      ERSEM, and the sea ice model, CICE, used in this study.
      Section~\ref{sec:EP} 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 \citet{Holt2014} is included in
      Sect.~\ref{sec:CS}. Section~\ref{sec:C:Chl} shows the ratio of carbon
      to chlorophyll in phytoplankton and particulate organic matter.
      Section~\ref{sec:POCPON} demonstrates how the model reproduces the
      ratio of particulate organic carbon and nitrogen as described by
      \citet{Redfield1934,Martiny2013}.  Section~\ref{sec:Sto} illustrates
      the internal stoichiometric relationships for ERSEM for each of the
      nutrient currencies modelled.  Finally, Sect.~\ref{sec:discussion}
      discusses the successes, potential and limitations of these methods.



\section{The ERSEM and NEMO models}%s2
\label{sec:NEMOERSEM}

      The European Regional Seas Ecosystem Model (ERSEM) is a~lower trophic
      level biogeochemical cycling, carbon based process model that uses the
      functional-group approach \citep{Baretta1995,Blackford2004,Butenschon2015}. 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. 
%       added:
      A~diagram showing the major organisms, nutrients, chemical systems,
      organic matter and fluxes modelled in ERSEM is shown in Fig.~\ref{fig:ERSEM}.
      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 remineralization 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 \citep{Garcia2010}. The
      initial iron concentrations were zonally averaged interpolations of
      the iron data from \citet{Tagliabue2012a}.  The iron dust surface
      deposition climatology was based on \citet{Mahowald2005}.  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 Nucleus for European Modelling of the
      Ocean, NEMO \citep{Madec2008}.  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 time scales. The version
      of NEMO used in this study was version 3.2 and the ocean circulation
      model was interfaced with the Los Alamos Sea Ice model, CICE
      \citep{CICE2013}.  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{\degree} 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
      dataset \citep{CORE2}.  The coupled model was run for
      117\,\unit{years}, from 1890 to 2007.  The initial conditions for the
      circulation model at 1890 were taken from a~60\,\unit{year} physics-only
      climatological spin up from a~still global ocean.  The first
      60\,\unit{years} (1890--1950) of coupled NEMO-ERSEM running were spun
      up using climatological CORE2 forcing. After 1950, the remaining
      57\,\unit{years} were run with the inter-annually 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
      \citep{Kwiatkowski2014}.  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 run under identical physical conditions, using the same
%       parameters and settings for NEMO, the same computing resource, the
%       same coding framework and the same initial conditions for
%       nutrients. 
      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.



% \section{Ecosystem functions in marine biogeochemical models}%s3

\section{Emergent relationships in marine biogeochemical models}%s3
\label{sec:EP}

      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. 
      Firstly, 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. 
      Secondly, 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. 
      
%       Furthermore, many different combinations of model choices and parameters may result in 
%       similar emergent properties.
%       In addition, a model may reproduce an emergent while failing to reproduce others.
      
      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 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 datasets, 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 can not 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 \citet{DeMora2013}.

      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,
      \citep{Redfield1934,Arrigo2005}.  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~causal relationship, then it may be possible to
%       use the model to discover a~causal relationship.
      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
%       a~strong indication that the model has a~appropriate representation of
%       the ecosystem functions that create the emergent 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
      \citet{Holt2014}.




\subsection{Phytoplankton community structure}%s3.1
\label{sec:CS}

      In modelling and observational marine biology, phytoplankton
      are often grouped into size or function based classifications \citep{WOODWARD2005}.
      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  
      \citep{Nelson1993,Stolte1995,Huete-Ortega2012,Zhang2012}.
      At the 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 
      \citep{Riegman1993,Finkel2010,Huete-Ortega2012}.
      In addition to size classes, phytoplankton function
      based classifications are also used in modelling and observational marine biology.
      Phytoplankton functional types allow 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 function based
      classification.  The phytoplankton classification in ERSEM follows
      both size and function based classifications. There are four plankton
      functional types (PFT) 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 \citet{Hirata2011,
      Devred2011,Brewin2012a,Brewin2014,Brewin2015a}.  Each of these
      fits is prepared using a~different dataset.  \citet{Hirata2011} used
      multiple HPLC databases from around the world.  \citet{Devred2011}
      used chlorophyll and absorption data from the Northwest Atlantic
      region that was collected between 1996 and 2003.  \citet{Brewin2012a}
      used HPLC data in the Indian ocean taken between 1995 and 2007.
      \citet{Brewin2014} used SFF from the Atlantic Meridianal Transect
      (AMT) cruises in the Atlantic ocean between 1996 and
      2012. \citet{Brewin2015a} used a~aggregation of 16 unique
      globally distributed databases.

%       Figure~\ref{fig:CSfits} shows the five fits of in situ community
%       structure and the least squared fit of ERSEM to the three-population
%       absorption model of \citet{Brewin2010}. 
      Figure~\ref{fig:CSfits} shows the five fits of in situ community
      structure \citep{Hirata2011,Devred2011,Brewin2012a,Brewin2014,Brewin2015a}. 
      The fit of the ERSEM simulation to the three-population absorption model
      of \citet{Brewin2010} 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.  The middle panel shows nanophytoplankton and the
      right panel shows picophytoplankton.  The dashed vertical line
      indicates a~typical detection limit of HPLC and SFF measurements of
      0.1\,\unit{mg\,Chl\,m^{-3}}.  The
      \cite{Devred2011,Brewin2012a,Brewin2014,Brewin2015a} are also
      fitted to the three-population absorption model, but the
      \cite{Hirata2011} 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 datasets and methodologies.

      Note that only the fits are shown in Fig.~\ref{fig:CSfits}, the in
      situ data itself is 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
      \citet{Hirata2011} has already been smoothed with a~5-point running
      mean, and that running mean is further smoothed to the fit line shown
      in Fig.~\ref{fig:CSfits}.

      Figure~\ref{fig:CS} 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
      \citet{Brewin2015a,Hirata2011}.  The fit to ERSEM, the
      Brewin 2015 and the Hirata 2011 lines are identical in
      Figs.~\ref{fig:CSfits} and~\ref{fig:CS}. For clarity, the other fits
      are not shown in Fig.~\ref{fig:CS}.  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 is 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 is not shown. 

      The model data in Fig.~\ref{fig:CS} was taken from the top
      40\,\unit{m} of the global ocean, excluding the shallow seas, inland seas,
      the Southern Ocean from 45\degree\ South 
      and the Arctic Ocean from 55\degree\ North.
      The Arctic and Southern
      oceans were excluded from this analysis for two reasons: firstly,
      there were no Arctic community structure data in four of the studies
      mentioned above, and \citet{Brewin2015a} only includes a~small number
      of Southern Ocean data. Secondly, they were excluded because this is
      a~region where the physics in this model run was particularly
      troublesome.  As mentioned in Sect.~\ref{sec:NEMOERSEM}, 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{\degree} global model
      can not 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 \citep{Brewin2010},
      and both sets of fit parameters are shown in Table~\ref{tab:CS}.

      The model reproduces the overall shape of the picophytoplankton
      community structure in Fig.~\ref{fig:CS}, in that the
      picophytoplankton contribution to community chlorophyll is higher at
      low chlorophyll concentrations.  However, the picophytoplankton
      community structures is concave in the model and convex in the
      three-population model.  At low chlorophyll concentrations,
      the model picophytoplankton community
      structure is more similar in shape to the \citet{Hirata2011} fit,
      which is the only in situ data fit from Fig.~\ref{fig:CSfits} not
      based on the three population model.

      The three-population absorption model of \citet{Brewin2010} can be
      summarised in the following equations.  The picophytoplankton
      component of the total community chlorophyll, $\text{chl}_\mathrm{p}$, is
      calculated as:
\begin{align}
&\text{chl}_\mathrm{p}=C^{\mathrm{m}}_\mathrm{p} \cdot (1 - \exp(- S_\mathrm{p} \cdot\,\text{chl}))
\label{eq:chlp}
\end{align}%e1
      where $\text{chl}$ is the total community chlorophyll,
      $C^{\mathrm{m}}_\mathrm{p}$ is the maximum picophytoplankton chlorophyll, and
      $S_\mathrm{p}$ is the initial slope of the exponential function for
      picophytoplankton.  Note that Figs.~\ref{fig:CSfits} and~\ref{fig:CS}
      show the phytoplankton functional type contribution as a~percentage of
      total community chlorophyll instead of as the chlorophyll
      concentration used in Eq.~(\ref{eq:chlp}).

      The combined picophytoplankton and nanophytoplankton chlorophyll, or
      ``piconano'' chlorophyll, $\text{chl}_{\text{p, n}}$, is calculated in
      Eq.~(\ref{eq:chlpn}).
\begin{align}
&\text{chl}_{\text{p, n}}=C^{\mathrm{m}}_{\text{p, n}} \cdot (1 - \exp (-S_{\text{p, n}} \cdot\,\text{chl}))
   \label{eq:chlpn}
\end{align}%e2
      where $C^{\mathrm{m}}_{\text{p, n}}$ is the maximum piconano chlorophyll and
      $S_{\text{p, n}}$ is the initial slope of the exponential function for
      piconano.  The $\text{chl}_{\text{p, n}}$ function is not shown
      explicitly in either Fig.~\ref{fig:CSfits} or Fig.~\ref{fig:CS}, but
      is used in the calculation of the nanophytoplankton, diatoms and large
      phytoplankton functional types community structure.

      Unlike $\text{chl}_\mathrm{p}$ and $\text{chl}_{\text{p, n}}$, the
      nanophytoplankton chlorophyll, $\text{chl}_\mathrm{n}$, and
      microphytoplankton chlorophyll, $\text{chl}_\mathrm{m}$, are not explicitly
      fitted. Instead, their contribution to total community chlorophyll are
      determined from a~combination of Eqs.~(\ref{eq:chlp})
      and (\ref{eq:chlpn}).  The nanophytoplankton chlorophyll,
      $\text{chl}_\mathrm{n}$, shown in Eq.~(\ref{eq:chln}), is the difference
      between the piconano group chlorophyll and the picophytoplankton
      chlorophyll.  This is shown in the middle pane of
      Figs.~\ref{fig:CSfits} and~\ref{fig:CS}.
\begin{align}
&\text{chl}_\mathrm{n} = \text{chl}_{\text{p, n}} - \text{chl}_\mathrm{p}
   \label{eq:chln}
\end{align}%e3

      The diatoms and large phytoplankton, (together also known as
      microphytoplankton) chlorophyll, $\text{chl}_\mathrm{m}$, are the remainder
      of the total chlorophyll that is not accounted for by the piconano
      component.  It is calculated is the difference between the piconano
      chlorophyll and the total chlorophyll. This is shown in the top pane
      of Fig.~\ref{fig:CS}.
\begin{align}
&\text{chl}_\mathrm{m}  = \text{chl}    - \text{chl}_{\text{p, n}}
   \label{eq:chlm}
\end{align}%e4

      Just to explicitly state the overarching assumption used here, the
      total chlorophyll, $\text{chl}$, is equal to the sum of the three components
      functional types chlorophyll.
\begin{align}
 &\text{chl} =  \text{chl}_\mathrm{m} + \text{chl}_\mathrm{n} + \text{chl}_\mathrm{p}
   \label{eq:chlt}
\end{align}%e5

      Note that the three fits of Fig.~\ref{fig:CS} 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.~(\ref{eq:chlp}),
      (\ref{eq:chln}) and (\ref{eq:chlm}). Equations~(\ref{eq:chlp})
      and (\ref{eq:chlpn}) are fitted to the simulated ERSEM data using in the
      least squares fit. 
%       In addition, the fits are most influenced by the
%       high data density regions in this figure. 
      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.

      




\subsection{The carbon to chlorophyll ratio in particulate organic matter and phytoplankton}%s3.2
\label{sec:C:Chl}

      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 \citep{Geider1987} and is affected by
      the cells response to irradiation, temperature and nutrient supply.
%       On the ecosystem scale, the carbon to chlorophyll ratio is influenced
%       by the food web, community structure and stoichiometric balance. 
      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  \citet{Sathyendranath2009,Geider1997}.
      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 \citep{Geider1987}.  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,
      \citep{MacIntyre2002,Polimene2014}. This behaviour is known as photo-acclimation.
      Similarly, it has long been known that the carbon to chlorophyll ratio
      increases with decreasing temperature \citep{Eppley1972}.  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 \citet{Geider1987}.  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 other particulate organic  matter \citet{Graff2015}.
      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 \citet{Sathyendranath2009}.
      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:
\begin{align}
&\label{eq:POCChl}
\text{POC} = m \cdot\,\text{chl}^{p}
\end{align}%e6
      where $\text{POC}$ is the particulate organic carbon, $\text{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, \citet{Sathyendranath2009} 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.~(\ref{eq:POCChl}):
\begin{align}
&\label{eq:CChl}
C = n \cdot\,\text{chl}^{q}
\end{align}%e7
      where $C$ is the phytoplankton carbon, $\text{chl}$ is the total
      community chlorophyll and $n$ and $q$ are the fitted parameters.

      ERSEM employs a~variant of the Geider model
      \citep{Geider1997,Geider1998} 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 a~direct consequence of
%       model choices, and can not be tuned with a~small number of
%       parameters. 
      For these reasons, the ratio of POC to community chlorophyll
      is not an obvious and purposefully prescribed consequence of model choices, 
      and can not be tuned with a~small number of parameters.       
      The ratio of the particulate organic carbon against the
      total chlorophyll is not explicitly restrained or parameterized in any
      way.

      Figure~\ref{fig:POCChl} shows the relationship between particulate
      organic carbon and total chlorophyll.  The model data is 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\,\unit{m} of the monthly climatology of the final ten years of the
      simulation, excluding shallow seas, inland seas, the Southern Ocean from 45\degree\ South,
      and the Arctic Ocean from 55\degree\ North.  A~fit of the models relationship between
      particulate organic carbon and total chlorophyll to
      Eq.~(\ref{eq:POCChl}) is also shown as a~full green line.  The two
      \citet{Sathyendranath2009} fits of POC to chlorophyll are also shown
      in 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~\ref{tab:POCChl}.

      The two 1\,{\%} quantile regressions fits are included in
      Fig.~\ref{fig:POCChl} 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\,{\%}
      \citet{Sathyendranath2009} 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\,\unit{m} deep,
      including the Arctic and Southern Oceans, and inland and shallow seas, the
      fraction of data becomes as high as 11\,{\%}.  
%       This is another
%       indication that either the model has not captured the behaviour of the
%       Polar regions, or the emergent property breaks down in these regions.
%       Unfortunately, \citet{Sathyendranath2009} did not include any data
%       from Polar regions in the winter that could be used to test this hypothesis.

      Figure~\ref{fig:PChl} 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.~(\ref{eq:CChl}).
      In a~modelling analysis, it is unnecessary to apply the 1\,{\%}
      quantile regression, as it is straightforward to extract phytoplankton
      carbon.  Figure~\ref{fig:PChl} also shows the two 1\,{\%} quantile
      regression fits from \citet{Sathyendranath2009} as dashed lines as in
      Fig.~\ref{fig:POCChl}.  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\,\unit{m} of the monthly
      climatology of the final ten years of the simulation, excluding
      shallow seas, inland seas and the Arctic and Southern Oceans.  Both
      the fits and the bulk of the model data distribution in this figure coincide
      with the two \citet{Sathyendranath2009} 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. 
      
%       In Fig.~\ref{fig:PChl}, there is a~region where the phytoplankton
%       carbon is much less than the \citet{Sathyendranath2009} fit, and this
%       group of data appears to have a~different slope.  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 dataset.  Similarly to
%       Sect.~\ref{sec:CS}, 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\,\unit{m} are included in the model, the number of
%       points included in these regions of the figure increases up to
%       approximately 10\,{\%} of the dataset.
% 
%       The carbon to chlorophyll relationship has many knock-on effects in
%       the model: it influences the entire carbon cycle and
%       has a~huge impact on the calculation of total global primary
%       production.  It also relates to the community structure shown in
%       Sect.~\ref{sec:CS}.  When the model successfully reproduces the carbon
%       to chlorophyll ratio in a~global ocean simulation, this is an
%       indicator that it has a~good 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 in a~range of environments.



\subsection{Particulate organic carbon and particulate organic nitrogen}%s3.3
\label{sec:POCPON}

      The ratio of carbon to nitrogen in the ocean has long been
      a~historically important measurement and was first identified
      by \citet{Redfield1934}. 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
      \citet{Martiny2013}.  Their paper used over forty thousands globally
      distributed POC-PON paired samples from 5383 unique stations in the
      upper 200\,\unit{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 dataset than to the mean or
      the mode.

      The POC\,:\,PON frequency distribution 
      from \citet{Martiny2013} was reproduced
      using simulated data in Fig.~\ref{fig:POCPON}.  This figure shows the
      POC\,:\,PON ratio calculated in the model, the POC\,:\,PON ratio of
      the in situ data from \citet{Martiny2013}, but also the canonical
      Redfield ratio.  Both the model and the in situ histograms in
      Fig.~\ref{fig:POCPON} were normalised to unity area.  A~summary of the
      statistical analysis of Fig.~\ref{fig:POCPON} is shown in
      Table~\ref{tab:Martiny}.  The mode shown in Table~\ref{tab:Martiny}
      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.~\ref{fig:POCPON}.

      The POC of the model in Fig.~\ref{fig:POCPON} 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.~\ref{sec:C:Chl}.  The PON was calculated as the sum of the
      nitrogen components of the same groups.  An artificial detection limit
      of 0.1\,\unit{{\mu}mol\,m^{-3}} was applied to the modelled PON
      component of the ratio.  The model data was a~monthly climatology of
      the final ten years of the simulation (1997--2007), excluding the
      Arctic ocean and all model data below 200\,\unit{m}.  Unlike the data
      selections of Sects.~\ref{sec:CS} and~\ref{sec:C:Chl}, the Antarctic
      domain, inland seas and shallow seas are included in this data.  This
      is because the \citet{Martiny2013} dataset 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, particulate organic matter (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.  
      %The three particulate detritus classes include carbon,
      %nitrogen and phosphorus currencies, and all three also have a~silicon
      %and/or an iron component. 
      
%       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.
% 
% %       In addition to those fields included in the POC and PON fields, ERSEM
% %       includes a~bacterial component.  The bacterial component of the model
% %       compete for nutrients with the phytoplankton, are predated by
% %       zooplankton, but do not contribute to particulate detritus.
% 
%       While the model captures the mode 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 in the tails below 5 or more than the Redfield
%       ratio.  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 \citet{Martiny2013} 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{\degree} by 1{\degree} 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{\degree} by 1{\degree} square of ocean.
%       Secondly, 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 standard error of a~distributions decreases with increasing number of data.  
% %       Although the POC\,:\,PON is
% %       clearly not a~Gaussian distribution, 
% %       The difference in sample size may partially explain the difference in standard deviation.
% %       While the Martiny dataset is composed of approximately 40k data,
% %       a~monthly climatology of a~1{\degree} model contains more than
% %       21\,million points even after removing land points, masked deep water
% %       and Arctic data points.
%       However, the range of stoichiometric
%       variability seen in measured POC\,:\,PON data is underestimated in the
%       model.
% 
%       Furthermore, there are precisely zero model data with a~POC\,:\,PON
%       ratio below 4.3 or above 16.5.  
%       The model  a~POC\,:\,PON distribution shows a tail on the right hand side of the 
%       distribution.
% %       a~better job of
% %       capturing the excess POC relative to PON than it does capturing the
% %       excessive PON. 
%       This behaviour 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.
% 
%       Despite these limitations, the ERSEM simulation was very successful at
%       reproducing the mode of the POC\,:\,PON ratio of \citet{Martiny2013}.
%       The means that a~the most common values in the model data are the same
%       as the most common values in the POC\,:\,PON ratio data.  The
%       reproduction of the mode of the dataset by the ERSEM model is a~strong
%       indication that the most common behaviour of the POC\,:\,PON
%       relationship is appropriately simulated.



\subsection{Intracellular elemental stoichiometry}%s3.4
\label{sec:Sto}

      Stoichiometry is the balance of each element in organisms and
      in the ecosystem \citet{Sterner2002}.  As mentioned previously,
      Redfield observed co-variability in the concentrations of dissolved
      nitrate and phosphate in seawater and in the composition of plankton,
      \citep{Redfield1934}.  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, \citep{Martiny2013}.
      The ratio of each element against carbon has been shown to vary
      significantly between region, taxa, ecosystem role and physiological
      status, \citep{Moore2013}.  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, co-limitation by
      two or more nutrients is still poorly understood, and appears
      infrequently in models, \citep{Moore2013}.  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 \citep{Droop1973}, 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, \citep{Johnson2011}.

      In their meta-study, \citet{Moore2013} compiled multiple datasets 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~\ref{fig:stoichio} shows a~version of this figure produced
      using the four ERSEM nutrients. This figure shows the model
      distribution with mean of the model data in circular markers, and the
      typical in situ value and observed range from \citet{Moore2013} in
      square markers with black bars.

      The organic component of Fig.~\ref{fig:stoichio} 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
      Sects.~\ref{sec:C:Chl} and~\ref{sec:POCPON}: 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 \citet{Moore2013}, 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 phytoplankton functional types, zooplankton
      functional types and detrital fields and their interactions need to be
      parameterized sensibly.  The overall particulate organic matter
      stoichiometry in ERSEM is not susceptible to tuning via a~small number
      of parameters.

%       Furthermore, the nutrient cycles of carbon, nitrogen, phosphorus and
%       iron are also influenced by the bacterial loop.  The bacterial biomass
%       does not contribute their biomass to the calculation of particulate
%       organic matter used here, but the bacterial functional type competes
%       with the phytoplankton for the inorganic nutrients except dissolved
%       iron and is an additional food source for zooplankton. This variable
%       stoichiometry bacterial functional type contributes to the flexibility
%       of ERSEM.  The bacteria also compete with the zooplankton to scavenge
%       non-silicon particulate detritus.  In addition, the bacterial
%       group only excretes to the dissolved organic matter detrital fields.

      The inorganic component of the model data was taken directly from the
      output of the simulation. The dissolved inorganic carbon (DIC) cycle
      in the model is described in \citet{Artioli2012}.  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 \citep{Aumont2003} 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.~\ref{fig:stoichio} are shown in Table~\ref{tab:Stoichio}.  This
      table includes the ``typical'' \citet{Moore2013} data mean and range,
      and the data from the model run.  \citet{Moore2013} did not include
      any measure of variability in the inorganic data; the inorganic range
      shown in Fig.~\ref{fig:stoichio} and Table~\ref{tab:Stoichio} 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 order 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 \citet{Moore2013}
      ``typical'' value as a 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 \citep{Garcia2010} for nitrogen,
      phosphorus and silicon, and GEOTRACES for Iron \citep{Henderson2007}.
      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.  The
      \citet{Moore2013} ``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, \citet{Moore2013} cite an
      observed maximum phytoplankton N\,:\,C quota of 0.169, 
      but the \citet{Martiny2013} dataset has observational data all the
      way down to their cut off point, which was a C\,:\,N ratio of 2\,:\,1.
      Note that \citet{Moore2013} used the  N\,:\,C ratio, but  \citet{Martiny2013} 
      used a C\,:\,N ratio. 
%       (Note that \citet{Martiny2013}
%       showed C\,:\,N instead of N\,:\,C, so their lower cut off value is
%       2.0 in Fig.~\ref{fig:stoichio}.)
      Also note that the model data does
      not show any histogram distribution, only presence/absence.  The
      presence/absence model data in Fig.~\ref{fig:stoichio} is extracted
      from the top 40\,\unit{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 was not further sub-sampled to match
      the sampling locations of the in situ data.  This data has monthly
      time resolution and is not climatologically averaged.

      In addition, it has been shown by \citet{Karl-awaitingReference} 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.
% 
%       The ratio of nitrogen to carbon is shown in blue in
%       Fig.~\ref{fig:stoichio}.  The model captured the mean organic ratio,
%       but had a~wider range of values than that quoted in the
%       \citet{Moore2013} work.  However, the maximum
%       value of the  N\,:\,C ratio has been extended from 0.169 to 0.5 after 
%       the \citet{Martiny2013} results have been 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 WOA dataset.  
% %       The inorganic nitrogen in the model never gets as
% %       depleted as the estimate of the minimum observed ratio.
% %       This is a~problem with the model
% %       parameterisations that has also been seen in Sect.~\ref{sec:POCPON}
% %       which will need to be addressed in future parameterisations.
% 
%       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 \citet{Moore2013} 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 lower limit of the dissolved
%       inorganic phosphorus to DIC ratio in the model is higher than
%       the lower limit of the same ratio observed in the estimate of the minimum observed ratio.       
%       The inorganic phosphorus in the model never gets
%       as depleted as the minimum observed ratio estimate. 
% %       However, the inorganic range shown
% %       in figure has no indication of the rarity or locations of such low
% %       P\,:\,C ratios.
%       However, the  range of the inorganic P\,:\,C ratio from the WOA estimates % shown in figure~\ref{fig:stoichio} 
%       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 extremely rare.
%       In addition, the model is the mean of a~1{\degree} by
%       1{\degree} patch of ocean, whereas the data originates from the mean of
%       a~one litre bottle, which would imply that we can expect fewer extreme values in the model.
%       The modelled P\,:\,C also illustrates the effect of
%       external resource limitation: as the inorganic  P\,:\,C ratio decreases, the
%       organic P\,:\,C ratio simultaneously decreases.
%       This figure does not 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.
% %       However, we can infer that a decrease in 
% %       However, it is likely that all three groups would see some of the decrease.
% %       modelled phytoplankton become nutrient stressed in low phosphorus environments, 
% %       which has a knock on effect of the total particulate organic matter.
% 
%       
%       The ERSEM model has four pelagic silicon fields: diatom silicon,
%       inorganic silicate, and medium and large detritus silicon.  The
%       C\,:\,Si ratio in \citet{Moore2013} and in Fig.~\ref{fig:stoichio} are
%       strictly limited to diatoms; there are no quotas associated with
%       silicon in other components of the ecosystem. The modelled ratio of
%       silicon to carbon, shown in purple in Fig.~\ref{fig:stoichio},
%       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, \citet{Garcia2010}, and 
%       is shown as a vertical purple dashed line in Fig.~\ref{fig:stoichio}.
%       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.  Also,
%       there is only a~very slim range of Si\,:\,C ratios allowed in the
%       modelled diatoms. 
%       This means that the nutrient stress effect seen in the phosphorus:carbon
%       ratio is not seen in the silicon:carbon ratio. 
%       Instead of becoming silicon-stressed at low inorganic Si\,:\,C ratios, 
%       the community structure changes to disfavour diatoms.
% 
%       The mean organic iron:carbon ratio in the model is lower than the 
%       same ratio in \citet{Moore2013}: 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 arguably 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 is ERSEM is
%       much less complex than that seen in a~nature. An example of a~more
%       natural iron model is \citet{Tagliabue2009}, which has three
%       bio-available forms of iron and two complexed forms of iron.  Despite
%       this, as the inorganic Fe\,:\,C ratio decreases, the organic ratio
%       also decreases, indicating that the phytoplankton 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, \citep{Paerl1997}.  Unfortunately, this model does
%       not include any anthropogenic nutrient loading, or indeed any flux of
%       riverine nutrients, and the dust deposition is forced with an annual
%       climatology.
% 
%       Overall, this figure informs us about the relationship between the
%       inorganic and organic component of the stoichiometric balance.
%       Effectively, it illustrates whether nutrient limitation and nutrient
%       stress are parameterized 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 datasets as the ocean is vastly under-sampled relative to the
%       model.  On the other hand, the model is the mean of a~1{\degree} by
%       1{\degree} patch of ocean, and the data is typically the mean of
%       a~one litre bottle, which would imply less variability in the model.
%       Furthermore, some of the in situ data may originate in coastal
%       datasets which have a~higher spatial variability than would be seen in
%       a~coarse global model.



\section{Discussions}%s4
\label{sec:discussion}

      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.  


      
% % % % % % % % % % % % % % % % % % % % 
% % %  Section 4.1
\subsection{Phytoplankton community structure}%s4.1
\label{dis:4.1}
% 
      % Move this into discussion?
      In the top panel of Fig.~\ref{fig:CS}, 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 increased in the polar
      regions \citep{Sosik2002}, we postulate that it is a~combination of
      multiple factors that creates this excess microphytoplankton. Firstly,
      in the polar regions there is an abundance of nutrients, and
      especially silicon, caused by excessive mixing in the physical model.
      Secondly, the model is parameterized to favour diatoms in low light
      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 three population model was performed to all the model data
      from the top 40\,\unit{m} of the surface.  The results of this fit are
      shown in the ERSEM (Top 40\,\unit{m}) column of Table~\ref{tab:CS}.       
      Relative to the \citet{Brewin2015a} 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.~\ref{fig:CSfits} and \ref{fig:CS} 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 was limited to the top 40\,\unit{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  \citep{Brewin2015a}, 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
      \citet{Hirata2011} than like the convex three-population model of
      \citet{Brewin2010} at low chlorophyll concentrations. 
      
% However, the picophytoplankton
%       community structures is concave in the model and convex in the
%       three-population model.      
      
%       Fig.~\ref{fig:CS} shows that the
%       model diatoms and large phytoplankton can sometimes dominate
%       the community structure at low chlorophyll concentrations in this simulation.
%       The model data was limited to the top 40\,\unit{m} of the
%       surface ocean, and the relationship breaks down in Arctic waters in
%       the model, and 
%       this effect is seen more prominnently when the Arctic and Southern Ocean data is included.
%       Hints of the breakdown of the community structure appear in
%       the in situ data \citet{Brewin2015a}, 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.
%       However, this effect is likely to be mitigated over most of the ocean by lower silicate
%       concentrations leading to slightly stronger silicate limitation for
%       diatoms.     
%       In addition, the picophytoplankton chlorophyll as a~fraction of the
%       community behaves much more like the empirical fit to data seen in
%       \citet{Hirata2011} than like the three-population model of
%       \citet{Brewin2010}.      
      
      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. 
      
%       Furthermore, the ability of the model to reproduce this emergent
%       property did not require a spatial comparison of model and data.
      
%       was not linked with any spatial 
%       Overall, it was shown in Figs.~\ref{fig:CSfits} and~\ref{fig:CS} that the model captures the
%       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 appropriate to bring out a~natural emergent community
%       structure.  

      

      
% % % % % % % % % % % % % % % % % % % % 
% % %  Section 4.2
\subsection{The carbon to chlorophyll ratio}%s4.2
      

      Together, the particulate organic carbon to chlorophyll and
      phytoplankton carbon to chlorophyll ratios from Sect.~\ref{sec:C:Chl}
      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.~\ref{fig:POCChl} and~\ref{fig:PChl}, there is a~region where the phytoplankton
      carbon is much less than the \citet{Sathyendranath2009} fit, and this
      group of data appears 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 dataset.  Similarly to
      Sect.~\ref{sec:CS}, 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\,\unit{m} are included in the analysis, the number of
      points included in these regions increases up to
      approximately 11\,{\%} of the dataset.  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, \citet{Sathyendranath2009} did not include any data
      from Polar regions in the winter that could be used to test this hypothesis.

%       However, the data below this line occur in the model only
%       on the edges of the Arctic and Antarctic oceans in the winter.  When
%       the entire model domain is included down to 40\,\unit{m} deep,
%       including the arctic and Antarctic and inland and shallow seas, the
%       fraction of data becomes as high as 11\,{\%}.  This is another
%       indication that either the model has not captured the behaviour of the
%       Polar regions, or the emergent property breaks down in these regions.
%       Unfortunately, \citet{Sathyendranath2009} 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 knock-on 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.
%       or selective grazing by zooplankton. 
      
%       While the model has been shown to reproduce the emergent property,
%       it has not shown that the model parameterisations of ecosystem function are perfect. 
   
      
%       Secondly, 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.       
      
%       Together, the particulate organic carbon to chlorophyll and
%       phytoplankton carbon to chlorophyll ratios from Sect.~\ref{sec:C:Chl}
%       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. 
% %       This figure also helped
% %       identify some unusual behaviour from the diatom and large
% %       phytoplankton functional types. 
%       As was the case for the community structure, %Sect.~\ref{sec:CS}, the
%       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 were suggested as causes of this discrepancy.     
      
      
% % % % % % % % % % % % % % % % % % % % 
% % %  Section 4.3    
\subsection{Particulate organic carbon and particulate organic nitrogen}%s4.3
\label{dis:POCPON}

      Figure~\ref{fig:POCPON} 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 \citet{Martiny2013} 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{\degree} by 1{\degree} 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{\degree} by 1{\degree} square of ocean.
      Secondly, 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 is 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 \citet{Martiny2013} data ranges from 2. 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 validators 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 \citet{Martiny2013}
      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 dataset 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 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.      
      
      
      
 % % % % % % % % % % % % % % % % % % % % 
% % %  Section 4.4     
\subsection{Intracellular elemental stoichiometry}%s4.4


      The stoichiometric variability of particulate organic matter against
      the ratio of inorganic nutrients to DIC was shown in Fig.~\ref{fig:stoichio}.
      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 at low as the estimates of this ratio, except iron.
      This is a~problem with the model
      parameterisations that has also been seen in Sect.~\ref{dis:POCPON}
      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.~\ref{fig:stoichio}.  The model captured the mean organic N\,:\,C ratio,
      but had a~wider range of values than that quoted by
      \citet{Moore2013}.  However, the maximum
      value of the  N\,:\,C ratio has been extended from 0.169 to 0.5 after 
      the \citet{Martiny2013} results have been 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 WOA dataset.  

      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 \citet{Moore2013} 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.       
%       the lower limit of the same ratio observed in the estimate of the minimum observed ratio.           
      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 is the mean of a~1{\degree} by
      1{\degree} patch of ocean, whereas the observational data originates from the mean of
      a~one 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 can not 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.~\ref{dis:POCPON},
      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 \citet{Moore2013} and in Fig.~\ref{fig:stoichio} 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.~\ref{fig:stoichio},
      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, \citet{Garcia2010}, and 
      is shown as a vertical purple dashed line in Fig.~\ref{fig:stoichio}.
      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  
      \citep{Butenschon2015}. 
      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 \citet{Moore2013}: 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 is ERSEM is
      much less complex than that seen in a~nature. An example of a~more
      natural iron model is \citet{Tagliabue2009}, which has three
      bio-available forms of iron and two complexed forms of iron.  Despite
      this, as the inorganic Fe\,:\,C ratio in our model decreases, 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, \citep{Paerl1997}.  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.~\ref{fig:stoichio} informs about the relationship between the
      inorganic and organic component of the stoichiometric balance.
      Effectively, it illustrates whether nutrient limitation and nutrient
      stress are parameterized 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 datasets as the ocean is vastly under-sampled relative to the
      model.  On the other hand, the model is the mean of a~1{\degree} by
      1{\degree} patch of ocean, and the data is typically the mean of
      a~one litre bottle, which would imply less variability in the model.
      Furthermore, some of the in situ data may originate in coastal
      datasets which have a~higher spatial variability than would be seen in
      a~coarse global model. 
      
      
      
    
      
      
% % % % % % % % % % % % % % % % % 
%    Summary of all them.
\subsection{The role of emergent properties}%s4.5
     
      Combined together, these relationships have informed 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 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 parameterized 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 \citet{Vichi2007}
      by demonstrating  using emergent properties that many
      of ERSEMs design choices and parameterisation are still appropriate in
      a~global context.  Furthermore, these emergent relationships were not
      explicitly parameterized in model development; all of them arise
      naturally out of a~combination of many other ecosystem functions.
%       relatively understood      natural behaviours.

      The combination of these well known phenomena have allowed 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 inform 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. % or bacteria growth efficiency.


      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 sections~\ref{sec:POCPON} or \ref{sec:Sto},
      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.
%       
%       In order to get a~better validation of the ecosystem functions of the
%       bacterial class, this would require more measurements of the ratio of
%       primary production to bacterial production.  Alternatively, the
%       bacterial growth efficiency could be used to gauge the model's
%       bacterial behaviour.      
%       This variable stoichiometry bacterial functional type contributes to the flexibility
%       of ERSEM, but is not currently verified by any emergent property.      
 
    
      

      The models grazers biomass were implicitly included in the
      POC\,:\,PON, POC\,:\,Chl and the stoichiometric
      relationships. However, there are no metric included here to study the
      zooplankton by themselves.  The authors are not aware of any stable global
      emergent property for describing the grazer community.
%       structure of grazers.

      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~point to 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.

%       Added:

      
      
% For future work, it would be good to see this sort of evaluation conducted in the context
% of testing hypotheses around bgc model structure. For instance, "X has suggested that
% bgc models need to include process Y to properly capture nutrient dynamics. If this is
% true, we would expect a model that did not include process Y to exhibit behaviour Z
% (e.g. systemically underestimating POC/DIN ratios)..."

\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.  They allow us to
      demonstrate how natural behaviour emerges from the model.  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.

      As ecosystem functions arise independently of local physical conditions of
      the ocean, they can be used to demonstrate model quality in the case
      when the physical features of a~sea are not co-located in the model
      and in nature.
      While these methods can not 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.
%       ANONUMYOUS 3 WANTS MORE DETAILS ABOUT THIS.

%       Most importantly, ecosystem functions are the only way to demonstrate
%       the models capacity to represent ecosystem function, as opposed to
%       quantitative metrics of absolute ecosystem state. 
      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's important to remember that each emergent property 
      that the model fails to reproduce is a new direction for future model development.

\section*{Code availability}

      The ERSEM model is available under the GNU Public License version 3
      (GPL3).  Interested parties are encouraged to register at the Shelf
      Seas Biogeochemistry home page, \url{www.shelfseasmodelling.org}, to
      register in order to download and use the model for their own
      purposes.

      The analysis toolkit used in this work is available from the author
      upon request under the terms of the Revised Berkeley Software
      Distribution (BSD) 3-clause license.



\begin{acknowledgements}
      This work is supported by the NERC National Capability in Modelling
      programme at Plymouth Marine Laboratory the Integrated Global
      Biogeochemical Modelling Network to support UK Earth System Research
      (iMarNet) project.

      We acknowledge use of the MONSooN system, a~collaborative facility
      supplied under the Joint Weather and Climate Research Programme, which
      is a~strategic partnership between the Met Office and the Natural
      Environment Research Council.

      Thanks to R. W. J. Brewin for early access to the most recent
      parameterisation of his three component community structure model.
\end{acknowledgements}



%bbl:
% \input{gmd-2015-135-discussions-typeset_source-version3.bbl}

\bibliographystyle{copernicus}
% %  \bibliographystyle{plainnat}
% %   \bibliographystyle{plain}
\bibliography{ldemoraPaper}


%T1
\begin{table}[t]
\caption{Parameters of the fits to the three population absorption model for
the \citet{Brewin2015a}, and for ERSEM. Two fits to ERSEM are shown:
the first is the fit to the ERSEM dataset excluding the polar, shallow and
inland seas, and the second includes all these regions to a depth of
40\,\unit{m}. The parameters are: $C^{\mathrm{m}}_{\text{p, n}}$: maximum
piconano chlorophyll, $S_{\text{p, n}}$ is the slope for piconano
chlorophyll, $C^{\mathrm{m}}_\mathrm{p}$ is the maximum picophytoplankton
chlorophyll and $S_\mathrm{p}$ is the slope for picophytoplankton.}
     \label{tab:CS}
%\scalebox{.7}[.7]{%
\begin{tabular}{lrrr}
\tophline
                              & \citet{Brewin2015a} & ERSEM & ERSEM (Top 40\,\unit{m})\\
\middlehline
$C^{\mathrm{m}}_{\text{p, n}}$ & 0.77                      & 0.345 & 0.742                   \\[4pt]
$S_{\text{p, n}}$             & 1.22                      & 2.147 & 0.799                   \\[4pt]
$C^{\mathrm{m}}_\mathrm{p}$    & 0.13                      & 0.061 & 0.289                   \\[4pt]
$S_\mathrm{p}$                & 6.16                      & 10.56 & 1.168                   \\
\bottomhline
\end{tabular}
%}
%\hack{\setlength\tabularwidth{0.7142\tabularwidth}}%
%\scalebox{.7}[.7]{%
%\belowtable{%
%\hack{\vspace*{2mm}}%        <-- aber NUR bei Verkleinerungsfaktor notwendig !!
%}
%}
\end{table}



%T2
\begin{table}[t]
\caption{The parameters of the particulate organic carbon to total chlorophyll and phytoplankton carbon to chlorophyll fits to the Eqs.~(\ref{eq:POCChl}) and~(\ref{eq:CChl}).}
\label{tab:POCChl}
%\scalebox{.7}[.7]{%
\begin{tabular}{lrrr}
\tophline
 & \multicolumn{2}{c} {Sathyendranath et~al.~(2009)}       & ERSEM \\
 & HPLC                                           & Turner & \\
\middlehline
\multicolumn{4}{c} {Particulate Organic Carbon}  \\
\cline{1-4}
$m$ & $180 \pm 2$      & $157 \pm 2$      & 145 \\
$p$ & $0.48 \pm 0.014$ & $0.45 \pm 0.013$ & 0.51\\
\cline{1-4}
\multicolumn{4}{c} {Phytoplankton Carbon}  \\
\cline{1-4}
$n$ & 79$^{\ast}$   & 64$^{\ast}$   & 66  \\
$q$ & 0.65$^{\ast}$ & 0.63$^{\ast}$ & 0.72\\
\bottomhline
\end{tabular}
%}
%\hack{\setlength\tabularwidth{0.7142\tabularwidth}}%
%\scalebox{.7}[.7]{%
\belowtable{%
%\hack{\vspace*{2mm}}%        <-- aber NUR bei Verkleinerungsfaktor notwendig !!
$^{\ast}$\,These values were calculated using a~1\,{\%} quantile regression.}
%}
\end{table}




%T3
\begin{table}[t]
\caption{Statistics to describe the POC\,:\,PON distribution as reported by \citet{Martiny2013}, compared to the result of this study and to the canonical Redfield Ratio.}
\label{tab:Martiny}
%\scalebox{.7}[.7]{%
\begin{tabular}{lrrr}
\tophline
                &\citet{Martiny2013} &Redfield &ERSEM    \\
\middlehline
      Mean      &7.06   &6.63  &5.93      \\
      Mode      &5.9    &&5.8      \\
      Median    &6.5    &&5.84      \\
      SD        &2.46$^{\ast}$   &&0.61  \\
\bottomhline
\end{tabular}
%}
%\hack{\setlength\tabularwidth{0.7142\tabularwidth}}%
%\scalebox{.7}[.7]{%
\belowtable{%
%\hack{\vspace*{2mm}}%        <-- aber NUR bei Verkleinerungsfaktor notwendig !!
$^{\ast}$\,This standard deviation was not included in the original publication; it was calculated based on their dataset.}
%}
\end{table}



%T4
\begin{table}[t]
\caption{Table showing the typical, minimum and maximum  organic and inorganic ratios against carbon for \citet{Moore2013} and for the ERSEM simulation.}
\label{tab:Stoichio}
%\scalebox{.7}[.7]{%
\begin{tabular}{llrr}
\tophline
         &&Observed   &ERSEM    \\
\middlehline
               &N\,:\,C     &0.13     &0.17     \\
     Mean      &P\,:\,C     &0.008     &0.02       \\
     Organic   &Si\,:\,C    &0.13      &0.14      \\
               &Fe\,:\,C    &6.1\,$\times10^{5}$     &5.4\,$\times10^{6}$     \\
\cline{1-4}
               &N\,:\,C    &[0.05, 0.17  (0.5$^{\ast}$)]     &[0.09,  0.21 ]      \\
       Range   &P\,:\,C    &[0.001, 0.026]   &[0.006, 0.039]      \\
       Organic &Si\,:\,C   &[0.08, 1.01]     &[0.13,  0.15]    \\
               &Fe\,:\,C   &[2.1\,$\times10^{6}$, 2.6\,$\times10^{4}$]      &[1.5\,$\times10^{6}$,3.5\,$\times10^{5}$]    \\
\cline{1-4}
\cline{1-4}
               &N\,:\,C     &0.013       &0.0008     \\
     Mean      &P\,:\,C     &0.00089     &0.0002     \\
     Inorganic &Si\,:\,C    &0.044       &0.0053     \\
               &Fe\,:\,C    &2.4\,$\times10^{7}$    &2.8\,$\times10^{7}$      \\
\cline{1-4}
               &N\,:\,C     &[8.9\,$\times10^{8}$, 0.02]$^{\dagger}$    &[5.0\,$\times10^{5}$, 0.012]    \\
       Range   &P\,:\,C     &[8.9\,$\times10^{8}$, 0.01]$^{\dagger}$    &[2.8\,$\times10^{6}$, 0.0008 ]     \\
   Inorganic   &Si\,:\,C    &[1.7\,$\times10^{7}$, 0.06]$^{\dagger}$    &[5.4\,$\times10^{6}$, 0.036] \\
               &Fe\,:\,C    &[8.7\,$\times10^{9}$, 2.6\,$\times10^{6}$]$^{\ddagger}$    &[6.8\,$\times10^{9}$, 2.4\,$\times10^{6}$ ]     \\
\bottomhline
\end{tabular}
%}
%\hack{\setlength\tabularwidth{0.7142\tabularwidth}}%
%\scalebox{.7}[.7]{%
\belowtable{%
%\hack{\vspace*{2mm}}%        <-- aber NUR bei Verkleinerungsfaktor notwendig !!
All Observation data is taken from \citet{Moore2013} except:
$^{\ast}$\,indicates values from \citet{Martiny2013} which were not taken
into account in \citet{Moore2013}, $^{\dagger}$\ indicates data calculated using World Ocean Atlas data \citep{Garcia2010},
and $^{\ddagger}$\ indicates data calculated using GEOTRACES data \citep{Henderson2007}.
}
\end{table}


\begin{figure}
\includegraphics[width=120mm]{ERSEM.pdf}
\caption{ 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.
	}
\label{fig:ERSEM}
\end{figure}

\begin{figure}
\includegraphics[width=120mm]{gmd-2015-135-discussions-f01.pdf}
\caption{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 \citep{Hirata2011,Devred2011,Brewin2012a,Brewin2014, Brewin2015a}.
      The dashed vertical line indicates a~typical detection limit of HPLC and SFF methods.}
\label{fig:CSfits}
\end{figure}



\begin{figure}
\includegraphics[width=90mm]{gmd-2015-135-discussions-f02.pdf}
 \caption{Phytoplankton community structure.
      The model data is 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 \citet{Brewin2010} is shown as a~full green line, and
      the fit of historic in situ data to the three-population absorption model from \citet{Brewin2015a} is shown in a~purple line.
      A~fit to data from \citet{Hirata2011} is shown in a~black line.
      The dashed vertical line indicates a~typical detection limited of HPLC and SFF methods.}
\label{fig:CS}
\end{figure}





\begin{figure}
\includegraphics[width=120mm]{gmd-2015-135-discussions-f03.pdf}
  \caption{The ratio of particulate organic carbon to total chlorophyll.
         The model data is shown as the logarithmically scaled two dimensional data histogram distribution in blue-scale.
         The full lines indicate the two \citet{Sathyendranath2009} fits to data, and a~fit of the model to Eq.~(\ref{eq:POCChl}).
         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.}
 \label{fig:POCChl}
\end{figure}




\begin{figure}
\includegraphics[width=120mm]{gmd-2015-135-discussions-f04.pdf}
 \caption{The ratio of phytoplankton carbon to total chlorophyll.
         The model data is shown as the logarithmically scaled two dimensional data histogram distribution in blue-scale.
         The dashed lines show the fit to data, a~fit of the model to Eq.~(\ref{eq:CChl}), and the results of the two 1\,{\%} quantile regression from \citet{Sathyendranath2009}.}
\label{fig:PChl}
\end{figure}




\begin{figure}
\includegraphics[width=85mm]{gmd-2015-135-discussions-f05.pdf}
  \caption{The ratio of Particulate Organic Carbon to Particulate Organic Nitrogen in the \citet{Martiny2013} in situ dataset and in the model.
      The Redfield ratio is also shown as a~red vertical line.
      The model data was taken from a~monthly climatology of the top 200\,\unit{m} of the final ten years of running, excluding the Arctic ocean.
      Both the ERSEM and the \citet{Martiny2013} histogram were normalised to unity area.}
\label{fig:POCPON}
\end{figure}





\begin{figure}
\includegraphics[width=90mm]{gmd-2015-135-discussions-f06.pdf}
\caption{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 \citet{Moore2013}
         are shown as square markers with horizontal bars.
         The vertical bars associated with the square markers were calculated from
         World Ocean Atlas data \citep{Garcia2010} and GEOTRACES \citep{Henderson2007}.}
\label{fig:stoichio}
\end{figure}

\end{document}\endinput
