Edinburgh Research Explorer Sensitivity analysis of an ocean carbon cycle model in the North Atlantic: an investigation of parameters affecting the air-sea CO2 flux, primary production and export of detritus

. The sensitivity of the biological parameters in a nutrient-phytoplankton-zooplankton-detritus (NPZD) model in the calculation of the air-sea CO 2 ﬂux, primary production and detrital export is analysed. We explore the effect on these outputs of variation in the values of the twenty parameters that control ocean ecosystem growth in a 1-D formu-lation of the UK Met Ofﬁce HadOCC NPZD model used in GCMs. We use and compare the results from one-at-a-time and all-at-a-time perturbations performed at three sites in the EuroSITES European Ocean Observatory Network: the Central Irminger Sea (60 ◦ N 40 ◦ W), the Porcupine Abyssal Plain (49 ◦ N 16 ◦ W) and the European Station for Time series in the Ocean Canary Islands (29 ◦ N 15 ◦ W). Reasonable changes to the values of key parameters are shown to have a large effect on the calculation of the air-sea CO 2 ﬂux, primary production, and export of biological detritus to the deep ocean. Changes in the values of key parameters have a greater effect in more productive regions than in less productive areas. The most sensitive parameters are generally found to be those controlling well-established ocean ecosystem parameterisations widely used in many NPZD-type models. The air-sea CO 2 ﬂux is most inﬂuenced by variation in the parameters that control phytoplankton growth, detrital sinking and carbonate production by phytoplankton (the rain ra-tio). Primary production is most sensitive to the parameters that deﬁne the shape of the photosynthesis-irradiance curve. control the of and the


Introduction
The ocean absorbs approximately 2 Pg per year of carbon from the atmosphere (Takahashi et al., 2002;Gruber et al., 2009) -around a third of current anthropogenic emissions. The exchange of CO 2 between the atmosphere and ocean is driven by the difference in CO 2 concentration between the air and surface water, which, in turn, is influenced by the uptake of carbon by photosynthesising phytoplankton in the euphotic zone. Phytoplankton and zooplankton mortality, undigested waste resulting from incomplete injestion of phytoplankton prey by grazing zooplankton, and excretion generate organic waste that sinks through the water column (export production). Sinking organic waste is largely remineralised in the water column with a tiny fraction (typically <1 % of export production) reaching the ocean bed where it is gradually broken down, either dissolving or remaining locked up in sediments, completing the biological carbon pump. A relatively simple way of modelling this process is to use a fourcompartment Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model. Such models are widely used in GCMs due to their computational efficiency and typically contain around twenty biological parameters. Many of these parameters represent bulk properties across the whole ocean, and are poorly constrained in value (Frenette et al., 1993;Fennel et al., 2001) leading to large uncertainties in model predictions. Understanding which of these parameters have the greatest influence on model output is important to understanding model results and developing improved models.
In this study we use the Hadley Centre Ocean Carbon Cycle (HadOCC) model of Palmer and Totterdell (2001) which is used to represent the role of the ocean ecosystem in the UK Met Office HadCM3 climate prediction GCM (Gordon et al., 2000) and its faster-running derivative the FAst Met Office and Universities Simulator (FAMOUS) GCM (Smith et al., 2008). HadOCC is an NPZD model that calculates  the flow of nitrogen and associated flows of carbon and alkalinity between the four model compartments. Twenty free parameters (detailed in Table 1) govern the biological processes of phytoplankton growth and mortality, zooplankton grazing and excretion, and detrital sinking and remineralisation. Analogous parameterisations are used in other NPZD and more complex ocean carbon cycle models (e.g. Popova et al., 2002;Waniek, 2003). For this work HadOCC has been coupled to the General Ocean Turbulence Model (GOTM) (Burchard et al., 1999;Kettle and Merchant, 2008) to create a one-dimensional (1-D) model (hereafter referred to as HadOCC-GOTM) that can be run at any location and is computationally efficient enough to perform the many model runs required for sensitivity analysis. The 1-D HadOCC-GOTM column is forced with meteorological data and uses a relaxation to a nutrient profile below the maximum depth at which production occurs to represent the horizontal transport of nutrient into the column -see Sect. 2.1 below and Appendix B.
Previous studies have assessed the importance of parameters in ocean ecosystem models at one specific site (e.g. Druon and Le Fèvre, 1999). Here, we analyse three sites that have been used for previous modelling studies (Kettle and Merchant, 2008) in different ocean regimes with different meteorological forcing and nutrient supplies to compare the regional importance of the model parameters. The Central Irminger Sea (hereafter referred to as the CIS) is located at 60 • N 40 • W between Greenland and Iceland in the North Atlantic subpolar gyre. Nutrients are brought up from  (2008)  16 Anderson and Pondaven (2003)  17 Waniek and Holiday (2006) single assimilation parameter 18 Druon and Le Fèvre (1999) single assimilation parameter 19 Anderson and Pondaven (2003) to DOM 20 Anderson et al. (2007) 21 Kawamiya et al. (2000) 22 Druon and Le Fèvre (1999) at all depths 23 Waniek and Holiday (2006) at all depths 24 Fujii and Chai (2007) 25 Fujii et al. (2005) 26 Phytoplankton concentration is assumed to be 1 mMol N m −3 for purposes of conversion 27 Cloern et al. (1995) deep water by winter storms extending the mixed layer depth allowing a phytoplankton bloom in May/June (Waniek and Holiday, 2006). We note that the CIS lies on the fringes of the seasonal ice zone and so due to the presence of ice algae the ecosystem at the CIS may in truth differ somewhat from the open ocean ecosystem that HadOCC represents (Deal et al., 2011 (Davenport et al., 2002;Zielin-ski et al., 2002). Further information on all these sites is available from the EuroSITES European Ocean Observatory Network (www.eurosites.info).
In this study we test the sensitivity of three model outputsnet annual air-sea CO 2 flux (mol CO 2 m −2 yr −1 ), annual primary production and annual deep export (mg C m −2 yr −1 ) to variation in the values of the model parameters. Air-sea CO 2 flux is used in the calculation of absorption of anthropogenic emissions by the ocean, and all three are key processes in climate models. Here we define deep export to be the annual accumulation of detrital material that sinks below the maximum mixed layer depth (MLD) of the water column and is therefore removed from the surface waters. As the mixed layer at the CIS extends to the deep ocean during the winter, returning sinking detrital material to the surface, the sensitivity of the calculated export to the input parameters at this site is not studied.
Sensitivity analysis (SA) is the study of how variation in model output can be apportioned to different sources in model input (e.g., Saltelli et al., 2001). SA is performed by varying the value of the parameters across their plausible ranges and recording the resultant change in the model output. Low-sensitivity parameters are those where plausible perturbation to their value does not have a significant effect on the model output, while perturbations of highly sensitive parameters have a large effect on model output. SA qualifies the relative importance of parameter uncertainty in order to better understand and improve the model. Four steps are required to perform SA. First, the uncertainty in each parameter is assessed and a reasonable range of values for each parameter developed. Second, the parameter ranges are sampled to generate a set of input parameter sets. Third, the model is run for each parameter set and the output recorded. Fourth, this information is used to perform a sensitivity assessment. At its simplest, SA is done One-At-a-Time (OAT), with each parameter varied individually while all others remain constant at a some reference value. However, this method assumes that parameter interactions have an insignificant effect on model output. Here we perform both OAT analysis and a more complex global analysis to explore the impact of each parameter across the full range of all the other parameters.

Model setup
We use a slightly modified version of HadOCC from that described by Palmer and Totterdell (2001) (see Appendix A for the model equations). HadOCC-GOTM is run at CIS, PAP and ESTOC in the North Atlantic ocean. Initial conditions at each location for alkalinity (Lee et al., 2006), dissolved inorganic carbon (DIC) (Key et al., 2004) are obtained from the Carbon Dioxide Information Analysis Center (CDIAC), and nitrate profiles (to which the model relaxes to resupply) from  (Garcia et al., 2006). The number of depth levels used by GOTM are set so that the total depth of the GOTM water column is as close as possible to the actual depth of the site. The heat flux method is that used as standard in GOTM (Burchard et al., 1999). To drive the model, meteorological data are taken from ECMWF 40-year Re-analysis (ERA-40) data (Uppala et al., 2005); these comprise air pressure, wind speed, relative humidity, air temperature and total cloud cover, at 6 hourly intervals. Appendix B gives further detail on the resupply of nutrient and meteorlogical forcing. For each location GOTM-HadOCC is run with one thousand different sets of parameters (as described below). The model is spun up for eight years with yearly repeating meteorological data from 2004 to remove sensitivity to the initial conditions. The three model outputs under investigation (airsea CO 2 flux, export and primary production) are then taken from the ninth year. This allows output differences between the runs for a given location to be attributed exclusively to the change in the parameter values.

Model parameters: ranges and sampling
We investigate twenty parameters in this study. The model equations in which these parameters are used are detailed in Appendix A. The parameters are listed in Table 1 with the current HadOCC value and values used in other comparable models. For each parameter, we estimate its likely range from expert opinions and the literature values in Table 1. The parameters are in two distinct categories -positive definite parameters (e.g., detrital sinking velocity), and those with a value confined in the interval 0 to 1 (e.g., the assimilation efficiency of zooplankton feeding).
The parameters that relate directly to phytoplankton growth (K nit , P s max , α, m 0 ) are generic to many NPZD and more complex models giving a good resource in determining their possible range. For these parameters the smallest and largest values found in the literature are used for the range (see Tables 1 and 3). Parameterisations of processes such as zooplankton grazing are however less universal so there is less information available to set a likely range. Following consultation with experts who use HadOCC (T. Anderson, R. Barciela, J. Hemmings, personal communication, 2008) we give these parameters a range of the nominal HadOCC value ±90 % of the difference between the nominal HadOCC value and the closest theoretical limit, in most cases either 0 or 1. To establish a range for the deep remineralisation rate Rm deep applied below 100 m depth 1 we use the maximum value for the shallow remineralisation rate Rm shall as the upper limit for Rm deep /138.9 m (the midpoint of the depth level below 100 m is at 138.9 m) from which the 90 % range is established as with other parameters above.
From these parameter value ranges one thousand parameter sets are selected using a maximin Latin Hypercube sampling method (Saltelli et al., 2001) which maximizes efficient coverage of the whole parameter space by ensuring selection of parameter values from their full range -see Appendix C. The use of the hypercube sampling method gives equivalent parameter space coverage to a randomly sampled selection of parameter value sets a factor of ten greater in size.
These parameter sets are then run in HadOCC-GOTM and the results recorded alongside the appropriate parameter set. The calculated outputs are screened to check for failed runs that have produced wholly unfeasible results created by parameter combinations causing numerical problems (e.g. division by zero). We use the Takahashi climatology (Takahashi et al., 2002) to screen the CO 2 flux results and set the condition that export and primary production must be greater than or equal to zero. No models runs were found to exceed these checks.

SA methods
Here, we start with the simplest method of SA -One-at-A-Time (OAT), then use the GEM-SA emulator (see below) to explore the possible effects of interactions between parameters. OAT analysis investigates the effect of changing the value of each parameter across its range while all others remain constant. OAT gives straightforward insight into the effect of parameter value variation on model outputs and the relative importance of parameters, but the amount of the full model parameter space covered is very limited as parameter interactions are not explored (Saltelli et al., 2001).

GEM-SA
In this work, the SA package Gaussian Emulation Machine for Sensitivity Analysis (GEM-SA) www.ctcd.group.shef.ac. uk/gem.html (O'Hagen, 2006), is used to perform and assess a global SA of HadOCC-GOTM. An issue with performing global SA is the computational expense of performing large numbers of model runs to generate the input-output data set for study. To reduce this computational cost, GEM-SA uses an input-output data set as a training set to build an emulation of the relationship between the inputs and output in the model, allowing the input-output set, and the corresponding cost of generation, to be smaller. This emulation is a statistical approximation to the actual model, which allows the SA to use information from the emulated regions between the points defined by the processed input-output data, as well as the points themselves. Because of this large computational benefit, statistical emulators like GEM-SA are increasingly being used to perform SA in both ecological modelling e.g. Petropoulos et al. (2009) and other research fields such as engineering e.g. Finley et al. (2009).
To emulate the relationship between the inputs and the output, GEM-SA interpolates between the points in the training set using a Gaussian process. At each of the points in the training set the emulator gives the same result as that found in the training set -at these points the uncertainty in the fit of the emulator is zero. Between these training set points the emulator gives a "best guess" to the true value of the output that would have been calculated if the corresponding parameter value had been explicitly used, with the uncertainty of this "best guess" as a normal distribution. This means that uncertainty in the goodness of fit of the emulator increases as the emulated parameter moves further away from the points of the training set. In practice, provided the relationship between the parameter and output is smooth, and a suitably sized training set is used, emulated data are virtually indistinguishable from genuine data.
SA is performed on the emulated model, assessing the contribution of the variance in the value of each parameter to the variance in the model output. Both the individual parameter contribution (without interactions with the other parameters), and the total parameter contribution which includes interactions with all other parameters are calculated.

OAT results
Figures 1, 2 and 3 show the one-at-a-time (OAT) results for each parameter at each of the sites. These are generated by running HadOCC-GOTM with fifty different equally spaced values of each parameter across the parameter ranges shown in Table 3 while all other parameters remain fixed at their original value.
Looking at the results in Figs. 1, 2 and 3, variation in some parameters has a noticeable effect on the output value calculated, while variation in others does not. For CO 2 flux ( Fig. 1) individual parameter value variation does not have a huge effect -variation in some parameters (e.g. the maximum rate of photosynthesis P s max and the detrital sinking rate v s ) has a slight effect on the flux calculated of around 2-4 mol CO 2 m − 2 yr−1 change in the CO 2 flux across the whole range of the parameter, while other parameters (e.g. F messy , F nmp and F zmort ) have a negligible effect. In the case of primary production (Fig. 2) individual parameter variation has a much stronger effect, most notably in the case of the maximum rate of photosynthesis P s max and the initial slope of the photosynthesis-irradiance curve α (primary production changes by between 0.7 and 2.8 mg C m − 2 yr−1 across the range of P s max and α), as is expected given these parameters role in governing the ability of the phytoplankton to photosynthesise and grow (see Appendix A). For deep export (Fig. 3)  OAT plots for CO 2 flux (y-axis) in mol CO 2 m −2 yr −1 for each of the sites. CIS green, PAP blue and ESTOC red. The y-axis scale is the same in all plots for easy comparison of the relative effect of variation in the value of each parameter (x-axis) on the CO 2 flux. The effect of individual parameter variation on the CO 2 flux is generally small, with changes to the CO 2 flux greatest for the parameters P s max , α, g max , v s and θ.
www.jn.net Journalname Fig. 1. OAT plots for CO 2 flux (y-axis) in mol CO 2 m −2 yr −1 for each of the sites. CIS green, PAP blue and ESTOC red. The y-axis scale is the same in all plots for easy comparison of the relative effect of variation in the value of each parameter (x-axis) on the CO 2 flux. The effect of individual parameter variation on the CO 2 flux is generally small, with changes to the CO 2 flux greatest for the parameters P s max , α, g max , v s and θ . role in controlling export and breakdown of detritus in the water column.
Overall, two trends are notable in these OAT results. First, that the effects of individual parameter value variation on the outputs are nearly all monotonic (change only in one direction), the sole exception being the case of the carbon to chlorophyll ratio of phytoplankton θ in the calculation of the primary production (Fig. 2). Second, that while as expected the output values at the different sites are different, the trends seen in the outputs at each site due to parameter value variation are all similar -i.e. if increase in the value of a parameter causes increase in the output, this effect is seen for all sites.
The clearest examples of this are perhaps to be seen in Fig. 2 where increase in the value of P s max and α increases the primary production.
These results are readily interpreted. Looking at the model equations (see Appendix A), changes in individual parameter values are to be expected to have a monotonic effect on model outputs. Similarly, while the model experiences different forcing at the different sites its overall behaviour in response to parameter variations will stay the same as the mechanisms of the model remain unaltered. Table 4 details the results of the model runs using the parameter inputs generated by the Latin hypercube sampling of the parameter space (see Sect. 2.2). Comparing these results (in particular the minimum and maximum output values calculated) with the OAT results in Figs. 1, 2 and 3 shows that changing the values of multiple parameters has a greater potential effect on the calculated output than changing the value of only one parameter. This is most marked in the case of the CO 2 flux where OAT changes (Fig. 1) have little effect on the flux, while large variation is seen in the CO 2 flux results for the global SA, with in the case of the PAP site a change from a minimum of ≈ −5 (outgassing) to a maximum of ≈ 10 (ingassing) mol CO 2 m −2 yr −1 . In the case of primary production and deep export some parameter combinations enable almost no primary production (e.g. low values of both P s max and α) and export (e.g. low primary production and slow sinking rate). This indicates that parameter interactions play a significant role in the output calculations. OAT plots for primary production (y-axis) in mgC m −2 yr −1 for each of the sites. CIS green, PAP blue and ESTOC red. The y-axis scale is the same in all plots for easy comparison of the relative effect of variation in the value of each parameter (x-axis) on the primary production. Individual parameter variation has a large effect on primary production in many cases with the greatest effect seen for P s max and α.

Journalname
www.jn.net OAT plots for primary production (y-axis) in mg C m −2 yr −1 for each of the sites. CIS green, PAP blue and ESTOC red. The y-axis scale is the same in all plots for easy comparison of the relative effect of variation in the value of each parameter (x-axis) on the primary production. Individual parameter variation has a large effect on primary production in many cases with the greatest effect seen for P s max and α. Table 4. Results of GEM-SA model runs for net annual air-sea CO 2 flux (mol CO 2 m −2 yr −1 ), primary production (mg C m −2 yr −1 ) and annual export below maximum MLD (mg C m −2 yr −1 ). MADM is the median absolute deviation from the median scaled such that for a normal distribution it is equal to the standard deviation σ .

Discussion of GEM-SA results
The GEM-SA results show several general trends. The parameters that have a high sensitivity for each output are mostly the same at the different sites -the most obvious exception being θ which has a very high sensitivity at ESTOC but not at the CIS or PAP sites. By extension, parameters that have little sensitivity for the output are mostly the same for the different sites. This is similar to the trend seen for the OAT results in Figs. 1, 2 and 3. Parameters that have a large overall sensitivity (individual and interaction contributions), also mostly have a large individual contribution. Looking at the GEM-SA results in detail, interactions between each individual possible combination of any two parameters (190 possible pairs of parameters) nearly all contribute to much less than 1 % of the output variance, with the majority contributing less than 0.01 %. For each parameter, there are 19 possible interactions with the remaining 19 parameters. The sum of these 19 interactions is shown by the difference between the overall contribution of the parameter and the individual contribution (visible section of the narrow bars in Figs. 4, 5 and 6). This shows that, while the sum of the parameter interactions is important to explaining the output variance, in general, interactions between individual parameter pairs are not. Table 5 details all the parameter interactions (9 in total) that contribute greater than 2 % to the output variance. These are mostly interactions that might be expected. For primary production, P s max and α jointly define the photosynthesisirradiance curve and a significant interaction between these parameters is to be expected. Similarly, as v s controls the speed at which detritus sinks to the deep, and Rm deep controls the rate at which that detritus breaks down, deep export is likely to be influenced by interaction between these parameters.
We now look at the differences in parameter sensitivities for each of the outputs between the different sites. Starting with CO 2 flux, a total of ten different parameters, six at the CIS, five at ESTOC and seven at the PAP have an overall contribution of greater than 5 % to the output variance (see and exporting it away from the surface through mortality. Ecosystem processes that increase DIC fixing (e.g. phytoplankton growth) encourage further CO 2 uptake from the atmosphere, while ecosystem processes that release organic carbon back to DIC in the upper ocean (e.g. remineralisation and carbonate formation) restrict CO 2 uptake.  At the CIS, nutrients and sunlight are in relative abundance compared to that available at ESTOC and the PAP , so the value of the maximum rate of photosynthesis P s max has a very large influence on the CO 2 flux as it controls the amount of DIC taken up by phytoplankton growth at the surface, and hence the ability of the surface ocean to uptake CO 2 from the atmosphere. The initial slope of the photosynthesisirradiance curve α is also highly influential for the same reason.
At the PAP, while P s max and α remain sensitive parameters for the reason described above, the sinking rate of detritus v s is by far the most influential. As shown by the OAT results in Fig. 1, increasing values of v s increase the CO 2 flux into the ocean, as faster detrital sinking reduces the amount 414 V. Scott et al.: Ocean carbon cycle model SA of detritus that breaks down to dissolved inorganic carbon (DIC) near the surface. This appears to be the key process influencing the CO 2 flux at the PAP.
At the ESTOC, the carbon to chlorophyll ratio of phytoplankton θ is found to be by a large margin the most influential parameter (for CO 2 flux and primary production θ individual contribution to output variance >50 %). The ESTOC site exhibits oligotrophic characteristics with limited nutrient availability restricting growth. As a result, changes in the growth rate parameters (P s max and α) have little effect on the DIC uptake by the phytoplankton, but changing the ratio of carbon to chlorophyll in the phytoplankton will change the DIC uptake. Similarly, the rain ratio ϒ c which controls the release of organic carbon to DIC by carbonate formation is seen to be influential. However, looking at Table 4, while these parameters are influential in relative terms, the effect of parameter variation on the CO 2 flux at the ESTOC is relatively small compared to that seen at the other sites (a range of 1.33 mol CO 2 m −2 yr −1 ). At a site where biological activity is restricted by low nutrient availability, variation in the parameters governing the ecosystem behaviour do not have a large effect on the uptake of CO 2 flux.
For primary production ( Fig. 5) a similar story is seen. At the more productive CIS and PAP sites the phytoplankton growth parameters P s max and α are the most influential -this can also be seen in the OAT results in Fig. 2. At the less-productive ESTOC, changing the carbon to chlorophyll ratio of phytoplankton θ has a much greater effect as, for increasing θ, the same restricted supply of nutrient can produce a greater mass of phytoplankton. Similarly, phytoplankton specific mortality m 0 at the ESTOC is more influential as nutrient limited growth is less able to replace the dead phytoplankton.
Lastly, we consider deep export. As in the case of CO 2 flux and primary production, θ proves highly influential at the ES-TOC due to limited nutrient supply. Otherwise, as expected from the OAT results (Fig. 3), the most influential parameters are the detrital sinking rate v s , and the deep remineralisation rate Rm deep .

Discussion
The results of the OAT analysis and GEM-SA analysis are generally in agreement and as expected from the model structure (see Appendix A). The CO 2 flux is most sensitive to parameters that influence the CO 2 partial pressure (pCO 2 ) of the ocean surface, and hence the air-sea CO 2 exchange, by altering the DIC content of the sea-surface level. The phytoplankton growth parameters P s max and α influence the ocean surface pCO 2 through controlling how efficiently phytoplankton fixes dissolved inorganic carbon. At the well-lit surface, P s max is more important in limiting phytoplankton growth, and so has greater influence on the calculation of CO 2 flux than α. The sinking rate of detritus v s controls the rate at which detritus is removed from the sea-surface level directly influencing the surface ocean DIC concentration and hence its pCO 2 . The work of Schneider et al. (2008) using the more complex PISCES model, corroborates the importance of the parameterisation of particulate sinking on the calculation of the surface pCO 2 and resulting CO 2 flux. The rain ratio ϒ c directly effects the DIC concentration by setting the amount of organic carbon released to DIC through carbonate formation. The strong influence of the carbon to chlorophyll ratio θ at the ESTOC arises due to the relatively low nutrient availability at the site, as discussed in Sect. 3.3 above.
Primary production is sensitive to the parameters that control phytoplankton growth: α, P s max and m 0 (the phytoplankton-specific mortality rate), as has been found in previous studies such as that by Druon and Le Fèvre (1999). Deep export is found to be most sensitive to the detrital sinking rate v s , and the deep remineralisation rate Rm deep . Again, high sensitivity to the carbon to chlorophyll ratio θ at the ESTOC arises due to low levels of growth (see Sect. 3.3).
As can be seen from the literature resources used to identify the parameter ranges in Tables 1 and 2, these influential parameters are widely used in many other NPZD type ocean ecosystem models as they are well established representations of observed behaviour. The less well-established parameterisations and corresponding parameters particular to HadOCC such as the zooplankton feeding parameters F injest , F messy , F nmp and F zmort are mostly found to have weak influence. This is encouraging for the comparison of the results of different NPZD models, as it indicates that provided the parameterisations of phytoplankton growth and detrital export are similar, differences in the parameterisation of other processes may not be as critical.
While the more complex GEM-SA approach reveals some greater detail about the sensitivity of the parameters, its results agree closely with those seen in the OAT analysis. As seen in the GEM-SA results (Figs. 4,5 and 6) individually important parameters remain the most important when parameter interactions are explored. While parameter interactions remain important in the calculation of the model outputs, changes to the value of individual parameters have a greater overall effect. Overall, these results show that the relatively simple form of an NPZD model is robust to assumptions about its behaviour based on its equations, and that parameter sensitivity can be assessed using simple methods.
For use in GCMs, NPZD models like HadOCC are tuned to match measured bulk properties of the ocean ecosystem. While these sensitivity analysis results show that only certain widely used parameters may need to be tuned, they also indicate that changes to these parameters that might arise from the alteration of the marine environment by anthropogenic activity can have a large effect on fundamental ocean biogeochemical processes. This is most clearly seen in the OAT results shown in Figs. 1, 2 and 3 where changing highly sensitive parameters over a reasonable range of values is seen to have a large effect on the model outputs, particularly at the more productive CIS and PAP sites.
The ocean ecosystem is subject to change in its environment due to the anthropogenic decrease in ocean pH arising from increased CO 2 absorption by the ocean (Key et al., 2004;Orr et al., 2005), and the warming of the ocean (Doney, 2006). The effects of these environmental changes on the ocean ecosystem and its biogeochemical behaviour remain uncertain (Cao and Caldeira, 2008;Doney et al., 2009b). Having only one generic phytoplankton, zooplankton, nutrient and detritus compartment, controlled by a limited set of basic parameters, NPZD models such as HadOCC will not capture the effect of environmental change on the ecosystems behaviour as currently formulated. This is perhaps most clearly apparent in the use of a fixed value for the rain ratio ϒ c . Ocean acidification is likely to have a detrimental effect on carbonate producing phytoplankton, possibly changing the global rain ratio (Doney et al., 2009a), and this effect is not captured in HadOCC. Models that contain greater flexibility have capacity for response to environmental changes.
The development of more advanced and flexible ocean biogeochemistry models is focussed on the use of phytoplankton functional types (PFTs) each with their own associated parameters to separately represent the behaviours of different types of phytoplankton (Hood et al., 2006;Le Quéré et al., 2005). There is much debate on the development of PFT models (e.g., Anderson, 2005;Flynn, 2006;Hood et al., 2006;and Le Quéré, 2006) due to issues of poorly understood ecology and lack of data. Selection of PFT groups and parameterisation of their characteristics is subject to the limitations of current knowledge, but does enable exploration of the possible effects of environmental change on ocean ecosystem behaviour. The changes seen in the model outputs arising from changes in the values of "bulk property" parameters demonstrate that more flexible representations of the ocean ecosystem are needed to be able to predict the ocean ecosystem's response to climate change and potential feedbacks.

Summary and conclusion
We have shown that reasonable variation in the values of the biological parameters used in the HadOCC NPZD ocean biogeochemistry model have a large effect on the calculation of three fundamental outputs for biogeochemical modelling and climate prediction (air-sea CO 2 flux, export and primary production). The parameters of greatest importance are generally found to be those that are well-established representations of observed behaviour and are widely used in many NPZD models. While parameter interactions are found to influence the output calculation, changes to the values of individual parameters have a much greater effect than interaction with other parameters. The results of a simple one-at-a-time sensitivity analysis are not overturned by the results of a more complex global sensitivity analysis.
Parameters that control phytoplankton growth (the maximum photosynthetic rate P s max and the initial slope of the photosynthesis-irradiance curve α), and the parameters that control the rate of sinking of detritus and the formation of carbonate (the rain ratio) have the greatest influence on the calculation of the air-sea CO 2 flux. The calculation of the primary production is most influenced by the values of the phytoplankton growth parameters P s max and α, and the specific mortality rate of the phytoplankton. Deep export (export to below the maximum annual mixed layer depth) is most influenced by the values of the sinking rate of detritus, and the remineralisation rate below the upper ocean. In a nutrient depleted oligiotrophic site, the ratio of carbon to chlorophyll of phytoplankton θ is the most influential parameter, but changes in its (and other parameter values) have much less absolute effect on the calculated outputs than in more productive regions.
The parameters of NPZD models used in GCMs are tuned to the large-scale properties of the ocean ecosystem. These large scale properties have the potential to change due to anthropogenic modification of the marine environment (Siegel and Franz, 2010). While NPZD models can be successfully tuned to reproduce the measured large-scale ocean ecosystem properties, their most influential parameters are not formulated to respond to large-scale changes in the marine environment. This potentially limits the power of NPZD models in GCMs to provide insight into the role of the ocean ecosystem in future climates.

Model equations
The NPZD model equations are described below. For more detail on these please refer to the original paper on HadOCC (Palmer and Totterdell, 2001). The model parameters are listed in Table 1.
where M P is the natural phytoplankton mortality. This is assumed to result entirely from viral infection, so the specific rate increases with population density such that M P = mP 2 where m = m 0 . However, if the phytoplankton concentration falls below the threshold of P ≤ 0.01 mMol N m −3 the natural mortality is switched off (m = 0) as the population is assumed to be too small to transmit infection. The nutrient limitation on phytoplankton concentration is given by The light limitation I lim is estimated using the dailyaveraged version of the spectrally-averaged parameterisation by Anderson (1993) such that where α max = 2.602α, the derivation of α # is detailed by Anderson (1993) and I is photosynthetically active radiation (W m −2 ).

Zooplankton graze on both phytoplankton and detritus.
If h is the grazing rate per unit food concentration then the losses to phytoplankton and detritus are H P = hP and H D = hD respectively, where and F = max(0,F tot − F th ), where F tot = B P P + B D D, F th = 0.1 and K F is the half saturation constant for grazing. The factors B P = 14.01+12.01 P 14.01+12.01 Red = 1 (where Red is the Redfield ratio), B Z = 14.01+12.01 Z 14.01+12.01 Red and B D = 14.01+12.01 D 14.01+12.01 Red are used to adjust for the different nitrogen content per unit biomass in zooplankton, phytoplankton and detritus. ∂C ∂t = P F nmp M P + P ηP + Z F zmort (A8) Appendix B

Model drivers: nutrient supply and meteorology
The replenishment of nutrient in the upper ocean by mixing depends partly on lateral advection from upwelling regions and estuarine regions. To represent this additional input from outside the 1-D GOTM column, the nutrient in HadOCC-GOTM is relaxed below the productive depth to a nutrient profile taken from Levitus et al. (1993). The productive depth is defined to be the greatest depth at which both light and nutrient are available in sufficient quantities for photosynthesis to take place -i.e. the shallower of the mixed layer depth or euphotic depth. The modelled nutrient profile (N) is relaxed to the fixed Levitus et al. (1993) profile (N Lev ) by reducing the difference between them by 1/60 per days worth of timesteps ( t d ) such that the replenished nutrient HadOCC-GOTM is forced with meteorological data for air pressure, wind speed, relative humidity, air temperature and total cloud cover at six hourly intervals. This data is taken from the ECMWF 40-year Re-analysis (ERA-40) dataset (Uppala et al., 2005). These are the meteorological processes that influence the rate of air-sea CO 2 transfer (wind, humidity, pressure and temperature and the amount of light (cloud cover) that reaches the ocean surface and is hence available for phytoplankton photosynthesis. The model uses the standard heat flux method in GOTM (Burchard et al., 1999).