The role of phytoplankton dynamics in the seasonal and interannual variability of carbon in the subpolar North Atlantic – A modeling study

. We developed an ecosystem/biogeochemical model system, which includes multiple

ranging from days to years. Here, we developed a coupled biogeochemical-physical model system, which includes multiple phytoplankton functional groups and carbon cycle dynamics, and applied it to assess the role of phytoplankton dynamics in the seasonal and interannual variability of carbon in the subpolar North Atlantic. The rationale for the site selection relates to a combination of environmental and biological factors, e.g., large range in mixed layer depth (MLD), low mean horizontal advection, clear seasonal succession of phytoplankton species, and bloom intensity. Spring phytoplankton blooms in the surface water of the North Atlantic Ocean and adjacent seas are known to cause a precipitous reduction of surface water partial pressure of carbon dioxide (pCO 2 ) and the concentrations of CO 2 and nutrients on a seasonal basis (Takahashi et al., 1993). In addition, research during the last decade has increased awareness of the relationship between key phytoplankton groups and their pivotal roles in the biogeochemical cycles of a range of elements (Boyd et al., 2010). For example, coccolithophores occupy a central role in the carbon cycle by the conversion of dissolved inorganic carbon (DIC) to both particulate organic carbon (POC) and particulate inorganic carbon (PIC) forms, albeit uncertainties exist over whether coccolithophore blooms are net sinks or sources of CO 2 to the atmosphere (Boyd and Trull, 2007). The ratio of PIC to POC in exported biogenic matter determines the relative strength of the biological carbon pump and consequently the flux of CO 2 across the surface ocean-atmosphere interface. Other functional groups, such as diatoms and dinoflagellates, are also important in the sequestration of carbon from the atmosphere to the deep ocean. The ability of some species of diatoms to form chains with built-in silica ballast and to produce large fast-sinking aggregates during the declining phase of the blooms means that they are key vectors in exporting and sequestering POC to the deep ocean (Lampitt, 1985). The competitive success of phytoplankton functional groups has important links to the environmental conditions. Although the predictive outcome of competition between two or more species has seldom been tested in the marine environment, some laboratory and mesocosm experiments reveal interesting results. Elaborate laboratory culture competition experiments with multiple marine phytoplankton species (coccolithophores excluded) (Sommer, 1994) and mesocosm experiments (Egge and Aksnes, 1992) show that diatoms dominate at high Si:N ratios, while diverse flagellates tend to dominate at low Si:N ratios. As Si is not regenerated above the MLD during a season, the growth period for diatoms fades out as the pycnoline is enforced and may be revoked every time a vertical mixing event occurs, as long as the available light is sufficient. Sommer (1994) showed no apparent effect of light intensities on the succession of functional groups, although it may exert some influence at the species level on the transition along the Si:N gradient. Nevertheless, light intensity has been shown to be very important on the blooms of certain species of coccolithophores. The large E. hux-leyi bloom south of Iceland in 1991 was modeled and high light was shown to be a positive cause (Tyrrell and Taylor, 1996). Under the intense bloom conditions ,mixing was not very deep and therefore average light intensities in the surface layer were high.
The physical characteristics of the North Atlantic subpolar waters also have an impact on the biogeochemical and carbon cycle variability. The waters surrounding Iceland are characterized by the cold polar water of the East Greenland Current and Arctic water of the East Icelandic Current from the north, and the warm North Atlantic water of the Irminger Current from the south (Gudmundsson, 1998). Figure 1a shows a map of the subpolar North Atlantic and Nordic Seas and the mean location of the Arctic Front (AF). Sea ice is present throughout the year along the east coast of Greenland, with the ice edge extending farthest offshore in winter-spring, and retreating in the fall. Ice melt in summer-fall freshens (32-34 psu) the surface waters north of the AF. South of the AF waters are saltier (∼35 psu) and warmer (8 to 12 • C), as a result of Atlantic water intrusions. These hydrographic characteristics have an impact on the seasonal vertical mixing, as shown in model simulations. During winter and spring, MLDs average 100 to 500 m south of the AF, with deeper values in the Irminger and Icelandic Basins; in summer, the vertical stratification is significant, with MLDs less than 20 m (de Boyer Montégut et al., 2004;de Boisséson et al., 2010;Carton et al., 2008). North of the AF, the MLDs are shallow (< 40 m) throughout the year due to the southward advection of fresher polar and Arctic waters, and ice melting during summer-fall. Some accounts of winter deep mixing are available for the region based on regular XBT surveys and CTD casts (G. Reverdin, personal communication, 2012). MLDs as deep as 500 m were found in the Irminger Current, in particular west of the Reykjanes Ridge, and occasionally close to Greenland and southwest of Iceland. In the central Irminger Sea, MLDs are not as deep, mostly ranging from 100 to 150 m. However, there are indications from Argos floats that deep convective chimneys occasionally occur in the central Irminger Sea (as in early 2008) down to nearly 1000 m (Vage et al., 2009). These large seasonal changes in stratification and vertical mixing play an important role in the euphotic zone nutrient renewal and on the onset and duration of the phytoplankton spring bloom (Henson et al., 2006). Variability in the intensity of primary production in general and of the timing of spring bloom in particular (Henson et al., 2009) affects the population dynamics of higher trophic levels, such as the commercially important Atlanto-Scandian herring (Jakobsson, 1978) in the region. The annual phytoplankton carbon production in the waters around Iceland is quite high. According to a recent estimate (Zhai et al., 2012), based on available regional measurements on the photosynthetic capacity and satellite biomass data, the annual primary production is 179 ± 36 and 238 ± 22 gC m −2 y −1 for the Arctic and the Atlantic waters, respectively. In addition, seasonal and interannual changes in phytoplankton production are tightly coupled to atmospheric CO 2 uptake and surface ocean pCO 2 variability (Corbiere et al., 2007), and therefore a major component of the carbon cycle.

Ice-ocean model (GPOMZ)
Physical forcing for the 1-D ecosystem/carbon model, ECO1D-2.0(location shown by the black triangle in Fig. 1a at 30 • W and 60 • N), is provided by a three-dimensional (3-D) coupled ice-ocean model of the Arctic and North Atlantic . The sigma-coordinate version has been previously used to simulate North Atlantic and Arctic variability (Hakkinen, 1999(Hakkinen, , 2001Hakkinen and Geiger 2000). Monthly temperature and salinity profiles obtained from GPOMZ are interpolated to the time steps and depths required by ECO1D-2.0. These profiles are assimilated at each time step by the 1-D model. A detailed description of GPOMZ is provided in Appendix A. In this subsection we present a few GPOMZ products to describe some of the physical properties of the region. In Appendix A we show three figures with physical properties derived from GPOMZ. Figure A1 shows the annual SST composites for 2009, derived from GPOMZ and MODIS-Aqua; Fig. A2 shows a comparison between GPOMZ and MODIS SST at ECO1D-2.0 simulation site, and Fig. A3 shows a comparison between the seasonal MLD from de Boyer Montégut et al. (2004) and GPOMZ at the ECO1D-2.0 site. Although there is a reasonable agreement between model and MODIS SST fields in Fig. A1, the Arctic Front is further away from the study site compared to MODIS data. However, as it will be shown later in section 4 and Appendix A, this has little or no effect on our results, since the bias-corrected SST at the modeling site agrees well with the data. The bias correction was done based on WOA05 monthly climatology and the resulting SST time series at the ECO1D-2.0 site was validated based on in situ observations and MODIS data (see caption in Fig. A2 for explanation). Figure 1b shows a smaller-scale map containing the locations of the in situ data available for model evaluation.
The 1-D model location lies on the western edge of the Reykjanes Ridge and within the Irminger Current, in an area where winter mixing is extremely vigorous (mixed layer depths approaching 600 m) (Bailey et al., 2005) and springsummer phytoplankton blooms are substantial (Fernandez et al., 1993;Holligan et al., 1993b;Weeks et al., 1993). The convergence of the polar and North Atlantic water masses form the Arctic Front, which varies slightly in location seasonally (Fig. 2) and is identified by strong gradients of sea surface temperature (SST) and salinity (SSS). The AF runs east-west, approximately along 66 • N to the west of Iceland, and north-south between 12 • W-8 • W, and then further in northerly direction from approximately 69 • N, northeast of Iceland. The above features are clearly shown in Fig. 2, which shows the seasonal SST and SSS fields for 2005 obtained from GPOMZ . GPOMZ products are available at monthly intervals for 1948-2009. We chose to show fields for 2005 as an example because it coincides with the year during which most of the surface ocean pCO 2 data are available (Chierici et al., 2009;Olsen et al., 2008).
The climatologic  seasonal surface currents from GPOMZ are shown in Fig. 3. Although there may be some differences between the model currents and what is known of the regional circulation patterns (Poulain et al., 1996;Gudmundsson, 1998), Fig. 3 shows a general cyclonic circulation in the Irminger Sea north of 59 • N, and a convergence south of the Denmark Strait between the cyclonic circulation and the East Greenland Current, in agreement with observations.

Ecosystem-carbon model (ECO1D-2.0)
The mixed layer component of ECO1D-2.0 originates from an existing one-dimensional physical-biogeochemical model (Signorini et al., 2001a;Signorini et al., 2001b), which utilizes a turbulence closure mixed layer scheme (Mellor and Yamada, 1982), identical to the scheme used for GPOMZ. ECO1D-2.0 has a vertical coordinate system that provides parameter values, including horizontal velocity components, temperature, salinity, and the vertical diffusivity coefficient, k v , at each time step and grid point. The surface boundary layer is resolved more accurately by using a stretched vertical coordinate with higher resolution near the surface.
The biogeochemical component of ECO1D-2.0 derives from a previous ecosystem model, ECO1D-1.0, configured for and applied at the Bermuda Atlantic Time Series (BATS) site, described in Signorini et al. (2003). ECO1D-2.0 includes additional conservation equations for diatoms, coccolithophores, calcium carbonate (CaCO 3 ), silicate (SiO 2 ) and alkalinity. Figure 4 shows a diagram illustrating the model components and their couplings. The details of the ecosystem model are described in Appendix B. Figure 5 (to be discussed later in Sect. 4) shows the satellite-derived net primary production (NPP) and calcification rate (P PCaCO3 ) for June 1998, during which the strongest coccolithophore bloom occurred during the SeaWiFS mission (Raitsos et al., 2006), with the location of the modeling site (black triangle at 30 • W and 60 • N). NPP was obtained from monthly CbPM files and the P PCaCO3 monthly composite was computed using the algorithm of Balch et al. (2007). The site lies at the western edge of the Reykjanes  Ridge in the Irminger Sea, a region of high phytoplankton productivity. The local depth is ∼1500 m, but the model vertical grid extends to 1000 m only, which accommodates the deepest MLD. ECO1D-2.0 can be forced by winds and other relevant atmospheric parameters to calculate heat and freshwater fluxes originating from NCEP-II Reanalysis products, which determined the 28-yr period of simulation . As an option, surface temperature and salinity can be specified at the surface instead of heat and freshwater fluxes. The physical component of ECO1D-2.0 is configured on a vertically-stretched logarithmic grid, while the biogeochemical component uses a uniform vertical grid with 1 m resolution. Deep water (z > 200 m) nutrient values are nudged within the bounds of the observed monthly climatology from the World Ocean Atlas 2005 (see details in Appendix B), and temperature and salinity values at all depths are nudged to values derived from the GPOMZ using the Newtonian relaxation method (see Appendix B for details). The Newtonian relaxation method (nudging) is a simple form of data assimilation. Surface salinity (SSS) from GPOMZ is imposed at the top layer of the 1-D model, which accounts for all processes that alter SSS, including ice melting-freezing and precipitation-evaporation. The Reynolds and Smith optimally interpolated (RSOI) SST is imposed at the model surface instead of heat flux.
Deep water dissolved inorganic carbon (DIC) is nudged (relaxation time of 2 days) to values obtained from a dataderived equation (DIC vs. T) and model temperature at each time step (0.5 h), plus has a superimposed DIC decadal trend  consistent with observations . The total alkalinity (TA) is nudged in deep water to values also obtained from a dataderived equation (TA vs. T and S), with a relaxation time scale of 10 days. The DIC and TA equations are based on CARINA observations, with details described in the Appendix B. The adopted method of forcing the one-dimensional biogeochemical model (ECO1D-2.0) with physical fields from GPOMZ relies on the assumption that vertical mixing pro-cesses, e.g., convective overturning and wind mixing, play a much more important role in the biological variability than the lateral advective processes. To better quantify and justify this assumption, an assessment of the nitrate transport via vertical mixing and horizontal advection was performed using a combination of monthly climatologic nitrate profiles from WOA05 and monthly climatologic  physical parameters (eddy vertical diffusivity, k v , and north and east velocity components, u and v) from GPOMZ. The calculation was done within ±2 • of the modeling site at 60 • N and 30 • W, and within the depth interval from 100 to 500 m, in the upper part of the water column. The minimum depth of 100 m was chosen to avoid the euphotic depth, within which nitrate is consumed by phytoplankton, while the maximum depth of 500 m is dictated by the availability of monthly WOA05 nitrate. Spatial and temporal averages were done to the resulting horizontal divergence and vertical mixing components of the nitrate transport, u∂NO 3 /∂x + v∂NO 3 /dy and ∂/∂z (k v ∂NO 3 /∂z), respectively. The results demonstrate that our assumption is a valid one. The annual mean local horizontal divergence is 0.005 mmol m −3 d −1 , while the vertical eddy diffusivity transport is 2.628 mmol m −3 d −1 , a factor of ∼550 times larger than the local horizontal divergence. The values for summer-fall (January-March and October-December) and spring-summer (April-September) are 0.007 and 0.003 mmol m −3 d −1 for horizontal divergence, and 5.124 and 0.133 mmol m −3 d −1 for vertical diffusion, respectively. Although the vertical diffusion of nitrate is significantly reduced in spring-summer, it is still ∼52 times larger than the horizontal divergence.

Satellite and in situ data sets
We rely on a combination of satellite data and field observations to provide an evaluation framework for the model. The satellite-derived data sets consist of SeaWiFS Chl a, Reynolds and Smith optimally interpolated (RSOI) SST (Reynolds and Smith, 1995) and sea ice concentration. These data were obtained from NOAA_OI_SST_V2 data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their web site at http://www.esrl.noaa.gov/psd/. We use primary production from the carbon-based productivity model (CbPM, (Behrenfeld et al., 2005)). The satellitederived PP data were obtained from the ocean productivity web site at Oregon State University (http://www.science. oregonstate.edu/ocean.productivity/). The inputs for CbPM are Chl a and backscatter at 443 nm ,estimated using the Garver-Siegel-Maritorena (GSM) semi-analytical algorithm (Maritorena et al., 2002) and monthly SeaWiFS water leaving radiances, SeaWiFS cloud-corrected PAR data, SeaW-iFS mixed layer light attenuation coefficients at 490 nm, advanced very high resolution (AVHRR) SST, and monthly mean regional mixed layer depths from the Fleet Numeric Meteorology and Oceanography Center (FNMOC) (Monterey, California). In addition, a limited amount of in situ PP data was obtained from C-14 incubations conducted during early August 2002 on a Marine Productivity (MarProd) cruise sponsored by the Natural Environment Research Council (UK) onboard the RRS Discovery. Taxonomic data (cell counts) were obtained from the Continuous Plankton Recorder (CPR) database (standard area B6 south of Iceland, http://www.sahfos.ac.uk/data-archive/standard-areas.aspx).
Surface ocean pCO 2 data for 2005 were acquired from onboard the container ship M/V Nuka Arctica (Chierici et al., 2009;Olsen et al., 2008). In addition, we used DIC, alkalinity, surface ocean pCO 2 , and nutrient in situ data from SURATLANTE (Corbiere et al., 2007;Metzl et al., 2010) The locations of all in situ data used to validate the 1-D ecosystem model are shown in Figure 1b. The atmospheric pCO 2 required to obtain the CO 2 flux at the atmosphereocean interface was obtained from GLOBALVIEW-CO 2 (NOAA, ftp.cmdl.noaa.gov). Note that these data sets used for model validation are from the surface or near-surface, which are not part of the deep water data utilized for model forcing and boundary conditions.

Ecosystem-carbon model evaluation
Using the data sets shown in Fig. 1b and discussed in Sect. 3, model versus data comparisons were compiled and analyzed. Due to the nature of the spatial distribution of the in situ data surrounding the model location, it is expected that some of the mismatches between model and observations are inherently related to the spatial variability and patchiness of measured quantities. For instance, the patchiness of primary production and calcite production can be readily seen in the satellite maps shown in Fig. 5. However, in spite of this spatial mismatch between the model single point simulation and available measurements in multiple locations around the model site, with few exceptions, the overall agreement between model and observations is quite good. It will be shown later in this section that the majority of model versus observed parameters are statistically correlated with p-values well below 0.05.
A quantitative approach based on metrics, assessment indices and skill scores is provided to evaluate model performance. For this purpose, we use a Taylor diagram (Jolliff et al., 2009) to conduct the model evaluation based on all available in situ and satellite observations. The normalized standard deviation (σ *) and the correlation coefficient (R) from the model (m) to reference field (r) comparisons may be displayed on a single Taylor diagram (Fig. 6). Note that σ * is displayed in both x and y coordinates of the polar diagram of Fig. 6. The Taylor diagram is a polar coordinate diagram that assigns the angular position to the inverse cosine of R (cos −1 (R)). A correlation of zero is thus 90 • away from a correlation of 1. The radial (along-axis) distance from the origin is assigned to the normalized standard deviation. The reference field point (black circle in Fig. 6), which is comprised of the statistics generated from a redundant reference to reference comparison, is indicated for the polar coordinate (1.0, 0.0). The model to reference comparison points may then be gauged by how close they fall to the reference point. This distance is proportional to the normalized unbiased root-mean-square difference, as defined by the equation Inspection of Fig. 6 reveals that the great majority of model versus reference comparisons fall within significant R values (0.7-0.9) and σ * close to 1. Except for alkalinity (TA), all other R-values are significant at the 95 % confidence level, with p-values well below 0.05. The TA reference data have a relatively low correlation with the model, probably as a result of the low local TA variability, even though the bias between model and observations is small. The Taylor diagram does not provide an estimate for the bias, but Table 2 provides surface ocean means and standard deviations (STDs) for all the variables shown in the Taylor diagram. The means and STDs for the model variables are very close to the observed values.  Despite having significant R values, the model has somewhat larger nitrate σ * than the observed reference. Time series of observed and model nutrients ( Fig. 9 in Sect. 6) reveal that the model underestimates all nutrients during 2005-2006, especially nitrate. However, the low winter nutrients are consistent with warmer winter SSTs, as given by both model and observations (Fig. 10 in Sect. 6), and consequently shallower MLDs. Since these are the only two years during which the model underestimates nutrients, it is possible that the monthly mean physical fields and monthly climatological nutrients used to derive the forcing for the model were of insufficient resolution for capturing the smaller scale variability responsible for nutrient enhancement that occurred during those two years. Surface ocean pCO 2 references are in good agreement with the model (R ∼ 0.7 and σ *∼1.0). Finally, the model PP provides good agreement with the PP estimates derived with satellite-based CbPM (Fig. 8 in Sect. 6) using SeaWiFS data at the precise location of the model simulation (r 2 = 0.68).

Sensitivity runs
In order to evaluate the effects of model forcing and biological uptake of carbon on the model results, a couple of sensitivity runs (Experiments B and C) were performed and compared with the baseline run (reference, Experiment A), which was conducted following the description given in 2.2. One of the sensitivity runs (Experiment B) was conducted without the relaxation of DIC and TA to the empirical values derived from the data. In addition, the deep water trend in DIC was turned off and deep water values of DIC and TA were maintained constant and equal to their initial values. The second sensitivity run (Experiment C) consisted of maintaining the DIC and TA forcing strategy used in the reference run but with the biological effects turned off, i.e., abiotic ocean case. Four surface ocean carbon-related variables were analyzed: DIC, TA, pCO 2 , and air-sea CO 2 flux. The climatologic  seasonal variability of these variables is shown in Fig. 7a. There are notable differences between the three different model experiments. Although Experiment B (without relaxation and DIC bottom trend) differs slightly from the reference, the largest differences are apparent in the abiotic run (Experiment C), since the biological sources and sinks that affect the DIC and TA seasonality were turned off in the model. These are the net community production, which is a sink for DIC (CO 2 uptake) and source for TA (nitrate uptake consumes hydrogen ions), dissolution (source for both DIC and TA) and CaCO 3 production (sink for both DIC and TA). With these terms turned off, only thermodynamic and vertical mixing effects are at work in the model and that is reflected in the seasonal changes of DIC, TA, pCO 2 , and CO 2 flux. Note that without the biology (Experiment C), the ocean becomes a source of CO 2 for most of the year, while the reference run (Experiment A), which agrees best with observations (see next section), shows an ocean that is a sink of CO 2 for most of the year and a small source during April-May. This result agrees with previous studies that show that biology is critical to these seasonal changes at this location and that vertical mixing also plays a major role (Takahashi et al., 1993(Takahashi et al., , 2002(Takahashi et al., , 2009Ullman et al., 2009;Bennington et al., 2009;Metzl et al., 2010 Table 1, which shows significant increase in the biases for the no relaxation (Experiment B) and abiotic (Experiment C) runs when compared to the reference run (Experiment A). For example, The biases for surface ocean pCO 2 are +3 µatm, −15 µatm, and +44 µatm for Experiments A, B, and C, respectively. The difference in interannual variability of surface DIC and ocean pCO 2 from two of the model experiments, one  including deep water DIC and TA relaxation and bottom DIC trend (Experiment A) and one without it (Experiment B), is shown in Fig. 7b. The available data are also shown for comparison. There is a striking difference between these two experiments. Experiment A shows a very good match with the available data, while the results from Experiment B have a very poor agreement with the data, including a significant reduction in the seasonal amplitude of DIC and pCO 2 . This result illustrates the degree to which the model results are di- rectly a function of the setting of DIC, TA and nutrients at >200 m, how much the deep reservoir impacts surface behavior, and more clearly that it is very important to be included in the model to appropriately represent the interannual variability driven by non-local processes. This ratio provides a measurement of how much of the daily light exposure the phytoplankton experiences during the growth season. Apparently, every year the most productive periods (when Zeu is shallowest) occur during highly stratified conditions (with shallow MLDs). That is an indication that light availability for photosynthesis is a major limiting factor, causing seasonal variability of phytoplankton growth and probably also affecting the species succession in the region (Nanninga and Tyrrell, 1996;Tyrrell and Taylor, 1996). Overall, the model agrees well with the corresponding satellite-derived products. In fact, the model reproduces the satellite PIC interannual variability rather well, albeit with the model showing a tendency to produce peak PIC concentration higher and later than the satellite observations. Notably, the interannual changes in peak PIC concentrations are correlated with light intensity, as they follow the variability of the Zeu/MLD ratio. The variability in this ratio is primarily driven by changes in MLD, since the standard deviation in MLD during the summer months (June-July-August, 1981-2008) is about 10 times larger than the standard deviation in Zeu (24.8 ± 10.8 m and 72.6 ± 1.4 m, respectively). The growth of calcite-forming coccolithophores in the model is, by design, three times more demanding on light exposure than the other two functional groups, since it has the lowest assigned P vs. I initial slope of 0.033 (Wm −2 ) −1 day −1 (see Table B2). The coherent model prediction of growth of coccolithophores and the variations in PIC concentration are consistent with the interannual changes in PIC observed by the satellite, which confirms the relevance of the Zeu/MLD ratio. Time series of model and in situ measurements of nutrients (NO 3 , PO 4 , and SiO 2 ) for the time period, for which there are data (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), are shown in Fig. 9; SST and surface ocean DIC and pCO 2 for the same time period are shown in Figure 10, also with corresponding in situ measurements. The GLOBALVIEW atmospheric pCO 2 , which was used to calculate pCO 2 and sea-air CO 2 flux, is shown in the bottom panel of Fig. 10 along with the surface ocean pCO 2 . The seasonal and interannual changes in nutrients and carbon are substantial and driven primarily by changes in biological production, vertical mixing, SST, and sea-air CO 2 flux. The observed and model seasonal partitioning of phytoplankton functional groups (diatoms, dinoflagellates, and coccolithophores) are shown in Fig. 11. The observed seasonality originates from the aggregated, monthly-averaged indices of phytoplankton groups from the CPR data (cell counts) in the standard area B6, south of Iceland (Fig. 11a). The size of standard area B6 is relatively large (59 • N to 64 • N; 19 • W to 31 • W), so some spatial smoothing is to be expected. The model shows (Fig. 11a) that diatoms are the most abundant functional group, with a large bloom in early spring and a secondary bloom in fall. The diatom bloom is followed by blooms of dinoflagellates and coccolithophores. The double diatom peak is also present in the model nitrogenbased concentrations of the functional groups (Fig. 11b). The abundance of coccolithophores in the CPR data is only a fraction of what one may expect, due to the small size of the organisms (5 µm) compared to the size of the mesh used (250 µm). The CPR data are however a relevant index of the relative annual changes and thus reveal the blooming of coccolithophores. The nitrogen-based concentration of coccolithophores in the model is thus much more substantial and in line with the findings of previous field work studies in the northeast North Atlantic, reporting large concentrations of coccolithophores during the summer bloom. For example, Fernandez, et al. (1993) reported a total estimated coccolithophore C biomass as large as 50 % of the total phytoplankton C biomass.

Biogeochemical response to physical forcing
The seasonal changes in biogeochemical properties as a result of physical forcing are given in Fig. 12, which shows the model seasonal cycles of chlorophyll, calcite (PIC), SST, PAR, the ratio of euphotic depth to MLD (Zeu/MLD), and nutrients (NO 3 , PO 4 , and SiO 2 ). The rise of Chl a starts in April, when the Zeu/MLD ratio is ∼0.5 and a drawdown of nutrients ensues. That is as one may expect, according to the Sverdrup (1953) theory on vernal blooming of phytoplankton, as the critical depths in oceanic environments are several times that of the Zeu.
The PIC concentration starts to rise in May when the Zeu/MLD has reached ∼2. The broad peak in the total Chl a from all three functional groups extends from May to September, associated with shallower MLDs (Zeu/MLD ratio > 2.0). Depletion of available light and nutrients for growth and an increase of grazing pressure reduce biomass significantly during the summer. The growth season cycle is completed during autumn, as wind stress increases the vertical mixing, promoting nutrient renewal. The coccolithophore bloom, indicated by the changes in PIC concentration due to calcite production, peaks in July-August, when light conditions (Zeu/MLD ratio > 4.0) provide favorable conditions for the bloom-forming coccolithophore Emiliania huxleyi, a species well-known to the area (Balch et al., 1992;Fernandez et al., 1993;Holligan et al., 1993a;Holligan et al., 1993b). As the MLD deepens and light (PAR) levels are significantly reduced, Chl a and PIC concentrations drop gradually after September. The model confirms established ideas of the interplay of light availability, vigorous winter mixing/summer restratification, and nutrient availability, which are central to the classical North Atlantic spring bloom (Nanninga and Tyrrell, 1996;Tyrrell and Taylor, 1996).

Impact of phytoplankton blooms on carbon uptake
To evaluate the impact of phytoplankton groups on the uptake of carbon, we conducted a model experiment with the coccolithophore components turned off (no carbonate pump) and compared the results with the baseline experiment (biological and carbonate pumps), including all three functional groups (diatoms, dinoflagellates, and coccolithophores). The results are shown in Figures 13 and 14 as seasonal 8-day climatologic averages for 1998-2008. In Figure 13 we see the results of seasonal changes in the upper 120 m verticallyintegrated phytoplankton biomass concentrations with and without the presence of coccolithophores ( Fig. 13a and b,  respectively), and the corresponding changes in net community production (Fig. 13c). In the presence of coccolithophores, the yearly-averaged integrated biomass for diatoms, dinoflagellates, and coccolithophores was 534, Geosci. Model Dev., 5, 683-707, 2012 www.geosci-model-dev.net/5/683/2012/ 318, and 184 mg C m −2 , respectively. The corresponding model-derived Np with all three functional groups was 134 mg C m −2 d −1 . With the exclusion of coccolithophores, the population of diatoms increased to an integrated value of 613 mg C m −2 and dinoflagellates to 341 mg C m −2 , respectively. However, there was a net decrease in Np to 104 mg C m −2 d −1 (or a 22.4 % decrease). At first glance, it would be expected that the diatoms and dinoflagellates would benefit from the absence of coccolithophores and have additional nutrients to draw from and increase their combined production. However, the explanation lies in the carbon uptake efficiency of coccolithophores when compared to the other two groups. In equation B35 of Appendix B, it can be seen that, for the baseline run, dinoflagellates and diatoms have a Redfield C:N (6.625), while coccolithophores have a higher C:N (9.4), which is an average, based on the range of 5.81 to 13.05 reported by Fernandez et al. (1993). The diatoms and dinoflagellates do pick up some of the slack but not entirely, because the coccolithophores uptake carbon 42 % more efficiently (100[1-9.4/6.625)). In addition, the peak demand of nutrients does not overlap among the three functional groups, as the maximum productive time periods of each group are not in phase. So the absence of coccolithophores is less influential on the availability of nutrients for the other two groups. This seems to be a consistent explanation of why the removal of coccolithophores from the model results in the reported differences in PP and pCO 2 drawdown. A couple of sensitivity tests were conducted to assess the effect of changing the C:N ratio for coccolithophores in the model. We used two values of C:N ratio for coccolithophores in the model, Redfield (R1 = 6.625), and twice the Redfield ratio (R2 = 13.25). The annual mean (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) results from R1 for the following key model parameters, pCO 2 , DIC, TA, air-sea CO 2 flux, net community production, CaCO 3 concentration, and PIC production are 358.8 µatm, 2114.6 µmol kg −1 , 2308.3 µmol kg −1 , 4.4 mmol C m −2 d −1 , 123.6 mg C m −2 d −1 , 3.3 mg m −3 , and 4.9 mg C m −2 d −1 . The equivalent values in the same units for experiment R2 are 353.7 (−5.1), 2111.3 (−3.3), 2308.0 (−0.3), 5.73 (+131 %), 146.1 (+118 %), 6.7 (+200 %), and 9.8 (+200 %). So the model is very sensitive to changes in the C:N ratio for coccolithophores.
The 22 % reduction in Np between the baseline run (C:N = 9.4 for coccolithophores) and the baseline run without coccolithophores has consequences in the biological uptake of carbon. Figure 13 shows that, without the net community production and calcification by coccolithophores, at the peak of the biological drawdown, the alkalinity (Fig. 14a) increases by up to 4 µmol kg −1 with a corresponding increase in DIC (Fig. 14b) of up to 16 µmol kg −1 . As a result, the net effect of the absence of the coccolithophores bloom is an increase in pCO 2 of more than 20 µatm (Fig. 14c). The largest differences occur from July to October, which translate into pCO 2 differences ranging from 24 % to 35 % during this period. The seasonal changes in the sea-air CO 2 flux are shown in Fig. 14d. At the peak of the bloom, the exclusion of coccolithophores causes a reduction of CO 2 uptake of more than 6 mmol m −2 d −1 . On average, the impact of coccolithophores on air-sea CO 2 flux is 27 %. This result highlights the importance of including all major functional groups in the modeling of carbon variability in the subpolar North Atlantic. Some sensitivity tests were conducted using different gas transfer coefficients (k). The baseline run uses a value of 0.31 (Wanninkhof, 1992), as indicated in Eq. (B39) in Appendix B. A value of k = 0.27 (Sweeney et al., 2009) was used to evaluate the model sensitivity of air-sea CO 2 flux to different gas transfer coefficients. A reduced CO 2 ocean uptake of ∼7 % is obtained with the experiment using k = 0.27, when compared with the experiment using k=0.31, with both experiments including coccolithophores. Model experiments using k = 0.27 with and without coccolithophores indicate that the exclusion of coccolithophores causes a reduced CO 2 ocean uptake of ∼28 %, a small difference (1 %) compared to the use of k = 0.31. Suykens et al. (2010) reported the carbon uptake impact of coccolithophore blooms in the northern Bay of Biscay based on DIC and alkalinity observations. They concluded that the decrease of DIC (and increase of pCO 2 ) due to net community calcification was overwhelmingly lower than the decrease of DIC (and decrease of pCO 2 ) due to net community production, a clear indication of the importance of coccolithophore blooms in the uptake of atmospheric CO 2 in the North Atlantic, where blooms of coccolithophores are the most intense and recurrent.
A spatial assessment on the reduction in uptake of atmospheric CO 2 under the exclusion of coccolithophores was made for the regional domain of the Irminger and Icelandic Basins. The region considered is bound by 44 • W -8 • W and 56 • N -65 • N, and by depths greater than 200 m to exclude coastal regions of Iceland and Greenland. The monthly areal (in m 2 ) extent of the bloom was determined using Sea-WiFS monthly composites within the above boundaries and considering only pixels with PIC concentrations larger than 5 mg m −3 , a value arbitrarily picked as the threshold of significant influence of coccolithophore blooms in the drawdown of atmospheric CO 2 (monthly concentrations range from 0 to 27 mg m −3 ). For comparison, a threshold of 2 mg m −3 in PIC concentration increases the flux estimate by about 12 % (or about 180 tonnes C yr −1 ). The monthly CO 2 flux difference was derived from the product of satellite area and the flux difference shown in Fig. 14d. The annual mean flux difference was estimated at nearly 1500 tonnes C yr −1 .
The seasonal drawdown of surface ocean pCO 2 is a result of two competing effects, i.e., temperature warming and biological uptake effects. Takahashi et al.(2002) developed a method to separate these two effects for the global oceans. The effect of biology (B e ) on the surface ocean pCO 2 in a given area is represented by the seasonal amplitude of pCO 2 , corrected to the mean annual temperature in that area. The effect of temperature changes (T e ) on the seasonal pCO 2 variations is represented by the seasonal amplitude of the mean annual pCO 2 , corrected to the range of observed temperatures. B e = ( pCO 2 ) bio = (pCO 2 at T mean ) max − (pCO 2 at T mean ) min T e = ( pCO 2 ) temp = (pCO 2 at T obs ) max − (pCO 2 at T obs ) min (1) where the subscripts "min" and "max" indicate the seasonal minimum and maximum values, T obs is the observed temperature for each pCO 2 value and T mean is the mean temperature over the entire seasonal cycle. The relative importance of the biology and temperature effects can then be expressed by the ratio B e /T e . Using climatologic  model seasonal surface ocean pCO 2 and SST, we calculate the ratio B e /T e as 1.92, which means that the biology effect is nearly twice the temperature effect in the shaping of the seasonal pCO 2 variability at the model location.
Although it is not straightforward to infer to what extent the model results from a single location in the subpolar North Atlantic are representative of the North Atlantic in general, previous studies have shown coherence of results in much broader regions of the North Atlantic. For example, McKinley et al. (2011) have shown that the North Atlantic can be separated in three large biogeographic regions ("biomes") in terms of surface ocean pCO 2 variability. These biomes were assigned on the basis of annual maximum mixed layer depth, annual mean SeaWiFS Chl a, and SST. Namely, the North Atlantic was divided into a northern seasonally stratified gyre (SP-SS) biome, a southern permanently stratified subtropical gyre (ST-PS) biome, and a seasonally stratified subtropical (ST-SS) biome. For a 29 yr-long observational period , the surface ocean pCO 2 trends converge with the atmospheric growth for all 3 biomes, which occupy 87 % of the total area of the North Atlantic. When a shorter time scale (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) is considered, the three biomes have different trends, but the SP-SS biome itself has a significant geographic extension, including a large region south of Iceland, and the Irminger, Icelandic, and Norwegian Seas. Our 1-D model location lies within the SP-SS biome and therefore is representative of surface ocean pCO 2 decadal variability for that region. However, for shorter time scales, more significant regional differences are likely to occur.

Summary and conclusions
A 1-D ecosystem model was developed to assess the relative contribution of phytoplankton functional types to atmospheric CO 2 uptake. Using a Taylor diagram, skill assessment of model versus field measurements reveals high scores for the majority of biogeochemical parameters, for which in situ data are available.
The seasonal patterns of phytoplankton concentrations are a response to the interplay between light availability, vigorous winter mixing/summer restratification, and nutrient availability, not unlike the classical North Atlantic spring bloom (Ducklow and Harris, 1993;Weeks et al., 1993). Functional groups compete seasonally for ideal growth conditions (Litchman et al., 2007;Sambrotto et al., 1993;Sieracki et al., 1993). Model results indicate that the springsummer bloom consists predominantly of diatoms, with still significant but less intense blooms of dinoflagellates and coccolithophores. The model shows that the diatom biomass peaks in May, with a secondary and less intense bloom in September. The dinoflagellates and coccolithophores peak in July through August, during which drawdown of surfaceocean CO 2 reaches its maximum value. The effect of biological changes in the surface ocean pCO 2 exceeds the temperature effect by a factor of almost 2, a clear indication of the importance of phytoplankton photosynthesis on the uptake of atmospheric CO 2 in the region, a result that is in agreement with previous studies (Takahashi et al., 1993(Takahashi et al., , 2002(Takahashi et al., , 2009Ullman et al., 2009;Bennington et al., 2009;Metzl et al., 2010).
Model experiments were conducted to investigate the seasonal changes in phytoplankton concentration with and without the presence of coccolithophores, and their impact on carbon uptake. Without the influence of coccolithophore blooms, the alkalinity increases by almost 4 µmol kg −1 and DIC is elevated by up to 16 µmol kg −1 . The net effect of coccolithophores blooms is a drawdown in pCO 2 of up to about 20 µatm during summer, with a corresponding increase of atmospheric CO 2 uptake of about 6 mmol m −2 d −1 , an indication of the importance of including all major phytoplankton functional groups when modeling the biological carbon pump variability in the subpolar North Atlantic.

Ice-ocean model description (GPOMZ)
The ocean model is hydrostatic and Boussinesq and uses a generalized vertical coordinate system, as described by  with a modified scalar advection scheme to avoid overshooting at sharp fronts (Mauritzen and Häkkinen, 1997). The equation of state is formulated in terms of in situ density (Mellor, 1991), expressed as a function of potential temperature, salinity, and pressure. The model's prognostic variables are the horizontal velocity components, potential temperature and salinity, and twice the kinetic energy, q 2 and q 2 l, where l is the turbulence macroscale (Mellor andYamada, 1974, 1982). These turbulence quantities, together with the vertical velocity shear and buoyancy, determine the vertical mixing coefficients for momentum and scalar variables. The dynamic-thermodynamic ice model is coupled to the ocean model via interfacial stresses and via salinity and heat fluxes through the ice-water interface. The thermodynamic component of the ice model follows Semtner's (1976) formulation but with modifications to account for leads. The  dynamics of the ice model is described as a continuum with a generalized viscous rheology (Häkkinen, 1987;Häkkinen and Mellor, 1992).
The ocean model was initialized with the winter average hydrographic climatology derived from the Polar Science Center Hydrographic Climatology (PHC 3.0) (Steele et al., 2001). A transport of 0.8 Sv was specified at Bering Strait and 0.8 Sv out at the southern boundary. At the northern and southern boundary, the salinities and temperatures are relaxed to monthly climatological values. Climatological monthly river discharge for 474 rivers within the model domain is specified with the global river database, described by Barron and Smedstad (2002). An additional river discharge of 700 km 3 yr −1 is uniformly distributed along the Eurasian margin of the Arctic. River runoff is treated as a virtual flux of salt. The model is forced with daily NCEP/NCAR Reanalysis sea level pressure, air temperature 2 m above sea level, surface momentum fluxes and relative humidity for 1948-2011. Climatological monthly mean cloudiness is based on the ISCCP D2 data set from 1984-2004, except north of ∼65 • N where a uniform monthly mean cloudiness is specified. Climatological monthly mean precipitation minus evaporation (P−E) is derived from the NCEP operational analysis (Rasmusson and Mo, 1996). An additional 1350 km 3 yr −1 of precipitation is added, and distributed with a peak value at ∼8 • N to simulate the Intertropical Convergence Zone.

Biogeochemical model description (ECO1D-2.0)
The model features multiple functional groups (diatoms, dinoflagellates, and coccolithophores), zooplankton, nutrients (NO 3 , PO 4 , NH 4 , SiO 2 , and Fe), POC, DIC, DOC, alkalinity (TA), calcite production, chlorophyll, complete carbonate chemistry, and air-sea CO 2 flux. The governing equations for the biogeochemical model are provided hereafter, where the subscripted index i = 1, 2, 3 represents diatoms, dinoflagellates, and coccolithophores, respectively. Iron (Fe) limitation, although included in the model, was not considered in this study. However, there is evidence (Nielsdóttir et al., 2009) of iron limitation of the post-bloom (July to early September) phytoplankton communities in the Iceland Basin, east of the model location, where high nutrient-low chlorophyll (HNLC) conditions may occur. Iron limitation studies in the Iceland Basin will be a topic of future studies using the same model.
The basic form of the model equations for a generic tracer of concentration C, differentiated with respect to time t and vertical coordinate z and balanced by the sources and sinks of the particular biogeochemical property, is: where w is the vertical velocity and k v is the vertical eddy diffusivity. We simplify the notation of the left hand side of all equations by substituting it with the single term dC dt .

B13 Irradiance model
The total (infrared plus visible) solar radiation is obtained using the Frouin model (Frouin et al., 1989). This model provides the total radiation and the photosynthetically available radiation (PAR). The infrared (I IR ) component is obtained by subtracting the PAR component from the total solar radiation. Using a spectral model for PAR (Gregg and Carder, 1990), the spectral PAR componentI PAR (λ)can be determined. The infrared component (for mixed layer model only) and the PAR component (mixed layer and biogeochemical model) of the penetrating irradiance are obtained from I PAR (λ, z) = I (λ, z − z) exp[−(a w (λ) + a ph (λ)) z] (B26) where a IR (3.75 m -1 ) is the attenuation coefficient for infrared radiation and a w (λ) and a ph (λ) are the wavelengthdependent light attenuation coefficients for water and phytoplankton, respectively. The water and chlorophyll-dependent attenuation coefficients from Morel (1988) were used in the model for this study. The dissolved matter attenuation coefficients, a dm (λ), are calculated by applying the IOP (inherent optical properties) model of Garver and Siegel (1997), which uses water leaving radiances from 6 SeaWiFS bands as input (level 3 binned monthly composites). The IOP model calculates the attenuation coefficient due to dissolved matter for the 443 nm wavelength. The attenuation coefficients for other wavelengths are obtained from a dm (λ) = a dm (443) exp [S(λ − 443)] where the exponential decay constant, S, is chosen to be 0.02061. A correction is applied to the downward irradiance pathway to account for seawater light refraction following Snell's law. After some algebraic manipulations, the correction is applied to zas follows: where α z is the solar zenith angle, and n s is the seawater refraction coefficient, which is expressed as a function of salinity and temperature adapted from Table 3.12 of Neumann and Pierson (1966) as The C:N ratio for coccolithophores (Coc) of 9.4 is the average from the reported range of 5.81 to 13.05 in Fernandez et al. (1993).

B17 Chlorophyll and Chl:N ratio
The model chlorophyll is calculated following the photoadaptation scheme for Chl:N ratio of Doney et al. (1996): The subscripts phy, zoo, and det refer to phytoplankton, zooplankton, and detritus, respectively. Table B1 defines the  model state variables and Table B2 provides the definition of the ECO1D-2.0 parameters and values used. Equivalent parameters for ECO1D-1.0, configured for and applied at the BATS site, are shown in parenthesis for comparison.

B18 Model forcing and relaxation approach
The terms δ(z)FCO 2 /ρ and δ(z)FO 2 /ρ, in (B15) and (B19), respectively, represent the CO 2 and O 2 sea-air fluxes at the surface. The Kroenecker delta (δ[z=0]=1; δ[z >0]=0) is used to denote that carbon dioxide and oxygen fluxes (FCO 2 and FO 2 , respectively) are only applied at the sea-air interface. The following formulations for the CO 2 and O 2 gas transfer were applied in the form of flux boundary conditions (FCO 2 and FO 2 in mmol m −2 yr −1 ) at the sea-air interface: where, K o is the gas transfer velocity, in m d −1 , which is a function of water temperature and wind speed (Wanninkhof, 1992), α is the CO 2 solubility in seawater (in mmol m −3 µatm), which is a function of temperature and salinity (Weiss, 1974), pCO 2 (in µatm) is the difference between sea and air pCO 2 , and O * 2 is the oxygen saturation concentration (in mmol m −3 ) in seawater, which is a function of temperature and atmospheric pressure (Weiss, 1970).
We adopt the following relationship between gas transfer and wind speed (W ) (Wanninkhof, 1992) using the NCEP 3 hourly winds: where Sc is the Schmidt number of CO 2 or O 2 (Wanninkhof, 1992).
To account for horizontal advective processes of heat and salt within deeper layers of the 1-D mixed layer model, temperature and salinity are assimilated from the 3-D model using a straightforward approach. The approach consists of relaxing the temperature and salinity profiles calculated by the 1-D mixed layer model to the pre-calculated values provided by the 3-D model at all depths. The assimilation of T and S is done using a Newtonian relaxation (nudging) method (Bauer and Wulfmeyer, 2009) with a relaxation time scale (τ ) of 10 days for T and 30 days for S. For model properties where the relaxation time scales for deep water are short (10 days or less), the term "nudging" is not applicable, as the properties are set to observed values.
A similar relaxation approach is used for nitrate, phosphate, and silicate, except that the nutrient values originate from T -dependent equations obtained from T , NO 3 , PO 4 , and SiO 2 climatologic monthly profiles (0-500 m) from the World Ocean Atlas 2005 at the model site. The total number of data points is 168 (N=14 depths×12 months).