Interactive comment on “ EMPOWER-1 . 0 : an Efficient Model of Planktonic ecOsystems WrittEn in R ”

This manuscripts explains the technical details of a simple NPZD model that runs in a two-layer vertical setup. The authors claim that using simple ecosystem models such the one described in the manuscript are still useful because one can run them very fast and then be able of evaluating how changes in equation formulation or parametrization affect the ecosystem dynamics. I can see the point of this argument and I somehow agree with it, although with some reservations. Personally I think that the dichotomy over "simple vs. complex" models is overstated and should not be a matter of too

elling testbed, EMPOWER-1.0, coded in the freely available language R. The testbed uses simple two-layer "slab" physics whereby a seasonally varying mixed layer which contains the planktonic marine ecosystem is positioned above a deep layer that contains only nutrient. As such, EMPOWER-1.0 provides a readily available and easy to use tool for evaluating model structure, formulations and parameterisation. The code is 10 transparent and modular such that modifications and changes to model formulation are easily implemented allowing users to investigate and familiarise themselves with the inner workings of their models. It can be used either for preliminary model testing to set the stage for further work, e.g., coupling the ecosystem model to 1-D or 3-D physics, or for undertaking front line research in its own right. EMPOWER-1.0 also serves as 15 an ideal teaching tool. In order to demonstrate the utility of EMPOWER-1.0, we carried out both a parameter tuning exercise and structural sensitivity analysis. Parameter tuning was demonstrated for four contrasting ocean sites, focusing on Station India in the North Atlantic (60 • N, 20 • W), highlighting both the utility of undertaking a planned sensitivity analysis for this purpose, yet also the subjectivity which nevertheless surrounds 20 the choice of which parameters to tune. Structural sensitivity tests were then performed comparing different equations for calculating daily depth-integrated photosynthesis, as well as mortality terms for both phytoplankton and zooplankton. Regarding the calculation of daily photosynthesis, for example, results indicated that the model was relatively insensitive to the choice of photosynthesis-irradiance curve, but markedly sensitive to

Introduction
Ecosystem models are ubiquitous in marine science today, used to study a range of compelling topics including ocean biogeochemistry and its response to changing climate, end-to-end links from physics to fish and associated trophic cascades, the impact of pollution on the formation of harmful algal blooms, etc. Models have become progressively elaborated in recent years, a consequence of both superior computing power and an expanding knowledge base from field studies and laboratory experiments. All manner of models have appeared in the published literature varying in terms of structure, equations and parameterisation. Anderson et al. (2014), for example, commented on the "enormous" diversity seen in chosen formulations for dissolved organic matter 10 (DOM) in the current generation of marine ecosystem models and asked whether reliable simulations can be expected given this diversity. This question applies not just to modelling DOM, but also to most processes and components considered in modern marine ecosystem modelling. A certain amount of variability among models is to be expected because of differing 15 objectives among modelling studies. A distinction can, for example, be made between models designed primarily for improving understanding of system dynamics, as opposed to those for out-and-out prediction (Anderson, 2010). Ultimately, however, much of the variability seen in model structure and equations is an outcome of personal choice on the part of the practitioner. Indeed, the art of modelling is in making de-In the early days of marine ecosystem modelling, it was necessary to resort to simple empirical approaches to deal with physics given the limited power of computers at the time. The so-called zero-dimensional "slab" models that came to the fore were the cornerstone of their discipline in the mid 20th century. Slab models have a simple physical structure consisting of two vertical layers. The depth of the upper (mixed) layer, 15 which can vary seasonally, was determined empirically from observations of vertical profiles of temperature or density. Containing the pelagic marine ecosystem, the upper layer was positioned above an essentially implicit (in that it is unchanging) bottom layer that contains a (typically fixed) nutrient concentration. Such slab models can be run quickly and straightforwardly, enabling both a multitude of runs and ease of analysing served them so well, allowing them to set up meaningful simulations from which they could so effectively draw conclusions and make progress in their field of study.
The need for preparation in terms of exploring sensitivity to ecosystem model formulations and parameterisation is no less in the modern era, indeed it is arguably greater given our deeper knowledge of the marine biota and a correspondingly larger multitude 10 of mathematical formulations to choose from. We propose that modellers can benefit from extensively "playing with" and testing their models and that the use of simple slab physics is an obvious choice in this regard, at least for ocean locations where the bulk of the biological activity occurs in the surface mixed layer. Experimentation of this kind may then be used to set the stage for the "serious" model runs that may follow, e.g. in 1-D or 3-D, although it is also entirely possible to undertake successful studies using only slab physics models. In addition, because they are straightforward to understand and do not require powerful computing resources to run, such simple models are ideal for use in teaching future generations of marine scientists about ecological structure and function. Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a fixed, rather than seasonally-varying, mixed layer depth. Applying the model to study the plankton of Fladen Ground and other regions in the northern North Sea, Steele demonstrated good agreement between the model and estimates of production from observations. Through work such as this, Steele emphasised that it is simplification that allows us to most easily address the controlling factors in marine ecosystems. One of 5 Steele's best-remembered findings, demonstrated again using simple models, is that the form of the zooplankton closure term has important consequences for ecosystem dynamics and export flux (Steele and Henderson, 1992). This finding remains relevant to modellers today and, indeed, we will examine model sensitivity to zooplankton mortality in Sect. 4.4. 10 It was Geoff Evans and John Parslow who would make the next major advance in the development of slab models with their "model of annual plankton cycles" (Evans and Parslow, 1985). Following Steele, they opted for an NPZ ecosystem embedded within the same two-layer framework with the marine ecosystem restricted to the upper layer and a fixed nutrient concentration in the lower. Evans and Parslow provided a more 15 complete representation of the interaction of the marine ecosystem with its physical environment by allowing the depth of the mixed layer to vary seasonally with direct impacts on the model state variables. As the mixed layer deepens, nutrients are entrained from below while phytoplankton density is diluted because their surface layer biomass is spread over a greater depth. Conversely, as the mixed layer shallows, the concen- 20 trations of nutrients and phytoplankton are unchanged although losses occur on a per unit area (m −2 ) basis. As many zooplankton can swim, Evans and Parslow assumed that they are able to avoid detrainment in a similar manner to the assumptions of prior models (e.g. Steele, 1958;Riley et al., 1949), as well as mixing, in which case their concentration increases as MLD decreases. 25 Evans and Parslow (1985) also took seasonal and daily irradiance forcing into consideration, in combination with depth integration of a non-linear P -I curve. As opposed to previous studies that had used observations, variation in light at the ocean surface was calculated from standard trigonometric/astronomical formulae (Brock, 1981) with transmission losses in the atmosphere as 70 % of cloud cover and photosynthetically active radiation (PAR) as 3/8 of total irradiance. Variation in light with time of day was assumed to be triangular (Steele, 1962), permitting analytic integration in time.
A notable contribution of Evans and Parslow's (1985) paper is the appendix which provides the equations required to construct a model subroutine to calculate daily depth- 5 integrated photosynthesis in a model layer as a function of noon irradiance (PAR entering the layer from above), day length, phytoplankton concentration, rate of light extinction (Beer's law) and parameters for maximum photosynthesis and initial slope that define the P -I curve.
In common with their predecessors, Evans and Parslow were interested in the fac- 10 tors controlling the initiation of the spring phytoplankton bloom, focussing on the role of vertical mixing. Bloom initiation, they concluded, required a low rate of primary production over winter, which is to be expected in the North Atlantic due to deep mixed layers at that time, and is also linked to coupling between phytoplankton and grazers. The simplicity of the slab model was key to their conclusions as articulated in their own 15 words: "It is worth emphasising the advantages of analysing simple models, and simplifying models until they can be analysed". The controls on phytoplankton dynamics in high-nutrient low-chlorophyll (HNLC) areas such as the Subarctic Pacific has remained a topical issue ever since, in large part because limitation by iron is also indicated (Martin et al., 1994;Coale et al., 1996), but the role of grazing and the link between 20 phytoplankton-zooplankton coupling and mixed layer depth remains firmly established as a key mechanism in these systems (Frost, 1987;Fasham, 1995;Chai et al., 2000;Smith and Lancelot, 2004). Perhaps the most famous slab modelling paper, published five years after Evans and Parslow (1985), is the study of nitrogen cycling in the Sargasso Sea by Fasham Introduction  1985), with seasonally varying mixed layer depth and irradiance forcing. The novel aspects of FDM90 were instead related to additional complexity of the ecosystem, expanding from a simple NPZ to explicitly separate new and regenerated production by including state variables for nitrate and ammonium (critical for calculating the f ratio; Eppley and Peterson, 1979), as well as having a simple microbial loop of dissolved 5 organic nitrogen and bacteria. Sinking detritus was also added as a state variable, facilitating the prediction of export flux. The success of this model was due to it being the first attempt to fully elucidate the processes involved in the recycling of nitrogen in the euphotic zone, as well as the complimentary roles of zooplankton and bacteria. The simplified physics of the model allowed it to be run on PCs of that era and Fasham 10 purportedly distributed the code on floppy disks, allowing other researchers to run the model on their PCs. The description of the marine ecosystem provided by FDM90 has largely served as the foundation for marine ecosystem modelling ever since. With the advent of increasing computer power, as well as increasing interest in the spatio-temporal behaviour of 15 plankton systems, most modelling studies are now undertaken in 1-D or 3-D physical frameworks. Nevertheless, many slab modelling studies have been published since FDM90 which follow the basic design described above, or slight modifications thereof (Table 1). A range of ecosystem models of varying complexity have been incorporated within slab physics and applied to contrasting sites throughout the world ocean. 20 The basic physical construction is similar in most cases consisting of a classic slab structure with a seasonal cycle of mixed layer depth specified from data and seasonal irradiance from standard trigonometric equations. Remarkably, Evans and Parslow's (1985) equations for calculating daily depth-integrated photosynthesis have prevailed and been used in most studies. A more sophisticated calculation method was devel-25 oped by Morel (1988Morel ( , 1991 and a simplified form of this (Anderson, 1993) is examined in Sect. 4.3. The models in Table 1 have been used for a diverse range of applications including studies of parameter optimisation (Spitz et al., 1998;Fennel et al., 2001;Schartau et al., 2001;Hemmings et al., 2004), parameter sensitivity analysis (Mitra, Mitra et al., 2007Mitra et al., , 2014, phytoplankton bloom dynamics (Findlay et al., 2006), nutrient cycling via organic and inorganic pathways (Llebot et al., 2010), primary production in HNLC systems (Kidston et al., 2013) and primary production and export flux in contrasting regions (Fasham, 1995;Onitsuka and Yanagi, 2005).
3 Model description 5 We demonstrate the use of EMPOWER-1.0 using a simple NPZD ecosystem model and forcing for four time series stations in the ocean. The code is readily adapted to incorporate other ecosystem models, including the relatively complex models of the modern era, and/or forcing for other ocean sites. 10 The model uses slab physics as per Evans and Parslow (1985), namely a seasonally varying surface mixed layer that contains the ecosystem positioned above a deep homogeneous layer containing unchanging nutrient and no plankton (Fig. 2). We have also included temperature dependencies for the physiological rates in the ecosystem model (see below).

68
• 25 E). These stations were chosen because of their contrasting environments, as illustrated by the differences in forcing variables: seasonally varying mixed layer depth (MLD), irradiance (I) and sea surface temperature (T ) (Fig. 3), as well as deep nitrate 20 (N 0 ; see below). Mixed layer depths were taken from World Ocean Atlas Locarnini et al., 2010). In common with most previous slab modelling studies, noon (peak daily) irradiance at the ocean surface for a given latitude as a function of time of year is calculated using standard trigonometric/astronomical equations. The effect of clouds on atmospheric transmission was calculated using the model of Reed  1977). The equations for irradiance forcing are not usually provided as part of published model descriptions but, for completeness, they are listed here in Appendix A. The bottom layer in most slab models is assumed to have a fixed concentration of nutrient, N 0 . There is in reality a gradient of nutrient with depth and this can be represented empirically in slab models using simple functions of nutrients vs. depth (Frost, 1987;Steele and Henderson, 1993;Fasham, 1995). We adopted this approach here for stations India and Biotrans, using simple linear relationships with depth:  15 The NPZD ecosystem model we have implemented in EMPOWER is presented in Fig. 4 with dissolved inorganic nitrogen (N; the sum of nitrate and ammonium), phytoplankton (P ), zooplankton (Z) and detritus (D) as state variables. It is a simplification of the marine ecosystem inspired by that of FDM90 (note that, because we focus here on Station India, the version of the FDM90 model applied to Station India, Fasham (1993) 20 provides the more pertinent foundation). Improved formulations are implemented for multiple-prey grazing, plankton mortality, regeneration and other detrital loss terms, as well as alterations to the parameterisation. The equations are described below; model parameterisation is described in Sect. 4.1. The phytoplankton equation is: where the terms are growth, grazing and non-grazing mortality (linear and quadratic), physical losses due to mixing across the bottom of the mixed layer, and dilution effects of entrainment. H(t) is mixed layer depth (m) at time t and H (t) denotes the rate of change of H when dH/dt is positive (dilution). As explained above, when dH/dt is negative the change in density due to detrainment of mass from the mixed layer is 5 exactly balanced by the increasing density due to decreases in volume, and therefore detrainment does not alter the concentration of remaining biomass. Variable µ P is the vertically-averaged temperature-dependent daily growth rate, defined as the product of a temperature-dependent maximum growth rate, µ max P (T ), and non-dimensional limitation terms for nutrients and light, L N (N) and L I (I(t, z)):

Ecosystem model description
Note that µ P is calculated on a daily basis averaging over the time of day (t) and depth (z). Temperature and nutrients are assumed to be uniformly distributed throughout the mixed layer, in which case µ P is: With the assumption of balanced growth, µ max P (T ) is equal to the equivalent maximum photosynthetic rate, V max P (T ). The temperature dependence of photosynthesis is from Eppley (1972): where k N is the half saturation constant. The calculation of L I is the most mathematically complicated aspect of characterising phytoplankton growth in this model as it takes into consideration both seasonal and 5 diurnal patterns of irradiance arriving at the ocean surface (I 0 ), attenuation of irradiance with depth and photosynthesis as a function of light intensity. Light is assumed to vary with depth according to Beer's law (I = I 0 exp(−k PAR z)) and photosynthesis calculated using a photosynthesis-irradiance (P -I) curve. The daily depth-average photosynthetic rate is calculated over the course of the day using an assumed daily variation 10 of light, from which the daily average is derived. The user of EMPOWER is provided with a choice between two photosynthesis-irradiance (P -I) curves, a Smith function (Eq. 7) and an exponential function (Eq. 8) (Fig. 5): 15 Integration with depth (inner integral of Eq. 4) can be calculated analytically for either of the two P -I curves; equations are provided in Appendix B. The default method of handling the diurnal variation in irradiance at the ocean surface (outer integral of Eq. 4) is to do a numeric integration. The user may choose between assuming either a sinusoidal (Platt et al., 1990) or triangular (Steele, 1962;Evans and Parlsow, 1985) 20 pattern of irradiance throughout each day, from sunrise to sunset and peaking at noon. Analytic depth integrals require a Beer's law attenuation of light within the water column characterised by a single attenuation coefficient, k PAR . The simplest assumption, provided as the first of two options in EMPOWER, is that k PAR is the sum of attenuation due to water and phytoplankton, parameters k w and k c respectively: Parameters k w and k c are often assigned values of 0.04 m −1 and 0.03 m 2 (mmol N) −1 respectively (e.g., FDM90); these values are used in EMPOWER.

5
The assumption of a single mixed layer value of k PAR is questionable because in reality the value of k PAR varies with depth as a result of the changing spectral properties of the irradiance field. Red light is mostly absorbed by water in the upper few meters while blue penetrates deepest, with relatively efficient absorption by chlorophyll at both wavelengths. Based on a complex treatment of submarine light (Morel, 1988), a piece-10 wise approach to light attenuation was developed by Anderson (1993) with different values, k PAR,i , with i = 1 for depth range 0-5 m, i = 2 for depth range 5-23 m and i = 3 for depths > 23 m (i = 1, 2, 3), in each case k PAR (i ) is related to pigment (chlorophyll) concentration, C:

15
This approach to light attenuation is provided as the default option for use in EM-POWER. The values of the polynomial coefficients are listed in Table 2. The diurnal variation in light at the ocean surface over the course of a day may be reasonably approximated by a sinusoidal function that is symmetric about noon irradiance (Platt, 1980). Further simplification is possible by use of a linear model, i.e., 20 triangular centred at noon (e.g. Steele, 1962;Evans and Parlsow, 1985) because this simplifies the time integration. It should be noted here that despite Evans and Parslow's (1985) claim that differences between the triangular and sinusoidal approximations are minimal if the area under the curve is the same, they did not make the "equivalent area" adjustment to their formula, nor is their statement generically true (i.e. it depends on the peak light intensity, the attenuation of light with depth and the nonlinear P -I relationship). In EMPOWER, the default method of handling the diurnal variation in irradiance at the ocean surface is to do a numeric integration. The user may choose between assuming either a sinusoidal (Platt et al., 1990) or triangular (Steele, 1962;Evans and Parlsow, 1985) pattern of irradiance throughout each day, from sunrise to sunset and peaking at noon. Undertaking a numerical time integral involves computa-5 tional cost and two empirical methods (Evans and Parslow, 1985;Anderson, 1993) have been published that provide analytic calculations (i.e pre-determined formulae) for daily depth-integrated photosynthesis in a water column. Both are provided as options for use in EMPOWER and have the advantage of faster run time. The first of the two EMPOWER options is the depth-averaged light-dependent calculation of growth of Evans and Parslow (1985) which assumes a triangular pattern of daily irradiance, Beer's law for light attenuation (Eq. 9) and a Smith function as the P -I curve (Eq. 7). It has been a popular choice in previous slab modelling studies (Table 1). The second option is from Anderson (1993), which was developed as an empirical approximation to the spectrally resolved model of light attenuation and photosynthesis of Morel (1988) 15 used in combination with the polynomial method of integrating daily photosynthesis of Platt et al. (1990). It assumes a sinusoidal pattern of irradiance through the day, a piecewise Beer's law light attenuation (Eq. 10) and an exponential function as the P -I curve (Eq. 8). Parameter α, the initial slope of the P -I curve, is also spectrally dependent. The method of Anderson (1993) calculates the variation of α with depth as 20 a function of chlorophyll in the water column. Daily photosynthesis is then calculated using a polynomial approximation. The methods for calculating daily depth-integrated photosynthesis of Evans and Parslow (1985) and Anderson (1993) are non-trivial and, for completeness, the equations are supplied in Appendix C.
Grazing by zooplankton is assumed to be on both phytoplankton and detritus. This 25 choice was made in part to illustrate how to implement ingestion on multiple prey types, as such functions are used for more complex models (e.g. when there are multiple phytoplankton size classes or functional types and/or omnivory by zooplankton  (Gentleman et al., 2003). For example, the multiple-prey grazing formula used in FDM90 and Fasham (1993) is classified as an active switching response (Gentleman et al., 2003) which can display anomalous behaviour such as sub-optimal feeding (i.e. ingestion rates decreasing when prey availability increases). We have therefore opted to improve upon Fasham's choice by using a different multiple-5 prey response, but one that is nevertheless commonplace in the literature. Specifically, we have adopted a passive switching response where density dependence of the prey preferences arises due to inherent differences in the single-prey responses (see Gentleman et al., 2003). This sigmoidal (or Holling Type 3) response characterised as ( Fig. 6): where the term in parentheses is the zooplankton specific ingestion rate. This Sigmoidal formulation implies that the single-prey response for both phytoplankton and detritus are each sigmoidal (Type 3). Parameter I max is the maximum specific graz- 15 ing rate, which is the same for both phytoplankton and detritus and equates to their single prey maximum ingestion rates. Although parameters ϕ P and ϕ D are often called preferences in the literature, the actual prey preferences associated with this response (i.e. relative amount in the diet as compared to the environment) are densitydependent, with the relative preference for phytoplankton to detritus is determined by half-saturation value in the literature, is actually an arbitrary parameter (i.e. this formulation is over-parameterized, see Gentleman et al., 2003) whose value determines the assumed single-prey half saturation constants based on choices for the ϕ parameters. The Sigmoidal response assumes an interference effect of alternative prey in that as detritus increases, ingestion of phytoplankton decreases (with the same interaction 5 for phytoplankton and ingestion of detritus). This interference effect is not so great as losing the benefit of generalism, i.e. total ingestion always increases for an increase in total prey density. The non-equal preferences reduce the interference effect for phytoplankton, i.e. the contours in the first panel of Fig. 6 are more vertical than for equal preferences. The corrollary effect is that the increased ingestion by consuming both 10 phytoplankton and detritus vs. just phytoplankton is reduced as compared to when prey have equal preferences.
Regarding non-phytoplankton grazing mortality, Fasham (1993) used a non-linear Michaelis-Menten saturating function although a linear mortality term is the usual choice (e.g., FDM90). We opted for the more flexible approach of using both linear 15 and nonlinear terms (Yool et al., 2011(Yool et al., , 2013a. The former may account for metabolic losses or natural mortality. The use of an additional nonlinear term represents densitydependent loss processes, notably mortality due to infection by viruses. The abundance of viruses is highly dependent on the density of potential host cells (e.g., Weinbauer, 2004) and, as reviewed by Danovaro et al. (2011), there is "compelling" evidence 20 that, at least in some instances, viruses are responsible for the demise of phytoplankton blooms based on observations of high proportions (10-50 %) of infected cells (e.g., Bratbak et al., 1993Bratbak et al., , 1996. A quadratic form was used for the nonlinear mortality term (e.g., Kawamiya et al., 1995;Oschlies and Schartau, 2005) and all phytoplankton nongrazing mortality losses were allocated to detritus. 25 The equation for rate of change of zooplankton density is: where the terms are growth, mortality (linear and quadratic) and losses due to mixing and changing MLD. Zooplankton growth can be described as the product of gross growth efficiency (GGE) and intake, where GGEs are typically between 0.2 and 0.3 (Straile, 1997). Gross growth efficiency is itself the product of absorption efficiency, β (more commonly, but incorrectly, known as assimilation efficiency; e.g. see Mayor 5 et al., 2011) and net production efficiency, k NZ . Splitting into these separate parameters permits three-way fractionation of intake between egestion (i.e. faecal pellet production, 1 − β), growth (β · k NZ = GGE; first term in Eq. 13) and excretion (β(1 − k NZ )). A variety of formulations exist in ecosystem models to describe zooplankton mortality and the appropriate functional form has been and continues to be a hotly debated topic (Steele and Henderson, 1992;Edwards and Yool, 2000;Mitra et al., 2014). Most common are the linear and quadratic terms, although some authors have chosen to employ other non-linear functions (e.g. Fasham, 1993 used a Michaelis-Menten relationship). As with phytoplankton, we used both linear and quadratic non-linear terms (Yool et al., 2011). The linear term represents density-independent natural mortality, 15 whereas the quadratic term is considered to be due to predation by carnivores (whose population tracks that of the zooplankton). The different sources of mortality result in different fates for these terms. Loss from natural mortality is allocated to modelled detritus, which implies a broader size-class of modeled particulates (and therefore higher sinking rates) than when just phytoplankton death contributes to this variable. 20 The fate of the predation-related mortality is less obvious because the metabolic activity of higher predators would result in ingested material being converted into dissolved nutrients as well as larger particulates (e.g. fecal pellets and death). Moreover, the higher predators may export material from the local region with migration. FDM90, along with a suite of follow-on models, therefore chose to allocate predation-related zooplankton mortality between nutrients (ammonium and DON, attributed to excretion by higher predators) and material that is immediately exported from the system (e.g. attributed to fast-sinking detritus generated by higher predators). Similarly, Steele and Henderson (1992) Salihoglu et al., 2008;Hinckley et al., 2009;Ye et al., 2012). We argue, however, that this is not necessarily realistic given that detrital particles related to higher-predators are larger and therefore even faster-sinking than that produced by the modelled plankton. We have therefore 5 here adopted to follow the sage approach of the model pioneers and assume that the predation-related mortality represented by our quadratic term is instantly exported and thereby entirely lost from the surface mixed layer of the model. As with phytoplankton, zooplankton are subject to changes in concentration via mixing and changes in MLD.
The equation for the rate of change of dissolved inorganic nitrogen (DIN) density is: DIN is taken up by phytoplankton (first term) and, via the food web, regenerated with terms 2 and 3 in Eq. (14) representing excretion by zooplankton and remineralisation of detritus respectively. The fourth term represents the net transport due to mixing (i.e. supply by the deep water and loss from the surface layer). The last term represents the 15 net effect of volume changes, i.e. increases in DIN density due to supply of deep water nutrients through entrainment and decreases in DIN density due to volume increases associated with entrainment. Finally, the detritus equation is: Detritus is produced by phytoplankton mortality, zooplankton natural mortality (linear term) and as zooplankton egestion (faecal pellet production). It is lost by zooplankton grazing and is also remineralised at a constant rate, m D . Detritus is mixed and subject 72 to changes via the seasonal cycle of MLD in the same manner as phytoplankton and zooplankton (terms 6 and 7), and also experiences losses due to gravitational sinking (last term). This occurs at rate v D (m d −1 ) and provides direct export of particulate organic matter to the layer below (where it is implicitly remineralised back to DIN). The first results Sects. 4.1 and 4.2 are devoted to parameterising the model for sta-5 tion India and a detailed description of values assigned to model parameters is provided therein.

Setup in R
We have chosen to code our model in the R programming language which can be readily downloaded for free over the internet. Input and output files are in ASCII text (.txt) format, avoiding the use of proprietary software. The structure of the code is designed to be transparent, where possible using conventional syntax common to different programming languages such as the use of loops, block IF statements, etc. As such, it can be relatively easily altered or translated into another programming language, if need be.
Where possible, we have followed what we consider to be best practice in developing iv. When a model run finishes, the summed annual fluxes associated with each term in the differential equations is displayed on the computer screen along with a re-5 port as to whether mass balance is achieved for each state variable (over the last year of simulation). Basic checking of mass balance is useful for ensuring that the model equations are error-free.
v. Regimented layout for clarity with extensive commenting throughout.
The R programming language is supported by various libraries that can be accessed 10 via the internet. One such library is for solving ordinary differential equations (Soetaert et al., 2010). Using this library has the advantage of minimising the length of the code and offers flexibility in terms of a range of numerical methods. On the other hand, its implementation requires that various conventions are adhered to and these can be restrictive when it comes to producing ancillary code, e.g., the formatting and export of 15 output files. As such, we opted to code the numerical solution of the ODEs manually within the core code of the model for several reasons: i. It offers full transparency for the interested user who wishes to see the method of integration.
ii. The use of manual code makes it considerably easier to export chosen variables 20 and fluxes to output files in desired formats and frequencies.
iii. In our case, the user is given the choice between two integration methods, Euler and fourth order Runge Kutta (RK4). These methods, particularly the latter, are entirely sufficient for the numerical task at hand and the coding of them is straightforward. iv. By using elementary syntax, the code can be easily altered or converted to other programming languages.
v. The code is stand alone and not subject to reformulation in the event of future changes in subroutine libraries.
The structure of the code is shown in Fig. 7. The functions come first, appearing 5 prior to the core code in R. The key function call is FNget_flux which contains the ecosystem model specification. The rate of change is calculated for each term in the differential equations and allocated to a 2-D array (flux no., state variable no.) which is then passed back to the core (permanent) code for processing. Other functions are: FNdaylcalc (calculates length of day), FNnoonparcalc (noon irradiance, PAR), FNLI-10 calcNum (undertakes numerical (over time) calculation of daily depth-integrated photosynthesis), FNLIcalcEP85 (calculates L I using the equations of Evans and Parslow, 1985), FNaphy (calculates chlorophyll absorption, effectively parameter α, in the water column after Anderson, 1993) and FNLIcalcA93 (calculates L I using the equations of Anderson, 1993). 15 Model setup comes next. Parameter values are read in from file NPZD_parms.txt. Simulation characteristics are then read in from file NPZDextra.txt. These include: i. Initial values for state variables.
ii. Run duration (years) and time step.
v. Choice of integration method: Euler or RK4.
vi. Choice of output characteristics: none, last year only or whole simulation, and a frequency of once per day or every time step. Model forcing for the chosen station of interest is then assigned. Monthly values of MLD and SST are read in and subject to linear interpolation in order to derive daily forcing. Other forcing variables are also set: latitude, deep nitrate (N 0 ; Eq. 1) and cloud fraction. At the end of the setup section there are a few lines of code that need to be altered if the ecosystem model is changed. These lines tell the computer how many 5 state variables the model has, the maximum number of flux terms associated with any one state variable and the maximum number of auxiliary variables to be stored for writing to output files. An advantage of this structure is that an initial section of customisable code is followed by a section of permanent code that does not require adjustment in the event of 10 changes to the equations that describe the ecosystem model, or indeed if a completely new ecosystem model is to be used. This code sets up a series of matrices to store fluxes and outputs and then integrates the model equations over time. State variables are updated and results exported to three output files: out_statevars.txt (state variables), out_aux.txt (chosen auxiliary variables) and out_fluxes.txt (all the terms in the 15 differential equations). These text files are readily imported to, for example, Microsoft Excel.
Results are plotted graphically on the computer screen at the completion of each simulation run. The graph plotting code is necessarily model specific and needs to be updated by the user as required. R is a user friendly programming language in this 20 regard and the code provided should be sufficient for the user to incorporate extra variables with ease.
Finally, a user guide is provided in Appendix D, outlining how to set up R, run the code, a summary of input and output files, and guidance on considerations when altering the ecosystem code and/or forcing.

Results
Model results are presented in four sections. First, a simulation is shown for station India using parameters taken from the literature (Sect. 4.1). Parameter tuning is then undertaken to fit all four ocean time series stations, India, Biotrans, Papa and Kerfix, to data for chlorophyll and nitrate at each site (Sect. 4.2). Moving on from the cali-5 bration of parameters, structural sensitivity analysis is then carried out by examining model sensitivity to equations for the calculation of daily depth-integrated photosynthesis (Sect. 4.3) and mortality of phytoplankton and zooplankton (Sect. 4.4).

Parameter initialisation: station India
Adjustment of parameters is a perennial problem for modellers. Parameters can be set 10 from the literature, sometimes directly on the basis of observation and experiment, but the usual starting point is to take values from previously published modelling studies. Almost inevitably, however, the resulting simulations will show mismatch with data and parameters are usually selected for adjustment (tuning) to improve the agreement with data. One option is to use objective tuning methods, such as the genetic algorithm or 15 adjoint method in which many or all of the model parameters are varied simultaneously in order to try and find a best fit solution to data (e.g., Friedrichs et al., 2007;Record et al., 2010;Ward et al., 2010;Xiao and Friedrichs, 2014). The advantage is objectivity, but difficulties include sloppy parameter sensitivities (parameters compensate for each other), different values of model parameters may be similarly consistent with the data 20 (the problem of identifiability), exploration of a huge parameter space may be required and local minima in misfit parameter space can make it difficult to find the true global minimum (Slezak et al., 2010). It is usually the case that models are underdetermined by data anyway ( Modellers more often than not carry out parameter adjustment by varying values of chosen parameters one at a time until satisfactory convergence with data is achieved. The skill is in deciding which parameters to vary. In principle, sensitivity analysis can be of help in this regard in that sensitive parameters can be identified and selected for adjustment if they can be justifiably altered (i.e., there is uncertainty regarding their 5 value). Here, we will demonstrate the use of EMPOWER for model calibration. Parameter sets will be derived for the four stations, India and Biotrans in the North Atlantic and the HNLC stations Papa (subarctic North Pacific) and Kerfix (Southern Ocean). The NPZD model we have presented is a new one and, as such, there is no readily available complete set of parameter values to draw upon. Using our experience, we chose 10 appropriate parameter values from the literature and adjusted others to give a good fit with the data for station India. This result is presented below along with a discussion of how we went about achieving this parameter set. Working from this parameter set, tuning of parameters is then undertaken to fit the other stations to the data.
Station India was previously modelled by Fasham (1993) and we used this publica- 15 tion as a starting point for the assignment of parameter values. Those parameters that differed from Fasham (1993) were otherwise parameterised from the literature where possible and/or selected as a best guess. The resulting parameter set, along with adjusted values (see below), is shown in Table 3. The maximum phytoplankton growth rate used by Fasham (1993) for station India 20 was 1.25 d −1 . The equivalent parameter in our model is the maximum rate of photosynthesis, V max P , which is usually expressed in units of g C (g Chl) −1 h −1 , requiring unit conversion. The Redfield C : N molar ratio of 6.625 is the obvious choice to convert between C and N. Carbon to chlorophyll ratios are more variable and a value of 50 g C (g Chl) −1 has previously been used in modelling studies (e.g., Fasham et al., 25 1990). However, C : Chl ratios are known to vary widely in response to ambient conditions. The recent study of Sathyendranath et al. (2009) found the North Atlantic ratio to typically vary between 50 and 100 g C (g Chl) −1 , so here we use an intermediate value of 75 g C (g Chl) −1 (parameter θ chl ). Converting units, V Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | lent to 3.9 g C (g Chl) −1 h −1 which is within a range of typical values for V P at ambient temperatures ranging between 1 and 5 g C (g Chl) −1 h −1 (e.g., Harrison and Platt, 1986;Cullen, 1990;Platt et al., 1990;Rey, 1991). We include temperature-dependence of this parameter and so, assuming that the rate of 3.9 d −1 occurs at a typical sea surface temperature during the bloom for station India of 10 Using the same conversion of units, Fasham's (1993) Yool et al., 2011Yool et al., , 2013a, and m P 2 , 0.025 (mmol N m −3 ) −1 d −1 (Oschlies and Garçon, 2005).
Zooplankton parameters I max and k Z were assigned directly from Fasham (1993) 10 with values of 1.0 d −1 and 1.0 mmol N m −3 , respectively. Note that assimilation efficiency as used by Fasham (1993) is in fact a growth efficiency whereas our use of absorption efficiency (parameter β) is more in keeping with contemporary zooplankton modelling (e.g., see Anderson et al., 2013) and refers to the fraction of material absorbed across the gut and is multiplied by a net production efficiency (parameter 15 k NZ ) to give growth efficiency. Values of 0.69 and 0.75 were assigned to parameters β and k NZ respectively (Anderson, 1994;Anderson and Hessen, 1995). In the model of Fasham (1993), zooplankton grazed on phytoplankton, bacteria and detritus. The model here has no bacteria and the relative ratio of grazing preferences for phytoplankton and detritus was maintained by assigning values for ϕ P and ϕ D of 0.67 and 0.33  , respectively. A detritus sinking rate of 1.0 m d −1 was used by Fasham (1993), a value at the low end of that typically used in models. Detritus is in reality composed of a range of sinking  Wilson et al., 2008), as well as slow-sinking material that is likely to be remineralised in the upper water column (Riley et al., 2012). A typical sinking rate used in ecosystem models is between 5 and 10 m d −1 (e.g, Fasham et al., 1990;Oschlies et al., 1999;Anderson and Pondaven, 2003;Llebot et al., 2010;Kidson et al., 2013). 5 We used a value for V D of 5.0 m d −1 , noting that results differed only slightly compared to using a sinking rate of 1.0 m d −1 . Note also that the detritus produced by quadratic zooplankton mortality is assumed to be very fast sinking and is instantly exported from the upper mixed layer. The remineralisation rate of detritus (parameter m D ) was set to 0.05 d −1 (Fasham, 1993). Finally, parameter w mix was set to 0.1 m d −1 (Fasham et al., 10 1990). Choices have to be made regarding the settings for calculating daily depth-integrated photosynthesis. A sinusoidal pattern of daily irradiance was set as default for this purpose, with a numeric integration over time of day. A Smith function was chosen as the P -I curve (Eq. 7) permitting a straightforward analytic depth integral for photosynthesis 15 (Appendix B). Photosynthesis at depth can be vertically integrated analytically, when light extinction in the water column is described by Beer's law with a constant coefficient. As default, we use the piecewise Beer's law treatment of Anderson (1993) in which the water column is divided into three depth zones (0-5, 5-23 and > 23 m) and a separate extinction coefficient calculated for each as a function of chlorophyll (Eq. 10). 20 Although this approach is more complicated that using a single extinction coefficient, it is easily justified a priori given the improved representation of light attenuation and its impact on predicted primary production (Anderson, 1993). Model sensitivity to these various assumptions regarding the calculation of light attenuation and photosynthesis will be examined in Sect. 4.3, including an assessment of the performance of the 25 algorithms of Evans and Parslow (1985) and Anderson (1993).
The model was run for three years, by which time it generates a repeating annual cycle of plankton dynamics. The chlorophyll data are SeaWiFS 8 day averages (O'Reilly et al., 1998). We had access to data from 1998 to 2013. Averaging data across years  Fig. 8. Nitrate (model DIN) is remarkably well predicted using these default parameter settings. Model chlorophyll shows a less good match with data. The timing of the spring bloom is too late although this could, at least in part, be due to the MLD forcing which was climatological, rather than for year 2006 (the chlorophyll data). Predicted chlorophyll also appears to be too high during the spring and summer 10 period. Parameter adjustment is therefore desirable in order to improve the fit with data.

Model calibration
Many modelers go about parameter adjustment on a trial-and-error basis, making ad hoc changes to parameters and observing the outcome. A more structured way of going about this is to undertake a systematic sensitivity analysis of parameters and 15 then, informed by this analysis, choose which parameters to vary. We use EMPOWER to demonstrate this practice here. Three variables were selected as simple measures of model mismatch with data: minimum DIN encountered during the seasonal cycle, N min , which is a logical choice because it is desirable to correctly predict DIN drawdown during the spring period, maximum chlorophyll at the peak of the spring bloom, chl max 20 and the average summer chlorophyll between days 200 and 300, chl av . Values of these three quantities, as outputs from the run shown in Fig. 8  where W S is the value of a given variable (in this case N min , chl max or chl av ) for the standard parameter set with parameter value p S , and W (p) is the value when the parameter is given value p. Results are shown in Table 4, ordered high to low for sensitivity of chl av . The chlorophyll data are too few in number to reliably infer the magnitude of the 5 spring bloom whereas there are many data points providing an average chlorophyll between days 200 and 300 of 0.29 mg m −3 . Looking at Table 4, chl av is sensitive to grazing parameters, notably k Z . As the first step to improving the model fit to data, k Z was decreased until predicted chl av was equal to 0.29 mg m −3 , resulting in a decrease in the value of this parameter from 1.0 to 0.52 mmol N m required N min is about 3.0 mmol N m −3 and in order to redress this mismatch with data parameter α was chosen for adjustment. This parameter shows high sensitivity for N min and relatively low sensitivity for chl av and chl max . Intuitively, α is a logical parameter to choose because nitrate drawdown occurs during rapid growth of phytoplankton at the onset of the spring bloom and increasing this parameter will therefore enhance draw- 20 down. An increase in α is also easily justified based on observational data (e.g., Rey et al., 1991). Increasing the value of α from 0.08 to 0.12 g C (g Chl) −1 h −1 (W m −2 ) −1 gave a predicted N min of 2.82 mmol N m −3 and an overall good fit to the data (Fig. 9). The only obvious mismatch is in the overwinter chlorophyll but extremely low values are a common feature of slab-type models. The mismatch can be improved by remov- 25 ing the linear phytoplankton mortality (see Sect. 4.4, and discussion therein). A further consideration is that phytoplankton may adjust their C : Chl ratio in winter to mitigate the effect of the low light intensities that they experience. We consider removing this mortality term unrealistic. It is no good getting the right result for the wrong reasons and 82 Introduction so chose to keep phytoplankton mortality as it is. There is also a hint that the timing of the bloom is a little late but, bearing in mind we used climatological cycle of annual mixed layer depth and light, whereas the data are for a single year, 2006, this is not particularly surprising. The associated seasonal cycles of P , Z and D, along with primary production, phy-5 toplankton grazing and mortality are shown in Fig. 10. The peak of Z lags seven days behind that of P , illustrating the decoupling of phytoplankton and zooplankton during the spring bloom period. Primary production remains relatively high over summer, tightly coupled to grazing, sufficient to keep phytoplankton biomass in check. It might be expected that Station Biotrans is simulated accurately with the same parameter values as those of Station India because of their relatively close proximity in the northern North Atlantic Ocean and this is indeed the case (Fig. 11).
The two HNLC stations can be expected to require alternative parameterisations to the two North Atlantic stations because the food web structure differs between the two types of system. In contrast to the diatom spring bloom in the northern North Atlantic, 15 iron-limited HNLC systems favour small phytoplankton which are tightly coupled to microzooplankton grazers (Landry et al., 1997(Landry et al., , 2011, "grazer controlled phytoplankton populations in an iron-limited ecosystem" (Price et al., 1994). Simulations for stations Papa, showing both the unfitted and fitted model, are shown below in in Fig. 12. The unfitted model solution corresponds to parameters as for the best-fit solution to Sta-20 tion India (Table 3). In common with the data, there is no predicted chlorophyll bloom. Predicted chlorophyll is however on the upper bound of the data and predicted drawdown of nitrate is too severe, suggesting that the modelled phytoplankton growth rate is too high. Low growth rate of phytoplankton may be expected relative to the North Atlantic because of iron limitation and so parameter V max P (0), acting as a proxy for iron 25 limitation, was halved in value to 1.0 g C (g Chl) −1 h −1 . With this parameter setting, the model does a remarkably good job at reproducing the data, without need for further parameter adjustment. A similar exercise was carried out for station Kerfix. Using the same parameter set as for station Papa, predicted chlorophyll was too high during the austral summer (Fig. 13). If grazing is dominated by microzooplankton, maximum grazing rate (parameter I max ) may be as high as 2.0 d −1 (Mongin et al., 2006). On this basis, I max was increased until predicted chl max (the maximum chlorophyll) equalled 0.35. A reasonable fit to the data 5 was achieved with I max equal to 1.4 d −1 .

Sensitivity to photosynthesis algorithm
Structural sensitivity analysis is performed to assess model sensitivity to the different assumptions for calculating daily depth-integrated photosynthesis. The best-fit simulation for Station India presented above ( Fig. 9) is used as the baseline for comparison. 10 Default settings in the baseline simulation were a numerical time integration (over the day), a Smith function for the P -I curve, and a sinusoidal pattern of daily irradiance and the piecewise application of Beer's law (Eq. 10; Anderson, 1993) for light attenuation in the water column. The first sensitivity test involved changing the P -I curve from a Smith function (Eq. 7) 15 to an exponential function (Eq. 8). Predicted seasonal cycles for chlorophyll and nitrate at station India are shown in Fig. 14. Results changed little with respect to the baseline simulation, with nitrate drawdown being slightly less when using the exponential P -I curve. Predicted chlorophyll was barely distinguishable between the two simulations. It is perhaps unsurprising that the model shows minimal sensitivity to choice of P -I curve 20 as the shapes of the two curves are similar. Slightly higher photosynthesis is predicted using the Smith function for mid-range irradiance (Fig. 5), consistent with higher drawdown of NO 3 . In a similar study by Anderson et al. (2010), however, remarkable sensitivity was seen to choice of the exact form of the zooplankton functional response.
Other studies have also shown "alarming" sensitivity to apparently small changes in 25 the specification of biological models (e.g. Wood and Thomas, 1999;Fussmann and Blasius, 2005). Reverting to the Smith function as the chosen P -I curve, model predictions were next compared for simulations using sinusoidal vs. triangular irradiance (Fig. 15). Once again, the difference between the two simulations is relatively minor, although predicted drawdown of nutrient was about 2 mmol m −3 less when using the triangular assumption. The triangular approximation underestimates the period of high light relative to 5 sinusoidal, for equivalent noon irradiance, with lower growth rate and associated drawdown of nutrient. It is worth noting, however, that the sensitivity shown here is at least as great as that for the choice of P -I curve, but has generally received much less attention in the literature. Model sensitivity of predicted primary production to the equations describing light at-10 tenuation in the water column was previously highlighted by Anderson (1993), although without extending to analysis using full ecosystem models. A marked difference was seen here when the piecewise Beer's law for calculating light attenuation (Eq. 10) was replaced with a simple Beer's law (Eq. 9) (Fig. 16). The difference between the simulations can be understood by comparing k PAR as a function of phytoplankton concentra-15 tion for the two algorithms (Fig. 17). The single Beer's law of Eq. (9) predicts a modest increase in k PAR from 0.04 m −1 at zero phytoplankton to 0.1 m −1 at P = 1 mmol N m −3 . The main difference with the piecewise Beer's law is the much greater light extinction in the upper 5 m of the water column, with k PAR of 0.13 m −1 at P = 0 mmol N m −3 , increasing to 0.23 m −1 at P = 1 mmol N m −3 . A lesser rate of light attenuation using the 20 simple Beer's law leads to greater penetration of light into the water column. The resulting higher photosynthesis over winter produced a larger spring bloom of phytoplankton and greater predicted drawdown of NO 3 . It is worth noting that the model sensitivity to this choice of light attenuation algorithm (both in terms of overestimating the spring bloom and the nutrient drawdown) is greater than that associated with the original pa-25 rameter adjustment exercise for station India, highlighting the importance of carefully selecting formulations for key processes prior to parameter tuning. Finally, there is the option to use the routines of Evans and Parslow (1985) and Anderson (1993)  using numerical integration over time. Evans and Parslow used a Smith function for photosynthesis in combination with a triangular pattern of daily irradiance. This corresponds exactly to the simulation in Fig. 15 above for triangular irradiance. Thus, running the model using the Evans and Parslow equations (Appendix C) produces a result indistinguishable from the numerical simulation. Matters are not so simple when using 5 the Anderson (1993) equations to calculate daily depth-integrated photosynthesis. The assumptions here are an exponential P -I curve and sinusoidal light, corresponding to the exponential P -I curve simulation in Fig. 14. But there is the additional assumption that parameter α, in addition to k PAR , is spectrally dependent and varies in the water column. Thus, running the model with both light attenuation and photosynthesis calculated as in Anderson (1993) gives rise to a different simulation (Fig. 18). It is noticeable that, when using the A93 method, primary production is higher over winter, a result of elevated α, giving rise to an earlier spring chlorophyll bloom and greater drawdown of nitrate. Nevertheless, the simulation is entirely credible and we can recommend the use of the Anderson (1993) for use in marine ecosystem models.

Mortality terms
The model includes two mortality terms, linear and quadratic, for each of phytoplankton and zooplankton. This approach has previously been used in other models (e.g., Yool et al., 2011Yool et al., , 2013a, giving maximum flexibility. The obvious question is whether all four terms are actually needed. As a simple structural sensitivity analysis, we removed each 20 of the four mortality terms in turn and show the impact on the predicted seasonal cycles of chlorophyll and nitrate, for Station India. Starting with the phytoplankton terms, setting m P or m P 2 to zero affected both the predicted timing and magnitude of the spring bloom (Fig. 19). One can argue that, although the predicted magnitude of the spring bloom looks a little low, removal of the linear term (setting m P = 0) improved the model 25 fit for chlorophyll, notably over winter. It seems hard to justify that loss rates should go to near zero at low population densities (the consequence of using a quadratic term only) because all organisms have metabolic requirements. Nearly all marine ecosys-86 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tem models do, therefore, include a linear term for phytoplankton mortality and, for our baseline simulation (Sect. 4.2), we chose to keep this term on a purely conceptual basis. Given deep mixing, it is surprising that phytoplankton biomass, as seen in the data, is maintained over winter in high latitude waters. The reasons why this is so remain a matter of conjecture with candidate theories including cyclic motion associated with 5 convective mixing (Huisman et al., 2002;Backhaus et al., 2003), and phytoplankton motility or buoyancy to remain near the ocean surface (see Ward and Waniek, 2007, and references therein). The slab model, and indeed 1-D and 3-D models, have difficulty dealing with this issue but there is no evidence that this seriously compromises results when it comes to the predicted timing and magnitude of the spring bloom and 10 associated ecosystem dynamics later in the year. In contrast to the representation of linear mortality, many models do not include a non-linear phytoplankton mortality term but it seemed to perform well here. When it was removed, the predicted phytoplankton spring bloom was rather too high. Results show that non-grazing phytoplankton mortality had a pronounced influence 15 on simulated phytoplankton biomass both prior to and during the initiation of spring bloom (Fig. 19). It is at this time of year correspond when grazing losses are minimal ( Fig. 10) such that phytoplankton dynamics are driven by the balance of growth and non-grazing mortality. Phytoplankton levels are low at the end of winter and hence removal of the quadratic mortality term had virtually no effect on pre-bloom phytoplankton 20 levels whereas removal of the linear term had a marked impact leading to a reduction in the peak of the bloom of about a third. This reduction can be explained by the fact that the higher phytoplankton density pre-bloom associated with removal of linear phytoplankton mortality enabled higher pre-bloom zooplankton grazing. In contrast, removal of the quadratic mortality term nearly doubled size of the bloom, as might be expected 25 based on the sensitivity analysis (Table 4). This strong effect on biomass indicates that it was the density-dependent (quadratic) mortality term that caused phytoplankton mortality to initially rival grazing (Fig. 10). Removing the zooplankton mortality terms in turn also significantly impacted on model predictions (Fig. 20). While changes in the linear mortality term had a noteworthy effect on both the bloom peak and minimum drawdown (as also shown in the sensitivity analysis Table 4), it was the quadratic zooplankton mortality term that had the most influence. Removal of quadratic mortality resulted in significantly lower 5 phytoplankton levels post-bloom (Fig. 20, Table 4) which is unsurprising since more zooplankton means more grazing. Perhaps less obvious is the result that removal of quadratic closure resulted in large changes in predicted post-bloom nitrate levels, even exceeding those arising from consideration of piecewise vs. simple light attenuation (Fig. 16). Predation-related losses, the quadratic term, were assumed to be instantly 10 exported and thereby lost from the surface mixed layer of the model. Thus, when these losses are set to zero (parameter m Z2 = 0), nitrate drawdown is significantly diminished because, instead of being instantly exported, zooplankton quadratic mortality is allocated to sinking detritus, part of which is remineralised in the mixed layer. Overall, the work highlights the need for careful consideration of the parameterisation of closure 15 in models, including the fate of material thereof.

Discussion
Simple models are all too often brushed aside in marine science today. When it comes to the representation of the marine ecosystem, complex models have come to the fore that have, for example, any number of plankton functional types, multiple nutrients, 20 dissolved organic matter and bacteria, etc. (e.g., Blackford et al., 2004;Moore et al., 2004;Le Quere et al., 2005). There is a similar trend with ocean physics toward large, computationally demanding models. Many publications in recent years have involved the use of 3-D models (e.g., Le Quéré et al., 2005;Wiggert et al., 2006;Follows et al., 2007;Hashioka et al., 2013;Yool et al., 2013b;Vallina et al., 2014), although 1-D 25 models are also well represented (e.g., Vallina et al., 2008;Kearney et al., 2012;Ward et al., 2013). Of course, the improved realism that is gained by using complex models is in general to be welcomed, with the caveat that improvements in prediction can only be achieved if the processes of interest can be adequately parameterised (Anderson, 2005). Despite the trend to complex ecosystem models embedded in advanced physical frameworks, there nevertheless remains a place today, we argue, for simple models. 5 Simple models are fast to run, transparent and easy to analyse. Marine ecosystem modelling can be somewhat of a black art in deciding what to include in terms of state variables, which formulations to apply for key processes such as photosynthesis, grazing and mortality, and in finding suitable parameter values. Simple models allow us to fully examine the subtle inner workings of models, assessing the merits of different choices for model specification. The pioneers of the field such as Riley, Steele and Fasham played extensively with their (simple) models, trying out different formulations and parameterisations, just to see what would happen. The simplicity afforded by using a zero-dimensional slab physics framework provides an ideal playground for familiarisation with ecosystem models, allowing for a multiplicity of runs and ease of analysis. 15 Here, we have presented an efficient plankton modelling testbed, EMPOWER, coded in the freely available language R. It provides a readily available and easy to use tool for evaluating model structure, formulations and parameterisations. EMPOWER has several advantages in that it is fast, easy to run, results are straightforward to analyse and, last but by no means least, the code is transparent and easily adapted to incorporate 20 new formulations and parameterisations. As such, the main purpose of EMPOWER is to provide an ecosystem model testbed that allows users to fully familiarise themselves with their models, allowing them to subsequently be incorporated with greater confidence into 1-D or 3-D models, as required. It may be that some amount of reparameterisation is required when transferring the model ecosystem between physical 25 codes, but this ought usually to be minimal in extent and will itself be greatly informed by the previous slab modelling work. Much better this approach, than starting out from scratch using computationally expensive and time-consuming 1-D or 3-D codes to undertake ecosystem model parameterisation. In order to demonstrate the utility of EMPOWER, we carried out both a parameter tuning exercise and a structural sensitivity analysis, the latter examining the equations for calculating daily depth-integrated photosynthesis, and mortality terms for both phytoplankton and zooplankton. In the parameter tuning exercise, a simple NPZD model, broadly based on the ecosystem model of Fasham (1993), was fitted to data (seasonal cycles) for chlorophyll and nitrate at four stations: India (60 • N, 20 • W), Biotrans (47 • N, 20 • W), Papa (50 • N, 145 • W) and Kerfix (50 • 40 S, 68 • 25 E). Formal parameter sensitivity analysis was carried out, highlighting which parameters phytoplankton stocks and nitrate drawdown are sensitive to. The model was successfully tuned to all four stations, the two HNLC stations (Papa and Kerfix) requiring different parameterisations, notably a halving of maximum photosynthetic rate (acting as a proxy for iron limitation) relative to the North Atlantic sites.
The parameterisation of the different stations highlighted the somewhat ad hoc process that most modellers go through when assigning parameter values. Some parameters may be set directly from the results of observation and experiment. More often 15 than not, however, the "path of least resistance" when assigning parameters is to simply select values from previously published modelling studies. Equations for processes such as photosynthesis, grazing and mortality can likewise be selected "on-the-shelf" from the published literature. Previous publication does not, of course, guarantee that equations or parameter values are necessarily best suited for a particular modelling ap-20 plication. Moreover, it is all too easy for less than ideal, even dysfunctional, formulations to become entrenched within the discipline and used in common practice (Anderson and Mitra, 2010). As a result, parameter tuning is almost inevitable in ecosystem modelling and we have shown how rigorous sensitivity analysis can help in this regard. Of course, even with a table of parameter sensitivities, there is still a considerable subjec-25 tive element to choosing which parameters to adjust. The most sensitive parameters should be selected, but the degree of uncertainty in parameter values is an additional consideration. It is no good tuning a sensitive parameter if its value is already well known from observation and experiment.  A necessary complement when ensuring that models show acceptable agreement with data is to remember that it is important that the theories and assumptions underlying the conceptual description of models are correct, or at least not incorrect (Rykiel, 1996). Indeed, it is the conceptual realisation of models that in many ways poses the greatest challenge, requiring expertise and practice to overcome observational or ex-5 perimental lacunae (Tsang, 1991). Subsequent to the parameter tuning exercise, we studied the sensitivity of the Station India simulation to chosen formulations for depthintegrated photosynthesis and both phytoplankton and zooplankton mortality. In the case of the photosynthesis calculation, some aspects showed relatively low sensitivity, namely the choice of P -I curve and whether to assume a triangular or sinusoidal pat-10 tern of irradiance throughout the day. In contrast, the way in which light attenuation in the water column is calculated showed marked sensitivity. Using a simple Beer's Law attenuation coefficient throughout the water column is clearly oversimplified because the spectral properties of irradiance vary with depth. Moving to a piecewise Beer's Law with separate attenuation coefficients for depth ranges 0-5, 5-23 and > 23 m (Anderson, 1993) led to more rapid light attenuation near the ocean surface. Depth-integrated photosynthesis declined accordingly, delaying the onset of the spring bloom and reducing its magnitude, along with drawdown of nutrient. The difference is of course in part due to parameter values, rather than the inherent difference in the equations. Additional sensitivity analysis and parameter tuning could be used to investigate this further but in 20 fact such an analysis was undertaken by Anderson (1993) who showed that no amount of parameter tuning can adequately account for the fact that attenuation will vary with depth, and cannot be assumed to be constant, because of the spectral properties of the irradiance field. In contrast to the sensitivity seen to equations for light attenuation, choice of P -I curve made only a negligible difference to model predictions. 25 When it comes to biogeochemical modelling studies in GCMs, it is possible that all manner of different methods are used to calculate light attenuation in the water column and resulting photosynthesis. Methodologies are often not reported in full within published texts, the assumption being that they are in some way routine and straightforward  , 2001). This divides light into two wavebands, "red" and "green-blue" that are attenuated separately by seawater, and a Smith function (Eq. 7) is used to calculate photosynthesis. But the published description omits a number of key details (although the model code was supplied), for instance that there is a 50 : 50 division of light between the two wavebands at the ocean surface, that the photosynthetically active fraction is 10 0.43 of total irradiance, that extinction coefficients for the two wavebands are a function of chlorophyll and that photosynthesis is calculated within each model layer (the model uses fixed layer depths, with 13 layers in the upper 100 m) as a function of average light within the layer. As a point of interest, we ran our model for station India again, this time using the 15 MEDUSA-2.0 method of light attenuation and a Smith function for the P -I curve (see Appendix E for details of the parameterisation of light attenuation). The calculation included replication the layer structure within the GCM in order to achieve a full comparison. Results (not shown) were almost identical to the baseline simulation for station India (Fig. 9) for all four state variables, with the minor exception that nitrate drawdown 20 was slightly greater (0.5 mmol N m −3 ) with the MEDUSA parameterisation. The similarity between the two simulations is because, remarkably, calculated light attenuation using the two red and green wavebands (MEDUSA) differs little from that using the Anderson (1993) piecewise Beer's law. Here, in a nutshell, is a classic example of the utility of EMPOWER. This result should alert GCM modellers to the fact that near iden- 25 tical results can be generated for light attenuation in the water column using these two contrasting sets of equations and a choice can be made as to which is most suitable for implementation based on computational efficiency. From a theoretical point of view, the result is also interesting. The equations of Anderson (1993) are an empirical approxi- mation of the full spectral model of Morel (1988) which divided PAR into 61 wavebands. It would appear that this model can be stripped down to just two wavebands, red and green, without serious degradation in accuracy when it comes to predicting light attenuation.
We also used EMPOWER to undertake an analysis of model sensitivity to the pres-5 ence/absence of linear and nonlinear mortality terms for phytoplankton and zooplankton. Whereas the use of linear phytoplankton mortality terms is commonplace in models, we investigated the performance of an additional quadratic phytoplankton mortality term. This term is intended to represent loss processes that scale with phytoplankton biomass that are not already accounted for in the model. Given that both self-shading 10 and grazing are explicitly modelled, we considered the quadratic term to represent mortality due to viruses. Model results were sensitive to this parameterisation, highlighting the potential importance of viruses in marine systems, which is consistent with field evidence (Bratbak, 1993(Bratbak, , 1996Danovaro et al., 2011). It has long been recognized that the parameterisation and functional form of zoo- 15 plankton mortality, the model closure term, can have a pronounced effect on modeled ecosystem dynamics (e.g. Steele and Henderson, 1981, 1992Murray and Parslow, 1999;Edwards and Yool, 2000;Neubert et al., 2004). Quadratic closure is a common choice, although other non-linear functional forms are also in use. While it is commonly stated that quadratic closure is dynamically stabilising, i.e., it prevents 20 both blooms and extinction of prey, there is a limit to this influence (Edwards and Yool, 2000) since other processes can come into play. In our case, it is obvious that quadratic closure had a stabilising effect on the model. Its removal caused the bloom peak to be higher and also post-bloom phytoplankton levels to decline to near-zero. In contrast to the community's broad recognition of the potential sensitivity to choice 25 of closure scheme, far less attention has been paid to model sensitivity regarding the fate of zooplankton mortality. There are likely various sources of zooplankton in reality including grazing by higher predators, starvation or disease. One can consider the grazing loss to be partitioned between an infinite series of higher predators (e.g.,  Fasham et al., 1990), with partitioning between detritus and dissolved nutrients in both organic and inorganic form. These fates will occur with time delays and potentially also with spatial separation due to migration of predators. Moreover, any detrital production by higher predators would comprise significantly larger "particles" than those due to plankton death, and would therefore be associated with much higher sinking 5 rates. Non-grazing mortality might lead to production of detritus in situ. There is no consensus on best practice, despite the fact that different approaches to partitioning of zooplankton losses can have a significant effect on modelled ecosystem function. Future structural sensitivity studies should be conducted to explore how the f ratio (the fraction of primary production fuelled by external nutrient) and e ratio (i.e. relative ex-10 port to total primary production) are affected by the various assumptions relating to zooplankton mortality and model closure.
We have described the utility of slab models as a testbed underpinning marine ecosystem modelling research. This is however by no means their only use. Slab models are ideal for teaching ecological modelling. They embrace the complex interplay 15 between primary production and the physical-chemical environment, combined with top-down control by zooplankton. Students often have difficulty grasping the relative significance of causal effects in ecosystems (Grotzer and Basca, 2003), e.g. the relative roles of bottom-up vs. top-down processes in structuring food webs. A certain amount of lecture material is of course needed, but there is no substitute for hands-on 20 modelling, providing an interactive approach whereby students can actively investigate ideas and interact between themselves and a teacher (Knapp and D'Avanzo, 2010). Insight can be gained by getting students to try simple things like switching grazing off, doubling phytoplankton growth rates, etc. The slab modelling framework provided herein is ideal for this purpose. The code is transparent, modular and readily adjusted 25 to include alternate parameterisations, it is easily set up for alternate ocean sites, the model runs fast with graphs of results appearing on the screen on completion, results are readily written to output files for more in depth analysis and, by coding in R, the models can be accessed and run without need for purchasing proprietary software. Finally, the great advances in marine ecology that the pioneers of plankton modelling achieved using slab models should not be forgotten. Riley, Steele and Fasham laid the foundations of today's marine ecosystem modelling using plankton models embedded within simple physics. Even in the modern arena, this use of simple physics cannot be dismissed as being too simple for practical application and there is no reason why 5 further scientific advances cannot be made on this basis. Models are, fundamentally, all about simplifying reality.

Appendix A: Irradiance calculations
Both the Evans and Parslow (1985) and Anderson (1993) subroutines for calculating daily photosynthesis require noon irradiance and daylength as inputs. When there are 10 data available, these data can be used as forcing for a model, akin to what is done for temperature. However, most typically light data is not available and so a light submodel must be used to prescribe the light forcing. A climatological approach is often used whereby these inputs are specified using trigonometric/astronomical equations. This task is not as straightforward as it might first appear. The basic equations are presented 15 in texts such as Brock (1981) and Iqbal (1983). Some adjustments were provided by Shine (1984) and we recommend the equation for short-wave irradiance at the ocean surface on a clear day published therein:  , Josey et al., 2003); the calculation of I clear is not sensitive to this parameter. The equation for R V is: where J is day of year (Julian day; i.e. 1 = 1 January). Solar zenith angle depends on latitude (ϕ), solar declination angle (δ) and on time of day (γ, where the Earth moves 5

15
• h −1 and γ is difference from noon): The cos(γ) term becomes irrelevant when considering noon irradiance. Solar declination angle is given by: δ = 23.45 sin(2π(284 + J)/365) (A4) 10 where h is hour angle which is the difference between the given time and noon (where 1 h is 15 • ). Note that δ is expressed in degrees in the above equation (1 radian = 180/π degrees). The flux of photosynthetically active solar radiation just below the ocean surface at noon, I noon , can now be calculated: where f PAR is the fraction of solar radiation that is PAR (λ between 400 and 700 nm), φ is ocean albedo and C FAC is the effect of clouds on atmospheric transmission. Parameters f PAR and φ are relatively invariant with typical values of 0.43 for f PAR and 0.04 for φ (e.g., Fasham et al., 1990). Dealing with the effects of clouds is a thorny issue for modellers. Simple empirical approaches have been developed, two of the most popular being those of Reed (1977) and Smith and Dobson (1984). We have opted for the former in which C FAC is a function of zenith angles (specified in degrees): where W is cloud fraction in oktas. A value of W = 6 was used for all four stations.

Appendix B: Analytic integrals for photosynthesis with depth
The average photosynthesis within a layer of depth H is: where V P is photosynthesis as a function of light intensity (specified as the P -I curve). Two P -I curves are provided with EMPOWER, a Smith function (Eq. 7) and exponential 10 function (Eq. 8). Analytic solutions to Eq. (B1) are provided here for each of these two P -I curves. In both cases a Beer's law attenuation with depth is assumed (parameter k PAR ), i.e., I(z) = I(0)e −k PAR z where I(0) is the irradiance entering the layer from above.

B1 Smith P -I curve
By performing a change of variables such that x = αI(z), the integral above becomes: This integral is solved analytically using a trigonometric transformation and then integration by parts, giving: B2 Exponential P -I curve 5 In order to integrate Eq. (B1) using an exponential P -I curve it is first useful to define (Platt et al., 1980): The integration over depth is then (see Platt et al., 1990):  Evans and Parslow (1985) provide an algorithm for calculating daily depth-integrated photosynthesis with the assumptions of a Smith P -I curve (Eq. 3), a triangular pattern of irradiance from sunrise to sunset and light extinction calculated with a single 5 Beer's law coefficient. The average daily rate of photosynthesis within the mixed layer is calculated as: where t, measured in days, is 0 at sunrise and τ at noon and H is layer depth. Assuming a triangular pattern of irradiance about noon, Eq. (C1) can be recast as (Evans and 10 Parslow, 1985): I noon is the photosynthetically active radiation (PAR) just below the ocean surface at noon. This integral solves as (Evans and Parslow, 1985):

C2 Anderson (1993) photosynthesis calculation
The subroutine of Anderson (1993) was developed as an empirical approximation to the spectrally resolved model of light attenuation and photosynthesis of Morel (1988) 5 used in combination with the polynomial method of integrating daily photosynthesis of Platt et al. (1990). It is based on an exponential P -I curve (Eq. 8), assumes a sinusoidal pattern of irradiance throughout the day and calculated light attenuation using a piecewise Beer's law (Eq. 10). The irradiance leaving the base of each layer is: 10 where I base,0 is the irradiance immediately below the ocean surface and z base,i is the depth of the base of the layer i (where z base,0 = 0). The subroutine of Anderson (1993) also takes account of the fact that, in reality, α depends on the spectral properties of light and therefore varies with depth in the water column. This parameter is the product of photosynthetic absorption cross section a c (λ), 15 which is spectrally dependent (λ denotes wavelength), and quantum yield ϕ A (Platt and Jassby, 1976;Platt, 1986): Ordinarily (e.g., Table 2), α is the initial slope of the P -I curve for white light (i.e., spectral distribution as for irradiance at the ocean surface). The corresponding value of α for the wavelength at which absorption is maximum, α max , is (Anderson, 1993): The value of α for any given wavelength of PAR, α(λ), is then: where a * (λ) is the dimensionless chlorophyll absorption cross section for wavelength λ. An additional complication, however, is that a * (λ) only applies when irradiance is specified as a scalar flux (Morel, 1991). Irradiance in the model is a downwelling flux and so Anderson (1993) converted between the two by defining a new version of the 10 chlrophyll absorption cross section (which can be used in Eq. (C9) in place of a * (λ), in combination with downwelling irradiance): Again using the piecewise three-layer scheme described above for k PAR , an average value of a # can be calculated for each layer by deriving an empirical approximation of 15 Morel's (1988) full spectral model. As a first step, a # at the ocean surface is calculated as: where the polynomial coefficients are given in Table C1. The a # at the base of each layer and the average a # in each layer are then calculated as: where a # calc,i is a lengthy empirical calculation: The coefficients, g x , are provided in Table C1. With irradiance assumed to vary sinusoidally through the day, the average rate of photosynthesis within a layer i is: where D is daylength (hours) and Ω j are the polynomial coefficients (Platt et al., 1990;  i. File NPZD_parms.txt. This file includes a single line header and then lists the value of each model parameter in turn, followed by a text string for the purpose of annotation. When changing the parameter list in the model, the corresponding section in the R code must be altered accordingly. 15 ii. File NPZD_extra.txt. This file holds initial values for state variables, additional parameters, and various flags: choice of station, choices for photosynthesis calculations (P -I curve, light attenuation, etc.) and grazing formulation. The user is at liberty to add to or remove from this list of flags as is desired. This file also contains flags for core model functions: run duration, time step, out-20 put type (none, last year, whole simulation), output frequency and integration method (Euler or Runge Kutta). These latter functions are required by the core code and should not be removed from this file.
iii. File stations_forcing.txt. This file has a header line for information, and then holds monthly values for forcing, in our case mixed layer depth and temper- 7. Graphical output. The model automatically generates graphical output on the computer screen on completion of each simulation. An advantage of R is that 5 the syntax for generating plots is straightforward and the user should have no problem, working from the plots provided, in generating extra graphs, as desired.