An intermediate complexity marine ecosystem model for the global domain

A new marine ecosystem model designed for the global domain is presented, and model output is compared with field data from nine different locations. Field data were collected as part of the international Joint Global Ocean Flux Study (JGOFS) program, and from historical time series stations. The field data include a wide variety of marine ecosystem types, including nitrogen- and iron-limited systems, and different physical environments from high latitudes to the mid-ocean gyres. Model output is generally in good agreement with field data from these diverse ecosystems. These results imply that the ecosystem model presented here can be reliably applied over the global domain. The model includes multiple potentially limiting nutrients that regulate phytoplankton growth rates. There are three phytoplankton classes, diatoms, diazotrophs, and a generic small phytoplankton class. Growth rates can be limited by available nitrogen, phosphorus, iron, and/or light levels. The diatoms can also be limited by silicon. The diazotrophs are capable of nitrogen fixation of N2 gas and cannot be nitrogen-limited. Calcification by phytoplankton is parameterized as a variable fraction of primary production by the small phytoplankton group. There is one zooplankton class that grazes the three phytoplankton groups and a large detrital pool. The large detrital pool sinks out of the mixed layer, while a smaller detrital pool, representing dissolved organic matter and very small particulates, does not sink. Remineralization of the detrital pools is parameterized with a temperature-dependent function. We explicitly model the dissolved iron cycle in marine surface waters including inputs of iron from subsurface sources and from atmospheric dust deposition. © 2001 Elsevier Science Ltd. All rights reserved.


Introduction
The decade-long series of field studies carried out as part of the Joint Global Ocean Flux Study (JGOFS) program has provided modelers with an unparalleled data set across a wide variety of marine ecosystem types. Currently, there are serious deficiencies with attempts to model biogeochemical cycles across diverse regimes; to incorporate multi-element nutrient limitation, trophic structure, and phytoplankton functional groups; and to predict marine ecosystem responses to climate change (Doney, 1999). The JGOFS data set in conjunction with other field programs and newly available satellite data (such as surface chlorophyll data from SeaWiFS, the Sea viewing Wide Field of View Sensor; McClain et al., 1998) should allow major improvements in global-scale modeling efforts (Doney, 1999;Evans, 1999). Here, we present a new marine ecosystem model designed for the global domain and compare model output with field data from nine JGOFS study sites. The model is an adaptation of the ecosystem model of Doney et al. (1996), expanded to allow for multielement limitation of phytoplankton growth (Fe, N, P, or Si) and to include three phytoplankton classes. The goal of this work is to develop an intermediate complexity marine ecosystem model suitable for modeling ocean biogeochemistry within global climate system models. In a companion paper, we present global fields of model output and examine iron cycling in surface waters of the world ocean, the sensitivity of the marine system to atmospheric iron deposition, and global-scale nutrient-limitation patterns (Moore et al., 2002, hereafter METB).
Over the last decade there has been an ever-growing awareness of the key role that the micronutrient iron plays in marine ecosystems. Bottle incubation experiments, field and satellite observations, and in situ fertilization experiments, such as the IronEx and SOIREE experiments, have conclusively shown that iron availability limits phytoplankton growth rates and determines ecosystem structure over large areas of the world ocean (Helbling et al., 1991;Martin et al., 1991;Martin, 1992;Price et al., 1994;de Baar et al., 1995;Coale et al., 1996;Landry et al., 1997;Gordon et al., 1998;Lindley and Barber, 1998;Behrenfeld and Kolber, 1999;Boyd and Harrison, 1999;Boyd et al., 2000).
The consistent pattern observed in these studies is that while the large and small phytoplankton species are iron-stressed in the high nitrate, low chlorophyll (HNLC) regions, it is the biomass of the larger diatoms that increases when available iron is increased. The smaller phytoplankton are efficiently grazed by microzooplankton, preventing large increases in biomass. Thus, the biomass of the smaller phytoplankton species is regulated by moderate iron limitation and strong grazing pressure, while diatom biomass is regulated largely by available dissolved iron concentrations (Price et al., 1994;Landry et al., 1997).  examined how iron availability determines ecosystem structure and strongly influences total and export production patterns in the Southern Ocean.
Global-scale ocean biogeochemical models must incorporate the effects of iron limitation on phytoplankton growth rates and on ecosystem structure to reliably capture the processes controlling total and export production. In recent years, iron limitation has been incorporated into marine ecosystem models in both simple and more complex ways Loukos et al., 1997;Leonard et al., 1999;Denman and Pe * na, 1999;Lancelot et al., 2000;Christian et al., 2002). Here we explicitly model phytoplankton iron limitation and the dissolved iron cycle in surface waters, including the atmospheric deposition of dust as an iron source.
It is well established that a strong preference for ammonium over nitrate in the iron-stressed HNLC regions is partially the reason for the observed high nitrate concentrations (Wheeler and Kokkinakis, 1990;Price et al., 1994;Armstrong, 1999a and references therein). There is good reason for a stronger ammonium preference in these regions, as valuable iron is required to reduce nitrate to a usable form (Armstrong, 1999a). Armstrong (1999a) examined the interactions between iron, ammonium, and nitrate use with a model that partitioned available iron in the cell between nitrate assimilation and photosynthetic carbon-related demands. Flynn and Hipkin (1999) also examined interactions between iron, light, ammonium, and nitrate use in a phytoplankton physiological model. They found an increased ammonium preference under iron stress, but they suggested the pattern might not be universal (Flynn and Hipkin, 1999). We incorporate this influence of iron stress on nitrogen uptake patterns, increasing the preference for ammonium over nitrate under iron-limiting conditions with a relatively simple parameterization.
Iron stress also has been observed to influence the relative uptake rates of silicate and nitrate by diatoms. Several recent papers have demonstrated that the ratio of silicate uptake to nitrate uptake increases with iron stress, from a value of B1 (mol/mol) under iron-replete conditions to a value of B2-3 (mol/mol) under iron-limiting conditions (Takeda, 1998;. One possible factor in this pattern is the stronger ammonium preference under iron limitation discussed above. However,  argue using several lines of reasoning that Fe-limited diatoms are more heavily silicified, that is, cellular Si/N and Si/C ratios are higher when cells are iron-stressed. Likely, both processes are important. We have incorporated an iron effect on silicon and nitrogen use by diatoms in the ecological model, such that the maximum silica cell quota (and thus Si/N and Si/C ratios) increases under iron stress.
In a recent paper, Fung et al. (2000) estimated 2-D spatial maps of the iron budget for the upper ocean. They compared atmospheric iron inputs with subsurface sources due to upwelling/mixing processes. The subsurface iron field was based on an extrapolation of the Moss Landing dissolved iron data set (Johnson et al., 1997). The Moss Landing data set was used to assign regional dissolved iron/nitrate ratios below the mixed layer (Fung et al., 2000). We adopt this approach to estimate subsurface iron concentrations with some modifications. We compare global scale iron budgets with the results from Fung et al. (2000) in METB.
In recent years, there has been a growing awareness that several key functional groups of phytoplankton play important roles in the oceanic carbon cycle (i.e., Doney, 1999). We have incoporated two of these groups explicitly in the model, the diatoms and the diazotrophs, and have implicitly included a third group, the coccolithophores. Diatoms are often a dominant component of the export flux from surface waters. Diazotrophs are phytoplankton capable of nitrogen fixation of N 2 gas into bioavailable nitrogen forms. We have modeled the diazotrophs based on available information for Trichodesmium spp., the dominant diazotroph in the world ocean. Diazotrophs represent a source of new nitrogen to surface waters in tropical regions and may influence primary/export production rates and whether an ecosystem is N-or P-limited (Karl et al., 1997;Falkowski, 1997;Capone et al., 1997). Another important functional group of phytoplankton are the coccolithophores, which incorporate CaCO 3 into their cells as outer platelets (coccoliths). This calcification can have significant effects on surface water carbon chemistry, pCO 2 concentrations, and the air-sea carbon flux (Holligan et al., 1993;Robertson et al., 1994). The bloom forming Emiliania huxleyi is one of the dominant coccolithophores globally (Holligan et al., 1993;Balch et al., 1992. We have parameterized calcification by phytoplankton in the ecosystem model by assuming that a variable percentage of the small phytoplankton group is coccolithophores. The percentage of coccolithophores varies as a function of ambient temperature and nutrient conditions and is also dependent on biomass levels of the small phytoplankton class. The model presented here is an eleven-component mixed-layer ecosystem model. We compare model output with data from nine JGOFS study sites: the Bermuda Atlantic Time Series (BATS), the Hawaiian Ocean Time series (HOT), the Subarctic North Pacific site at Station P (STNP), the Arabian Sea (ARAB), the Equatorial Pacific (EQPAC), the North Atlantic Bloom Experiment (NABE), the Kerfix time series (KERFIX), the Ross Sea (ROSS), and the Antarctic Polar Front region (PFR).
The ecosystem model has been developed for incorporation into the NCAR 3-D ocean circulation model. However, the results presented here are from a simplified global mixed-layer grid that captures the local processes of turbulent mixing, vertical velocity at the base of the mixed layer, and seasonal mixed-layer entrainment and detrainment. Horizontal advection is not included. Thus, the model is run independently at each grid point with no lateral exchange. Other forcings include sea-surface temperature, percent sea-ice cover, surface radiation, and the atmospheric deposition of iron and silicon. For the implementation used here we neglect vertical structure and focus strictly on mean values for the surface mixed layer.

Ecosystem model
The ecosystem model is adapted from Doney et al. (1996) and consists of 11 main compartments, zooplankton, small phytoplankton, diatoms, diazotrophs, two detrital classes, and the nutrients nitrate, ammonium, phosphate, iron, and silicate. The small phytoplankton size class is meant to represent nano-and pico-sized phytoplankton, and may be Fe-, N-, P-, and/or light-limited. The larger phytoplankton class is explicitly modeled as diatoms and may be limited by the above factors as well as Si. Growth rates of the diazotrophs may be limited by Fe, P, and/or light. Calcium carbonate production by phytoplankton is parameterized as a variable percentage of primary production by the small phytoplankton group. Many of the biotic and detrital compartments contain multiple elemental pools as we track carbon, nitrogen, phosphorus, iron, silicon, and calcium carbonate through the ecosystem. A schematic of the model including the elemental fields is shown in Fig. 1. A full listing of model terms, parameters and equations is provided in Appendix A. Computer code for the global mixed-layer grid and the full ecosystem model are available through the US JGOFS Synthesis and Modeling Project website (http://usjgofs.whoi.edu/mzweb/syn-mod.htm).
Phytoplankton growth rates are determined by available light and nutrients using a modified form of the growth model of Geider et al. (1998, hereafter GD98). The GD98 model was modified to allow for two nitrogen sources (ammonium and nitrate) and to allow for growth limitation by the other nutrients. Iron, phosphorus, and silicon limitation are modeled in the same manner as nitrogen in GD98. Phytoplankton cell quotas of Chl, N, P, Fe, and Si are measured relative to C ratios. The maximum and minimum iron quotas we have used encompass the mean value of 5.0 mmol/mol estimated by Johnson et al. (1997). However, these Fe/C ratios are generally low compared to coastal phytoplankton species and even some open ocean species (i.e. Sunda and Huntsman, 1995). The species that thrive in the HNLC areas will have adapted to low iron concentrations, and we thus feel that these are reasonable iron quotas for use in a global model. They are also in very good agreement with an iron requirement, B8-fold lower than the diazotrophs, and the available evidence regarding diazotroph iron quotas (see Discussion and References below). We examine the model sensitivity to higher maximum Fe cell quotas in METB.
If phytoplankton are unable to attain their maximum cell quota through uptake, carbon fixation (growth rate) is reduced proportionately (Relations (106), (107), (111), (112), (116) and 117)). Thus, if the iron quota (the cellular Fe/C ratio) is at half of its maximum value, the maximum rate of carbon fixation (growth rate) is reduced by half. Whichever nutrient is currently most-limiting (lowest cell quota relative to the maximum quotas) modifies the carbon fixation rate. There is good laboratory evidence for a relationship between cell quotas (measured as nutrient/C ratios) and specific growth rates (Sunda and Huntsman, 1995;Geider et al., 1998). For example, Sunda (1997) noted that based on laboratory measurements at a Fe/C ratio of 2 mmol/ mol would reduce the specific growth rate of the diatom Thalassiosira oceanica by 65%. Using the quota scenario outlined here, growth rate would be reduced by 71% at an Fe/C ratio of 2 mmol/ mol. Geider et al. (1998) present evidence for a similar effect due to changing N/C ratios. Some of the carbon fixed by phytoplankton (proportional to nitrogen uptake) is used for biosynthesis and thus is not added to the phytoplankton carbon pool (lambda term in Relations (110), (115) and (119), see GD98). This biosynthesis cost is higher when cells are growing on nitrate than when ammonium is the N source (Geider, 1992). We thus modify this term based on the current F-ratio to account for the lower cost of growth on ammonium (Relations (109) and (114)).
There are two detrital pools in the model, a non-sinking pool, which largely represents dissolved organic matter (DOM) but would also include colloidal and small non-sinking particulate organic matter, and a large particulate detrital pool, which represents particulate organic matter (POM) and sinks at a rate of 20.0 m/day (Appendix A). Biomass incorporated into small phytoplankton and diazotroph pools is largely recycled within the mixed layer through the small detrital pool, while biomass incorporated into the diatom pool is largely exported through the sinking largedetrital pool (see following text and Appendix A), an approach similar to that of Laws et al. (2000). We model only the labile/semi-labile fractions for both detrital pools. Both detrital pools are remineralized at a temperature dependent rate with a maximum rate of 10%/day at 301C (Relations (41), (42) and (155)-(164)). The nutrients nitrogen, phosphorus, and iron all remineralize at this temperature dependent rate. Carbon remineralizes at a slightly slower rate (95%) than the nutrients. Similarly remineralization of silica and calcium carbonate is set at even slower rates (50% and 1% of the value for the other nutrient pools, respectively). The length scale of remineralization for the sinking POM pool varies from B200-1000 m depending on temperature. The time scale for the remineralization of nutrients in the DOM pool ranges from B10-50 days depending on temperature. Length and time scales for remineralization of C, Si, and CaCO 3 would be proportionally longer.
Nitrogen uptake is partitioned between ammonium and nitrate according to the substitutive equations of O'Neill et al. (1989) based on Fasham (1995). These equations are based on Michaelis-Menten kinetics modified by the half-saturation coefficients (Relations (62), (63), (70) and (71)). They allow for increasing ammonium preference as ambient ammonium levels increase, without reducing total nitrogen uptake. We have set the half-saturation constants (Ks) for ammonium relatively low, and those for nitrate relatively high to give phytoplankton a reasonably strong preference for ammonium (Appendix A). Nitrification occurs only at low light levels corresponding to winter conditions at high latitudes (radiation o4.0 W À2 averaged over the mixed layer (Relations (48) and (49)). The small phytoplankton have lower half-saturation constants than the larger diatoms for all nutrients, while the diazotroph halfsaturation constants are intermediate (2 times that of the small phytoplankton for phosphate) (Appendix A). Diazotrophs (more specifically Trichodesmium spp.) often occur as large colonies (i.e. Karl 1996, 1998). Surface to volume considerations imply higher half-saturation constants for these colonies relative to the small phytoplankton class. It is also known that Trichodesmium spp. are not very efficient at inorganic nutrient uptake (i.e., McCarthy and Carpenter, 1979).
Under iron stress, we increase the preference for ammonium over nitrate by modifying the halfsaturation constants for N uptake (Relations (52)-(57)). Half-saturation values for nitrate increase with increasing iron stress (decreasing Fe cell quota) up to a maximum of 150% of their original value. Likewise, half-saturation constants for ammonium decrease with increasing iron stress to a minimum of one-half the original value. These modifications of half-saturation constants in conjunction relative uptake equations are a simple way to capture the observed stronger ammonium preference in iron-limited systems. It is likely that phytoplankton in ironlimited systems adapt to become more efficient at ammonium uptake even in the presence of relatively high nitrate concentrations (see Armstrong, 1999a). Our parameterization attempts to include this iron effect on ammonium preference while maintaining a relatively constant total nitrogen uptake.
Relative uptake rates of Fe, P, and Si are modeled using Michaelis-Menten kinetics (Appendix A). Fitzwater et al. (1996) estimated a community half-saturation constant for iron of 0.12 nM in the equatorial Pacific, which includes uptake by the diatoms whose biomass increased with iron addition and smaller cells which were efficiently grazed. The small phytoplankton species that dominate in HNLC regions would be expected to have a lower half-saturation constant than the larger bloom forming diatoms. Price et al. (1994) estimated halfsaturation constants of 0.035 nM for natural populations at the equator and of 0.22 nM at 151N. We use a relatively low value for the iron half-saturation constant of 0.08 nM for the diazotrophs and small phytoplankton, and a higher value of 0.2 nM for the larger diatoms. We assume that the diazotrophs have adapted to relatively efficient iron uptake due to their larger iron demands. In addition, Trichodesmium colonies may be able to extract additional iron from dust particles (Reuter et al., 1992), a process not accounted for in our model.
Half-saturation constants for silicate uptake can vary considerably in different regions (Nelson and Tr! eguer, 1992). Brzezinski and Nelson (1989) found values as low as 0.5 mM within Gulf Stream warm-core rings. Nelson and Tr! eguer (1992) found values of 1.1-4.6 mM in the Southern Ocean ice-edge zone. Much higher values have sometimes been reported in Southern Ocean waters (i.e. Nelson et al., 2001). Half-saturation Si constants for growth are lower than the value for Si uptake (Nelson et al., 2001). Here we adopt a value of 1.2 mM for a globally applicable value (Appendix A). We examine the sensitivity to a higher half-saturation constant for silicate uptake in METB. For all nutrients, uptake rates are near their maximum value when cell quotas are low, and decline towards zero as the maximum cell quota is approached (GD98). We have modified the maximum uptake equation from GD98 to give a somewhat more gradual taper off of uptake rates as cell quota approaches its maximum value (i.e. see Relation (78)).
Photoadaptation is modeled according to Geider et al. (1996Geider et al. ( , 1998 (Appendix A, Relations (121)-(126). Chlorophyll synthesis is regulated by the balance between photosynthetic carbon fixation and light absorption (the ratio of energy assimilated to energy absorbed) (Geider et al., 1996(Geider et al., , 1998. Chlorophyll synthesis is assumed to be proportional to nitrogen uptake, reflecting the need for the synthesis of proteins used in light harvesting complexes and elsewhere in the photosynthetic system (GD98). A maximum Chl/N ratio is one of the input parameters to the model (Appendix A). We have used a value of 3.0 mg Chl/mmol N.
To model diazotrophs in our ecosystem model we have drawn largely on recent models efforts for Trichodesmium spp. (Fennel et al., 2002;Hood et al., 2001) as well as field data, laboratory studies, and global distributions of Trichodesmium spp. Karl, 1996, 1998;Karl et al., 1997;Capone et al., 1997;Kustka et al., 2002). It is known from these studies (and references therein) that Trichodesmium thrives in warm tropical/subtropical waters under well-stratified, high light conditions. Trichodesmium is typically observed in waters warmer than 201C and blooms usually occur in waters warmer than 251C (Capone et al., 1997). We modify the maximum growth rate for the diazotroph as a function of temperature (as with the other phytoplankton) and assume that production and biomass of the diazotrophs are negligible when sea-surface temperatures are below 161C (Relations (120) and (151)). Trichodesmium have very slow natural growth rates, with doubling times typically between 3 and 5 days (Capone et al., 1997). The known affinity of Trichodesmium for relatively high light levels is thought to be a function of high respiration costs and a low initial slope of the production vs. irradiance curve (alpha) (see Fennel et al., 2002, and references therein). We have set the non-grazing mortality (which includes respiration) at 15% per day for the diazotrophs compared with 10% per day for the other phytoplankton. We also use a much lower alpha value and a much slower maximum growth rate for the diazotrophs compared with the other phytoplankton groups (Appendix A). Thus, there is a strong light constraint on diazotroph growth in our model such that high biomass levels can only occur during prolonged periods of relatively shallow mixed-layer depths. Fennel et al. (2002) assumed respiration losses of 17% per day in addition to a second quadratic mortality term. Little is known about the grazing losses of Trichodesmium, Fennel et al. (2002) did not include an explicit grazing loss term, while Hood et al. (2001) set a relatively high grazing loss rate similar to the other phytoplankton but a much lower respiration loss. We have set the maximum grazing rate much lower for the diazotrophs than for the other phytoplankton (Appendix A). This allows for blooms of the diazotrophs under optimum environmental conditions. Trichodesmium spp. have been shown to exude up to 50% of the nitrogen fixed from N 2 gas as dissolved organic nitrogen (DON) (Capone et al., 1994;Glibert and Bronk, 1994). We assume that 30% of fixed N is released as DON (enters the small detrital pool) and that enough N is fixed to maintain near maximum N cell quotas for the diazotrophs despite this loss (Relations (76)-(83)). Furthermore, following Fennel et al. (2002) we assume that all N requirements for the diazotrophs are met by fixing N 2 gas and that the supply of N 2 is unlimited. While Trichodesmium are capable of uptake of nitrate, ammonium, and amino acids, it is thought that this uptake is a small fraction of N fixation in the field. Orcutt et al. (2001) found that uptake of ammonium and amino acids during summer months at BATS was typically o1% of nitrogen fixation rates (although some significant uptake of these other N forms was observed during winter and early spring months, when biomass was generally low).
The elemental composition of diazotrophs differs from other phytoplankton which tend to have the traditional Redfield ratios.  reported a mean N/P ratio for Trichodesmium of 45 mol/mol and a mean C/N molar ratio of 6.3 close to the Redfield value of 6.6. Mean Trichodesmium C/N ratios at BATS are also approximately 6.6 mol/mol (Orcutt et al., 2001). Fennel et al. (2002) set a fixed P/N value of 45 mol/mol for Trichodesmium in their model. We set the optimum C/N ratio at near Redfield values for the diazotrophs in the same manner as the other phytoplankton (i.e. the same minimum and maximum N/C ratios, see Appendix A). We assume an optimum N/P molar ratio for the diazotrophs of 45.0 based on the above references but let this ratio vary like all the other model elemental ratios. In practice this is accomplished by setting the minimum and maximum P/C ratios lower for the diazotrophs than the other phytoplankton (Appendix A).
It has been argued that iron requirements for diazotrophs are much higher than for other phytoplankton due to the iron requirements of N fixation (Raven, 1988;Reuter et al., 1992). Kustka et al. (2002) recently estimated that the iron requirement for Trichodesmium is B8 times that of eukaryotic phytoplankton based on theoretical calculations in conjunction with field measurements. Based on this work we assume that the iron requirements of the diazotrophs are 8 times higher than the other phytoplankton and thus the minimum and maximum Fe/C quotas of diazotrophs are increased to 8.0 and 56.0 mmol/mol (Appendix A). These iron requirements are in good agreement with recent laboratory experiments where Trichodesmium growth rates and nitrogen fixation rates approached maximum values at cellular Fe/C ratios of B50 mmol/mol (Ilana Berman-Frank, Jay Cullen, Yaeila Hareli and Paul Falkowski, personal comm.). In addition, Orcutt et al. (2001) noted that molar N : Fe ratios for Trichodesmium colonies in the North Atlantic approach a value of 3200 during winter months and increase during summer months. Assuming a Redfield C/N ratio (Orcutt et al., 2001) for these colonies, this translates to Fe/C ratio of 47 mmol/mol. Molar N/Fe ratios at BATS during summer months range from 346-5337 ( Fig. 4 of Orcutt et al., 2001). This translates to a broad range of Fe/C ratios between 28 and 438 mmol/mol. Dust particles stuck to colonies as described by Reuter et al. (1992) could have resulted in an overestimation of iron content in some samples. There was no significant correlation between iron content and nitrogen fixation rates (Orcutt et al., 2001), implying that the Trichodesmium colonies were generally iron-replete even at the wintertime Fe/C ratios of B47 mmol/mol.
To parameterize calcification rates by the phytoplankton we made several assumptions based on the known distributions of calcification in the oceans as described by Milliman (1993). Calcification rates are quite low at the highest latitudes (lowest sea-surface temperatures), and also generally lower in the mid-ocean gyres. Rates are higher in mid-latitude temperate waters and in coastal upwelling systems (Milliman, 1993). Blooms of E. huxleyi are most common in temperate waters, particularly in the North Atlantic (i.e. Holligan et al., 1993). In the equatorial Pacific, Balch and Kilpatrick (1996) reported calcification rates that were highly variable but on average about 10% of photosynthetic carbon fixation. Balch et al. (2000) found that coccolithophore calcification was typically between 1% and 5% of total photosynthesis rates in the Arabian Sea.
Calcification rates in E. huxleyi vary as a function of environmental factors such as light levels and available nutrients, and also with growth phase (Balch et al., 1992. Calcification by phytoplankton is a complex process. Our calcification parameterization is an admittedly crude attempt to improve on previous modeling efforts which often assumed a constant calcification/ production ratio, as in the Ocean Carbon Model Intercomparison Project II (OCMIP) (Doney et al., 2002a;Najjar et al., 2002).
In our parameterization of calcification we begin with a base calcification rate equal to 5% of the photosynthetic carbon fixation of the small phytoplankton group (Relation (127)). This base rate is then modified by several factors. The base rate is reduced as a function of the current minimum nutrient quota squared (Relation (128)). Thus, calcification rates are lower under extreme nutrient limitation. Under these conditions (as is often the case in the mid-ocean gyres) we assume that the coccolithophores would be largely outcompeted by smaller picoplankton who would have an advantage due to their larger cell surface/volume ratios. The base calcification rate is similarly progressively reduced as sea-surface temperatures fall below 51C, and calcification is assumed to be negligible at sea-surface temperatures below 01C (Relations (129)-(130)). Lastly we progressively increase the calcification/photosynthesis ratio when the small phytoplankton class blooms (i.e. when biomass exceeds 2.0 mmol C/m 3 , Relation (131)). This is justified because coccolithophores are frequently a dominant component of non-diatom phytoplankton blooms.
Based on Armstrong (1999b) we have only one zooplankton compartment, which grazes on the large detrital pool and all three phytoplankton classes. This ''single'' grazer system is more stable than one with multiple classes of zooplankton (Armstrong, 1999b;Lima et al., 2002), and encompasses the actions of both the microzooplankton and larger zooplankton in a highly parameterized fashion. Maximum growth rates for the zooplankton varies with the food source. When small phytoplankton are the food source (microzooplankton grazing) the maximum growth rate is set higher than when diatoms, large detritus or diazotrophs are the food source (larger zooplankton grazing) (Relations (132)-(135)). Thus, under near optimum growth conditions, the diatoms are able to escape grazing control and bloom until nutrients become limiting. It is more difficult for the small phytoplankton and the diazotrophs to bloom, but notably, not impossible for either group.
The grazing equations contain a quadratic density dependent term such that grazing rates decline at low prey biomass. Lessard and Murrell (1998) recently found evidence for such a threshold effect on grazing working with natural microzooplankton assemblages in the Sargasso Sea. Grazing rates declined sharply at chlorophyll concentrations o0.1 mg/m 3 (Lessard and Murrell, 1998). Fitting curves to the experimentally determined grazing rates, they suggested the threshold where grazing rates approached zero was at chlorophyll concentrations of B0.024-0.035 mg/m 3 . Since chlorophyll values at or below this level are rarely seen in situ, even in the most oligotrophic gyre regions, they suggested that this grazing threshold effect may be an important factor in determining the consistently low phytoplankton biomass levels seen in many regions of the world ocean (Lessard and Murrell, 1998). The grazing threshold for the diatoms and the large detritus is set lower than for the small phytoplankton and the diazotrophs (Relations (132)-(135)). The lower thresholds are justified as the large zooplankton that feed on diatoms, and large detrital particles are typically capable of mobility and of actively filtering the water in search of prey.
Zooplankton incorporate 30% of grazed matter into new zooplankton biomass (Straile, 1997). The remaining 70% of the grazed material (termed z slop, see Appendix A) is lost to sloppy feeding, remineralization within the gut, to excretion, or to faecal pellets and is routed in the following manner. Half of this 70% (35% of total) grazed phytoplankton C, N, Fe, P, Si, and CaCO 3 is remineralized by grazing processes and enters directly the appropriate nutrient pool. Remineralized C and CaCO 3 are not currently tracked in the model. The other half enters the small detrital pool (for grazing on small phytoplankton) or the large detrital pool (for grazing on diatoms and large detritus). The half not remineralized is split evenly between the detrital pools when the diazotrophs are the food source. The CaCO 3 associated with the small phytoplankton class that is not remineralized (50%, based on Harris, 1994) is added to the large detrital pool. All grazed silica not remineralized is added to the large detrital pool (65% of total as none is incorporated by zooplankton).
There are three types of phytoplankton mortality/loss in the model. In addition to grazing losses there are losses due to a non-grazing mortality and due to aggregation/sinking. The nongrazing mortality is set at a constant value of 10% per day for the diatoms and small phytoplankton and at 15% per day for the diazotrophs. This term would include losses due to viral lysis (Fuhrman, 1992), as well as internal respiration/degradation, and excretion. The GD98 model included a term Rref to account for internal respiration/degradation losses. We have set Rref equal to 0.0, (as did GD98) but retain the term in our equations for completeness (Appendix A). These losses are instead accounted for in our non-grazing mortality term. Cell aggregation losses for the diatoms and small phytoplankton are according to quadratic equations (Relations (141)-(149)). Thus, aggregation losses are minimal when biomass is low but become significant under bloom conditions (Alldredge and Jackson, 1995, and references therein). A minimum aggregation loss of 5% per day is set for the diatoms to account for direct sinking losses. A maximum aggregation loss of 70% per day is imposed on both phytoplankton groups. There is no aggregation loss for the diazotrophs, as we assume that the diazotrophs cannot sink out of surface waters due to their gas vacuoles (see Letelier and Karl, 1998).
Aggregation losses for both phytoplankton classes enter the sinking detrital pool (Appendix A). Non-grazing mortality losses for the diatoms are routed in the following manner. Seventy-five percent of the C, N, Fe, and P enters the small detrital pool, representing excretion and release of dissolved and small particulate matter. The remaining 25% (plus all of the silica) enters the large detrital pool, representing organic matter bound to the frustule and larger particulate matter. Most of the small phytoplankton and all of the diazotroph non-grazing mortality biomass enters the small detrital pool. A typically small fraction of the small phytoplankton biomass associated with the CaCO 3 enters the POM pool. The fraction routed to the POM pool is equal to 25% of the C, N, Fe, and P associated with CaCO 3 , assuming a mean cellular inorganic C/organic C ratio of 0.5 (Appendix A, Relations (27)-(36)). This ratio can vary considerably as a function of growth stage and other environmental factors (Balch et al., 1992(Balch et al., , 1993. Balch et al. (1992) estimated a calcification/organic carbon production ratio of B0.37 in a bloom of E. huxleyi in the Gulf of Maine. Other phytoplankton contribute to the organic carbon production so this value is likely an underestimate for the coccolithophores themselves (Balch et al., 1992).  found that daily integrated calcification/photosynthesis ratios varied from 0.2 to 0.7 in lightlimited cultures of E. huxleyi. Balch et al. (1992) found ratios greater than 1.0 for logarithmic growth phase cultures of E. huxleyi.
We set a minimum phytoplankton biomass (0.001 mmol C/m 3 , for the small phytoplankton) using the Pprime term where the non-grazing and aggregation mortality loss terms are set equal to zero (Relations (141)- (152)). The minimum threshold is needed to preserve viable seed populations of diatoms and small phytoplankton through the winter at high latitudes, and to ensure a viable seed population for the slow growing diazotrophs in tropical waters. The minimum biomass is set at 0.005 mmol C/m 3 for the diatoms and at 0.03 mmol C/m 3 for the diazotrophs. A higher threshold for the diatoms is justified because some diatoms can overwinter within the sea ice and thus begin the growing season at relatively high biomass levels. Other diatoms form resting cysts that also can provide seed stock in the spring. The diazotrophs have much slower maximum growth rates than the other phytoplankton (see Appendix A), and a higher minimum biomass is needed to ensure a viable seed population throughout tropical/ subtropical areas. Fennel et al. (2002) also set a relatively high background biomass for Trichodesmium in their ecosystem model (0.01 mmol N/m 3 ). The minimum biomass for diazotrophs is set much lower (0.001 mmol C/m 3 ) in waters where sea-surface temperature is less than 16.01C (Relation (151)). We assume negligible biomass and production by diazotrophs in these cold waters. In the strongly nutrient-limited mid-ocean gyre regions with a typical biomass of B0.5 mmol C/m 3 , the mortality of the small phytoplankton group is reduced by only B0.2% by the use of the Pprime term. For the diatoms in these areas, at a typical biomass of 0.1 mmol C/ m 3 , mortality is reduced by B5% by the use of the Pprime term. We also set a minimum biomass for zooplankton (0.01 mmol C/m 3 ) to retain a seed population through the winter at high latitudes (Relations (153) and (154)).
Zooplankton mortality closes the ecosystem model at the upper end of the food chain. Zooplankton morality consists of a linear (6% per day) and a quadratic biomass-dependent morality loss (see Appendix A parameter list and Relations (153) and (154)), which parameterizes zooplankton losses due to respiration and predation by higher trophic levels (Steele and Henderson, 1992). These loss terms effectively set a strong cap on zooplankton biomass and when the phytoplankton escape grazing control, large blooms occur if there is sufficient light and nutrients available. Upon mortality zooplankton biomass is partitioned between the sinking and non-sinking detrital pools based on food source. A fraction of zooplankton mortality ( f zoo detr, Relation (136)) consisting of 80% of the current grazing on large detritus and diatoms plus 30% of the current grazing on small phytoplankton plus 50% of the grazing on diazotrophs divided by the total current grazing is routed to the sinking detrital pool (Relations (27)-(30)). The remainder (1.0Àf zoo detr) enters the non-sinking detrital pool (Relations (33)-(36)). Thus, most of the effectively ''larger'' zooplankton enter the sinking pool upon death, while a smaller fraction of the microzooplankton grazing is passed up the food chain and enters the sinking detrital pool.
There are still many unanswered questions concerning the cycling of iron in marine surface waters, including the solubility of iron in dust particles, the portion that is bioavailable, conversion between dissolved and particulate forms, rates of biological uptake and particle scavenging, and just what are the absolute dissolved iron concentrations (i.e. Bruland et al., 1994;Wells et al., 1995;Johnson et al., 1997;Measures and Vink, 1999;Fung et al., 2000;Jickells and Spokes, 2002). We have incorporated current understanding and theories to model the dissolved iron cycle in surface waters. This necessarily involved a number of assumptions and simplifications. Undoubtedly the model will need to be refined as our understanding of the marine iron cycle improves.
We assume that all of the dissolved iron is bioavailable and that particulate iron is not bioavailable (and therefore neglected in the model). Little is known about abiotic scavenging rates for iron. It is known that at low iron concentrations (oB0.6 nM) most (>99%) of the dissolved iron is bound to organic ligands, and are not strongly particle reactive (van den Berg, 1995;Bruland, 1995, 1997). Bruland et al. (1994) estimated a lifetime for iron in the deep ocean of 70-140 years, a relatively short period. Johnson et al. (1997) argued that scavenging was negligible due to ligand binding at iron lower concentrations (oB0.6 nM). However, this would result in a substantially longer lifetime than 70-140 years. Lef" evre and Watson (1999) and Archer and Johnson (2000) examine several scavenging scenarios including with and without ligand ''protection''.
We have adopted a simple particle scavenging loss term for dissolved iron. We assume that at low iron concentrations (oB0.6 nM) most dissolved iron is bound to organic ligands, but is still weakly particle reactive, and that at iron concentrations above 0.6 nM, particle scavenging rates increase rapidly with increasing iron concentrations. Thus, at iron concentrations below 0.6 nM we include a weak scavenging loss for dissolved iron of 1% per year (Relation (50)). This corresponds to a 100-year lifetime neglecting biological uptake. As iron concentrations exceed this threshold, scavenging increases with increasing iron concentrations according to Michaelis-Menton type kinetics up to a maximum rate of 2.74% per day with a half-saturation constant of 2.5 nM (Relation (51)). There are several reasons for increasing the scavenging loss rate at higher iron concentrations. As iron exceeds the concentrations of organic ligands the dissolved iron would become more particle reactive. As iron concentrations increase, precipitation to form iron hydroxides also would result in a loss of dissolved iron. In the model, iron only exceeds B2.0 nM in areas of very high atmospheric dust deposition. In these regions the portion of iron that enters the dissolved phase may be lower than in other regions (Jickells and Spokes, 2002). In addition there would be more potentially scavenging particles present in the these high dust areas. Our increased scavenging rate accounts for these effects. Scavenged iron is added to the large detrital pool and may be remineralized during sinking.

Model grid and forcings
The model is run on a 2-D global grid (100 Â 116 gridpoints) that corresponds to the top layer of the coarse resolution NCAR Climate Ocean Model (NCOM) (Large et al., 1997;Gent et al., 1998;Doney et al., 1998Doney et al., , 2002b. Longitudinal resolution is 3.61 and latitudinal resolution varies from 1-21, with higher resolution near the equator. There is no lateral advection or exchange between grid points; thus the model is run independently at each gridpoint. Model forcings at each location include surface shortwave radiation, sea-surface temperature, mixed-layer depth, the deposition of atmospheric iron and silicon, vertical velocity at the base of the mixed layer, turbulent mixing at the base of the mixed layer, and percent sea-ice cover (Appendix A). The model is run for 2 years as spin-up, and then monthly data are summed and output from a third year.
Surface shortwave radiation is from the ISCCP cloud-cover-corrected data set (Bishop and Rossow, 1991;Rossow and Schiffer, 1991). Photosynthetically available radiation (PAR) is assumed to be 45% of the surface shortwave flux, and is averaged over the mixed layer using an absorption coefficient, which depends on water absorption and chlorophyll concentration (Relations (37) and (38)). Sea-surface temperature data are from the World Ocean Atlas 1998 (hereafter WOA98) data set (Conkright et al., 1998) remapped to the NCOM grid. Climatological mixed-layer depths (based on the 0.125 sigmaÀy difference criteria) are used to force the model (Monterey and Levitus, 1997). Mixed-layer depths were modified so that the minimum mixed-layer depth is 25 m. This was done so that a significant fraction of total primary production would be accounted for in the model and is roughly consistent with the wind-mixed-layer depth. At several of the JGOFS sites, the climatological mixed-layer depths were further modified, either to remove anomalous points (i.e. at KERFIX the August mixedlayer depth was o5 m), or to make mixed-layer depths more similar to the observed data at the JGOFS sites during the period when the in situ data with which we compare was collected (see Fig. 3 for in situ measurements, original climatology, and modified climatological mixed-layer depths).
Vertical velocity at the base of the mixed layer was derived from a run of the NCAR 3-D ocean model and then smoothed with a 3 Â 3 boxcar filter to account somewhat for the effects of lateral advection. Thus instead of very strong upwelling at the equator combined with lateral advection, we have strong upwelling on the equator and moderate upwelling in adjacent regions. Downwelling has no effect on mixed-layer concentrations, while upwelling is assumed to displace an equal volume of water laterally. We do not track this laterally exported water in the model, and it is thus a form of ''export'' from the surface mixed layer. This is a weakness of the simple global grid used here. Physical dynamics in areas where lateral advection is important, such as at the equator and in coastal upwelling regions will be better represented in a 3-D ocean model implementation. Turbulent mixing at the base of the mixed layer is set at a constant value of 0.15 m/day. This is somewhat higher than the value of 0.1 m/day used in many modeling studies, and has been increased to account for nutrient inputs from physical processes not included in the model (i.e. storm passage, eddy events, etc. y).
We have used the dust deposition model study of Fung (1994, 1995) to estimate inputs of iron from the atmosphere. Iron was assumed to comprise 3.5% (by weight) of atmospheric dust based on a recent study of mineral dust by Zhu et al. (1997). This is also the fraction of iron within dust assumed by Duce and Tindale (1991). We further assume that 2% of the deposited iron enters the dissolved iron pool upon deposition, based on the review of Jickells and Spokes (2002, see also Fung et al., 2000). In the companion paper, we examine how varying total dust/iron flux and the soluble iron fraction affects marine ecosystem dynamics at the global scale (METB). We also include atmospheric deposition of Si by assuming a constant Si/Fe ratio in the dust. Based on Duce and Tindale (1991) we assume that atmospheric dust is a constant 30.8% Si by weight and that 7.5% of this Si is soluble and thus enters our dissolved silicate pool upon deposition.
Global monthly fields of percent sea-ice cover (NASA Team algorithm) from the DMSP-F13 satellite sensor were obtained from the EOSDIS NSIDC Distributed Active Archive Center for the period October 1998-September 1999(National Snow and Ice Center, 2000. Sea-ice cover acts to reduce short wave radiation levels at the sea surface in direct proportion to %coverage (i.e. 50% sea-ice cover reduces surface radiation flux by 50%). Monthly fields of each forcing parameter are input to the model, and a cubic spline curve fit to the monthly data are then used to interpolate the forcings at each time step.

Boundary conditions
For each of the nutrients nitrate, phosphate, silicate, and iron, the bottom-boundary condition (or concentration below the mixed layer) was computed as a function of a surface intercept value and a slope value that increases linearly with the depth of the mixed layer. The slope was computed as the deep-water value (value at maximum mixed-layer depth or 200 m depth, whichever was shallower)Fthe surface value, divided by maximum mixed-layer depth (or 200 m). We assume that nutrients are not seasonally depleted by the biology below 200 m. Nutrient concentrations below the mixed layer were then computed as the surface value+(slope * mixedlayer depth) at each timestep in the model. For all other ecosystem variables, the concentration just below the mixed layer was assumed to be similar to mixed-layer values at shallow mixed-layer depths (75% of mixed-layer values at 25 m mixed-layer depth), and to decrease linearly with increasing mixed-layer depth to a value of 0.0 at 100 m.
Bottom-boundary conditions for silicate, phosphate, and nitrate were derived from WOA98 (Conkright et al., 1998). The surface intercept value was determined using summer season data (for each hemisphere), which represents a depleted annual minimum value. These surface fields were smoothed with a 5 Â 5 gridpoints boxcar filter to remove small-scale variability. The concentration at maximum mixed-layer depth (or 200 m) also was taken from WOA98 (annual fields). The surface intercept value for iron used in computing the bottom boundary condition was set at 50 pM based on Johnson et al. (1997). Mean surface concentration was somewhat higher in that study (68 pM), but because many of their samples were below detection limits, Johnson et al. (1997) suggested the mean was overestimated. We have used the iron/nitrate ratios estimated by Fung et al. (2000) in conjunction with the WOA98 nitrate data to determine iron concentrations at the maximum mixed-layer depth (or at 200 m) in most areas. We modified the iron/nitrate ratios of Fung et al. (2000) in the North Atlantic where we set a higher iron/nitrate ratio of 40 (mmol/mol) and in coastal waters (defined as areas where depth o600 m on our coarse resolution global grid, rather than areas within 200 km of the coastline) where we set the slope such that iron concentrations reach 1.05 nM at 100 m depth. Of the JGOFS locations we compare here, only the ROSS site falls into this ''coastal'' category. In some cases the surface nutrient value (iron, phosphate, nitrate, or silicate) exceeded the deep-water concentration, in these instances the surface value was reduced to 70% of the deep-water value. Maximum bottom-boundary conditions for each of these four nutrients were also imposed so that unrealistically high nutrient levels would not be introduced where mixed-layer depths exceeded 200 m (32.0 mM for nitrate, 120 mM for silicate, 2.0 mM for phosphate, and 2.0 nM for iron).
We increased the subsurface iron concentrations in the North Atlantic for several reasons. Fung et al. (2000) found in their analysis that the North Atlantic should be an iron-stressed region comparable to the Subarctic North Pacific. Yet field measurements and bottle incubation studies suggest that this region is not strongly iron-limited (Martin et al., 1993). The high Fe/C ratios estimated by Sunda (1997) for the North Atlantic also argue that this region is relatively ironreplete. Fung et al. (2000) suggested that iron inputs from coastal upwelling, rivers, and shallow shelf areas (all not included in their analysis) may be important in this region. Likewise, Martin et al. (1993) and Luther and Wu (1997) suggest that some of the iron requirements in this region are likely met by lateral transport from the continental margins. Therefore, we have increased subsurface iron concentrations to account for this advective source of iron.

JGOFS and satellite data sets
We compare data collected during the nine JGOFS studies with model output from the corresponding grid point on the global grid described above. Our approach has been to average available mixed-layer measurements to compare with the model output. One exception is dissolved iron concentrations where in situ data were quite sparse. We plot all iron measurements made at depths of 50 m or less. For primary production data, we integrate available measurements (typically 24 h incubations) over the mixed layer, and then divide by mixed-layer depth to estimate mean production per cubic meter.
Location of the JGOFS sites are shown over annual mean chlorophyll concentrations from SeaWiFS in Fig. 2. Note the wide diversity in ecosystems types among the different stations. HOT and BATS are located in the mid-ocean gyres with typically low chlorophyll concentrations. Several high-latitude sites experience deep winter mixing (NABE, KERFIX, PFR, and ROSS) and attain moderate to high chlorophyll concentrations, at least seasonally (Fig. 2). The STNP site has perpetually moderate chlorophyll concentrations. While the ARAB site has moderate chlorophyll concentrations most of the year, it is influenced seasonally by coastal upwelling associated with the Southwest Monsoon in late summer.
In some of the JGOFS studies, transects were made over large areas and across different ecosystem types. Thus, we focus on specific areas from these studies. For the Arabian sea study we have extracted data from station S7 of the 1994-1996 Arabian Sea Expedition . The US Southern Ocean JGOFS study, also known as Antarctic Environment and Southern Ocean Process Study (AESOPS), focused on two regions, the Ross Sea continental shelf region (ROSS), and the second was in the vicinity of the Antarctic Polar Front (PF) region along 1701W (PFR) . For the ROSS site we compare with all data from the transects along B761S. Due to the latitudinal gradients in nutrient concentrations across the PF, we restrict our comparison to waters at the PF and just to the south (168-1721W, 60.5-62.51S, see Smith et al., 2000). At the NABE site, we compare with available in situ data from the region 17-211W by 46-481N. At all other sites we compare with data collected within 11 of latitude and longitude of the locations given in the figures. Ocean Station Papa (STNP) at 1451W, 501N was part of the recent Canadian JGOFS study in the NE subarctic Pacific (see Boyd and Harrison, 1999). We will compare with some literature values from this study; however, the in situ nutrient and chlorophyll data we have plotted for comparison are from earlier historical data collected at STNP. The in situ data from the recent Canadian JGOFS study were not yet available. Surface chlorophyll measurements (as opposed to mixed-layer averages at other sites) are plotted from the period (1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976) (Anderson et al., 1977;Stephens, 1977). Historical nutrient data from the period (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) also are compared with model output from this site.
Mixed-layer depth criteria for the in situ data differ at some sites, with most regions using the criteria of a density difference from the surface of 0.125 (sigmaÀy units). At all three Southern Ocean sites a sigmaÀy difference from the surface of 0.02 was used. The entire mixed-layer data set along with the details of data sources and our processing is available through the US JGOFS Synthesis and Modeling Project website (http://usjgofs.whoi.edu/mzweb/syn-mod.htm) and as a NCAR Technical Report (Kleypas and Doney, 2001).
We also compare model estimates of chlorophyll concentration with SeaWiFS satellite-based estimates for each JGOFS site. Level-3 standard mapped images of global surface chlorophyll (Version 2.0) covering the period October 1998 through September 1999 (post-El Ni * no period) were obtained from the Goddard Distributed Active Archive Center (McClain et al., 1998). The chlorophyll data were log-transformed before averaging and then remapped to our global grid. We compare model predicted dissolved iron concentrations with data from the Moss Landing database (Johnson et al., 1997), and with more recent data collected as part of the JGOFS studies (Measures and Vink (1999) for the ARAB site; and measurements from AESOPS for the ROSS and PFR sites, obtained from the US JGOFS database; Fitzwater et al., 2000;Measures and Vink, 2001;Coale et al., 1996).

Results
Climatological mixed-layer depths from the Monterey and Levitus (1997) data set are compared with in situ measurements from the JGOFS sites in Fig. 3. Also, shown in Fig. 3 are the modified climatological mixed-layer depths used in this study. Modifications to the climatological depths were made at the ARAB site (to more closely match the deep mixing periods) and at the ROSS, PFR, and KERFIX sites to more closely match the observed spring stratification in the in situ data (Fig. 3). We wanted the model forcings to be similar to the in situ data to enable a more direct comparison. The seasonal cycles are in good agreement at most sites. The climatological mixed-layer depths never exceed 60 m at HOT; thus, nutrient inputs at this site are quite low. The EQPAC in situ mixed-layer depths are quite variable over short timescales (diurnal to several days). This short-term variability is not present in our monthly forcing data.
Seasonal patterns of phytoplankton and zooplankton biomass (carbon units) are shown in Fig. 4. In the spring, there is a mixed assemblage of diatoms and smaller phytoplankton at BATS, transitioning to a small phytoplankton dominated regime during the summer. A similar pattern exists at HOT all year (Fig. 4). In subsequent sections, we will show that strong nitrogen limitation of the diatoms is driving this pattern. At both BATS and HOT the biomass of the diazotrophs increases during summer months and peaks in late summer (Fig. 4). This is consistent with field observations at both sites (Karl et al., 1997;Orcutt et al., 2001). There is a low, relatively constant diazotroph biomass at EQPAC (B0.04 mmol C/m 3 ) and moderate biomass during the March-July period at ARAB (peaking at B0.13 mmol C/m 3 during May) (Fig. 4). Diazotroph biomass at the other sites is negligible all year (Fig. 4). At the high latitude HNLC Southern Ocean sites (PFR and KERFIX), the small phytoplankton dominate the assemblage with only modest diatom biomass (Fig. 4). We will show that the pattern is driven by strong iron limitation of the diatoms. At STNP, the small phytoplankton dominate the assemblage during winter and spring months, and a mixed assemblage of small phytoplankton and diatoms is seen in late summer and fall. At EQPAC, biomass remains relatively low and constant all year (Fig. 4). Subsequent figures will show that the diatoms experience substantial iron stress at EQPAC and Fig. 3. Measurements of mixed-layer depths from the nine JGOFS sites. In most cases a 0.125 sigmaÀy difference from the surface was the criteria used to determine mixed-layer depth (see text for details). Climatological mixed-layer depths are shown as squares (Monterey and Levitus, 1997). Modified climatological mixed-layer depths used in this study are shown as triangles.
that modest iron stress and strong grazing pressure maintain low biomass levels for the small phytoplankton. At NABE, there is a strong spring bloom initially of diatoms, then by the small phytoplankton (Fig. 4). There is also a strong phytoplankton bloom in the Ross Sea from December through February and diatom blooms at the ARAB site during fall and winter months (Figs. 3 and 4). At all of the mid-to high-latitude locations zooplankton biomass peaks during summer months (Fig. 4). We will examine the processes driving these seasonal patterns in more detail in the context of the following figures.
There are relatively few estimates of carbon biomass in the scientific literature of the components as shown in Fig. 4. Caron et al. (1995) estimated the particulate carbon in different classes of microorganisms from two cruises (spring March-April, and late summer in August) in the Sargasso Sea. Mean phytoplankton biomass was estimated as 0.92 and 0.43 mmol C/m 3 for the spring and summer cruises, respectively (from Tables 4 and 6 of Caron et al., 1995). Adding the model small phytoplankton, diazotroph, and diatom biomass from BATS gives comparable values of 0.71 mmol C/m 3 for March, 0.74 mmol C/m 3 for April, and 0.83 mmol C/m 3 during August (Fig. 4). Similarly, microzooplankton biomass (including heterotrophic nano-and microplankton) was estimated to be 0.42 and 0.35 mmol C/m 3 for the spring and summer cruises (from Tables 4 and 6 of Caron et al., 1995). The model estimates for BATS zooplankton biomass are 0.40 mmol C/m 3 for March, 0.39 mmol C/m 3 for April, and 0.41 mmol C/m 3 for August. Considering the uncertainties in the field estimate conversion factors (see Caron et al., 1995) the agreement between the model and observations is reasonable. Roman et al. (1995) estimated total zooplankton biomass in surface waters during EQPAC in October to be B0.43 mmol C/m 3 during the day and B0.63 mmol C/m 3 during the nighttime (their Fig. 6). The model estimates for EQPAC zooplankton biomass are B0.48 mmol C/m 3 , in good agreement with the field measurements ( Fig. 4).
At the ARAB site, Garrison et al. (1998) estimated carbon biomass within several ecosystem components during the SW Monsoon period (August-September). At station S7 mixed-layer autotrophic nanoplankton biomass was estimated to be 0.47 mmol C/m 3 and autotrophic microplankton biomass (mostly diatoms) was 1.2 mmol C/m 3 (Garrison et al., 1998). Model estimates for this period for diatom biomass were similar at 0.59 mmol C/m 3 during August and 3.7 mmol C/m 3 during September (Fig. 4). The model biomass of the small phytoplankton plus the diazotrophs was 0.37 mmol C/m 3 during August and 0.40 mmol C/m 3 during September in good agreement with the field estimates. Model zooplankton biomass was 0.63 mmol C/m 3 during August and 0.84 mmol C/m 3 during September. Garrison et al. (1998) estimated heterotrophic dinoflagellate abundance at station S7 to be >2.0 mg C/l (and o2.3 mg C/l) and heterotrophic nanoflagellates to be >1.4 mg C/l (and o2.0 mg C/l) in surface waters, and ciliate abundance at >3 mg C/l (their Fig. 2). If mean ciliate abundance was 4 mg C/l, dinoflagellate was 2.2 mg C/l and nanoflagellate biomass were 1.7 mg C/l, this would total to B0.66 mmol C/m 3 , a number close to our total zooplankton biomass.
Model estimates of surface mixed-layer nitrate concentrations are compared with the in situ data in Fig. 5. The seasonal patterns are an excellent fit to the in situ data at most of the JGOFS sites. Nitrate is fully depleted at HOT and seasonally depleted at BATS, ARAB, and NABE (Fig. 5). At the EQPAC site, it must be recalled the spring in situ data (February-April) were collected during El Ni * no conditions (weak upwelling), while our model is only forced with the La Ni * na (stronger upwelling) conditions. Model nitrate is depleted at the ARAB site by the winter diatom bloom (Figs. 4 and 5). At the ROSS site, initial spring nitrate concentrations are too low, but the seasonal drawdown by the biology is of a similar magnitude to the in situ data (Fig. 5). The model spring bloom at the PFR site does not deplete nitrate as strongly as observed during AESOPS (Fig. 5).
Ammonium concentrations from the model run are compared with in situ data from the JGOFS sites in Fig. 6. Ammonium concentrations are quite low all year at HOT and BATS (subsequent figures will show strong nitrogen limitation at these sites). At the high-latitude locations, ammonium concentrations increase at the end of the summer and in post-bloom conditions when heterotrophic processes exceed autotrophic ones (seen in the model output and in situ data, Figs. 4 and 6). Wheeler and Kokkinakis (1990) did not find any ammonium concentrations above 0.4 mM during the spring/summer in the Subarctic North Pacific, with most values below 0.25 mM. Similarly, Varela and Harrison (1999) measured ammonium concentrations between 0.1 and 0.5 mg-at/l at Ocean Station P. Model-predicted ammonium concentrations are too low at PFR during austral fall and too high at ARAB in January and February when compared with the in situ data (Fig. 6).
If we do not include the iron effect on ammonium preference (Relations (54)-(57)), maximum ammonium concentrations increase at KERFIX, PFR, EQPAC, and STNP by roughly a factor of two (not shown). Subsequent figures will show these regions to be iron-limited. In the other regions ammonium concentrations are only modestly different from the standard run shown in Fig. 6. Seasonal patterns of phytoplankton biomass and nitrate concentration were similar in the alternate run without the iron effect on ammonium preference. Modeled f-ratios were higher at KERFIX, PFR, and STNP without the iron effect on ammonium preference, and little changed elsewhere. Thus, the main result from including this iron effect is to lower ambient ammonium concentrations with a modest decrease in f-ratios (B10-20%) in the iron-limited regions.
Monthly averaged total phytoplankton f-ratios are displayed in Fig. 7. Square symbols indicate f-ratios calculated as nitrate uptake divided by total nitrogen uptake by the small phytoplankton and diatoms; triangle symbols indicate the f-ratio calculated as total nitrate uptake plus total nitrogen fixation divided by total N uptake by all phytoplankton. Comparing the two f-ratios it can be seen that nitrogen fixation accounts for more than 50% of new nitrogen during summer months at HOT and roughly 50% at EQPAC (Fig. 7). The data at HOT are in good agreement with field estimates that nitrogen fixation accounts for at least 27% (and up to 50%) of new production at HOT on an annual basis Karl et al., 1997). We are unaware of any nitrogen fixation data from the equatorial Pacific, and it seems likely that diazotroph biomass and nitrogen fixation are overestimated by the model in this region (even though total N fixation is relatively low, see Fig. 17). The model predicts total f-ratios (>0.4) at the EQPAC site (Fig. 7). Much of this is due to the influence of nitrogen fixation, but even ignoring this factor, f-ratios are too high (Fig. 7), as observations would suggest f-ratios of o0.12  (Landry et al., 1997). Nitrogen fixation is also significant in terms of new production at the ARAB site and BATS during summer months and has little effect at the other sites (Fig. 7). The lowest f-ratios of o0.1 are seen seasonally at HOT and at BATS, while the highest values are seen at the ROSS site (Fig. 7). At NABE, the f-ratio is 0.3-0.45 during the spring bloom and then declines during summer months. Martin et al. (1993) estimated an f-ratio of 0.45 during the spring bloom period at NABE based on carbon fluxes. Despite high available nitrate concentrations, f-ratios remain relatively low in the iron-stressed regions of the subarctic North Pacific (STNP), the equatorial Pacific (EQPAC), and in the Southern Ocean (PFR and KERFIX). In contrast, high f-ratios are seen in the Ross Sea, where iron is less limiting (Fig. 7). In the Southern Ocean, high f-ratios (>0.5) are associated with phytoplankton bloom periods (most often seen in areas of elevated iron inputs, i.e. coastal waters and following the retreating ice edge), while low productivity open-ocean waters typically have low f-ratios (o0.3) . Waldron et al. (1995) found f-ratios up to 0.6 in a front-associated phytoplankton bloom, and lower values (0.28) away from the front in the Bellingshausen Sea. Values in excess of 0.85 have been reported during spring in the Weddell Sea (Kristiansen et al., 1992). In the subarctic North Pacific, Varela and Harrison (1999) found depth-integrated f-ratios to range from 0.05 to 0.37 (mean of 0.21), which is in good agreement with the model results (Fig. 7). In the Arabian Sea, McCarthy et al. (1999) found mean f-ratios of 0.15 (November) and 0.13 (January) during the NE Monsoon period; the model predicted f-ratios at the ARAB site are considerably higher, especially when N fixation is included (Fig. 7).
Model-predicted phosphate concentrations are shown in comparison with in situ data in Fig. 8. At most sites, the seasonal drawdown is of comparable magnitude between the model and in situ data. There are offsets at HOT and the ROSS site (Fig. 8). Our coarse resolution grid likely overestimates (underestimates) the subsurface concentrations at these sites. The magnitude of the seasonal drawdown is in good agreement with the in situ data at ROSS (Fig. 8). Again, the spring drawdown at PFR is lower than was observed during AESOPS (Fig. 8). Phosphate is strongly depleted at BATS and NABE during summer months (Fig. 8).
Model-predicted silicate concentrations are compared with the in situ data in Fig. 9. Seasonal patterns are generally in good agreement with the in situ data. There are offsets at the HOT and KERFIX sites. Both locations lie within strong gradients in silicate in the WOA98 climatologies and it is likely that averaging within our coarse resolution grid overestimates the subsurface source. The seasonal drawdown at both sites is of similar magnitude in the model and in situ data (Fig. 9). The drawdown at the PFR site is much weaker than was observed during AESOPS. The spring bloom at NABE depletes silicate from surface waters, a phenomenon also observed in the Fig. 7. Monthly f-ratio values from the model defined as nitrate uptake divided by (nitrate+ammonium uptake (squares; include only the small phytoplankton and diatoms) and as (nitrate uptake+nitrogen fixation) divided by total nitrogen uptake (triangles; all phytoplankton groups).
in situ data (Fig. 9). Silicate concentrations are low at ARAB during winter due to strong depletion by the diatom bloom (Figs. 4 and 9). In a separate model run (not shown), we did not include the iron effect on diatom Si quotas (Relation (58)). Without the iron effect, ambient silicate concentrations are increased by B1.0 mM at EQPAC and by several mM during summer months at the STNP and ROSS sites. Silicate concentrations also decreased substantially during summer months at ARAB as there was a midsummer diatom bloom during June-July not seen in our standard run (subsequent figures will show this site to be silicon limited). Despite strong iron limitation at the PFR and KERFIX sites (see following figures and Discussion), silicate concentrations in the two runs are quite similar. This is because low ambient iron concentrations limit diatom production to low levels at both sites (see Fig. 4, and following sections). Seasonal patterns of diatom biomass were nearly identical in the alternate and standard model runs except at the ARAB site (standard run data shown in Fig. 4). The significantly lower ambient silicate concentrations at EQPAC when the iron effect on silicate uptake is included, imply that there could be silicon or iron/silicon co-limitation in this region, particularly off the equator away from the strong upwelling zone (and in other areas with only moderately high ambient silicate concentrations, like the Arabian Sea and Subantarctic waters) (see Dugdale and Wilkerson, 1998). We return to this topic of co-limitation in METB.
Dissolved-iron concentrations are shown in Fig. 10. The highest model iron concentrations are seen at the ARAB site (large atmospheric dust inputs) and at the ROSS site (coastal subsurface source), while the lowest average iron concentrations are seen at EQPAC, KERFIX, and PFR (Fig. 10). Model iron concentrations are also relatively low at HOT, STNP, and seasonally at NABE. Summertime diazotroph production results in marked drawdown of dissolved iron at HOT (Figs. 4 and 10). In Fig. 10, we compare with in situ data from the Moss Landing group (shown as asterisks, from the Johnson et al., 1997 database and AESOPS measurements, Fitzwater et al., 2000;Coale et al., 1996; and shown as diamonds from Measures and Vink (1999, for the ARAB site; and AESOPS measurements at PFR, Measures and Vink, 2001). It is apparent that there is a discrepancy between the measurements at the PFR site. The two groups use different methods for measuring dissolved iron (see Measures et al., 1995;Johnson et al., 1997), and at the PFR site Measures and Vink had more near surface measurements where iron concentrations may have been higher due to inputs from melting sea ice (see Discussion section). The discrepancy highlights the need for iron standard solutions used for calibration and for intercomparison studies among groups measuring dissolved iron. Since our subsurface iron boundary conditions are based on measurements by the Moss Landing group, our model predicted concentrations are, not surprisingly, closer to their measurements. The detection limit for the Moss Landing group was 0.03-0.04 nM; thus iron concentrations plotted at this level actually may have been lower.
Iron is strongly depleted by mid-summer at all of the Southern Ocean sites in the model output (Fig. 10). Iron concentrations are also quite low at NABE during summer months. Deep mixing during winter months increases iron concentrations at all of the high-latitude sites (Fig. 10). Note that our increase of the dissolved iron/nitrate ratio in the subsurface waters of the North Atlantic does not lead to unrealistically high iron concentrations being introduced by winter convection. Values as high as 0.6 nM have been measured in North Atlantic surface waters Luther, 1994, 1996). Martin et al. (1993) measured deep-water iron concentrations at the NABE site of 0.36 nM (600 m) and 0.4 nM (700 m). Given persistent winter mixed-layer depths of 500 m in the model (and in situ data, see Fig. 3), our springtime iron concentration of B0.6 nM is not unreasonably high. The ROSS site begins the growing season with high dissolved-iron concentrations but is strongly depleted by the phytoplankton bloom, reaching very low concentrations by January, in good agreement with the in situ measurements (Fig. 10).
In an alternate model run (not shown) without any scavenging loss for dissolved iron in the surface layer, the seasonal patterns were very similar to those shown in Fig. 10. Peak ambient iron concentrations in the low iron areas (all sites except ARAB and ROSS) were little affected with no scavenging loss (differences less than 1%). Maximum mean monthly dissolved iron concentrations at ARAB and ROSS were approximately 2% and 7% higher, respectively, without the scavenging loss. At the ROSS site, dissolved iron concentrations remained slightly above 1.2 nM through Fig. 10. Model mixed-layer dissolved iron concentrations compared with in situ measurements from the JGOFS sites. Star symbols are from Johnson et al. (1997), and from Fitzwater et al. (2000) and from unpublished AESOPS data for the PFR. Diamond symbols are from Measures and Vink (1999) at the ARAB site and Measures and Vink (2001) at the PFR site.
November but were still depleted by January. Seasonal patterns of phytoplankton biomass at all sites were nearly identical in the alternate run, with no iron scavenging except at the ARAB site, where a mid-summer diatom bloom developed during June-July (not shown). This pattern was quite similar to the bloom that developed when the iron effect on cellular Si quotas was not included.
Next we examine how these available nutrient concentrations (Figs. 5,6,8,9 and 10) affect phytoplankton growth rates (by influencing nutrient cell quotas). In Figs. 11 and 12, monthly mean relative cell nutrient quotas for the diatoms and small phytoplankton are displayed. Cell quota relative to the maximum cell quota is shown for all nutrients. Phytoplankton growth is limited by whatever nutrient is most limiting (lowest cell quota relative to the maximum quota, see Methods section and Appendix A). Phytoplankton are relatively nutrient-replete (at or near maximum cell quotas) during winter months at the high latitude sites (Figs. 11 and 12). Recall that cell quotas are measured relative to cellular carbon in the model. During winter, when carbon fixation (growth rates) are minimal due to strong light limitation, it is easy to meet the cellular quota even if uptake rates are very low. We can now return to the seasonal patterns in phytoplankton biomass displayed in Fig. 4 and explain them in terms of light and nutrient limitation, as well as grazing pressure from the zooplankton. Fig. 11. Relative nutrient cell quotas of the diatoms. Cell quota relative to the maximum cell quota is displayed (see text for details). Whichever nutrient is most-limiting (lowest relative cell quota) limits phytoplankton growth rate. Nutrient pools depicted include nitrogen (square), iron (triangle), phosphorus (diamond) and silica (star).
At BATS, deep mixed layers during the winter lead to low light levels within the mixed layer (Fig. 4). Low carbon fixation due to light limitation, combined with nutrient entrainment from below, allows both size classes of phytoplankton to attain near maximum cell nutrient quotas, although the diatoms are still slightly nitrogen-stressed (Figs. 11 and 12). In the summer, both phytoplankton classes become nitrogen limited, with the diatoms more strongly limited due to their higher half-saturation constants for nitrate and ammonium uptake and the very low ambient nitrate and ammonium concentrations (Figs. 5, 6, 11 and 12). Thus, small phytoplankton dominate the assemblage at BATS during the summer (Fig. 4). The diatoms also experience some phosphorus stress during summer months at BATS (Fig. 11). A similar dominance of the assemblage by the small phytoplankton is observed at HOT all year, where a lack of deep mixing in the winter strongly limits nutrient inputs to the surface mixed layer (Figs. 3, 4, 11 and 12). Note that nitrogen stress is decreased for both the diatoms and the small phytoplankton at HOT during late summer/fall months (Fig. 12). This is due to new nitrogen being added to the ecosystem due to nitrogen fixation by the diazotrophs (Figs. 4 and 17).
At STNP, EQPAC, KERFIX, and PFR, a common pattern is observed with strong iron limitation of diatom growth rates and moderate iron limitation of the small phytoplankton (Figs. 11 and 12). The iron limitation of both phytoplankton classes is less severe at STNP relative to the other sites due to higher ambient iron concentrations (Figs. 10-12). The small Fig. 12. Relative nutrient cell quotas of the small phytoplankton. Cell quota relative to the maximum cell quota is displayed (see text for details). Whichever nutrient is most-limiting (lowest relative cell quota) limits phytoplankton growth rate. Nutrient pools depicted include nitrogen (square), iron (triangle), and phosphorus (diamond). phytoplankton tend to dominate the assemblage at these sites, with significant diatom biomass seen in late summer at STNP (Fig. 4). Strong grazing pressure plays a key role in keeping the small phytoplankton from blooming at these sites.
Comparing Figs. 4,5,9 and 11, it is apparent that silicon limitation leads to the demise of the diatom bloom at the NABE site. The seasonal patterns are in excellent agreement with the field data from NABE, where silicon limitation terminated the initial diatom bloom, which was then followed by a bloom of the smaller phytoplankton (Lochte et al., 1993;Sieracki et al., 1993). Diatoms become strongly Si-limited during June (Fig. 11). Then as the smaller phytoplankton bloom and diatom biomass remains low, the diatoms become N-limited during latter summer months (Figs. 4 and 11). It is interesting to note that all of the key nutrients are strongly depleted by the blooms and high export at NABE, and, thus, all cell quotas fall well below their optimum values (Figs. 5,6 and 8,9,10,11 and 12). Thus small differences in nutrient inputs/exports or recycling efficiencies could lead to the summer phytoplankton assemblage in this region being limited by Si,N,P,or Fe (Figs. 11 and 12). This is in contrast to the other regions where one nutrient is clearly limiting.
The smaller phytoplankton size class is relatively nutrient replete all year at the ARAB and ROSS sites (Fig. 12). The very different patterns in biomass of the small phytoplankton at these two sites (Fig. 4) are due to differences in grazing pressure. At the ROSS site, zooplankton biomass is quite low in the spring and is unable to prevent the spring bloom of both the diatoms and small phytoplankton (Fig. 4). This is consistent with findings from AESOPS that losses due to grazing were very low in this region, particularly early in the season . In contrast, zooplankton biomass is relatively high all year at the ARAB site, and the small phytoplankton are never able to escape grazing control (Fig. 4). The diatoms experience silicon limitation for much of the year at the ARAB site (Fig. 11).
Nutrient limitation patterns for the diazotroph group were significantly different than for the other phytoplankton (cell quotas for diazotrophs not shown). The diazotrophs were nutrient replete at all of the JGOFS sites all year except at BATS, HOT and to a lesser extent EQPAC. That is, relative cell quotas for all nutrients were B1.0 all year, and the primary controls on diazotroph biomass were light (a function of mixed-layer depth) and temperature, with a modest contribution by grazing losses. At BATS, the diazotrophs became strongly P-limited during latter summer months with mean cell quotas for phosphorus of 80%, 59%, and 79% of the maximum value during July, August, and September, respectively, a pattern similar to that observed for the diatom P quota (Fig. 11). McCarthy and Carpenter (1979) suggested P limitation of Trichodesmium spp. was a distinct possibility in the central North Atlantic. At HOT, the diazotrophs became Fe limited in late summer, with the Fe cell quotas falling to 93%, 86%, 81%, and 92% for the months July through October. The diazotrophs become slightly Fe-stressed at EQPAC during October-November.
We compare model chlorophyll concentrations with in situ data and satellite-based estimates from SeaWiFS in Fig. 13. In general, the model does an excellent job of reproducing the seasonal patterns of chlorophyll concentrations observed at the JGOFS sites (Fig. 13). Chlorophyll concentrations are on the high end of observed in situ values during the summer at BATS (Fig. 13). Persistently moderate chlorophyll concentrations are seen at the STNP, EQPAC, PFR, and KERFIX sites (Fig. 13). There is a modest spring bloom at PFR, with model chlorophyll concentrations reaching B0.7 mg/m 3 . Chlorophyll concentrations at EQPAC are somewhat low relative to the in situ measurements but are in good agreement with the SeaWiFS data (Fig. 13). Chlorophyll values during November, December, and January are lower than the in situ and SeaWiFS data at KERFIX. The SeaWiFS data are undoubtedly influenced by the nearby high chlorophyll values associated with the Kerguelen Islands and the Kerguelen Plateau (see  due to our coarse resolution grid. SeaWiFS values during December and January therefore are considerably higher than the model and the in situ measurements. Depending on circulation patterns, the in situ data also could be influenced by this nearby shallow water iron source, as the Kerfix site is located near relatively shallow waters (see Kleypas and Doney, 2001). Some of the in situ data at this site in excess of 0.8 mg/m 3 are much higher than is typically observed in open ocean areas of the Southern Ocean . Model chlorophyll values at NABE are in good agreement with the SeaWiFS data, peaking slightly later in the season than the in situ data (Fig. 13). The model does reproduce well the seasonal peaks in chlorophyll concentrations seen in the SeaWiFS and the in situ data during August (SW Monsoon) and in the SeaWiFS data during winter months at the ARAB site (Fig. 13). The wintertime increase is overestimated compared with the SeaWiFS data (Fig. 13).
Mean mixed-layer primary production values (mg C/m 3 /day) are compared with in situ measurements in Fig. 14. Also shown in Fig. 14 is the daily total grazing rate by the zooplankton (includes consumption of all three phytoplankton classes and large detritus). We first examine the production rate estimates. Several early spring in situ values at BATS were higher than 40 mg C/ m 3 /day (thus off the scale of Fig. 14). The high production rates occasionally seen in late winter/ early spring in situ data from BATS are usually associated with a rapid shallowing of the mixed layer. Our model is forced with consistently deep mixed layers over the January-March period (see Fig. 3), leading to light limitation of photosynthesis. Summertime production is consistently lower than the in situ data at the HOT and EQPAC sites (Fig. 14). At HOT and BATS, incubations are performed only during the light period and thus substantially overestimate daily production rates by not accounting for respiration losses during the dark period (see Geider, 1992). At the other sites, in situ results from 24-h incubations are shown. Biomass was somewhat lower than observed in situ at the EQPAC site (see Fig. 13), which may account for some of the discrepancy in production estimates. Model estimates of primary production are in good agreement with the field data at the NABE, PFR, and ROSS sites (Fig. 14). Boyd and Harrison (1999) reported primary production rates in surface waters at STNP ranging from low values during winter (o10 mg C/m 3 /day) to peak values of 30-40 mg C/m 3 /day during late summer. These are surface values, mixed-layer averages of primary production would be lower, although still likely higher than our model predictions, particularly during winter months (Fig. 14). Welschmeyer et al. (1993) reported mean water column production (May-September) of 707 mg C/m 2 /day over a mean depth of 73.4 m (which corresponds to 9.6 mg C/m 3 /day) during the SUPER program in the subarctic North Pacific. Fig. 14. Mean mixed-layer daily primary production rates (squares) compared with in situ data from the JGOFS sites. Also shown are daily total grazing rates (triangles) for each site.
Comparing primary production and grazing rates, it is apparent that grazing consumes a large fraction of the daily production at most of the JGOFS sites (Fig. 14). Under bloom conditions at NABE, ROSS, and to a lesser extent the ARAB site, a much smaller percentage of daily production is consumed by grazing (Fig. 14). These results are consistent with the in situ observations, which indicate that microzooplankton grazing consumes most production at low biomass locations such as the BATS site (Lessard and Murrell, 1998). In the equatorial Pacific, microzooplankton consume most of the production by the smaller phytoplankton species (Landry et al., 1995(Landry et al., , 1997Verity et al., 1996). Verity et al. (1996) concluded that microzooplankton at the EQPAC site consumed roughly half of the production by larger diatoms, with grazing by larger zooplankton and sinking/aggregation processes accounting for the residual.
In Fig. 15, we compare total mixed-layer primary production with the sinking of carbon out of the mixed layer within the large detrital pool (the particulate organic carbon, or POC flux). In the ecosystem model, there is considerable export due to mixing/detrainment processes as well as the sinking flux (see METB, this issue). Comparing Figs. 4 and 15, it can be seen that in areas with large diatom blooms (NABE, ARAB, and ROSS sites), the sinking carbon flux is high. In contrast, in areas where the small phytoplankton dominate the assemblage the sinking export flux is quite small (most notably at HOT, during summer months at BATS, and late summer at PFR, Figs. 4 and 15). There are smaller seasonal increases in the export flux at the BATS, STNP, Fig. 15. Model-predicted total mixed-layer primary production (squares) and the sinking particulate organic carbon flux (triangle) at each JGOFS location. KERFIX, and PFR sites, typically associated with increases in diatom biomass (compare Figs. 4 and 15). At EQPAC, the export flux is moderate and relatively constant, despite seasonal changes in total production (Fig. 15). At BATS, the high production values during spring are due to the deep mixed layers at this time, not to high production rates per unit volume (compare Figs. 4, 14  and 15).
The mixed-layer primary production and sinking POC flux can be compared with measurements of primary production and particulate organic carbon export from the JGOFS sites. At BATS, primary production (integrated from 0 to 150 m) peaks during the spring, with values ranging from 400 to 1100 mg C/m 2 /day (with a lot of interannual variability); summer values of typically 200-300 mg C/m 2 /day (Michaels and Knap, 1996). Our model estimates for the surface mixed layer are in very good agreement with these values, especially considering that significant production would occur below the surface mixed layer during summer months (Fig. 15). Sediment trap carbon flux at 150 m at BATS has typical summertime values of B20 mg C/m 2 /day, with peaks generally during the spring of B40-50 mg C/m 2 /day (Michaels and Knap, 1996). This springtime export peak is reproduced by the model, while summer export is quite low (Fig. 15). During stratified periods, much of the export production at mid-ocean gyre locations such as HOT and BATS may come from diatom production occurring below the surface mixed layer, deeper in the water column. Very little export production is produced in the shallow surface mixed layer. Export during the summer at BATS and all year at HOT in our model output is quite low, reflecting the low diatom production in surface waters (Fig. 15). Nelson and Brzezinski (1997) estimated that diatoms account for 15-25% of primary production at BATS. The diatoms contributed 24% of total production in our model.
At HOT, average primary production (integrated from 0 to 200 m) is 463 mg C/m 2 /day, and the particulate carbon flux at 150 m depth averages 29 mg C/m 2 /day . Much of the primary production occurs below the surface mixed layer, with significant production during summer months at 100-200 m .  estimated that the upper 25 m of the water column accounted for on average 30% of primary production. Taking these depth variations into account, our mixed-layer averages seem reasonable, although spring values may be somewhat low (Fig. 15). Recall our climatological mixed-layer depths remain fairly shallow all year at this site (Fig. 3). Thus, sporadic wintertime deep mixing and eddy upwelling events are not present in the forcings and nutrient inputs from below may be underestimated during this season (see Letelier et al., 2000). Export of sinking biogenic carbon at in the model peaks at 29 mg C/m 2 /day during October, lagging maximum nitrogen fixation by 1 month (see following sections and Figs. 15 and 17).
At STNP Boyd and Harrison (1999) estimate total water-column primary production to be 300 mg C/m 2 /day during winter months and 850 mg C/m 2 /day during spring/summer months. They estimate mixed-layer primary production at STNP to be 173 mg C/m 2 /day during February 1994 and 322 mg C/m 2 /day during May 1994 (Boyd and Harrison, 1999). Model results are in very good agreement with these mixed-layer values (Fig. 15). Welschmeyer et al. (1993) measured water-column mean spring/summer production at 707 mg C/m 2 /day. Later in the summer, when the model mixed layers are shallow (see Fig. 3), mixed-layer production is lower, but there would be substantial production beneath the mixed layer at this time (Fig. 15). Wong et al. (1999) found the POC flux to sediment traps at 200 m depth at STNP was elevated during summer months with peak values in May of B38 mg C/m 2 /day. The model predicted export at STNP peaks in April at 81 mg C/m 2 /day at the base of the mixed layer (Fig. 15). Significant portions of this flux would remineralize at depths less than 200 m.
At the EQPAC site in the fall total water-column primary production was estimated to be 1213 mg C/m 2 /day , and estimates of the sinking flux of particulate organic carbon range from 36 mg C/m 2 /day (Bacon et al., 1996) to 144 mg C/m 2 /day (Murray et al., 1996). Quay (1997) argued that the higher POC export value was needed to balance the carbon budget at this site. Our mixed-layer export is between these two estimates (mean export flux of 61 mg C/m 2 / day, Fig. 15). Total production is likely underestimated by the model due to a lower phytoplankton biomass than observed in situ and due to production occurring below the surface mixed layer (Figs. 13-15). Diatoms accounted for 29% of primary production at this site in good agreement with the estimate of 34% by Blain et al. (1997) for the equatorial Pacific.
Mean water column productivity during the latter half of May 1989 at the NABE site was 1086 mg C/m 2 /day (Martin et al., 1993), higher than our model mean for May but only slightly more than the model production rates for June and July (Fig. 15). Martin et al. (1993) estimated the export of carbon from the upper 35 m over the entire NABE study period (April 24-June 1) to be 468 mg C/m 2 /day, in reasonable agreement with the model estimates for May and June (Fig. 15). Particulate organic carbon flux at 150 m depth at NABE was B118 mg C/m 2 /day (Lochte et al., 1993).
In the Ross Sea, Smith et al. (1996) estimated total primary production of 2630 mg C/m 2 /day during January 1990, and 780 mg C/m 2 /day during February of 1992. We calculate that primary production integrated over the mixed layer during summer months averaged 1217 mg C/m 2 /day for the ROSS site and 493 mg C/m 2 /day at the PFR during AESOPS. Model results for the surface mixed layer at the ROSS and PFR sites are of comparable magnitudes (Fig. 15). A deep sediment trap (MS-3 at 1000 m) near the PFR site recorded peak organic carbon fluxes of B19.5 mg C/m 2 / day during December (Honjo et al., 2001). This lags our peak in POC export from surface waters by 1 month (Fig. 15).
Seasonal variations in the size of the two detrital pools are shown in Fig. 16. The large (sinking) detrital pool is typically smaller than the non-sinking pool because it sinks out of the mixed layer rapidly once formed (see Fig. 15). The sinking pool reaches maximum values in the presence of large diatom blooms at the ARAB and ROSS sites (compare Figs. 4 and 16). At NABE, particulate sinking export is high during the spring bloom (Fig. 15), but large detritus is rapidly exported from the shallow mixed layers, and, thus, mixed-layer concentrations remain fairly low (Figs. 3 and 16). This is in contrast to BATS, where deep mixed layers allow significant amounts of large detritus to accumulate even though export is moderate (Figs. 15 and 16). The size of the sinking detrital pool is closely tied to diatom production in the model. Thus, there is a small buildup in the spring at BATS, KERFIX, and the PFR sites when some diatom production is occurring. The large detrital pool at each of these sites is very small later in the summer season when diatom biomass is very low (compare Figs. 4 and 16). A similar pattern is observed all year at HOTS, with some increased diatom production and export in late summer associated with increased nitrogen fixation (compare Figs. 4,16 and 17). The small detrital pool increases during spring/summer months at most sites as production exceeds remineralization (BATS, HOT, STNP, NABE, KERFIX, and PFR, Fig. 16). This pool then decreases during fall/winter as remineralization exceeds production. Both detrital pools are nearly completely remineralized during winter at the high-latitude sites (STNP, ROSS, NABE, KERFIX, and PFR, Fig. 16).
The small detrital pool is meant to include the labile and semi-labile portions of the dissolved and suspended organic matter. Thus, a rough comparison between our model small detritus and measurements of dissolved organic carbon (DOC) is possible. Carlson et al. (1994) reported that DOC builds up during the spring season at BATS and then is partially consumed during the summer/fall period. The seasonal buildup ranged from B1-7 mM carbon over winter values in surface waters, with a larger seasonal increase beneath the surface mixed layer (Carlson et al., 1994). The seasonal pattern of small detrital matter at BATS is in agreement with these measurements (Fig. 16). At the EQPAC site, Carlson and Ducklow (1995) reported surface water DOC values approximately 20 mM higher than deep-water samples, with little seasonal variability. Our small detrital carbon values for this region are considerably lower. In the Southern Ocean along 61W, K. ahler et al. (1997) reported DOC values ranging from B5-15 mM above deep-water values in the upper 100 m of the water column. Similar values have been reported for the Ross Sea (Carlson et al., 1998). The model predicted small detrital pool for the KERFIX and PFR sites is of comparable magnitude. During AESOPS, DOC values in surface waters as high as 25 mM above winter values were observed (compare with Fig. 16).
We next examine the role of nitrogen fixation in the model. Monthly mean nitrogen fixation rates at the JGOFS sites are shown in Fig. 17. Highest rates are observed at HOT where annual N fixation is 34.9 mmol N/m 2 /yr. This total increased to 71.8 mmol N/m 2 /yr in an alternate Fig. 16. Biogenic material in the sinking (square) and non-sinking (triangle) detrital pools at monthly intervals in carbon units. iron-replete model run (see METB). Thus, Fe limitation has a significant impact on N fixation rates at HOT in the model. Estimates of total water column N fixation at HOT range from 31-51 mmol N/m 2 /yr (Karl et al., 1997), with significant fixation below the surface mixed layer . The seasonality of N fixation and diazotroph biomass is also consistent with field observations, with peak values in mid-to late summer and low values during winter and early spring months (Karl et al., 1997). The release of DON by diazotrophs in the model at HOT clearly has a stimulatory effect on the other phytoplankton groups. The N cell quotas of diatoms and the small phytoplankton as well as the biomass of both classes increase (Figs. 4,11 and 12). Sinking export and primary production also increase in conjunction with increased N fixation (Figs. 14,15 and 17). This stimulatory role of diazotrophs at this location is well documented from in situ data Karl, 1996, 1998;Karl et al., 1997). We quantify this stimulatory effect of the diazotrophs at the global scale in METB.  estimated that, on average, Trichodesmium spp. accounted directly for 4% of primary production and at least 27% of new production at the HOT site, with considerable interannual variability. Nitrogen fixation is the dominant source of new nitrogen during summer months in the model (Fig. 7). The diazotrophs account for 2.5% of total primary production in the model output on an annual basis.
The other sites in the model with significant N fixation are BATS, ARAB, and EQPAC, with negligible N fixation and diazotroph biomass at the other locations (Figs. 4 and 17). The seasonal cycle at BATS is similar to HOT, with a peak in mid-to late summer after a period of prolonged shallow mixed layers (Figs. 3 and 17). Maximum N fixation rates at BATS are significantly lower than HOT and annual N fixation was also lower at 13.3 mmol N/m 2 /yr. Orcutt et al. (2001) estimated mean N fixation at BATS by Trichodesmium spp. to be 15.0 mmol N/m 2 /yr, with considerable interannual variations, peaking during summer months as predicted by the model (Fig. 17). Unlike at HOT, Trichodesmium biomass and N fixation maxima were typically in surface waters (Orcutt et al., 2001).
Annual N fixation at the ARAB site was 15.3 mmol N/m 2 /yr, with rates peaking during the stratified period of April-June (compare Figs. 3 and 17). At mixed layers deepened during summer months, N fixation was decreased due to light limitation (Figs. 3 and 17). Capone et al. (1998) documented an extensive bloom of Trichodesmium erythraeum in the central Arabian Sea and noted that the bloom occurred during periods of low wind speed and high atmospheric iron flux. The timing of this bloom coincides with our seasonal increase in diazotroph biomass during the spring intermonsoon. Capone et al. (1998) estimated an areal N fixation rate of 0.129 mmol N/ m 2 /day during the bloom, in excellent agreement with the model predictions for May and June (Fig. 17). They also estimated a background (non-bloom) N fixation rate for other seasons of 0.04 mmol N/m 2 /day (Capone et al., 1998), only slightly higher than the model estimates during late summer/fall months.
At EQPAC there was a low relatively constant amount N fixation all year with a total of 8.2 mmol N/m 2 /yr. Although a relatively low number, this N fixation accounted for a large portion of new production, as nitrate uptake by the diatoms and small phytoplankton was relatively low due to a strong ammonium preference. The in situ data showed a large variability in mixed-layer depths not captured in our monthly climatological forcings (Fig. 3). It is likely that the variability in light levels caused by this short-term variability in mixed-layer depths would suppress N fixation in situ, and thus our model likely overestimates diazotroph biomass and N fixation at this location.
Mean monthly calcification rates and the sinking export of CaCO 3 within the large detrital pool predicted by the model are displayed in Fig. 18. Where grazing is the dominant loss term for the small phytoplankton group, approximately 50% of the CaCO 3 produced by phytoplankton is remineralized in the surface layer (Fig. 18). There is increasing field and laboratory evidence for shallow depth remineralization of CaCO 3 , likely due to grazing processes (Harris, 1994;Milliman et al., 1999;Balch et al., 2000). Under bloom conditions when aggregation/sinking losses become important, more than half of the phytoplankton CaCO 3 produced is exported (i.e. the summer small phytoplankton bloom at NABE,Figs. 4 and 18). Calcification rates are negligible at the ROSS site and quite low at the KERFIX and PFR sites due to the effects of low sea-surface temperatures in our parameterization (Fig. 18, Appendix A). Similarly, calcification is depressed at HOT and at BATS during summer months due to the effect of strong nutrient limitation in our parameterization. There are few in situ measurements of calcification to compare with these model results. Balch and Kilpatrick (1996) estimated water column calcification rates along 1401W that ranged between 2 and 10 mmol C/m 2 /day near the equator. Balch et al. (2000) measured calcification rates in the Arabian Sea that were again highly variable but on average about 5% of integrated primary production (or a water column calcification of 6.2 mmol C/m 2 / day based on mean production of 1485 mg C/m 2 /day, Table 2 of Balch et al. (2000)). This would be for the whole water column, and thus our model mixed-layer results are lower but in reasonable agreement (Fig. 18). Note that calcification is less affected by light levels than photosynthesis and often does not decrease as much as photosynthesis at low light levels, sometimes leading to subsurface maximum in calcification/photosynthesis ratios (Balch et al., 2000; see also Balch et al., 1992). Highest calcification rates were observed during the SW Monsoon in the coastal upwelling regions (Balch et al., 2000). In the model, calcification averaged B3% of primary production during spring months, with lower values later in the season. Balch et al. (2000) also noted that mean calcification rates were 3 times higher during the July 17-August 14 cruise (SW Monsoon) compared with the cruise from October 24-November 25 during the NE Monsoon. Encouragingly, a similar seasonal pattern is seen in our model output (Fig. 18). Perhaps the best way to evaluate our parameterization of calcification is not a site-to-site comparison but an examination of the resulting basin to global scale patterns as we do in the companion paper (METB).

Discussion
Overall the marine ecosystem model described here does a good job in reproducing patterns observed in the field data from nine diverse oceanic locations. Early spring nutrient concentrations are in most cases quite similar to the in situ data (Figs. 5,6,8,9 and 10). This indicates that our bottom boundary conditions are introducing roughly correct concentrations to the surface mixed layer. The seasonal drawdown patterns of ambient nutrients is also in good agreement with the field data, suggesting biological uptake, export, and regeneration are being correctly estimated by the model (Figs. 5,6,8,9 and 10). Seasonal patterns of primary production, f-ratio, N fixation, phytoplankton biomass, chlorophyll concentration, and nutrient limitation seen in the model output are also generally in good agreement with the in situ data (several exceptions are discussed below). The timing and magnitude of the sinking flux of biogenic carbon also largely agree with field measurements (Fig. 15). Comparisons with sediment trap data are complicated by the fact that significant portions of the POC flux may be produced below the surface mixed layer (i.e. Goldman, 1993). Perhaps most importantly, phytoplankton blooms occur or do not occur at each station in agreement with the field data, and for the correct reasons. Thus, N limitation is the key mechanism in the central ocean gyres (BATS and HOT), while Fe limitation of the large diatoms determines bloom dynamics in the HNLC regions (STNP, EQPAC, KERFIX, and PFR). This ecosystem model thus can be applied with some confidence to the global domain (METB).
The spring diatom bloom at the PFR site observed during AESOPS was considerably stronger than predicted by our model. Silicate and nitrate drawdown were significantly stronger during December in the in situ data relative to the model (Figs. 5 and 9). Diatom growth was strongly Felimited in the model during this period (Fig. 11). The discrepancy may be related to physical forcing and nutrient inputs in the model. A meandering jet such as the Antarctic Polar Front increases the flux of nutrients from subsurface waters through localized areas of upwelling (Moore et al., 1999. These mesoscale physical dynamics are not resolved in our global grid; thus, the flux of iron from subsurface waters is probably underestimated. Similarly, melting sea ice can be a source of iron to surface waters (Sedwick and DiTullio,1997;Measures and Vink, 2001). This iron source also is not included in our model. The seasonal sea-ice cover nearly reaches the PF along 1701W and in many years it likely does reach the PF just upstream of the JGOFS study site (Moore et al., 1999. Thus, iron inputs during spring at the PFR site likely have been underestimated in the model, preventing the observed diatom bloom from developing. Primary production and phytoplankton biomass indices tended to be lower than observed at the EQPAC site. At the EQPAC site, our smoothly seasonally varying mixed-layer depths do not capture the short-term variations seen in the in situ data (Fig. 3). These short-term variations in mixed-layer depths would introduce additional nutrients to the mixed layer. Our coarse resolution physical grid likely underpredicts the flux of nutrients into the surface mixed layer due to upwelling. We smoothed the vertical velocities at the base of the mixed layer (see Methods section). Lastly, our mixed-layer formulation does include lateral advection, which is an important process in this region.
Inclusion of our simple Fe-stress effect on ammonium preference improved the fit to the JGOFS data. In a large number of model runs, the model fits to observed ammonium, nitrate, and f-ratio data were improved by including this increased preference for ammonium over nitrate under Felimited conditions. Our inclusion of an iron effect on cell silica quotas, and thus on the ratio of silicate uptake/nitrogen uptake, was less conclusive in regard to the in situ data (Fig. 9). We return to this topic of silicon/iron co-limitation and the iron influence on Si/N uptake ratios in the companion paper (METB).
Our inclusion of the atmospheric source for dissolved silicate had a generally negligible effect on ecosystem parameters, including ambient silicate concentrations, except at the NABE and ARAB sites, regions where Si limited diatom growth (Fig. 11). The wintertime diatom bloom at ARAB (see Fig. 4) was restricted to January in an alternate model run with no atmospheric Si deposition. The effect at NABE was smaller; the diatom spring bloom at peaked somewhat later when the atmospheric Si source was not included. Mean diatom biomass was 25% lower during May and 7% higher during June without the atmospheric Si inputs. Monthly mean dissolved-silicate concentrations in runs with and without the atmospheric Si input were nearly identical during these 2 months at NABE. Thus, the atmospheric Si input can significantly affect ecosystem dynamics in some regions, particularly areas with high dust flux inputs like the Arabian Sea. However, at the global scale atmospheric Si input has little impact on primary production or spatial patterns of nutrient limitation (see METB).
In contrast the atmospheric iron source had a large impact on ecosystem dynamics at nearly all the sites. Without the atmospheric iron source seasonal patterns of nutrient limitation and phytoplankton biomass were drastically different at BATS, HOT, STNP, and the ARAB site. The EQPAC, KERFIX, ROSS, NABE, and PFR sites were less affected (although phytoplankton biomass was lower) as most iron comes from subsurface sources at these locations. We examine the relative importance of the atmospheric iron source in more detail in METB.