Factors controlling interannual variability of vertical organic matter export and phytoplankton bloom dynamics – a numerical case-study for the NW Mediterranean Sea

Mid-latitude spring blooms of phytoplankton show considerable year-to-year variability in timing, spatial extent and intensity. It is still unclear to what degree the bloom variability is connected to the magnitude of the vertical flux of organic matter. A coupled three-dimensional hydrodynamic-biogeochemical model is used to relate interannual variability in phytoplankton spring-bloom dynamics to variability in the vertical export of organic matter in the NW Mediterranean Sea. Simulation results from 2001 to 2010, validated against remote-sensing chlorophyll, show marked interannual variability in both timing and shape of the bloom. Model results show a tendency for the bloom to start later after cold and windy winters. However, the onset of the bloom occurs often when the mixed layer is still several hundred metres deep while the heat flux is already approaching zero and turbulent mixing is low. Frequency and intensity of wind episodes control both the timing and development of the bloom and the consequent export flux of organic matter. The wintertime flux is greater than zero and shows relatively low interannual variability. The magnitude of the interannual variability is mainly determined in March when the frequency of windy days positively correlates with the export flux. Frequent wind-driven mixing episodes act to increase the export flux and, at the same time, to interrupt the bloom. Perhaps counterintuitively, our analysis shows that years with discontinuous, low-chlorophyll blooms are likely to have higher export flux than years with intense uninterrupted blooms. The NW Mediterranean shows strong analogy with the North Atlantic section within the same latitude range. Hence, our results may also be applicable to this quantitatively more important area of the world ocean.


Introduction
The dynamics regulating the vertical flux of organic matter in the ocean determine the partitioning of carbon between surface and deep ocean and the transfer of organic matter to higher trophic levels.These dynamics affect both climate and the ocean's ability to sustain fisheries.Primary production of organic matter occurs all year round, with a seasonal bloom at all mid-latitude oceans.These blooms are an important component of the total CO 2 ocean uptake (Takahashi et al., 2009) which is achieved through the export of organic matter to the deep ocean.Blooms show considerable variability in time, space and intensity, and this variability is related to physical processes that affect vertical mixing.
Timing and intensity of the bloom show latitudinal and interannual variability mainly determined by the variability in atmospheric forcing (Henson et al., 2009;Ueyama and Monger, 2005;Waniek, 2003).Henson et al. (2009) showed that the onset of the bloom in the North Atlantic, between 40 • N and 45 • N, can vary from year to year by as much as 20 weeks.This area represents the transition between subpolar light-limited and subtropical nutrient-limited environments (Dutkiewicz et al., 2001).The NW Mediterranean Sea is enclosed between these latitudes and its seasonal bloom shows high variability concurrent with its latitudinal counterpart in the North Atlantic Ocean.The NW Mediterranean bloom has been studied using remote-sensing chlorophyll.
The depth of the mixed layer is commonly related to the onset of blooms according to the "critical depth hypothesis" (Sverdrup, 1953).Since the hypothesis assumes phytoplankton to be homogeneously distributed over the mixed layer, the mixed layer regulates phytoplankton mean exposure to light.The bloom develops as soon as the mixed layer, shoals and becomes shallower than a critical depth, such that vertically integrated phytoplankton growth wins over phytoplankton losses.According to this theory, the interannual variability in the timing of the bloom would be the result of the interannual variability in the timing of re-stratification.For the Irminger basin (NE Atlantic), Henson et al. (2006) described a preconditioning effect of the winter atmospheric forcing on the bloom timing: a deeper mixed layer would take longer to shoal up to the critical depth, resulting in a later start date for the bloom.
The critical-depth hypothesis has been questioned as a predictor of the onset of the bloom, after growing evidence of blooms taking place in deep mixed layers (Townsend et al., 1992;Eilertsen, 1993;Dale et al., 1999;Koertzinger et al., 2008).Huisman et al. (1999) used a turbulent diffusion model to show that a bloom can develop if turbulent mixing is less than some critical value, regardless of the depth of the mixed layer.Recently, Taylor and Ferrari (2011) related this critical turbulent diffusivity to the atmospheric forcing, showing that when cooling subsides (heat flux is close to zero) turbulent mixing becomes weak, increasing phytoplankton residence time in the euphotic layer and allowing blooms to develop even in the absence of stratification.The authors focused their analysis on thermally driven convection pointing out, however, that their results can be extended to scenarios with turbulence generated by wind forcing and evaporation.
Several authors have pointed out the importance of windinduced mixing during the bloom period as the key forc-ing process contributing to the interannual variability in both timing and intensity of the bloom (Ueyama and Monger, 2005;Waniek, 2003).Nevertheless, to which extent the bloom variability is connected to the variability of the vertical flux of organic matter remains under discussion.Windinduced mixing episodes are expected to have an impact on the export flux as high biomass surface water is mixed with low-biomass underneath water.The effectiveness of vertical mixing in exporting organic matter to the depths has been described before.Ho and Marra (1994) ascribed a significant part of the Northeast Atlantic export of primary production to intermittent early spring vertical mixing.Bishop et al. (1986) showed how variations of the mixed layer in a warm-core ring were able to remove up to 67 % of primary production.Gardner et al. (1995) described a day-night "mixed-layer pump" as an important mechanism to sustain new primary production and to remove particles from surface waters.Koeve et al. (2002) showed how storm-induced variations of the mixed layer in the Northeast Atlantic were able to interrupt the spring bloom and transport particles to depth through convective mixing.
Based on numerical simulations validated against field data, we relate the interannual variability in time and intensity of the bloom to the interannual variability in the export flux of organic matter in the NW Mediterranean Sea.We hypothesize that wind-induced mixing episodes during the bloom are responsible for shaping both the bloom development and the export flux by effectively redistributing phytoplankton at depth.The interannual variability in the frequency and intensity of wind forcing during the bloom would, thus, account for a significant part of the interannual variability in the bloom characteristics and in the vertical flux of organic matter.We evaluate the effectiveness of the closeto-zero heat flux, as proposed by Taylor and Ferrari (2011), as predictor for the onset of the bloom in our numerical simulation and stress the role of wind mixing in driving both the heat flux and the distribution of phytoplankton at depth.The present numerical study is based on a newly implemented hydrodynamic model setting for the Western Mediterranean Sea, coupled to a newly developed biogeochemical model (Bernardello, 2010).Model validation is performed with reconstructed mean dynamic topography, argo float data and remote-sensing chlorophyll.

The model
The hydrodynamic component of this model is the Stony Brook Parallel Ocean Model (Jordi and Wang, 2012), a parallelised version of the Princeton Ocean Model (POM) (Blumberg and Mellor, 1987).The model domain includes the Western Mediterranean Sea between the Atlantic side of the strait of Gibraltar and the Sicily channel.The horizontal resolution is 1/20 • so that the mesh size is constant in longitude (5560 m) and decreases northwards (from 4456 m to 3964 m).In the vertical dimension the grid is resolved by 52 sigma-layers unevenly distributed with higher resolution at the surface and at the bottom.The bottom topography is obtained from the ETOPO1 1 arc-minute global relief model (Amante and Eakins, 2009), after bilinear interpolation.The atmospheric forcing is prescribed by using archived forecast analysis data provided by the European Centre for Medium Range Weather Forecast (ECMWF) with a spatial resolution of 0.25 • and a time-step of 6 h.These data are used as inputs to a set of bulk-formulas used to represent air-sea boundary processes (see Estournel et al., 2009, as an example).Daily outputs of temperature, salinity and velocities from the Mediterranean Forecasting System (MFS) high resolution model are used as lateral boundary conditions.Data come from the implementation MFS1671 with a horizontal resolution of 1/16 • and 72 unevenly spaced vertical levels, based on the numerical code OPA8.1 (Tonani et al., 2008(Tonani et al., , 2009;;Pinardi and Coppini, 2010).The MFS1671 fields are used for fluxes directed into our domain while the model computes the outward fluxes.
The model is initialised at rest using MEDAR-MEDATLAS climatology for temperature and salinity.A spin-up of 13 yr is performed.During this phase the model is forced by repeating the same year for the atmospheric forcing and the open lateral boundaries.Both datasets are obtained by averaging ten years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) of ECMWF and MFS1671 data, respectively.At the end of the 13th yr the interannual simulation starts by prescribing the interannual forcing for the period 2001-2010.In order to avoid drift in the simulation, during the spin-up and the period 2001-2010, temperature and salinity fields are restored to the climatology values with a timescale of 30 days.This timescale is short enough to control the drift of the simulation for the period studied and is long enough to leave most of the surface interannual signal untouched.
A newly developed biogeochemical component is linked to the above physical module.This is an aggregate-type model based on previous work (Fasham et al., 1990;Varela et al., 1992;Bahamón and Cruzado, 2003).It consists of different compartments representing nitrate, ammonium, phytoplankton, bacteria, zooplankton, detritic matter and dissolved organic matter and uses nitrogen as currency.In order to better represent the vertical flux of organic matter, particulate matter is split into small and large detritus and dissolved organic matter into labile and semi-labile fractions.Large and small detritus and phytoplankton are vertically redistributed at each time-step assuming sinking rates of 210 m, 1.5 m and 0.8 m per day, respectively.Some degree of complexity is added by allowing variable C:N ratios in both detritus and dissolved organic matter.Bacterial processes are modified from Fasham et al. (1990) by introducing a bacterial stoichiometric sub-model (Anderson, 1992) that depends on the variable C : N ratio of labile dissolved organic matter.Also, a temperature dependent limitation is introduced for zooplankton growth following Simonot et al. (1988).
The biogeochemical component is initialised using nitrate and chlorophyll vertical profiles compiled by Manca et al. (2004) using EU/MEDAR/MEDATLAS II and MATER databases.The rest of the variables are initialised at low values.The biogeochemical spin-up starts on the 9th yr of the hydrodynamic spin-up and runs for four years.During this phase and the following interannual simulation, the open lateral boundary conditions are prescribed from seasonal climatological values for nitrate and chlorophyll.

The area of the bloom
In order to avoid coastal dynamics, we limit the analysis to the area of the bloom defined as the set of pixels where chlorophyll concentration is higher than 0.5 mg m −3 , in the spring MODIS-Aqua climatology, and bottom depth is deeper than 200 m in the model topography.This area (Fig. 1) is similar to the bloom area identified by D'Ortenzio and d'Alcala ( 2009) from K-means cluster analysis of SeaW-iFS remote-sensing chlorophyll, confirming the consistency of our criteria.
The bloom area roughly coincides with the center of the general cyclonic circulation that characterises the NW Mediterranean Sea (Millot, 1999).On the northern boundary, the northern current flows along the slope in the SW direction.This current forms from the joining of the Western and the Eastern Corsican currents flowing northwards on both sides of Corsica (Taupiere- Letage and Millot, 1986).The southern boundary of the gyre is less well defined because of the intense mesoscale activity that characterises the North Balearic front.This front separates the northern cold, salty and older Modified Atlantic Water (MAW) from the southern warm, fresher and younger MAW.
The area selected is known as a region where deep wintertime mixing and deep-water formation take place.The cyclonic circulation leads to an uplift of the isopycnals forming a dome centred on the south of the Gulf of Lions at 42 • N 5 • E. During winter the Mistral from the Rhone Valley and the Tramontana from the north side of the Pyrenees blow over the area.These winds are cold and dry and typically occur in strong bursts, lasting for a few days, that are able to erode the near-surface stratification and expose the weakly stratified waters underneath (Leaman and Schott, 1991;Mertens and Schott, 1998).Deep convection can then occur with a consequent nutrient enrichment of the upper layers that will fuel the next spring bloom.However, the deep-water formation process is very irregular and can be completely absent during some years (Mertens and Schott, 1998).This variability is likely to be reflected also in bloom dynamics.

Data treatment
The present analysis is focused on the period from 1 December to 31 May from 2001 to 2010 (9 yr), as we are interested in the interannual variability of spring bloom dynamics.However, as the Aqua-MODIS data series starts in June 2002, this 9-yr period is reduced to 8 yr for the purpose of the satellite-model data comparison.

Hydrodynamic model results are validated by comparison between modelled Mean Dynamic Topography (MDT)
with associated geostrophic circulation and data-model reconstruction.MDT is the sea elevation due to the mean oceanic circulation.The only MDT available for the Mediterranean Sea (hereafter called RioMDT) was reconstructed by combining oceanic observations from altimetry and in situ measurements and outputs from an ocean general circulation model with no data assimilation for the period 1993-1999 by Rio et al. (2007).The bloom area defined for this study is superimposed on both MDTs as a reference for the discussion (Fig. 1).
Table 1.Mean Mixed Layer Depth (MLD), mean daily vertically integrated nitrogen phytoplankton biomass (0-75 m, PHY) and primary production (0-75 m, PP), mean daily organic nitrogen export flux at 75 m depth (EF), mean daily neat heat flux (HF), mean daily wind speed (WS), day of the maximum MLD, week of the bloom onset for the model and for MODIS. 2002MODIS. 2003MODIS. 2004MODIS. 2005MODIS. 2006MODIS. 2007MODIS. 2008  To validate the model estimate of the Mixed Layer Depth (MLD) we use Argo floats data obtained from the Coriolis Data Assembly Center.The float data are part of the MedArgo Programme (Poulain et al., 2007) started in 2003.We select float profiles for our specific bloom area and for the period of the simulation.The MLD is calculated (in both data and model) as the depth at which potential temperature or potential density vary by more than T = 0.2 • C or σ θ = 0.03 kg m −3 relative to their values at 10 m depth as in de Boyer Montegut et al. (2004).The model output is sampled at the location and day of each valid Argo-float profile to allow model-data comparison (Fig. 2).
We use chlorophyll level-3 8-day composite maps from sensor Aqua-MODIS (2009.1 Ocean Colour Reprocessing) obtained from the NASA Ocean Colour Home Page.Images are interpolated to the model grid because their resolution (∼1/16 pixel km −2 ) is higher than that of the model (∼1/23 pixel km −2 ).The data series obtained spans from June 2002 to December 2010 for a total of 391 maps.The surface chlorophyll estimated by the model is averaged at 8-day intervals to match the MODIS data series.Finally, MODIS data and model estimate are averaged over the area of the bloom and compared (Fig. 3).
We calculate the maximum range of variability in the timing of the bloom following Henson et al. (2009).First, for each grid element, the onset of the bloom is determined for each year and for both series (model and MODIS).To this end, the method proposed by Siegel et al. (2002) is modified by considering the onset of the bloom as the first week of the year when the chlorophyll concentration is 10 % (instead of 5 %) above the median value for the period January-May (instead of the whole year).We obtain a map for each year (and both series) representing the onset date of the bloom.For each pixel the maximum and minimum values across the eight years are selected to define the maximum range of variability in the bloom onset date over the period 2003-2010.These maps are shown in Fig. 4 with the area of the bloom superimposed for comparison.In Table 1 we present a synthesis of model results that are used throughout the analysis.MLD is obtained from daily averages of temperature and salinity fields.The net heat flux (HF) represents the sum of shortwave and longwave radiation, latent and sensible heat assuming the positive sign for heat gained by the ocean.Wind speed (WS) is the magnitude of the resultant of the meridional and zonal components of the wind forcing applied to the model.MLD, HF and WS are presented as the average over the bloom area of daily mean values for each year.We present also the day of the maximum MLD and the week of the bloom onset (referred to 1 January) for both the model and MODIS data.These dates are calculated also from the spatial average.
Model data for primary production (PP), export flux (EF) and phytoplankton biomass (PHY) are reported in terms of nitrogen as daily total (PP and EF) and daily mean (PHY).PP and PHY are calculated as depth-integrated values in the layer between the surface and 75 m depth while EF is considered at 75 m depth.EF includes contributions from zooplankton, detritus, dissolved organic matter, phytoplankton and bacteria.It takes into account the vertical transport as a result of advection, diffusion and, for phytoplankton and detritus, gravitational sinking.Sinking can contribute only positively to the export flux while advection and diffusion can operate also upwards, giving rise to negative contributions.In this sense, EF is a measure of the daily net vertical export flux of organic nitrogen at a fixed depth.We choose 75 m depth to characterise the export from the euphotic zone.
The MLD during winter is almost always deeper than 75 m, thus, resuspension of organic nitrogen is very likely.However, contributions to EF from resuspension are taken into account in the definition of EF itself.Another reason for the choice of the EF reference depth is that we are interested in mixing events that operate at daily timescale.Such mixing events are primarily wind-induced and their typical lengthscale in the area during winter is ∼80 m according to the Obukhov length relation.

Results
The two MDT reproduce the main cyclonic circulation in the northern sector, roughly coinciding with the bloom area defined for this study (Fig. 1).In both model and data, the northern current originates in the NE sector after the joining of the eastern and western Corsican currents as described by Taupiere-Letage and Millot (1986).In the RioMDT the northern current is clearly visible along the northern coast from Italy up to the Ibiza channel while in the model MDT it appears somewhat more intense on the Italian coast and starts decreasing along the Spanish coast, in agreement with current measurements of ∼5 cm s −1 off the Ebro delta (Font et al., 1995).From climatological studies Font et al. (1988) described the deflection of one branch of the northern current that would then return cyclonically to the northeast to form the Balearic current.Both our model and data represent this deflection and the resultant Balearic current flowing towards the center of the Algero-Provenc ¸al Basin.Here, between the Balearic Islands and Corsica the mean flow is not well defined as this area is characterised by strong mesoscale circulation across the North Balearic frontal zone.The model depicts a dominant NE flux that eventually joins the northward eastern Corsican current and closes the cyclonic circulation around the bloom area.
Model MLD is reasonably validated by ARGO float data.An example of this validation is shown in Fig. 2, for years 2005 and 2008.The model tends to overestimate the MLD throughout the winter in both years.Argo data and model estimates are significantly, though weakly, correlated in both cases (2005, r = 0.33, p < 0.05; 2008, r = 0.23, p < 0.05).There is no clear difference in the mixed layer evolution between the two years for both data and model.However, Argo-float data are sparse in space and the drift between subsequent profiles (up to 30 km) means they rarely profile the same water column twice (Smith et al., 2008).Furthermore, Argo floats drift with the prevailing currents and are found primarily in the southern and northern part of the bloom area during 2005 and 2008, respectively, along the path of the main cyclonic circulation.The northern sector is more directly exposed to the wind forcing than the southern side, suggesting that the representation of the mean mixed layer for the whole bloom area obtained from these locations is likely to be negatively/positively biased respectively in 2005/2008.
The absolute values of the MLD predicted by the model are in good agreement with those observed in past studies of deep water formation.In particular, the absolute magnitude of the marked interannual variations predicted by our model in both MLD and HF (Table 1) agree with those reconstructed by observations and numerical modelling by Mertens and Schott (1998) (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010).The bloom onset date is calculated as the first day that meets the chosen criteria for the bloom onset within the week identified as the first week of the bloom (see Sect. 3, "Data treatment").
Gulf of Lions that roughly coincides with an area of well documented deep-water formation, often referred to as the MEDOC area (MEDOC, 1970).In this and other successive studies (Gascard, 1978;Marshall and Schott, 1999)  Surface chlorophyll estimated by the model and by MODIS averaged over the bloom area is shown in Fig. 3.The timing and intensity of the spring bloom is well captured by the model (r = 0.76, p < 0.01) although the yearly maximum concentration is systematically underestimated.The interannual variability that determines the differences in the shape of the bloom seems to be reasonably captured by the model as the modest discontinuous blooms of 2003, 2007 and 2008 are clearly distinguishable from the shorter, more intense blooms of 2005, 2006 and 2010.In general, the model tends to overestimate surface chlorophyll during December, simulating a fall peak that is not observed in MODIS data.
The maximum range of variability in the timing of the bloom is shown for MODIS and model chlorophyll in Fig. 4. The bloom area is superimposed in both maps and shows, in the case of MODIS data, a strict coincidence with an area characterised by variability higher than 7-8 weeks.The model depicts a similar pattern, though more heterogeneous, with higher values (∼10 weeks) at the NE and SW sides and lower values (∼5-6 weeks) in the centre of the bloom area.When chlorophyll is averaged over the bloom area, the maximum variability in the timing of the bloom is lower for both series as can be observed in Table 1.In this case, the agreement between MODIS and simulated data is generally good except for 2009 when a transient increase in MODIS chlorophyll, on the third week, is wrongly interpreted as the onset of the bloom by our chosen criteria.
MLD and HF show evident interannual variability with years characterised by a less severe heat loss having shallower MLD and vice versa (Table 1, Fig. 5).Cold years are characterised by a late bloom onset and low mean phytoplankton biomass (PHY) and primary production (PP) with respect to warm years.These magnitudes are significantly correlated with the mean heat flux while the export flux (EF) does not show the same behaviour as can be observed in Fig. 5. Higher than average values of EF occur during both warm (i.e., 2008) and cold (i.e., 2006) years as well as for lower than average EF values (i.e., 2007 and 2005).We example of interannual variability.The relative differences between the two years in the mean MLD and heat flux are 57 % and 56 %, respectively.These two years also show very different bloom dynamics (Fig. 3).
The PP is fairly constant throughout the winter and tends to increase during the spring-bloom as can be observed for years 2005 and 2008 in Fig. 6.In both years (2005 and 2008) the peaks of EF are always associated to a sudden deepening of the mixed layer caused by a wind episode and the associated heat loss (Figs. 6 and 7).The vertical mixing redistributes the upper water column properties at depth resulting in a net downward transport of organic matter.This seems to affect the PP only partially, as PP never goes below a background value close to the winter average.The model estimates a mean winter PP (December-February) ranging from 0 Morán and Estrada (2005) in winter for the NW Mediterranean.
The simulated export flux are comparable to sediment trap data collected in the bloom area.Lee et al. (2009) reported a mean daily export of organic carbon of ∼29 mg C m −2 d −1 for the period March-May 2003 at 238 m depth at the DY-FAMED station in the northeastern sector of the bloom area.The model estimates a mean EF of particulate organic carbon of ∼40 mg C m −2 d −1 at the same depth, location and for the same period.Data from a neighbouring area, but from sediment traps deployed just below the euphotic zone (80 m) were presented by Miquel et al. (1994) for the pe-riod between mid April 1987 and mid November 1988.The authors reported for the second half of April 1987 a mean export of ∼4 mg N m −2 d −1 .The model estimate for the same period averaged over the 10 yr of the simulation gives In Fig. 8 we show the cumulative export flux for years 2001 to 2010, the 10 yr mean, and the maximum range of interannual variability calculated as the difference between the year with the highest value and the year with the lowest value (at each day).The cumulative EF grows throughout the winter and stabilises in April as vertical mixing decreases.The final maximum range of interannual variability is around 800 mmol N m −2 (∼66 % of the final mean accumulated EF), half of which is accumulated during December and the other half during March.Although the cumulative EF keeps growing during January and February, the maximum range of interannual variability is almost constant during this period (dotted line, Fig. 8).We focus the discussion on the export flux, excluding the contribution from December because this is likely to be affected by the model overestimation of phytoplankton biomass during this period.

Discussion
Cold years tend to have a deeper mixed layer and a later bloom onset date.This relationship seems at first to agree with the winter preconditioning effect postulated by Henson  (2006) in the sense that a deeper mixed layer would take longer to re-stratify and thus would result in a later bloom.However, pulses of high primary production have often been observed in absence of thermal stratification (Townsend et al., 1992;Eilertsen, 1993;Dale et al., 1999;Koertzinger et al., 2008) demonstrating that a shallow mixed layer is not a necessary condition for a bloom to start.Huisman et al. (1999) and more recently Taylor and Ferrari (2011) related the beginning of the bloom to a critical threshold of turbulent diffusivity coinciding with the decrease of the heat loss to the atmosphere.Crucial to their analysis is the distinction between a mixed layer with a uniform density and a mixing layer with a uniform density and active turbulence (Brainerd and Gregg, 1995).In other words, there is a timelag between the shutting off of the heat loss from the sea surface and the re-stratification.During this time the low active turbulence allows an increase in the average time of exposure of phytoplankton to light even if the mixed layer is still deep.Our results confirm that a condition of close-to-zero heat flux is indeed a better estimator than the mixed layer depth for the onset of the bloom (Fig. 9).Also the onset of the bloom is clearly associated with a significant drop in the vertical diffusivity which depends on the turbulent kinetic energy, as posited by Taylor and Ferrari (2011).Although we do not explicitly calculate the Sverdrup's critical depth, Fig. 9 shows how the bloom can start even when the mixed layer is still several hundred metres deep as long as the heat flux approaches zero and the vertical diffusivity is low.Since the timing of the bloom is independent of the depth of the mixed  decrease in wind together with an increase in solar radiation result in a sign inversion in the heat flux and a marked drop in the turbulent diffusivity.The same mechanism occurs in 2008 but 3 weeks earlier.
In our example, winter 2005 is colder than winter 2008 and a higher portion of the winter is characterised by intense wind forcing, higher heat loss and hence a deeper mixed layer.
In winter 2008 a less intense storminess both decreases the mean MLD and increases the probability of an early onset of the bloom.We suggest then that a deeper mixed layer and a later bloom are both consequences of a long and intense winter forcing, but are not directly related as the maximum MLD has no direct influence on the timing of the bloom.
Our model results suggest that the frequency of wind episodes regulates the intensity of blooms (Figs. 6 and 7).In general, weather conditions tend to improve with the transition from winter to spring.However, as pointed out by Waniek (2003), this is not a smooth transition and passing weather systems may interrupt the development of the bloom by mixing the phytoplankton at depth.For example, in 2008 (Fig. 6) the spring bloom is interrupted by two main wind episodes at the beginning and at the end of March.These two events act as a control on the maximum biomass accumulation at the surface, by physically transporting phytoplankton at depth and by decreasing their average light exposure, thus, limiting growth.By contrast, a relatively calm weather period following the onset of the 2005 spring bloom allows phytoplankton to grow undisturbed and reach higher concentrations than in 2008.In conclusion, we posit that wind-induced mixing during the bloom period is a key forcing agent contributing to interannual variability in both the intensity and timing of the blooms, in agreement with the satellite chlorophyll analysis performed by Ueyama and Monger (2005) for the North Atlantic.
Our model results indicate that export flux is driven by a combination of available surface biomass and vertical mixing.During January and February, EF stays at a background level, very close to PP indicating that mixed conditions are favourable to an efficient transport of organic matter at depth.No interannual variability is added during these two months because those years characterised by more intense winter mixing (i.e., 2005) tend to have lower surface biomass than those years characterised by a shallower winter mixed layer (i.e., 2008).Since EF is driven by a combination between biomass availability at the surface and vertical mixing, the result is a compensation that leads to an almost year-to-year constant winter EF.This compensation partly explains why EF is not correlated with HF in Fig. 5.
About half of the maximum range of interannual variability in EF is accumulated in March (Fig. 8).In particular, interruptions of the bloom as observed in March 2008 are associated to the peaks in EF that are responsible for this half of the interannual variability.We have found that the portion of interannual variability in EF, accumulated in March, is determined by the interannual variability in wind forcing during this month.The frequency of windy days during March, calculated as the percentage of days with average wind above the 10-yr average mean daily value for the same month, correlates well (r = 0.9, p < 0.01) with the March mean daily EF (Fig. 10).
Therefore, March is a key period during which a significant part of the interannual variability in the total EF is determined.This is because the range of variability in the date of the bloom onset is centred on the first week of this month (week 8.4, average blooming week in Table 1) amplifying the effect of the interannual variability of wind-induced mixing episodes on EF.In other words, the bloom has already started during, at least, a part of March, every year.This means that there is abundance of surface biomass to be exported at depth by eventual wind-driven mixing episodes.
Importantly, our results suggest that satellite chlorophyll could be used to identify years with higher EF.Counterintuitively, years with a bloom of intermittent nature characterised by low maximum chlorophyll concentrations are likely to also be years of high EF while years characterised by an intense uninterrupted bloom would have lower EF.
The bloom area considered in this study lies between 40.5 • N and 44 • N. Henson et al. (2009) defined a transition zone in the North Atlantic between 40 • N and 45 • N, based on the significance of the temporal correlation between weekly SeaWiFS chlorophyll concentration and MLD estimated from optimally interpolated Argo float data.North/south of this transition zone the correlation was negative/positive as expected for a subpolar (light limited)/subtropical (nutrient limited) environment.The same correlation for our model results (not shown) shows a variable (∼0-0.5),primarily negative correlation in the bloom area.This suggests that this area encompasses a variety of behaviours ranging from the southern subpolar up to the tran-sitional.We, therefore, emphasise the relevance of our results also for an important portion of the North Atlantic where high interannual variability has been observed in the timing and intensity of the spring bloom (Henson et al., 2009;Ueyama and Monger, 2005).
Shifts in storm tracks and modes of atmospheric circulation are responsible for the interannual variability in the passage of weather systems.The transition to springtime is a critical moment when the passage of these weather systems is able to modulate bloom dynamics and consequently the vertical export of organic matter.

Conclusions
Model results show a clear connection between the year-toyear variability of the spring phytoplankton bloom and the year-to-year variability of the vertical export flux of organic matter.Late blooms and deep mixed layers are both the result of intense wind forcing, but are not necessarily related to each other.A condition of close-to-zero surface heat flux as postulated by Taylor and Ferrari (2011) seems to be a better predictor for the bloom timing than the onset of stratification.The bloom onset is centred, on average, on the first week of March with earlier blooms starting during low-wind periods in February.The frequency of wind episodes is a key factor controlling the onset and the following development of the bloom.The passage of weather systems in the transition from winter to spring can interrupt the development of the bloom by actively mixing phytoplankton at depth.This results in a net vertical export of organic matter, driven by a combination of surface biomass availability and vertical mixing.Since the bloom has already started during at least some part of March each year, the frequency of wind-driven mixing episodes during this month has a strong impact on the export flux accumulated in winter-spring.Our simulation shows that at least half of the maximum interannual variability in this flux is determined during March.Importantly and perhaps counterintuitively, our results show that years with less intense and discontinuous blooms are likely to have higher export flux than years with intense uninterrupted blooms.

Fig. 2 .
Fig. 2. Comparison between the mixed layer depth estimated from Argo floats data (circles) and that estimated by the model (asterisks) for the years 2005 and 2008.In the above panels the positions of the floats are displayed.

Fig. 3 .
Fig. 3. MODIS and model chlorophyll concentration averaged over the area of the bloom.Only the period from 1 December to 31 May is shown for each year.As a consequence, the two series are not continuous and the gray vertical bars represent the interruptions (from 31 May to 1 December).

Fig. 5 .
Fig. 5. Scatter plots of modelled heat fluxes vs. mixed layer depth (upper left), bloom onset date (upper right), primary production (lower left) and export flux (lower right).Each point represents the December through May average for an individual year(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010).The bloom onset date is calculated as the first day that meets the chosen criteria for the bloom onset within the week identified as the first week of the bloom (see Sect. 3, "Data treatment").
the authors observed deep convection and a complete homogenization of the water column down to the bottom (>2000 m) over an area of 50-100 km in diameter.The model simulates deep homogenisation of the water column also in the Ligurian subbasin, on the northeastern side of the bloom area.Recently, Smith et al. (2008) showed evidence of deep convection in the Western Ligurian subbasin during winter 2006 and in the Catalan Sea during winter 2005, both locations outside of the MEDOC area.Moreover, Hong et al. (2007) reported results from a numerical simulation showing how deep convection can extend into the Ligurian subbasin.

Fig. 6 .
Fig. 6.Area averaged mixed layer depth (black continuous, left axis), vertically integrated phytoplankton nitrogen biomass (0-75 m, red continuous) and primary production (circles) and organic nitrogen export flux (asterisks) at 75 m depth.The period from December to May is shown for 2005 (above) and 2008 (below).The vertical bars represent the day of the bloom onset (thin) and the day of the maximum mixed layer depth.

Fig. 7 .
Fig. 7. Model daily neat heat flux and mean wind speed averaged over the area of the bloom for the years 2005 (above) and 2008 (below).The vertical lines represent the day of the maximum MLD (thick line) and the spring bloom onset date (thin line).

Fig. 8 .Fig. 9 .
Fig.8.Cumulative EF for each year (thin lines), their average (dots) and the maximum range of variability (thick line -i.e., the highest year minus the lowest at each day).Circles represent the spring bloom onset date for each year.

Fig. 10 .
Fig. 10.Mean daily March EF (on the abscissa) and frequency of days during March with wind above the daily March 10-yr average (on the ordinate).
Maximum range of bloom start date from MODIS data (left) and model estimate (right).