Abstract
The Gulf of Thailand is ideal for studying eustatic sea level fluctuations in Southeast Asia due to its shallow basin and tectonic stability. However, our understanding of how this region’s relative sea level (RSL) has fluctuated over the Holocene epoch is far from complete. In this study, we used lithostratigraphy, loss on ignition, grain size, and pollen analyses to reconstruct the environmental changes in the Sam Roi Yot wetland, which was significantly influenced by seawater intrusion, driven by fluctuations in RSL in the Gulf of Thailand. Therefore, the analyzed pollen records of the sediment core from the wetland reflected variabilities in the RSL in the Gulf of Thailand. Subsequently, we found that after a sea level highstand prior to 4000 cal y BP, the RSL gradually fell with two significant regressions at c. 2950 and 1850–1450 cal y BP before rising at 1450–1050 cal y BP and declining after that. The inconsistency between RSL reconstruction based on our results and the global sea level changes simulated by the Glacial Isostatic Adjustment (GIA) model further suggests that long–term El Niño Southern Oscillations (ENSO) variabilities may have played a significant role in sea level changes in the Gulf of Thailand over the Late-Holocene period. Thus, during extended El Niño or La Niña conditions, the sea level would have been consistently lower or higher than expected from eustatic and isostatic processes alone. Overall, this study emphasizes the importance of considering regional factors such as ENSO to understand sea level changes in Southeast Asia.
Introduction
Rising sea levels due to climate change have been intensely studied (Bindoff et al., 2007; Cazenave et al., 2014; Dangendorf et al., 2019; Mimura, 2013; Rahmstorf, 2010). The melting of the Earth’s glaciers has led to global isostatic adjustments leading to either a rise or subsidence of landmass (Whitehouse, 2018; Wu et al., 2010). There is also the human-induced subsidence that has resulted in the disproportionate increase, three times the global average rate, in sea level rise along the coastal regions, particularly in Southeast (SE) Asia (Nicholls, 2021; Nicholls et al., 2021). This is attributed to increasing groundwater extraction due to sprawling megacities along the coast, leading to compaction of sediments within the aquifers and thereby causing the above land to sink (Nicholls, 2021; Nicholls et al., 2021). Despite our current understanding, how the relative sea level (RSL) has fluctuated in SE Asia, particularly over the Late-Holocene, remains far from understood. The Gulf of Thailand, a shallow sea bordering the eastern part of the Southeast Asian Mainland, is an ideal location to study RSL fluctuations (Figure 1). The Gulf of Thailand is considered a far–field site and tectonically stable over the Holocene epoch and will thus not require glacial isostatic and tectonic adjustments (e.g. Chabangborn and Wohlfarth, 2014; Chabangborn et al., 2020; Hanebuth et al., 2011; Horton et al., 2005; Oliver and Terry, 2019).

Location of the study area on (a) the western coast of the Gulf of Thailand with (b) the study in the Sam Roi Yot wetland in a red square and the sites included in discussions (further information on each site is in Supplementary 5, available online). The wetland covers the green area of approximately 37 km2 (c) comprised of (d) densely emergent plants enclosing the open water bodies. (e) Location of the coring points in the Sam Roi Yot wetland. Blue, green, and red lines are the correlation for the eastern, western, and NW-SE transects, respectively. Core CP4 is a red circle.
Several studies from this region have reported a general RSL rise during the early Holocene and reached its highstand at the Mid-Holocene (e.g. Bird et al., 2007, 2010; Chabangborn et al., 2020; Hanebuth et al., 2011; Oliver and Terry, 2019; Scoffin and Tissier, 1998). However, the RSL changes during the Late-Holocene are under debate. For example, Sinsakul et al. (1985), and Sinsakul (1992) suggested a gradual RSL fall, punctuated by two significant transgressions at about c. 3500 and 2000 cal y BP (Supplementary 1, available online). Similarly, Tjia (1992), based on marine shells and peat, identified c. 6000, c. 4000, and c. 2700 cal y BP as Mid-/Late-Holocene highstands and attributed it to periods of short regression and transgression (Supplementary 1, available online). Others, on the contrary, have argued for a gradual decline in RSL after the Late-Holocene highstand, although the timing differs (e.g. Choowong et al., 2004; Chua et al., 2021; Horton et al., 2005; Hutangkura, 2014; Nimnate et al., 2015; Scoffin and Tissier, 1998; Surakiatchai et al., 2018) (Supplementary 1, available online). Scoffin and Tissier (1998), based on the study of intertidal reef–flat corals (microatolls), identified a single Holocene highstand with a constant rate of RSL fall from c. 6000 cal y BP. On the other hand, Horton et al. (2005) identified the Mid-Holocene highstand of ∼5 m above present sea level (APSL) from ∼4900 to 4500 cal y BP, based on pollen records (Supplementary 1, available online). These discrepancies could be attributed to the lack of understanding of the competing factors influencing the RSL changes in the region, resulting in the differences in the timing and magnitude of RSL fluctuations reported (e.g. Chabangborn et al., 2020; Choowong et al., 2009; Jiwarungrueangkul et al., 2022; Oliver and Terry, 2019; Sainakum et al., 2021; Surakiatchai et al., 2018).
Considering the societal impact of global sea level rise, understanding how the RSL has fluctuated over the Holocene epoch is necessary to predict the risks of rising sea levels on the coasts of SE Asia. This study, therefore, seeks to reconstruct the environmental history of Sam Roi Yot wetland, located in the Khao Sam Roi Yot National Park on the west coast of the Gulf of Thailand (Figure 1). Since this area is directly affected by seawater intrusion, the sedimentary sequence from the Sam Roi Yot wetland has been suggested to record the RSL changes in the Gulf of Thailand (Koskelainen, 2014; Sritipsak, 2014; Surakiatchai et al., 2018). The environmental reconstruction was based on lithostratigraphy, loss on ignition (LOI), grain size, and pollen analyses and constrained by nine radiocarbon dating results. Our study does not only contribute to the understanding of the temporal variability of the Late-Holocene RSL changes in the Gulf of Thailand. Still, it could also be analogous to the predicted sea level rise driven by global warming.
Regional setting
The Sam Roi Yot wetland (12.24°N, 99.93°E) is a freshwater marsh of about 1–1.5 m water depth (Koskelainen, 2014) located in the Khao Sam Roi Yot National Park (12°04′–12°21′N and 99°51′–100°02′E), Prachuap Khiri Khan Province, on the west coast of the Gulf of Thailand (Figure 1b). The mean sea level is 1.61 m above the tide datum of the mean lowest low water, and the mean tidal range is 1.03 m (Williams et al., 2016). The park covers an area of approximately 98.8 km2 and can be divided into wetland, limestone mountain range, and coastal plain from west to east based on their geomorphology (Koskelainen, 2014; Surakiatchai et al., 2018).
The freshwater wetland is approximately 37 km2, located in the western part of the national park (Figure 1c), and is covered by water all year round. Precipitation and the small streams from the limestone mountain ranges are the central inflow feeding into the wetland (Sritipsak, 2014). The west canals originate from the Tenasserim Range and partly supply water to the wetland (Niyomdecha, 2019). In addition, it is occasionally affected by seawater through a tidal creek, called Klong Khao Daeng, which connects the Gulf of Thailand to the Sam Roi Yot wetland (Parr et al., 1993; Sritipsak, 2014; Surakiatchai et al., 2018; Figure 1b). The wetland is densely vegetated by emergent plants, for example, Phragmites karka, Eleocharis dulcis, Arundo donax, and Typha angustifolia (Sritipsak, 2014), while submerged plants, for example, Utricularia aurea, Najas graminea, and Hydrilla verticillata, can be found in the open water bodies that are scattered in the wetland (Sritipsak, 2014) (Figure 1d).
The center of the park has a series of limestone hills along a north–south direction (Figure 1c and d), which are mainly composed of gray to bluish–gray limestone interbedded by brown feldspathic and calcareous sandstone that were deposited during the Permian (Surakiatchai et al., 2018). The thin soil and the barren rocks allow only dwarf evergreen and deciduous trees and shrubs, for example, Dracaena cochinchinensis, Opuntia cochenillifera, and Wrightia lanceolata, to thrive on the mountain (Koskelainen, 2014; Niyomdecha, 2019).
The coastal plain consists of beach ridge plains in the southeastern and northeastern parts of the national park and a tidal flat in the eastern part (Surakiatchai et al., 2018, 2019). The southeast beach ridge plains comprise a series of NE–SW ridges and swales, which arrange in the NW–SE direction in the northeast plain. The tidal flat is dissected by a tidal creek network (Surakiatchai et al., 2018) (Figure 1c). Mangrove plants, for example, Rhizophora apiculata, R. mucronata, Lumnitzera littorea, and Suaeda maritima, can be found in the tidal–flat area (Niyomdecha, 2019).
Since Prachuap Khiri Khan Province is in the rain shadow region caused by the obstruction of the southwest monsoon by the Tenasserim Range, the mean annual precipitation is 1112 mm, slightly lower than the mean annual precipitation in Thailand (Amatayakul and Chomtha, 2017). However, the increase in rainfall in the early winter season (October–November) corresponds with other areas along the west coast of the Gulf of Thailand and is explained by humidity transportation from the Gulf of Thailand by the northeast monsoon (Amatayakul and Chomtha, 2017; Chabangborn et al., 2020). The mean annual temperature is 27.1°C (Amatayakul and Chomtha, 2017).
Methodology
Sediment cores were taken from nine points in the Sam Roi Yot wetland in July 2018 on a coring platform using a modified Russian corer (7.5 cm diameter, 1 m length) (Figure 1e). The sediment cores were collected at a 50 cm overlap depth with each coring point to archive a continuous sequence. Since core CP4 was taken from the deepest water depth (1.8 m), it was used for further analyzes of the loss on ignition (LOI), pollen composition, grain size distribution, and radiocarbon dating. At the other core sites, the water depth was 1.7 m at sites CP1, CP7, and CP8, 1.5 m at sites CP5 and CP9, 1.4 m at site CP6, 1.3 m at site CP2, and 1.2 m at site CP3 (Figure 2a–c). The shallow water north of site CP8 and dense vegetation in the west and south of CP1/CP2 made it difficult to move further by the inflatable boat (Figure 1e). Each sediment core was carefully wrapped in plastic, placed in PVC tubes, and transported to the Department of Geology, Chulalongkorn University, for further analysis.

Lithostratigraphic correlations for the (a) east, (b) west, and (c) NW-SE transects (the locations of coring points are presented in Figure 1e).
The sediment cores were described following Wohlfarth (2014) and correlated between the overlapping cores to construct the composited depth and lithostratigraphic column for individual coring points. The lithostratigraphic columns were correlated in two parallel N–S and one NW–SE transects (Figure 1e). The N–S transects were in the eastern (i.e. cores CP1, CP4, CP7, CP9, and CP8) and western (i.e. cores CP2, CP3, BW2011 (which was taken by the Asian monsoon project at the Stockholm University), and CP6) parts of the wetland (Figures 1e, 2a and b). The NW–SE transect included cores CP4, CP5, BW2011, and CP6 (Figures 1e and 2c).
The samples for 14C dating obtained from core CP4 were prepared following Björck and Wohlfarth (2001). First, each sample was sieved under running tap water using a mesh size of 0.5 mm diameter. The sieved remains were cleaned several times with deionized water and carefully identified under a binocular light microscope. Next, the plant remains were picked out, for example, leaves, charcoal, small twigs, and wood fragments. The selected samples were cleaned multiple times with deionized water and dried at 60°C for 12 h. Next, the dried samples were pretreated following the acid-base-acid method (de Vries and Barendsen, 1952) using (a) 1 M hydrochloric acid (HCl) and (b) 1 M sodium hydroxide to remove (a) carbonate and (b) fluvic and humic acids, respectively. The pretreated samples were cleaned with deionized water and dried at 60°C for 12 h. Finally, the samples were sent to the DirectAMS for radiocarbon analysis using an accelerator mass spectrometry (AMS). The 14C age ± 1 standard deviation (SD) was calculated following the conventions of Stuiver and Polach (1977) using a Libby half-life of 5568 years and a fractionation correction based on δ13C measured on the AMS.
The radiocarbon dating results were calibrated with IntCal20 (Reimer et al., 2020). The age–depth model was constructed using Bacon, a Bayesian statistics-based routine, from 3.645 to 1.80 m depth below the water surface (DBWS) (Blaauw and Christen, 2011). The lowermost part of core CP4 (layer 1: 3.645–3.905 m DBWS) found several angular to sub-angular gravels (approximately 5%) which had conceivably collapsed from the limestone mountain. Consequently, it was excluded from further analysis hereafter.
For the LOI analysis, consecutive 2 cm samples were dried at 105°C for 12 h and then combusted at 550°C for 6 h and 950°C for 3 h to assess the organic matter (OM) and carbonate contents, respectively (Heiri et al., 2001; Wohlfarth et al., 2012). The LOI was calculated as the percentage weight loss of the dried sample.
Grain size distribution was analyzed for each consecutive 4 cm depth. The samples were prepared following Rowell (1994). In brief, each sample was treated with 20 mL of 10% (v/v) HCl and 20–10 mL of 30% (v/v) H2O2 to remove the carbonate and OM. To disperse the particles, 5% (w/v) sodium hexametaphosphate (Calgon) was added to the prepared sample. The samples were then sent to the Scientific and Technological Research Equipment Center, Chulalongkorn University, for the grain size distribution analysis by laser diffraction using a Malvern Mastersizer 3000. The results were further analyzed using the GRADISTAT program (Blott and Pye, 2001).
Subsamples (about 1 cm3) for pollen analysis were taken at 5 cm intervals. The subsamples were sieved with a mesh size of 0.2 mm, and the sieved remains were treated by the acetolysis method (Ellison, 2008). Pollen grains were extracted using the heavy liquid approach with sodium polytungstate (3Na2WO4•9WO3•H2O) (Englong et al., 2019). Afterward, an aliquot of each treated subsample was dried onto a coverslip and then mounted onto a glass slide. Pollen grains from each depth interval were identified by comparing to modern pollen references (Punwong, 2007), and counted under light microscopy at the Faculty of Environmental and Resource Studies, Mahidol University, and the Department of Geology, Faculty of Science, Chulalongkorn University. In these samples, pollen grains were enumerated to achieve a count of 300 pollen grains, but the pollen content of some samples was very sparse and insufficient to achieve this target (i.e. less than 300 pollen grains/cm3). Pollen taxa were categorized into ecological assemblages: mangroves, back mangroves, non-mangrove taxa, and unknown. Mangroves and back mangroves were grouped according to the distribution of trees and shrubs in the mangrove formations in Thailand by Santisuk (1983).
Principal component analysis (PCA) was performed using SPSS version 28 to gain insight into the relationship between LOI, sediment components, and palynological compositions.
Results
Lithostratigraphic correlations from SRY wetland
The sedimentary sequences were divided into six units, from bottom to top: units A, B, C, D1, D2, and E, based on the physical properties of sediments and dating results from core BW2011 and CP4. Unit A is characterized by dark gray clayey silt layers, common in the lowermost part of the west transect and core CP5. The top of unit A is comparable across cores CP2, CP3, BW2011, and CP6 at 2.65 m DBWS, while it is slightly higher at 2.5 m DBWS in core CP5 (Figures 1e, 2b and c). Unit A is overlaid by unit B, which is composed of brown clayey silt with a gradual lower boundary. Unit B is present in all cores, including the west transect, but becomes indeterminable in core CP5 (Figures 1e, 2b and c). Unit B is approximately 0.15 m thick, and the top of this unit roughly resembles the top of unit A in core CP5. Unit C comprised dark gray clayey silt, analogous to the sediments in unit A, and was typical in the bottom of cores in the east transect. The exception is core CP8, in which the lowermost layer was brown sand, and unit C was indeterminable (Figures 1e and 2a). In core CP4, unit C was overlaid on the bedrock preventing deeper penetration. The top of unit C was comparable at 2.65 m DBWS. Unit D is a combination of brown clayey silt and brown silty clay (unit D1) overlaid by light or dark gray or dark brown clayey silt (unit D2). Unit D1 is 0.7 m thick in the east transect cores but becomes thinner in core CP8 (Figures 1e and 2).
Unit D2 is present in the west transect cores, with varying thicknesses (Figures 1e and 2b). The transition between units D1 and D2 is gradual. Unit D2 is present in cores CP8 and CP9 in the east transect but undetermined in others (Figures 1e and 2a). Unit D2 is composed of dark brown clayey silt overlaid by dark gray clayey silt with a thickness of 0.25–0.3 m in cores CP3, CP6, and BW2011 (Figures 1e and 2c). Unit E was the uppermost unit and mainly consisted of dark gray gyttja clay overlaid by gyttja layers, found in all cores in the west transect except for core CP2, in which it was interbedded between gyttja clay and gyttja (Figures 1e and 2). The thickness of unit E varied among cores. The sedimentary units were identified by their characteristics and thickness, which varied among the cores.
Lithostratigraphic description of core CP4
Core CP4 comprises three sedimentary units, namely units C, D1, and E, further divided into 16 sedimentary layers. The lowermost part of core CP4 is composed of dark gray clayey silt with several angular to sub-angular gravels (about 5%) (layer 1). The layers above it comprise dark gray clayey silt (layers 2–4 and 6–9), separated by a layer of dark brown silty clay (layer 5). Unit C is gradually replaced by the brown clayey silt layers of unit D1 (layers 10–14), while Unit E is characterized by brown gyttja clay overlaid by clayey gyttja (layers 15 and 16). More detailed information can be found in Supplementary 2, available online.
Chronology of core CP4
Core CP4 was dated using nine radiocarbon measurements of charcoal, with a range of c. 3750–60 cal y BP (Figure 3 and Supplementary 3, available online). Two samples (D–AMS 034856 and 032441) from near the boundary between unit C and D1 at 2.58 and 2.64 m DBWS did not align with the general pattern, with one being broadly consistent with the age–depth model (D–AMS 034856) and the other being excluded (D–AMS 032441). This discrepancy may be due to gradual changes and mixing of sediments in the 6 cm thick samples used for dating. The sequential dates from unit C indicate a relatively fast sedimentation rate of 0.09–0.15 cm/year from 3.675 to 2.635 m DBWS, but this rate suddenly drops to 0.04–0.02 cm/year at the top of unit D1. The sedimentation rate in unit E is estimated to be 0.06 cm/year (Figure 4).

(a) Lithostratigraphy and (b) age-depth model constructed by Bacon (Blaauw and Christen, 2011) for core CP4 (see Supplementary 3, available online for details on the 14C dates).

Lithostratigraphy, sedimentation rate, LOI at 550°C and 950°C, the percentage of clay, silt, and sand fractions, grain size distribution, and mangroves, back mangroves, and mangrove associates pollen accumulation rate (PAR).
Loss on ignition (LOI)
The loss on ignition (LOI) at 550°C generally showed an opposite trend compared to the LOI at 950°C in the depth range between 3.645 and 1.8 m DBWS (Figure 4). The LOI at 550°C increased from 5% to 10%. In contrast, at 950°C, it decreased from 7% to 2.5% in this depth range (Figure 4). The LOI at 550°C increased to 20% at the top of unit D1 (layer 14: 2.045–1.930 m DBWS) and remained relatively stable in the lower layer of unit E (layer 15: 1.930–1.865 m DBWS) as detailed in Supplementary 2, available online and shown in Figure 4. The LOI at 550°C abruptly increased and reached its maximum at 30% in layer 16 (1.865–1.800 m DBWS). The LOI at 950°C slightly increased from 2.5% to 4% between 2.045 and 1.8 m DBWS (Figure 4).
Grain size analysis
Forty-two samples were analyzed using laser diffraction to determine their grain size distribution. Sedimentary unit C had a composition of around 70–80% silt and 20–30% clay (Figure 4). The sand content varied from 2% to 10%, and the variations in the sand fraction had a reverse relationship with the silt component in unit C (Figure 4). In the lower part of unit D1 (2.635–2.45 m DBWS), the percentages of silt and clay fractions were similar to those in the underlying layer of unit C. However, the sand fraction in this depth gradually increased from 5% to 8%, while the clay content significantly increased from around 20% to 40% (Figure 4). In contrast, the percentage of silt and sand significantly decreased from 80% and 8% to 50% and 1%, respectively, at 2.45 m DBWS. Above this depth, the variations in the percentage of clay and silt were inversely proportional, with the clay fraction progressively increasing from 40% to 50% and the silt component decreasing from 50% to 40%, resulting in a minor sand content above 2.45 m DBWS.
For unit E, the percentage of clay decreased significantly to 30% in layer 15, congruent with an increase in the silt content to 70% (Figure 4). In contrast, the grain size analysis in layer 16 demonstrated an increased clay component and reduced silt fraction. The volumetric mean diameters increased from 12 to 20 μm from 3.645 to 2.45 m DBWS, but suddenly declined to 5 µm at 2.45 m DBWS and only slightly changed after that (Figure 4).
Pollen analysis
The sedimentary sequence CP4 was divided into four pollen zones (Figure 5), and analysis was performed by dividing the pollen taxa into mangroves, back mangroves, and non-mangrove taxa.

Pollen percentage diagram of sediment sequence CP4, divided into four pollen zones based on the percentages of the major pollen taxa.
Zone 1: 3.645–2.25 m DBWS
Pollen zone 1 was primarily composed of mangrove pollen taxa. This zone was further divided into zones 1a and 1b based on the number of pollen grains found, with ranges of 2–300 and 200–300 grains/cm3, respectively. Zone 1a, which extended from 3.645 to 2.65 m DBWS (layers 3–9 and the lower part of layer 10) (Figure 5), was primarily composed of Rhizophora pollen at 80–90%, followed by Bruguiera pollen at around 5–10%. Small amounts of Avicennia and Sonneratia mangrove pollen taxa were also present. Pollen from Suaeda, Acrostichum, Stenochlaena, Poaceae, and Cyperaceae, which are back mangroves and non-mangrove taxa, respectively, were also found in small quantities. Although the percentages of Rhizophora pollen and Acrostichum spores were 50% at 3.00–3.01 m DBWS, only two pollen grains were found in this sample (Figure 5). Pollen was more abundant in zone 1b, which extended from 2.65 to 2.25 m DBWS and was mainly made up of Rhizophora and Bruguiera pollen at around 80%–90% and 5%–10%, respectively (Figure 5). Avicennia pollen was present at around 3%. In contrast, Suaeda pollen slightly increased from 1.5% to 4%, while the Rhizophora and Bruguiera pollen percentage declined at the top of this pollen zone (Figure 5).
Zone 2: 2.25–2.15 m DBWS
Pollen zone 2 consisted of the upper part of layer 12 and the lower part of layer 13, spanning from 2.25 to 2.15 m DBWS. The number of pollen grains in this zone increased from 7 to 300 grains/cm3 from the bottom to the top of this zone. In this zone, Rhizophora pollen, the most abundant pollen taxa in the preceding zone, decreased in proportion from around 30% to 3%. This zone was mainly occupied by Avicennia pollen, which comprised approximately 60%. The other mangrove pollen taxa, Bruguiera and Stenochlaena, were present in small amounts. Suaeda pollen declined from around 60% to 20%. However, the non-mangrove taxa were also present in small quantities, at around 5% (Figure 5).
Zone 3: 2.15–2.045 m DBWS
In pollen zone 3, located in the upper part of layer 13 (2.25–2.045 m DBWS), the pollen count suddenly dropped to 10–20 grains/cm3. The Avicennia pollen decreased to 40%, while other mangrove pollen taxa were unidentifiable (Figure 5). However, Suaeda pollen increased to around 40% at 2.05 m DBWS before reaching a maximum of 70% at 2 m DBWS. Acrostichum, a back mangrove taxon, made up 7% of the pollen at 2.05 m DBWS. On the other hand, Poaceae pollen accounted for approximately 20% of the pollen in zone 3.
Zone 4: 2.045–1.80 m DBWS
Pollen zone 4 was identified in sedimentary layers 14–16, which range from 2.045 to 1.8 m DBWS. Despite the low pollen count of fewer than 5 grains/cm3, a concentration of around 70 grains/cm3 was found at 1.9 m DBWS in layer 15. Non-mangrove and back mangrove taxa became more prevalent in this zone, as shown in Figure 5. Rhizophora and other mangrove pollen, such as Avicennia and Bruguiera, were present in only small amounts in some samples. In contrast, Suaeda pollen, the most dominant back mangrove pollen assemblage, accounted for approximately 20% of the pollen. The most prevalent non-mangrove taxa were Poaceae and Cyperaceae, making up about 30% and 20%, respectively (Figure 5).
Principal component analysis (PCA)
The Principal Component Analysis (PCA) included data on LOI, sediment components, and the most common pollen taxa (Rhizophora, Avicennia, Bruguiera, Suaeda, Acrostichum, Oncosperma, Cyperaceae, Poaceae, and Stenochlaena) as shown in Figure 6. The first two components of the PCA scores explained 60% of the cumulative variance. In comparison, the first four components explained approximately 77% of the total variance (Figure 6a). The eigenvalues of axes 1 and 2 were 6 and 2.5, respectively (Figure 6b), but for axes 3 and 4, it was around 1. Therefore, only the PCA of axes 1 and 2 were considered to identify the drivers of environmental changes in the Sam Roi Yot wetland over the last 4000 years (as detailed in Supplementary 4, available online).

LOI, sediment components, and the percentage of most frequent pollen taxa obtained from core CP4 were included in the principal component analysis (PCA). The PCA showing (a) the cumulative variance of the first two components is 60%, and of the first four components is 77% of the total variance and (b) the eigenvalues of axes 1 and 2 are the most dominant, and therefore (c) axes 1 and 2 were identified to play a role in the SRY wetland. Groups A and B (the light gray areas) and Groups C and D (in the white lines) are the data analyzed on axes 1 and 2, respectively.
Discussion
Environmental reconstruction in the SRY wetland
Based on the factor loadings, the pollen data could be categorized into four groups: A and B along axis 1, and C and D along axis 2 (Figure 6c). Group A mainly consisted of Poaceae, Cyperaceae, organic matter (OM) (from the LOI at 550°C), and the clay-fraction. Rhizophora was the only member of group B (Figure 6c). This is significant because it suggests that the dominant vegetation in this area at that time was mainly composed of Poaceae, Cyperaceae, and Rhizophora, similar to the present day. A 60% score for groups A and B explaining the cumulative variance supports this interpretation.
The wetland today is a freshwater marsh with dense Poaceae (Phragmites karka and Arundo donax) and Cyperaceae (Eleocharis dulcis) growth, with the eastern tidal flat dominated by Rhizophora – mangrove forest (Koskelainen, 2014; Niyomdecha, 2019; Parr et al., 1993; Sritipsak, 2014). The water level in the Sam Roi Yot wetland is mainly controlled by rainfall and small streamlines from the limestone mountain, runoff from Tenasserim Range, and saline water intrusion through the tidal creek (Koskelainen, 2014; Niyomdecha, 2019; Parr et al., 1993; Sritipsak, 2014). We, therefore, suggest that the RSL fluctuations influence the salinity of the wetland hydrology. This indicates that the corresponding environmental shift between a freshwater marsh and mangrove forest in the Sam Roi Yot wetland was dictated by the variability of RSL over longer timescales.
Group C was composed of carbonate (from the LOI at 950°C) (e.g. Heiri et al., 2001), silt-fraction, and a proportion of Rhizophora pollen (Figure 6c). Suaeda and Avicennia pollen and a fraction of clay clustered in Group D. Group C and D account for a total variance of 17%. This variance is likely related to precipitation and runoff because the carbonate content, the most dominant factor loading on axis 2, is primarily sourced from the limestone mountain ranges close to the wetland. The correlation between inorganic carbon input and the percentage of the silt and the inverse correlation between inorganic carbon and clay supports this interpretation. Additionally, the inverse relationship between the inorganic carbon input with the relative abundance of Suaeda, a hypersaline herb (Farooqui and Nautiyal, 2016; Hait and Behling, 2009, 2019), indicates a dry period. Dryer conditions will lower the water level, expanding the saline wetland fringe, thus favoring the relative increase of Suaeda, which is more adapted to hypersaline conditions (Farooqui and Nautiyal, 2016). Therefore, the variability of Suaeda pollen in the sedimentary record indicates an occasionally dry period caused by water level lowering and increasing salinity on the wetland.
Usually, typhoon strikes are the possible depositional mechanism for the coarse sediment in this area (Kongsen et al., 2021; Phantuwongraj et al., 2013; Williams et al., 2016). However, the absence of sharp contact with enclosing sediments, which are the critical criteria of storm surge sedimentation, rules out this possibility (Phantuwongraj et al., 2013; Williams et al., 2016). Consequently, based on PCA analysis, we suggest that the environmental shifts in the Sam Roi Yot wetland are significantly influenced by seawater intrusion, which is dictated by sea level changes in the Gulf of Thailand. Thus, the pollen records of the Sam Roi Yot wetland’s sediment core reflect the variabilities in the RSL in the Gulf of Thailand.
RSL fluctuation in the Gulf of Thailand during the Late-Holocene
After the Late-Holocene highstand at 4 ka, the mangrove pollen taxa, especially Rhizophora, were continuously abundant. In contrast, the back mangrove and non-mangrove taxa were very low, suggesting a continuous inflow of saline water to the wetland at 3750–1850 cal y BP (Figure 5). The variations in the sedimentation rate, sand fraction, and volumetric mean diameter were approximately congruent and demonstrated a tidal inundation fluctuation at c. 3750–1850 cal y BP (Figure 4). The mangrove pollen accumulation rate (PAR) during c. 3750–2950 cal y BP was significantly lower than at c. 2950–1850 cal y BP (Figure 5). This low mangrove pollen influx suggested that the mangrove pollen were possibly transported to the site, which was characterized by lagoonal conditions and a relative sea level of c. 2.6 m above present sea level (APSL) during c. 3750–2950 cal y BP (Surakiatchai et al., 2018) (Figure 7a and b). The gradual increases in the sedimentation rate, sand fraction, and grain size distribution indicated the increased RSL at c. 3750–3200 and 3050–2950 cal y BP, which was intervened by a lowering in the tidal inundation, possibly reflected regression at 3200–3050 cal y BP (Figure 4).

Coastal evolutions based on the Digital Elevation Model (DEM) at the RSL of (a) 3, (b) 2, and (c) 1 m APSL in Sam Roi Yot area. The location of core CP4 is the red circle.
The shift from unit C to unit D1 coincided with a significant decline in the sedimentation rate and a concurrent decrease in RSL of approximately 1.5 m APSL (Surakiatchai et al., 2018). According to a beach ridge study in Sam Roi Yot, the decrease in RSL of about 1.5 m led to the shrinkage of the wetland around 2950 cal y BP (Surakiatchai et al., 2018, Figures 2c, 4, 7b and c). This is consistent with the expansion of mangrove swamp forest at the coring site, as indicated by a significant increase in mangrove pollen influx, particularly of Rhizophora and Bruguiera, from 1000 to 5700 grains/cm2/year at 2950 cal y BP (Hicks and Hyvärinen, 1999) (Figure 4). Regression in Sam Roi Yot was slightly earlier than that obtained by the beach ridge study and the coastal dune formation at Chumphon (#5), which possibly began c. 3000 years ago (Miocic et al., 2022; Nimnate et al., 2015; Figure 1b and Supplementary 5, available online). However, it was not recognized in peat layers from the Chao Phraya Delta (#3) and the Straits of Malacca (#10), pollen records from Thale Noi (#8) and Kuantan (#9), and radiocarbon dating of oyster bands attached to the sea-notch at Ko Mae Koh (#6) and Ko Sichang (#2) (Geyh et al., 1979; Hesp et al., 1998; Horton et al., 2005; Oliver and Terry, 2019; Sinsakul, 1992; Somboon and Thiramongkol, 1992; Zhang et al., 2021; Figure 1b and Supplementary 5, available online). These discrepancies in the RSL reconstruction were possibly caused by the inability to capture this excursion with the relatively large dating uncertainties.
The enlargements in the sand and silt fractions coincided with the gradual rise in the mangrove PAR, which implied the continued intrusion of saline water accompanied by mangrove forest expansion at c. 2950–2450 cal y BP (Figure 4). The replacement of the silt with the clay component was possibly related to a reduced tidal inundation and RSL changes and environmental shifts caused by the thriving mangrove forests in the wetland at c. 2450–1850 cal y BP. These observations together supported the gradual shrinkage of the lagoons, possibly available until 1850 cal y BP. The increase in Suaeda – a back mangrove taxon – was concurrent with the decline in Rhizophora, thus reflecting a fall in the RSL at c. 1850 cal y BP (Figure 5). The meager pollen grain count seemed to be caused by wetland exposure and oxidation, implying the continued fall in the RSL until 1450 cal y BP. However, since the significant change was undeterminable in the sediment compositions at c. 1850–1450 cal y BP, we assumed that the wetland was possibly exposed seasonally when the RSL was about 1 m APSL, regarded to beach ridge at Sam Roi Yot and dating of the oyster band at Ko Phaluai (#7) (Oliver and Terry, 2019; Surakiatchai et al., 2018; Figures 1b, 5 and 7c, and Supplementary 5, available online).
Mangrove pollen taxa, especially Avicennia, became dominant from 1450 to 1050 cal y BP (Figure 5). Although Avicennia is a salt-tolerant plant that can grow over a wide range of Watson’s inundation classes, its increased pollen dominance in these samples implied an increased tidal inundation reflecting an excursion of transgression at 1450–1050 cal y BP (Santisuk, 1983; Somboon, 1990). Suaeda, Avicenna, and non-mangrove taxa (Poaceae and Cyperaceae) gradually increased after that (Figure 5). The timing of the RSL fall is consistent with the age of the rock shelter paintings in the Sam Roi Yot Mountain (The Fine Arts Department, 2004). These RSL fluctuations overlapped with those obtained from the pollen records from a mangrove forest in Koh Chang (#1) that demonstrated the RSL fall at c. 1500–1300 cal y BP and a rise at c. 1300–500 cal y BP (Englong et al., 2019) (Figure 1b and Supplementary 5, available online). The different environmental conditions between Ko Chang and the Sam Roi Yot wetland could explain the inconsistency in the timing of the sea-level changes. However, by c. 500 cal y BP, the percentage of Cyperaceae and Poaceae pollen gradually increased, suggesting an evolution of the wetland to a freshwater marsh, similar to the present (Figure 5) (e.g. Koskelainen, 2014; Niyomdecha, 2019; Parr et al., 1993; Sritipsak, 2014).
Impact of ENSO on the RSL fluctuations in the Gulf of Thailand
The Holocene sea level fluctuations in Southeast Asia have been primarily attributed to hydro–isostasy, which is the adjustment of the Earth’s crust to changes in the volume of ice sheets and glaciers (Horton et al., 2005). However, the spatial variability of the magnitude of sea level fluctuations varies depending on the geographical location within the region (Chabangborn et al., 2020; Choowong et al., 2009; Jiwarungrueangkul et al., 2022; Oliver and Terry, 2019; Sainakum et al., 2021). For the Late-Holocene, the global sea level has been reported to decrease gradually, as simulated by the ICE–6E_C Glacial Isostatic Adjustment (GIA) model (Peltier et al., 2015), or show no evidence of fluctuation (Lambeck et al., 2014). However, the predominance of Rhizophora, a mangrove species, at the bottom of core CP4 suggests that the sea level rose prior to c. 4000 cal y BP, concurrent with the reconstruction of the sea level highstand at Thale Noi, Pattalung (#8) (Horton et al., 2005) (Figures 1b and 5, and Supplementary 5, available online). The discrepancies between the GIA model simulation and the sea level reconstruction allow us to examine the other possible mechanisms interacting with hydro-isostasy in the Gulf of Thailand (Peltier et al., 2015).
Due to the region’s tectonic stability, tectonic adjustment is unlikely to contribute to the sea level changes in the Gulf of Thailand. However, the transgression in the Gulf of Thailand coincides with the abrupt climate shift at the 4-ka event, which was related to the cold phase increasing ice rafted debris in the North Atlantic (Berkelhammer et al., 2012; Bond et al., 1997; Wang et al., 2016), and potential connection to the El Niño Southern Oscillations (ENSO) in the Equatorial Pacific (Dang et al., 2020; Toth and Aronson, 2019).
Previous studies have shown that ENSO dominates rainfall dynamics in Thailand over multiple timescales (Buckley et al., 2014; Yamoah et al., 2016, 2021). In addition, the crucial role in the inflow from the South China Sea to the Gulf of Thailand is the easterly wind (Sojisuporn et al., 2010), which is directly related to ENSO. Therefore, the large-scale atmosphere–ocean mechanisms of ENSO variability, manifested by alternating warm El Niño and cold La Niña conditions, will significantly influence the volume of water within the Gulf of Thailand. Indeed, this is consistent with modern observations based on sea level changes at Ko Lak, Thailand (Figure 1b). By comparing the mean annual sea surface temperature (SST) anomaly at Niño 3.4 region (5°N–5°S, 120°–170°W) – a widely used marker of El Niño variability – with sea level changes at Ko Lak, it was observed that the sea level in the Gulf of Thailand gradually fell during prolonged warm SST in 1990–1995 and 1997 whereas the La Niña years in 1988–1989, 1995–1996, and 1998–2000 was characterized by higher sea levels (National Centers for Environmental Information, 2023; WOCE Sea Level, 2006; Figure 8). This means that ENSO can significantly impact the sea level changes in the Gulf of Thailand and that during times of extended El Niño or La Niña conditions in the past, the sea level would have been consistently higher or lower than expected from eustatic and isostatic processes alone during the Holocene period. Therefore, we hypothesized that the long-term fluctuations of ENSO possibly contributed to sea level fluctuations in the Gulf of Thailand during the Late-Holocene.

Comparison between (a) the sea surface temperature (SST) anomaly at Niño 3.4 region (5°N–5°S, 120°–170°W) (National Centers for Environmental Information, 2023) and (b) the monthly sea level from coastal tide gauge station Ko Lak, the Gulf of Thailand (WOCE Sea Level 2006).
Conclusion
Our study suggests that environmental changes in the Sam Roi Yot wetland are significantly influenced by seawater intrusion and that the sedimentary record can reflect fluctuations in the RSL in the Gulf of Thailand. The predominance of Rhizophora in the bottom of core CP4 suggests that the RSL in the Gulf of Thailand rose prior to around 4000 years before the present. This coincides with the abrupt climate shift at the 4-ka event, which was related to the cold phase in the North Atlantic, and ENSO in the Equatorial Pacific. Although RSL fluctuations in Southeast Asia have been attributed to hydro-isostasy, the magnitude of RSL fluctuations reported differs based on the geographical location within the region, indicating other possible localized mechanisms such as ENSO. We, therefore, hypothesize that the long-term fluctuations of ENSO possibly contributed to RSL fluctuations in the Gulf of Thailand during the Mid-/Late-Holocene high stand. After the Late-Holocene highstand at 4 ka, mangrove pollen, especially Rhizophora, were continuously abundant. At the same time, back-mangrove and non-mangrove taxa were very low, indicating a continuous inflow of saline waters to the wetland. The variations in the sedimentation rate, sand fraction, and volumetric mean diameter demonstrated a tidal inundation fluctuation. Between c. 3750 and c. 2950 cal y BP, the mangrove pollen accumulation rate was significantly lower compared with the period from c. 2950 to 1850 cal y BP indicating that the mangrove pollen were possibly transported to the site. The lowering in RSL gradually shrunk the lagoon after c. 2950 cal y BP. Regression in the Sam Roi Yot area was slightly earlier than that obtained by the beach ridge study and the coastal dune formation at Chumphon.
Supplemental Material
sj-pdf-1-hol-10.1177_09596836231197745 – Supplemental material for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene
Supplemental material, sj-pdf-1-hol-10.1177_09596836231197745 for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene by Akkaneewut Jirapinyakul, Worakamon Nudnara, Paramita Punwong, Robin Nohall, Apichaya Englong, Chawisa Phujareanchaiwon, Kweku Afrifa Yamoah and Montri Choowong in The Holocene
Supplemental Material
sj-pdf-2-hol-10.1177_09596836231197745 – Supplemental material for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene
Supplemental material, sj-pdf-2-hol-10.1177_09596836231197745 for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene by Akkaneewut Jirapinyakul, Worakamon Nudnara, Paramita Punwong, Robin Nohall, Apichaya Englong, Chawisa Phujareanchaiwon, Kweku Afrifa Yamoah and Montri Choowong in The Holocene
Supplemental Material
sj-pdf-3-hol-10.1177_09596836231197745 – Supplemental material for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene
Supplemental material, sj-pdf-3-hol-10.1177_09596836231197745 for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene by Akkaneewut Jirapinyakul, Worakamon Nudnara, Paramita Punwong, Robin Nohall, Apichaya Englong, Chawisa Phujareanchaiwon, Kweku Afrifa Yamoah and Montri Choowong in The Holocene
Supplemental Material
sj-pdf-4-hol-10.1177_09596836231197745 – Supplemental material for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene
Supplemental material, sj-pdf-4-hol-10.1177_09596836231197745 for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene by Akkaneewut Jirapinyakul, Worakamon Nudnara, Paramita Punwong, Robin Nohall, Apichaya Englong, Chawisa Phujareanchaiwon, Kweku Afrifa Yamoah and Montri Choowong in The Holocene
Supplemental Material
sj-pdf-5-hol-10.1177_09596836231197745 – Supplemental material for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene
Supplemental material, sj-pdf-5-hol-10.1177_09596836231197745 for ENSO may have contributed to sea level changes in the Gulf of Thailand during the Late-Holocene by Akkaneewut Jirapinyakul, Worakamon Nudnara, Paramita Punwong, Robin Nohall, Apichaya Englong, Chawisa Phujareanchaiwon, Kweku Afrifa Yamoah and Montri Choowong in The Holocene
Footnotes
Acknowledgements
The authors thank the Department of National Parks, Wildlife, and Plant Conservation for collecting samples, Robert Butcher for proofreading and assistance during manuscript preparation, and all those who helped in diverse ways to make this project a success.
CRediT author statement
Akkaneewut Jirapinyakul: Conceptualization, Investigation, Formal analysis, Visualization, Writing original draft preparation, Writing – review & editing, Project Administration, Funding acquisition. Worakamon Nudnara: Investigation, Formal analysis, Writing original draft preparation, Writing – review & editing. Paramita Punwong: Conceptualization, Investigation, Formal analysis, Writing original draft preparation, Writing – review & editing. Robin Nohall: Investigation and Formal analysis, Writing -review & editing. Apichaya Englong: Investigation and Formal analysis. Chawisa Phujareanchaiwon: Investigation and Formal analysis. Kweku Afrifa Yamoah: Resources, Conceptualization, Writing – review & editing. Montri Choowong: Resources, Conceptualization, Writing – review & editing.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was financed through the Fundamental Fund (2021), Chulalongkorn University (CUFRB65_dis(11)_099_23_29). Worakamon Nudnara acknowledges financial support from the Faculty of Science, Chulalongkorn University, under the research assistant grant (2018).
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
