Ecohydrogeochemistry of the Slikken van Flakkee, a former tidal wetland in the Netherlands

The Slikken van Flakkee, a former tidal wetland in the Grevelingen, the Netherlands, has been changing under freshening conditions since 1971; in that year, the salty Grevelingen estuary was closed from the North Sea as part of the Dutch Delta Works. Our research aimed to analyse vegetation development over 1972–2016 and relates succession to quantified ecohydrogeochemical processes in the Slikken van Flakkee. We analysed succession and physico-chemical conditions at 63 monitoring locations. Based on extensive fieldwork in 2016 and 2017, we quantified the identified geochemical processes via inverse modelling and applied Canonical Correspondence Analysis (CCA) to search for interrelations with vegetation. Nine per cent of the 229 observed plant species are rare and threatened or vulnerable, according to the Dutch Red List. Succession series showed a development from saline pioneer communities towards grassland and shrub vegetation dependent on management measures (grazing and mowing versus spontaneous development), desalinisation and geochemical processes, of which dissolution of carbonates appeared to be the dominant one. During desalinisation, the water-type shifted from NaCl towards CaHCO3. The CCA confirmed the impact of geochemical and texture variables, elevation and management on vegetation. Grazing and mowing led to species-rich grasslands and forb vegetation harbouring twice as many species than the spontaneously developed area after 50 years, demonstrating the importance of active nature management. We argue that hydrological-geochemical modelling is necessary to support policy-making before deciding on the partial re-establishment of the tide in Lake Grevelingen.


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 saltfresh 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 CaCO 3 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 CO 2 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 CaHCO 3 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 reestablishment 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 km 2 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).
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.

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 1 3 the 2016 samples, we subsampled the porewater while conserving field redox conditions, and the 2017 samples we extracted with H 2 O. Porewater subsamples, extracts and soil subsamples were chemically analysed, and the data were processed for the composition of porewater (PW), reconstructed porewater (PW re ) 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 Δ ), and the inverse modelling. Both calculation and inverse modelling are based on the end members presented in Table 1.

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).

Vegetation development
Vegetation types and succession series Succession is still continuing in some PQs. Figure  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 subareas. 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.   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 permeabilitydependent 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 (

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 SO 4 2− concentrations; reduction of sulphate to decreased SO 4 2− concentrations (Fig. 3b). The dissolution of carbonates prohibited acidification (Fig. 3c), even in the A-horizon with low CaCO 3 contents (Fig. 3d). Cation exchange was revealed by Na + Δ , K + Δ and Mg 2+ Δ . The accumulated organic matter contains carbon and nitrogen in a depth-dependent ratio. We described these processes more extensively in Online Resource 9. 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 saltwatertype NaCl (right in the diamond) to the freshwater-type CaHCO 3 (left in the diamond). The Figure also shows the diverging from the water-type NaCl and converging to the water-type CaHCO 3 : 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. Figure 5 illustrates TDS concentrations in PW and PW re and those by conservative mixing of the end members. TDS Δ 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.

Inverse modelling
In the modelling, all I Δ 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 CaCO 3 . Table 3 also presents all phase transfers' averages per method and horizon. The table demonstrates that the dissolution of CaCO 3 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 PW re in the A-horizon are enhanced, probably because of oxidation since sampling. Because the sum of net transfers is almost equal to TDS Δ , we conclude that the inverse modelling simulated the geochemical processes rather well.

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 CaCO 3 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.

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  (229)      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. 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 SO 4 2− 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)

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 timedependent 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 CaHCO 3 ; 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 subarea ('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 CaCO 3 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 CaHCO 3 water-type.
Diving deeper into how vegetation is related to environmental conditions and management, we observed that vegetation types clearly relate to TDS concentration, CaCO 3 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: CaCO 3 -A, silt-A, TDS Δ -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 CaCO 3 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.  High  45  10  17  Middle  49  58  22  Low  7 3  5 9  6 5   Table 3 Average phase transfers for occurring geochemical processes since reclamation per method and horizon 1) Averages per method and per horizon (Porewater-A: n = 9; Porewater-B: n = 9; H 2 O extract-A: n = 54; H 2 O extract-B: n = 57) 2) The exchange complex X remains part of the solid phase and, therefore, does not contribute to the sum of the phase transfers 3) The sum of all cation exchange expressed in meq/l is zero, but because of the different mass and charges of the involved cations, their sum in g/l is not zero 4) The sum of the modelled absolute phase transfers, i.e. the amount that has entered the porewater (by dissolution and desorption) plus the amount that left the porewater (by precipitation and adsorption) 5) The net sum of the modelled phase transfers, i.e. the net modelled amount that has entered the porewater and has increased the modelled amount of dissolved solids in the porewater 6) The increase of the measured total dissolved solids (TDS) concentration relative to the TDS concentration by conservative mixing of original porewater and local rainwater 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  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 elevation gradients, using freshwater seepage, and applying low-intensity mowing and grazing. Applying our method in comparable restoration projects elsewhere, e.q. the Niles'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 Significance codes 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1 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 CaHCO 3 . On average, TDS Δ 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).