Increased oceanic dimethyl sulfide emissions in areas of sea ice retreat inferred from a Greenland ice core

Ocean phytoplankton are an important source of dimethyl sulfide, which influences marine cloud formation. Model studies suggest that declines in Arctic sea ice may lead to increased dimethyl sulfide emissions, however observational support is lacking. Here, we present a 55-year high-resolution ice core record of methane sulfonic acid flux, an oxidation product of dimethyl sulfide, from the southeast Greenland Ice Sheet. We infer temporal variations in ocean dimethyl sulfide emissions and find that springtime (April–June) fluxes of methane sulfonic acid correlate well with satellite-derived chlorophyll-a concentration in the Irminger Sea. Summertime (July–September) methane sulfonic acid fluxes were 3 to 6 times higher between 2002–2014 than 1972–2001. We attribute this to sea ice retreat day becoming earlier and a coincident increase in chlorophyll-a concentration in the adjacent open coastal waters. Summer sea ice retreat in the high latitude North Atlantic since 2000 is associated with increases in oceanic emissions of dimethyl sulfide, according to a record of methane sulfonic acid fluxes from an ice core in southeast Greenland.

D imethyl sulfide (DMS), which originates from ocean phytoplankton, is oxidized to non-sea-salt sulfate (nssSO 4 2− ) and methane sulfonic acid (MSA) after emission into the atmosphere 1 . These aerosols contribute to the formation of cloud condensation nuclei in the marine atmosphere and participate in cloud radiative forcing and direct aerosol radiative forcing 2,3 . As the oceanic flux of DMS accounts for 10-40% of the total global sulfur flux 4 and impacts the formation of cloud condensation nuclei in the marine atmosphere, DMS is an important factor used to calculate the radiative forcing from the climate model 5,6 .
DMS emissions from the ocean to the atmosphere are determined by the gas transfer coefficient from the ocean to the atmosphere, the solubility of gaseous DMS, and the DMS concentrations in water and air 7,8 . The gas transfer coefficient shows a linear dependence on mean horizontal wind speed to 11 m s −1 , but the relationship between the gas transfer coefficients and wind speed weakens at higher wind speeds 9,10 . The DMS in the ocean is converted from dimethylsulfoniopropionate (DMSP), a product of phytoplankton, by the enzyme system of bacteria 1 . Therefore, the concentration of DMS in the ocean is estimated from the chlorophyll-a (Chl-a) concentration. Gali et al. 11 estimated the DMS concentration at the sea surface at global and regional scales with a new remote sensing algorithm using remotely sensed Chl-a along with light penetration and climatological mixed-layer depth. The estimation of the DMS concentration at the sea surface by the new algorithm was improved at the global scale, but the estimation disagreed with the regional scale and seasonal variation. Gali et al. 12 noted that the shortcomings of the new algorithm were due to sparse and biased in situ DMS data.
Arctic temperatures are rising twice as fast as the global average, and the summer seasonal sea ice extent has declined sharply in recent decades 13,14 . In the Arctic region, the decline in seasonal sea ice cover increases the light penetration on the surface and subsurface of the ocean and promotes the primary production of phytoplankton 12,14-17 . In the satellite and modelbased estimation, the DMS emissions increased in the last decade because of the decline in the seasonal sea ice cover in the pan-Arctic 12,17 . However, there is little observational evidence for the recent increase in DMS emission in the pan-Arctic. To estimate the changes of the impact of DMS emission on the global climate, long-term and continuous aerosol monitorings are needed.
When DMS is released from the ocean to the atmosphere, DMS is oxidized to MSA and non-sea-salt sulfate. The MSA in aerosols is used as an indicator of oceanic DMS emission because the natural origin of MSA is almost exclusively limited to the oxidation of DMS. Sharma et al. 18 showed that aerosol MSA concentrations at three stations in the Arctic (Alert-Nunavut, Barrow-Alaska, and Ny-Ålesund-Svalbard) increased since 2000 with respect to 1990, and suggested that the increase in MSA resulted from the northward migration of the marginal ice edge zone. Except for the 3 stations, there are no long-term continuous monitoring data for aerosol MSA in the pan-Arctic region 19 . Ice cores drilled from ice sheets and glaciers are useful archives for reconstructing past aerosol history 20,21 . Aerosol MSA is preserved in ice sheets and glaciers after wet or dry deposition from the atmosphere to the snow surface 22 , but the MSA in the snowpack at low accumulation areas is released from the snow surface to the atmosphere, called postdepositional loss 23,24 . Therefore, drilling sites must be selected at high accumulation and cold areas of ice sheets and glaciers.
The southeastern dome of the Greenland Ice Sheet (SE-Dome; 67.18°N, 36.37°W, 3170 m above sea level, Fig. 1) has the distinct characteristic of having a high accumulation rate (1.01 ± 0.22 m water equivalent yr −1 ) and low annual mean temperature (−20.9°C) 21,25,26 . The accumulation rate is the highest among the domes on the polar ice sheets, at 3-4 times that of typical domes of Greenland Ice Sheet 25,26 . Therefore, even nitrate, which is a volatile species, was preserved in the SE-Dome ice core without any effects of postdepositional loss 21 , the dating of the SE-Dome ice core was determined by a comparison of δ 18 O records from the ice core and a record simulated by isotopeenabled climate models with an accuracy of 1 month in most period 26 (Method). Consequently, the ice core obtained from the SE-Dome allows us to reconstruct the seasonal variations in MSA from 1960 to 2014 21,26 . The objective of this study is to reveal the temporal variations in DMS emissions from the ocean to the atmosphere around the seasonal sea ice region using the hightime-resolution MSA profile in the SE-Dome ice core.

Results and discussion
Decadal and seasonal variation in the MSA flux. The mean annual flux of MSA (MSA flux ) (Method) from 2002-2014 was substantially higher than that from the complete study period (1960-2014) ( Fig. 2a and Table 1). On the other hand, the mean of annual MSA flux from 1972-2001 was substantially lower than the mean of the entire period ( Table 1). The interannual trend of annual MSA flux substantially decreased from 1960 to 2001 (−0.0061 mg m −2 yr −1 , p < 0.05). The decrease trend in the MSA flux in the SE-Dome ice core from 1960-2001 is consistent with the variations in MSA concentration in ice cores from the inland Greenland Ice Sheet, which showed declining trends over the late nineteenth and twentieth centuries [27][28][29] . However, the increase in the MSA flux in the SE-Dome ice core after 2002 has not been found in other ice cores. We divided the entire period into three periods (1960-1971, 1972-2001, and 2002-2014) to investigate the behavior of MSA for the following discussion.   (Fig. 3, black). The July to September MSA flux from 2002-2014 was 3-6 times higher than that from 1972-2001 (Fig. 3). In the northern part of the North Atlantic, phytoplankton blooms start during April-June offshore of the north of 45°N and during July-August in the seasonal sea ice area on the southeastern coast of Greenland 30 . Based on the seasonal variation in the MSA flux and the phytoplankton bloom, we defined April-June as "spring", July-September as "summer", October-December as "autumn", and January-March as "winter".
Most MSA (76%) is deposited on the snow surface by wet deposition 31 . Therefore, a high annual snow accumulation results in a high amount of MSA deposition. The reconstructed seasonal accumulation rate showed no seasonal difference (spring: 0.25 ± 0.07 m, summer: 0.25 ± 0.09 m, autumn: 0.27 ± 0.09 m and winter: 0.25 ± 0.07 m) 26 . Moreover, the accumulation rates in spring and summer did not change from 1960 to 2014 26 . Consequently, we suggest that the recent increase in summer MSA flux was caused by the increase in MSA concentration in the atmosphere from the changes in the seasonality of phytoplankton blooms and/or DMS emissions in the northern part of the North Atlantic, and not by the changes in precipitation seasonality. Next, we discuss the process of the MSA flux variation in the spring (April-June), summer (July-September), autumn (October-December), and winter (January-March).
Spring MSA flux. The spring MSA flux did not substantially change from 1960-2014 (+0.0002 mg m −2 yr −1 , p = 0.8888) (Fig. 2b). Important factors controlling DMS emissions from the ocean to the atmosphere include the DMS concentration at the ocean surface and wind stress on the ocean 7,8 . High concentrations of DMS at the ocean surface enhance the emission of DMS from the ocean to the atmosphere 8,32 . The gas transfer coefficient for DMS emission is exhibit the highest values at intermediate wind speeds (8-11 m s −1 ) 9,10 . We investigated the contributions of these factors to the variation in MSA flux in the SE-Dome ice core described below. The correlation coefficient between the aerial average of Chl-a (Methods) in the Irminger Sea (55°N-62°N, 22°W-38°W, dashed line box in Fig. 4a, c) and spring MSA flux was 0.69 (p < 0.01) (Fig. 4a, b). On the other hand, the spring MSA flux did not correlate with the frequency of intermediate wind (Methods) in the Irminger Sea (r = 0.04, p = 0.87) (Fig. 4c, d). A backward trajectory analysis (Methods) showed that transportation of air masses from the possible source region of MSA to the SE-Dome has been unchanged (Supplementary Fig. 1). Consequently, the spring MSA flux in the SE-Dome ice core was dominated by the Chl-a concentration in the Irminger Sea.
Summer MSA flux. The summer MSA flux substantially increased from 1960 to 2014 (+0.0037 mg m −2 yr −1 , p < 0.05) (Fig. 2c). Moreover, the summer MSA flux from 2002-2014 was 3-6 times higher than that from 1972-2001 (Fig. 3). The spatial distribution of air mass transportation was stable for the entire period (Supplementary Fig. 2). We can assume that the recent high values of summer MSA flux were not caused by changing transportation processes but by changes in phytoplankton production, sea ice conditions, and wind stress along the southeastern coast of Greenland. However, distribution maps of the correlation coefficient between the summer MSA flux and Chl-a concentration, frequency of intermediate wind, and sea ice concentration show that these indices did not correlate with the summer MSA flux (Supplementary Fig. 3). Although the spring MSA flux was quantitatively correlated with the Chl-a concentration in the Irminger Sea (Fig. 4a, b), the summer MSA flux did not correlate with any factors.
We focused on a sea ice retreat day for each period (1960-1971, 1972-2001, 2002-2014 Fig. 5b). The appearance frequency of higher Chl-a concentrations in CDF40 region has clearly increased in the 2002-2014 (Fig. 6). Therefore, we conclude that the summer DMS emissions due to the high production of phytoplankton have been enhanced and thus preserved as MSA in the SE-Dome ice core.
We can speculate why the summer MSA flux from 2002-2014 increased 3-6 times compared to that from 1972-2001 as follows.   5 and Supplementary Fig. 5a).
The lack of sea ice during July increases the light penetration at the surface and subsurface, and then enhances the production of phytoplankton and lengthens the bloom period. The summer MSA flux in the ice core results from the early retreat of the sea ice in CDF40 followed by an increase in DMS emissions from enhanced phytoplankton production.   (Fig. 3). Commonly, phytoplankton growth in mid-and high-latitude ocean in winter is limited because of deep mixing layer and low solar radiation. However, recently widespread winter phytoplankton blooms were observed in a large part of the North Atlantic sub-polar gyre triggered by intermittent restratification of the mixed layer when mixed-layer eddies led to a horizontal transport of lighter water over denser layers 34 . We speculate that the recent increase in MSA flux during winter could result from the recent observed phytoplankton bloom in winter in the North Atlantic.
Implication. The results of SE-Dome ice core analysis showed the 3-6 times increase in summer MSA flux since 2002. We suggested that the reason for the increase in summer MSA flux was the early sea ice retreat followed by an increase in DMS emissions from enhanced phytoplankton production. The similar early retreat of sea ice and the increase in Chl-a concentration also observed in offshore of Ilulissat in Baffin Bay after 2002 ( Fig. 5 and Supplementary Fig. 5). Additionally, process-based studies indicated that the accelerated decline of the sea ice in the pan-Arctic region since the 2000s could be associated with upward trends in oceanic DMS flux in the pan-Arctic 12,17 . Therefore, the increase in oceanic DMS emission in summer, which was the new observational truth by the SE-Dome ice core, could already occur in other seasonal sea ice areas. Under the future climate warming, we expect that the ocean areas increasing in the summer DMS emission were expanded. This suggests that negative cloud radiative forcing may be enhanced and mitigate the positive icealbedo feedback in the pan-Arctic because of increasing low clouds in the marine atmosphere. Although the spring bloom has been pointed as the most dominant factor for the oceanic natural sulfur emission until now, our study suggests that the process of the summer DMS emission and the direct and indirect effect of the radiation balance in accordance with the summer DMS emission would be more important for the estimation of the future radiation balance.

Conclusions
We reconstructed the annual and seasonal MSA flux from 1960 to 2014 with monthly resolution from a high-resolution ice core obtained from the SE-Dome, southeastern Greenland Ice Sheet. The annual MSA flux decreased from 1960 to 2001 and markedly increased after 2002 (Fig. 2). We investigated the most dominant factor for the recent increase in the annual MSA flux by defining April-June as "spring", July-September as "summer" October-December as "autumn", and January-March as "winter". The spring MSA flux was significantly correlated with the Chl-a concentration in the Irminger Sea (r = 0.69, p < 0.01) (Fig. 4). The summer MSA flux from 2002-2014 was 3-6 times higher than that from 1972-2001 (Fig. 3). This was the most dominant reason for the increase of annual MSA flux from 2002. We suggested that the reason for the increased summer MSA flux was the early sea ice retreat followed by the increase in oceanic DMS flux from the increase in phytoplankton production at the ocean surface and subsurface. The autumn MSA flux was more detected from 2002-2014 (Figs. 2d and 3). The winter MSA flux substantially increased from 1960 to 2014 (Figs. 2e and 3). The winter MSA in the SE-Dome ice core could be sourced by the winter phytoplankton bloom in the North Atlantic sub-polar gyre.
This ice core study concludes that ocean phytoplankton productivity and oceanic DMS emissions have increased since the early 2000s in seasonal sea ice areas. Data suggest that such increase in the oceanic DMS emissions could already occur in other regions where the early retreat of seasonal sea ice occurs. This suggests that negative cloud radiative forcing may be enhanced because of increasing low clouds in the marine atmosphere.

Methods
Dating of the SE-Dome ice core. The SE-Dome ice core was dated based on pattern matching of the δ 18 O variations between the ice core record and a simulated profile in Furukawa et al. 26 . The simulated profile of δ 18 O (δ 18 O model ) was constructed by a regional climate model with isotope diagnostics included in the hydrological cycle (REMO-iso) 35 and an atmospheric general circulation model (GCM), into which stable water isotopes are incorporated (iso-GSM) 36 . The SE-Dome δ 18 O and δ 18 O model variations were matched by selecting manually 170 tie points from 0.8 to 86.5 m depth using AnalySeries software 37 . This matching provides a relationship between the depth (ice core) and date (model). Then, the SE-Dome core age scale based on the isotope composition was constructed by linearly interpolating between the points. To obtain the annual accumulation rate, the depth was resampled at 1-year intervals based on the age scale. The depth of a 1-year interval indicates the annual accumulation rate in snow equivalents. Then, the snow accumulation rate in water equivalent depth was calculated by multiplying snow density. In the same way, the seasonal accumulation rates were obtained based on seasonal boundaries of 1 March, 1 June, 1 September, and 1 December. The varying uncertainty from section to section was calculated by using an algorithm 38 based on a Monte Carlo simulation fitting ensembles of straight lines to subsets of the age data. An age determination error at each tie point should be assigned to calculate a confidence limit. We assumed that the error at each tie point is ±1 month because most of the interannual δ 18 O model peaks, selected for tie points, consisted of three data points. In most periods, the 95% confidence limit was around 1 month (average ± 0.9 months). However, there were larger uncertainties where the number of tie points was small. The largest uncertainty was found in October 2004 (±2.4 months).
MSA analysis and processing in the SE-Dome ice core. We analyzed a 90.45 m depth ice core drilled at the SE-Dome site 39 . For analyses of chemical species in the ice core, we divided the ice core into 100 mm depth sections in a cold room. The samples were decontaminated using a clean ceramic knife in a cold clean room (class 10,000), put into a cleaned polyethylene bottle, and then melted in the bottle at room temperature in a clean room. The methane-sulfonate ion (MS -) was measured by ion chromatography (Thermo Scientific, ICS-2100). We used a Dionex AS-14A column with 23 mM KOH gradient eluent for MS -(we describe MSas MSA). The injection volume of the sample was 1 mL. The number of samples was 395 from 1960 to 2014, corresponding to~7 samples per year on average. Annual and seasonal fluxes of MSA (MSA flux ) were calculated by multiplying an MSA concentration by a water equivalent value of the vertical thickness of the sample for each season (annual: January-December, winter: January-March, spring: April-June, summer: July-September, autumn: October-December).
Uncertainty due to the ice core dating. We evaluated the uncertainty of interannual and seasonal variations in MSA flux resulting from the age determination error in the SE-Dome ice core. The age determination error of each sample was assigned a random number based on a standard normal distribution, which was derived from the 95 % confidence limit of the age of each sample in Furukawa et al. 26 . The age of each sample was re-assigned from the age determination error. We re-caluculated the interannual and seasonal variations in MSA flux using the reassigned age scale. We repeated this processes in 1000 times and obtained 1000 different interannual and seasonal variations in MSA flux to eliminate the bias based on the random number generation. We calculated mean values and standard deviations with these 1000 profiles of MSA flux . As the result, the interannual and seasonal variations in simulated MSA flux were not different from these in original MSA flux even considering the age determination error in the MSA flux (Supplementary Figs. 6 and 7).
Backward trajectory method. To investigate the source area of MSA contained in ice cores at the SE-Dome, we analyzed air mass position along the backward trajectory from the SE-Dome during the past 7 days using the National Oceanographic and Atmospheric Administration (NOAA) Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model 40 and National Centers for Environmental Prediction (NCEP) reanalysis data. The initial positions of the air mass were set at 10, 200, 400, 600, 800, and 1000 m above ground level over the SE-Dome site. The initial date and time were every 6 h from 1960 to 2014. We calculated the probability of the existence of an air mass with a 1°resolution. Considering the MSA supply from the ocean surface, we excluded air masses over 1000 m above ground level. Most MSA is deposited on the snow surface by wet deposition 31 . The existence probability was weighted by the daily amount of precipitation when the air mass arrived at the core site. The daily amount of precipitation was extracted from the ERA5 reanalysis dataset 41 . We calculated the cumulative distribution function (CDF) of the existence probability. We described the sea regions where the CDF was 100-60%, 100-40%, and 100-20% as CDF40, CDF60, and CDF80, respectively.
Satellite-derived chlorophyll-a datasets. The precursor of DMS in the ocean is dimethylsulfoniopropionate (DMSP), which is produced in phytoplankton. We used the merged multiple satellite Chl-a concentration products provided by the European Space Agency GlobColour project (http://www.globcolour.info/) as chlorophyll concentrations for the following discussion. The product is produced by a semianalytical ocean color merging model that estimates Chl-a concentrations, combined dissolved and detrital absorption coefficients, and particle backscattering coefficients from normalized water-leaving radiances retrieved by multiple sensors (SeaWiFS, MERIS, MODIS-Aqua, VIIRS) 42,43 . The temporal coverage and resolution are improved compared to a product of a single sensor 44,45 . The temporal resolution was 8-day and monthly. The spatial resolution was 4 km.
Frequency of intermediate wind and sea ice conditions. The DMS emission is promoted when the wind speed is 8-11 m s −1 9 . The frequency of intermediate wind was defined as the ratio of the total hours when the wind speed at 10 m above sea level was from 8-11 m s −1 to the total hours of the concerned period, such as season or year. We used the hourly product of the wind speed and sea ice concentration from the ERA5 reanalysis dataset 41 . The spatial resolution was 0.25°. Data were discarded if the sea ice concentration in the grid was greater than 10% because the DMS emission from the ocean to the atmosphere was restricted by sea ice covering on the sea surface.
The sea ice retreat day was defined as the late day when the daily sea ice concentration was less than 10% during April-September in a given year. The daily sea ice concentration was calculated by an hourly sea ice concentration product of ERA5 41 .
The sea ice extent was calculated by multiplying the sea ice concentration by the size of the grid cell. The size of the grid cells was calculated according to the Earth ellipsoid.
Reporting summary. Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Code availability
The Fortran90 code was used for data analyses. To read NetCDF datasets, we used a netCDF Fortran (version 4.5.4) and a netCDF C (version 4.8.1) libraries developed by UCAR/Unidata (https://doi.org/10.5065/D6H70CW6). Additional data analyses were done using the Interactive Data Language (IDL) (version 8.7.1) developed by Exelis Visual Information Solutions, Boulder, Colorado. Map figures were generated using the Generic-Mapping Tools (GMT) (version 6.1.0) 47 (https://www.generic-mapping-tools.org/). The statistical tests and graph figures in this study were done using the R software (version 4.1.0) 48 (https://www.r-project.org/).