Introduction

Environmental site conditions of vegetation depend on multiple driving factors. In landscape ecology, these drivers are categorised into a hierarchy of decreasing dominance in which climate comes first, followed by relief, hydrology, minerals, and organics (Kemmers et al. 2011). Nearby areas will experience a similar climate, but the other drivers may be very different, depending on the landscape setting and their developments over time. In deltas, gradients over small distances may cause considerable differences in site conditions, leading to different habitats where plants can colonize and succeed each other (Beeftink et al. 1985).

Since deltas are attractive to people who settled for agricultural use and fisheries, deltas have also undergone human interference to control floods: dikes were built for protection, consolidating reclamation. These dikes prevented most floods but not all. For example, in Zeeland, the southwestern province of the Netherlands, reclaimed areas were flooded due to a heavy storm surge in February 1953, causing 1836 casualties and damaging buildings and infrastructure (Nienhuis 2008). To prevent new floods, most of the estuaries in the Dutch delta were closed by high sea-walls constructed in 1958–1986 (Bannink et al. 1984), nowadays known as the Dutch Delta Works.

Traditionally, Dutch reclaimed areas were put to agricultural use, but the former estuaries’ reclaimed areas were left unused. These areas became former tidal wetlands and were designated as nature areas. They all have their specific year of reclamation and management common at that time (Slim and Oosterveld 1985), making them exciting objects for studying succession. Initially, management was trial and error based in these reserves, but monitoring programmes were conducted to understand hydrology-driven processes such as freshening and acidification. Unfortunately, monitoring stopped after some years in most areas. A notable exception was the wetland Slikken van Flakkee in the Grevelingen, where Saeijs (1973), building on the work of Beeftink (1965), initiated annual vegetation monitoring. This monitoring was continued by De Jong and Van der Pluijm (2003) and Van Haperen (2009), resulting in a unique time series of c. 50 years of vegetation monitoring. Slager et al. (1986) and Slager and Visser (1990) monitored the water and salt movement, and Van Dam (1997) was the first to model these movements hydrologically. The generated research results were regularly evaluated to identify the interrelations, e.g., De Jong and Van der Pluijm (2003) revealed that the salt-fresh gradient is necessary for valuable botanic vegetation. Van Haperen also added some chemical data, allowing him to relate physico-chemical conditions to vegetation (while focusing on Salix and Hippophae rhamnoides thickets).

Joenje and Verhoeven (1993) summarized the hydrological and soil processes in the Dutch former tidal wetlands and related them to vegetation succession. These authors stated: “The patterns of salt and water movements provide the predominant factors for the development of the vegetation.” They mentioned the Slikken van Flakkee as an example that underpinned this statement and referred to Joenje (1978) for the prediction of successional trends related to the diminishing CaCO3 content of the upper soil layer, the lowering of the pH and the accumulation of humus and nutrients. A comparable prediction was made by Van Noordwijk-Puijk et al. (1979) for the Middelplaten in Lake Veere, a former tidal wetland c. 30 km southwest of the Slikken van Flakkee.

It is obvious that geochemical processes drive the mentioned soil processes referred to by the authors mentioned above. Appelo and Postma (2005) pointed out which geochemical processes can be expected in freshening (desalinating) soils, also demonstrating the advantages of inverse modelling with the program PHREEQC (Parkhurst and Appelo 2013). Based on a column experiment with soil from a nearby comparable area, Appelo et al. (1998) modelled these processes and found pyrite oxidation to trigger calcite dissolution, cation and proton exchange, and CO2 sorption. Stuyfzand et al. (2014) also encountered redox reactions, dissolution of calcite and cation exchange as dominant processes during freshening based on their extensive hydrological and hydrogeochemical research for quantifying physico-chemical conditions related to Liparis loeselli on the Veermansplaaat, an island next to the Slikken van Flakkee. They suggested that mixing (of rainwater and salt Lake Grevelingen water) and the mentioned processes led to CaHCO3 as the dominant water-type.

Of course, these processes also occur in tidal wetlands elsewhere. Often sulphur and lime content dictate the occurring processes as they determine if acidic conditions can develop, e.g., in tidal lowlands in Indonesia where acidic conditions developed after reclamation, hampering agricultural use (Suryadi 1996). In Australia, remediation of these conditions by re-establishing marine tidal inundation (Johnston et al. 2011) and adding lime (Johnston et al. 2012) were quantified, using PHREEQC for speciation calculations. Frei (2012) applied PHREEQC more advanced for studying the interaction between hydrology and biogeochemistry within riparian wetlands in Germany by virtual modelling. The inverse modelling with PHREEQC facilitates reconstructing the evolution of groundwater composition and thus understanding and quantifying the occurring geochemical processes (Belkhiri et al. 2010; Karlsen et al. 2012). For vegetation studies, the most effortless assessment of environmental conditions is using plant species as indicators (Ellenberg et al. 1992), but more common is to assess them directly with some geochemical data and to reveal the interrelations with statistical analysis, e.g., by Abdelaal et al. (2018), Sánchez-Higueredo et al. (2020), and Anderson et al. (2022).

To our knowledge, vegetation studies combined with inverse modelling based on complete geochemical data have not been applied yet for former tidal wetlands, and a full understanding and quantification of the occurring processes are still lacking. Therefore, we see a unique opportunity to apply inverse modelling for the Slikken van Flakkee in combination with long-term vegetation monitoring.

Our research aim is to quantify the occurring ecohydrogeochemical processes in the Slikken van Flakkee in response to the tide-stop. By analysing the available vegetation data and determining the present hydrological, geochemical and nature management situation, we expect to identify and quantify these processes. Accordingly, we have three research questions: (1) How did the vegetation develop since the tide-stop? (2) How did the hydrological and geochemical conditions develop? (3) What is the relationship between the vegetation and the physico-chemical conditions in the Slikken van Flakkee? Moreover, we will compare this case study to other published literature about former tidal wetlands and discuss how the geochemistry differs between vegetation types and relates to management.

Such research would help understand natural processes in the Grevelingen relevant for policy-making. Notably, a partial re-establishment of the tide to improve water quality in the lake is being considered for Lake Grevelingen (Postma and Lammers 2014). Ultimately, such partial re-establishment of the tide has to anticipate to sea-level rising and will then affect the Slikken van Flakkee, at least by partial flooding and an upward shift of habitats along the elevation gradient.

We named our multidisciplinary research ecohydrogeochemistry to express ecohydrology’s enhancement with geochemistry.

Site description

The former tidal wetland Slikken van Flakkee is located in the southwest of the Netherlands (Fig. 1a) at the border of Lake Grevelingen (Fig. 1b). The area has a temperate sea climate with a mean annual precipitation surplus of c. 200 mm; mean temperature in January of 4 °C and in July of 17.5 °C (KNMI 2017). The area is bordered at one side by the lake and at the land side by dikes. The area is c. 16 km2 and includes former marshes, mudflats and interconnected creeks. The sediment consists of fine sand, locally with peat and sandy clay underneath at some m below ground surface (bgs) (only next to the former gullies: a central large former gully and smaller ones in the northern and southern part; Fig. 1c).

Fig. 1
figure 1

Topographic maps of a former estuary Grevelingen, the Netherlands, b the Slikken van Flakkee next to Lake Grevelingen, c elevation contours, top of impermeable layer and sub-areas, d monitoring locations. All locations sampled in 2017 are labelled with their PQ number, and those both sampled in 2016 and 2017 are also labelled with their symbol indicating their generalised position (North/South I/South II; High/Middle/Low)

Full size image

In 1971 tidal Grevelingen was blocked from the sea by the Brouwersdam, which led to the falling dry of the upper half of the former intertidal zone of the Slikken van Flakkee. To limit erosion in the new landscape, mixtures of cereals and perennial grasses were sown (Zuidema et al. 1973). The water in Lake Grevelingen is still saline with an aimed Cl− concentration of ≥ 16,000 mg/l (Wetsteijn 2011), and during strong winds and storms, the Slikken van Flakkee is partly flooded with this water (Slager et al. 1986; Slager and Visser 1990).

Since 1971, the groundwater level in the Slikken van Flakkee has been controlled by infiltration, evapotranspiration, flow in the soil directly influenced by its permeability, and seepage. Consequently, variation in groundwater levels and desalinisation depend on the position relative to the gullies (Fig. 1c): relatively stable water levels and completed desalinisation at the former gullies; more dynamic conditions and higher salt contents in the rest of the area.

Nature management in the Slikken van Flakkee is differentiated: in ‘North’, vegetation is developing spontaneously without active nature management; ‘South’ is permanently grazed by cattle and horses in low densities (0.1–0.2 animal per hectare) and is partly mown triennial.

In Online Resource 1, more specific information about the Slikken van Flakkee is presented: geology, soil type, geomorphology, water level of Lake Grevelingen, desalinisation, vegetation monitoring, and policy-making for partial re-establishment of the tide in the Grevelingen.

Methods

Vegetation

For monitoring plant species, permanent quadrats (PQs) were used (see Fig. 1d). Species numbers and cover percentages were registered following Braun-Blanquet and processed in TURBOVEG (Version 2.118, Hennekens and Schaminée 2001). Ellenberg indicator values (Ellenberg values) were calculated (Ellenberg et al. 1992) to obtain qualitative insight into abiotic conditions' development over the monitoring period. Subsequently, we grouped the values according to the gradients present. In JUICE (version 7.0.84, Tichý 2002), we clustered our data with TWINSPAN (Hill 1979) and constructed a synoptic table enabling discernment of succession series. Online Resource 2 gives additional details on monitoring and data processing.

Hydrology

We assessed the main hydrological conditions at the Slikken van Flakkee for an impression of water flow as it directly influences water composition. These main conditions are the extent of floodings or salt spray by Lake Grevelingen water and the permeability of the sediment; the former can raise or even reset Cl− in the porewater, and the latter influences flow rate and, consequently, the rate of de- or resalinisation. Although we know an example in which the flooding extent was mapped (Slager et al. 1986), in most cases, they are unknown because the surface-water level at the Slikken van Flakkee is not recorded anymore, and the continuous recordings of the water level on Lake Grevelingen underestimate the surface-water level at the Slikken van Flakkee (Hissink 2018). As an alternative, we used Cl− concentrations and surface elevations at the PQs as a proxy to estimate the extent of floodings or salt spray. Large-scale permeability was deduced from the map ‘Top impermeable layer’ (Slager et al. 1986), and permeability values were calculated from measured grain size distributions in samples from all PQs.

Geochemistry

We took 36 soil samples in 2016 at 9 representative PQs and 124 samples in 2017 at all still existing 63 PQs. From the 2016 samples, we subsampled the porewater while conserving field redox conditions, and the 2017 samples we extracted with H2O. Porewater subsamples, extracts and soil subsamples were chemically analysed, and the data were processed for the composition of porewater (PW), reconstructed porewater (PWre) and soil. The 2017 soil samples were also analysed on grain size distributions. As a stepwise approach for quantifying the geochemical processes, firstly, we assessed the occurring geochemical processes (see “Geochemical processes”), secondly, we determined the water-types of the porewater compositions (see “Porewater composition”), and lastly, we inversely modelled the water compositions with PHREEQC (Parkhurst and Appelo 2013) to quantify the geochemical processes that led to these water compositions (see “Inverse modelling”). Online Resource 3 gives a detailed description of the fieldwork and its two campaigns, the laboratory work (including the applied analytical methods as listed in Table S2), the data processing (including the calculation of difference concentrations, symbolised with \(\Delta\)), and the inverse modelling. Both calculation and inverse modelling are based on the end members presented in Table 1.

Table 1 Water compositions for inverse modelling
Full size table

Interrelations

In addition to the grouped Ellenberg values to express interrelations, we calculated average concentrations of porewater and soil contents per present vegetation type. We also carried out a Canonical Correspondence Analysis (CCA) (Legendre and Legendre 2012) with the vegetation data of 2016 (or, if not available, of the last year for which data were available) that were constrained with the geochemical data of 2017. The vegetation data were transformed by calculating each plant species’ cover as the fraction of the sum of all plant species’ cover (leading to a number between 0 and 1), and environmental variables were ln-transformed. For the ordination, the geochemical and texture variables were forward selected by ordistep and tested by anova, both functions within the vegan package in R (Oksanen et al. 2019). Ellenberg values, management, elevation and soil composition (as a proxy for permeability of the soil profile) were fitted to the ordination by the function envfit. Allocation of the variables (selecting or fitting) was based on the type of variables [direct (geochemical and texture) or indirect (all other), respectively] following Zelený and Schaffers (2012).

Results

Vegetation

Synoptic table

The result of our clustering of the c. 30,000 observations (2047 relevés; 229 plant species) is 9 clusters (vegetation communities) and 18 sub-clusters (vegetation types). Table 2 presents these clusters in which 43 diagnostic plant species are used (19% of all plant species; 62% of all observations). The clusters are (1) Halophytic pioneer community, (2) Transition community with Puccinellia maritima and Epilobium hirsutum, (3) Transition community with Tripolium pannonicum and Agrostis stolonifera, (4) Grassland communities with Juncus gerardii and Agrostis stolonifera, (5) Grassland vegetation dominated by sown grasses, (6) Grassland communities related to cluster 7, (7) Transition community with Calamagrostis epigejos and Hippophae rhamnoides, (8) Communities of tall-growing forbs and grasses on the former salt marshes, (9) Shrubland with Hippophae rhamnoides and Salix species. Online Resource 4 describes these clusters, their sub-clusters and their syntaxonomical position.

Table 2 Synoptic table of vegetation types in Slikken van Flakkee 1972–2016
Full size table

Vegetation development

Vegetation types and succession series

Succession is still continuing in some PQs. Figure 2a expresses the average change rate of the vegetation development along the succession series. Since 2000, this rate has stabilised at 0.05–0.10, i.e., in 5–10% of all PQs, vegetation type changes each year. Vegetation development is along 9 identified succession series as shown in Online Resource 4: (A) Regular inundated saline shore as pioneer zone, (B) Low-lying saline areas, (C) Shore with deposited organic material, (D) Largely desalinated grasslands, (E) Diverse grasslands and vegetations with Salix repens, (F) Thickets with Hippophae rhamnoides, (G) Thickets and bushes with Salix cinerea, (H) Thickets and bushes with Salix caprea, (I) Former salt marshes. The photos in Online Resource 5 show the landscape and all present vegetation types.

Fig. 2
figure 2

a Moving averages of change rate vegetation type versus time. Change rate = ‘sum changed vegetation types’/‘number PQs’ in monitoring year. The 21 years’ moving averages were allocated to the central year. b Average number of plant species per PQ per area versus time for Slikken van Flakkee and sub-areas. Vertical lines express standard deviation multiplied by 0.2 to minimise overlap

Full size image
Plant species richness

Species richness is plotted versus time for all PQs (Fig. 2b), showing a rapid increase in plant species numbers since eliminating the tide in 1971. After increasing the first three years to an average of c. 13 species per PQ, numbers remained relatively constant; however, we observed slight differences between the sub-areas. In ‘North’, plant species richness increased more steeply in the first 3 years to an average of c. 15 species per PQ, subsequently continued increasing until the average reached a maximum of 18 species per PQ in 1981, and, since then, has decreased gradually to an average of c. 13 species per PQ. The standard deviation shows a decreasing trend. In ‘South’, the plant species numbers have gradually increased to c. 22 species on average per PQ from 1971 onwards, and standard deviation remained relatively constant. These developments reflect the spontaneous development in ‘North’ versus the grazing and mowing in ‘South’, the latter creating opportunities for higher species diversity through more favourable light conditions and low nutrient conditions.

Red list plant species

Of the 229 occurring plant species in the Slikken van Flakkee, 36 are rare, 29 are threatened or vulnerable, and 24 are rare and threatened or vulnerable, according to the Red List. Table 2 presents only some; the complete list is in Online Resource 6. The total number of annual observations of rare and threatened or vulnerable plant species fluctuated between 3 and 18% of all species and was, on average, c. 9%, and per PQ in ‘South’ 1.35 times higher than in ‘North’.

Ellenberg values

Online Resource 7 shows how Ellenberg values vary per sub-area and elevation class for all PQs over time. The development of each Ellenberg value is briefly described in this Resource. Of all Ellenberg values, the salinity value is most clearly dependent on the surface elevation and secondary on soil composition as a proxy for permeability of the soil profile (see Online Resource 8).

Hydrology

The flooding border was mapped for the severe storm on 26 November 1983: between 0.5 and 1.0 m m.s.l. Figure 3a shows that the high Cl− concentrations in the porewater occur till c. 0.4 m m.s.l., raised only below that elevation by floodings or salt sprays from the salty Lake Grevelingen. If the high Cl− concentrations are caused by evaporation, the extent of the floodings or salt sprays is even smaller.

Fig. 3
figure 3

a Cl− concentration in porewater (PW) from 2016 samples at 4 depths and H2O extract (reconstructed porewater; PWre) from 2017 samples at 2 depths versus surface elevation. Cl− concentration range is zoomed to 0–20,000 mg/l. b SO42− concentration versus Cl− concentration in PW and PWre. Blue line represents the conservative mixing of original porewater and local rainwater. Both axes are logarithmic. Cl− concentration range is zoomed to 1–20,000 mg/l, and SO42− concentration range is zoomed to 1–10,000 mg/l. c pH in porewater and H2O extract versus CaCO3 content in soil samples. d CaCO3 content versus organic matter content in soil samples.

Full size image

For the A-horizon, the mean grain size ranges from fine (125–250 µm) to medium sand (250–500 µm), and for the B-horizon, from very fine (62–125 µm) to fine sand (125–250 µm). The organic matter content in the A-horizon (on average c. 14% DS) shifted the grain size distributions towards the coarser class; in the B-horizon, organic matter is almost absent (on average c. 1% DS). The sorting of the grains in the B-horizon ranges from very poorly sorted to well-sorted, and its average is moderately sorted. The sludge (< 16 µm) fraction in the B-horizon ranges from 0 to 5.5%, and its average is c. 1%.

Calculated permeability in the B-horizon decreases with elevation: near-shore, it is on average c. 13 m/day, and near the dike, on average c. 7 m/day, classified as permeable and semi-permeable, respectively (Bear 1988). We interpret this permeability-elevation relation as imposed by the tide-related current of floodwater before the closure of the Grevelingen in 1971. This current decreased in strength from shore to dike (because of the increased resistance of the soil surface by the vegetation present) and deposited the grains accordingly. These calculated permeabilities are higher than those reported by Stafleu and Gunnink (2016) for the lithoclasses in the Walcheren Member that outcrop in the Slikken van Flakkee (TNO-GSN 2021): 0.02–5 m/day. The lithoclass peat of the Hollandveen Member, present at some m bgs, has the lowest permeability: 0.001 m/day, matching the impermeable layer shown in Fig. 1c. The permeability of the Wormer Member underneath ranges 0.02–5 m/day.

Consequently, we regard the hydrological system as a locally thin unconfined aquifer, being infiltrated with rainwater (that replaces the original porewater at a permeability-dependent velocity), being flooded near-shore incidentally (with decreasing frequency) or salt sprayed or evaporated.

When the water level management in Lake Grevelingen becomes partial tidal re-established and climate-adapted (based on the present most likely scenario), the average water level will elevate stepwise (≤ 2029: – 0.20 m m.s.l.; 2030–2046: – 0.14 m m.s.l.; ≥ 2066: – 0.04 m m.s.l.) (Boeters 2022), leading in 2066 to the reduction of c. 2.5 km2 of the area of the Slikken van Flakkee (being c. 15% of the present total) and consequently loss of Natura 2000 habitat types and species.

Geochemistry

Geochemical processes

The identified geochemical processes during desalinisation are redox processes, dissolution of carbonates, cation exchange and accumulation and decomposition of organic matter. Our results show that desalinisation progress in the A-horizon is, on average, 88%, and in the B-horizon, on average, 78%. The oxidation of sulphides in the soil led to increased SO42− concentrations; reduction of sulphate to decreased SO42− concentrations (Fig. 3b). The dissolution of carbonates prohibited acidification (Fig. 3c), even in the A-horizon with low CaCO3 contents (Fig. 3d). Cation exchange was revealed by \({\mathrm{Na}}_{\Delta }^{+}\), \({\mathrm{K}}_{\Delta }^{+}\) and \({\mathrm{Mg}}_{\Delta }^{2+}\). The accumulated organic matter contains carbon and nitrogen in a depth-dependent ratio. We described these processes more extensively in Online Resource 9.

Porewater composition

Figure 4 visualises all water compositions grouped per vegetation type in a Piper diagram. Based on this figure, we interpret that in the Slikken van Flakkee, the water composition develops from the original saltwater-type NaCl (right in the diamond) to the freshwater-type CaHCO3 (left in the diamond). The Figure also shows the diverging from the water-type NaCl and converging to the water-type CaHCO3: Vegetation types 1, 2, 6 and 11 show the slightest variation in water composition, the vegetation types with the intermediate compositions (7, 12, 16, 17 and 18) have a larger variation, while vegetation type 15 has a specific water composition. Vegetation type 12 shows the most extensive variation, indicating that management is more decisive for this vegetation type than the water composition. Online Resource 10 shows the Piper diagram with the individual points from which the convex hulls for each vegetation type group were constructed.

Fig. 4
figure 4

Piper diagram with porewater compositions grouped per vegetation type. The main anions are presented in the right triangle, the main cations in the left triangle, while the triangles’ positions are combined in the diamond. The concentrations are in meq/l%, and the encircled numbers refer to the vegetation types.

Full size image

Inverse modelling

Figure 5 illustrates TDS concentrations in PW and PWre and those by conservative mixing of the end members. \(\text{TDS}_{\Delta }\) is the net result of the occurring geochemical processes during desalinisation; the assumed processes and their phases used in the inverse modelling are specified in Table 3.

Fig. 5
figure 5

Total dissolved solids (TDS) concentration in PW and PWre versus fraction rainwater. Blue line represents the conservative mixing of original porewater and local rainwater. The dotted line was drawn manually through the TDS concentrations in porewater. Distance between both lines represents \({\text{TDS}}_\Delta\). TDS range is zoomed to 0–40,000 mg/l

Full size image
Table 3 Average phase transfers for occurring geochemical processes since reclamation per method and horizon
Full size table

In the modelling, all \({I}_{\Delta }\) were allocated to phase transfers. Online Resource 11 presents a selection of these phase transfers versus Cl−. The figures show that phase transfers decrease during desalinisation; transfer maxima are c. – 3 g/l for Fe(OH)3 and c. 6 g/l for CaCO3. Table 3 also presents all phase transfers’ averages per method and horizon. The table demonstrates that the dissolution of CaCO3 is enhanced by the entering of gasses, dissolution of sulphides, precipitation of iron hydroxides and cation exchange; averages of the sum of net phase transfers vary between 0.80 and 2.49 g/l. This range means that the extent of phase transfers in the soil depends on the sampled horizon and the method; especially the phase transfers of PWre in the A-horizon are enhanced, probably because of oxidation since sampling. Because the sum of net transfers is almost equal to \(\text{TDS}_{\Delta }\), we conclude that the inverse modelling simulated the geochemical processes rather well.

Interrelations

Chemical characterisation of vegetation types

The interrelations found in the Piper diagram (Fig. 4) motivated us to characterise the present vegetation types chemically, as shown in Fig. 6a by the diminishing TDS concentration trend (clear and consistent for both campaigns and horizons), in Fig. 6b as the opposite trend with increasing CaCO3 phase transfer, and in Fig. 6c as varying organic matter content (revealing that the highest contents occur in vegetation type 11 and 12 and thus coincide with management). Online Resource 12 contains the complete chemical characterisation of the vegetation types.

Fig. 6
figure 6

Chemical characterisation of vegetation types as averages per vegetation type for a TDS concentration, b CaCO3 phase transfer, and c organic matter content, and with d belonging map of vegetation types. Symbol size in the map represents the number of PQs along a transect with the same vegetation type

Full size image

CCA

In addition to this chemical characterisation, Fig. 7 shows the CCA model’s plots with the selected geochemical and texture variables (that significantly relate to vegetation) and the variables fitted to this model. The positions of the sample location (clustered by the vegetation types) correspond rather well with the two main and opposite gradients in the field: elevation and saline conditions in the root zone (represented by TDS-A and Cl-A), and with management that is perpendicular to these gradients.

Fig. 7
figure 7

a Canonical Correspondence Analysis plot (CCA[vegan]-R) showing relations between cover of plant species (2016) constrained by selected geochemical and texture variables from samples (taken in 2017 at 2.5 and 40 cm bgs) from 63 locations in the Slikken van Flakkee. To model fitted variables are also shown. Eigenvalue axis 1 = 0.7711, axis 2 = 0.7189; 43.7% variance accounted for by these two axes. Selected model variables and variables fitted to the model are described and quantified in Table 4. Sample locations belonging to a specific vegetation type are indicated with black polygons with a number referring to the vegetation type. The other vegetation types (3–5, 8–10, 13 and 14) are absent because they did not occur in 2016. b CCA plot with Ellenberg values fitted to the same model

Full size image

The plot also reveals that the grain size fraction ‘very fine sand’ (V_F_Sand) correlates positively with elevation, while the grain size fraction ‘fine sand’ (F_Sand) correlates negatively with elevation. Cl− and SO42− correspond with halophytic vegetation types 1 and 2. The CCA model explains c. 44% of the total variance of the plant species, and 14 variables are significant (Table 4). TDS in the 0–5 cm depth samples explains the largest proportion of variance. Variables measured in the shallow samples explain c. 26% of the variance and those measured at 35–45 cm depth c. 17%, indicating that environmental conditions at both depths seem relevant for vegetation. Ellenberg light and salinity values and elevation strongly correlate to the CCA model.

Table 4 Forward selected CCA model variables and variables fitted to the model
Full size table

The vegetation types can be characterised as follows:

  • 1 and 2: Near-shore locations with the highest TDS-A and \(\text{TDS}_{\Delta}\text{-B}\) concentrations, highest Ellenberg light, salinity and moisture values combined with the lowest elevation and Ellenberg nutrient values.

  • 6, 7, 15–18: Locations along decreasing gradients of TDS-A, fine sand-A, Ellenberg light, salinity and moisture values, and increasing gradients of CaCO3-A, silt-A, \(\text{TDS}_{\Delta}\text{-B}\), clay-A, elevation and Ellenberg nutrient values.

  • 11, 12: Locations all present in ‘South’ and thus related to management (grazing and mowing).

Discussion and conclusions

Discussion

After Van Haperen (2009) evaluated the long-term vegetation monitoring in the former tidal wetland Slikken van Flakkee, we took up the unique possibility of actualising vegetation monitoring and complementing it with inverse modelling of the occurring geochemical processes. This modelling also enabled us to verify hypothesised processes during desalinating (Appelo and Postma 2005). By quantifying ecohydrogeochemical processes, we assess how processes evolve and are interrelated, and what impact management has.

We distinguished 18 vegetation types that formed specific succession series from 1972 to 2016 in the Slikken van Flakkee. A total of 229 individual plant species were found, and their numbers have changed over time, dependent on the sub-area. In ‘North’ (spontaneous development), the average number of plant species stabilised at c. 14 plant species per PQ, while in ‘South’ (with management consisting of grazing and additional mowing), the number of species is still increasing. In ‘South’ in 2016, plant species reached an average number of 22 per PQ. Similarly, Esselink et al. (2002) also found an increased average number of species per PQ (from 8 to 10) as the effect of management in a tidal salt marsh in the Dollard, the Netherlands, although in that case, reduction of both grazing intensity and drainage were confounding factors. Rare and threatened plant species constituted c. 9% of all observed plant species in our area; in the actively managed part, the number was 1.35 times higher than in the spontaneously developed part, illustrating how effective active nature management contributes to plant diversity and protection of vulnerable species.

Not unexpectedly, Ellenberg values reflect the effect of developing hydrological and geochemical conditions and management measures since the tide-stop in 1971 (see Online Resources 7 and 8). Salinity is the most time-dependent and thus clearly changing gradient in contrast to elevation and management, both constant in time. The other Ellenberg values also changed over time but within smaller ranges.

The Piper diagram clearly expresses the water-type shift from NaCl towards CaHCO3; the largest variation in water composition is observed for vegetation type 12. Cation exchange is limited, unlike other Dutch former tidal wetlands that contain clay (Hissink 1938). We also observed organic matter accumulation, which is in line with the humus accumulation as found by Slager and Visser (1990), who also considered peat-forming in ‘North’ a possibility. However, we found higher organic matter content in the managed sub-area (‘South’) than in the spontaneously developed sub-area (‘North’), suggesting that the vegetation-driven accumulation is caused directly by management. More likely, it reflects the different distribution of the organic matter in the soil profiles: in ‘South’ concentrated at shallow depth; in ‘North’ more dispersed because of the deeper roots.

The inverse modelling clarified that desalinisation was followed by oxygen and carbon dioxide entrance, oxidation of sulphides, precipitation of iron hydroxides, dissolution of carbonates, cation exchange and organic matter oxidation. Because of the excess carbonates present, the conditions remain alkaline to neutral, i.e., no acidic conditions occur, even in the A-horizon, where the accumulated organic matter influences the abiotic conditions greatly. These identified processes were quantified by the inverse modelling, which revealed that during desalinisation, the dissolution of CaCO3 is the most dominant process and is enhanced by the other processes. On average, c. 1–3 g/l enters the solution, ultimately leading to a clear CaHCO3 water-type.

Diving deeper into how vegetation is related to environmental conditions and management, we observed that vegetation types clearly relate to TDS concentration, CaCO3 phase transfer and organic matter content. We may also conclude from the CCA that elevation and total dissolved solids (TDS) were identified as two major-related gradients, management as a distinctive factor and several geochemical variables as derivative gradients for the occurrence of the vegetation types in the Slikken van Flakkee. The CCA showed their position along these gradients and management groups. The early successional stages (vegetation types 1 and 2 = ‘1–2’) are near-shore locations with low elevation and the highest TDS concentrations. All other successional stages with spontaneous development (grassland, i.e., ‘6–7’; forbs, i.e., ‘15’; shrubland, i.e., ‘16–18’) occur along the gradients. These gradients are increasing (elevation and in order of variance: CaCO3-A, silt-A, \(\text{TDS}_{\Delta}\text{-B}\) and clay-A) and decreasing (TDS-A and fine sand-A). Vegetation types ‘11–12’ only occur in ‘South’ and thus are related to management. Clearly, vegetation on elevated sites experiences desalinisation. Our hydrogeochemical data and inverse modelling confirm several processes already qualitatively described by Joenje and Verhoeven (1993): At places where freshwater lenses develop, desalinisation and dissolution of CaCO3 occur.

Furthermore, our long-term time series shows that vegetation succession still develops towards species-rich meadow vegetation when managed by grazing and mowing: in ‘South’, only vegetation types 11 and 12 developed (see Online Resource 4), limited expressing hydrological-geochemical gradients present. This management of periodic biomass removal is crucial in most management schemes applied in nature reserves to enhance biodiversity, e.g., as concluded already by Van de Rijt et al. (1996) for the adjacent former estuary Haringvliet and by Joenje and Verhoeven (1993) for the older Dutch former tidal wetland Piamer Kooiwaard. Nevertheless, based on similar wetlands worldwide, e.g., Bantilan-Smith et al. (2009), Garde et al. (2004), Karberg et al. (2018), Wu et al. (2020), we expected that hydrological-geochemical gradients would be the most dominant factors for developing vegetation types in the Slikken van Flakkee. Initially, this was indeed the case in our area, but gradually management gained more influence on vegetation development and appeared to be a dominant factor for the vegetation succession’s direction in addition to the gradients present.

In general, the soil of the Slikken van Flakkee can be characterised as calcareous clay-less fine-sandy sediment which, since the tide-stop, gradually desalinates, dissolves lime, and gains organic matter from the developing vegetation depending on management. Although these are universal ecohydrogeochemical processes, they have led to specific conditions that made this former tidal wetland unique worldwide.

A climate-adapted plan, in which partial re-establishment of the tide is meant to improve the water quality in Lake Grevelingen (Wijsman et al. 2022), will when implemented, lead to higher lake water levels in time. We expect this will shift the habitats towards higher elevation at the same pace as the water level rises and shrink the land area, which is the main adverse effect. With progressing sea-level rise, this seems inevitable. Therefore, we suggest compensation measures to be taken in advance since the area is protected by law. These measures can now be optimised by making use of our results, e.g., by creating new land or re-allocating land currently used for agriculture to nature and realising favourable conditions for vegetation development by measures, e.g., avoiding acidification, creating salinity and elevation gradients, using freshwater seepage, and applying low-intensity mowing and grazing. Applying our method in comparable restoration projects elsewhere, e.q. the Ni-les’tun restoration project on the southern Oregon coast (Janousek et al. 2022), may be beneficial for quantifying ecohydrogeological processes and help identify potential succession lines.

Conclusions

We described vegetation development and related succession to quantified ecohydrogeochemical processes in the Slikken van Flakkee. The availability of the unique vegetation time series appeared to be very informative. In c. 50 years, the area developed from tidal mudflats, salt marshes and interconnected creeks towards grassland, tall forbs and thickets under the influence of changing physico-chemical conditions where desalination is the dominant and ongoing process causing a shifting salinity gradient. Desalination is a function of the hydrological conditions in the area dictated by the elevation gradient. During desalinisation, carbonates dissolve, leading to a water-type shift from NaCl towards CaHCO3. On average, \(\text{TDS}_{\Delta }\) equals c. 1–3 g/l, as quantified by inverse modelling. The applied management measures (grazing and mowing) appeared to be very important for successional development towards grassland and species-rich forb vegetation harbouring many rare species. Identifying the rate and extent of long-term ecohydrogeochemical changes appeared to be relevant not only for our case area but also for comparable former tidal wetlands after tide-stop, thus expanding our knowledge of the functioning of these wetlands and the impact of management.

Recommendations

Because we expect that the vegetation in the Slikken van Flakkee will keep developing, we recommend continuing monitoring to follow the development closely. To support further the policy-making for the possible partial re-establishment of the tide in Lake Grevelingen by predicting the vegetation development and the related hydrological-geochemical conditions, we also recommend modelling these conditions spatio-temporal at different water level scenarios with MODFLOW6 (2022) and PHREEQC, similar to Karlsen et al. (2012).