A Multiscale Framework for Sustainable Management of Tufa-Forming Watercourses: A Case Study of National Park “Krka”, Croatia

: Tufa sedimentary systems are sensitive ﬂuvial landscapes subject to various external disturbances. Tufa landscape degradation reﬂected in negative hydrological changes and a decrease in the intensity of the tufa formation process have been detected in National Park Krka (Croatia). The main causes were recognized in the uncontrolled spread of invasive vegetation ( Ailanthus altissima ) and increased anthropogenic inﬂuence. Therefore, the Park administration launched the project, Management and Maintenance of Macro-Vegetation at Skradinski Buk (SB)—Development of a Multicriteria Model for Sustainable Management. The methodological framework was divided into three scales of research. The macro-scale research comprised a set of activities aimed at selecting the most suitable test surface within a wider area of the Skradinski Buk (SB) waterfall. The meso-scale research involved mapping the reference and ﬁnal state of the vegetation and hydrological network after the removal of invasive vegetation and mitigation of negative anthropogenic impact. At the micro-scale, a monitoring system was established to track the quality of the tufa sedimentary system. Special emphasis was placed on the measurement of tufa formation dynamics (TFD) on limestone plates using a new methodological approach based on structure from motion (SfM) photogrammetry. Implementation of the proposed multiscale framework resulted in reactivation of tufa-forming watercourses, prevention of invasive vegetation regeneration and achievement of sustainable conditions for the tufa formation process. In reactivated watercourses, the average tufa growth rate was 4.267 mm a − 1 (n = 18). Potential users of this framework include local authorities and administrators of protected areas.


Introduction and Background
Landscapes by their very nature necessitate a multiscale approach in their monitoring, modelling and management [1,2]. One of the most complex, vulnerable and fragile systems in the world is the karst landscape [3,4]. Due to the importance of natural resource preservation and the possibility of their exploitation, interest in the management of karst landscapes is constantly growing [4,5]. One of the most spectacular depositional forms in karst landscapes are tufa [6] and travertine barrages and cascades [7] developed within river and lake systems [8]. They are geomorphological forms of universal scientific and aesthetic values [9] and are often included in the UNESCO World Heritage  Tufa landscape degradation in NPK has been observed in the wider area of SB. Therefore, the macro-level of research refers to the wider area (50 ha) of SB ( Figure 1C). The SB waterfalls are regarded as one of the most valuable and beautiful tufa waterfalls in Europe [36]. This area is influenced by a moderately warm Mediterranean rainy climate (Csa) with dry and hot periods in summer. Most of the annual precipitation (approximately 1200 mm) falls during October-February. At the meso-level of research, the case study area refers to the selected, most suitable, test surface (0.8 ha) within the wider are of SB. It is located on the northeast side near the Marasovića Lake ( Figure 1D). The selection of test surface was performed using multi-criteria GIS analyzes (MCDA-GIS). Finally, the micro-level Water 2020, 12, 3096 4 of 34 of research refers to the previously active and reactivated tufa-forming watercourses, respectively a specific limestone plates (PLs) installed in it.

General Structure of the Methodological Framework
The methodology is based on an interdisciplinary research. This implies the treatment of the perceived problem (tufa landscape degradation) as a whole, taking into account the knowledge from different disciplines. This methodological framework can be regarded as the result of the project implementation. The methodology was divided into macro-, meso-and micro-scale of research following the principle of multi-resolution (hierarchical) landscape modeling ( Figure 2). Understanding the changes in the tufa landscape system requires a multiscale approach because analyses over a uniform scale are inapplicable in this case. Landscape hierarchy theory states that understanding of a studied element comes from the lower hierarchical level and significance from a higher hierarchical level [43]. Therefore, the dynamics of the focal element evolution (tufa formation) are influenced by the sets of elements from lower (e.g., biohydrochemical parameters) and higher (e.g., vegetation density, soil morphology) hierarchical levels. Thus, examining such phenomena requires studying the change in elements at different scales [3,44,45].
Water 2020, 12, x FOR PEER REVIEW 4 of 36

General Structure of the Methodological Framework
The methodology is based on an interdisciplinary research. This implies the treatment of the perceived problem (tufa landscape degradation) as a whole, taking into account the knowledge from different disciplines. This methodological framework can be regarded as the result of the project implementation. The methodology was divided into macro-, meso-and micro-scale of research following the principle of multi-resolution (hierarchical) landscape modeling ( Figure 2). Understanding the changes in the tufa landscape system requires a multiscale approach because analyses over a uniform scale are inapplicable in this case. Landscape hierarchy theory states that understanding of a studied element comes from the lower hierarchical level and significance from a higher hierarchical level [43]. Therefore, the dynamics of the focal element evolution (tufa formation) are influenced by the sets of elements from lower (e.g., biohydrochemical parameters) and higher (e.g., vegetation density, soil morphology) hierarchical levels. Thus, examining such phenomena requires studying the change in elements at different scales [3,44,45]. Therefore, the first (macro-) scale of research includes a set of activities aimed at selecting the most suitable smaller surface area within the wider area of Skradinski Buk (Figure 3), where specific management measures could be implemented. The second (meso-) scale of research includes mapping of the initial and final state of tufa-forming watercourses (active and dried up) and vegetation cover at defined test surface (Figure 3). At that scale, the implementation of specific management measures (selective removal of vegetation and mitigation of negative anthropogenic effects) were done with the aim of dried-up watercourse reactivation. At the last (micro-) scale of research, a system for monitoring the state of the tufa landscape was established. It consisted of four components. The first included tufa formation dynamic (TFD) monitoring on the installed artificial limestone PLs (Figure 3). The second component included systematic monitoring of the flow rate in initially active and reactivated tufa-forming watercourses. The third referred to the interval measurement of key physico-chemical parameters of water. The last, fourth component, included the monitoring of the invasive vegetation regeneration using the cyclic photography method. Basic scheme of landscape hierarchy theory in the context of tufa landscape degradation problem. Therefore, the first (macro-) scale of research includes a set of activities aimed at selecting the most suitable smaller surface area within the wider area of Skradinski Buk (Figure 3), where specific management measures could be implemented. The second (meso-) scale of research includes mapping of the initial and final state of tufa-forming watercourses (active and dried up) and vegetation cover at defined test surface (Figure 3). At that scale, the implementation of specific management measures (selective removal of vegetation and mitigation of negative anthropogenic effects) were done with the aim of dried-up watercourse reactivation. At the last (micro-) scale of research, a system for monitoring the state of the tufa landscape was established. It consisted of four components. The first included tufa formation dynamic (TFD) monitoring on the installed artificial limestone PLs ( Figure 3). The second component included systematic monitoring of the flow rate in initially active and reactivated tufa-forming watercourses. The third referred to the interval measurement of key physico-chemical parameters of water. The last, fourth component, included the monitoring of the invasive vegetation regeneration using the cyclic photography method.

Macro-Scale Research
At the macro-scale, the main research objective was to select the most suitable test surface within the wider area of SB on which specific management measures would be implemented. The most suitable test surface means that this area has a detected tufa landscape degradation problem or that there exists a high probability of this occurrence. Furthermore, this area must be accessible for performing various data acquisition methods (terrestrial LiDAR, UAV photogrammetry, installation of PLs, etc.). Implementation of management measures on the entire SB area was not necessary because tufa landscape degradation did not occur everywhere. Additionally, since these measures have not been applied in other protected areas, it was necessary to test their effectiveness on a smaller test surface. The logistical and financial capabilities of the project were also limiting factors.

Field Surveys
The first step involved a preliminary field survey conducted by an interdisciplinary science team with the members of the Park administration. In the one-month period, several field surveys were performed within the wider area of SB. Potential locations of tufa landscape degradation were photographed and mapped using RTK-GPS the Stonex S10 (Paderno Dugnano, MI, Italy). The most suitable test surface in the wider area of SB was initially identified. However, it was concluded that

Macro-Scale Research
At the macro-scale, the main research objective was to select the most suitable test surface within the wider area of SB on which specific management measures would be implemented. The most suitable test surface means that this area has a detected tufa landscape degradation problem or that there exists a high probability of this occurrence. Furthermore, this area must be accessible for performing various data acquisition methods (terrestrial LiDAR, UAV photogrammetry, installation of PLs, etc.). Implementation of management measures on the entire SB area was not necessary because tufa landscape degradation did not occur everywhere. Additionally, since these measures have not been applied in other protected areas, it was necessary to test their effectiveness on a smaller test surface. The logistical and financial capabilities of the project were also limiting factors.

Field Surveys
The first step involved a preliminary field survey conducted by an interdisciplinary science team with the members of the Park administration. In the one-month period, several field surveys were performed within the wider area of SB. Potential locations of tufa landscape degradation were photographed and mapped using RTK-GPS the Stonex S10 (Paderno Dugnano, MI, Italy). The most suitable test surface in the wider area of SB was initially identified. However, it was concluded that due to the complex terrain morphometry, high vegetation density and available historical data, it was necessary to use GIS-based multicriteria decision analysis (GIS-MCDA) [46] to select the most suitable test surface.

Aerial LiDAR Scanning
Aerial LiDAR scanning of the six waterfalls and both banks of the Krka River was performed for the purposes of generating the criteria for GIS-MCDA, production of referent digital surface model (DSM), mapping of high vegetation and other needs of the Park administration. Aerial LiDAR scanning was performed using a full-waveform laser measuring system Riegl LMS-Q780. The system was connected to a Eurocopter EC120B helicopter with an integrated Hasselblad H39 camera and a Novatel OEV/OEM4 internal GPS receiver. Data were collected in the NW-SE direction of flight. The laser beam stabilization or deviation accuracy was approximately 0.25 mrad, corresponding to a positional error of approximately 25 cm in the field at a relative flight altitude of 1000 m. The flight altitude was 500-600 m, while the flight speed was 112.12 km/h. The point density was 24.7 points/m 2 (precision 20 mm). Data processing was performed using the software packages Microstation (v2004) and Terrasolid, while Grafnav and IGI Aerooffice were used to process the GPS data. The process involved input and point classification, manual correction of automatic point classification (pulse time), data transformation and ultimately the creation of a digital terrain model (DTM). The entire scanned area contained a total of 350 dense point clouds with more than 781 million ground points (last pulse). Classified data were saved in las format (LAS format) and converted into the LAS Dataset structure from which detailed statistics of a classified point cloud were generated. After point classification, the average sampling density dropped to 10-15 points per m 2 , depending on the density of the vegetation cover. A DTM was created using the LAS Dataset to Raster tool in ArcMap software ( Figure 4).
Water 2020, 12, x FOR PEER REVIEW 6 of 36 due to the complex terrain morphometry, high vegetation density and available historical data, it was necessary to use GIS-based multicriteria decision analysis (GIS-MCDA) [46] to select the most suitable test surface.

Aerial LiDAR Scanning
Aerial LiDAR scanning of the six waterfalls and both banks of the Krka River was performed for the purposes of generating the criteria for GIS-MCDA, production of referent digital surface model (DSM), mapping of high vegetation and other needs of the Park administration. Aerial LiDAR scanning was performed using a full-waveform laser measuring system Riegl LMS-Q780. The system was connected to a Eurocopter EC120B helicopter with an integrated Hasselblad H39 camera and a Novatel OEV/OEM4 internal GPS receiver. Data were collected in the NW-SE direction of flight. The laser beam stabilization or deviation accuracy was approximately 0.25 mrad, corresponding to a positional error of approximately 25 cm in the field at a relative flight altitude of 1000 m. The flight altitude was 500-600 m, while the flight speed was 112.12 km/h. The point density was 24.7 points/m 2 (precision 20 mm). Data processing was performed using the software packages Microstation (v2004) and Terrasolid, while Grafnav and IGI Aerooffice were used to process the GPS data. The process involved input and point classification, manual correction of automatic point classification (pulse time), data transformation and ultimately the creation of a digital terrain model (DTM). The entire scanned area contained a total of 350 dense point clouds with more than 781 million ground points (last pulse). Classified data were saved in las format (LAS format) and converted into the LAS Dataset structure from which detailed statistics of a classified point cloud were generated. After point classification, the average sampling density dropped to 10-15 points per m 2 , depending on the density of the vegetation cover. A DTM was created using the LAS Dataset to Raster tool in ArcMap software ( Figure 4).

Application of GIS-MCDA in the Selection of the Test Surface
The GIS-MCDA was used for the selection of the most suitable test surface. Two concepts were used: Boolean and continuous values [47]. The Boolean approach, which serves to limit the alternatives under consideration, relies on the binary statement (true/false) of the criteria for a considered goal or decision [48]. Four exclusionary criteria were used ( Figure 5):

Application of GIS-MCDA in the Selection of the Test Surface
The GIS-MCDA was used for the selection of the most suitable test surface. Two concepts were used: Boolean and continuous values [47]. The Boolean approach, which serves to limit the alternatives under consideration, relies on the binary statement (true/false) of the criteria for a considered goal or decision [48]. Four exclusionary criteria were used ( Figure 5):
Built-up areas;  Three criteria, which are recognized as key in the terrestrial laser scanning process, were used in the continuous value approach: 1. Terrain ruggedness; 2. Slope; 3. Viewshed/visibility Then, the standardization process was performed. Each criterion was ranked into (1-5) classes using the Jenks classification method, where (5) represents very high convenience, (4) high, (3) medium, (2) low and (1) very low convenience ( Figure 6). Pairwise comparison of criteria using the analytical hierarchical process (AHP) was performed to determine the criteria weights. Then, the consistency ratio (CR) was calculated. The determination of weight coefficients is valid if the CR value is < 0.1 or if the error is less than 10%. In this case, the CR was 0.01 (Table 1). This measure show how consistent the mutual comparisons of the criteria are.  Three criteria, which are recognized as key in the terrestrial laser scanning process, were used in the continuous value approach:
Viewshed/visibility Then, the standardization process was performed. Each criterion was ranked into (1-5) classes using the Jenks classification method, where (5) represents very high convenience, (4) high, (3) medium, (2) low and (1) very low convenience ( Figure 6).  Three criteria, which are recognized as key in the terrestrial laser scanning process, were used in the continuous value approach: 1. Terrain ruggedness; 2. Slope; 3. Viewshed/visibility Then, the standardization process was performed. Each criterion was ranked into (1-5) classes using the Jenks classification method, where (5) represents very high convenience, (4) high, (3) medium, (2) low and (1) very low convenience ( Figure 6). Pairwise comparison of criteria using the analytical hierarchical process (AHP) was performed to determine the criteria weights. Then, the consistency ratio (CR) was calculated. The determination of weight coefficients is valid if the CR value is < 0.1 or if the error is less than 10%. In this case, the CR was 0.01 (Table 1). This measure show how consistent the mutual comparisons of the criteria are.  Pairwise comparison of criteria using the analytical hierarchical process (AHP) was performed to determine the criteria weights. Then, the consistency ratio (CR) was calculated. The determination of weight coefficients is valid if the CR value is < 0.1 or if the error is less than 10%. In this case, the CR was 0.01 (Table 1). This measure show how consistent the mutual comparisons of the criteria are. Analysis of historical data about tufa-forming watercourses on wider SB area may indicate the intensity of hydrological and geomorphological changes. The Austro-Hungarian topographic map of the wider area of SB was georeferenced and analyzed ( Figure 7a). The map was created between 1881 and 1884 [49] at a scale of 1:720. The full name of the map was Situations-plan des Krkafall-gebietes bei Scardona (Figure 7b). Tufa-forming watercourses on the map are represented by arrow symbols (Figure 7c). They were mapped with a manual vectorization method (Figure 7c). The collected data served as a basis for the detection of dried tufa-forming watercourses in the wider SB area. Analysis of historical data about tufa-forming watercourses on wider SB area may indicate the intensity of hydrological and geomorphological changes. The Austro-Hungarian topographic map of the wider area of SB was georeferenced and analyzed ( Figure 7a). The map was created between 1881 and 1884 [49] at a scale of 1:720. The full name of the map was Situations-plan des Krkafall-gebietes bei Scardona (Figure 7b). Tufa-forming watercourses on the map are represented by arrow symbols (Figure 7c). They were mapped with a manual vectorization method (Figure 7c). The collected data served as a basis for the detection of dried tufa-forming watercourses in the wider SB area.

Meso-Scale Research
Meso-scale research was conducted on the test surface selected through the applied GIS-MCDA, where the activities include: 1. Mapping the referent state of vegetation and tufa-forming watercourses using terrestrial LiDAR scanning (Z + F Imager5600i), UAV photogrammetry (Phantom 4 Pro) and field survey (RTK-GPS Stonex S10).

Meso-Scale Research
Meso-scale research was conducted on the test surface selected through the applied GIS-MCDA, where the activities include:

2.
Manual and mechanical removal of invasive vegetation and implementation of management measures to mitigate the negative anthropogenic impact on the selected test surface.

Terrestrial LiDAR Scanning
Terrestrial LiDAR scanning of the vegetation and tufa-forming watercourses referent state was performed on 5-7 April 2016. Scanning was performed with the Z+F Imager 5010 TLS (Figure 8a). The scanning resolution was set to high, while the quality was normal. A total of 244,736,910 points ( Figure 8) were recorded at 156 stops (over 1.5 million points per stop). The distribution of scanning stops was the result of the test surface accessibility and its morphometric characteristics. In between two scans, at least two marks were left at the same locations. The marks were placed within the scanning range and moved through the scanning process. The coordinates of the six marks were collected by the RTK-GPS Stonex S10 receiver and used for dense cloud georeferencing.
Terrestrial LiDAR scanning of the vegetation and tufa-forming watercourses referent state was performed on 5-7 April 2016. Scanning was performed with the Z+F Imager 5010 TLS (Figure 8a). The scanning resolution was set to high, while the quality was normal. A total of 244 736 910 points ( Figure  8) were recorded at 156 stops (over 1.5 million points per stop). The distribution of scanning stops was the result of the test surface accessibility and its morphometric characteristics. In between two scans, at least two marks were left at the same locations. The marks were placed within the scanning range and moved through the scanning process. The coordinates of the six marks were collected by the RTK-GPS Stonex S10 receiver and used for dense cloud georeferencing. Terrestrial LiDAR scanning of the final state of vegetation and tufa-forming watercourses was performed two years after the implementation of management measures. It was performed on 19 December 2019. Scanning was performed with the Faro Focus M70 TLS ( Figure 9). There were a total of 22 scanning stops in the final scan. The number of stops was lower than in the first scanning since the Faro Focus M70 TLS was more advanced than the Z+F Imager 5010 and since the vegetation density was lower after selective vegetation removal. In the scanning process, 1/4 resolution (with a point distance of approximately 6 mm at 10 m distance) and 2× quality was used. Additionally, the HDR texture parameter was set. The scanning per stop took 6-7 min.
The collected scans were processed in Faro Scene software. The software allows automatic detection of reference targets and register related scans when the targets are placed within the suggested range. TLS scanning stops were set at a distance of 20-30 ms to enable target registration. Target registration is the process of aligning multiple scans in a parent coordinate system using specific markers (spheres, checker boards and manual targets) with reference positions common between scans. Terrestrial LiDAR scanning of the final state of vegetation and tufa-forming watercourses was performed two years after the implementation of management measures. It was performed on 19 December 2019. Scanning was performed with the Faro Focus M70 TLS ( Figure 9). There were a total of 22 scanning stops in the final scan. The number of stops was lower than in the first scanning since the Faro Focus M70 TLS was more advanced than the Z+F Imager 5010 and since the vegetation density was lower after selective vegetation removal. In the scanning process, 1/4 resolution (with a point distance of approximately 6 mm at 10 m distance) and 2× quality was used. Additionally, the HDR texture parameter was set. The scanning per stop took 6-7 min.

Cyclic UAV Photogrammetry
The objective of cyclic UAV photogrammetry was to create digital orthophotos (DoPs) on which changes in vegetation, watercourses and lakes before and after the conduction of specific management measures would be visible. Cyclic UAV photogrammetry was performed using a UAV Phantom 4 ( Figure 10a). The first aerial recording was performed in March 2017. The second was done in March 2018, shortly after the selective removal of vegetation and placement of the artificial substrates in the reactivated watercourses. Before performing the flight, the IMU system and UAV compass were calibrated ( Figure 10b). Ground control points (GCPs) were collected using the RTK-GPS Stonex S10 The collected scans were processed in Faro Scene software. The software allows automatic detection of reference targets and register related scans when the targets are placed within the suggested range. TLS scanning stops were set at a distance of 20-30 ms to enable target registration. Target registration is the process of aligning multiple scans in a parent coordinate system using specific markers (spheres, checker boards and manual targets) with reference positions common between scans.

Cyclic UAV Photogrammetry
The objective of cyclic UAV photogrammetry was to create digital orthophotos (DoPs) on which changes in vegetation, watercourses and lakes before and after the conduction of specific management measures would be visible. Cyclic UAV photogrammetry was performed using a UAV Phantom 4 ( Figure 10a). The first aerial recording was performed in March 2017. The second was done in March 2018, shortly after the selective removal of vegetation and placement of the artificial substrates in the reactivated watercourses. Before performing the flight, the IMU system and UAV compass were calibrated ( Figure 10b). Ground control points (GCPs) were collected using the RTK-GPS Stonex S10 ( Figure 10c). The test surface was recorded with several Double Grid missions at different heights using Pix4Dcapture software. The image workflow process was completed in Agisoft Metashape 1.5.1 software. UAV compass calibration and (c) data acquisition of orientation points using the GPS Stonex S10.

Mapping of Vegetation Cover and Tufa-Forming Watercourses
In March 2017, a mapping of vegetation cover and tufa-forming watercourses referent state was carried out. The objective was to quantify these landscape elements to develop a plan for invasive vegetation removal and to track hydrological changes after the implementation of management measures. Registered and georeferenced dense point clouds of the test surface, acquired in initial and final terrestrial laser scanning, were used in the vegetation and watercourse mapping process ( Figure  11). Based on registered dense point clouds and acquired TLS images, it was possible to recognize each tree type and active watercourse on the test surface. Trees were manually mapped and classified from LAS dataset data in ArcScene software. The whole process was performed using three software packages: Cloud Compare, Faro Scene software and ArcScene. (b) UAV compass calibration and (c) data acquisition of orientation points using the GPS Stonex S10.

Mapping of Vegetation Cover and Tufa-Forming Watercourses
In March 2017, a mapping of vegetation cover and tufa-forming watercourses referent state was carried out. The objective was to quantify these landscape elements to develop a plan for invasive vegetation removal and to track hydrological changes after the implementation of management measures. Registered and georeferenced dense point clouds of the test surface, acquired in initial and final terrestrial laser scanning, were used in the vegetation and watercourse mapping process ( Figure 11). Based on registered dense point clouds and acquired TLS images, it was possible to recognize each tree type and active watercourse on the test surface. Trees were manually mapped and classified from LAS dataset data in ArcScene software. The whole process was performed using three software packages: Cloud Compare, Faro Scene software and ArcScene.
The dried-up tufa-forming watercourses, whose stream beds were filled with invasive vegetation, could not be detected in this way. They were detected through several steps ( Figure 12). In the ArcMap, a surface runoff network was generated based on hybrid DTM (aerial LiDAR data + terrestrial laser scanning data). Then, the surface runoff network was overlaid with tufa-forming watercourses mapped on an Austro-Hungarian topographic map, providing good insight into where dried tufa-forming watercourses could potentially be found. Ultimately, the derived data (referent state of watercourses and vegetation) were checked in the field. Field research was performed using RTK-GPS Stonex S10 and an analogue map (A0 format) divided into regular coded polygons. On the map, extracted surface runoff network, tufa-forming watercourses mapped on an Austro-Hungarian topographic map and vegetation cover were marked. Detection of dried watercourses was also facilitated by the discovery of (old) "dead" tufa at locations of tufa-forming watercourses mapped on an Austro-Hungarian topographic map. Then, the surface runoff network data were modified and an attribute (active or inactive) was added for each watercourse ( Figure 12). After the removal of the invasive vegetation and mitigation of negative anthropogenic impact, the mapping of the vegetation and the tufa-forming watercourses final state was performed following the same methodology. UAV compass calibration and (c) data acquisition of orientation points using the GPS Stonex S10.

Mapping of Vegetation Cover and Tufa-Forming Watercourses
In March 2017, a mapping of vegetation cover and tufa-forming watercourses referent state was carried out. The objective was to quantify these landscape elements to develop a plan for invasive vegetation removal and to track hydrological changes after the implementation of management measures. Registered and georeferenced dense point clouds of the test surface, acquired in initial and final terrestrial laser scanning, were used in the vegetation and watercourse mapping process ( Figure  11). Based on registered dense point clouds and acquired TLS images, it was possible to recognize each tree type and active watercourse on the test surface. Trees were manually mapped and classified from LAS dataset data in ArcScene software. The whole process was performed using three software packages: Cloud Compare, Faro Scene software and ArcScene. The dried-up tufa-forming watercourses, whose stream beds were filled with invasive vegetation, could not be detected in this way. They were detected through several steps ( Figure 12). In the ArcMap, a surface runoff network was generated based on hybrid DTM (aerial LiDAR data + terrestrial laser scanning data). Then, the surface runoff network was overlaid with tufa-forming watercourses mapped on an Austro-Hungarian topographic map, providing good insight into where dried tufa-forming watercourses could potentially be found. Ultimately, the derived data (referent state of watercourses and vegetation) were checked in the field. Field research was performed using RTK-GPS Stonex S10 and an analogue map (A0 format) divided into regular coded polygons. On the map, extracted surface runoff network, tufa-forming watercourses mapped on an Austro-Hungarian topographic map and vegetation cover were marked. Detection of dried watercourses was also facilitated by the discovery of (old) "dead" tufa at locations of tufa-forming watercourses mapped on an Austro-Hungarian topographic map. Then, the surface runoff network data were modified and an attribute (active or inactive) was added for each watercourse ( Figure 12). After the removal of the invasive vegetation and mitigation of negative anthropogenic impact, the mapping of the vegetation and the tufa-forming watercourses final state was performed following the same methodology. Various methods of invasive vegetation removal have been applied in many protected areas [50][51][52]. Prior to the development of a plan, it was necessary to obtain a permit from the Ministry of the Environment to experimentally remove the vegetation of the invasive A. altissima. This species aggressively reduces habitat biodiversity, reducing the value of natural ecosystems and causing

Development of a Plan for Invasive Vegetation Removal
Various methods of invasive vegetation removal have been applied in many protected areas [50][51][52]. Prior to the development of a plan, it was necessary to obtain a permit from the Ministry of the Environment to experimentally remove the vegetation of the invasive A. altissima. This species aggressively reduces habitat biodiversity, reducing the value of natural ecosystems and causing immeasurable damage to karst ecosystems, nature parks and tourism [53]. Due to the sensitive and protected tufa ecosystem of NPK, it was decided that a manual and mechanical method would be applied in the removal of invasive vegetation. This type of manual and mechanical removal promotes the growth of shoots from the trunk and root shoots. Research has shown that A. altissima cannot withstand long-lived exposure to moist and flooded soils; however, once settled, it can survive on dry soils and tolerate drought conditions [54,55]. With continuous removal of invasive vegetation and reactivation of tufa-forming watercourses, unfavorable conditions for the regeneration and survival of this resistant plant were sought.

Micro-Scale Research
At the micro-level of the research, a monitoring system was established to track the quality of the reactivated tufa-forming watercourses. It consists of four components.
Monitoring of invasive vegetation regeneration after manual and mechanical removal on the test surface (cyclic photography-interval camera).

Measurement of TFD
The TFD on the selected test surface was measured at 14 locations on 28 artificial PLs. The design of PLs was explained in Reference [56]. Sixteen PLs at 8 locations ( Figure 13) were installed in dried tufa-forming watercourses that were reactivated after the execution of specific management measures (removal of invasive vegetation and mitigation of negative anthropogenic impacts) at the meso scale of research.
Water 2020, 12, x FOR PEER REVIEW 13 of 36 reactivation of tufa-forming watercourses, unfavorable conditions for the regeneration and survival of this resistant plant were sought.

Micro-Scale Research
At the micro-level of the research, a monitoring system was established to track the quality of the reactivated tufa-forming watercourses. It consists of four components.
1. Measurement of TFD on artificial PLs using structure from motion (SfM) photogrammetry (NIKON D5300 + macro lens LAOWA) 2. Measurement of water flow velocity in reactivated tufa-forming watercourses (current meter digitizer, Rickly Hydrological Co., Columbus, OH, USA) 3. Measurement of physico-chemical parameters in reactivated tufa-forming watercourses (YSI EXO2 multiparameter) 4. Monitoring of invasive vegetation regeneration after manual and mechanical removal on the test surface (cyclic photography-interval camera).

Measurement of TFD
The TFD on the selected test surface was measured at 14 locations on 28 artificial PLs. The design of PLs was explained in Reference [56]. Sixteen PLs at 8 locations ( Figure 13) were installed in dried tufa-forming watercourses that were reactivated after the execution of specific management measures (removal of invasive vegetation and mitigation of negative anthropogenic impacts) at the meso scale of research. Each location was given a name and a specific code engraved on each base of the plate. The plates were placed in different fluvial environments. Table 2 gives detailed information about the locations. Measurement of TFD based on SfM photogrammetry was performed with a newly designed photogrammetric measurement system called the coordinate measuring macrophotogrammetry device (CMD). A detailed description of the CMD construction, components of the device, data collection, image workflow process and calculation of tufa growth and erosion rates are  Each location was given a name and a specific code engraved on each base of the plate. The plates were placed in different fluvial environments. Table 2 gives detailed information about the locations. Measurement of TFD based on SfM photogrammetry was performed with a newly designed photogrammetric measurement system called the coordinate measuring macro-photogrammetry device (CMD). A detailed description of the CMD construction, components of the device, data collection, image workflow process and calculation of tufa growth and erosion rates are given in Reference [56]. However, TFD can be measured using others direct methods of measurement as modified micro-erosion meter (MEM), mass increment method, accretion pins and so forth. After the implementation of management measures, the water flow velocity was measured at seven locations. They were selected to cover all major water inflows to the test surface (P1-P4), key surface water flows (P5-P6) and runoff water from the surface (P7) to the old water well at the dried waterfall ( Figure 14). Flow velocity was measured with a current meter digitizer, Rickly Hydrological Co., Although it would be ideal to measure the water flow velocity rate at individual plate locations, device logistical limitations precluded this operation. Before the implementation of management measures, there was no water at sites P1, P5, P6 and P7. These were dried-up tufa-forming watercourses. The water flow velocity was measured several times a month.

Measurement of Water Quality Using the YSI EXO2 Multiparameter Probe
Systematic monitoring of water quality helps in identifying specific natural processes (e.g., natural eutrophication) and human impacts (e.g., the presence of pollutants or heavy metals) on this sensitive ecosystem. Near the P2 location, a multiparameter YSI EXO2 probe was installed on 18 May 2018. It is designed for water quality monitoring and characterized by high-precision sensors with built-in memory, wireless connectivity and easy integration into marine, freshwater and groundwater systems. The probe was positioned below the walkway and continuously measured the quality of the water coming to the test surface. To reduce the risk of theft, the probe was placed in a metal cage of nonoxidizing iron at a depth of approximately 40 cm at location P3 ( Figure 15). The probe measurement interval was set to 15 min. Ninety-six data logs were collected daily. Data processing was performed in three steps ( Figure 16). First, the data were grouped by day and the mean, maximum, minimum values and standard deviations (SDs) were calculated for each day. Then, derived data by day (mean value) were grouped for a specific month. Then, the maximum, mean, minimum

Measurement of Water Quality Using the YSI EXO2 Multiparameter Probe
Systematic monitoring of water quality helps in identifying specific natural processes (e.g., natural eutrophication) and human impacts (e.g., the presence of pollutants or heavy metals) on this sensitive ecosystem. Near the P2 location, a multiparameter YSI EXO2 probe was installed on 18 May 2018. It is designed for water quality monitoring and characterized by high-precision sensors with built-in memory, wireless connectivity and easy integration into marine, freshwater and groundwater systems. The probe was positioned below the walkway and continuously measured the quality of the water coming to the test surface. To reduce the risk of theft, the probe was placed in a metal cage of non-oxidizing iron at a depth of approximately 40 cm at location P3 ( Figure 15). Systematic monitoring of water quality helps in identifying specific natural processes (e.g., natural eutrophication) and human impacts (e.g., the presence of pollutants or heavy metals) on this sensitive ecosystem. Near the P2 location, a multiparameter YSI EXO2 probe was installed on 18 May 2018. It is designed for water quality monitoring and characterized by high-precision sensors with built-in memory, wireless connectivity and easy integration into marine, freshwater and groundwater systems. The probe was positioned below the walkway and continuously measured the quality of the water coming to the test surface. To reduce the risk of theft, the probe was placed in a metal cage of nonoxidizing iron at a depth of approximately 40 cm at location P3 ( Figure 15). The probe measurement interval was set to 15 min. Ninety-six data logs were collected daily. Data processing was performed in three steps ( Figure 16). First, the data were grouped by day and the mean, maximum, minimum values and standard deviations (SDs) were calculated for each day. Then, derived data by day (mean value) were grouped for a specific month. Then, the maximum, mean, minimum  The probe measurement interval was set to 15 min. Ninety-six data logs were collected daily. Data processing was performed in three steps ( Figure 16). First, the data were grouped by day and the mean, maximum, minimum values and standard deviations (SDs) were calculated for each day. Then, derived data by day (mean value) were grouped for a specific month. Then, the maximum, mean, minimum values and SD of each parameter were calculated for each month. In the third step, the annual trends of the mean values of the selected physico-chemical data, grouped by months, were determined.
Water 2020, 12, x FOR PEER REVIEW 16 of 36 values and SD of each parameter were calculated for each month. In the third step, the annual trends of the mean values of the selected physico-chemical data, grouped by months, were determined. Figure 16. Schematic representation of physico-chemical parameter data processing.

Monitoring of Invasive Vegetation Regeneration
An interval HD hunting camera was mounted on the test surface after the removal of invasive vegetation and reactivation of the tufa-forming watercourses. The camera was hidden and fixed to a tree. It was set up in a location that captures one of the reactivated tufa-forming watercourses (Vrtnjak) and the area where the A. altissima was highly represented before invasive vegetation removal ( Figure 17). Interval image capture was set to a time resolution of 15 min for a period of one year. From interval images, video was created with the aim of invasive vegetation regeneration monitoring.

Monitoring of Invasive Vegetation Regeneration
An interval HD hunting camera was mounted on the test surface after the removal of invasive vegetation and reactivation of the tufa-forming watercourses. The camera was hidden and fixed to a tree. It was set up in a location that captures one of the reactivated tufa-forming watercourses (Vrtnjak) and the area where the A. altissima was highly represented before invasive vegetation removal ( Figure 17). Interval image capture was set to a time resolution of 15 min for a period of one year. From interval images, video was created with the aim of invasive vegetation regeneration monitoring. values and SD of each parameter were calculated for each month. In the third step, the annual trends of the mean values of the selected physico-chemical data, grouped by months, were determined. Figure 16. Schematic representation of physico-chemical parameter data processing.

Monitoring of Invasive Vegetation Regeneration
An interval HD hunting camera was mounted on the test surface after the removal of invasive vegetation and reactivation of the tufa-forming watercourses. The camera was hidden and fixed to a tree. It was set up in a location that captures one of the reactivated tufa-forming watercourses (Vrtnjak) and the area where the A. altissima was highly represented before invasive vegetation removal ( Figure 17). Interval image capture was set to a time resolution of 15 min for a period of one year. From interval images, video was created with the aim of invasive vegetation regeneration monitoring.

Field Survey Results
Field surveys conducted in the wider area of SB with the guidance of the official expert staff of NPK have found that invasive vegetation, more specifically, A. altissima, has spread over a large area of the NPK (Figure 18a,b). [57] were the first to mention this problem.
Water 2020, 12, x FOR PEER REVIEW 17 of 36

Field Survey Results
Field surveys conducted in the wider area of SB with the guidance of the official expert staff of NPK have found that invasive vegetation, more specifically, A. altissima, has spread over a large area of the NPK (Figure 18a,b). [57] were the first to mention this problem. Furthermore, the locations where the tufa formation process was stopped (old/dead tufa) ( Figure  19b) or its intensity was reduced, along with dried-up tufa-forming watercourses and waterfalls, were detected. Additionally, locations with a change in direction and a slowdown in the flow velocity were detected (Figure 19a).  Furthermore, the locations where the tufa formation process was stopped (old/dead tufa) (Figure 19b) or its intensity was reduced, along with dried-up tufa-forming watercourses and waterfalls, were detected. Additionally, locations with a change in direction and a slowdown in the flow velocity were detected (Figure 19a).

Field Survey Results
Field surveys conducted in the wider area of SB with the guidance of the official expert staff of NPK have found that invasive vegetation, more specifically, A. altissima, has spread over a large area of the NPK (Figure 18a,b). [57] were the first to mention this problem. Furthermore, the locations where the tufa formation process was stopped (old/dead tufa) ( Figure  19b) or its intensity was reduced, along with dried-up tufa-forming watercourses and waterfalls, were detected. Additionally, locations with a change in direction and a slowdown in the flow velocity were detected (Figure 19a).  Additionally, many potential sites have been identified that may indicate a possible natural eutrophication process (Figure 19c). In addition, a smaller number of trees were detected that fell into the tufa-forming watercourses due to weather conditions (storms) partially blocking the water flow in these watercourses (Figure 20). Additionally, many potential sites have been identified that may indicate a possible natural eutrophication process (Figure 19c). In addition, a smaller number of trees were detected that fell into the tufa-forming watercourses due to weather conditions (storms) partially blocking the water flow in these watercourses ( Figure 20).   At some locations, negative hydrodynamic phenomena were also observed as a result of anthropogenic impact (trail and resort construction, tourist influence, gravel accumulation under trail/bridges) ( Figure 21). Additionally, many potential sites have been identified that may indicate a possible natural eutrophication process (Figure 19c). In addition, a smaller number of trees were detected that fell into the tufa-forming watercourses due to weather conditions (storms) partially blocking the water flow in these watercourses ( Figure 20). At some locations, negative hydrodynamic phenomena were also observed as a result of anthropogenic impact (trail and resort construction, tourist influence, gravel accumulation under trail/bridges) ( Figure 21).

Aerial LiDAR Scanning Results
Aerial scanning was performed with the Riegl LMS-Q780 full-wave laser measuring system. The average point density was 24.7 points/m 2 . The process of data classification was performed based on reflection time. The ground class included 350 dense point clouds for the entire NPK with a total of just over 781 million points. The spatial resolution of the created DTM was 30 cm. From the DTM of the whole NPK, a model of the wider SB area (approximately 50 ha) was extracted. The extracted DTM served ( Figure 22) as the basis for deriving the morphometric criteria used in the GIS-MCDA.
Water 2020, 12, x FOR PEER REVIEW 19 of 36

Aerial LiDAR Scanning Results
Aerial scanning was performed with the Riegl LMS-Q780 full-wave laser measuring system. The average point density was 24.7 points/m 2 . The process of data classification was performed based on reflection time. The ground class included 350 dense point clouds for the entire NPK with a total of just over 781 million points. The spatial resolution of the created DTM was 30 cm. From the DTM of the whole NPK, a model of the wider SB area (approximately 50 ha) was extracted. The extracted DTM served (Figure 22) as the basis for deriving the morphometric criteria used in the GIS-MCDA.

Tufa-Forming Watercourses on the Austro-Hungarian Topographic Map
The Austro-Hungarian topographic map (Situations-plan des Krkafall-gebietes bei Scardon) was georeferenced. Watercourses in the wider area of SB ten years prior to the construction of the hydropower plant (HE) Jaruga were mapped. Figure 23 shows the georeferenced Austro-Hungarian topographic map with mapped water bodies (watercourses and lakes). This area is characterized by a very dense tufa-forming watercourse network that connects larger active lakes (Šupuk Lake, Golubinka, Veliki vir, Veliko Lake, Draga, Dedinac, etc.).

Tufa-Forming Watercourses on the Austro-Hungarian Topographic Map
The Austro-Hungarian topographic map (Situations-plan des Krkafall-gebietes bei Scardon) was georeferenced. Watercourses in the wider area of SB ten years prior to the construction of the hydropower plant (HE) Jaruga were mapped. Figure 23 shows the georeferenced Austro-Hungarian topographic map with mapped water bodies (watercourses and lakes). This area is characterized by a very dense tufa-forming watercourse network that connects larger active lakes (Šupuk Lake, Golubinka, Veliki vir, Veliko Lake, Draga, Dedinac, etc.). The suitability index for test surface selection was generated using GIS-MCDA. The index value ranged from 0 (least suitable area) to 1 (most suitable area). Values were classified by the equal interval method in the 5 classes. For each class, the descriptive attribute was assigned as follows: (1) very low suitability, (2) low, (3) medium, (4) high and (5) very high suitability. From the derived data,

Selected Test Surface Using GIS-MCDA
The suitability index for test surface selection was generated using GIS-MCDA. The index value ranged from 0 (least suitable area) to 1 (most suitable area). Values were classified by the equal interval method in the 5 classes. For each class, the descriptive attribute was assigned as follows: (1) very low suitability, (2) low, (3) medium, (4) high and (5) very high suitability. From the derived data, a potential test surface was detected ( Figure 24). The test surface was proposed for the following two reasons.
(1) The area was surrounded by a path that makes it easier to perform TLS and (2) the selected surface (approximately 0.8 ha) fell almost entirely in classes of high (4) or very high suitability (5). In addition, the vast majority of examples of the changes in tufa-forming watercourses due to the uncontrolled growth of invasive vegetation and anthropogenic influences listed in Section 3.1.1 were located in that area. However, after the selection of the test surface, further field surveys were conducted to detect tufa-forming watercourses vectorized on the Austro-Hungarian topographic map. However, the detected tufa-forming watercourses were inactive on the test surface. Therefore, taking into account field surveys, GIS-MCDA and mapped tufa-forming watercourses on the Austro-Hungarian map, a test surface north of Marasović Lake was selected as the most suitable for further activities (mesoand micro-scale research).

Initial State of Vegetation Cover and Tufa-Forming Watercourses
The surface runoff model was generated using a hydrological analysis from hybrid DTM (aeroLiDAR + TLS). It was then corrected and modified by field survey. Figure 25 shows the reference state of the surface runoff model before the intervention (removal of invasive vegetation and mitigation of the negative anthropogenic effects). The active flows are located mostly on the southeast side of the test surface ( Figure 25). Many inactive, dried-up tufa-forming watercourses were identified as active streams on the Austro-Hungarian topographic map. However, in an interview with the expert park's staff, it was revealed that some of the tufa-forming watercourses and the waterfall below the NE part of the test surface had dried up several years earlier.
On the test surface, 358 trees were mapped from 15 different species (Figure 25). Along with black hornbeam (Ostrya carpinifolia), the most represented species on the test surface was the tree of heaven (A. altissima) with 20.11%. If all categories of the tree of heaven by age are considered, this species becomes the most represented vegetation on the test surface. For the purpose of restoring the sustainable tufa formation condition, it was necessary to implement specific management measures (removal of invasive vegetation and mitigation of negative anthropogenic impact).

Initial State of Vegetation Cover and Tufa-Forming Watercourses
The surface runoff model was generated using a hydrological analysis from hybrid DTM (aeroLiDAR + TLS). It was then corrected and modified by field survey. Figure 25 shows the reference state of the surface runoff model before the intervention (removal of invasive vegetation and mitigation of the negative anthropogenic effects). The active flows are located mostly on the southeast side of the test surface ( Figure 25). Many inactive, dried-up tufa-forming watercourses were identified as active streams on the Austro-Hungarian topographic map. However, in an interview with the expert park's staff, it was revealed that some of the tufa-forming watercourses and the waterfall below the NE part of the test surface had dried up several years earlier.

Removal of Invasive Vegetation and Mitigation of Anthropogenic Impact
The removal of invasive vegetation was carried out by the expert staff of NPK in June 2017, which, along with July, represents the best period since root food reserves are used to grow shoots and leaves [58]. The following approaches were used: All young seeds and A. altissima trees were removed from the test surface. Specifically, 72 trees with a diameter larger than 15 cm were removed from the test surface. The removal of invasive vegetation and fallen accumulated plant fragments from the dried-up flows resulted in micromorphological changes on the test surface. Furthermore, the removal of artificially accumulated

Removal of Invasive Vegetation and Mitigation of Anthropogenic Impact
The removal of invasive vegetation was carried out by the expert staff of NPK in June 2017, which, along with July, represents the best period since root food reserves are used to grow shoots and leaves [58]. The following approaches were used:

Removal of Invasive Vegetation and Mitigation of Anthropogenic Impact
The removal of invasive vegetation was carried out by the expert staff of NPK in June 2017, which, along with July, represents the best period since root food reserves are used to grow shoots and leaves [58]. The following approaches were used: All young seeds and A. altissima trees were removed from the test surface. Specifically, 72 trees with a diameter larger than 15 cm were removed from the test surface. The removal of invasive vegetation and fallen accumulated plant fragments from the dried-up flows resulted in micromorphological changes on the test surface. Furthermore, the removal of artificially accumulated All young seeds and A. altissima trees were removed from the test surface. Specifically, 72 trees with a diameter larger than 15 cm were removed from the test surface. The removal of invasive vegetation and fallen accumulated plant fragments from the dried-up flows resulted in micromorphological changes on the test surface. Furthermore, the removal of artificially accumulated gravel was performed at those locations where its influence on the change in the flow outflow direction occurred (Figure 27).
Water 2020, 12, x FOR PEER REVIEW 23 of 36 gravel was performed at those locations where its influence on the change in the flow outflow direction occurred (Figure 27). Additionally, the path was raised in those locations where it slowed down the flow and influenced the change in its direction ( Figure 28). Given the small ruggedness of the terrain, these changes directly influenced the alteration of the hydrological network.  Additionally, the path was raised in those locations where it slowed down the flow and influenced the change in its direction ( Figure 28). Given the small ruggedness of the terrain, these changes directly influenced the alteration of the hydrological network.
Water 2020, 12, x FOR PEER REVIEW 23 of 36 gravel was performed at those locations where its influence on the change in the flow outflow direction occurred (Figure 27). Additionally, the path was raised in those locations where it slowed down the flow and influenced the change in its direction ( Figure 28). Given the small ruggedness of the terrain, these changes directly influenced the alteration of the hydrological network.  The combination of intense rainfall with increased water flow (gravel removed under the trail/bridge) caused the process of washing away plant fragments and reactivating the dried tufa, forming watercourses (Video S1) ( Figure 29). These watercourses were initially mapped on the 19th-century Austro-Hungarian topographic map. The combination of intense rainfall with increased water flow (gravel removed under the trail/bridge) caused the process of washing away plant fragments and reactivating the dried tufa, forming watercourses (Video S1) ( Figure 29). These watercourses were initially mapped on the 19thcentury Austro-Hungarian topographic map. Figure 29. Screen shot from Video S1, which shows reactivation of tufa-forming watercourses.

Final State of Vegetation Cover and Tufa-Forming Watercourses
After the implementation of management measures (removal of invasive vegetation, mitigation of negative anthropogenic effects), a new mapping of vegetation and tufa-forming watercourses was performed. Figure 30 shows the final state of vegetation cover and tufa-forming watercourses. All watercourses (old-active and new-reactivated) were given names-Vrtnjak, Crni grab, Ivy, Matica and Smokovik. On the test surface, one dried-up lake (Lake of six trees) was reactivated. During the colder part of the year, water spills from reactivated tufa-forming watercourses and lakes, all over the test surface, were recorded.  . Screen shot from Video S1, which shows reactivation of tufa-forming watercourses.

Final State of Vegetation Cover and Tufa-Forming Watercourses
After the implementation of management measures (removal of invasive vegetation, mitigation of negative anthropogenic effects), a new mapping of vegetation and tufa-forming watercourses was performed. Figure 30 shows the final state of vegetation cover and tufa-forming watercourses. All watercourses (old-active and new-reactivated) were given names-Vrtnjak, Crni grab, Ivy, Matica and Smokovik. On the test surface, one dried-up lake (Lake of six trees) was reactivated. During the colder part of the year, water spills from reactivated tufa-forming watercourses and lakes, all over the test surface, were recorded. The combination of intense rainfall with increased water flow (gravel removed under the trail/bridge) caused the process of washing away plant fragments and reactivating the dried tufa, forming watercourses (Video S1) ( Figure 29). These watercourses were initially mapped on the 19thcentury Austro-Hungarian topographic map. Figure 29. Screen shot from Video S1, which shows reactivation of tufa-forming watercourses.

Final State of Vegetation Cover and Tufa-Forming Watercourses
After the implementation of management measures (removal of invasive vegetation, mitigation of negative anthropogenic effects), a new mapping of vegetation and tufa-forming watercourses was performed. Figure 30 shows the final state of vegetation cover and tufa-forming watercourses. All watercourses (old-active and new-reactivated) were given names-Vrtnjak, Crni grab, Ivy, Matica and Smokovik. On the test surface, one dried-up lake (Lake of six trees) was reactivated. During the colder part of the year, water spills from reactivated tufa-forming watercourses and lakes, all over the test surface, were recorded.

Tufa Growth in Reactivated Tufa-Forming Watercourses
The production of digital submillimeter resolution (0.0236 mm) models of tufa surfaces is explained in Reference [56]. On all PLs installed in reactivated tufa-forming watercourses during a one-year period, the occurrence of tufa was measured ( Figure 31).
Water 2020, 12, x FOR PEER REVIEW 25 of 36

Tufa Growth in Reactivated Tufa-Forming Watercourses
The production of digital submillimeter resolution (0.0236 mm) models of tufa surfaces is explained in Reference [56]. On all PLs installed in reactivated tufa-forming watercourses during a one-year period, the occurrence of tufa was measured ( Figure 31). The tufa growth rate (TGR), as expected, varied significantly by location. The average TGR per PL (n = 16) and by location (n = 8) was calculated. In both cases, the TGR was 4.267 mm a −1 . The highest TGR was recorded at location Near young ash (6.770 mm a −1 ) and on PL12 (6.790 mm a −1 ) ( Figure 32). The SDs of 2.076 mm a −1 by plates and 2.016 mm a −1 by locations were measured. As expected, the lowest TGR was measured at the Lake of six trees (0.327 mm a −1 ). Prior to the implementation of specific management measures, this was a smaller drained lake. According to the defined fluvial classification, this was a location of stagnant water (SW) where plates were placed to a depth > 100 cm; splashing and water turbulence was absent and thus no significant mechanical release of CO2 occurred. In contrast, the largest TGR was recorded at the Near young ash (6.770 mm a −1 ) and Hard tufa (6.659 mm a −1 ) sites. The sites were located in the reactivated tufa-forming stream Vrtnjak in fastflowing water (FFW), where significant water turbulence and splashing facilitated the mechanical release of CO2 and led to CaCO3 precipitation. Furthermore, at measuring surface of PL13 (Hard tufa), large number of tufa-forming organisms have been detected ( Figure 33A,C). On the measuring surface of PL10 (Near young ash) incrusted leaf was detected which is probable reason for the higher tufa growth rate at that PL ( Figure 33B). The tufa growth rate (TGR), as expected, varied significantly by location. The average TGR per PL (n = 16) and by location (n = 8) was calculated. In both cases, the TGR was 4.267 mm a −1 . The highest TGR was recorded at location Near young ash (6.770 mm a −1 ) and on PL12 (6.790 mm a −1 ) ( Figure 32).

Tufa Growth in Reactivated Tufa-Forming Watercourses
The production of digital submillimeter resolution (0.0236 mm) models of tufa surfaces is explained in Reference [56]. On all PLs installed in reactivated tufa-forming watercourses during a one-year period, the occurrence of tufa was measured ( Figure 31). The tufa growth rate (TGR), as expected, varied significantly by location. The average TGR per PL (n = 16) and by location (n = 8) was calculated. In both cases, the TGR was 4.267 mm a −1 . The highest TGR was recorded at location Near young ash (6.770 mm a −1 ) and on PL12 (6.790 mm a −1 ) ( Figure 32). The SDs of 2.076 mm a −1 by plates and 2.016 mm a −1 by locations were measured. As expected, the lowest TGR was measured at the Lake of six trees (0.327 mm a −1 ). Prior to the implementation of specific management measures, this was a smaller drained lake. According to the defined fluvial classification, this was a location of stagnant water (SW) where plates were placed to a depth > 100 cm; splashing and water turbulence was absent and thus no significant mechanical release of CO2 occurred. In contrast, the largest TGR was recorded at the Near young ash (6.770 mm a −1 ) and Hard tufa (6.659 mm a −1 ) sites. The sites were located in the reactivated tufa-forming stream Vrtnjak in fastflowing water (FFW), where significant water turbulence and splashing facilitated the mechanical release of CO2 and led to CaCO3 precipitation. Furthermore, at measuring surface of PL13 (Hard tufa), large number of tufa-forming organisms have been detected ( Figure 33A,C). On the measuring surface of PL10 (Near young ash) incrusted leaf was detected which is probable reason for the higher tufa growth rate at that PL ( Figure 33B). The SDs of 2.076 mm a −1 by plates and 2.016 mm a −1 by locations were measured. As expected, the lowest TGR was measured at the Lake of six trees (0.327 mm a −1 ). Prior to the implementation of specific management measures, this was a smaller drained lake. According to the defined fluvial classification, this was a location of stagnant water (SW) where plates were placed to a depth > 100 cm; splashing and water turbulence was absent and thus no significant mechanical release of CO 2 occurred. In contrast, the largest TGR was recorded at the Near young ash (6.770 mm a −1 ) and Hard tufa (6.659 mm a −1 ) sites. The sites were located in the reactivated tufa-forming stream Vrtnjak in fast-flowing water (FFW), where significant water turbulence and splashing facilitated the mechanical release of CO 2 and led to CaCO 3 precipitation. Furthermore, at measuring surface of PL13 (Hard tufa), large number of tufa-forming organisms have been detected ( Figure 33A,C). On the measuring surface of PL10 (Near young ash) incrusted leaf was detected which is probable reason for the higher tufa growth rate at that PL ( Figure 33B). The average difference in TGRs between adjacent PLs in reactivated tufa-forming watercourses was 0.313 mm a −1 .Furthermore, on all PLs placed in previously active tufa-forming watercourses during a one-year period, the occurrence of tufa formation was measured (Figure 34). At these locations, the TGR varied more than that measured in reactivated tufa-forming watercourses. Specifically, SDs of 8,108 mm a −1 by plates and 6.198 mm a −1 by locations were measured. The highest TGR was recorded at the location Ducks waterfall (19.302 mm a −1 ) and on PL07 (32.293 mm a −1 ) ( Figure  34). As expected, the lowest TGR was measured at the Cascade (1.474 mm a −1 ). This was the location of slow flow velocity (SFW) where PLs were placed to a depth of approximately 30 cm, splashing and water turbulence was absent and thus no significant mechanical release of CO2 occurred. The largest TGR was measured at the Ducks waterfall site, where significant variability in TGR between adjacent PLs (25.981 mm a −1 ) was found. The average difference in TGR between adjacent PLs in previously active tufa-forming watercourses was 4.755 mm a −1 . The average difference in TGRs between adjacent PLs in reactivated tufa-forming watercourses was 0.313 mm a −1 .Furthermore, on all PLs placed in previously active tufa-forming watercourses during a one-year period, the occurrence of tufa formation was measured (Figure 34). At these locations, the TGR varied more than that measured in reactivated tufa-forming watercourses. Specifically, SDs of 8108 mm a −1 by plates and 6.198 mm a −1 by locations were measured. The highest TGR was recorded at the location Ducks waterfall (19.302 mm a −1 ) and on PL07 (32.293 mm a −1 ) ( Figure 34). As expected, the lowest TGR was measured at the Cascade (1.474 mm a −1 ). This was the location of slow flow velocity (SFW) where PLs were placed to a depth of approximately 30 cm, splashing and water turbulence was absent and thus no significant mechanical release of CO 2 occurred. The largest TGR was measured at the Ducks waterfall site, where significant variability in TGR between adjacent PLs (25.981 mm a −1 ) was found. The average difference in TGR between adjacent PLs in previously active tufa-forming watercourses was 4.755 mm a −1 .

Regeneration of Invasive Vegetation
No regeneration of the invasive vegetation on the test surface was observed from the interval photographs collected during one-year recording with a 15 min interval. Furthermore, no regeneration was observed through regular monthly field surveys of the test surface. A likely reason for this is a reactivation of the tufa-forming watercourses with an intensified inflow of water beneath the raised path. This caused occasional flooding of larger areas within the test surface (Figure 35). Such a flooded environment did not favor the regeneration of invasive vegetation since A. altissima, does not tolerate long exposure to moist and flooded soils [54,55].

Flow Velocity in Reactivated Tufa-Forming Watercourses
The flow velocity of the test surface was monitored at seven locations (P1-P7) ( Figure 36). Sites P1, P5, P6 and P7 before the intervention had no water (dried-up tufa-forming watercourses). In the one-year period, the fastest water flows were (P1) 0.897 m/s (major water inflow to the test surface) and (P7) 0.808 m/s (runoff water from the test surface). In October and December 2018, the flow velocity could not be measured due to device malfunction. During the one-year period, there was no

Regeneration of Invasive Vegetation
No regeneration of the invasive vegetation on the test surface was observed from the interval photographs collected during one-year recording with a 15 min interval. Furthermore, no regeneration was observed through regular monthly field surveys of the test surface. A likely reason for this is a reactivation of the tufa-forming watercourses with an intensified inflow of water beneath the raised path. This caused occasional flooding of larger areas within the test surface ( Figure 35). Such a flooded environment did not favor the regeneration of invasive vegetation since A. altissima, does not tolerate long exposure to moist and flooded soils [54,55].

Regeneration of Invasive Vegetation
No regeneration of the invasive vegetation on the test surface was observed from the interval photographs collected during one-year recording with a 15 min interval. Furthermore, no regeneration was observed through regular monthly field surveys of the test surface. A likely reason for this is a reactivation of the tufa-forming watercourses with an intensified inflow of water beneath the raised path. This caused occasional flooding of larger areas within the test surface ( Figure 35). Such a flooded environment did not favor the regeneration of invasive vegetation since A. altissima, does not tolerate long exposure to moist and flooded soils [54,55].

Flow Velocity in Reactivated Tufa-Forming Watercourses
The flow velocity of the test surface was monitored at seven locations (P1-P7) ( Figure 36). Sites P1, P5, P6 and P7 before the intervention had no water (dried-up tufa-forming watercourses). In the one-year period, the fastest water flows were (P1) 0.897 m/s (major water inflow to the test surface) and (P7) 0.808 m/s (runoff water from the test surface). In October and December 2018, the flow velocity could not be measured due to device malfunction. During the one-year period, there was no

Flow Velocity in Reactivated Tufa-Forming Watercourses
The flow velocity of the test surface was monitored at seven locations (P1-P7) ( Figure 36). Sites P1, P5, P6 and P7 before the intervention had no water (dried-up tufa-forming watercourses). In the one-year period, the fastest water flows were (P1) 0.897 m/s (major water inflow to the test surface) and (P7) 0.808 m/s (runoff water from the test surface). In October and December 2018, the flow velocity could not be measured due to device malfunction. During the one-year period, there was no drying up of measurement sites (P1, P5, P6 and P7). At site P1, the highest flow velocity was recorded (1.51 m/s) in April 2018. The lowest values were recorded at sites P4 (0.103 m/s), P3 (0.128 m/s) and P2 (0.323 m/s), which are tufa-forming watercourses coming from Marasović Lake at the test surface. Note that at the initial field surveys (in the macro-level research), it was already recognized that there is a real danger of these streams drying out due to uncontrolled vegetation growth and consequently a change in the direction of runoff (Figure 17a). As expected, at all locations, the lowest flow velocities were measured in July and August. Note that at the initial field surveys (in the macro-level research), it was already recognized that there is a real danger of these streams drying out due to uncontrolled vegetation growth and consequently a change in the direction of runoff ( Figure 17a). As expected, at all locations, the lowest flow velocities were measured in July and August.  Table 3 shows the grouped data for 14 parameters that were collected during a one-year period (May 2018-May 2019). In the analysis, only four parameters were used (percentage of dissolved oxygen at saturation, ODO%, in mg/L; temp in °C; pH; and oxidation-reduction potential, ORP, in mV). The reason was tufa formation around the probe sensors, which resulted in a large variety of data for the following parameters-SpCond in µS/cm, turbidity FNU, chlorophyll RFU, chlorophyll in µg/L, BGA-PC RFU, BGA-PC in µg/L, fDOM RFU and fDOM QSU. To avoid possible errors that might have occurred while interpreting the results, the mentioned parameters were not analyzed.   Table 3 shows the grouped data for 14 parameters that were collected during a one-year period (May 2018-May 2019). In the analysis, only four parameters were used (percentage of dissolved oxygen at saturation, ODO%, in mg/L; temp in • C; pH; and oxidation-reduction potential, ORP, in mV). The reason was tufa formation around the probe sensors, which resulted in a large variety of data for the following parameters-SpCond in µS/cm, turbidity FNU, chlorophyll RFU, chlorophyll in µg/L, BGA-PC RFU, BGA-PC in µg/L, fDOM RFU and fDOM QSU. To avoid possible errors that might have occurred while interpreting the results, the mentioned parameters were not analyzed. In the observed period, the ODO% averaged 84.04%. The maximum value was recorded in May 2018 (94.38%) and the minimum in October (72.80%). Generally, a percentage of 80% to 120% is considered very good, while higher (>125) and lower (<50) values can endanger the water biome [59]. The amount of dissolved oxygen (ODO, in mg/L) in a one-year period averaged 8.34 mg/L. The highest level was recorded in February 2019 (10.53 mg/L) and the lowest in August and September of the same year (6.49 mg/L). Compared to other streams in karst areas [60,61], slightly lower values of dissolved oxygen concentration and saturation were recorded. These lower values can be associated with higher organic matter content, which can cause eutrophic conditions. This is expected given the location of the probe. Specifically, the probe was placed under the path into the Mjernik stream, which extends from Marasovića Lake towards the test surface. This location was chosen for security and logistical reasons. However, this location is characterized by extremely rich and dense aquatic vegetation, which is why in the initial field surveys, it was recognized as a potential dried-up watercourse. Therefore, the lower values measured during August and September may be due to slow and calm flow [62], dense aquatic vegetation and fast metabolism of O 2 -consuming organisms. The concentration of dissolved oxygen is inversely related to water temperature, so a negative Pearson correlation coefficient (−0.75) between water temperatures and dissolved oxygen concentrations was expected.

Physico-Chemical Characteristics of Tufa Forming Watercourses at Test Surface
The average annual water temperature at the selected location was 16.32 • C. The maximum water temperature was measured in August 2018 (25.59 • C) and the minimum in January 2019 (7.06 • C). The average pH measured at the selected location was 8.65. The measured pH values were similar in other active tufa-forming systems [22,[62][63][64][65]. No significant seasonal fluctuations in pH were observed. Nevertheless, the analysis of pH values showed a clear upward trend (Table 3) as a result of the abundant higher aquatic plants that surround this location. Specifically, underwater photosynthesis during the day is expected to increase the pH because CO 2 is extracted from the water. Then, at night, photosynthesis eventually stops, the pH drops and respiring organisms add CO 2 to the water [66]. However, when plants or algae are growing rapidly, more CO 2 is removed each day by photosynthesis than is added each night by respiration. As a result, potential alkalization can occur in natural waters brought about by higher aquatic plants [67].
The average annual value of the ORP (mV) was 360.6. The highest value was recorded in January 2019 (450.9) and the lowest was recorded in May 2018 (206.7). Higher ORP values are an indicator of the presence of a greater amount of dissolved O 2 in water. The measured values meet the criterion that the ORP ranges from 300 to 500 mV in healthy waters [68].

Discussion
The uncontrolled growth of invasive vegetation and anthropogenic influences (bridge construction, gravel accumulation, etc.) have been recognized as a key problem causing negative hydrological changes and thus endangering the sustainable conditions of tufa formation in the wider SB area. The park management has recognized these threats and allowed the removal of invasive vegetation within the selected test surface and mitigation of the detected anthropogenic impacts. The removal of A. altissima and clusters of gravel underneath the path, which was raised, caused smaller micromorphological and larger hydrological changes on the test surface. Increased water flow caused the accumulated plant material to wash away, resulting in the reactivation of dried-up tufa-forming watercourses, lakes and smaller waterfalls. Although, in a one-year period, no regeneration of A. altissima on the test surface was observed, it is necessary to implement active and frequent control and mapping practices [69] as a basic component of the NPK management measures. However, it is necessary to emphasize that the removal of large amounts of macrophytes and invasive vegetation can lead to an increase in the erosive effect of the flow for a short period of time.
The reactivation of the dried tufa-forming watercourses is supported by the following evidence. First, the video was taken after the implementation of the abovementioned management measures. It shows water flowing through previously inactive tufa-forming watercourses. Furthermore, the reactivation is evidenced by the flow velocity data. Now, the flow velocity measured at the defined locations is in the range of values that are recognized as one of the basic conditions for tufa formation (0.5-3.5 m/s) [61,70]. Furthermore, reactivated tufa-forming watercourses do not have pronounced erosive properties, although minor erosive events have been measured [56], especially on tufa-forming organisms. Specifically, some authors state that a flow velocity faster than 1.5 m/s [71] or even 0.8 m/s can erode the tufa surface. However, flow velocities greater than 1.5 m/s were recorded only in April 2018 at the P1 location. Nevertheless, the values of the flow velocity at P1 (average value 0.897 m/s) and P7 (average value 0.808 m/s) indicate that the Vrtnjak tufa-forming watercourses at some locations had potential erosive properties. Although none of the sites of flow velocity measurement dried up, it is important to point out an example of how the non-removal of naturally fallen vegetation can influence the change in water runoff. Specifically, during the removal of the PLs in the colder part of the year, a change in the direction of water runoff within tufa-forming watercourses Crni grab was observed. This was the result of the fallen large tree and the accumulated plant material that created a natural dam (Figure 37). After the removal of the fallen tree, the watercourse was reactivated.
Water 2020, 12, x FOR PEER REVIEW 30 of 36 (0.5-3.5 m/s) [61,70]. Furthermore, reactivated tufa-forming watercourses do not have pronounced erosive properties, although minor erosive events have been measured [56], especially on tufaforming organisms. Specifically, some authors state that a flow velocity faster than 1.5 m/s [71] or even 0.8 m/s can erode the tufa surface. However, flow velocities greater than 1.5 m/s were recorded only in April 2018 at the P1 location. Nevertheless, the values of the flow velocity at P1 (average value 0.897 m/s) and P7 (average value 0.808 m/s) indicate that the Vrtnjak tufa-forming watercourses at some locations had potential erosive properties. Although none of the sites of flow velocity measurement dried up, it is important to point out an example of how the non-removal of naturally fallen vegetation can influence the change in water runoff. Specifically, during the removal of the PLs in the colder part of the year, a change in the direction of water runoff within tufa-forming watercourses Crni grab was observed. This was the result of the fallen large tree and the accumulated plant material that created a natural dam ( Figure 37). After the removal of the fallen tree, the watercourse was reactivated. Therefore, the downstream part of the stream dried out and the PLs at the location Near the memorial ran out of water. This resulted in a lower TGR (3.298 mm a −1 ) than in other locations (In the ivy and No name) in the same fluvial environment (MFW) and a higher TGR in the warm period (2.832 mm a −1 ) than in the cold period (0.209 mm a −1 ). This and other examples ( Figure 38) clearly demonstrate how non-maintenance of tufa-forming watercourses (removal of fallen and accumulated vegetation) can compromise sustainable sedimentation conditions at some micro-locations.
The reactivation of tufa-forming watercourses is also evidenced from the results of the TGRs. Specifically, on all PLs initially placed in inactive (dried-up) flows, the TGR was measured in the range of from 0.327 to 6.770 mm a −1 . The average TGR per PL in reactivated tufa-forming watercourses was 4.267 mm a −1 , which is 3.440 mm a −1 lower than the average TGR per PL (7.707 mm a −1 ) measured in previously active tufa-forming watercourses. Thus, in reactivated tufa-forming watercourses, a lower TGR (by a factor of 1.8) was measured compared to PLs that were placed in watercourses before the implementation of management measures. This was due to a combination of two factors. The first involves the specific micro-location conditions of the PLs at the Duck waterfall location. Specifically, PL07 lay in the immediate vicinity of the moss cluster, while the adjacent PL25 was farther away. The moss cluster overgrew the measurement surface of PL07, which gradually increased the nucleation surface and thus artefactually increased the TGR. This is expected. For example, Reference [72] stated that certain mosses can form spongy layers with a growth rate of up to 140 mm a −1 . If PL07 were excluded from the analyses, the average TGR on PLs placed in initially active tufa-forming watercourses would be 5.471 mm a −1 , which is only 1.28 times higher than the Therefore, the downstream part of the stream dried out and the PLs at the location Near the memorial ran out of water. This resulted in a lower TGR (3.298 mm a −1 ) than in other locations (In the ivy and No name) in the same fluvial environment (MFW) and a higher TGR in the warm period (2.832 mm a −1 ) than in the cold period (0.209 mm a −1 ). This and other examples ( Figure 38) clearly demonstrate how non-maintenance of tufa-forming watercourses (removal of fallen and accumulated vegetation) can compromise sustainable sedimentation conditions at some micro-locations.
The reactivation of tufa-forming watercourses is also evidenced from the results of the TGRs. Specifically, on all PLs initially placed in inactive (dried-up) flows, the TGR was measured in the range of from 0.327 to 6.770 mm a −1 . The average TGR per PL in reactivated tufa-forming watercourses was 4.267 mm a −1 , which is 3.440 mm a −1 lower than the average TGR per PL (7.707 mm a −1 ) measured in previously active tufa-forming watercourses. Thus, in reactivated tufa-forming watercourses, a lower TGR (by a factor of 1.8) was measured compared to PLs that were placed in watercourses before the implementation of management measures. This was due to a combination of two factors. The first involves the specific micro-location conditions of the PLs at the Duck waterfall location. Specifically, PL07 lay in the immediate vicinity of the moss cluster, while the adjacent PL25 was farther away. The moss cluster overgrew the measurement surface of PL07, which gradually increased the nucleation surface and thus artefactually increased the TGR. This is expected. For example, Reference [72] stated that certain mosses can form spongy layers with a growth rate of up to 140 mm a −1 . If PL07 were excluded from the analyses, the average TGR on PLs placed in initially active tufa-forming watercourses would be 5.471 mm a −1 , which is only 1.28 times higher than the TGR in reactivated watercourses. The second factor refers to lower periphyton abundance in reactivated tufa-forming watercourses compared to initially active watercourses [25]. In other words, microfauna have just begun to colonize reactivated tufa, forming watercourses. However, the periphyton community at these locations can have higher diversity, presumably due to lower competition and predatory rates [73]. In conclusion, the measured TGR, environmental data and protozoan assemblage patterns that were analyzed by Reference [73] show that reactivated tufa-forming watercourses are now representing favorable habitats for microfauna colonization.
Water 2020, 12, x FOR PEER REVIEW 31 of 36 TGR in reactivated watercourses. The second factor refers to lower periphyton abundance in reactivated tufa-forming watercourses compared to initially active watercourses [25]. In other words, microfauna have just begun to colonize reactivated tufa, forming watercourses. However, the periphyton community at these locations can have higher diversity, presumably due to lower competition and predatory rates [73]. In conclusion, the measured TGR, environmental data and protozoan assemblage patterns that were analyzed by Reference [73] show that reactivated tufaforming watercourses are now representing favorable habitats for microfauna colonization. In future research, it is desirable that a complete monitoring system of the tufa landscape includes additional components. These would include monitoring the dynamics of microfauna colonization, that is, analyzing the periphyton community, measuring the sun intensity of each site at which the dynamics of tufa formation are monitored and isotopically analyzing of newly formed tufa samples. It would be desirable to measure all parameters within each component at the same X, Y and Z locations. In other words, each site at which TFD is monitored should have specific data about flow velocity, physico-chemical parameters, sun intensity, periphyton, isotope data and so forth. This would allow the development of accurate prediction models of TFD and identify key criteria that influence this process.

Conclusions
In this research, a multiscale framework for the sustainable management of tufa-forming watercourses is proposed, with the main objective of achieving and maintaining sustainable conditions for tufa formation. This can provide the administration of protected tufa landscape guidelines for sustainable management of sensitive tufa area. The development of this methodological framework featured a range of activities at different levels of research. Since these experimental measures has not yet been applied in other protected areas; their effectiveness was tested within smaller test surfaces. Selective removal of invasive vegetation was performed on In future research, it is desirable that a complete monitoring system of the tufa landscape includes additional components. These would include monitoring the dynamics of microfauna colonization, that is, analyzing the periphyton community, measuring the sun intensity of each site at which the dynamics of tufa formation are monitored and isotopically analyzing of newly formed tufa samples. It would be desirable to measure all parameters within each component at the same X, Y and Z locations. In other words, each site at which TFD is monitored should have specific data about flow velocity, physico-chemical parameters, sun intensity, periphyton, isotope data and so forth. This would allow the development of accurate prediction models of TFD and identify key criteria that influence this process.

Conclusions
In this research, a multiscale framework for the sustainable management of tufa-forming watercourses is proposed, with the main objective of achieving and maintaining sustainable conditions for tufa formation. This can provide the administration of protected tufa landscape guidelines for sustainable management of sensitive tufa area. The development of this methodological framework featured a range of activities at different levels of research. Since these experimental measures has not yet been applied in other protected areas; their effectiveness was tested within smaller test surfaces. Selective removal of invasive vegetation was performed on selected test surfaces. A system for monitoring the quality of the tufa sedimentary environment with four components was established. Implementation of the proposed methodological framework on the selected test surface suppressed the proliferation of invasive vegetation and enabled reactivation of dried-up tufa-forming watercourses, thus achieving sustainable tufa-forming conditions. This was confirmed with measured TGRs, the water flow velocity at reactivated tufa-forming watercourses, physico-chemical parameters, the occurrence of tufa-forming organisms at PLs and the absence of invasive vegetation regeneration. Potential beneficiaries of the proposed methodological framework for this dissertation include national, regional and local authorities, as well as governing boards of protected and endangered natural areas. The results obtained in this study may be useful to scholars involved in studying the dynamics of tufa formation through various aspects. The derived data serve as a reference state when conducting spatiotemporal analyses in TFD research.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: