Response and recovery of a Sphagnum peatland from long-term human-induced alkalinisation

Northern peatlands are significant terrestrial carbon stores but are increasingly threatened by human activities. Ombrotrophic peatlands, being naturally acidic, are particularly vulnerable to alkaline pollution. Despite their importance, the effects of alkalinisation on peatlands remain insufficiently studied. In Estonia, alkaline pollution from a cement industry and oil shale power plant emissions have degraded several peatlands since the 19th century. Although some sites have recovered in recent decades, more severely impacted areas remain in poor condition.

We investigated the effects of alkalinisation on Varudi peatland, a forested site in northeast Estonia, which was exposed to 125 years of alkaline emissions from a nearby cement factory. Using a multi-proxy, high-resolution palaeoecological approach combined with a precise and reliable age-depth model, we reconstructed changes in environmental, chemical, botanical, and hydrological conditions over the past millennium. Our findings revealed three successional phases: during the mid-13th century CE, land clearance and increased mineral deposition caused the site to transition from a bog to a poor fen phase between approximately 1250–1570 CE; and while the cement factory operated without efficient filters, the site became a pine-dominated fen between 1871–1995.

After the installation of filters in 1996, peatland pH returned to near pre-disturbance levels, and some recovery was observed. However, the site remains degraded. Our results indicate that alkalinisation significantly disrupts peatland functioning, reducing carbon storage and altering vegetation communities. These effects can persist for decades even after the source of contamination is removed, underscoring the need for more comprehensive monitoring of peatlands impacted by alkaline pollution globally.

1 Introduction

Despite only covering ca. 3 % of the Earth's surface, northern peatlands contain >500 Gt of carbon (Bridgham et al., 2008; Yu et al., 2010; Yu, 2012). Their capacity to accumulate and store carbon results from the waterlogged and acidic nature of their soils (Clymo et al., 1998). These conditions preserve organic material, which accumulates and may be stored indefinitely (Harenda et al., 2018). Since their initiation, peatlands have slowly removed carbon from the atmosphere, imparting a weak but persistent cooling effect upon global climate over millennial timescales (Frolking et al., 2006).

Despite being recognised as a valuable tool for climate change mitigation, peatlands still receive little protection, regionally or nationally (Rawlins and Morris, 2010). As of 2018, approximately 10 % of the remaining peatlands worldwide are in a degraded state (Leifeld and Menichetti, 2018), while in Europe this rises to 25 % (Tanneberger et al., 2021). Such disturbance can disrupt the fragile hydrological balance that maintains the carbon sink function of peatlands and may cause them to shift from sinks to sources of atmospheric carbon, exacerbating climate change (Leifeld and Menichetti, 2018).

Estonia is one of Europe's most peat-rich countries, with peatlands covering ca. 22.5 % of its land area (Orru and Orru, 2008). Due to this abundance, peat is a significant natural resource for Estonia, and has been heavily exploited, particularly after the Industrial Revolution (Paal et al., 2010; Łuców et al., 2022). The rise in anthropogenic pollution since this time has caused substantial changes in global geochemistry such that this era is informally termed “the Anthropocene” (Fiałkiewicz-Kozieł et al., 2018; Waters et al., 2023). During this time in Estonia, emissions from industrial sources were characterised by high levels of calcium-rich particulate matter, with most of the emissions concentrated in the northeastern industrial region of the country (Liblik et al., 1995; Karofeld, 1996). Ombrotrophic peat bogs, which are the dominant type of peatlands in Estonia, being naturally acidic and nutrient-poor ecosystems are particularly sensitive to alkaline atmospheric pollution (Paal et al., 2010). These emissions caused significant changes in the geochemical and botanical composition of bogs adjacent to pollution sources, which resulted in dramatic increases in pore-water pH and losses of bog-specific vegetation (including Sphagnum mosses) near affected sites (Paal et al., 2010; Vellak et al., 2014). This was often followed by encroachment of Pinus sylvestris (Scots pine) onto polluted sites alongside other species typical of nutrient-rich alkaline environments (Pensa et al., 2004, 2007; Ots and Reisner, 2006; Kaasik et al., 2008; Kask et al., 2008).

By the 1990s, industrial emissions in Estonia began to fall, following a decline in power generation and improved filtration systems in factories (Liiv and Kaasik, 2004). Following these reductions, polluted peatland sites began to show signs of recovery, with acidic conditions and bog-specific vegetation returning (Karofeld, 1996; Kaasik et al., 2008; Paal et al., 2010). However, in more heavily polluted sites, this recovery has been slow, and the impact of past alkaline pollution persists to this day in some areas (Ots and Reisner, 2006). It remains unclear whether current levels of atmospheric pollution are sufficiently low to permit their full recovery in the future, or how long this process will take (Paal et al., 2010). Despite growing concerns over alkaline pollution and its potential future effect on peatlands, particularly concerning their role as carbon reservoirs, research exploring the effects of alkalinisation upon peatland ecosystems and their subsequent recovery has been limited.

Atmospheric pollution remains a significant threat to peatland ecosystem functioning (Bobbink et al., 1998; Turetsky and St Louis, 2006; Osborne et al., 2024). The effects of alkaline pollution upon peatlands have been relatively overlooked relative to those of acid rain due to its effects being more localised (Vellak et al., 2014; Sutton et al., 2020). However, nearly two billion tonnes of alkaline residues are emitted into the atmosphere each year (Gomes et al., 2016). Despite environmental standards curbing emissions in recent decades, in some areas these regulations are not consistently enforced or are merely declarative (Abril et al., 2014; Ivanov et al., 2018). Following global reductions in acid rain since the 1980s, the relative proportion of alkaline pollutants in airborne particulate matter has increased in the UK, much of Europe, North America and China since 1986 (Turetsky and St Louis, 2006; Sutton et al., 2020). Additionally, climate change may exacerbate the effects of increased alkalinity in boreal regions, as permafrost thawing may cause the expansion of areas of open water, increasing surface runoff and infiltration in some regions (Walvoord and Kurylyk, 2016). This may allow for longer contact times between surface water and carbonate bedrock, accelerating weathering and raising the pH of surface waters which may then enter peatlands in the surrounding catchment (Schindler, 1997; Osterkamp et al., 2000; Lehmann et al., 2023).

In this study, we focus on the effects of alkalinisation resulting from over 125 years of intense emissions from a nearby cement factory upon Varudi soo (bog), a formerly ombrotrophic peatland in northeastern Estonia. By employing a high-resolution, multi-proxy palaeoecological framework, we reconstruct changes in the chemical, botanical, hydrological and environmental conditions of the site over the past millennium to address the following questions:

1. What is the current state of a heavily polluted raised bog almost 30 years after the reduction in alkaline pollution?

2. How has alkaline pollution altered the ecosystem functioning of the site and how does this compare with pre-disturbance conditions?

3. To what extent has this ecosystem function recovered 30 years after removing the point source of pollution?

4. Can we identify critical transitions that can be broadly applied to assess peatland condition and recovery following alkaline pollution?

2 Methods and materials

2.1 Study area

Varudi bog (59°26^′19^′′ N, 26°35^′13^′′ E) is located in Lääne-Virumaa, northeastern Estonia, consisting of fen-bog habitats. The site is approximately 10 km south of the coast of the Gulf of Finland and Baltic Klint (Fig. 1) and spans ca. 12.6 km². The site is primarily a forested Sphagnum bog interspersed by numerous bog pools and hummocks, with an overstory of Pinus sylvestris. Varudi peatland receives 478 mm of rainfall per year, has a mean annual temperature of 7.3 °C, and prevailing winds are from the southwest and south. The underlying bedrock is composed of Cambrian and Lower Ordovician siliciclastic sedimentary rocks and Middle Ordovician limestones (Sibul et al., 2017) covered with a relatively thin layer of glacial and post-glacial sediments.

Figure 1 Map of study locations. (A) Modern-day Estonia (in Gray). The red shaded square in map (A) indicates the area mapped in map (B). (B) Locations of sites relevant to this study (Kunda cement factory: Red, Aru-Lõuna Quarry: Yellow, Varudi peatland: Blue). (C) Satellite image of Varudi bog (© Microsoft) showing the extent of peat cutting and drainage that has taken place over the past century. The red shaded area indicates the area mapped in map (D). (D) Close-up of the location from where core VAR1 was recovered (red star). (E) Coring location.

Figure 2 Photographs of the Kunda cement factory and locale taken during the 1980s and 1990s, showing (A) Cement dust emissions from the chimneys. (B) Photograph of chimneys and factory surroundings. (C) Photograph of the factory workings, showing substantial cement dust deposition in the surrounding area. (D) Photograph showing cement dust deposition upon house near the factory. Photography by: (A, D) Atko Heinsalu taken in the early 1990s; (B) Estonian National Archives Photo Database (code EFA.204.0.268452) August 1994 Albert Truuväärt. (C) Estonian National Archives Photo Database (code EFA.204.0.168924) August 1989 Tiit Veermäe.

Varudi was selected for study due to its proximity (7.5 km NW) to the Kunda Nordic Tsement Factory (henceforth Kunda cement factory) (Figs. 1 and 2), which has been operational since the 1870s. During the late 1970s, cement production peaked at ca. 1.2 million t yr⁻¹ (Fig. 3). The dust emissions from the cement plant have fluctuated between 45 000 and 99 000 t yr⁻¹ during the last decades of the 20th century, with the highest dust emissions recorded in 1991 (Ots and Mandre, 2012). In addition to emissions from the cement factory, the site has also received emissions from nearby alkaline generating industries and oil shale power plants, including the Balti power plant, located approximately 100 km west, and the Aru-Lõuna limestone quarry, situated 8 km to the east (Karofeld, 1996, Fig. 1) which supplies raw materials for cement production at Kunda. However, significant emission reductions occurred after 1996, following the installation of pollution control filters, which lowered emissions from 14 000 t in 1996 to 530 t in 2000 and 8 t in 2020 (Fig. 3).

The input of alkaline cement dust emissions generated drastic increases in porewater pH. During 1996–1997, the pH levels taken from peat pools at Varudi varied from 7.6 to 8.5, whereas in natural conditions bog pools have a pH in the range of 3–4. Disturbances at the site are compounded by drainage and peat harvesting for horticulture, which continues to the present day, with harvesting affecting around 40 % of the former centre of the site (Fig. 1C). There is limited documentary information available related to the timing of drainage and peat cutting at Varudi, although historic maps (EtoMesto, 2025; Maa-ja Rummiamet, 2025) indicate that most of the extensive drainage took place before 1977 and after 1944, probably relating to the large-scale drainage projects undertaken in Estonia in the 1960s under Soviet rule (Paavilainen and Päivänen, 1995). Peat extraction matching the areal extent of present-day works at the site by 1988. Prior to this, there is documentary evidence for drainage in the northern, eastern and western margins of the site having taken place around 1866–1911, likely indicative of the systematic forest drainage that took place in Estonia around this time (Paavilainen and Päivänen, 1995).

Figure 3 Cement annual production and emissions by the Cement Factory in Kunda. Data adapted from Trumm et al. (2010), digitised from a figure at page 209 using WebPlotDigitiser: https://automeris.io/WebPlotDigitizer/ (last access: 18 November 2024) Emissions data (available from 1985 to 2023) sourced from Raukas (1993); Pärtma (2023) and Heidelberg Materials Kunda AS (2023). Note the logarithmic scale for emission values.

The chemical composition of the cement dust emissions from Kunda consists primarily of CaO (12 %–17 %), SiO₂ (6 %–9 %), and several other trace metals, including lead (Pb, ca. 60 mg kg⁻¹), cadmium (Cd, ca. 0.9 mg kg⁻¹), and zinc (Zn, 129 mg kg⁻¹) (Mandre and Ots, 1999). The pH of the cement dust in the water solution ranged between pH 12.3 to 12.6 (Mandre and Korsjukov, 2007).

2.2 Coring method

In August 2022, an 86 cm peat core (VAR1) was extracted from Varudi peatland using a 1 m-long Wardenaar peat corer (Wardenaar, 1987, coring location: 59°26^′23^′′ N, 26°35^′55^′′ E, Fig. 1D). The core site was located near an actively harvested peat area but is within an intact section of the site in what was originally the bog's central raised dome. Vegetation was characterised by Pinus sylvestris, with the ground cover dominated by Eriophorum sp., Menyanthes trifoliata, Vaccinium sp. and scattered Betula nana. Both Sphagnum and brown mosses were also present. Despite its proximity to a drainage ditch, the coring location was representative of the overall condition of the site.

Following recovery, the core was wrapped in plastic and transported to the Faculty of Geographical and Geological Sciences at Adam Mickiewicz University, Poznań, Poland for analysis. The core was stored at 4 °C prior to sub-sampling.

2.3 Dating methods and age-depth model

Identifiable above-ground plant macrofossils were picked from 1 cm thick sub-samples taken at various depths throughout the core, following frameworks by Piotrowska et al. (2011) and Nilsson et al. (2001). Sphagnum stems, branches and leaves were preferentially used for dating where present. Where these were not available, above-ground remains of ericaceous plants (leaves, stems, seeds) were used instead. Initial samples were taken at 20 cm intervals throughout the core to establish a baseline chronology, which was then used to determine where additional samples would be selected. A total of 12 samples were sent for radiocarbon analysis. Each sample was pre-treated using the acid-base-acid approach and analysed by accelerator mass spectrometry (AMS) at the Poznań Radiocarbon Laboratory, Poland.

To provide a reliable chronology for recently accumulated peat, ²¹⁰Pb and ^239+240Pu analyses were used on material from the upper 40 cm of the core at 1 cm contiguous resolution, following methods outlined by Appleby (1998). Peaks in the activity of ^239+240Pu, which are linked to nuclear fallout events (e.g., 1950s atmospheric nuclear tests), were also measured as independent time-markers to validate and supplement the age-depth model (Mróz et al., 2017; Cwanek et al., 2021). Samples were processed at the Polish Academy of Sciences' Institute of Nuclear Physics, Kraków, Poland using an AlphaAnalyst™ 7200 spectrometer (Mirion Technologies). Quality control and accuracy was ensured by measuring blanks and certified reference material (IAEA 447, IAEA 385) alongside the core samples (results are provided in Table S1. The activity concentration of ²¹⁰Pb ($T\mathrm{1}/\mathrm{2}=\mathrm{22.3}$ years) was estimated by measurement of its decay product – ²¹⁰Po ($T\mathrm{1}/\mathrm{2}=\mathrm{138.4}$ d), while ^239+240Pu was measured directly. Activity concentrations of ²¹⁰Pb and ^239+240Pu are reported in units of Bq kg⁻¹.

The age-depth model was constructed by integrating all ¹⁴C, ²¹⁰Pb, and ^239+240Pu data within a Bayesian framework in R using the package “rplum” (rplum package, Aquino-López et al., 2018; R Core Team, 2023). rplum generates maximum age probabilities at user-defined intervals (here every 1 cm), together with maximum and minimum ages based upon calculated 95 % credible intervals. This method allows for the integration of ¹⁴C dates with the ²¹⁰Pb dates without the need for re-modelling (Aquino-López et al., 2018, 2020). Radiocarbon dates were calibrated in rplum using the INTCAL20 curve, with post-1950 samples using the BOMB1 curve for the Northern Hemisphere (Reimer et al., 2020; Uno et al., 2013). In this study, the resultant ages are expressed as calendar years (cal) CE, with 0 BP equal to 1950 cal CE.

2.4 Palaeoenvironmental proxies

2.4.1 Testate Amoebae

Samples were processed following a modified version of protocols by Hendon and Charman (1997). Samples were placed into 50 mL centrifuge tubes filled with deionised water and agitated for ca. 10 min. These were sieved through a 300 µm mesh and the smaller fraction was retained. Sieved samples were centrifuged at >3000 rpm for 5 min and a sub-sample of the resultant material was transferred to a microscope slide for identification at 400× magnification. Samples were not heated or micro-sieved, to retain small but ecologically sensitive species, as recommended by Avel and Pensa (2013). A minimum count of 100 tests for each sample was considered statistically significant (Payne and Mitchell, 2009). Tests were identified to the species level where possible, with reference to Siemensma (2023) and Mazei and Tsyganov (2006) and were later pooled into taxonomic groups defined by Amesbury et al. (2016). The relative abundance (%) of each taxa count was calculated for each sample. Water table depth was reconstructed using the pan-European tolerance down weighted with inverse de-shrinking transfer function model by Amesbury et al. (2016), based on a training set of 1302 samples spanning 35° of latitude and 55° of longitude. Small species (<10–25 µm, broadly oval-shaped) not included in these groupings but present in the core from Varudi were grouped under the “Cryptodifflugia oviformis” group. We also reconstructed changes in pH using a transfer function developed by Šímová et al. (2022), based on over 250 samples from east-Central Europe. The calibration dataset includes samples from peatland sites ranging from acidic bogs to calcareous fens, with a pH range of 3.5–7.9 making it suitable for reconstructing pH in sites that have likely undergone substantial changes in pH over time, such as Varudi. The transfer function was constructed using maximum likelihood calibration of Gaussian response curves (MLRC) with squared chord distance as the dissimilarity measure. The full dataset transfer function (mineral and ombrotrophic sites combined) was tested by leave-one-out cross-validation, achieving a RMSEP₁₀₀ of 0.45 pH units and an R² of 0.86. Both reconstructions were performed using the “Rioja” package in R (Juggins, 2019). Comparisons between the modern calibration and fossil datasets using the “compare.datasets” function in Rioja indicate good levels of overlap, with no non-analogue assemblages identified in the fossil data.

Stratigraphically Constrained Cluster Analysis (CONISS) was used to quantitatively define stratigraphic zones in the sub-fossil testate amoeba data (Grimm, 1987). This method is used to determine statistically significant zones, reflecting changes in testate amoebae community composition. The data were square root transformed prior to applying CONISS. The numbers of statistically significant zones were determined using Broken-Stick modelling.

We applied the framework outlined by Burge et al. (2023) to identify significant ecological transitions in the testate amoeba communities due to disturbance. Sub-fossil data were Hellinger transformed and analysed using the “prcurve” function in the “analogue” package in R (Simpson and Oksanen, 2025). Principal response curves (PrC), which reduce multi-dimensional community data to a single-dimensional curve, calculating the (dis)similarity between sample scores indicative of the difference between samples (De'ath, 1999; Van den Brink and Braak, 1999) provided the best fit for our data. A generalised additive model (GAM), an approach effective in capturing rapid and non-linear changes in palaeoecological studies, was applied to the PrC data to account for temporal autocorrelation (Auber et al., 2017; Beck et al., 2018; Burge et al., 2023). Given the abrupt changes in our data, we also applied an adaptive spline GAM following Burge et al. (2023). However, this method does not account for temporal autocorrelation and thus remains incompatible with the GAM framework as outlined by Simpson (2018).

2.4.2 Plant macrofossils

Plant macrofossil analysis followed procedures adapted from Mauquoy et al. (2010). Samples of approximately 5 cm³ were sieved through a 200 µm test sieve and the larger fraction was retained. Botanical composition was estimated as percentages under a low-powered microscope at between 10–100 times magnification, using a 10×10 grid eyepiece graticule to aid quantification of plant remains. Seeds, fruits, spindles, leaves, and wood were also identified at the species level where possible and counted as individual counts. Identification was aided by identification guides (Katz et al., 1965, 1977; Grosse-Brauckmann, 1972, 1974; Tobolski, 2000; Mauquoy and van Geel, 2007; Bojnanský and Fargašová, 2007). Sphagnum remains were identified to sub-generic sections, with 100 leaves examined per sample, where possible, to calculate the relative abundance of each sub-section as a percentage of the total leaves identified.

2.4.3 Pollen, non-pollen palynomorphs and microscopic charcoal analyses

Past changes in vegetation cover at the landscape scale were assessed using pollen analysis (Seppä and Bennett, 2003). A total of 22 samples were prepared following the laboratory procedures outlined by Berglund and Ralska-Jasiewiczowa (1986). Each 1 cm thick sample, measuring 1 cm³, was sub-sampled at 5 cm intervals throughout the core. The samples were treated with 10 % potassium hydroxide (KOH) to remove humic compounds before acetolysis. A Lycopodium tablet (batch no: 280521291, 18 407 spores per tablet; Manufacturer: Lund University) was added to each sample to calculate pollen concentrations, following methods by Stockmarr (1971). Samples were transferred to microscope slides and mounted in glycerine jelly for analysis. Pollen, spores, and selected non-pollen palynomorphs (NPPs) were identified and counted using a high-powered stereo microscope. Identification was based on established atlases and keys (Pollen: Moore et al., 1991; Beug, 2009; NPPs: van Geel, 1978; van Geel and Aptroot, 2006; Miola, 2012). Although a target of 500 terrestrial pollen grains was aimed for per sample, this count was not always achievable due to low pollen concentrations in some core sections. The relative abundance of spores and NPPs was calculated as a proportion of the terrestrial pollen sum (TPS), which includes both arboreal (AP) and non-arboreal (NAP) pollen, excluding aquatic and wetland plant spores, Ericaceous pollen, and NPPs.

In addition to pollen counts, microscopic charcoal particles (10–100 µm) were counted from pollen slides as past fire activity, both natural and anthropogenic (Finsinger and Tinner, 2005), while spheroidal carbonaceous particles (SCPs) counted as indicators of industrial activity (Patterson et al., 1987; Swindles, 2010). Microcharcoal concentrations per cm³ were calculated by dividing the number of particles counted with the number of Lycopodium spores and multiplying this by the total number of particles counted, and accumulation rates (reported as particles cm⁻³ yr⁻¹) were calculated by dividing the microcharcoal concentration with the rPlum-derived age increments per cm slice.

2.4.4 Apparent rates of carbon accumulation

To measure apparent rates of carbon accumulation (aCAR), contiguous 1 cm-thick sub-samples were taken throughout VAR1. The volume of each wet sample was determined by water displacement, and then the samples were dried in an oven at 105 °C until no further weight loss occurred. Dry bulk density was calculated by dividing the wet volume by the dry mass of each sample. Organic matter content (derived from LOI) was determined by ashing the samples at 550 °C for 6 h, following the method of Chambers et al. (2011).

The carbon content of each sample was estimated indirectly by multiplying the LOI content by 0.52, based on the average ratio of organic carbon (OC) and LOI in ombrotrophic peat from multiple studies (Ball, 1964; Dean, 1974; Gorham, 1991; Clymo et al., 1998). Carbon density was calculated by multiplying the dry bulk density (g cm⁻³) by the percentage of carbon content, as described by Chambers et al. (2011). Apparent carbon accumulation rates (aCAR) were then calculated by dividing the carbon density from each peat slice by the sedimentation rate, determined from the age-depth model (Young et al., 2019, 2021).

2.4.5 µXRF-Core Scanning (ITRAX)

To identify the section of the core affected by cement dust pollution, we followed methods similar to those used by Varvas and Punning (1993) to assess pollution histories from Estonian lake sediments. They identified rapid increases in micro-element concentrations associated with alkaline fly-ash emissions, accompanied by a decrease in organic matter, indicating the presence of particulate emissions from oil-shale combustion by power plants.

In this study, the concentrations of geochemical elements throughout the core were measured using an ITRAX µXRF core scanner equipped with a molybdenum X-ray tube. Element concentrations were quantified as counts per second, based on the number of secondary fluorescence detected for each element over a given period. Measurements were taken at 5 mm intervals (30 kV, 50 mA, exposure time: 30 s per step) at the Institute of Geography, University of Bremen, Germany. The scanner identified the activity of the following elements: Al, Si, P, S, Ar, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ga, Ge, Se, Br, Rb, Sr, Y, Ba, Pt, Pb, Bi, and Fe.

The mean concentrations of each element were calculated for each 1 cm slice of the core. To account for variability in element concentrations due to sedimentation throughout the core, the results were normalised by dividing the total counts by the sum of the coherent and incoherent peaks, following the approach of Orme et al. (2015) as recommended by Longman et al. (2019). The chemical signature of cement dust was identified based on the composition of clinker emissions from the Kunda cement factory, as detailed by Klõšeiko et al. (2011) (Clinker emissions being the emissions produced by the production of clinker, the primary ingredient of cement, created by the high-temperature heating of limestone and other materials). Due to the uneven and unconsolidated nature of the top 6 cm of the core, this section was not scanned and was removed prior to analysis.

Principal Component Analysis (PCA) was used to summarise patterns of variation in the geochemical data, using the “vegan” package in R (Oksanen et al., 2019). To mitigate scaling effects, the data were standardised to z scores. The analysis was conducted with varimax rotation in correlation mode, to explore correlations between elements and organic matter percentage (Silva-Sánchez et al., 2014). The number of components to retain for analysis was determined using a Broken-Stick model.

3 Results

3.1 Dating and age-depth model

A total of 13 radiocarbon dates were analysed from core VAR1, four of which were post-bomb dates extending to a depth of 19.5 cm. The oldest dated section of the core produced a date of >1000 years. Radiocarbon and calibrated dates are shown in Table 1.

Table 1 Uncalibrated and calibrated radiocarbon dates from core VAR1, including depths and materials used for dating.

PMC: Percent modern carbon.

Figure 4 Comparison between ²¹⁰Pb activity concentrations throughout core VAR1 and cement production. Values interpolated from a GAM fitted to annual cement production data from Kunda cement factory (Fig. 2; Trumm et al., 2010) at the same age frequency as the VAR1 age depth model (Fig. 5) using a GAMM. Age ranges are based upon the median ages from the age-depth model.

The ²¹⁰Pb activity profile throughout core VAR1 (Fig. 4) declined throughout the core, although not in a typical monotonic pattern, with a substantial increase from 12.5 to 19.5 cm, indicating that the rate of ²¹⁰Pb accumulation was not constant throughout the core. This spike is due to the enrichment of ²¹⁰Pb from fly-ash fallout from the cement factory and other sources (Vaasma et al., 2014, 2017). Comparison of the ²¹⁰Pb activities and cement production rates applied to the age-depth model support this (Fig. 4) with a significant correlation (τ_b=0.487, $p=<\mathrm{0.0001}$), indicating that ²¹⁰Pb accumulation in the core cannot be solely attributed to precipitation. From 32 cm core depth, ²¹⁰Pb activities achieve low levels, although they do not reach the background activity. Unsupported ²¹⁰Pb was calculated using linear regression of the last 5 samples, showing that the samples assumed background activities following the final measured samples.

The ^239+240Pu activity profile reflects the history of atmospheric deposition at the site (Fig. 5). A clear peak at 28.5 cm corresponds to the onset of nuclear testing in 1945, followed by a second, larger peak at 22.5 cm, likely associated with the peak in global fallout from bomb testing in 1963, before the signing of the Partial Test Ban Treaty (Cwanek et al., 2021). Smaller peaks at 16.5 cm suggest possibly the 1986 Chernobyl fallout (Ketterer et al., 2004), and at 5.5 cm potentially originate from the 2011 Fukushima disaster (Bossew, 2013).

Figure 5 ^239+240Pu activity with depth in core VAR1. Selected median dates from the rPlum-derived age-depth model show good correspondence between known fallout events. Each peak may be tentatively related to a known nuclear fallout event, as shown on the right y axis.

The age-depth relationship for core VAR1 (Fig. 6) was calculated by aligning ²¹⁰Pb and ^239+240Pu data with the calibrated radiocarbon dates. Model diagnostics are shown in Fig. S1. Considering the enrichment of ²¹⁰Pb from industrial fallout, the assumption of a constant unsupported ²¹⁰Pb supply was violated. Despite this, anchoring the model with the known 1963 peak in ^239+240Pu activity and adding a constant but uncertain reservoir effect of ca. 15 years for radiocarbon dates improved the model alignment with the peaks in ^239+240Pu. A radiocarbon reservoir effect is possible in Sphagnum peatlands, due to the recycling of “old” gaseous carbon by mycorrhizal fungi associated with ericaceous plants near the peat surface. However, this does not usually occur when individual plant remains are dated, as was the case for our study (Piotrowska et al., 2011). A more likely explanation could be the uptake of the uptake of old carbon via dissolved carbonates', particularly in the more recent samples (Madeja and Latowski, 2008) where old carbon may derive from the buried cement dust.

The resulting age-depth model (Fig. 6) indicates that peat accumulation rates were stable from the base of the core until around 1255 cal CE, after which accumulation slowed to around 0.22 yr cm⁻¹, remaining low until around 1940 cal CE, where there was a sharp acceleration to ca. 2 yr cm⁻¹.

Figure 6 Age-depth model for core VAR1, including calibrated radiocarbon dates and uncertainty (green) ²¹⁰Pb activity (Bq kg⁻¹) in blue and the calendar date for the 1963 nuclear treaty peak (orange). Uncertainties for the model are shown as the shaded area.

3.2 Peat physical and chemical properties

Figure 7 shows that activities of elements associated with clinker dust pollution began to increase around ca. 1873 cal CE, marking the beginning of cement production in Kunda. There is evidence of increasing lithogenic dust deposition throughout the core, beginning as early as the mid-12th century. There is a clear negative correlation between the counts of lithogenic elements associated with cement dust (Ti, Pb, Ca, K, Fe, Cr) and organic matter content, most noticeably between ca. 1942–2006 cal CE. Organic matter content falls sharply until ca. 1988 cal CE, suggesting the accumulation of cement dust occurred during this period, mirroring the results of studies in NE Estonia lakes impacted by oil shale powerplants (Varvas and Punning, 1993; Koff et al., 2016).

The aCAR for the Varudi core is illustrated in Fig. 7, which also shows the trends in LOI %, bulk density, and carbon accumulation throughout the core. Over the past millennium, the average aCAR for the entire core is 148.7 g OC m² yr⁻¹. This is especially high in the upper section of the core, where carbon accumulation rates peak at 342.2±231.7 g OC m² yr⁻¹. For most of the record (ca. 1045 to 1910 cal CE), the mean aCAR is 72.9±28.9 g OC m² yr⁻¹, aligning more closely with average values reported for Northern peatlands (Roulet et al., 2007). The rate of carbon accumulation is relatively stable from the base of the core until ca. 1280 cal CE, with aCAR averaging 94.5±19.8 g OC m² yr⁻¹. After this point, aCAR falls to 51.3±19.7 g OC m² yr⁻¹, remaining low until ca. 1840 cal CE. After this date, aCAR increases significantly, reaching an average of 256.6±221.2 g OC m² yr⁻¹.

Figure 7 Comparison of selected elements detected by the µXRF core scanner, as indicated by the PCA analysis, for the detection of dust fallout from the Kunda cement factory. In addition to peat physical property parameters: Loss on ignition (LOI %), bulk density and apparent carbon accumulation rates (aCAR). µXRF Data are presented as counts per second and normalised by dividing the sum of incoherent (inc) and coherent (coh) activities. The shaded area represents the section of core where most of the cement dust is concentrated. The uppermost aCAR sample was removed from the figure to aid interpretation.

The broken stick model shows that the two first components together explain a significant proportion (48.3 %) of the total variance in the peat's chemical composition (PCA1: 25.5 %, PCA2: 22.8 %). These components highlighted elements associated with clinker dust deposition from the Kunda cement factory and are significantly negatively correlated with peat LOI content (Fig. 8).

Figure 8 Biplot of first and second principal components showing factor loadings of individual elements scanned across core VAR1. Positively correlated variables point to the same side of the plot, while negatively correlated elements point to the opposite sides. Peat organic matter % (Derived from Loss on Ignition) content is shown in blue. The colour of the lines represents the sum factor loading of each variable for both axes, representing how strongly each variable contributes to the principal component.

3.3 Palaeoecological reconstructions

3.3.1 Plant Macrofossils

The results of plant macrofossil analysis are illustrated in Fig. 9. The plant macrofossil data show three major phases of vegetation change, corresponding broadly with shifts observed in the testate amoeba and palynological records (Figs. 10, 11 and 12).

Figure 9 Plant macrofossil stratigraphic diagram illustrating the changes in the botanical composition throughout core VAR1. Note the mixed data types used, percentages and total counts.

From ca. 1000 cal CE, the site is dominated by Sphagnum mosses, particularly those of the subgenus Sphagnum, although by the beginning of the 12th century, Sphagnum section Acutifolia becomes the most abundant. By ca. 1450 cal CE, there was an increasing abundance of monocots, likely from Eriophorum species and Cyperaceae, owing to the presence of spindles and fruits identified to these species. During this period, Sphagnum gradually declined, eventually disappearing entirely by the start of the next phase.

The latter portion of the record, starting ca. 1970 cal CE, is characterised by a shift towards more “woody” vegetation. Shrub-type taxa increase initially, followed by a rise in ligneous remains, particularly those of Pinus sylvestris. Mycorrhizal roots, bark fragments, and pine needles become significant components of the peat's botanical composition in the upper section. Betula is represented in the uppermost samples by its characteristic catkin scales and fruits, along with a Betula nana leaf recovered from the surface sample. This late phase is also characterised by a high percentage of unidentified organic material. At the top of the core, remains of the brown moss Tomentypnum nitens, characteristic of calcareous fens, were identified (Hájek et al., 2021).

3.3.2 Testate amoeba sub-fossil communities

Throughout core VAR1, a total of 105 distinct testate amoeba taxa were identified. The results of the testate amoeba analysis and reconstructions are presented in Fig. 10. Three distinct zones were identified throughout the record by CONISS.

The first zone, extending from the base of the core to 44.5 cm (ca. 1600 cal CE) is characterised by a mix of proteinaceous and mixotrophic species, including Archerella flavum and Hyalosphenia papilio. Other abundant species in this zone include Difflugia pulex, Phryganella acropodia, Assulina muscorum and Assulina seminulum. The second zone, ending at 19.5 cm (ca. 1980 cal CE), shows a replacement of the dominant species by Galeripora discoides, Amphitrema stenostoma, and xenosomic species such as Cyclopyxis arcelloiodes. By 35.5 cm (ca. 1821 cal CE), mixotrophic species have almost disappeared, replaced by Cyclopyxis eurystoma, Centropyxis aculeata, Difflugia ampullula, Difflugia brevicolla and Pseudodifflugia gracilis. The final zone, at the top of the core, is dominated by Centropyxis elongata, Centropyxis sylvatica, Plagiopyxis sp., Euglypha rotunda type, and Euglypha leavis.

Figure 10 Stratigraphic plot showing changes in the relative abundance (%) of testate amoebae taxa identified throughout core VAR1, as well as reconstructed peat pore water pH and water table depths (WTD) and uncertainties. The results of CONISS are illustrated to the right of the figure, with horizontal red lines representing the zone boundaries defined by CONISS. Only species with maximum abundancies above 5 % are illustrated. The full dataset is available in the Supplement.

The relative abundance of testate amoeba taxa was used to reconstruct changes in water table depth and peat pore water pH over time. Both reconstructions exhibit notable trends throughout the past ca. 1000 years. From the base of the core until ca. 1330 cal CE, water table depth was relatively constant, averaging 7.0±2.9 cm. Conditions became progressively wetter after this time, reaching a minimum of −2.5 cm by ca. 1540 cal CE, possibly indicating a period of open water at the coring location. Water table depth began to decrease gradually, accelerating after ca. 1960 cal CE and reaching a maximum by ca. 2012 cal CE.

Reconstructed pH values were relatively stable throughout much of the record (4.7±0.6). Around ca. 1850 cal CE, pH began to rise, reaching a maximum of 7.2±0.6 by ca. 1985 cal CE, remaining high until the uppermost samples, where pH returned to near pre-disturbance levels (4.8±0.5), this change occurring after ca. 1998 cal CE.

3.3.3 Pollen, spores, and non-pollen palynomorphs

The results of the palynological analysis are illustrated in Fig. 11 (Pollen and spores) and Fig. 12 (Non-Pollen Palynomorphs – NPPs). The pollen sequence from the core VAR1 is dominated throughout by arboreal taxa, particularly Pinus sylvestris and Betula sp., with lesser contributions from shrubs such as Calluna vulgaris and Vaccinium type. The high dominance of these arboreal taxa and a low sampling resolution make detecting subtle shifts in human activity in the palynological record difficult to infer (Favre et al., 2008). Other palynological studies from Estonia also report low variation in pollen assemblages and dominance of P. sylvestris and Betula across different regions, even during periods of significant land-use changes (e.g., Poska et al., 2004; Veski et al., 2005; Łuców et al., 2022).

Figure 11 Stratigraphic plot illustrating the relative abundance of selected pollen taxa from core VAR1, as well as microcharcoal (expressed as the total number of counted fragments). Relative abundances of arboreal taxa are shaded green, disturbance indicators in pink, crops in orange, herbaceous taxa in brown and shrubs in teal. Grey shading represents 5× exaggeration. Only taxa with a minimum abundance of 2 % are shown.

Figure 12 Stratigraphic plot illustrating the abundance of selected non-pollen palynomorphs from core VAR1, as well as spheroidal carbonaceous particles (SCPs) and Pinus stomata, both expressed as counts. Relative abundances of spores of aquatic taxa are shaded blue, algae: light green, fungi: red, invertebrate remains grey and testate amoebae in dark green. Grey shading represents 5× exaggeration. Only taxa with a minimum abundance of 2 % are shown.

The oldest section of the core (ca. 1000–1250 cal CE) provides little evidence of anthropogenic disturbance, indicating a largely forested landscape. This period is dominated by Pinus sylvestris, Betula sp. and Picea abies, suggesting a stable, predominantly woodland environment. After ca. 1250 cal CE, there is a slight increase in pollen from grasses (Poaceae) from around 3 % to ca. 8 % of the total assemblage until ca. 1680 cal CE, as well as the appearance in low concentrations of plants indicative of disturbance, such as Plantago lanceolata, Chenopodiaceae and Ranunculus acris type and a general increase in Secale pollen which together may indicate agricultural activity, although wild populations of Secale are believed to have existed in Estonia (Veski, 1998; Poska, 2003). Despite this, arboreal pollen remains the dominant component of the record, with declines mainly affecting Pinus sylvestris (49 % at ca. 1310 cal CE to 34 % by ca. 1680 cal CE) and Picea abies (16 % at ca. 1220 cal CE to 5 % by ca. 1680 cal CE). Tree species associated with early succession, such as Betula, remain largely unchanged, while Alnus increases during this period from ca. 6 % at the start of the record to ca. 10 % by ca. 1680 cal CE. By ca. 1680 cal CE, arboreal pollen has declined to its lowest relative abundance (around 82 % of the total pollen sum). This decline is driven mainly by the changes in arboreal species, particularly declines in P. sylvestris, which becomes progressively less common through the core, a trend that continues until ca. 1965 cal CE. By ca. 1680 cal CE a peak in cereals occurs, including Hordeum, although these are only present in low concentrations. The pollen record is stable after this period, with only slight (<1 %) increases in Calluna vulgaris and Vaccinium type taxa occurring. These changes correspond with a rapid increase in microcharcoal accumulation rates, rising from 311 to 1027 particles cm⁻³ yr⁻¹ following the opening of the cement factory, peaking at 9837 particles cm⁻³ yr⁻¹ by ca. 1970 cal CE. After this, Pinus sylvestris begins to recover, but the relative abundance of Betula sp. continues to rise. By ca. 2010 cal CE, Pinus sylvestris becomes dominant once again, while Betula declines throughout the most recent samples, and microcharcoal accumulation reaches a maximum at ca. 2010 cal CE, of 172 237 particles cm⁻³ yr⁻¹.

The NPP record provides a more detailed picture of changes in the local environment than the pollen record. In the earliest portion of the record (ca. 1000–1220 cal CE), fungi such as HdV-90 and HdV-13 (cf. Entophlyctis lobata) are common, and HdV-27 (Bryophytomyces sphagni) is present throughout. These fungi are typically associated with oligotrophic and ombrotrophic conditions, although HdV-90 can also thrive in more minerotrophic or poor-fen environments (van Geel, 1978; Kuhry, 1985).

Between ca. 1220–1460 cal CE, HdV-13 increases in abundance. Around 1850 cal CE, taxa associated with the earlier section of the record start to decline, being replaced by the microalgae Botryococcus taxa, indicative of aquatic conditions in addition to several fungi including HdV-55A (Sordaria type), HdV-112 (Cercophora type). The most notable species in the upper section of the core during this period is HdV-201 (cf. Xylomyces chlamydosporis), a wood-inhabiting fungus linked to freshwater environments or pool vegetation (Goh et al., 1997; Kuhry, 1997). This species is especially abundant ca. 2006 cal CE, comprising around 34 % of the assemblage. After ca. 1950 cal CE, the diversity of NPPs declines. A key species identified in this section is Botryococcus braunii, a green alga typically found in environments with high levels of inorganic phosphorus, which peaks around 1985 cal CE (Órpez et al., 2009) suggesting a shift toward more nutrient-enriched conditions. Overall, the record suggests a transition from nutrient-poor, Sphagnum-dominated peat towards an increasingly nutrient-enriched system, followed by a change towards decomposers of ligneous material.

3.3.4 Rate of change analysis

The principal response curves (PrC) explained 79 % of the variance in the sub-fossil testate amoeba data from Varudi and are illustrated in Fig. 13. The poor fit of GAM models to the PrCs prompted us to use adaptive splines with a GAM, which are recommended for data exhibiting abrupt changes (Burge et al., 2023). However, adaptive spline GAMs cannot yet be used within the GAM framework as described by Simpson (2018). One outlier was removed due to the difference in testate amoebae composition owing to a high proportion of Cryptodifflugia angusta (33.5 cm).

The PrC revealed two separate periods of change at Varudi, as well as highlighting periods of rapid change associated with the onset of cement production in Kunda in 1871, and an earlier change indicated from the palynological record to coincide with increasing human activity in the landscape (Fig. 11). Testate amoeba community composition was relatively consistent from ca. 1040 cal CE until ca. 1330 cal CE and again between ca. 1645 to ca. 1886 cal CE. Following ca. 1998 cal CE, the direction of change in the PrC curve reverses, coinciding chronologically with the installation of filters in the cement factory in 1996. By the end of the record, testate amoeba community composition was characterised by the return of species such as Cyclopyxis arcelloides, Cryptodifflugia oviformis, and Assulina muscorum, which were not abundant since before ca. 1970 cal CE.

Figure 13 Changes in principal response curve scores derived from sub-fossil testate amoebae data, modelled using an adaptive spline generalised additive model GAM (black line). Solid red lines indicate periods of rapid change, identified where the modelled confidence interval of the slope does not include zero. Vertical dashed lines denote key chronological events potentially influencing environmental conditions at Varudi: The onset of human activity recorded in the pollen data (ca. 1250 cal CE, orange), the beginning of cement production in Kunda (1871, blue) and the installation of pollution control filters at the cement factory (1996, black).

4 Discussion

4.1 Palaeoecological assessment of ecosystem and functional changes before, during, and after intensive alkaline dust pollution

The palaeoecological record of the Varudi peatland reveals three distinct phases occurring over the past 1000 years: an ombrotrophic bog phase (ca. 1000–1250 cal CE), a poor fen phase (ca. 1250–1570 cal CE), and a wetland forest phase that continues to the present day. Most contemporary research on peatland recovery has been conducted on historically drained, mined, or agricultural sites (Roderfeld, 1993; Graf et al., 2008; Wagner et al., 2008; Paal et al., 2010). In contrast to these studies, an advantage of our palaeoecological approach is that it provides a unique pre-disturbance context for understanding the effects of atmospheric pollution upon this site.

The different phases of development exhibited by the site in contrast with the different periods of anthropogenic activity and disturbance are summarised in Fig. 14. During the early ombrotrophic phase, low palynological richness, minimal microcharcoal counts and abundant peat organic matter indicate a largely forested landscape with limited human impact near the site, consistent with other palaeoecological records throughout Estonia during this period (Veski et al., 2005). The dominance of Sphagnum mosses and mixotrophic testate amoeba species such as Archerella flavum and Hyalosphenia papilio reflect the relatively undisturbed and productive nature of the bog during this phase (Marcisz et al., 2020; Łuców et al., 2022). During this period, water table depths were high, and pH levels were low, consistent with natural conditions for ombrotrophic bogs in the Northern Hemisphere (Warner and Asada, 2006).

Figure 14 Summary of key trends in the Varudi peatland record. Shown are changes in reconstructed pH and water table depth, percentages of non-arboreal pollen (serving as a proxy for forest clearance related to agricultural activity), apparent carbon accumulation rates (aCAR), periods of significant change (as identified by principal response curves with Generalised Additive Models, GAMs), and a general timeline of environmental shifts. Dashed horizontal lines represent key historical events relevant to the Varudi record: the yellow line marks the approximate onset of human activity ca. 1250 cal CE, possibly linked to the Livonian Crusades; the blue line indicates the establishment of the Kunda cement factory in 1871; the black line denotes the installation of effective filters at the Kunda cement factory; and the grey-shaded area highlights the section of the core with substantial cement dust, identified through LOI and µXRF core scanning.

The beginning of agricultural activity around ca. 1250 cal CE is broadly coincident with the German-Danish crusaders occupation of Estonia starting in 1227 CE (Veski et al., 2005). From ca. 1250 to 1570 cal CE, increased Poaceae, Secale cereale, and field weed pollen indicates a shift toward more open, meadow-like conditions (Poska and Saarse, 2002; Niinemets and Saarse, 2009), which may relate to the environmental effects resulting from economic changes resulting from the Teutonic order crusades, similar to those identified by Brown and Pluskowski (2011) in a peatland record from northern Poland.

Declines in Pinus sylvestris and Picea abies, alongside increases in early successional species such as Betula and Alnus, and minute increases in microcharcoal accumulation rates (Fig. 11) suggest land clearance using slash-and-burn methods in the surrounding landscape (Jääts et al., 2010). While arboreal pollen remained dominant, indicating continued forest cover around the Varudi peatland, mineral soil enrichment of the peat suggests overland flow or aeolian deposition resulting from land clearance. Similar patterns of early human-induced soil erosion have been recorded in peatland records from central Estonia (Heinsalu and Veski, 2010). This is reflected by decreasing organic matter content and increasing concentrations of lithogenic elements such as Ca and Ti, patterns like those observed in Estonian lake sediment records following land clearance events (Koff et al., 2016; Vandel and Vaasma, 2018). Ti, a conservative lithogenic element, indicates soil erosion from the surrounding catchment (Boyle, 2001; Boës et al., 2011).

By ca. 1570 cal CE, the site had transitioned into a wetter, more minerotrophic fen, indicated by the gradual replacement of Sphagnum by monocots such as Carex sp. and Eriophorum sp., increases in Sphagnum sect. Cuspidata and subsequent declines in Sphagnum sect. Acutifolia and Sphagnum in the plant macrofossil record (Fig. 9). This shift towards wetter conditions is also supported by the testate amoebae reconstruction, possibly indicating permanent or seasonal standing water at this time (Fig. 10), and the presence of Scheuchzeria palustris pollen between ca. 1500–1800 cal CE (Fig. 11), a flood-tolerant species (Tallis and Birks, 1965). This interpretation is consistent with observations from contemporary bog pools, where Sphagnum is often rare or absent in shallow or temporary water bodies on a peatland's surface, while monocots such as Eriophorum species are abundant, as exemplified by the British National Vegetation Classification (NCV) bog pool community classification (M3) (Rodwell, 1998). Nutrient enrichment due to soil erosion resulting from forest clearance and agricultural activities in the landscape at this time likely contributed to this shift (Hölzer and Hölzer, 1998).

The establishment of the Kunda cement factory in 1871 brought significant ecological changes to the site. Substantial increases in Ti, Ca, Fe, and other elements associated with cement dust fallout occurred along with rapidly increasing pH and declining organic matter content, indicating the chemical influence of cement dust pollution upon the peat (Figs. 7 and 8). This pollution appears to have driven rapid changes in testate amoeba communities by ca. 1890 cal CE, as seen in the PrC record (Figs. 13 and 14). These results are consistent with previous palaeoecological studies on alkaline pollution impacts in Estonian peatlands and lakes (Varvas and Punning, 1993; Koff et al., 2016; Vandel and Vaasma, 2018; Vellak et al., 2014). The enrichment of the peat by cement dust is reported to have had a fertiliser-like effect on ligneous and ericaceous species, increasing tree cover, size, and stand density within the industrial region of Estonia at this time (Pensa et al., 2004; Ots et al., 2011).

The concurrent hydrological changes experienced by the site, characterised by persistent water table lowering after ca. 1960, were likely driven by drainage and enhanced transpiration rates following tree and shrub expansion onto the site (Stelling et al., 2023). The simultaneous compounding effect of hydrological change upon the site confounds interpretation of the effects of alkaline pollution at Varudi peatland, as it is likely that the changes that occurred at this time were driven by the combined effects of both disturbances.

From 1955 CE onward, cement production increased, reflected in peaks of Ca, Ti, and K, alongside rising bromine (Br) levels, coincident with increased microcharcoal accumulation rates (increasing from ca. 11 to 58 particles yr⁻¹. The halogen Br is often linked to changes in storminess in peat palaeoecological studies due to being commonly sourced from marine aerosols (e.g., Shotyk, 1997; Turner et al., 2014). Here, elevated Br levels correlate strongly with cement dust indicators, suggesting that increased Br in the core corresponds with enhanced plant decay and oxidation of abiotic organic matter, resulting in the formation of methyl bromide (Lee-Taylor and Holland, 2000; Keppler et al., 2000). The increased presence of fungi in addition to microalgae species Botryococcus braunii in the palynological and plant macrofossil records further suggests that the soil decomposition rate increased during this period, coincident with periods of increased seasonal wetness (Barthelmes et al., 2012; Defrenne et al., 2023; Buttler et al., 2023; Thormann et al., 2003) (Fig. 12).

Following the restoration of Estonia's independence in 1991, more rigorous environmental regulations resulted in a reduction of industrial emissions (Kask et al., 2008). The effect of this decline is evident in our palaeoecological record, with declining pollution markers and pH beginning to return towards pre-disturbance levels after 1996 (Figs. 7 and 10), as well as the apparent (albeit limited) recovery of testate amoeba communities shown by the PrC curve (Fig. 13). Despite these improvements, the site has remained densely forested and enriched, with a persistently lowered water table. We find no evidence of botanical succession toward pre-disturbance conditions in the plant macrofossil record (Fig. 9). The presence of the calciphilous moss Tomentypnum nitens at the core surface likely reflects the legacy of alkaline pollution at the site (Maimer et al., 1992; Hájek et al., 2006; Apolinarska et al., 2024). These cement deposits may still influence surface water chemistry and vegetation through the upward movement of water enriched by the buried dust. Nonetheless, Sphagnum has re-established at the site since its disappearance in the 1960s, suggesting a degree of recovery. Identification and monitoring of the present-day flora could provide additional insights into the extent of this recovery.

4.2 Ecosystem Recovery After Reduction in Atmospheric Alkaline Pollution

Our study found evidence of recovery in PrC of the sub-fossil testate amoeba assemblages (Fig. 13). The community composition of these environmentally sensitive proxies showed signs of turnover nearly 30 years after installing filters at the Kunda cement factory. However, the extent of recovery indicated by the testate amoebae is limited, as the community composition in the top sediment layer more closely resembles conditions from ca. 1977 cal CE, the peak period of cement production at Varudi (Figs. 3 and 13).

We used adaptive spline GAMs to analyse our data, which while suitable for datasets exhibiting abrupt changes such as the example from Varudi (Burge et al., 2023), they cannot yet be incorporated into the Generalized Additive Mixed Model (GAMM) framework described by Simpson (2018) and are more susceptible to issues such as temporal autocorrelation and boundary effects (Simpson, 2018; Burge et al., 2023). While the model best fit our data, more evidence is needed to draw conclusive interpretations about recovery than these results alone. More recently accumulated samples may be more influenced by temporal autocorrelation (Burge et al., 2023). Some testate amoeba species may infiltrate the sub-fossil community composition downcore, making near-surface communities appear more similar (Liu et al., 2024). Another limitation of this study is that our results are based upon reconstructions taken from a single core and as such do not account for the spatial heterogeneity characteristic of peatlands. Although palaeoecological trends tend to be well replicated across multiple cores within a given site (e.g., Barber et al., 1999; Hendon et al., 2001), it may be that the trends reconstructed from VAR1 represent conditions specific to the sampling site, rather than the whole site (Bacon et al., 2017). Paal et al. (2010) examined peatland vegetation communities in eleven bogs in northeastern Estonia at varying distances from sources of atmospheric alkaline pollution, finding evidence of Sphagnum re-establishment and the recovery of other bog-specific plant species at several sites as surface waters re-equilibrated to pre-pollution conditions, particularly within microforms that are less affected by contaminated groundwater. Therefore, recovery at Varudi may be spatially heterogenous, with some areas recovering more quickly depending upon their contact with polluted peat layers.

Studies investigating the recovery of vegetation communities following alkaline pollution show that ecosystem recovery may be limited due to the present, degraded condition of the ecosystem. For example, Vellak et al. (2014) studied bryophyte recovery following reductions in atmospheric fly-ash pollution across Estonia and northwestern Russian sites. They found that bryophyte species growing in more heavily affected sites tended to be adapted for growing in low light conditions, due to competing with the larger and more dense vascular plants that encroached on these sites, in response to the more enriched and alkaline conditions for which they are better adapted than bog-specific vegetation (Pärtel et al., 2004; Zvereva et al., 2008a, 2008b). Increased nutrient and litter supply and root exudates, coupled with the faster growth of vascular plants and trees may further hinder the reestablishment of bog and fen communities, delaying or preventing recovery (Konings et al., 2019). Gunnarsson et al. (2002) demonstrated that nutrient enrichment, often accompanied by higher pH, can give vascular plants a competitive advantage over bog-specific vegetation (Limpens et al., 2003; Dieleman et al., 2015).

4.3 How has the pollution impacted apparent rates of carbon accumulation at the site?

Each of the phases in ecosystem development in the palaeoecological record corresponds with substantial changes in the apparent carbon accumulation rate (aCAR) throughout the Varudi record (Fig. 7). Between ca. 1000 to 1250 cal CE, when the site was in a relatively pristine, undisturbed condition, aCAR rates were higher (32.0±10.7 g C m² yr⁻¹) than the average for northern hemisphere ombrotrophic peatlands (23±2 g C m² yr⁻¹) (Korhola et al., 1995; Yu et al., 2009). Throughout the transition from bog to poor fen, from ca. 1360 to 1570 cal CE, average aCAR rates most likely initially increased (38.5±34.3 g C m² yr⁻¹). However, they were highly variable, likely reflecting the incorporation of vascular plant roots from the overlying forested phase, introducing younger carbon into deeper sediments. Between ca. 1630 and 1800 cal CE, rates fell as low as 19.0±5.4 g C m² yr⁻¹, aligning with reported values for Finnish fens (Korhola et al., 1995). Overall, minerotrophic fens typically exhibit lower carbon accumulation rates than ombrotrophic bogs (Loisel and Bunsen, 2020; Yang et al., 2023).

Following the succession from poor fen to forested fen from ca. 1871 cal CE to the present, there is an apparent increase in aCAR, especially in the most recently accumulated peat (280.7±211.3 g C m² yr⁻¹). Due to the incomplete decomposition of labile organic matter at the peat surface, care must be taken when interpreting recently accumulated carbon from peat core records (Young et al., 2019, 2021). This increase is likely due to high litter deposition from trees and shrubs, and the rapid decline in aCAR downcore suggests that this labile material is rapidly cycled back into the atmosphere rather than sequestered in the soil. Our results suggest that the current vegetation composition, dominated by trees and shrubs, is less effective for long-term carbon sequestration compared to earlier phases. While the water table remains reduced, the site shows limited potential for recovery to its original ombrotrophic condition due to ongoing drainage. At the same time, nutrient inputs and evapotranspiration from the overlying trees compound this limiting factor.

4.4 Indicators of critical thresholds to assess peatland condition and recovery due to alkaline pollution

By comparing our pH reconstruction with significant successional shifts at the site, we may infer the thresholds at which transitions in the steady state of the ecosystem occurred in response to changes in alkalinity at Varudi bog. Defining tipping points that may be broadly applied to other sites allows for predicting ecosystem shifts in response to global change, allowing for prediction of such shifts in the future (Munson et al., 2018). To our knowledge, our study represents the first step towards defining thresholds for tipping points in peatland ecosystems in response to changes in alkalinity.

Around ca. 1910 cal CE, the reconstructed pH for Varudi showed an increasing trend from the previous, relatively stable conditions, reaching a maximum of 7.2±0.6 by ca. 1985 cal CE (Fig. 10). This value aligns with previous pH measurements for forest litter in the zone surrounding the Kunda cement factory of 7.1–7.4, supporting the reliability of our reconstruction (Paal and Degtjarenko, 2015). We find that the transition from an ombrotrophic bog to fen-like conditions during the mid 13th century, driven by mineral soil enrichment resulting from land clearance for agriculture in the surrounding landscape, was associated with an increase of average pH levels from ca. pH 3.9 to pH 4.1 which represents (due to the logarithmic scale of pH) more than 1.5 times increase in alkalinity. In contrast, the subsequent shift from fen-like to forested conditions corresponded with a small average increase in pH to ca. 5.7, over 40 times more alkaline than the previous state and more than 60 times more alkaline than during the bog phase. Our results suggest that while substantial increases in pH can drive rapid and substantial ecosystem regime shifts, relatively small increases in pH (increases in pH of 0.2–0.3) may also drive critical ecosystem transitions that may be more gradual but can have equally as significant and long-lasting impacts upon peatland ecosystem functioning, such as the capacity for carbon accumulation. However, we stress that the threshold values defined here remain uncertain, owing to the relatively significant uncertainties associated with transfer function models in general (Amesbury et al., 2016). Higher pH values measured from bog pools at the site in 1996–1997 of 7.6–8.5 are slightly higher than our estimates for this period (6.1–6.6), and by 2022 pH had only fallen to 6.3–6.8 (Pärtma, 2023) although our reconstruction indicates pH had fallen to 4.8 by this time, suggesting that our reconstruction, despite having a strong predictive power, underestimates pH levels experienced by the site. Alternatively, this may simply reflect local variation in pH within Varudi, as the pH measurements were taken from bog pools in the NW margin of the bog, not adjacent to the coring location. This may warrant future re-evaluation of our data using transfer functions specifically developed for reconstructing pH in polluted peatlands. In addition, we cannot rule out the possibility that other impacts of cement deposition (e.g., reduced photosynthesis rates due to dust deposition) and, significantly – drainage; may have also played a role in driving local environmental changes at the site.

Contemporary field and lab experiments that manipulate the alkalinity of peatland soils directly may provide more precise values for ecosystem thresholds in response to alkaline pollution. However, to date few studies have done this for alkalinisation: one example is Kang et al. (2018) who conducted a series of field and laboratory pH manipulation experiments across seven peatlands in the UK, Japan, Indonesia and South Korea, finding that increases in pH resulted in significant changes in microbial composition, resulting in increased phenol oxidase activity and enhanced DOC releases. Another is the long-term (active since 2002) Whim Bog experiment near Edinburg, Scotland (Sheppard et al., 2011; Levy et al., 2019), where the effects of enrichment upon a peatland by gaseous ammonia (NH₃) and wet-deposited ammonium (NH${}_{\mathrm{4}}^{+}$) and nitrate (NO${}_{\mathrm{3}}^{-}$) are compared. Analysis by Sutton et al. (2020) indicates that vegetation declined three times more quickly when exposed to gaseous ammonia and three times for ammonium than for a similar dose of nitrate.

A significant limitation of these studies in defining tipping points is that the long timescale necessary for critical transitions to occur in some cases typically exceeds the lifespan of most observational and experimental studies (Taranu et al., 2018; Lamentowicz et al., 2019). Therefore, further palaeoecological work across a more extensive range of sites or with multiple cores from within a single site may also advance our understanding of ecosystem tipping points in response to alkalinisation. Furthermore, disentangling the effects of alkalinisation from those of drainage is challenging, as both factors may influence plant community composition, carbon accumulation rates, overall ecosystem functioning (Word et al., 2022) and possibly resilience to alkalinisation in peatlands. Since significant drainage and cement production appear to have begun at roughly the same period in the site's history, some of the ecological variability attributed here to cement dust pollution may instead be the result of drainage impacts. However, it is notable that despite reconstructed pH values decreasing following the cessation of significant emissions from Kunda, water table depths have not recovered. This indicates that drainage is unlikely to have influenced the apparent modest ecosystem recovery observed at Varudi in recent decades.

5 Conclusion

This study is the first to investigate the long-term impacts of alkaline emissions on a peatland over centennial timescales. It establishes the first attempt to define threshold indicator values for ecosystem tipping points in response to alkalinisation. Our findings demonstrate that alkaline pollution, along with later drainage, has had a profound influence upon ecosystem development at Varudi for more than 750 years, with ecosystem succession following both low-level, sustained mineral soil enrichment due to agricultural activities and intensive fly-ash fallout sustained over 160 years. We find that an increase in pH of 0.2 to 0.3 (corresponding to, approximately, a two to threefold increase in alkalinity) is sufficient to trigger a critical ecosystem shift, while an increase of 1.5 pH units (a $\sim \mathrm{40}\times$ increase in alkalinity) will cause rapid ecosystem transition, resulting in the loss of Sphagnum mosses and changes in ecosystem functioning. This can lead to a long-term decline in carbon accumulation over long timescales, and such changes may be slow to recover or permanent, even if the point source of pollution is eliminated. Our results, while insightful, have limitations that underscore the need for additional experimental and palaeoecological research to assess peatland responses and resilience to alkalinisation across a range of spatial and temporal scales. This would allow a better understanding of the timescales required for peatland recovery and how these ecosystem transitions influence greenhouse gas dynamics from affected sites. Alkalinisation poses a growing threat to peatlands worldwide and is a developing challenge for the 21st century (Sutton et al., 2020). As carbon-rich ecosystems, peatlands are important for regulating future atmospheric greenhouse gas concentrations. Understanding how peatlands will respond to future alkalisation is essential for predicting the role of peatlands in climate change mitigation.