The potential natural vegetation of large river floodplains – From dynamic to static equilibrium

The potential natural vegetation (PNV) is a useful benchmark for the restoration of large river floodplains be- cause very few natural reference reaches exist. Expert-based approaches and different types of ecological models (static and dynamic) are commonly used for its estimation despite the conceptual differences they imply. For natural floodplains a static concept of PNV is not reasonable, as natural disturbances cause a constant resetting of succession. However, various forms of river regulation have disrupted the natural dynamics of most large European rivers for centuries. Therefore, we asked whether the consideration of succession dynamics and time dependent habitat turnover are still relevant factors for the reconstruction of the PNV. To answer this we compared the results of a simulation of the vegetation succession (1872–2016) of a segment of the upper Rhine river after regulation (damming, straightening and bank protection) to different statistic and expert-based modelling approaches for PNV reconstruction. The validation of the different PNV estimation methods against a set of independent reference plots and the direct comparison of their results revealed very similar performances. We therefore conclude that due to a lack of large disturbances, the vegetation of regulated large rivers has reached a near-equilibrium state with the altered hydrologic regime and that a static perception of its PNV may be justified. Consequently, statistical models seem to be the best option for its reconstruction since they need relatively few resources (data, time, expert knowledge) and are reproducible.


Introduction
River floodplains are amongst the most species-rich and productive ecosystems (Naiman and Decamps, 1997). At the same time, these ecosystems are one of the most threatened and modified worldwide (Tockner and Stanford, 2002), highlighting the need for conservation and restoration efforts (Buijse et al., 2002;Myers et al., 2000). Such efforts, however, are challenged by a lack of natural reference sites for orientation (Whited et al., 2007). Indeed, 90% of river floodplains in Europe and North America are used for agriculture or forestry and no longer harbor natural vegetation communities (Tockner and Stanford, 2002). Where no natural references exist, a potential natural vegetation (PNV) is often reconstructed and used as benchmark (Carranza et al., 2003;Hickler et al., 2012;Justice et al., 2017;Klimas et al., 2009;Schleupner and Schneider, 2013;Shi et al., 2016).
But the concept of PNV and the methods for its reconstruction are highly controversial (Chiarucci et al., 2010;Loidi and Fernández-González, 2012;Mucina, 2010;Somodi et al., 2012). PNV was first defined by Tuxen (1956) as the vegetation that would develop under present site conditions if human influences were excluded completely, and succession would reach its climax stage at once. It has often been thought of as an historic, pre-human reference condition (Hall and McGlone, 2006;Willis and Birks, 2006). However this idea has provoked disagreement (Carrión and Fernández, 2009;Mitchell, 2005) because it ignores that environmental conditions have changed since pre-human times (Dotterweich, 2008;Nilsson et al., 2005). Therefore, it has been argued that an estimation of the natural vegetation based on the assessment of present-day natural vegetation remnants is more reliable (Kowarik, 1987). But in areas with historically high levels of land use transformation this assessment is also prone to uncertainties (Zerbe, 1998). Another much discussed issue is the consideration of ecosystem dynamics and stochasticity (Chiarucci et al., 2010;Härdtle, 1995), which are especially relevant in naturally disturbed areas (Jackson, 2013;Leuschner, 1997;Zerbe, 1998). The traditional PNV estimation method is expert-based and follows a floristic-sociological approach (Westhoff and Van Der Maarel, 1978). It relies on the fieldworker's assessment and understanding of ecology to extrapolate present-day natural vegetation remnants to similar environments (Kowarik, 1987;Moravec, 1998). This approach, however, lacks transparency and reproducibility. Furthermore, its implied static perception of PNV makes predictions in dynamic systems questionable (Mucina, 2010). More comprehensive and also widely used for PNV estimations are ecological models based on the relationship between vegetation and environmental variables (Somodi et al., 2012;Zerbe, 1998). Two types can be differentiated that imply fundamental conceptual differences: static models and dynamic models (Hannon and Ruth, 1998). Static models describe a phenomenon at a given point in time and assume equilibrium between the vegetation and its environment. Dynamic models (e.g. process-based models or mechanistic models) are based on ecological processes and differ from static models by explicitly incorporating timedependent changes in the system state. Therefor they are able to capture the transient response of vegetation to a changing environment (Hannon and Ruth, 2014). Little attention, however, has been given to issues surrounding PNV in fluvial contexts. The first modelling approaches of PNV in floodplains were static and tried to explain the vegetation patterns along vertical (height above channel) and lateral (distance away from channel) gradients (Bowman and Mcdonough, 1991;Ellenberg, 1996;Glavac et al., 1992;Hosner and Minckler, 1963;Hughes, 1988;Nixon et al., 1977;Roberts and Ludwig, 2016;Robertson et al., 1978;Ward and Stanford, 1995). Later more advanced static models emerged that relate the vegetation distribution to a set of hydrologic variables based on expert rules (Aggenbach and Pelsma, 2003;Baptist et al., 2004;Fuchs et al., 2012;Jungwirth et al., 2002;Lenders et al., 2001;Pieterse et al., 1998;Runhaar, 2003) or statistical analyses (Auble et al., 1994;Franz and Bazzaz, 1977;Menuz, 2011). Dynamic floodplain vegetation models combine a simulation of the hydro-dynamics with the modelling of ecological processes (e.g. growth and mortality, recruitment, succession/retrogression, competition) and are used to describe the vegetation development over time (Benjankar et al., 2011;Camporeale and Ridolfi, 2006;García-Arias and Francés, 2016;Pearlstine et al., 1985). More recently, the dynamic feedbacks between vegetation and hydrogeomorphological processes have also been incorporated (Camporeale et al., 2013;van Oorschot et al., 2016). We argue that for floodplains of unregulated rivers the original static concept of PNV based on a climax stage of vegetation is not reasonable because successional sequences are repeatedly rejuvenated and reset by hydro-geomorphological disturbances (Pringle et al., 1988). It has been theorized that at the appropriate scale the proportion of successional phases would remain constant when processes of regression are compensated by progression, a dynamic equilibrium referred to as shifting steady-state mosaic (Bormann and Likens, 1979;Geerling et al., 2006;Stanford et al., 2005). While methodological advances reflect the growing recognition in the importance of allowing dynamic change in PNV estimation, various forms of river regulation have disrupted the natural dynamics in most large rivers (Buijse et al., 2002;Dynesius and Nilsson, 1994). River damming and bank stabilization are the main reason for the impediment of dynamic geomorphological processes (e.g. avulsion, meandering, braiding) and decrease of hydrodynamic variability and disturbance in the riparian zone (Church, 1995;Magilligan and Nislow, 2005;Nilsson and Berggren, 2000;Petts and Gurnell, 2005).
In this context, we investigated whether the consideration of succession dynamics and habitat turnover are still relevant factors for the model-based reconstruction of the PNV of regulated large river floodplains. Our hypothesis is that they can be neglected because riparian vegetation of regulated large rivers has reached a stable equilibrium due to the loss of natural disturbance dynamics. To test this idea we compared the results of a simulation of the succession dynamics (1872-2016) of the floodplain vegetation of a segment of the heavily regulated upper Rhine River to different static approaches for the estimation of its PNV a) a statistical model based on hydrologic predictors and the geomorphological age of site b) a statistical model only based on hydrologic predictors and c) a gradient approach only based on the distance to the mean water level.

Study area
The Rhine River is one of the largest rivers in Central Europe, with a length of approximately 1230 km and a catchment area of approximately 185300 km 2 . Our study area lies in the Upper Rhine region where the nival discharge regime is strongly influenced by snow-melt in the Alps (Belz and Frauenfelber-Kääb, 2007). Until the beginning of the 19th century the Upper Rhine could still be considered in natural condition (Gallusser and Schenker, 1992) and was classified as an highly dynamic, island-dominated, anabranching river system (Gurnell and Petts, 2002;Herget et al., 2005). During the course of the 19th century, however, the Upper Rhine River was transformed into a singlethread channel by cutting off meander bends and building groins and bank revetments (Bernhardt, 2000). In the 20th century river regulation intensified through the construction of 10 hydropower plants in the main channel or in artificial side channels (Dister et al., 1990). During this time industry and settlements also expanded in the study area (Habersack and Piégay, 2007). The study area is the "Raststatter Rheinaue", a nature reserve on the eastern, German side of the floodplain that includes a 9-km segment of the Upper Rhine River downstream from the Iffezheim dam to the confluence of the river Murg (Rhine km 335.8-345, 114-110 m a.s.l.). It is only limited by flood dykes towards the east and is still regularly flooded. The study area covers approximately 645 ha (including water bodies).

Historic maps
Because the Upper Rhine has been the border between France and Germany, detailed maps were produced for the planning of the river straightening in the beginning of the 19th century. These indicate the location of water bodies, islands and gravel/sand bars within the aquatic area, as well as land uses in the floodplain (grasslands, forests, croplands and settlements). Our work is based primarily on four historical maps from that time (1816,1838,1852,1872) that were georeferenced and classified in natural (natural water body, gravel/sand bar, grassland and forest) or anthropic (artificial water body, cropland, settlement and industry) habitat categories (for details see Table C1 and Diaz-Redondo et al., 2017).

Discharge data and hydrodynamic model
The analysis of the flow regime of the study area for the whole simulation period is based on the Maxau gauging station (Rhine km 362.3) which has the longest continuous record of daily discharge (1921 -today) and also records of the annual low, mean and high discharges for the period 1872-1921(Diaz-Redondo et al., 2017. A two-dimensional hydrodynamic model (SRH-2D; Lai, 2008) of the Rhine river and its eastern floodplain was set up. The model bathymetry is based on a high-resolution (1 m) DEM (Wasserstrassen-und Schifffahrtsverwaltung des Bundes (WSV), 2016) supplemented by longitudinal profiles and cross sections through the main waterbodies in the study area (Díaz-Redondo et al., 2018). The model mesh consists of 149,900 nodes with an average distance of 20 m in the main channel, 10 m in the floodplain and down to 2 m in the river bank and dam zones. Break lines were integrated manually. Water surface elevations (WSE) for 6 flood events with return periods between 1 and 100 years provided by the German Federal Agency for Hydrology (BfG) were used for model calibration and setting of the lower boundary condition. Manning roughness coefficients were first appointed to the model elements based on different land use and lie around 0.083 for the floodplain forest and around 0.026 for the side channels and 0.037 for the main river channel. Calibration was performed by adjusting Manning's roughness coefficients to minimize the difference between modelled and measured WSEs. Mean WSE errors were between 1 cm for flood events with short return periods and 20 cm for higher return periods.

Calibration and validation data
We used an expert-based PNV map (Ochs et al., 2019) and analyzed the historical land-use to delineate likely reference areas for four main vegetation types: Reeds, softwood forest, transition forest and hardwood. Within these areas a total of 130 random sampling plots (radius = 5 m) were distributed with a minimum distance between them of 50 m. The PNV-type of the plots was verified during several field visits (Föll and Egger, 2017). The verification was guided by indicator species from the herb and shrub layers (see Table A1). Reeds could be confirmed in 8 plots, softwood forest in 36 plots, transition forest in 40 plots and 37 plots could be clearly identified as hardwood forest. To increase the independency of the assessment of the predictive performance and comparison between the different modelling approaches we split the study area geographically perpendicular to the river axis (Wenger and Olden, 2012). The downstream part that represents between 30% and 40% of the reference plots of each vegetation type was used for validation ( Fig. 1).

Dynamic succession model (DM)
The dynamic floodplain vegetation model CASIMIR (Benjankar et al., 2011;Egger et al., 2013) was used to predict the PNV by simulating the succession of the floodplain vegetation from 1872 to 2016. The time period was chosen because by 1872 the study area had already suffered the main hydro-morphological impacts through river straightening and channelization (Bernhardt, 2000). In the model, the riparian vegetation is represented in succession lines and their respective succession phases (Table C2). The dynamic modules are: Recruitment, controlled by the spring mean water level as described by the recruitment box model (Mahoney and Rood, 1998) and Succession (Progression/Retrogression), controlled by the disturbance indicators "flood duration" and "shear stress" (Formann et al., 2014). Each year the recruitment module checks for bare soils in the bank and floodplain zone as well as the water levels that allow seedling survival (Table C3) and the disturbance module checks whether the critical values of the disturbance indicators are surpassed (Table C4 and Table C5). The result of one simulated year will be used as input for the next year. The parametrization of the model was based on analyses of historic maps and historic discharge data. The model was calibrated against an expertbased PNV estimation of the upstream part of the study area so that the reference plots for validation can also be considered independent (for a detailed description of the model functioning and calibration/validation see Ochs et al., 2019 as well as Appendix B and C).
For comparability the succession phases of the final year were aggregated to match the main PNV types in the study area (Table C2).

Statistic models (SM1 and SM2)
The statistic modelling approach for the classification of the main PNV types was based on the random forest algorithm (Breiman, 2001). Random forest selects random bootstrap samples from a given dataset to build a set of decision trees. The final prediction is based on the majority vote from the individually developed trees.
We chose three predictors representing the hydrological control factors of riparian vegetation: flood duration, water depth and shear stress. For the flood duration raster, the average flood duration of each grid cell during the growing periods between 1921 and 2016 was calculated (see Appendix B). Maps of water depth and shear stress were calculated for HQ10 (4100 m 3 /s) so that the whole floodplain could be represented. In addition, we tested the influence of the habitat age. The geomorphological age of different areas of the floodplain was derived through the analyses of the changes from water surfaces to sand and gravel bars on historic maps (Diaz-Redondo et al., 2017).
We built two different models: SM1 was based on the hydrological predictors and geomorphological age, SM2 considered only the hydrological predictors. The models were set to grow 1000 trees based on bootstrap samples from the calibration plots. The sample was balanced to compensate for the overrepresentation of softwood forest in the reference plots which according to an expert-based PNV map of 2017 covered around 10% of the study area. The statistical modelling was done in the R environment using the "randomForest" package (Liaw and Wiener, 2002).

Gradient model (GM)
The German Federal Institute of Hydrology (BfG) developed a gradient model for the large-scale assessment of the main floodplain vegetation types for the free-flowing parts of River Rhine and Elbe in Germany. It is based on the field-observation that the occurrence of Salix alba at a site correlates with the relative height to the mean water level and the mean annual flood duration (Schleuter, 2014). The mean annual flood duration of a grid cell is calculated as follows (Schleuter, 2016): Ln X 70.599 ( 0.50) 88.711 F = mean annual flood duration X = relative height to mean water level (m) The PNV types are then assigned based on expert knowledge (Table 1).

Model validation and comparison
For validation of the predictive performance all models were tested against the same set of geographically separated reference plots. Based on a confusion matrix we calculated the global metrics overall accuracy (OA) and kappa coefficient (K) (Cohen, 1960), which corrects the OA for chance agreement. In addition, we calculated Sensitivity and Specificity for each PNV class. For comparison between the models all area wide predictions were directly compared to each other by calculating the metrics Kappa (K), Kappa Location (KLoc) and Kappa histogram (Khist). KLoc describes the similarity of spatial allocation of categories of the two compared maps, and Khist describes the quantitative similarity (Pontius, 2000). The following rating system was applied: values greater than 0.75 indicate very good-to-excellent agreement, values between 0.40 and 0.75 indicate fair-to-good agreement, and values of 0.40 or less indicate poor agreement (Landis and Koch, 1977).

Results
The results of the different approaches to reconstruct the PNV of our study area are shown in Fig. 2. The overall agreement between the models was good. All approaches predicted hardwood forests to be the dominant vegetation class followed by transition forests, softwood forests and reeds (Table 2). Along the same sequence the agreement of the predictions between the approaches diminished (Table 3). Hardwood forests were predicted for about 50% of the study area by all models and the agreement (spatial and quantitative) was excellent to very good. Transition forests were estimated to cover around 35% by DM, SM1 and SM2 but only 26% by GM. The similarity of the predictions of DM and SM1/2 was very good but only fair when compared to GM. Softwood forest were predicted on only 7% of the area by the DM but nearly twice as much by the other models. The agreement between the DM and SM1/2 still can be considered fair but showed high discrepancies to GM. Reeds presented poor agreement between all models, especially spatially.

Validation
Overall the DM, SM1 and SM2 showed good predictive performance, and that of the GM was fair (see Table 4). Notably, SM1 and SM2 performed identically. All models were unable to detect reeds (sensitivity = 0). Softwood forest and hardwood forests were predicted with very good accuracy. But SM1, SM2 and GM only identified about 50% of transition forest reference plots correctly.
The "Mean Decrease Accuracy" and "Mean Decrease Gini" measures of the random forest models both revealed flood duration to be the most Table 1 Main PNV-types, relative height to the mean water level and mean annual flood duration according to the GM (Schleuter, 2016 K. Ochs, et al. Journal of Hydro-environment Research xxx (xxxx) xxx-xxx important predictor. "Habitat age" and "shear stress" were the least important ones.

Discussion
The equally good performances of the dynamic and static modelling approaches in predicting the PNV of our study area support the hypothesis that due to the loss of natural disturbance dynamics the riparian vegetation in our study area has reached a stable equilibrium with the hydrological control factors.
Sensitivity analyses of the statistic model and the DM (Ochs et al., 2019) revealed that the PNV is mainly determined by "flood duration". But we show that the resulting pattern of softwood, transition and hardwood forest is explained equally well by a static average as a reconstruction of the temporal dynamics of the flood regime (DM). Even more, the fair results of the gradient approach show that the relative height to the mean water level also captures most of the influencing factor of riparian vegetation. With the predictor "habitat age", we wanted to include a time dimension to the static modelling approach as indication of a possible successional progression. However, roughly 150 years after geomorphological changes have been impeded habitat age proved to have no influence on the present vegetation communities. The transition of the large-scale dynamic equilibrium of natural floodplain ecosystem to a more mature and stable state after river regulation has also been recognized by other studies (Diaz-Redondo et al., 2017;Hohensinner et al., 2004;Ollero, 2010;Tockner and Stanford, 2002) and has been mainly attributed to an impediment of morphologic dynamics through bank stabilization and flow regulation (Florsheim et al., 2008;Hohensinner et al., 2014).
To allow for the comparison of the model predictions we validated the results against a geographically separated holdout sample. The spatial blocking strategy increases the independency of the sample and allow an effective test of a models transferability (Roberts et al., 2017;Wenger and Olden, 2012). It meant however, a trade-off with the sample size used for calibration of the statistical model which already had to be considered small (Wisz et al., 2008). Nevertheless, random forests are recognized as one of the most accurate species distribution modelling techniques (Cutler et al., 2007;Elith et al., 2006). Also, some confidence about our results can be drawn from the good agreement between the modelling approaches themselves. The reference plots were identified based on indicator species from the herb and shrub layers (see Table A1) that usually develop without direct human manipulation (Gilliam, 2007;Metzger and Schultz, 1984). As opposed to area-wide expert based assessments of PNV that are often used to validate ecological models in areas of high anthropogenic transformation (Hickler et al., 2012;Somodi et al., 2017) the reference plots are not extrapolated and therefore more comprehensive and less prone to uncertainties.
All tested modelling approaches simplify the complex floodplain ecosystem and are based on several assumptions. They assume that in our study area the hydrological control factors are most relevant and neglect other factors that are known to influence plant communities in floodplains. Regarding the occurrence of reeds this seems to be an oversimplification since no model was able to detect it. The proliferation and dominance of Phragmites australis can be linked to nutrient competition and allelopathy (Hazelton et al., 2014;Uddin and Robinson, 2018). The fixed topographic input and disregard of the complex hydro-morphological processes normally occurring within the floodplain (Gurnell, 2016) can be justified in part by river regulation measures and artificially stabilized banks (a further in-depth discussion of the uncertainties regarding the dynamic model Casimir vegetation can be found here: Benjankar et al., 2011;Ochs et al., 2019). Another obvious source of prediction bias for both the statistic and dynamic model are possible errors in the hydrological model and the historic maps that were used for parameterization.
In addition to the predictive performance, other important criteria when choosing a model are the required resources and deployment time. The simulation of succession dynamics for nearly 150 years was only possible with access to data of high spatiotemporal resolution and a high level of expert knowledge as well as a laborious calibration process. The static models on the other hand needed less data, knowhow and time. Especially the very simple gradient model (GM) which still showed fair agreement with the other modelling approaches doesn't even need hydraulic simulations since it is only based on the relative distance to the mean water level.  Table 3 Agreement metrics between the four modelling approaches (green = good-toexcellent; yellow = fair-to-good; grey = poor).

Table 4
Accuracy measures of the four modelling approaches (green = good-to-excellent; yellow = fair-to-good; grey = poor).

Conclusion
The high degree of transformation of large river floodplains through forestry and agriculture makes PNV a valuable concept, particularly as a benchmark for conservation measures. Although the conceptual and methodological issues around PNV are much discussed (Chiarucci et al., 2010;Loidi and Fernández-González, 2012;Somodi et al., 2012) the specific challenges for its reconstruction in river flood plains have gained little attention. Because natural floodplains are a disturbancedriven ecosystem, the classical, static PNV definition is not reasonable. However, through the direct comparison of process-based and statistic modelling approaches for PNV we showed that after 150 years of river regulation and impediment of geomorphological dynamics the riparian vegetation has reached a stable equilibrium state with its hydrologic control factors. A static perception of its PNV seems justified. Consequently, statistical models are the best option for its reconstruction, since they need relatively few resources (data, time, expert knowledge) and are reproducible.

Acknowledgements
The first author is part of the FLUVIO Doctoral Program, supported by a grant from the Fundação para a Ciência e Tecnologia (PD/BD/ 114354/2016).
We acknowledge the collaboration of the Karlsruhe Institute of Technology, Department of Wetland Ecology in Rastatt (Germany).

Appendix B. Calculation of flood duration
Map representations of flood duration were needed as input for the DM and static models. For the DM five representative years (dry, medium wet, wet, very wet and extreme wet) from the period 1921-2016 were selected based on their maximum, mean and minimum discharge and the representativeness of the flow duration curve (Table B1). For the SM the average flood duration of the vegetation period (1921-2016) was considered. To reduce the calculational efforts for the hydrodynamic modelling, 14 discharges were selected to best represent the average flow duration

PNV-Type
Indicator species and forest layer (TL -tree layer, SL -shrub layer, HL -herb layer) Reeds curve of the vegetation period (Table B2). Using the two-dimensional hydrodynamic model, the water surface elevations (WSE) for these 14 discharges were calculated for the whole study area. In a first step the results from the irregular hydrodynamic model mesh were transferred into raster (regular grids). For the calculation of the flood duration raster of each representative year we attributed to the WSE of the 14 modeled discharges the number of days that they were exceeded during the growing period through analyzing the respective hydrographs (Table B2). The final raster was then composed through superimposing the WSE of each modelled discharge and their number of days exceeded. The flood duration of the grid cell located between the water edge lines of two neighboring WSE was calculated according to the relative vertical position of the grid cell between the two calculated water surface elevation, as shown in Fig. B1 183  183  183  183  183  183  1257  MQ  39  112  135  183  147  109  1600  Other  14  63  96  149  111  59  2000  Other  2  8  46  99  66  18  2200  Other  0  5  20  67  53  11  2450  Other  0  2  12  40  39  5  2724  HQ1  0  0  8  13  35  2  2850  HQ2  0  0  7  7  29  2  3000  Other  0  0  6  4  21  1  3150  Other  0  0  5  3  13  1  3594  HQ5  0  0  3  0  7  0  4100  HQ10  0  0  1  0  3  0  4900 HQ50 0 0 0 0 0 0 5300 HQ100 0 0 0 0 0 0 Fig. B1. Illustration for the calculation of the flood duration of a grid cell located between the water edges of two modeled discharges (Eq. (B1)).
K. Ochs, et al. Journal of Hydro-environment Research xxx (xxxx) xxx-xxx Appendix C. CASiMiR model parameters A detailed description of the model functioning as well as a description of the parametrization, calibration and validation of the model can be found in Ochs et al., 2019.

Table C1
Deduction of the starting condition derived from the analysis of the chronological sequence of forest, gravel and water surfaces on historical maps.

Table C2
Succession lines, phases and age spans of the DM. The succession lines and phases were aggregated to match the main PNV types in the study area: IP, PP and HP to "Reeds", SP and ESWP to "Softwood forest", LSWP as "Transition forest" and EFP and TS to "Hardwood forest".

Table C3
Water levels for the bank and floodplain zone that allow recruitment. The maximum water levels of the siltation lines are chosen so high that potential outliers will be included.