Holocene millennial-scale variability of coastal environments on the southern coast of Korea and its controlling factors

Abstract

Coastal evolution is influenced by past sea-level changes and resultant shifts from fluvial- to marine-dominant environments and the accompanying significant geochemical and isotopic changes in the water mass and sediments. We investigated the elemental and isotopic features of coastal sedimentary cores (27 m in length) from a small paleo-bay located on the southern coast of Korea to determine such geochemical variability and specify past changes in the bay environment and anoxic conditions and possible links to past climate changes. We analyzed total organic carbon (TOC), total sulfur (TS), their isotopes (δ¹³C_TOC and δ³⁴S_TS), and pyrite. The δ¹³C_TOC values ranging from −25 to −19‰ (a proxy for terrestrial influence) were lower than average (−22.5‰) before 8300 cal yr BP and since 500 cal yr BP, while the intervening Early to Late Holocene showed higher δ¹³C_TOC values, indicating a shallow coastal environment. The δ³⁴S_TS values fluctuating between −35 and +5‰ resembled sedimentation rate change. Based on the changes in the ratios of TOC to TS (C/S ratios), sedimentation rate, and δ³⁴S_TS, we found five possible periods with higher salinity and intensified anoxic conditions at millennial timescales: 8900–8200, 7950–6500, 5200–4300, 3500–2600, and 2000–1100 cal yr BP. These intensified anoxic conditions seem to have been influenced by increased air temperature and sea-surface temperature conditions, which could have intensified the intensity of thermal stratification (less ventilation and mixing) between surface and bottom waters and resultant anoxic conditions.

INTRODUCTION

The most widespread severe anthropogenic influences on estuarine and marine ecosystems are hypoxia and anoxia. The number of dead zones (or hypoxic areas with dissolved oxygen below 2 mg/L) in the coastal oceans has been increasing since the 1960s, mainly due to increases in primary production and consequent coastal eutrophication fueled by riverine runoff of fertilizers and sewage (Diaz and Rosenberg, 2008; Irby et al., 2018). With regard to anticipating global warming, deoxygenation driven by anthropogenic influences is expected to have greater impacts as ocean warming and acidification increase, leading to greater challenges for aquaculture and fisheries from hypoxia. The expected hypoxia and anoxia in coastal areas are mainly associated with a semi-enclosed hydro-geomorphology, hindering water exchange and resulting in water column stratification (Diaz and Rosenberg, 2008; Lee et al., 2018). However, there have been limited studies on the past coastal hypoxia or anoxia variability at centennial to millennial timescales and their controlling factors from the viewpoint of natural climate change and geomorphological evolution, which makes it difficult to constrain future predictions.

Past coastal evolution can be understood in terms of changes in long-term morphology and seawater-influenced bottom-water conditions (hypoxic, anoxic, euxinic, and so forth). Many previous studies have suggested that current coastal morphological evolution is the result of local or regional responses to past climate and sea-level changes since the last glacial maximum and deglacial periods (Stanley and Warne, 1994; Kim and Kennett, 1998; Dellwig et al., 2001; Chen et al., 2004; Hori et al., 2004; Mackie et al., 2005; Lamb et al., 2007; Shin et al., 2007; Nahm et al., 2008; Yang et al., 2008; Yu et al., 2010, 2011; Kim et al., 2012; Ishihara et al., 2012). These studies clearly demonstrated Holocene transgression in various regions in terms of millennial timescales, but there has been little study regarding coastal bottom water and sediment conditions due to the lack of a suitable index and compatible sedimentary records with high-resolution dating ages. Recent geochemical and isotopic data from coastal sediments from the past have shown potential to provide information about such coastal environmental changes (e.g., Pasquier et al., 2017; Lim et al., 2019; Liu et al., 2019).

δ¹³C_TOC values in coastal sediments

Transgression represents a remarkable coastal change during the Holocene, although its timing was asynchronous in each present coastal area due to different geomorphological settings and elevations. This coastal change, including proximity to coasts responsible for the first-order riverine organic input, has been reconstructed using organic carbon isotope records preserved in estuary and coastal sediments (Mackie et al., 2005; Lamb et al., 2006, 2007; Yu et al., 2010, 2011; Williams et al., 2014; Lim et al., 2015). In general, marine aquatic algae and plants have δ¹³C values between −17.3 and −21.7‰, whereas freshwater algae and terrestrial plants (e.g., C₃ plants) have δ¹³C values between −24.9 and −32.5‰ (e.g., Meyers, 1994; Mackie et al., 2005; Lamb et al., 2006, 2007), forming two groups that can be clearly distinguished by their different δ¹³C values. This distinct distribution has been confirmed in different environmental settings. For example, Yu et al. (2010) tested the applicability of δ¹³C values as indicators for sources of organic matter in deltaic and estuarine sediments based on the geochemical information from organic matter from source to sink in the Pearl River catchment, delta, and estuary. They found that the average δ¹³C of estuarine surface sediment increases from −25.0‰ ± 1.3‰ in the freshwater environment to −21.0‰ ± 0.2‰ in the marine environment. Based on these variances in the δ¹³C value, past relative sea-level changes and terrigenous input changes driven by freshwater (e.g., summer monsoons) have been reconstructed in various areas (Mackie et al., 2005; Lamb et al., 2006, 2007; Zong et al., 2006, 2012; Yu et al., 2011; Lim et al., 2015).

C/S ratio and δ³⁴S value in coastal sediments

Pyrite is generally considered a final product of the early diagenesis of sulfur in marine sediments through metabolic processes starting from seawater sulfate (SO₄²⁻) with +20‰ of δ³⁴S value, including microbial sulfate reduction and sulfur isotopic fractionation up to 70‰ (e.g., Canfield et al., 2010; Sim et al., 2011). Pyrite can be formed under the two different processes (diagenetic and syngenetic) (Raiswell and Berner, 1985; Werne et al., 2003; Lyons et al., 2009). In normal marine environments (oxygenated bottom water) the H₂S needed for pyrite formation is produced diagenetically (after deposition) by sulfate-reducing bacteria only below the sediment–water interface. In euxinic environments (anoxic with H₂S in the water column), pyrite forms above the sediment–water interface syngenetically (before burial) as well as in the sediments. It is clear that the pyrite from H₂S is enriched in the light sulfur isotope (³²S), resulting in enrichment of ³⁴S in the parent sulfate (SO₄²⁻). In the euxinic marine setting, the δ³⁴S value is determined by the combination of fractionation through bacterial sulfate reduction (more than −50‰) and reservoir effect, which can suppress net ³⁴S depletion in H₂S and pyrite due to limited sulfate supply in the local reservoir. This reservoir effect will be amplified in the syngenetic process. In euxinic environments, pyrite forms in the water column isolated from the upper oxygenated seawater. As pyrite increases, the heavy ³⁴S isotope in the euxinic water column gradually becomes enriched, resulting in a decrease in the net ³⁴S depletion in H₂S and pyrite due to limited sulfate supply (Lyons et al., 2009).

Compared with the use of carbon isotope analysis for the reconstruction of past sea level–related coastal changes, sulfur isotope analysis of coastal sediments has not been widely used. But several past studies have pointed out the potential for reconstructing alternating marine and nonmarine environments through sedimentary sulfur isotope values. For example, Middelburg (1991) suggested that freshwater sediments underlying the marine sediments can be recognized on the basis of their δ³⁴S values. The Holocene marine sediments in Kau Bay, Indonesia, showed isotopically lighter values (δ³⁴S = −20‰) but underlying freshwater sediments had a δ³⁴S value of +15‰. Similarly, Wilkin and Arthur (2001) reported different δ³⁴S values in each different depositing environment. The δ³⁴S values observed in shallow-water sedimentary cores from the Black Sea covering the past 15 ka ranged widely between −38.0‰ in the stratified anoxic–sulfidic water column and +11‰ in a freshwater lake. The δ³⁴S value in coastal sediments has been reported to vary considerably between inorganic (predominantly pyrite) and organic components. For organic sulfur, δ³⁴S values from −2.0 to +10.0‰ have been reported from the late Pleistocene and Holocene coastal plain sediments in Taiwan (Chen et al., 2004). For pyrite, which is the main form of inorganic sulfur typically formed during bacterial sulfate reduction, much lower values have been reported, varying from +15.2 to −27.6‰ in coastal plain sediments in Taiwan (Chen et al., 2004). Dellwig et al. (2001) investigated marshland sediment in northwest Germany covering the entire Holocene and showed that the peat layers were characterized by significant enrichment of pyrite. The bulk δ³⁴S values of the peat layer sediments ranged from −2.6 to −26.7‰ and were attributed to microbial sulfate reduction following the input of seawater sulfate. Furthermore Chen et al. (2004) found a close relationship between the δ³⁴S of pyrite and the late Pleistocene and Holocene sea-level curve in the coastal plain sediments in Taiwan and suggested that sulfur analysis is a useful technique for reconstructing alternating marine and nonmarine environments.

Furthermore, the ratio of organic carbon to pyrite sulfur (C/S ratio) is thought to be affected by the salinity at the time of deposition, as diagenetic pyrite forms more readily in marine than freshwater sediments due to the relatively high availability of dissolved sulfate in seawater (Berner, 1984; Berner and Raiswell, 1984; Morse and Berner, 1995). For example, C/S ratios were 17–34 at salinities <1‰, whereas at salinities of 19–21‰, C/S ratios were 1.4–1.8 (Berner and Raiswell, 1984 and references therein). Based on the relationship between salinity and C/S ratio, it was suggested that freshwater and marine sediments have high (>10) and low C/S values (0.5–5), respectively, providing two endmembers (Berner, 1984; Berner and Raiswell, 1984; Woolfe et al., 1995, and references therein). Furthermore, it has been suggested that the C/S ratio is 2.8 ± 0.8 in normal marine sediments. In addition, deviations toward higher C/S ratios indicate a freshwater environment, while deviations toward lower C/S ratios suggest burial under a euxinic environment (Berner and Raiswell, 1984).

Recently, several studies have tested the applicability of C/S ratios by measuring the ratios from present river and inner-continental sediments (e.g., Liu et al., 2021, 2022; Chang et al., 2022). For example, Liu et al. (2021) established a new C/S ratio data set from surface sediments of five rivers (Aojiang, Feiyunjiang, Jiaojiang, Oujiang, and Qiantangjiang Rivers), as well as core sediments offshore the Oujiang Estuary. Based on the C/S ratios from the various environmental conditions, they suggested that a C/S ratio of 2.8 can effectively separate freshwater environments (river sediments) from marine environments (Liu et al., 2021). But these studies suggested that C/S ratios should be interpreted with other geochemical indicators, because sedimentation processes, including influence of sea-level change, can modulate the diagenetic path in the marine sediments (Chang et al., 2022; Liu et al., 2022). Thus, the combined information of C/S ratios, the δ³⁴S values, and δ¹³C values in coastal sediments can be used to diagnose past anoxia and test possible controlling factors.

Here, we report and discuss the content and isotopic composition of organic carbon and total sulfur in southern coastal sediments of South Korea. We trace past changes in C/S ratios, the δ³⁴S values, and δ¹³C values in coastal sedimentary sequences and interpret these changes in terms of salinity and anoxic condition changes during the Holocene at multicentennial to millennial timescales.

STUDY AREA

The study area, a small paleo-bay with an area of ~3 km × 3 km, is located where the Seomjin River drains into the northern reaches of Gwangyang Bay, forming an estuarine environment, as shown in Figure 1. There are elongated sandbars in the present river mouth on the northwest of Gwangyang Bay. Subaqueous parts of the sandbars are migrating seaward, forming sandy-silt tidal–subtidal flats on the south inlet of the study area, resulting in a semi-isolated bay environment. The present tide at the Gwangyang Bay and its surrounding coastal areas is semidiurnal with an amplitude of 3 m. This coastal area is characterized by low water exchange between the bay and the outside coastal waters due to subaqueous sills and narrow inlets. The surface sediments in the bay mainly consist of mud and sandy mud, and the benthic environment is partly anoxic, as evidenced by the C/S ratios and hydrogen sulfide content (1–367 ppm) (Hyun et al., 2003).

Figure 1. Study areas (Gwangyang Bay and Goheung Bay) and coring sites on the southern coast of Korea. (A–D) Coring sites: STP17-14 (this study), STP16-20 (Lim et al., 2019), STP17-13 (Park et al., 2019), and STP18-05 (Lim et al., 2022). (E and F) Testing past geomorphological changes between 1920 and 2013. A geomorphological map made in the 1920s showed that study area including a coring site (STP17-14) was a small semi-closed bay characterized by low-lying mountain-fringed coastal lines.

METHODS

Part of the study area has been reclaimed as rice fields. Based on a geomorphological map made in the 1920s showing a small semi-closed bay characterized by low-lying mountain-fringed coastal lines, we selected a coring site (STP17–14) in the central part of the bay. We recovered sedimentary cores 27 m in length in 2017 using a rotary corer, which returns a core sample 1 m in length and 50 mm in diameter in a plastic liner (Fig. 2). After taking pictures, we subsampled the cores at 5 cm intervals. During this process, plant fragments for radiocarbon dating (n = 16) were selected. To test the sedimentation rate (SR) change in core STP17–14 in terms of regional perspective, we collected additional plant fragments (n = 18) from the core STP17–13 located at the mouth of Seomjin River (Park et al., 2019). After pretreatment, including application of the acid–alkali–acid method and graphitization, age dating was performed using the accelerator mass spectrometry (AMS) facility at the Korea Institute of Geoscience and Mineral Resources (KIGAM).

Figure 2. (A) Lithological features, (B) radiocarbon dates and sedimentation rate (SR) change, (C) age–depth model generated from the ¹⁴C ages using the R code package clam (Blaauw, 2010), and (D) photograph of core STP17-14 from Gwangyang Bay on the southern coast of Korea.

To analyze the total organic carbon (TOC) content and its isotope (δ¹³C_TOC), bulk subsamples of ~500 mg were treated with 1 N HCl at ~100°C for 1 h to remove carbonate and then rinsed with distilled water. Samples of approximately 3–5 mg HCl-treated subsamples were loaded into tin combustion cups, and TOC contents were determined using a carbon, nitrogen, and sulfur (CNS) elemental analyzer (vario MICRO Cube, Elementar, Langenselbold, Germany). Simultaneously, stable carbon isotope (δ¹³C) analyses were performed using a continuous-flow isotope ratio mass spectrometer (Isoprime100, G.V. Instruments, Manchester, UK) coupled to the CNS elemental analyzer. The results are expressed in delta (δ) notation relative to the Vienna Pee Dee belemnite standard. The reference material was CH-6 (sucrose, δ¹³C = −10.45‰ ± 0.033‰), which was obtained from the International Atomic Energy Agency (IAEA). Replicated analyses had precision >0.2‰. Analyses of the total sulfur (TS) content and its isotope (δ³⁴S_TS) were performed with bulk subsamples using the same continuous-flow isotope ratio mass spectrometer coupled to the CNS elemental analyzer. The results are expressed in delta (δ) notation relative to the Vienna Canyon Diablo Troilite. The reference material used was NBS127 (seawater sulfate, δ³⁴S = 20.3‰ ± 0.4‰) obtained from the IAEA. Replicated analyses had precision >0.6‰.

The pyrite content of sediments (n = 22) was determined at 3.4, 4.0, 4.8, 5.6, 6.2, 6.8, 8.2, 9.2, 9.8, 10.8, 11.6, 12.2, 13.2, 13.8, 14.2, 14.8, 16.8, 18.2, 19.0, 19.6, 20.2, and 20.4 m, using finely ground samples with a Panalytical X'Pert³ powder diffractometer (Cu Kα radiation; 45 kV; 40 mA). SIROQUANT (v. 3.0, Sietronics, Canberra, Australia) software was used to quantify mineralogy.

RESULTS

Age dating and lithological features

The sedimentary cores from the paleo-bay were formed mainly during the Holocene period based on the ¹⁴C dating results (Tables 1 and 2, Figs. 2 and 3). The sedimentary core STP17-14 can be divided into six units based on lithological features: Unit 1, gray silty sand layer overlaid on the weathered bedrock, located at 26.5–25.8 m; Unit 2, gray coarse-grained sand layer with little structure, located at 25.8–24.0 m; Unit 3, gray silty sand unit with clear sand layers corresponding to ca. 9600–8900 cal yr BP and located at 24.0–19.6 m; Unit 4, green-gray silty mud grains with no clear layers corresponding to ca. 8900–500 cal yr BP and located at 19.6–3.2 m; Unit 5, silty sand sediment, corresponding to ca. 500 cal yr BP to present, located between 3.2 and 1.1 m; and Unit 6, the uppermost layer, consisting of reworked gravel and sandy deposits formed by reclamation.

Table 1. Results of accelerator mass spectrometry (AMS) ¹⁴C dating and calibrated dates for core STP17-14.

^a Calibrated with Radiocarbon Calibration Program (CalPal), http://c14.arch.ox.ac.uk/embed.php?File=oxcal.html.

^b Excluded from SR calculation.

Table 2. Results of accelerator mass spectrometry (AMS) ¹⁴C dating and calibrated dates for core STP17-13.

^a Calibrated with Radiocarbon Calibration Program (CalPal), http://c14.arch.ox.ac.uk/embed.php?File=oxcal.html.

^b From Park et al. (2019).

Figure 3. Radiocarbon dates, high sedimentation rate (SR) events, and photograph of core STP17-13 recovered in the Seomjin River estuary on the southern coast of Korea (modified from Park et al., 2019). The age–depth model for core STP17-13 was generated from the ¹⁴C ages by using the R code package clam (Blaauw, 2010).

This study revealed significant SR changes during the Holocene based on the dating results (n = 16). As shown in Figure 2, in the Early Holocene, SR was not stable and fluctuated with higher SR periods of 9400–8900 cal yr BP and 8300–7500 cal yr BP. In the Middle Holocene, SR showed a decreasing trend, with one SR peak at 5800–5500 cal yr BP. The gradual decrease in SR during the Late Holocene was disturbed by a recent SR increase at 520–150 cal yr BP. The age–depth model for core STP17-14 was generated from the ¹⁴C ages by using the R code package clam (Blaauw, 2010; Fig. 2).

The use of a combination of additional ages (n = 18) dated in this study and previously reported ages in sedimentary core STP17-13 (Park et al., 2019) indicated rapid SR changes during the Early Holocene, including SR events centered at 9800, 9100, 8200, and 7900 cal yr BP (Table 2 and Fig. 3). These SR changes were compared with those of core STP17–14 to test whether they were local to the Seomjin River estuary or at least regional.

Geochemical/isotopic analyses of TOC, TS, and their relationships

No information is available regarding TOC% in Units 1 and 2, because we analyzed samples from Units 3–5 based on age dating results (Fig. 4). In Unit 3, which consisted of a silty sand layer, TOC% showed relatively low values of about 0.3–0.7%. In Unit 4, TOC% increased, with a long-term trend from 0.4% at 19.6 m to 1.6% at 5.7 m, and then decreased. The decreasing trend in the upper part of Unit 4 was continued in Unit 5, reaching 0.15% at 1.4 m. TS% showed a slight increasing trend in Unit 3, and the lower part of Unit 4, reaching 0.55% at 14.4 m. An abrupt increase in TS% was found at 14.2 m with a value of 1.8%, corresponding to the maximum value during the Holocene, and then a gradual decreasing trend was dominant in the rest of Unit 4 and Unit 5. The δ¹³C_TOC values showed rather simple patterns following the features of each unit. In Unit 3, δ¹³C_TOC values were relatively low, fluctuating between −25 and −23‰. At the starting point of Unit 4, the δ¹³C_TOC value increased rapidly to −22‰ and then remained stable, although there were slight fluctuations of less than 1‰, except for one anomaly of −19.8‰ at 5.8 m. Unit 5 showed decreasing δ¹³C_TOC values from −22.3‰ at 3.2 m to −24.5‰ at 1.4 m.

Figure 4. Lithological unit features, radiocarbon dates, and results of geochemical (TOC, total organic carbon; TS, total sulfur), isotopic (δ¹³C_TOC and δ³⁴S_TS), and mineral analyses of sedimentary core STP17-14 from Gwangyang Bay on the southern coast of Korea. The depth profile was divided into five stages according to the geochemical and isotopic features, mainly based on their averages.

As shown in Figure 4, δ³⁴S_TS values were characterized by significant fluctuations. The δ³⁴S_TS values showed a wide range of variation between −32 and 5‰. In Unit 3, there were three remarkable peaks with δ³⁴S_TS values of ca. −11 to −8‰. Unit 4 showed clustering of higher δ³⁴S_TS values up to 6‰ located at depths of 14.2 to 13.2 m. Then, a gradual decreasing trend was observed in order of increasing depth, reaching a minimum δ³⁴S_TS value of −32‰. In Unit 5, δ³⁴S_TS values increased to −7‰. C/S ratios in the depth profile changed within a narrow range of 0.3–4.8, showing relatively low values in Unit 4 compared with Units 3 and 5. These elemental and isotopic data showed a long-term change consisting of five stages that reveal geochemical features of the deviation from the averages of the sedimentary δ¹³C_TOC values, TS%, δ³⁴S_TS values, and C/S ratios (Fig. 4).

Pyrite% was measured at irregular intervals to check its content in the sediments of low or high δ³⁴S_TS values. For example, the sediments at 3.4 m with the minimum δ³⁴S_TS value of −32‰ had a pyrite content of 3.4%, and at 14.2 m with the maximum δ³⁴S_TS value of 5.2‰, the pyrite content was 3.1%, with values of 1.8–3.8% in Unit 3 and the upper part of Unit 4.

DISCUSSION

Long-term coastal response to Holocene sea-level change

As indicated in many previous studies, it is clear there were significant sea-level changes of up to 20 m during the Holocene (Tanigawa et al., 2013; Lambeck et al., 2014). To specify past long-term coastal environmental changes, it is essential to constrain water depth in relation to the sea level at that time. There are few reconstruction records of Holocene sea-level changes on the southern coast of Korea, including the present study area. We estimated past water-depth changes in the coring site (Fig. 5) based on the reconstructed sea-level changes from the Yellow Sea and East China Sea (Lambeck et al., 2014 and references therein). Based on these estimated water depths, the present tidal range (~3 m), and geochemical information, reconstructions of past coastal environmental changes in the paleo-bay located in the northern part of Gwangyang Bay were created and tested in terms of spatial coastal changes by comparison with other coastal changes on the southern coast of Korea.

Figure 5. Changes in total organic carbon (TOC), total sulfur (TS), C/S ratios, their isotopes (δ¹³C_TOC and δ³⁴S_TS), pyrite, and sedimentation rate (SR) of core STP17-14 from Gwangyang Bay on the southern coast of Korea during the Holocene. Past water depth in the coring site was estimated based on the reconstructed past sea levels from the Yellow Sea and East China Sea (Lambeck et al., 2014 and references therein). The Holocene environmental change in the study area was divided into five stages according to the geochemical and isotopic features.

Water depth at the coring site during the Holocene was controlled by both sea-level change and SR (or sediment accumulation rate) (Bird et al., 2007). Water depth can be approximated by subtracting the elevation of the sediment–water interface (or seafloor) in the core from the elevation of reconstructed sea-level curve (Lambeck et al., 2014 and references therein). As shown in Figure 5, the reconstructed sea level at ~9600 cal yr BP was −18 m, and the sediment–water interface elevation in the core site was −24 m, suggesting a water depth of ~6 m at that time. Considering the present tidal range with an amplitude of ~3 m in Gwangyang Bay, the coring site seems to have been included in a subtidal environment, indicating that the coring site experienced transgression before 9600 cal yr BP. This water depth increased according to the subsequent sea-level rise and reached up to 17 m at 8200 cal yr BP, suggesting that the coring site became a small semi-closed bay. After that, it is likely that the water depth decreased, and the study site then had an intertidal to subtidal environment at 2200 cal yr BP. This long-term coastal evolution traced by simple comparison with the reconstructed sea-level changes can be supported by geochemical indicators, including sedimentary δ¹³C_TOC values, δ³⁴S values, C/S ratios, and SRs.

The present surficial δ¹³C_TOC values along the river and coastal areas in China and Korea show the endmember to trace past coastal environments (e.g., Zhan et al., 2011; Williams et al., 2014). For example, Zhan et al. (2011) reported different δ¹³C values according to modern sedimentary environments from the lower Yangtze River to the East China Sea. In the river, the δ¹³C values from suspended particulate matter and sediments range from −28 to −24.4‰, and values in the mid- to lower river mouth decrease to −24.3 and −20.5‰. Samples from shallow-marine deposits show values between −22.7 and −20‰. A similar pattern was reported from modern surficial sedimentary δ¹³C values along the Yeongsan Estuary, Korea (Williams et al., 2014). The δ¹³C of modern lake surficial sediments located in the upper part of the Yeongsan Estuary varied from −26 to −24‰. Inner estuary sediments ranged between −23.5 and −18‰, whereas outer estuary and coastal sediments had the highest values, varying between −18 and −16.5‰ (Williams et al., 2014), suggesting clear spatial differences in the δ¹³C values.

Sedimentary δ¹³C_TOC values in this study varied between −25 and −21.5‰ with an average of −22.5‰ (except for a sample with −19.0‰) and showed a long-term change consisting of five stages that reveal geochemical features of the deviation from the average, for example, the δ¹³C_TOC values, TS%, δ³⁴S_TS values, and C/S ratios (Fig. 5). Based on the present surficial δ¹³C_TOC values along the river and coastal areas in China and Korea (e.g., Zhan et al., 2011; Williams et al., 2014), during the period between 9600 and 8400 cal yr BP corresponding to Stages 1 and 2, the sedimentary δ¹³C_TOC values fluctuated between −25 and −23‰, suggesting a significant terrestrial organic input. In Stage 1, there were significant changes in SR and δ³⁴S_TS values. Stage 1, which covers the period 9600–8900 cal yr BP, showed a relative increase in SR up to 7.1 mm/yr and enriched peaks of δ³⁴S_TS values at 8900, 9200, and 9500 cal yr BP. In Stage 2, which corresponds to the period 8900–8400 cal yr BP, both SR and δ³⁴S_TS values decreased. This stage was followed by a clear shift to a bay environment with a water depth of 16 m in response to the Early Holocene sea-level rise, as shown by the increase in δ¹³C_TOC value from −24 to −22‰. This bay environment of Stage 3 in the period 8400–7500 cal yr BP was characterized by a remarkable increase in SR up to 20.4 mm/yr and enrichment of δ³⁴S_TS up to +5‰ coupled to the minimum C/S ratio (0.4) and maximum pyrite content (3.8%), suggesting the frequent occurrence of significant anoxic conditions. In Stage 4, the stable δ¹³C values during the Middle to Late Holocene indicated relatively stable coastal environments in the study area. During this stage, SR and δ³⁴S_TS values decreased gradually, reaching 0.7 mm/yr and −30‰, respectively. Since 1500 cal yr BP (late Stage 4), the shallow coastal environments shifted to tidal-influenced environments, as shown by the deceased δ¹³C_TOC value of −24‰ and water depth of ~3 m. In the relatively recent period from 500 cal yr BP to the present (Stage 5), SR increased to 5.4 mm/yr and δ³⁴S_TS increased to +10‰.

Changes in the δ³⁴S values and sedimentation rate during the Early Holocene

The sedimentary isotope records in this study clearly indicated changes in the δ³⁴S_TS values and related coastal environmental changes in a semi-isolated bay during the Early Holocene. As pointed out by Chen et al. (2004), the δ³⁴S value in coastal sediments shows different ranges between inorganic (predominantly pyrite) and organic components. For example, organic δ³⁴S values vary from −2.0 to +10.0‰ and inorganic sulfur values vary from +15.2 to −27.6‰ (Chen et al., 2004). Compared with these values, the δ³⁴S_TS values in core STP17-14 showing much depleted variation between −32 and +5‰ may have been affected mainly by inorganic components, although the influence of organic sulfur cannot be excluded. This can be supported by the high pyrite content (2–3.5%) in Figure 5, suggesting that TS mainly consisted of pyrite. Thus, the δ³⁴S_TS values may provide information mainly about pyrite formation in the coastal sediments and this assumption will be tested in the following discussion.

Isotopic fractionation during the metabolic processes of pyrite formation can be influenced by a variety of biochemical and environmental conditions, including limitation of sulfate, organic matter, iron minerals, and depositional environment (e.g., SR) (e.g., Berner, 1984; Dellwig et al., 2001; Werne et al., 2003; Pasquier et al., 2017; Liu et al., 2019). The δ³⁴S values in the sedimentary cores in this study were mainly formed after transgression, indicating that there was sufficient sulfate in the paleo-bay provided by seawater. As shown in Figure 5, there was no significant covariance between TOC content and δ³⁴S value, suggesting that organic matter input was not the main factor regulating the δ³⁴S values in this study. The coring site in this study is located near the river mouth, and iron was available during the process of pyrite formation.

Recently, a possible linkage between the δ³⁴S values of pyrite and SR change was reported from a study using a 300 m drill core of Mediterranean sediments deposited over the past 500,000 yr (Pasquier et al., 2017). This study revealed that the fluctuations of the δ³⁴S values of pyrite in glacial–interglacial timescales have been affected by sea-level and temperature changes controlling SR change in the study area. Pasquier et al. (2017) suggested that the change in the δ³⁴S values of pyrite could be influenced by local depositional conditions, especially changes in SR that control connectivity with the overlying water column. Related to a possible millennial-timescale linkage is the report by Liu et al. (2019), who suggested that the δ³⁴S values of pyrite ranging between −38.2 and +15‰ during the last 12,500 cal yr BP were influenced by SR changes based on geochemical and isotopic studies using the 60 m drilled core sediments from the inner shelf of the East China Sea. This study clearly showed the δ³⁴S values have been controlled by SR change in the inner shelf linked to the winter monsoon–driven coastal current intensity. For example, the δ³⁴S values are markedly depleted during the intervals with lower SR (δ³⁴S values of pyrite < −30‰) relative to that observed during the intervals with higher SR (δ³⁴S values of pyrite >0‰). These results were attributed to limited connectivity between porewaters and overlying seawater, in which biological fractionation is oppressed, resulting in enriched δ³⁴S values.

As shown in Figures 5 and 6, the SR change in the study area was closely correlated with the long-term sedimentary δ³⁴S values during the Early Holocene. The remarkable variability superimposed on the long-term change indicated fluctuation of the δ³⁴S_TS values on multi-centennial timescales (indicated by the letters c–e in Fig. 6). These occurrences of higher δ³⁴S_TS values were consistent with the higher SR and increased pyrite deposition rate of 3.8%. To examine whether the observed synchronous changes in SR and δ³⁴S_TS values were local or regional, we performed additional high-resolution age dating using sedimentary cores recovered from the Seomjin River estuary (Park et al., 2019) and calculated SR changes in the sedimentary core STP17–13 (Fig. 3). As shown in Figure 6, SR in core STP17–13 was similar to the fluctuations of SR and δ³⁴S_TS values of core STP17–14. For example, the enriched δ³⁴S_TS events at 8900, 8100, and 7900 cal yr BP in core STP17–14 corresponded to increased SR in core STP17–13 within a dating error of ~100 yr. Similarly high SR has been found in a river mouth located on the southern coast of Korea at around 8000 cal yr BP and attributed to seawater spilling onto the site as a result of sea-level rise (Lim et al., 2015). This suggests that the SR changes responsible for the enriched δ³⁴S_TS events in core STP17–14 were not driven by local geomorphological features in the paleo-bay. The similarity between sedimentary δ³⁴S values and SR changes identified in this study from the southern coast of Korea supports the recent results showing SR-controlled δ³⁴S values in ocean-bottom sediments, resulting in open or closed diagenetic processes (Pasquier et al., 2017; Liu et al., 2019).

Figure 6. Comparison of total sulfur (TS), C/S ratios, δ¹³C_TOC and δ³⁴S_TS, and sedimentation rate (SR) of core STP17-14 with SR change in core STP17-13 on the southern coast of Korea during the Early to Middle Holocene. Past water depth in the coring site was estimated based on the reconstructed past sea levels from the Yellow Sea and East China Sea (Lambeck et al., 2014 and references therein).

Furthermore, significantly increased (more than ~10 times) SR could have supplied increased organic matter to the water column, encouraging oxygen consumption, although relative organic matter content (TOC%) was not high. Under this marine setting, the δ³⁴S value may have been influenced by the combination of fractionation through bacterial sulfate reduction (> −50‰) and a reservoir effect from the water column, which can suppress net ³⁴S depletion in H₂S and pyrite due to limited sulfate supply in the local reservoir. Given that the study site is characterized by its semi-enclosed geomorphology and limited water circulation resulting in strong seasonal stratification (Lee et al., 2018), this reservoir effect coupled with significantly increased SR could have been amplified in the syngenetic process, although it is difficult to specify relative contribution of pyrites formed syngenetically (before burial) from dominant ones formed diagenetically (after burial) in this study.

Centennial- to millennial-timescale variability of salinity and anoxic conditions during the Holocene

As shown in Figures 5 and 6, the time series of the δ³⁴S_TS values and C/S ratios revealed significant fluctuations during the Holocene, suggesting possible centennial to millennial variability in salinity and anoxic conditions.

To examine whether these fluctuations occurred locally or regionally, we compared the time series of C/S ratios in this study with those from Goheung Bay on the southern coast of Korea (Lim et al., 2019) and from the previously seawater-filled Daesan Basin (paleo-Daesan Bay) located in the middle reach of the present Nakdong River, Korea (Lim et al., 2022). Based on the relationship between salinity and C/S ratio (Berner, 1984; Berner and Raiswell, 1984; Woolfe et al., 1995 and references therein) and differences from present river and inner-continental sediments (e.g., Liu et al., 2021, 2022; Chang et al., 2022), the C/S ratios at the three sites (cores STP16-20 [Lim et al., 2019], STP17-14 [this study], and STP18-06 [Lim et al., 2022]) revealed frequent coastal environmental changes in terms of anoxic conditions (Fig. 7). Given that core STP18-06 is located in the middle reach of Nakdong River, the similar variability found among these cores suggests simultaneous responses to seawater-related environmental changes at that time. Based on the average values, these fluctuations can be roughly divided into five periods with reduced C/S ratios, suggesting higher salinity and more intensified anoxic conditions: 8900–8200, 7950–6500, 5200–4300, 3500–2600, and 2000–1100 cal yr BP.

Figure 7. (A–C) Comparison of changes in the C/S ratio of core STP17-14 (this study) with those of core STP16-20 from Goheung Bay on the southern coast of Korea (Lim et al., 2019) and core STP18-06 from Daesan Plain (paleo-Daesan Bay), located at the middle reach of the Nakdong River, Korea (Lim et al., 2022). See Fig. 1 for details. (D) Air temperature during the Holocene reconstructed by using argon and nitrogen isotopes within trapped air in a Greenland ice core (Kobashi et al., 2017). (E) Holocene Mg/Ca-derived sea-surface temperature (SST) from foraminifers in the sediment cores in the western tropical Pacific Ocean (Stott et al., 2004). (F) Asian summer monsoon index (Dykoski et al., 2005).

At present, the hypoxic event on the southern coast of Korea is a chronic problem that occurs every summer (Lim et al., 2006). For example, Jinhae Bay, where physical energy is low due to its semi-enclosed geomorphology and limited water circulation, experiences strong seasonal stratification (Lee et al., 2018). During March to May, thermal stratification is formed due to solar heating, decreasing bottom-water dissolved oxygen. These hypoxic conditions are greatly intensified by thermal stratification due to summer solar heating. Furthermore, this water column stability is intermittently strengthened by freshwater discharges. From the end of summer in late August, the stratification is destroyed by active vertical mixing between surface and bottom waters due to thermal cooling and elevated wind forcing (Lee et al., 2018). Past anoxic conditions in the study area could have been influenced mainly by temperature, and the factors responsible for controlling possible anoxic conditions in the southern coastal areas in Korea may have been limited water circulation and thermal stratification linked to temperature.

To examine the possible influence of temperature on past anoxic condition changes in the study area, we performed comparisons with air temperatures reconstructed by using argon and nitrogen isotopes within trapped air in a Greenland ice core (Kobashi et al., 2017). Interestingly, the fluctuation of the C/S ratios was quite similar to the temperature changes in Greenland Summit during the Holocene in terms of centennial to millennial timescales (Fig. 7E). As shown in Figure 7, the periods of higher temperature in Greenland correspond to decreased C/S ratios. For example, the proposed periods with more intensified anoxic conditions seem to have occurred during the synchronous high-temperature periods (7950–6500, 5200–4300, 3400–2600, and 2000–1000 cal yr BP), suggesting possible linkage between them. Kobashi et al. (2017) reported that the reconstructed temperature of Greenland was correlated with significant atmospheric circulation patterns over Greenland and monsoon activity in Oman and China, indicating that the Greenland temperature changes were hemispheric signals. Based on this suggestion, the similarity observed in this study suggests that the fluctuations in C/S ratios may have been influenced by global temperature and atmospheric circulation changes. At present, the hypoxic conditions on the southern coast of Korea are mainly controlled by the intensity and duration of thermal stratification due to solar heating according to seasonal change. It is likely that the higher air temperature conditions linked to higher temperatures in Greenland could have intensified the thermal stratification (less ventilation and mixing) and resultant anoxic conditions, and the higher temperatures could have extended the anoxic conditions in the following autumn by delaying vertical mixing due to higher temperatures and reduced wind forcing.

With regard to the possible influence of sea-surface temperature (SST) on coastal anoxic conditions, the Mg/Ca-derived SST from foraminifers in the sediment cores in the western tropical Pacific Ocean (Stott et al., 2004; Fig. 7F) seem to have affected the anoxic conditions to some extent. For example, the significantly increased anoxic conditions at 8800–8200, 7950–6500, 5200–4300, 4000–3600, 3400–2600, and 2000–1000 cal yr BP corresponded to higher SST periods partly associated with millennial-timescale fluctuations (Fig. 7). Isono et al. (2009) reported clear millennial-timescale periodicity (e.g., 1470 yr) in SST changes in the northwestern Pacific off central Japan during the Holocene. They suggested a climatic link between the North Pacific gyre system responsible for changes in the SST and the pathway of the Kuroshio Current and the high-latitude North Atlantic thermohaline circulation. The similarity in fluctuations between C/S ratios and SST in the western tropical Pacific suggest that the millennial-timescale changes in anoxic conditions in coastal areas of East Asia may have been linked to SST in the western tropical Pacific through atmospheric–oceanic circulation changes. In addition, the increased SST signals shown in Figure 7 may have been transferred to the southern coast of Korea by the Kuroshio Current and Tsushima Current, which are responsible for transporting heat energy in the western Pacific Ocean, influencing the chances of thermal stratification between the surface and bottom waters in the study area.

The C/S ratios in the bay in this study could have been affected by freshwater input mainly driven by the East Asian summer monsoon. Lim et al. (2015) suggested that sedimentary C/S ratios in the river mouth and surrounding areas could have been controlled by past freshwater input events during the Holocene. To test this, we compared the C/S ratios with the change in the Asian summer monsoon intensity reconstructed from the δ¹⁸O values of a speleothem in Dongge Cave, South China (Dykoski et al., 2005; Fig. 7F). The two data sets are quite different in terms of the long-term change and their millennial-timescale variability. This suggests that the influence of summer monsoonal rainfall on the C/S ratios in the study area (a semi-closed bay environment) was very weak and supports possible strong influences of atmospheric temperature and SST in the western tropical Pacific on the past coastal hypoxic and anoxic events.

CONCLUSIONS

To test possible change in coastal anoxic conditions driven by natural climate change, we reconstructed past coastal evolution and anoxic conditions of bottom water based on TOC%, TS%, C/S ratios, and isotopic values of δ¹³C_TOC and δ³⁴S_TS of coastal sedimentary cores. We identified five possible periods with intensified anoxic condition variability at millennial timescales characterized by enriched δ³⁴S_TS values, increased SR, and decreased C/S ratio: 8900–8200, 7950–6500, 5200–4300, 3500–2600, and 2000–1100 cal yr BP. These intensified anoxic conditions seem to have been influenced mainly by global air temperature and SST conditions at these times, which could have intensified thermal stratification and resultant anoxic conditions. Our results demonstrate the possible utility of δ¹³C_TOC, δ³⁴S_TS, C/S ratios, and high-resolution SR records for tracing past coastal environment changes, including anoxic conditions. With regard to anticipation of global warming and coastal responses, temperature increases in the atmosphere and ocean will exacerbate the bottom-water conditions and increase hypoxic conditions, causing frequent anoxic conditions. Further studies in various areas using proxies for coastal environments (e.g., δ¹³C_TOC, δ³⁴S_TS, and C/S ratio) are required in terms of regional- and global-scale responses to climate change.