Modelling the 13C and 12C isotopes of inorganic and organic carbon in the Baltic Sea

Journal of Marine Systems 148 (2015) 122–130

Contents lists available at ScienceDirect

Journal of Marine Systems journal homepage: www.elsevier.com/locate/jmarsys

Modelling the 13C and 12C isotopes of inorganic and organic carbon in the Baltic Sea Erik Gustafsson a,⁎, Carl-Magnus Mörth a,b, Christoph Humborg a,c, Bo G. Gustafsson a a b c

Baltic Nest Institute, Baltic Sea Centre, Stockholm University, Stockholm 106 91, Sweden Department of Geological Sciences, Stockholm University, Stockholm 106 91, Sweden Department of Applied Environmental Science, Stockholm University, Stockholm 106 91, Sweden

a r t i c l e

i n f o

Article history: Received 28 November 2014 Received in revised form 19 February 2015 Accepted 23 February 2015 Available online 3 March 2015 Keywords: Baltic Sea Modelling Carbon isotopes

a b s t r a c t In this study, 12C and 13C contents of all carbon containing state variables (dissolved inorganic and organic carbon, detrital carbon, and the carbon content of autotrophs and heterotrophs) have for the first time been explicitly included in a coupled physical–biogeochemical Baltic Sea model. Different processes in the carbon cycling have distinct fractionation values, resulting in specific isotopic fingerprints. Thus, in addition to simulating concentrations of different tracers, our new model formulation improves the possibility to constrain the rates of processes such as CO2 assimilation, mineralization, and air–sea exchange. We demonstrate that phytoplankton production and respiration, and the related air–sea CO2 fluxes, are to a large degree controlling the isotopic composition of organic and inorganic carbon in the system. The isotopic composition is further, but to a lesser extent, influenced by river loads and deep water inflows as well as transformation of terrestrial organic carbon within the system. Changes in the isotopic composition over the 20th century have been dominated by two processes — the preferential release of 12C to the atmosphere in association with fossil fuel burning, and the eutrophication of the Baltic Sea related to increased nutrient loads under the second half of the century. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The isotopic composition of carbon in coastal areas is influenced by properties of external carbon sources (e.g. river loads and atmospheric deposition), isotope fractionation during processes such as primary production and air–sea CO2 exchange, as well as interannual changes and long-term trends in e.g. atmospheric CO2 concentration and isotopic value, and production/mineralization of organic carbon. Riverine dissolved and particulate organic carbon (DOC and POC) is typically well depleted in 13C (~− 28‰ to − 27‰ compared to the PDB standard (Deutsch et al., 2012)). Depending on properties of the catchment area, the isotopic composition of dissolved inorganic carbon (DIC) in rivers can vary within a comparatively large range, but is generally less depleted in 13 C than the terrestrial organic carbon. Primary producers are preferentially assimilating 12C, resulting in 13 C enrichment in surface water DIC during the productive season. Mineralization of organic carbon again decreases the δ13C of DIC. Changes in nutrient loads affect CO2 assimilation, POC deposition and burial and by that also the isotopic composition of inorganic and organic carbon in water column and sediments. Allochthonous and autochthonous organic carbon typically exhibit different isotopic compositions, resulting in distinct isotopic signals related to respiration of organic ⁎ Corresponding author. E-mail address: [email protected] (E. Gustafsson).

http://dx.doi.org/10.1016/j.jmarsys.2015.02.008 0924-7963/© 2015 Elsevier B.V. All rights reserved.

carbon from terrestrial and marine sources respectively (e.g. Alling et al., 2012). The δ13C of atmospheric CO2 is continuously decreasing as a result of CO2 evasion from fossil fuel burning (the Suess effect). Thus, the oceans are because of the rising atmospheric CO2 partial pressure (pCO2) increasingly absorbing CO2, and the atmospheric CO2 is gradually becoming more and more depleted in 13C. In this study, carbon isotopes are included in the BALTSEM model, enabling calculations of carbon isotope fractionations related to air– sea exchange, primary production, and mineralization of organic carbon from both marine and terrestrial sources. Our main purposes are to 1. Describe the model parameterizations. 2. Estimate the model sensitivity to isotopic composition of e.g. riverine DIC and DOC, the fractionation during primary production and air–sea CO2 fluxes, etc. 3. Examine the sensitivity of δ13C to long-term changes in DIC and organic carbon related to the Suess effect and eutrophication. The study begins with a model description in Section 2. Results are presented in Section 3, followed by a sensitivity study and discussion in Section 4, and finally concluding remarks in Section 5. 2. Material and methods 2.1. Model description BALTSEM is a coupled physical–biogeochemical model where the Baltic Sea is divided into 13 interconnected sub-basins (Figure S1).

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

Each basin is described as horizontally homogeneous but with a high vertical resolution and depth dependent area distribution. A hydrodynamical module (Gustafsson, 2000) simulates transports of state variables between and within basins. Nutrient and plankton dynamics are simulated by a biogeochemical module that closely follows Savchuk's (2002). The model was recently expanded by the addition of state variables for the marine carbonate system as well as particulate and dissolved organic carbon of both marine and terrestrial origin (Gustafsson et al., 2014a). 2.1.1. Hydrodynamical module Exchange of water, salt, heat, and dissolved and particulate constituents between sub-basins is forced by wind, fluctuating sea level and horizontal density differences between basins. Flows through the straits that connect basins are controlled by frictional resistance and dynamical flow contraction due to Bernoulli and Coriolis effects (Gustafsson, 2000, 2003). Within each basin, the vertical stratification is resolved by a variable number of layers shaped by water inflows of differing densities. Fusion of layers keeps the total number of layers below a prescribed maximum value (Gustafsson, 2000). Vertical mixing is described by a mixed-layer model for the surface layer (Stigebrandt, 1985) together with a wind and stratification dependent parameterization for deepwater mixing (Axell, 1998; Stigebrandt, 1987). Further, the entrainment of surrounding water into dense gravity currents is simulated according to Stigebrandt (1987). The model includes a sea-ice model (Björk, 1997; Nohr et al., 2009) and bulk formulas for heating/cooling and evaporation at the sea surface (Björk, 1997; Gustafsson, 2003). At the open boundary between the Baltic and North Seas, time dependent profiles for all state variables comprise the lateral boundary conditions for BALTSEM. 2.1.2. Biogeochemical module The biogeochemical module includes oxygen, hydrogen sulfide and four inorganic nutrients; ammonium, nitrate, phosphate and silicate. DIC and total alkalinity (TA) are included as state variables as well, enabling calculations of pH, pCO2, and air–sea CO2 exchange. DIC is composed of dissolved CO2, carbonic acid, bicarbonate and carbonate. Autotrophs are described as three functional groups: diatoms, flagellates and cyanobacteria (Savchuk, 2002; Savchuk et al., 2012). The cyanobacteria group has the capacity to fix atmospheric nitrogen gas and is consequently not limited by dissolved nitrogen in the productive layer. One functional group for heterotrophs represents all organisms that feed on phytoplankton and detritus. Dead organic matter is included in the forms of dissolved and particulate organic carbon, nitrogen and phosphorus as well as biogenic silica (diatom shells). Detrital carbon is described by two state variables, separating carbon from terrestrial and marine sources respectively (Gustafsson et al., 2014a). The same distinction is made for DOC; autochthonous and allochthonous DOC are then further separated into refractory and bioavailable fractions — in all, four state variables describe DOC in the model. Allochthonous DOC is imported by rivers and atmospheric deposition, autochthonous DOC is produced by phytoplankton exudations, zooplankton excretions as well as solubilization of detrital carbon. A fraction of the refractory DOC is phototransformed into labile DOC. Further, a fraction of the DOC is transformed into POC by flocculation (Gustafsson et al., 2014a). Particulate organic matter sinks through the water column and is added to the depth dependent sediment pools of organic nitrogen and phosphorus, biogenic silica, as well as autochthonous and allochthonous organic carbon. Temperature dependent mineralization of organic matter occurs both in the water column and in the sediments. Resuspension, the following near-bottom transports and subsequent resettlement is parameterized as a downwards transport of organic matter (Savchuk et al., 2012). A fraction of the settled organic matter is permanently buried in the sediments.

123

An oxygen dependent fraction of the phosphate resulting from benthic mineralization of organic phosphorus is retained in the sediments. The remaining fraction is released into the overlying water. Further, an oxygen dependent fraction of the nitrate produced by mineralization and nitrification in the sediments is released to the water column whereas the remainder is denitrified. Under anaerobic conditions all mineralized nitrogen is released as ammonium. If the oxygen concentration approaches zero, first nitrate and then sulfate will be used to oxidize the organic matter (denitrification and sulfate reduction). Mineralization of organic carbon produces CO 2 which then equilibrates with bicarbonate and carbonate — depending on e.g. temperature and TA. Most biogeochemical processes result in either a production or consumption of TA. Internal TA sources and sinks are calculated according to Wolf-Gladrow et al. (2007). Additional TA generation resulting from unresolved processes in the sediment is included according to Gustafsson et al. (2014b). Simulated nutrient concentrations and fluxes have been discussed in detail and validated against observations by Savchuk et al. (2012). Simulated parameters of the carbonate system (TA, DIC, pH and pCO2) as well as the concentrations of autochthonous and allochthonous DOC respectively were validated against observations by Gustafsson et al. (2014a). 2.1.3. Carbon isotopes In this study, the model has been further expanded by adding one extra state variable for each carbon containing state variable (DIC, all plankton groups, terrestrial and marine detrital and dissolved organic carbon; cf. Gustafsson et al., 2014a). The carbon fractionation resulting from air–sea fluxes, biogeochemical processes and external loads are described in Sections 2.2–2.4. 2.1.4. Model forcing The hydrodynamical module is forced by 3-hourly meteorological data: wind, air temperature, humidity, cloudiness, pressure and precipitation. Further, daily average sea level in the Kattegat and stratification at the Skagerrak boundary are needed as well as monthly average runoff data. Monthly average river loads and atmospheric deposition of inorganic and organic nutrients and carbon as well as DIC and TA are together with profiles at the open model boundary needed for the biogeochemical modelling. Historical forcing data, covering the period 1850–2006, are described in detail by Gustafsson et al. (2012). River loads of DIC, DOC and TA are based on all available observations in 1996–2000 (cf. Gustafsson et al., 2014b). These data were used to calculate average concentrations of riverine DIC, DOC and TA. The monthly loads thus depend on monthly freshwater input multiplied by the average concentrations. Model forcing files have been updated to present day (December 2014): Meteorological forcing files for 1999–2014 were created from MESAN data; these data were downloaded from databases provided by the Swedish Meteorological and Hydrological Institute, SMHI (http:// opendata-download-grid-archive.smhi.se/explore/?modeltype=4). Further, the average daily Kattegat water level was updated to 2014. Data were downloaded from databases provided by the SMHI (http://opendata-download-ocobs.smhi.se/explore/?parameter= 4#). River runoff and nutrient loads have been based on PLC data (Pollution Load Compilation; cf. Svendsen et al., 2013) been updated to year 2012. Runoff data for 2013–2014 were not available, for that reason we repeated the monthly loads from 2012 also in the years 2013–2014. 2.2. Air–sea CO2 exchange Gas exchange between sea surface and atmosphere is controlled by the air–sea difference in CO2 partial pressure and by the wind speed

124

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

dependent transfer velocity. The exchange (mol kg−1 cm h−1) can be written as follows (e.g. Wanninkhof et al., 2009):

Freeman and Hayes, 1992; Hinga et al., 1994; Hollander and McKenzie, 1991). The isotopic value for aqueous CO2 is calculated from

FðCO2 Þ ¼ kwCO2 K0 ðpCO2a −pCO2w Þ:

δ CCO2ðaqÞ ¼ δ CDIC −f HCO3 εHCO3‐g −f CO3 εCO3‐g

ð1Þ

Here, K0 (mol kg−1 atm−1) is the CO2 gas solubility (Weiss, 1974), whereas pCO2a and pCO2w (atm) are CO2 partial pressures in air and surface water respectively. The parameterization for transfer velocity follows that by Weiss et al. (2007) as follows: kwCO2

rffiffiffiffiffiffiffiffiffi 660 ¼ k660 : Sc

ð2Þ

Here, k660 is the normalized transfer velocity (Wanninkhof et al., 2009), and Sc is a temperature dependent Schmidt number (Wanninkhof, 1992). The fractionation calculations associated with CO2 exchange between air and sea follow Zhang et al. (1995). Air–sea exchange of 13 C-CO2 can thus be written as follows: F

13

C−CO2 Þ ¼ kwCO2 K0 αk αaq‐g pCO2a Ra −pCO2w

! RDIC : αDIC‐g

ð3Þ

The fractionation factors α in Eq. (3) are related to per mil fractionations ε according to: αk ¼

εk þ1 1000

αaq‐g ¼

ð4Þ

εaq‐g þ1 1000

αDIC‐g ¼

ð5Þ

εDIC‐g þ 1: 1000

ð6Þ

These ε-values are adopted from Zhang et al. (1995) (cf. Table S1). RDIC and Ra are 13C/12C ratios for DIC and atmospheric CO2 respectively. Due to fossil fuel burning, the isotopic value of atmospheric CO2 is time dependent and has undergone a decrease from a preindustrial value of approximately −6.4‰ to the present-day value of −9‰ (as of 2014). Based on the data by Francey et al. (1999), we calculate this value as a linear function of atmospheric CO2 partial pressure as follows: 13

6

δ CCO2a ¼ −0:0222  pCO2a  10 :

ð7Þ

In Figure S2, reconstructed past atmospheric pCO2 values are displayed together with the corresponding isotopic values according to Eq. (7). 2.3. Production and mineralization A Rayleigh process is used to describe the fractionation during primary production (cf. Figure S3). 13

13

δ CPOCðinst:Þ ¼ εprod þ δ CCO2ðaqÞ 13

RPOCðinst:Þ ¼

δ CPOCðinst:Þ 1000

13

ð10Þ

where fHCO3 and fCO3 are the instantaneous bicarbonate and carbonate fractions of DIC. The isotopic value for DIC is 13

δ CDIC ¼

  RDIC −1  1000 RStd

ð11Þ

where Rstd is the PDB standard. δ13CDIC is continuously undergoing changes because of processes such as production, mineralization, air– sea exchange, diapycnal mixing, river loads, etc. Silicon fractionation has been observed during dissolution of biogenic silica (Sun et al., 2014). Here it is assumed that there is no carbon fractionation associated with grazing or mineralization of organic carbon. Hence, organic carbon is consumed/respired in accordance with its isotopic composition. Changes in 12C-DIC and 13C-DIC because of pelagic processes can be written as follows:   d 12 C−DIC dt   13 d C−DIC dt

¼ wDOCm þ wDOCt þ eZOO −aPHY

ð12Þ

¼ wDOCm  RDOCm þ wDOCt  RDOCt þ eZOO  RZOO −aPHY  RPOCðinst:Þ :

ð13Þ

Here, wDOCm and wDOCt indicate mineralization of marine and terrestrial DOC respectively, eZOO is DIC excretion by heterotrophs and aPHY is DIC assimilation by autotrophs. The 13C/12C-ratios for DOCm, DOCt, heterotrophs and autotrophs are denoted by RDOCm, RDOCt, RZOO, and RPOC(inst.) respectively. 2.4. External loads For the deep water inflows of DOCm and DOCt at the Skagerrak– Kattegat boundary we prescribe typical isotopic values for marine and terrestrial DOC respectively: δ 13 C DOCm-SK = − 21‰ and δ13CDOCt-SK = − 28‰ (Deutsch et al., 2012). The present-day isotopic value of DIC at the Skagerrak–Kattegat boundary, δ13 CDIC-SK, is assumed to be 1.5‰ (cf. Quay et al., 2007). The rate of change for δ13 CDIC-SK due to the Suess effect is set to − 0.26‰ per decade (Körtzinger et al., 2003). Isotopic values of riverine DIC can vary considerably depending on properties of the catchment, but also varies within a catchment as well as in individual rivers over a season (Giesler et al., 2013). Isotopic values for riverine DOC exhibit comparatively smaller variations and are typically found to be in the range from −28‰ to −27‰ (Deutsch et al., 2012). The assumed isotopic values for riverine DIC and DOC (δ13CDICt and δ13CDOCt) in different sub-basins are listed in Table 1. We do not know the temporal development of isotopic values for riverine DIC and DOC, although it appears reasonable that at least the isotopic

ð8Þ

! þ 1 RStd

13

ð9Þ

Here, δ13CPOC(inst.) and RPOC(inst.) are the instantaneous isotopic value and ratio respectively during primary production. The fractionation between CO2 and the instantaneous product is for simplicity assumed to be constant (εprod = −13‰, cf. Alling et al., 2012) although fractionation by phytoplankton can depend on phytoplankton species, ambient CO2 concentration, pH, and temperature (Fogel and Cifuentes, 1993;

Table 1 Isotopic values for riverine DIC and DOC supplied to the different sub-basins (own estimates). Sub-basin

δ13CDICt (‰)

δ13CDOCt (‰)

Kattegat Danish Straits Baltic Proper Bothnian Sea Bothnian Bay Gulf of Riga Gulf of Finland

−10 −10 −10 −12 −14 −8 −10

−28 −28 −28 −28 −28 −28 −28

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

values for riverine DIC should to some extent be affected by the decreasing isotopic value for atmospheric CO2. 2.5. Data DIC concentrations are calculated by using observed pH, TA, temperature, salinity etc. The measured parameters were extracted from the SHARK data base provided by the SMHI. In addition, our own measurements of δ13CDIC (cf., Table S2) were used to validate model results (Section 4.1). The sampling and preservation techniques were similar to those described in Alling et al. (2012). Niskin bottles attached to a Seabird CTD were used to collect the samples. After filtration, the samples were directly injected with a syringe into 12 ml septum-seal glass vials (Labco Limited) which had been flushed with argon gas (75 ml min−1) for 5 min. For each sample, duplicate aliquots (1–4 ml) were taken. 100 μl of 85.5% phosphoric acid was added to each vial to act as a preservative and to transform all the bicarbonate and carbonate to CO2(g). Both air entering the vial and CO2 leaving were avoided by means of transfer by needle injection through the septa. The samples were stored under cold (+ 4 °C) and dark conditions until analysis. Stable carbon isotopic compositions were determined using a Gasbench II extraction line coupled to a Finnigan MAT 252 mass spectrometer. 3. Results In Fig. 1, simulated DIC concentrations in the Gotland Sea (sub-basin 9, Figure S1) are displayed together with corresponding δ13CDIC values. DIC assimilation in the surface layer during the productive season is accompanied by increasing isotopic values in accordance with the Rayleigh process (Figure S3). In the deep water, extended stagnation periods are characterized by increasing DIC concentrations and decreasing isotopic values due to mineralization of organic carbon. DIC vs. δ13CDIC values are presented in Fig. 2. Mineralization is the main process affecting isotopic values in the deep, and there is a clear linear relationship between the two properties. The slope of the line depends on the isotopic composition of the organic carbon that is mineralized. Thus, the larger the contribution of terrestrial organic carbon, the steeper the slope becomes (the isotopic value decreases faster during a stagnation period). In the surface layer the isotopic composition is in addition to mineralization affected by river loads, primary production and air– sea CO2 exchange. The development of DIC vs. δ13CDIC values in surface water over one productive season is further explored in Figure S4. DOC concentrations and δ13CDOC values are presented in Fig. 3. The seasonal cycle in surface water DOC is coupled to respiration, phytoplankton exudation, as well as zooplankton excretion and sloppy feeding. The seasonality of the isotopic value is connected to the corresponding δ13CDIC value but less pronounced since a large fraction of the DOC is of terrestrial origin (Deutsch et al., 2012; Gustafsson et al., 2014a) which is assumed to have a constant isotopic value (−28‰).

125

Surface water POC in the model is dominated by phytoplankton bloom events, although there is also a small contribution of POC from terrestrial sources. In summer, the δ13CPOC value is close to the isotopic value for phytoplankton, whereas in winter the isotopic value approaches that of allochthonous POC due to the insignificant production of autochthonous POC (Fig. 4). The isotopic value for primary producers varies considerably over the productive season because of the Rayleigh process. The fractionation between CO2 and phytoplankton has in this study been assumed to be constant (εprod = − 13‰). The fractionation between DIC and phytoplankton, εPHY-DIC, on the other hand exhibits seasonal variations (Figure S5) due to the temperature dependent fractionation between CO2 and DIC (cf. Sections 2.2–2.3). 4. Discussion This new model formulation has the potential to improve our understanding of e.g. CO2 assimilation and mineralization rates, air–sea CO2 fluxes, and the impact of allochthonous DOC mineralization in different Baltic Sea sub-basins. However, because of a general shortage of observational data, our simulated isotopic values cannot as yet be validated in any detail. A few measured δ13CDIC profiles are available from 2012 (own data, Table S2) — a comparison between the simulated and available measured δ13CDIC values is presented in Section 4.1. In Section 4.2, the sensitivity of δ13CDIC to adjustments of various parameters is explored. Implications on a long-term scale are discussed in Section 4.3. 4.1. Measured vs. modelled isotopic values Measured DIC concentrations from 2012 in the Gotland Sea are somewhat higher than simulated values in 2012 from the same area, especially in the deep water (Fig. 5). The measured isotopic values lower on the other hand are lower (except at the surface). The deep water bias can partly be explained by underestimated deep water salinity and TA in the model (not shown). Thus, the simulated input of highly saline and alkaline water from the North Sea was insufficient in this period, resulting in a discrepancy also in the DIC concentrations and isotopic values. The simulated DIC concentrations in the Bothnian Sea are overall slightly underestimated, whereas the δ13CDIC values are within the observed range (Fig. 6). In the Bothnian Bay on the other hand, we find that although there is no clear bias between measured and simulated DIC concentrations, the simulated δ13CDIC values exceed measured ones by approximately 1.5‰ (Fig. 7). It has earlier been recognized that Baltic Sea physical–biogeochemical models fail to reproduce nutrient and carbon cycling in the Bothnian Bay, resulting in an underestimated productivity by at least a factor two (e.g. Eilola et al., 2011; Gustafsson et al., 2014a; Savchuk et al., 2012). If POC production is underestimated, then the export and deep water mineralization of organic carbon are

Fig. 1. Simulated surface and deep water DIC concentrations (μmol kg−1) and isotopic values (‰) in the Gotland Sea.

126

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

Fig. 2. Simulated surface and deep water DIC (μmol kg−1) vs. δ13CDIC (‰) in the Gotland Sea over the 2000–2014 period.

underestimated as well. In a model experiment, the carbon assimilation in the Bothnian Bay was amplified by a factor two by increasing the recycling rate of phosphorus. The resulting increased uptake of atmospheric CO2 as well as POC transport to deeper layers had a comparatively small effect on the DIC concentrations, but on the other hand a very large effect on the δ13CDIC value (Fig. 7, dashed line). 4.2. Sensitivity study In this section we explore the model sensitivity to a number of different adjustments in isotopic values as well as river loads of carbon and nutrients. The different model adjustments performed in the sensitivity tests are listed in Table 2. Sensitivity to modifications of isotopic values. • In the present study, we have assumed that the DOC supplied from rivers (DOCt) has a constant isotopic value of −28‰. In Section 4.2.1, we examine the sensitivity of the system if this value is adjusted by ±2‰. • Further, the sensitivity to assumed isotopic values for riverine DIC (Table 1) is examined. In Section 4.2.2, the isotopic values for riverine DIC are adjusted by ±2‰. • In the 1960–2006 period, the δ13C for atmospheric CO2 went from about − 7.0‰ to − 8.5‰ due to fossil fuel burning (Figure S2). In Section 4.2.3 we examine the impact on δ13CDIC values in the Gotland Sea if the atmospheric value is adjusted by ±2‰. • The δ13CDIC for inflowing North Sea water is assumed to have a contemporary value of 1.5‰, and is further assumed to decline by 0.23‰ per decade (cf. Section 2.4). In Section 4.2.4 we examine the effect on δ13CDIC if the value is modified by ±2‰. • The fractionation between CO2 and organic carbon during primary production (εprod) has been assumed to have a constant value of − 13‰, although the fractionation can be a function of e.g. the CO2 concentration itself or surface water pH (Freeman and Hayes,

1992; Hinga et al., 1994; Hollander and McKenzie, 1991). In Section 4.2.5, this value is adjusted by ±2‰.

Sensitivity to alterations in river loads. • River loads of DIC and nutrients (which affect production) have an impact on the CO2 partial pressure and thus air–sea CO2 exchange. Changes in these loads will thus in addition affect the isotopic composition of inorganic and organic carbon. In Section 4.2.6 we examine the sensitivity to DIC loads by either increasing or decreasing the loads by 25%. • In Section 4.2.7 we instead either increase or decrease the loads of dissolved inorganic nitrogen and phosphorus (DIN and DIP) by 25%.

4.2.1. Sensitivity to isotopic values of riverine DOC As expected, if the isotopic value DOCt is increased or decreased, the isotopic value for DIC increases/decreases as well because of DOCt mineralization (tests 1–2, Fig. 8a). The difference compared to the standard model run is less than 0.1‰ throughout the water column. An adjusted isotopic composition of DOCt thus has a rather small effect according to our model simulation. The reason is that mineralized DOCt contributes only by a marginal fraction (0.52 Tg C y−1, Table 3) to the total DIC pool, which in the period 2000–2014 amounts to an average 240 Tg C in the Gotland Sea according to our model calculations. 4.2.2. Sensitivity to isotopic values of riverine DIC By increasing or decreasing the isotopic value of riverine DIC by 2‰ in all rivers (Table 2), the delta value for DIC increases/decreases throughout the water column (tests 3–4, Fig. 8a). The increase/decrease in δ13CDIC is slightly larger than in the former case (tests 1–2). The reason is that the contribution from riverine DIC to the DIC pool

Fig. 3. Simulated surface and deep water DOC concentrations (μmol kg−1) and isotopic values (‰) in the Gotland Sea.

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

Fig. 4. Upper panel: Simulated surface water POC concentration (μmol kg−1). Lower panel: isotopic values (‰) for POC (black line) and phytoplankton (grey line) in the Gotland Sea.

127

Fig. 6. Simulated average profiles from 2012 for DIC and δ13CDIC in the Bothnian Sea. DIC and δ13CDIC observations from 2012 (full circles) are included.

(1.9 Tg C y−1, Table 3) is larger than the contribution from mineralized DOCt.

4.2.3. Sensitivity to the isotopic value of atmospheric CO2 A higher or lower isotopic value of atmospheric CO2 has a major impact on δ13CDIC in the system. Thus, when the isotopic value of atmospheric CO2 is adjusted by ±2‰ compared to the standard model run, we see a similar change of δ13CDIC in the water column (tests 5–6, Fig. 8b). The long term impact of the Suess effect is further explored in Section 4.3.

4.2.4. Sensitivity to the isotopic value of inflowing North Sea water By increasing or decreasing the isotopic value of DIC in inflowing North Sea water by 2‰, the delta value for DIC again increases/decreases throughout the water column (tests 7–8, Fig. 8c). The effect is as expected more pronounced in the deep water where the contribution from Skagerrak water is larger than in the surface water. A long-term influence is anticipated as well, as the isotopic value of inflowing Skagerrak water gradually decreases because of the Suess effect (cf. Section 2.4).

Fig. 5. Simulated average profiles from 2012 for DIC and δ13CDIC in the Gotland Sea. DIC and δ13CDIC observations from 2012 (full circles) are included.

4.2.5. Sensitivity to fractionation during primary production Compared to observations of isotopic values for phytoplankton (cf. Rolff, 2000; Zohary et al., 1994), the simulated range (δ13C ~ −23‰ to − 20‰; Fig. 4) is small. This could be explained by the fact that a constant fractionation between CO2 and phytoplankton POC is used during primary production. Production/mineralization of organic carbon is associated with large carbon flows (Table 3). Thus, δ13CDIC is sensitive to the fractionation during primary production. By increasing εprod from − 13‰ to − 11‰ (test 9), we expect the δ13CDIC of surface water to decline since the preferential 12C uptake — and removal via sedimentation — is less pronounced than in the standard model run. The decreased 12C removal is however to some extent compensated by air–sea CO2 fluxes. In the deep water on the other hand, there are no compensating processes. Thus, mineralization of the comparatively isotopically heavier organic carbon in the deep water sediments results

Fig. 7. Simulated average profiles from 2012 for DIC and δ13CDIC in the Bothnian Sea. DIC and δ13CDIC observations from 2012 (full circles) are included. For comparison, the dashed lines indicate the resulting average DIC and δ13CDIC when CO2 assimilation was amplified by a factor two.

128

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

Table 2 Sensitivity experiments. Test

Parameter

Adjusted isotopic value

Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 Test 9 Test 10 Test 11 Test 12 Test 13 Test 14

δ-Value of riverine DOC δ-Value of riverine DOC δ-Value of riverine DIC δ-Value of riverine DIC δ-Value of atmospheric CO2 δ-Value of atmospheric CO2 δ-Value of North Sea DIC δ-Value of North Sea DIC Fractionation during primary production Fractionation during primary production DIC loads × 1.25 DIC loads × 0.75 DIN & DIP loads × 1.25 DIN & DIP loads × 0.75

δ13CDOCt +2‰ δ13CDOCt −2‰ δ13CDICt +2‰ δ13CDICt −2‰ δ13CCO2a +2‰ δ13CCO2a −2‰ δ13CDIC-SK +2‰ δ13CDIC-SK −2‰ εprod +2‰ εprod −2‰ – – – –

in a clear increase in δ13CDIC compared to the standard model run (test 9, Fig. 8c). If εprod is instead decreased by 2‰, the effect on δ13CDIC values is the opposite. 4.2.6. Sensitivity to changes in DIC loads Riverine DIC typically has considerably lower isotopic values than surface water in the Baltic Sea (Table 1, Fig. 1). Thus, if the DIC loads are increased by 25%, this adjustment alone would result in decreased δ13CDIC values in the system. However, an increased DIC load also has the effect that the CO2 partial pressure in surface water increases. This in turn results in a diminished absorption/enhanced outgassing of CO2, and it is preferentially 12C that is exchanged between air and sea. The combined effect of these two processes is a slight increase of δ13CDIC values in the water column. If DIC loads instead are decreased by 25%, the effect is the opposite (tests 11–12, Fig. 8d). 4.2.7. Sensitivity to changes in DIN and DIP loads If DIN and DIP loads are increased by 25%, this results in an enhanced primary production and thus increased CO2 assimilation and POC export to deeper water and sediments. Since 12C is preferentially assimilated, the immediate effect of enhanced primary production is that δ13CDIC increases (cf. Fig. 1, Figure S3). By increasing the nutrient loads compared to the standard model simulation, the preferential 12C assimilation, and by that also the δ13CDIC values, are enhanced. Another effect of increased primary production is that during the productive season the surface

Table 3 Processes affecting the DIC pool and δ13CDIC values in the Gotland Sea. Average fluxes over the 2000–2014 period (Tg C y−1). Process

12

13

Riverine DIC load Net DIC export Net uptake of atmospheric CO2 CO2 assimilation by phytoplankton Pelagic mineralization of autochthonous organic carbon Pelagic mineralization of allochthonous organic carbon Sediment release of DIC (marine origin) Sediment release of DIC (terrestrial origin)

1.9 3.4 1.2 21 17 0.52 4.3 0.23

0.021 0.039 0.014 0.23 0.18 0.0056 0.048 0.0025

C flux

C flux

water CO2 partial pressure decreases (more than in the standard model run) and the absorption of atmospheric CO2 thus increases. According to our results the net effect of these two processes is that in the surface water, increased absorption of atmospheric CO2 largely compensates the loss of 12C associated with enhanced primary production and POC sedimentation. In the deep water on the other hand, there is no compensating process and the increased POC input and mineralization thus results in decreased δ13CDIC values. If the nutrient loads are instead decreased by 25%, the result is the opposite (tests 13–14, Fig. 8d). 4.3. Implications Although several processes influence both DIC concentration and isotopic composition, we found that δ13CDIC values in the water column are mainly controlled by primary production/mineralization and the associated CO2 exchange between air and sea due to changes in CO2 partial pressure (cf. Broecker and Maier-Reimer, 1992). The reason is that the carbon flows coupled to production and mineralization by far exceed other carbon flows such as river loads and DOCt mineralization (Table 3). We further found that the isotopic value of atmospheric CO2 has a major impact on the system. To further explore the long-term effects of changing productivity and the Suess effect, the model sensitivity in terms of isotopic composition will be evaluated in this section in relation to eutrophication as well as the long-term changes in δ13C of CO2 in the atmosphere and DIC in inflowing North Sea water. Two different model simulations covering the years 1940–2006 will be compared. In the first simulation (standard loads) we use reconstructed nutrient loads (cf. Gustafsson et al., 2012)

Fig. 8. Difference in δ13CDIC between the different experimental model runs (Table 2) and the standard model simulation. Average profiles from the Gotland Sea in 2000–2014 are displayed. (a) Full black and grey lines indicate results from tests 1 and 2, whereas dashed black and grey lines indicate results from tests 3 and 4 respectively. (b) Full black and grey lines indicate results from tests 5 and 6 respectively. (c) Full black and grey lines indicate results from tests 7 and 8, whereas dashed black and grey lines indicate results from tests 9 and 10 respectively. (d) Full black and grey lines indicate results from tests 11 and 12, whereas dashed black and grey lines indicate results from tests 13 and 14 respectively.

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

whereas in the second simulation nutrient loads are assumed to be at early 20th century levels throughout the simulation period (low loads). River loads of DIC, DOC and TA are identical in the two simulations. Fig. 9 illustrates the simulated development of DIC and δ13CDIC respectively in Gotland Sea surface and deep water. Variations in DIC concentrations prior to ~ 1970 depend mainly on changing runoff. In the case where standard loads were used, substantially increased DIC concentrations are obtained in both surface and deep water. Deep water DIC increases due to an amplified POC deposition (Fig. 10) and mineralization coupled to increased nutrient loads. The resulting oxygen depletion is associated with a strong internal TA generation (e.g. Krumins et al., 2013). The subsequent increasing TA concentration in the water column enhances the buffer capacity and thus the absorption of atmospheric CO2. Thus, surface water DIC increases in spite of the increased DIC assimilation and POC sedimentation coupled to eutrophication. The isotopic value of surface water DIC declines in both simulations because of the Suess effect (cf. Figure S2). In the standard case this decrease is however somewhat dampened as a result of the enhanced DIC assimilation, where 12C is preferentially removed from the productive layer. In the deep water on the other hand the isotopic value decreases more rapidly in the standard case because of the enhanced mineralization of organic carbon (which is depleted in 13C). In Fig. 10, simulated POC deposition fluxes in the Gotland Sea are shown from the two simulations. Results from the standard simulation clearly demonstrate the impact of increasing nutrient loads — POC deposition more than doubled in 1940–2006. The corresponding isotopic values behave differently in the two simulations. In the low loads simulation, the isotopic value of deposited POC decreases since the isotopic values of primary producers decline — which in turn is related to a decreasing isotopic value for DIC because of the Suess effect (Figure S2). In the standard simulation on the other hand, the isotopic value of deposited POC is compensated as the deposited material to a higher and higher degree is dominated by autochthonous POC which is isotopically heavier than allochthonous POC. 4.4. Future outlook The lack of measured δ13C of DIC, DOC and POC limits the possibility to validate model results — and thus limits the usability of the model. In order to validate the model more thoroughly, seasonally resolved measurements from the major Baltic Sea sub-basins would be required. I.e., monthly profiles of δ13C of DIC, DOC and POC from all the BALTSEM basins (Figure S1), covering a period of at least one year, would be highly useful. In addition, monthly measurements in the major rivers entering the system would further reduce uncertainties in the model simulations. We would thus highly recommend a measurement program focused on the seasonal dynamics of δ13C in different parts of the Baltic Sea system.

129

Fig. 10. Annual mean POC deposition (Tg C y−1) and corresponding isotopic values (‰) for deposited POC in the Gotland Sea from two simulations: standard nutrient loads (black lines) and low loads (grey lines).

5. Conclusions • BALTSEM is a coupled physical–biogeochemical model that recently was expanded to include DIC as well as DOC and POC from both terrestrial and marine sources (Gustafsson et al., 2014a). This version of BALTSEM has in the present study been developed further and now explicitly includes the 12C and 13C fractions of all carbon containing state variables. • Different processes in the transformations of organic and inorganic carbon have specific fractionation values. Thus, depending on the original composition as well as the dominant processes transforming the involved carbon species, particular isotopic fingerprints are formed. • The model can be used to calculate not only e.g. by DIC concentrations, but also by the help of simulated vs. measured isotopic values. It is further possible to estimate the dominant processes in the ecosystem that resulted in a certain concentration. • For the time being, only a few measurements of carbon isotopes are available, and these measurements (from 2012) are not covered by the simulation period which ends in 2006 due to a lack of contemporary forcing files. Nonetheless, our simulated time series for δ13C of DIC, DOC and POC appear reasonable. • In a sensitivity study we show that production and mineralization and the related air–sea CO2 fluxes are processes that largely control the isotopic composition in the system, since they also represent the main carbon flows. • The gradually decreasing isotopic value of atmospheric CO2 associated with fossil fuel burning has had a major impact on the isotopic composition of inorganic and organic carbon in the Baltic Sea under the 20th century — the δ13C of e.g. DIC has steadily decreased throughout the water column. In the deep water, this decrease has been further enhanced as a result of eutrophication. An amplified productivity

Fig. 9. Annual mean surface and deep water DIC concentrations (μmol kg−1) and isotopic values (‰) in the Gotland Sea from two simulations: standard nutrient loads (black lines) and low loads (grey lines).

130

E. Gustafsson et al. / Journal of Marine Systems 148 (2015) 122–130

and POC export to deeper layers result in an increased deep water mineralization of comparatively light organic carbon. • In order to validate the model more thoroughly and fully take advantage of the model's potential, further measurements of δ13C values in the system would be required. Monthly profiles from the different sub-basins, covering a period of at least one year, as well as measurements from the major rivers would highly improve the usability of the model.

Acknowledgements We acknowledge Heike Siegmund from Stable Isotope Laboratory (Stockholm University) for the analysis of stable isotopes of DIC. Development of BALTSEM is a part of the core activities of the Baltic Nest Institute supported by the Swedish Agency for Marine and Water Management. This project was in addition partly funded by FORMAS under the 215-2009-813 contract and by the strategic marine initiative Baltic Ecosystem Adaptive Management (BEAM) at Stockholm University. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.jmarsys.2015.02.008. References Alling, V., Porcelli, D., Mörth, C.-M., Anderson, L.G., Sanchez-Garcia, L., Gustafsson, Ö., Andersson, P.S., Humborg, C., 2012. Degradation of terrestrial organic carbon, primary production and out-gassing of CO2 in the Laptev and East Siberian Seas as inferred from δ13C values of DIC. Geochim. Cosmochim. Acta 95, 143–159. http://dx.doi.org/ 10.1016/j.gca.2012.07.028. Axell, L.B., 1998. On the variability of Baltic Sea deepwater mixing. J. Geophys. Res. 103, 21667–21682. http://dx.doi.org/10.1029/98JC01714. Björk, G., 1997. The relation between ice deformation, oceanic heat flux, and the ice thickness distribution in the Arctic Ocean. J. Geophys. Res. 102, 18681–18698. http://dx. doi.org/10.1029/97JC00789. Broecker, W.S., Maier-Reimer, E., 1992. The influence of air and sea exchange on the carbon isotope distribution in the sea. Global Biogeochem. Cycles 6, 315–320. http://dx.doi.org/10.1029/92GB01672. Deutsch, B., Alling, V., Humborg, C., Korth, F., Mörth, C.M., 2012. Tracing inputs of terrestrial high molecular weight dissolved organic matter within the Baltic Sea ecosystem. Biogeosciences 9, 4465–4475. http://dx.doi.org/10.5194/bg-9-4465-2012. Eilola, K., Gustafsson, B.G., Kuznetsov, I., Meier, H.E.M., Neumann, T., Savchuk, O.P., 2011. Evaluation of biogeochemical cycles in an ensemble of three state-of-the-art numerical models of the Baltic Sea. J. Mar. Syst. 88, 267–284. http://dx.doi.org/10.1016/j. jmarsys.2011.05.004. Fogel, M.L., Cifuentes, L.A., 1993. Isotope fractionation during primary production. Organic Geochemistry. Plenum, New York, pp. 73–98. Francey, R.J., Allison, C.E., Etheridge, D.M., Trudinger, C.M., Enting, I.G., Leuenberger, M., Langenfelds, R.L., Michel, E., Steele, L.P., 1999. A 1000-year high precision record of δ13C in atmospheric CO2. Tellus B 51, 170–193. http://dx.doi.org/10.1034/j.16000889.1999.t01-1-00005.x. Freeman, K.H., Hayes, J.M., 1992. Fractionation of carbon isotopes by phytoplankton and estimates of ancient CO2 levels. Global Biogeochem. Cycles 6, 185–198. http://dx. doi.org/10.1029/92GB00190. Giesler, R., Mörth, C.-M., Karlsson, J., Lundin, E.J., Lyon, S.W., Humborg, C., 2013. Spatiotemporal variations of pCO2 and δ13C-DIC in subarctic streams in northern Sweden. Global Biogeochem. Cycles 27, 176–186. http://dx.doi.org/10.1002/gbc.20024. Gustafsson, B.G., 2000. Time-dependent modeling of the Baltic entrance area. 1. Quantification of circulation and residence times in the Kattegat and the straits of the Baltic sill. Estuaries 23, 231–252. Gustafsson, B.G., 2003. A Time-dependent Coupled-basin Model of the Baltic Sea. C47. Earth Sciences Centre, Göteborg University, Göteborg. Gustafsson, B.G., Schenk, F., Blenckner, T., Eilola, K., Meier, H.E.M., Müller-Karulis, B., Neumann, T., Ruoho-Airola, T., Savchuk, O.P., Zorita, E., 2012. Reconstructing the de-

velopment of Baltic Sea eutrophication 1850–2006. Ambio 41, 534–548. http://dx. doi.org/10.1007/s13280-012-0318-x. Gustafsson, E., Deutsch, B., Gustafsson, B.G., Humborg, C., Mörth, C.-M., 2014a. Carbon cycling in the Baltic Sea — the fate of allochthonous organic carbon and its impact on air–sea CO2 exchange. J. Mar. Syst. 129, 289–302. http://dx.doi.org/10.1016/j. jmarsys.2013.07.005. Gustafsson, E., Wällstedt, T., Humborg, C., Mörth, C.-M., Gustafsson, B.G., 2014b. External total alkalinity loads versus internal generation: the influence of nonriverine alkalinity sources in the Baltic Sea. Global Biogeochem. Cycles 28, 1358–1370. http://dx.doi. org/10.1002/2014GB004888. Hinga, K.R., Arthur, M.A., Pilson, M.E.Q., Whitaker, D., 1994. Carbon isotope fractionation by marine phytoplankton in culture: the effects of CO2 concentration, pH, temperature, and species. Global Biogeochem. Cycles 8, 91–102. http://dx.doi.org/10.1029/ 93GB03393. Hollander, D.J., McKenzie, J.A., 1991. CO2 control on carbon-isotope fractionation during aqueous photosynthesis: a paleo-pCO2 barometer. Geology 19, 929. http://dx.doi. org/10.1130/0091-7613(1991)019b0929:CCOCIFN2.3.CO;2. Körtzinger, A., Quay, P.D., Sonnerup, R.E., 2003. Relationship between anthropogenic CO2 and the 13C Suess effect in the North Atlantic Ocean. Global Biogeochem. Cycles 17. http://dx.doi.org/10.1029/2001GB001427. Krumins, V., Gehlen, M., Arndt, S., Van Cappellen, P., Regnier, P., 2013. Dissolved inorganic carbon and alkalinity fluxes from coastal marine sediments: model estimates for different shelf environments and sensitivity to global change. Biogeosciences 10, 371–398. http://dx.doi.org/10.5194/bg-10-371-2013. Nohr, C., Björk, G., Gustafsson, B.G., 2009. A dynamic sea ice model based on the formation direction of leads. Cold Reg. Sci. Technol. 58, 36–46. http://dx.doi.org/10.1016/j. coldregions.2009.04.005. Quay, P., Sonnerup, R., Stutsman, J., Maurer, J., Körtzinger, A., Padin, X.A., Robinson, C., 2007. Anthropogenic CO2 accumulation rates in the North Atlantic Ocean from changes in the 13C/12C of dissolved inorganic carbon. Global Biogeochem. Cycles 21. http://dx.doi.org/10.1029/2006GB002761. Rolff, C., 2000. Seasonal variation in δ13C and δ15N of size-fractionated plankton at a coastal station in the northern Baltic proper. Mar. Ecol. Prog. Ser. 203, 47–65. http://dx.doi. org/10.3354/meps203047. Savchuk, O.P., 2002. Nutrient biogeochemical cycles in the Gulf of Riga: scaling up field studies with a mathematical model. J. Mar. Syst. 32, 253–280. http://dx.doi.org/10. 1016/S0924-7963(02)00039-8. Savchuk, O.P., Gustafsson, B.G., Müller-Karulis, B., 2012. BALTSEM — a marine model for the decision support within the Baltic Sea region (Technical Report No. 7). BNI Technical Report Series. Stigebrandt, A., 1985. A model for the seasonal Pycnocline in rotating systems with application to the Baltic proper. J. Phys. Oceanogr. 15, 1392–1404. http://dx.doi.org/10. 1175/1520-0485(1985)015b1392:AMFTSPN2.0.CO;2. Stigebrandt, A., 1987. A model for the vertical circulation of the Baltic deep water. J. Phys. Oceanogr. 17, 1772–1785. http://dx.doi.org/10.1175/1520-0485(1987)017b1772: AMFTVCN2.0.CO;2. Sun, X., Olofsson, M., Andersson, P.S., Fry, B., Legrand, C., Humborg, C., Mörth, C.-M., 2014. Effects of growth and dissolution on the fractionation of silicon isotopes by estuarine diatoms. Geochim. Cosmochim. Acta 130, 156–166. http://dx.doi.org/10.1016/j.gca. 2014.01.024. Svendsen, L.M., Staaf, H., Gustafsson, B., Pyhälä, M., Kotilainen, P., Bartnicki, J., Knuuttila, S., Durkin, M., 2013. Review of the Fifth Baltic Sea Pollution Load Compilation for the 2013 HELCOM Ministerial Meeting (Baltic Sea Environment Proceedings No. 141). Helsinki Commission—HELCOM. Wanninkhof, R., 1992. Relationship between wind speed and gas exchange over the ocean. J. Geophys. Res. 97, 7373–7382. http://dx.doi.org/10.1029/92JC00188. Wanninkhof, R., Asher, W.E., Ho, D.T., Sweeney, C., McGillis, W.R., 2009. Advances in quantifying air–sea gas exchange and environmental forcing. Annu. Rev. Mar. Sci. 1, 213–244. http://dx.doi.org/10.1146/annurev.marine.010908.163742. Weiss, R.F., 1974. Carbon dioxide in water and seawater: the solubility of a non-ideal gas. Mar. Chem. 2, 203–215. http://dx.doi.org/10.1016/0304-4203(74)90015-2. Weiss, A., Kuss, J., Peters, G., Schneider, B., 2007. Evaluating transfer velocity–wind speed relationship using a long-term series of direct eddy correlation CO2 flux measurements. J. Mar. Syst. 66, 130–139. http://dx.doi.org/10.1016/j.jmarsys.2006.04.011. Wolf-Gladrow, D.A., Zeebe, R.E., Klaas, C., Körtzinger, A., Dickson, A.G., 2007. Total alkalinity: the explicit conservative expression and its application to biogeochemical processes. Mar. Chem. 106, 287–300. http://dx.doi.org/10.1016/j.marchem.2007.01.006. Zhang, J., Quay, P.D., Wilbur, D.O., 1995. Carbon isotope fractionation during gas–water exchange and dissolution of CO2. Geochim. Cosmochim. Acta 59, 107–114. http:// dx.doi.org/10.1016/0016-7037(95)91550-D. Zohary, T., Erez, J., Gophen, M., Berman-Frank, I., Stiller, M., 1994. Seasonality of stable carbon isotopes within the pelagic food web of Lake Kinneret. Limnol. Oceanogr. 39, 1030–1043.