Marine productivity response to Heinrich events: a model-data comparison

. Marine sediments records suggest large changes in marine productivity during glacial periods, with abrupt variations especially during the Heinrich events. Here, we study the response of marine biogeochemistry to such an event by using a biogeochemical model of the global ocean (PISCES) coupled to an ocean-atmosphere general circulation model (IPSL-CM4). We conduct a 400-yr-long transient simulation under glacial climate conditions with a freshwater forcing of 0.1 Sv applied to the North Atlantic to mimic a Hein-rich event, alongside a glacial control simulation. To evaluate our numerical results, we have compiled the available marine productivity records covering Heinrich events. We ﬁnd that simulated primary productivity and organic carbon export decrease globally (by 16 % for both) during a Hein-rich event, albeit with large regional variations. In our experiments, the North Atlantic displays a signiﬁcant decrease, whereas the Southern Ocean shows an increase, in agreement with paleo-productivity reconstructions. In the Equatorial Pa-ciﬁc, the model simulates an increase in organic matter export production but decreased biogenic silica export. This antagonistic behaviour results from changes in relative uptake of carbon and silicic acid by diatoms. Reasonable agreement between model and data for the large-scale response to Hein-rich events gives conﬁdence in models used to predict future centennial changes in marine production. In addition, our model allows us to investigate the mechanisms behind the observed changes in the response to Heinrich events.


Introduction
Marine primary productivity (PP) is a key component of climate-active biogeochemical cycles such as the carbon cycle.It also sustains upper trophic levels and marine resources (Pauly and Christensen, 1995) and is the first level of marine food web impacted by climate change.The response of PP to future climate change is largely uncertain (e.g.Taucher and Oschlies, 2011) and natural variability hampers the detection of unequivocal trends from the 15-yr ocean colour satellite record (Henson et al., 2010).Coupled climate -marine biogeochemical models are used to simulate the evolution of marine PP over the historical period and under future scenarios.Four models (Steinacher et al., 2010) show a global decrease in PP and in the export production of organic carbon at the base of the euphotic layer (EXP) to deeper waters of between 2 and 20 % by 2100, relative to preindustrial conditions, notwithstanding regional variability in the response.A reduced input of nutrients to the euphotic zone from subsurface waters due to increased stratification decreases PP in the North Atlantic and tropical regions, whereas lower light limitation increases Southern Ocean PP (Bopp et al., 2001;Steinacher et al., 2010).Another model (Schmittner et al., 2008) shows an increase in PP and a decrease in EXP due to enhanced remineralisation in response to increasing temperature.Nevertheless, evaluation of these models on such decadal-to-centennial time scales is still difficult due to sparse data covering these time scales (Schneider et al., 2008) and relatively moderate climatic change.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Turning to the geologic past, the last 70 000 yr of paleorecords may permit the evaluation of such climate-marine biogeochemical models on centennial time scales.These records show large and rapid variations linked to abrupt events such as Heinrich events (HEs) (Heinrich, 1988).During these events, massive iceberg discharges occur in the North Atlantic Ocean (Broecker et al., 1992), affecting the global ocean circulation through a collapse of the Atlantic meridional overturning circulation (AMOC) (Mc-Manus et al., 2004).At the same time, there is a cooling over the North Atlantic Ocean both seen in data (Voelker, 2002;Bard et al., 2000;Cacho et al., 1999) and reproduced in climate model experiments (Kageyama et al., 2010;Swingedouw et al., 2009;Ganopolski and Rahmstorf, 2001).Alongside the climate and ocean circulation changes, there is ample evidence of abrupt variations in biogeochemical variables.From ice cores, we know that atmospheric N 2 O concentrations decreased during HEs (Flückiger et al., 2004;Stocker and Schilt, 2008).Schmittner and Galbraith (2008) explain this decrease by a global reduction in primary productivity, especially in the Indo-Pacific region, which induces increased subsurface oxygen concentrations and then reduced N 2 O production.We can also find biogeochemical variations such as EXP variations in marine sediments (e.g.Anderson et al., 2009;Nave et al., 2007).While marine paleo-records give a picture of regional processes during HEs, such as a decrease in North Atlantic EXP and an increase in the Southern Ocean, they cannot give an integrated view of the global response of marine biogeochemistry to such events due to their sparse distribution.
Biogeochemical models can provide information concerning centennial-scale marine carbon cycle dynamics.Until now, most "paleo" setting studies aimed at understanding glacial-interglacial variations of atmospheric pCO 2 observed in ice core data (Monnin et al., 2001).To explain the lower pCO 2 during glacial times compared to interglacial times, different mechanisms have been investigated, such as Southern Hemisphere westerly-wind modification (Toggweiler et al., 2006), sinking of brines (Bouttes et al., 2010), marine biology enhancement through iron fertilisation (Bopp et al., 2003) or larger nutrient availability (Matsumoto et al., 2002).Alongside past CO 2 concentrations, recent models have also been used to understand changes in marine productivity between glacial and interglacial times (Kageyama et al., 2012;Oka et al., 2011;Tagliabue et al., 2009;Bopp et al., 2003).In agreement with the marine productivity compilation of Kohfeld et al. (2005), they all find a dipole effect in the Southern Ocean with enhanced EXP in the middle latitudes due to increased iron deposition and decreased EXP in the high latitudes due to increased light limitation following enhanced glacial sea-ice.The response in the other regions of the ocean are model dependent Nonetheless, on submillennial time scales, three model studies using Earth system models of intermediate complexity (EMICs) (Schmittner, 2005;Menviel et al., 2008;Schmittner and Galbraith, 2008) and two model studies using atmosphere-ocean general circulation models (AOGCMs) (Obata, 2007;Bozbiyik et al., 2011) have investigated the response of marine biogeochemistry to Heinrich-like events.Schmittner (2005), Obata (2007) and Bozbiyik et al. (2011) studies have been performed under preindustrial background climate, but as seen in Menviel et al. (2008) the structure of the marine productivity changes obtained under preindustrial and LGM conditions is fairly similar, so they can be compared to marine records of HEs.Bozbiyik et al. (2011) do not show changes in EXP, but the other studies do.They all simulate a global decrease of marine productivity and a common response in certain regions (North Atlantic Ocean, Benguela coast, Mauritanian coast) matching the data, whereas in other regions they do not agree with the marine response seen in data (Eastern Equatorial Pacific, Southern Ocean, see Table 1).These regional differences between model results and marine sediment records suggest that physical or biogeochemical processes might be missing or underestimated in these models.Such model-data mismatches need to be investigated further to understand the main mechanisms controlling key ocean regions.
Moreover, in the future, the Greenland ice sheet may melt, releasing an amount of freshwater that might resemble a HE, albeit smaller and released more to the North.The impact of such a release on the AMOC and marine biogeochemistry implies different opposing mechanisms with large uncertaintities (Swingedouw et al., 2007).Validating climate models with marine biogeochemistry data from the past is clearly a key to properly evaluate the response of marine system to future freshwater input.Importantly, it is still rare that identical marine biogeochemical models are employed and tested in "paleo" settings, as well as being used for predictions of our future climate.
In this study, we investigate the global and regional response of marine biogeochemistry to a Heinrich-like event with the state of the art biogeochemical model PISCES, which has also been used for future climate projections (Steinacher et al., 2010).In a first section we present the data compilation gathered, the experimental design and the comparison method performed between model and data.Then a second section details the results of the simulations both globally and regionally.A third section both discusses the discrepancies between our results and previous model studies and give some insights on potential effect of future Greenland ice sheet melting on marine biology.

Data compilation
In order to evaluate our numerical results, we have compiled existing marine sediment cores documenting millennialscale paleoproductivity changes during the past 70 000 yr.
Table 1.Data compilation of paleoproductivity changes in response to Heinrich events (HEs).For each record and each event, we summarise the trend of the proxies for paleoproductivity by a certain number of "−" or "+" signs.Each "−" (resp."+") corresponds to a sensor indicating a significant decrease (resp.increase) in EXP.We discriminate between significant and not significant trend in sensors by applying a simple test on each time series.We compare the trend of the first 500 yr following a dated HE to the mean and variance of the 3000 yr preceding this dated HE (resp.2000 yr for the Younger Dryas (YD)).If the value of the first 500 yr is higher (resp.lower) than the addition of the mean and the variance (resp.the difference between the mean and the variance) of the last 3000 yr, we consider it as a significant increase (resp.decrease).When there is not such a significant trend, we display an "x".Some records, while having a sufficient time resolution (less than 500 yr between two measurements), had no dated HEs or YD.Would we have discarded them, only 52 out of the 74 studies would have remained.We decided to keep them and to date the HEs and YD ourselves on the time series, taking as an onset for these events −13 kyr for YD, −18 kyr for HE1, −26.5 kyr for HE2, −32.7 kyr for HE3, −40.2 kyr for HE4, −50 kyr for HE5 and −63.2 kyr for HE6 (following Sanchez Goñi and Harrison (2010)).As there are time dating error bars issues on paleo-data, we consider these results as exploratory hypotheses, in the absence of a real timing for these events.To differentiate them from the well dated results we added a "?" in front of the concerned results.For each record, we summarise the results of all HEs and YD in one "value": when we add all the "-" and "+", if we obtain one "-" (resp."+"), we write one "-" (resp."+"), if we obtain more than two "-" (resp."+"), we write "--" (resp."++") and finally if we have an equal number of "-" and "+" (or only "x"), we write an "x" because we consider that the records does not give a significant trend.Following Kohfeld et al. (2005) for the choice in paleoproductivity proxies, we consider either products of direct biogenic origin such as organic carbon, biogenic opal and alkenones, or indicators of biological activity such as the ratio Ba/Al or the δ 13 C in foraminifers.Like all the species measured in marine sediments, we need to keep in mind they are subject to degradation in the water column and in the sediments, so the initial process they relate to might be deteriorated.

Core
Our compilation (Table 1) encompasses the last six HEs and the Younger Dryas (YD).We selected only high resolution cores (with less than 500 yr between two measurements) in order to capture the hypothetical change in biogeochemical variables induced by a HE and get a reasonable comparison with simulations at the centennial time scale.
Because paleoproductivity is reconstructed from different proxies that are not quantitatively comparable, our method is deliberately qualitative.We associate with each record a sign for productivity changes during HEs and a degree of confidence, which is deduced from a combination of the trend of the proxies, the number of proxies and the number of documented events.The details of our method are explained within the legend of Table 1.We end up with 74 data points that are gathered in Fig. 2 and will be discussed later.

The biogeochemical model PISCES
The Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) marine biogeochemical model includes a simple representation of marine ecosystem and of main biogeochemical cycles and is composed of 24 pools in total.Among them, there are two size classes of phytoplankton (nanophytoplankton and diatoms), two size classes of zooplankton, five pools of nutrients (phosphate, ammonium, nitrate, silicic acid, iron), small and large particulate organic carbon, one of dissolved organic carbon and another one of dissolved inorganic carbon.These pools interact with each other following the main biogeochemical processes such as photosynthesis, respiration, grazing, particle aggregation, particle sinking, remineralisation and sedimentation (for more details on the model, see Aumont and Bopp, 2006).Phytoplankton growth in the model is limited by five nutrients (nitrate, phosphate, ammonium, silicic acid, iron) and by light availability (a function of photosynthetically active radiation reaching the ocean surface, the optical properties of the water and the mixed-layer depth).The ratios for C/N/P are kept constant following the Redfield ratios.However, the ratio of silica and iron to carbon in phytoplankton biomass varies in function of nutrient availability and light.For instance, following culture experiments (Hutchins and Bruland, 1998;Takeda, 1998), the ratio of silica to carbon in diatom cells is modulated by the degree of iron limitation making diatoms more silicified under more severe iron limitation.
The model has already been evaluated under glacial conditions (Tagliabue et al., 2009;Bopp et al., 2003)    ) averaged for the simulated years 350-399 (filled field), alongside paleoproductivity changes during Heinrich events compared to glacial mean state reconstructed from our compilation (points).Dark and light blue points represent significantly lower (SL) and slightly significantly lower (SSL) export production respectively (equivalent to "--" and "-" or "?-" in Table 1).Dark and light red points represent significantly higher (SH) and slightly significantly higher (SSH) export production respectively (equivalent to "++" and "+" or "?+").
The grey points represent no significant change (NSC) (equivalent to "x" or "?x").
reproduces roughly the paleoproductivity reconstruction of Kohfeld et al. (2005).One of the more consistent patterns is a dipole effect in the Southern Ocean with enhanced EXP in the middle latitudes due to increased iron deposition and decreased EXP in the high latitudes due to increased light limitation following enhanced glacial sea ice.

Experimental design
PISCES is forced offline by the atmosphere-ocean general circulation model IPSL-CM4 (Marti et al., 2010), which includes the ocean dynamical model NEMO with 2 • × 2 • -0.5 • horizontal resolution and 31 vertical levels, 10 being located in the first 100 m.
Two 400-yr-experiments have been performed under Last Glacial Maximum (LGM) conditions, with the orbital parameters, greenhouse gas concentrations and ice sheets from 21 000 yr before present (see Kageyama et al., 2009 for a detailed presentation of the climate setup).The biogeochemical simulations based on these experiments use a constant atmospheric CO 2 concentration fixed at LGM level (190 ppm) and constant LGM dust deposition distribution (Mahowald et al., 2006), which is important for aeolian iron interaction with marine biology.The first experiment is an equilibrated glacial run (GLA) used as a reference run.The second experiment is a hosing experiment (FWF), starting from year 100 of the reference run.In this experiment, the and in the Arctic Ocean, which mimics the icebergs melting during an HE.In FWF simulation, the AMOC collapses in around 250 yr.This is a relatively long time response compared to the simulation time but a relatively short time response compared to the resolution of most of the marine records.Indeed these records do not allow to distinguish if the changes in the AMOC happened within 10 or 400 yr.One limitation of this experimental set-up is that it uses a full LGM reference state rather than a reference state closer to the MIS3 conditions, i.e. with smaller ice-sheets, less depleted atmospheric greenhouse gases and orbital parameters favouring more insolation in the summer hemisphere (e.g.Van Meerbeeck et al., 2009).We chose an LGM reference state because it was closer to MIS3 conditions than a preindustrial run and because we obtained an abrupt collapse of the AMOC for a rather small amount of freshwater hosing, which is relevant for the study of mechanisms occurring during Heinrich events.This is partly due to the AMOC of the reference state being rather strong and hence closer to an interstadial than to the full LGM.Besides, Heinrich events 1 and 2 did happen actually relatively close to the LGM so these conditions are also relevant for this reason.
For this study, we focus on the biogeochemical results of these simulations.Figure 1 gives an overview of the main dynamical fields (sea surface temperature, mixed-layer depth, wind stress and Ekman pumping anomalies, sea-ice extent) useful in the interpretation of the biogeochemical results.Further ocean atmosphere dynamics details of these simulations can be found in Kageyama et al. (2009) where our simulations GLA and FWF correspond respectively to the simulations LGMb and LGMc.
We first study EXP as it is clearly more relevant to be compared with sediment core observations.But, as Taucher and Oschlies (2011) have pointed out for projection simulations, there might be cases where the response of PP and EXP to climate can be decoupled, so we will examine this potential decoupling once the fidelity of the model performance has been assessed.
In the following, we define a typical Heinrich-like event in the model through the difference between the FWF and GLA simulations averaged over the last 50 yr of the 400-yr simulations.We compare such a signature with our data compilation, where the mean response of the 7 events is considered to represent a typical HE in the observations.Nonetheless, we need to keep in mind that those reconstructions span time periods that were very different in solar irradiance due to the changes in the Earth's orbit, which may have affected the response.The comparison between 50-yr averaged simulations and 500-yr resolution data is not ideal, but because of calculation time constraints, we could only run two simulations of 400 yr each.Schmittner (2005) and Schmittner and Galbraith (2008) highlight the time dependence of the response in EXP to HEs.For example, full reduction in PP in the Indian and Pacific oceans is only expressed more than 1000 yr after the AMOC shutdown in their simulations.It would be interesting to consider the time component in longer simulations than what was possible in this study.

Statistical match between model and data
In total, we found 74 marine records studies capturing at least one HE (or YD).In short, 35 record a decrease in EXP during this period, while 21 show an increase and 18 do not provide a significant trend (Table 1).If we consider the 56 records that display a significant trend, the model outputs match 35 of these cores (Fig. 2).Among the 21 modeldata mismatches, 5 of the data points (ODP Site 1144, Core 17950-2, GeoB3359-3, MD97-2120 and NBP9802-6PC, see Table 1 for precise location) are located on a model front area.PISCES does represent the main features of nutrients distribution, but some fronts can be shifted, due to the biases of the model, so this could explain part of the mismatches.Six of the data points (MV99 GC31/PC10, MD02-2508, MD02-2519, MD02-2524, MD02-2529 and TN057-13-4PC) are located in an area where the model does not simulate a significant change in EXP (between −1 and 1 gC m −2 yr −1 ).Three data points (PL07-57PC, EW0408-11JC and M23415) are displaying an increase in EXP, while other data points (respectively M35003-4 and PL07-39PC for PL07-57PC; EW0408-85JC and EW0408-66JC for EW0408-11JC; SU9039 and BOFS-5K for M23415) within the same area display a decrease in EXP just as the model does.These last data points raise the problem of comparing local data records to regionally averaged model cells.Even if we consider that the signal captured by marine cores has been well preserved within both the water column and the sediments, each data record is nonetheless the addition of a regional signal and a local signal.The model simulates only the regional part of the signal, and we would need several data points for one cell of the model to be entirely confident with the statistical comparison.Five points (MD04-2805 CQ, GeoB5546-2, GeoB7925-2, GeoB9526-5 and GeoB9527-5) are located on the Mauritanian coast and are clearly in contradiction with the model results.We will discuss in more details the results for this region in the following section.Finally, the two remaining non-matching data points (GeoB3912-1 and ODP Site 1233) are both located on the coast, so there might have been continental influences on the results that are not included in our simulations.
In conclusion, if we exclude the cores located on the model front areas and the ones that do not agree regionally with other nearby records, we match the sign of the response of EXP in 73 % of the locations (35 over 48).Hence, our model  seems to be able to simulate the main response of EXP to HEs correctly, or at least its first order mechanisms.

Regional analysis
In order to use our model-data comparison to better understand the response of EXP in key regions, we define 7 box regions (NATL, BEN, EEP, IND, MAU, NZL and SEP, acronyms defined in Table 2 and regions visible in Fig. 2).We first focus on regions where proxy data is available: either model and data are in agreement as in NATL, BEN, EEP and IND, or there is a clear contradiction as in MAU.In addition, we also specifically focus on the Southern Ocean, because of its meridionally diverse response of EXP to HEs.

North Atlantic
The North Atlantic (NATL) is the region where most of the marine records are located (19 in total, 13 with significant trends, see Table 1) and our modelled Heinrich-like event matches the sign of all the significant records except for one (M23415, Weinelt et al., 2003) which records a slightly significant increase in EXP while two other nearby cores (SU9039 and BOFS-5K) record a significant decrease in EXP.In the model, a shallowing winter mixed-layer depth (by more than 150 m in February, Fig. 4a) reduces the nutrient flux to the upper ocean (not shown).This was already shown in Schmittner (2005), who also gets a similar percentage decrease in EXP.Moreover, an increasing sea-ice cover (Fig. 4b) reduces light availability.Both these processes cause a decrease in EXP of 44 % (Table 2).Gil et al. (2009) find an increase in EXP at core OCE326GGC6, located south of the present day Gulf Stream.Their data did not pass our test for significant data because the duration of their record before HE1 was less than 3000 yr, limiting the comparison to a glacial state.Nonethe- less, they explain this increase in EXP by an iceberg migration to the subtropics inducing an isolated environment involving turbulent mixing, upwelled water and nutrient-rich meltwater supporting productivity.This hypothesis has not been tested in our model set-up and could not be confirmed or dismissed because (1) we do not have enough horizontal resolution to account for meso-scale processes and (2) we do not take into account nutrient input accompanying freshwater discharge.
Overall, our model matches the reduction in EXP during HEs in the North Atlantic in most cases and we suggest that the response results from greater limitation of PP and thus also of EXP by both nutrients and light.

Southern Ocean
In the Southern Ocean, there are 4 records indicating an increase of EXP (MD97-2120, TN057-13-14PC, E27-23, NBP9802-6PC).Our model shows significant zonal variability in this area (see Fig. 2).Nonetheless, these records are located in or close to simulated increased EXP regions.We focus here on two regions of interest that illustrate the major trends: an area in the southeast of New Zealand (NZL) and an area in the southeastern Pacific (SEP).
In the NZL area, south of the polar front, our model simulates a 6.4 % increase in EXP (Table 2) in agreement with two different sediment cores from this same area (E27-23, NBP9802-6PC).The model suggests that a deepening of the Austral winter mixed layer (Fig. 5d), combined with a strengthening of the upwelling (see positive Ekman pumping anomalies in Fig. 5c), act together to increase the nutrient flux to the euphotic zone and stimulate PP and EXP.In particular, we note an increase in silicic acid concentrations in surface waters in FWF compared to GLA (see Fig.  b) as suggested by Anderson et al. (2009).Silica export positive anomalies in the NZL region (Fig 3) confirm this last finding.In the model, enhanced simulated wind stress intensifying the upwelling and retreat of Austral winter sea ice poleward (Fig. 5d) both reduce the insulation effect and increase the momentum flux between the atmosphere and the ocean, which a deepening of the mixed layer.The retreat of sea ice is itself linked to increased SST (see-saw effect, Fig. 1a).
In the SEP area, EXP decrease by 17 % in the model (Table 2), in contrast to the NZL region.Our model simulates increasing Austral winter sea ice that isolates the surface waters from the atmosphere in this area.This "barrier effect" of sea ice means that winds cannot deepen the winter mixed layer enough to introduce nutrients in the upper ocean, which would be available for consumption by phytoplankton in spring.Moreover, greater sea ice also reduces the springtime light availability.Both these processes (reduced winter mixing and less springtime irradiance) induce a decrease in EXP in this area.As far as we know, no paleo-productivity proxy with millennial-scale resolution has been published in this area, so we cannot compare our results with any reconstructions.Nonetheless, we found one sediment core (E11-2 in Mashiotta et al., 1999) displaying SST reconstructions for the last deglaciation and HE1.This core shows that SST increases from the LGM to the Holocene in this region, so in particular during HE1.From our model results, we would have expected a decrease of SST during HE1 in this area, to explain the simulated increased sea-ice extent.Some global climate models also find this SST anomaly in the SEP, but others do not (Kageyama et al., 2010;Merkel et al., 2010;Otto-Bliesner and Brady, 2010;Kageyama et al., 2009), thus this pattern could be model dependent.We need further high resolution data to discriminate between model results and better understand the main processes involved in this area.
Overall, our model matches the increased EXP noted in the sediment cores in the Southern Ocean, but we highlight a large degree of spatial variability in the response of EXP to different patterns of winter mixing and sea ice.

East Equatorial Pacific
In the East Equatorial Pacific (EEP), 5 cores studies (3 with a significant trend) are available at the same location (they appear as one red point in Fig. 2) and show an increased EXP.In this area, the model indicates a 2 % increase in EXP (Table 2), which is the result of greater equatorial upwelling due to stronger trade winds (see Fig. 1c).This mechanism increases the nutrient availability in this area and agree with concomitant decreased reconstructed SST and increased organic carbon content in sediment (Kienast et al., 2006).Although the model simulates an increase in EXP (Fig. 6d), it shows a simultaneous decrease in the export of silicate (Fig. 6c), a biomarker for diatoms, which might appear counter-intuitive.However, laboratory data has shown that when diatoms are iron limited, they adapt to this new environment consuming more silicate relative to C and N (Hutchins and Bruland, 1998;Takeda, 1998) thereby allowing the export of Si and C to become decoupled; this process is included in our model.In the context of HEs, the exact opposite happens, with diatoms experiencing greater iron availability, which decreases their relative silicic acid uptake and decreases the Si/C ratio for diatoms (Fig. 6b).This process drives a decrease in the export of silicate (see also Fig. 3, EEP region), even if EXP increases.This process has already been pointed out by Pichevin et al. (2009) at ODP Site 1240 (see Table 1) for glacial/interglacial time-scales and our simulations reveal that it seems to be a critical process on submillennial time-scales as well (see Fig. 6).The EEP is in fact a high-nutrient, low-chlorophyll (HLNC) region, which is iron-limited, so variation in the input of iron can induce high variations in EXP.As for the origin of iron, Leduc et al. (2007)   that conveys more airborne iron.In our model, the airborne iron flux is kept constant at glacial levels between the two simulations, so we cannot capture this effect.Nevertheless, the model simulates greater iron supply to EEP surface waters due to enhanced vertical supply.The subsurface ocean is clearly to be considered as a potential source for iron during HEs in this region.Even if the simulated increase (2 %) in EXP is fairly small, it is encouraging that our model reproduces the trend of EXP recorded in proxies, as well as capturing the decoupling between the export of carbon and silica noted in the geologic record in EEP.

Coastal regions
Coastal regions are not the best areas to test our model results as the model cannot capture specific coastal processes because of its coarse resolution.These areas are nevertheless regions where most data are available because of an important sedimentation rate, which allows a high-resolution analysis, so we endeavour to examine them.Our assumption is that the main signature found in coastal area is related to large-scale changes.Four proxy-based studies are available on the Mauritanian coast (MAU); they all find increased EXP in response to HEs linked to an enhanced upwelling (Penaud et al., 2010;Holzwarth et al., 2010;Zarriess and Mackensen, 2010;Romero et al., 2008).However, in contrast to these records, our runs simulate a 58 % decrease in EXP (Table 2).We do observe an enhanced upwelling in our simulations (see positive Ekman pumping anomalies in Fig. 1c) but it is completely offset by the overall thinning of the mixed layer (Fig. 1b) induced by the freshwater forcing in North Atlantic.It is plausible that our idealised freshwater forcing might be too strong compared to real HEs or that the zone of freshwater input does not exactly correspond to the region where icebergs melt.
On the Benguela coast (BEN), EXP decreases both in data (GeoB1706-2, GeoB1711-4, GeoB3606-1) and in the model (by 66 %, Table 2).In the mean glacial state, there is an important upwelling in this area that decreases significantly during our HE experiment.Therefore, most of nutrient enriched sub-surface waters do not reach the euphotic layer, which reduces EXP by increasing nutrient limitation.
In the Indian Ocean (IND), more precisely in the Arabian Sea, the model simulates decreased EXP (by 47 %, Table 2), in good agreement with high resolution data (Ivanochko et al., 2005;Higginson et al., 2004;Altabet et al., 2002;Schulte and Müller, 2001).In this region, EXP is primarily controlled by upwellings, themselves induced by monsoon Westerlies (Bassinot et al., 2011).Kageyama et al. (2009) found a weaker monsoon in the Heinrich-like event simulation we use.This weaker monsoon thus induces weaker upwellings and then decreased EXP.Overall, despite the potential problems in comparing global model results to coastal sediment cores, our model succeeds in reproducing the observed trends in the BEN and IND regions.We suppose that the mismatch in the MAU region could be due to the highly idealised way in which we simulate HEs.

Global analysis
When globally integrated, EXP decreases by 16 %, or 1.5 Pg (Table 2 and Fig. 7b) by the end of the simulations, in response to a HE.EXP strongly depends on nutrients concentrations in the illuminated waters, and the nutrient availability itself is highly dependent on nutrient supply through ocean ventilation and mixing.The constant freshwater flux that we use to approximate a HE induces a decrease of the thermohaline circulation (shown by a reduction in the strength of the AMOC by 87 % (or 13 Sv), see Fig. 7a) associated to a strong ocean stratification in the North Atlantic.This leads to a reduced upwelling of deep water (as already shown in Schmittner, 2005) and a decrease in the ventilation of the subsurface ocean, which induce a decrease of nutrient supply (cf. the de-crease of almost 13 % of global nitrate concentration at the surface in Fig. 7c).This explains the global decrease in EXP.We note that whereas the AMOC stabilises at around 2 Sv by the end of our simulation, EXP continues to decrease linearly over the entire 400 yr period.
As explained in the experimental design section, we have chosen to compare our modeled EXP to available marine productivity data, making the hypothesis that PP and EXP are varying in the same direction.However, Taucher and Oschlies ( 2011) recently showed that it is not always the case in response to climate variability.When the temperature increases, metabolic effects cause an increase in both PP and remineralisation of organic matter by bacteria.Greater remineralisation can reduce EXP, but may also yield positive feedback on PP via the subsequent increase in available renewed nutrients due to greater heterotrophy.In order to investigate if PP and EXP respond similarly in our experiment, we plot the ratio of the comparative change in these two quantities in Fig. 8.For the areas we focused on here (Fig. 2) PP and EXP vary in the same direction, thus our modeled EXP is a correct "proxy" of PP in these areas.Alternatively, there are some regions where the model simulates opposite responses for PP and EXP.Most of these areas are located in the boundary between an increased and a decreased EXP so they may be due to horizontal advection effects.Other PP-EXP decoupled areas, like the ones south of New Zealand correspond to regions showing lower EXP, higher PP and warmer temperatures (Fig. 8).Taucher and Oschlies (2011)'s hypothesis can thus partly explain the differences observed between PP and EXP (the parameterization of PP and remineralisation is indeed dependent on temperature in PISCES).

Comparison with other model studies
Other model studies have examined the response of EXP to freshwater forcing with models of intermediate complexity (Menviel et al., 2008;Schmittner, 2005) or atmosphereocean general circulation models (Obata, 2007).The freshen-Fig.8. Comparative change in PP (primary productivity) and EXP (export production) (with X = X(FWF)−X(GLA) X(GLA) , X averaged on simulated years 350-399) (filled field) and sea surface temperature anomaly (isolines each 0.5 degrees Celsius) for the simulated years 350-399.
ing scenarios applied were different from those employed in this study, so we cannot compare the results quantitatively, as explained by Bouttes et al. (2012) who discuss how different hosing scenarios modulate Heinrich events and their impact on the ocean carbon cycle.We can however qualitatively examine the patterns in the EXP response between models.
In general, the reduced EXP in NATL and BEN are consistent between the models suggesting that these are robust responses to freshwater fluxes input in the North Atlantic.However, there are regions where the models have a different response to hosing.We will focus on EEP, NZL, MAU and IND.These are regions where EXP is mainly controlled by upwelling.As coarse-resolution climate models usually have difficulties to simulate upwellings correctly, we first check that upwellings are well represented in the modern version of each model before drawing any conclusion on glacial centennial-scale changes in these regions.The IPSL-CM4 model represents all these upwellings (Bassinot et al., 2011;Steinacher et al., 2010;Lenton et al., 2009;Schneider et al., 2008) though it tends to underestimate their intensity both in terms of upwelled water flux and surface productivity.The LOVECLIM model represents the EEP and NZL upwellings but not the MAU and IND ones (see Menviel et al., 2008, Figs. 5 and 6).The UVic model represents the EEP, NZL and IND upwellings, but not the MAU one (see Schmittner, 2005, Fig. 2).The MRI model represents all the upwellings (see Obata and Kitamura, 2003, Fig. 3a).We need also to point out that while Menviel et al. (2008)'s, Obata (2007)'s and our simulations are performed with interactive winds, Schmittner (2005)'s simulations use prescribed modern winds, so his simulations cannot capture wind-driven changes in upwellings.In EEP and most of the Southern Ocean, we find an increased EXP which is in contrast to the three previous model studies.Note however that when forced by changes in insolation and CO 2 for HE1, LOVECLIM also simulates an increase in EXP for the Southern Ocean (Menviel et al., 2011).For both regions, we attribute this difference to the winddriven increased upwelling.We explain the difference with Schmittner (2005)'s study by his non representation of winds changes.For Menviel et al. (2008)'s andObata (2007)'s studies that do represent these upwelling areas in modern times and computes wind changes, we need to find other explanations.We hypothesise that the discrepancies between Menviel et al. (2008)'s model and ours are probably due to two factors.On one hand, our model has an increased atmospheric resolution, so we hypothesise it captures better the wind changes in upwelling areas.On the other hand, our model has a parameterization of Si/C as a function of temperature and iron availability that is not implemented in LOVE-CLIM and we have shown in the results section that this parameterization was key to simulate the EEP region changes.For the discrepancies between Obata (2007)'s model and ours, the explanation is different.The MRI-CGCM2 model has a comparable resolution to the IPSL-CM4 one, so the discrepancies cannot come from a problem of low resolution winds representation.In Obata ( 2007)'s Heinrich-like simulation, the sea-ice cover increases in the Southern Ocean, which should induce a shallower mixed layer and also a decreased upwelling.This can explain why the MRI-CGCM2 model simulates a decreased EXP in this area.For the EEP region, the MRI-CGCM2 model does not have either the parameterization of Si/C implemented in the PISCES model, so the model is not able to simulate an increase in EXP.Nonetheless, the EEP region is fairly narrow in our simulations and is bounded to the north and south by negative EXP anomalies, so we would need further marine records to document this area and draw robust conclusions.

Clim
In IND, our study finds a decrease in EXP due to a weaker upwelling as do Schmittner (2005) and Obata (2007), whereas Menviel et al. (2008) simulate an increase in EXP.As Schmittner (2005)'s simulations use prescribed preindustrial winds, the decrease in upwelling seen in his simulations is thermohaline-driven.As explained above, the modern simulations of LOVECLIM do not represent the upwelling in this area, so it cannot obviously simulate changes of upwelling regimes.
In MAU, the decreased EXP we simulate is at odds with both data and previous model studies, which present an increased EXP due to an enhanced upwelling.We need to explain why our model, supposed to be able to capture changes in wind-driven upwellings, is not able to capture the signal of increased productivity seen in data.We make the hypothesis that it is due to the location of our hosing.The hosing is applied between 40 • N and 90 • N in our experiment, which is south enough to make the freshwater being advected through the subtropical gyre directly to the Maurita-nian coast.On the contrary, in the other studies, the hosing is applied more northward (between 50 • N and 65 • N for Menviel et al. (2008), between 50 • N and 70 • N for Obata (2007) and between 45 • N and 65 • N for Schmittner (2005)), which makes the freshwater flow northward through the North Atlantic Current.Hence, the surface waters are considerably mixed through the subpolar gyre before they arrive on the Mauritanian coast.We do have an increase of upward vertical velocities in the upper ocean of this region (see positive Ekman pumping in the Mauritanian coast in Fig. 1c), but this increased upwelling is completely balanced by the freswater lid that does not allow the nutrients to spread out from the subsurface waters.In addition, the modern simulations previously performed with our model (Steinacher et al., 2010) display a really weak upwelling in this area compared to data, so this underestimation of vertical velocities can explain why it cannot balance the lid effect.Nonetheless, the other studies simulate an increase in productivity in MAU.This is understandable for Obata (2007) because the MRI-CGCM2 model represents the Mauritanian upwelling and the simulations have interactive winds.As seen before, the two other studies do not represent properly the Mauritanian upwelling in preindustrial settings.The greater EXP obtained in Menviel et al. (2008) off the Mauritanian coast might be due to a greater upwelling and/or an increase in nutrient content in the source water.Indeed, the stratification generated in the North Atlantic by the freshwater input leads to a subsurface positive nutrient anomaly, which is advected to lower latitudes.In Schmittner (2005), as the surface winds are constant, the EXP increase in MAU is mainly due to a greater thermohaline driven upwelling and/or a greater nutrient content in the source water.To summarise, in MAU, according to data, there was an enhanced upwelling during HEs inducing increased EXP, but our model cannot reproduce it because of (1) a too weak intensity of the upwelling in modern times and (2) a freshwater forcing too close to the upwelling area.
In conclusion, the differences between model results come mostly from atmospheric resolution and biogeochemical parameterization differences but also from the different locations of freshwater forcing and the way the models advect them from the area of freshening.Overall, except for the Mauritanian coast, our study correctly simulates the response of EXP to HEs, and points out the importance of the Si/C ratio parameterization in marine biogeochemical models.

What do these results tell us for 21st century projections?
We do not yet have any certainty on the sign of the global evolution of marine PP with global warming (Taucher and Oschlies, 2011).Investigating the response of marine biogeochemical models to HEs could be of use in examining the predicted impact of climate change on marine biogeochemistry.Valdes (2011)  The problem with the Mauritanian upwelling system needs though to be fixed in the future for more reliable projections, as this upwelling system is very important for future food (e.g.fish) supply.Finally, we need to point out that while our model simulates a global decrease of PP by 16 % with HEs, we still do not have any global constraint to validate this result.New isotopic methods using for instance the triple isotopes (Landais et al., 2007) or the Dole effect (Landais et al., 2010) should be used to investigate this result.Simulating HEs could be a benchmark for coupled climate carbon cycle models to test their ability to simulate abrupt transient climate changes, and our compilation of paleoproductivity proxies could be compared to the results of the models.
Studying the response of EXP to HEs can also give insights of mechanisms that may affect EXP under global warming.The radiative effect of increased CO 2 and subsequent warming on marine biology has already been tested.Steinacher et al. (2010) performed a model inter-comparison between four coupled climate-carbon cycle models for future climate.Significant regional differences between the models in the response of EXP to climate change appear but there are shared patterns like a decrease of EXP in NATL and an increase in the Southern Ocean.Global warming is also accompanied by a melting of the Greenland ice sheet, which is not taken into account in most of the actual coupled models.Our study could represent an analogue to this future Greenland melting, which is implemented in the study of Steinacher et al. (2010) only for one model (IPSL-CM4) out of four models.Of course our study starts with a glacial climate background so this could induce differences in the intensity of the response of the system compared to the same freshwater forcing with an interglacial climate background, so we need to be careful when comparing our HEs simulations with global warming freshwater forcing projections.Nonetheless, we can point out some significant trends, like a Fig. 9. FWF-GLA time response (in years) of EXP (export production): defined as the time after which half of the final signal (averaged on the simulated years 350-399) in EXP has been reached.We only considered regions where the final anomaly of EXP was more than 1 gC m −2 yr −1 in absolute value: the white areas in the ocean do not have such an important anomaly in the end.decrease of EXP in the North Atlantic, a region that is already projected to undergo a decrease of EXP due to global warming.This result is in agreement with the fact that the IPSL model in Steinacher et al. (2010) study was already simulating a more important decrease in this region.Hence, actual projections may underestimate the decrease of EXP in this region.Swartz et al. (2010) have shown that actual increased fishing induces a higher percentage of required PP to sustain global fish populations, with a special increase in North Atlantic Ocean.Projections including a freshwater forcing may be of use to help constraining EXP response in the future, especially for the areas of actual intense fishing that could be strongly affected in the coming decades.

Conclusions
This study first aimed at evaluating the response of a marine biogeochemistry model (PISCES) to centennial-scale events in glacial times, using marine cores for comparison.Despite the fact that we used full glacial boundary conditions rather than more realistic MIS3 conditions, the model results regarding the response of marine biology to Heinrich events are most of the time qualitatively consistent with paleo-data, which is encouraging for its ability to simulate future climate impacts on primary productivity and especially abrupt centennial-scale changes.The data compilation for paleoproductivity we gathered and used to test our model results can be used as a tool to evaluate other coupled biogeochemicalclimate model responses to Heinrich events and their ability to simulate a centennial-scale climate change.Our work also highlights the importance of the Si/C ratio parameterization in models as a key mechanism to simulate certain regions ecosystem, here in the Eastern Equatorial Pacific.This study also points out the importance of multi-proxy analysis to interpret paleoproductivity in the sediments.
The second aim of this study was to use the model to more accurately understand the global and regional response of marine productivity to Heinrich events.We simulate a global decrease of primary productivity of 16 % following the freshwater forcing, with some regional differences.According to our data-model intercomparison, it is very likely that the North Atlantic Ocean, the southwestern coast of Africa and the Indian Ocean experienced a decrease in primary productivity, whereas the Southern Ocean and the Eastern Equatorial Pacific experienced an increase during Heinrich events.This study also gives us an insight of what could be the contribution of a melting of Greenland ice sheet in the coming century: an accentuated decrease of organic matter export in the North Atlantic Ocean.

Fig. 1 .
Fig. 1.Changes in some physical fields between the FWF and GLA simulations, averaged for the simulated years 350-399.(a) Sea surface temperatures (SST) anomalies (shaded area, in • C) and sea-ice extent (each contour represents the isoline of 10 % of sea-ice cover, respectively in green for the GLA run and in red for the FWF run), (b) Mixed-layer depth (MLD) anomalies (shaded area, in m) and sea-ice extent and (c) Ekman pumping (m yr −1 ) and wind stress (Pa) anomalies.The Ekman pumping is masked between 10 • S and 10 • N because it is not defined (Coriolis parameter f is close to 0).

Fig. 5 .
Fig. 5. Silicic acid concentration (in µmol l −1 ) for (a) GLA and (b) FWF runs (transect averaged longitudinally over NZL area (see longitudes details in Table 2; (c) Ekman pumping (shaded area, in m yr −1 ) and wind stress (vectors, in Pa) anomalies; (d) mixed layer depth (shaded area, in m) and sea-ice extent (each contour represents the isoline of 10 % of sea-ice cover, respectively in green for the GLA run and in red for the FWF run).All the fields plotted here are averaged on months August and September for the simulated years 350-399.

Fig. 6 .
Fig. 6.Average on the EEP area (simulated years 350-399) of (a) diatom concentration (in molC m −2 averaged on the first 200 m of the water column), (b) Si/C ratio (also averaged on the first 200 m of the water column), (c) silica export (ExpSi) at 100 m (in molSi m −2 ) and (d) export production (EXP) at 100 m (in molC m −2 ).

Table 2 .
Definition and characteristics (latitudes, longitudes, area)of the regions chosen for the data-model comparison.The column EXP(GLA) gives the amount of carbon exported to the deep waters for each region defined.This amount is an average over the last 50 yr of simulation GLA.The column EXP gives the difference in carbon exported to the deep waters between the simulations FWF and GLA.For example, "−0.27 (−44)" means that in the FWF simulation, the NATL region exports 0.27 PgC yr −1 less than in the GLA simulation, which corresponds to a 44 % decrease.
showed that during HEs, the southward shift of the Intertropical Convergence Zone (ITCZ) can induce drier air www.clim-past.net/8/1581/2012/Clim.Past, 8, 1581-1598, 2012 Beaufort et al. (1997)) changes like the actual global warming and that they need to be tested on past abrupt climate changes.This study shows that the model IPSL-CM4 including PISCES is able to represent the main features of EXP response to a HE, except for the Mauritanian region which needs to be further investigated.With the sediment data currently available, we are not able to test if the model answers quantitatively well.We need further more attempts of PP or EXP calibration from different biomarkers as it has been done bySalgueiro et al. (2010),Voelker et al. (2009)orBeaufort et al. (1997)to have a more direct comparison with our model outputs.Nonetheless, the model can simulate an EXP response to HEs qualitatively consistent with available data, and more importantly with a fast time response in certain regions, like in the Atlantic Ocean which responds strongly within a hundred years (see Fig.9).Even with an AMOC decreasing on slower timescales ( 250 yr), PISCES model forced by IPSL-CM4 output seems to be able to simulate transient climate changes on centennial scales.This is encouraging as it has been used for climate projections.