Using historical data to access the surface subsidence in the vegetable belt of the Three Lakes Region, Switzerland

The ‘Grosses Moos’ in the Three Lakes Region (Berner Seeland) was formed by the Wallis glacier resulting in a peat bog-dominated landscape. During the last two centuries this area was drained by a complex channel system to enable highly fertile agricultural activity. As such, the region is the vegetable belt of Switzerland. However, as the peat has degraded due to overexploitation, the situation for agriculture production has become critical. Consequently, measures are needed to prevent soil degradation and water accumulation. The extent of surface changes for the last ~ 100 years is, however, not precisely known and is assumed to have varied spatially. To understand the historical evolution of this area, we used a unique map from 1920 to compare the surface height with the new-est digital terrain models (DTM) in order to estimate surface subsidence, and thus soil degradation, for the last 100 years. More than 44,000 single measurement points drawn on the historical map were digitized to derive the DTM that served as a basis for further analysis. The in-depth investigation of the observation methods and the accuracy assessment allows us to conclude that differences in the DTMs of up to 2.4 m (i.e., averaged rate of 2.4 cm yr −1 ) can be attributed to large-scale soil subsidence and that degradation was heterogeneously distributed over the investigated area. The data provide the basis for further soil restoration efforts in the area.


Introduction
Areas around rivers and lakes were historically difficult to populate, as well as cultivate, due to frequent flooding and uncontrolled river systems.Moreover, they were often dominated by peat bogs and swamps, which are also a source of potential diseases.During industrialization, empowered by technological developments, such lakedominated landscapes became the focus of new drainage systems and water correction efforts (Nast, 2006).One of these areas is the "Grosses Moos" in the western part of Switzerland between the lakes of Murten, Biel, and Neuchâtel.This "Three Lakes Region" was made accessible to agricultural use in the late nineteenth century.Drained through a complex channel system and complemented by the optimization of the Swiss Jura water levels, it became the primary vegetable-producing area in Switzerland and is commonly called "the vegetable belt" (Nast, 2006).
Peatlands, which are the characteristics of this area, are highly fertile for agriculture but are also highly sensitive to soil degradation and peat decomposition, especially when the water table decreases (Egli et al., 2021;Hutchinson, 1980).The degradation of peatlands, resulting from the loss of organic matter (i.e., biochemical degradation), shrinkage, and consolidation when drained, as well as the intensive use of mires for agriculture, contributes significantly to soil degradation.The estimated average C-losses in drained agricultural land often totals several tons per hectare (e.g., 4.9 t ha −1 in Egli et al., 2021) and can cause significant soil subsidence that is important to monitor.The calculation of such soil subsidence when it occurs to a few meters maximum is usually done using single profiles, monitoring plots (e.g., Evans et al., 2022), or transects and either by observing them over a given time or by retrospectively deriving the occurred soil losses (van Asselen et al., 2018).Using this approach, van Asselen et al. (2018) were even able to derive time-dependent subsidence rates in the Rhine-Meuse delta over the last 1000 years.The indirect and qualitative derivation of ground subsidence caused by more intrusive activities (e.g., mining), has also been studied using remote sensing techniques (Landsat monitoring; Yi et al., 2021), and Ikkala et al., (2021) quantified subsidence and corresponding rates using historical maps and LiDAR (Light Detection And Ranging) measurements.Furthermore, Hoogland et al. (2012) modeled degradation rates based on 1923 historical point locations (for an area of about 12 km 2 ) combined with recent LiDAR measurements of a polder in the Netherlands.Soil subsidence due to peat degradation is also relevant and observed in subtropical regions and has been modeled using historical input data (water table, biomass input, peat thickness, etc.) from over ~ 100 years  by Rodriguez et al. (2020) in the Florida Everglades.
In the case of the "Grosses Moos", agricultural overexploitation gave rise to a very thin soil layer in the lower decimeter range (Zihlmann et al., 2019) that endangers future agricultural use.To allow continued sustainable use for agriculture, as well as to support the restoration of critical areas, data about the initial state, historical degradation, and related subsidence rates are needed.With this paper we present an analysis of soil subsidence since the early twentieth century.A historical map from 1920 was available for analysis and was compared with modern surface terrain models, the digital terrain models (DTM), based on different observation techniques.The paper includes an in-depth analysis of expected measurement accuracy and the various error sources alongside a quality assessment of the different data sets.We interpret the differences between the historical and recent DTMs to provide insights into the heterogeneously distributed rates of degradation over a time span of approximately 100 years in the Grosses Moos area.

History and overview of the study area
The study focuses on the Grosses Moos (the great marsh), which is the most important agricultural area for Swiss vegetables (Nast, 2006).The region lies in the Western part of Switzerland as part of the Swiss plateau and borders the Jura mountains to the northwest.The area is surrounded by the three lakes of Murten, Biel, and Neuchâtel and lies between the cities of Ins, Kerzers, and Aarberg (Fig. 1).The area of interest is very flat, having only 15-20 m difference in altitude over 15 km of horizontal distance, and it is cultivated with different crop types and vegetables.
The topographic characteristics were shaped by the Wallis glacier system and experienced regular coverage by glaciers during the Pleistocene (e.g., Egli et al., 2021;Wohlfarth-Meyer, 1987a).Owing to several ice ages and advances of the glaciers, the bedrock was eroded and filled again with sediments before finally becoming covered by mires.The three lakes have a predominantly glacial origin (via moraines), and the landscape of Grosses Moos was further shaped by the Aare River until about 5000 years ago (Wohlfarth-Meyer, 1987a).
Between 1868 and 1897, the first infrastructural intervention to regulate the level of the three lakes and the Aare River was performed resulting in better living conditions for the local population as productive agriculture was made possible (Vischer & Feldmann, 2005).Nevertheless, the area continued to experience flooding and strong surface subsidence that impacted agricultural activities.Thus, a second optimization was necessary between 1962 and 1973 (Nast, 2006), which helped lower the flood risk.Nevertheless, the degradation of the mires, and in particular the compression and the mineralization of organic matter, has resulted in large-scale subsidence that puts the Swiss vegetable belt under severe risk.

Historical topographical maps
This analysis is based on historical maps found in an archive that was made accessible for digitization and investigation.The original maps from 1920 were available on paper and drawn using the plane table technology (Wainwright, 1915).With this technique, the surrounding surface of the surveyor was observed by measuring angles using a ruler with a telescopic sight mounted on a horizontal table for drawing the map.The height difference from the surveyor's position to all visible points around the table was measured in relation to a known reference point, resulting in a very detailed map (Wainwright, 1915).The topography of the Grosses Moos was mapped using a grid of around 20-30 m between the measurement points spacing with additional measurement where needed (e.g., along rivers or roads).Afterwards, these measurements were interpolated with contour lines and complemented with topographical information.Terrain information (i.e., contour lines and escarpment) is drawn in brown, while water bodies are represented in blue and paths in black.Figure 2 shows a section of the available map.Each of the observation points is marked with a dot and a corresponding number between 0.0 m and 9.9 m.Two digits are included per dot to reduce the writing effort and for a better arrangement.The number is calculated as the positive difference to a local height in decimeters.The heights of the region shown in Fig. 2 are thus relative to 430.0 m rounded down to the decimeter from the reference point (i.e., the triangle marked with "431.62"m).The triangle indicates the connection to the higher-order Swiss observation network.Due to the given digits of the reference height and the individual measurement points in combination with the very flat terrain, we can expect a measurement accuracy in the lowest decimeter range.
An overview of available maps for the region Berner Seeland is shown in Fig. 3 covering the agricultural area from the Lake of Neuchâtel in direction to Ins until Kerzers and Finsterhennen.The background map used in Fig. 3 is based on a swisstopo (Federal Office of Topography, Switzerland) topographic map from 2019 and does not show the historical landscape from 1920.The map areas that are considered within the presented work are labeled with 4, 5, 6, 7, 8, 9, 11 and were digitized for the analysis of the area of the Grosses Moos.
The historical maps were scanned with a dedicated map scanner and imported into Quantum GIS (QGIS) for georeferencing (plugin Georeferencer GDAL, QGIS 3.2.2).Georeferencing was performed using the coordinate system frame at the edge of the maps and transformed using a Helmert transformation accounting for translation, rotation, and scale into the Swiss coordinate system.Thanks to the high-quality scans, a single scale was sufficient to account for the distortion.The coordinate grid used for the historical maps is based on the land survey of 1903 (LV03); however, it still uses non-military coordinates.The maps were, thus, labeled with the coordinates relative to the coordinate system origin in Berne (0N/0E).After 1920, all coordinates changed to the military system adding an offset of 200,000 m and 600,000 m for north and east coordinates, respectively.By adding these offsets, interchange latitude and longitude can be avoided (Swisstopo, 2016).
After georeferencing, the maps were digitized by point registration and by manually entering the related height (Fig. 2).A total of 44,319 points were digitized covering an area of ~ 56 km 2 .The number of points distributed within the different tiles can be found in Table 1.The table also shows that maps were generated by different surveyors in the same year.Looking at the total number of more than 44,000 recorded points, creating these highresolution maps with such high quality, must have been an enormous work.We highly acknowledge the work of the surveyors as well as the persons digitizing the points that made this analysis possible.
For further processing, the points were merged into one DTM modeled with the Triangulated Irregular Network (TIN) interpolation.TIN was used because the point network is irregular but dense (about every 20 m a point) and almost equally spaced (e.g., Mitas & Mitasova, 1999).The DTM is solely derived from single point height observations as the object information such as the rivers and streets were not digitized as break-lines.Blunder detection was performed based on the visual inspection of the DTM.Thanks to the very flat investigation area, blunders can easily be detected and corrected.

Recent data
For comparison with the present-day situation, surface information available through the Swiss Federal Office of Topography (swisstopo) was used.The vector data set swissTLM3D (Swisstopo, 2022b) includes surface information as individual geo-objects such as water bodies, streets, and buildings.This data was used for comparing the horizontal alignment of specific objects from historical and new maps.The results are described in Fig. 4 and the section Validation of height system.
The surface change analysis is based on the newest DTM swissALTI3D available from swisstopo from 2016 (southern part) and 2018 (northern part).The swiss\ ALTI3D is Switzerland's high-resolution digital terrain model where vegetation and development (such as houses) were subtracted.The swissALTI3D DTMs used for the area of interest have a resolution of 2 m pixel spacing with an accuracy of 0.5 m and are based on the LiDAR technology, a technology measuring the distance between the sensor (e.g., an aircraft) and the surface by sending out laser beams (Swisstopo, 2022a).Thus, the resolution of the swissALTI3D is around 5× higher than the 1920 maps and includes also breaklines (i.e., rapid changes in the topography).The contours and/or break-lines drawn on the 1920 map were not digitized within this project.
In addition, the so-called "DTM-AV" of the year 2000 with data also measured by a LiDAR and interpolated to a 2 m pixel grid (derived from a point cloud having a density of 1 point per 2 m 2 ) was used as the second most recently available data set.The DTM is designated with an accuracy in height of ± 0.5 m for nonvegetated areas and ± 1.5 m for vegetated areas.As our area of interest is dominated by agricultural fields (i.e., changing vegetated area), the additional influence of the vegetated and potentially changing surface needs to be considered in the analysis.The vegetated areas could be visible in a noisy surface model and/or visible individual fields in the difference map.
The difference in height and quality between the two recent data sets lies mainly in the lower pixel density of the older data set resulting in a less detailed surface model.
All the data sets presented are available in the Swiss reference coordinate system LV03 (Swisstopo, 2016) and the Swiss reference height system LN02 (Hilfiker, 1902;Marti & Schlatter, 2002).Using the Swiss reference system in horizontal location from 1903 (LV03) and in height from 1902 (LN02) allows for a direct comparison of the different data sets without transformation.

Historical DTM: co-registration and interpolation
Several factors contributed to the overall residual error and the accuracy of the historical DTM.The sources of total uncertainty in georeferencing include the scanning process (e.g., scanner distortion), the state of the individual sheets (e.g., with crinkles), and the quality of the drawn map with the line width.The line width also contributes to precision limitations, along with the operator's choice of points and the digitalization process itself.Fig. 4 Comparison between the historical underlaid map (map 9) and the vector data set from swisstopo swissTLM3D (2022).Present day streams (dark red line) and paths (dashed purple lines) mostly still exist and are well co-registered, some buildings also still exist such as the Ziegelei in the Eastern part with the new buildings shown in orange and the historical ones in black.Surface Subsidence in the Three Lakes Region, Switzerland The results of the georeferencing procedure additionally depend on the original size and format of the map influencing the transformation residuals of the reference points.They differ after the transformation between 7 pixels and 17.5 pixels, with an averaged root mean square error (RMSE) of 11.43 pixels (Table 2).The point precision derived for the co-registration is in the expected range because the drawn line of the coordinate cross is around 2.5 pixels in thickness corresponding to ~ 1 m that is then resulting in the RMSE of 4.6 m for the total co-registration.The highest uncertainties are measured for Map 11 and Map 4. Map 4 is a particularly large sheet, strongly elongated in the East/West direction where more distortion and uncertainty are expected due to its shape, and Map 11 displays a rather small section located in the North of the area of interest.The effect of a horizontal shift of 7 m (in case of Map 4) would result in 0.009 m error on the height measurements if we calculate with a height difference of 20 m for 15 km horizontal distance.The effect of a horizontal misalignment in the range of the identified residuals can thus be neglected for the large-scale analysis.
Visual inspection between the georeferenced maps and the swissTLM3D vector data set (Swisstopo, 2022b) with a focus on streets, water bodies, and buildings revealed a very good alignment of the historical and recent data sets with differences in the lower meter range (Fig. 4).Many paths still exist at the same location as well as the major streams despite the 100 years between the two data sources.However, some paths might also have been moved or were not accurately drawn on the map.The horizontal alignment of the historical and new data can be assumed to be good enough so that a comparison of the DTMs can be performed without further validation of the horizontal alignment (also thanks to the very flat terrain).

Validation of height system
To check for unexpected biases between the historical and recent data sets, measurement points from 1920 were compared with the corresponding height extracted from the DTM 2016/2018 at the same coordinate.Areas expected to have remained stable over the course of 100 years were sought after.Finding such areas was challenging, especially since most of the area comprises agricultural fields, with only a few urban areas included in the mapping.In some areas (e.g., north and east of Müntschemier), roads and channels were newly structured, while as shown in Fig. 4 other areas largely retained the 1920 structure.The comparison of the points presented in Table 3 is used to estimate if large biases (i.e., in the decimeter range) might exist across the entire area, which could indicate systematic inconsistencies (for example, in the reference height systems).Differences between these reference points also give an indication in which range the absolute accuracy lies.Table 3 lists the differences in height for selected points.The distribution of these validation points including the two examples (Kerzers and Brühlteilen/Finsterhennen) are shown in Fig. 5.Additional validation points listed in Table 3 are shown in Appendix A. The height of these points as written on the 1920 map was directly compared to the height extracted from the 2016/2018 DTM at the same location.

Comparison of historical and recent DTM
On the 1920-DTM, the surface structure is visible (Fig. 6 left) with the agricultural fields surrounded by water channels or pathways and some river meanders.When comparing it with the recent 2016/2018-DTM (Fig. 6 right), the different resolutions of the data sets are easily visible with many more details such as creeks and paths in the recent model.Well discernible in both DTMs are the Broye channel in the southern-western part (connecting the two lakes) and the old Aare meanders from around 5000 years ago or older (Wohlfarth-Meyer, 1987b).As already shown in Fig. 4, the historical path and creeks still exist nowadays; however, most of them are missing in the historical DTM as only the height measurements were digitized, and break-lines such as the small creeks and paths smooth out when generating the DTM using the TIN interpolation method.The individual topographic objects such as streets are, as a result, missing.We were mainly interested in the large-scale surface changes, and individual topographic objects were omitted in the analysis.However, the differences in resolution need to be considered for the interpretation of the height values (Fig. 7).
The difference between the DTM from 1920 and 2016/2018 is shown in Fig. 7, where red represents land subsidence, while green or blue colors mark areas with higher surface height in 2016/2018 than in 1920.The clearly recognizable river meanders correspond to the old Aare riverbed.Parts of these old meanders are filled up with standing water or ponds and are thus particularly well recognizable in the recent DTM.The highest subsidence was measured in the Kerzersmoos/Längenmatten area with up to 2.4 m in height loss since 1920.Another area with large subsidence can be seen in the region of Birkenhof and Lindenhof (south and south-east of Ins) with height differences of up to 1.6-1.8m, and in the region of Finsterhennen with similar values.These areas can be considered as hot-spots.Close to the old Aare riverbed, however, lower or almost no subsidence was detectable (white area in Fig. 7).Areas with positive surface changes can mostly be attributed to new urban areas (e.g., new buildings), forest areas (Seewald/Tannenhof ), and lake areas with different water levels.The largest region of positive surface changes was in the range of 0.20-0.40m near the Broye channel in the southwestern part of the map (dark red line in Fig. 7).
When looking at the statistics using the histogram in Fig. 8B, we observe a mean and median subsidence over the entire region of − 0.70 m with a standard deviation of 0.67 m.The area of investigation underwent irregular subsidence detectable with two different local maxima at − 0.31 m and − 1.11 m, with a local minimum at − 0.84 m.Based on these two peaks, the investigated area can be separated into the area surrounding the old Aare River meander with lower subsidence rates (< 0.84 m) than the red areas with higher rates (> 0.84 m) using the local minima between the two peaks for distinguishing the two magnitudes of subsidence (Fig. 8B).
The area of the Grosses Moos has been used for cultivating vegetables since the nineteenth century (Nast, 2006), so land degradation is not new.Agricultural use intensified in the twentieth century because a higher yield was needed, and the soils were accordingly fertilized and plowed.Consequently, the recent evolution of surface subsidence is of interest.For this reason, subsidence rates were also checked using the DTM from swisstopo of the year 2000.Figure 9A depicts this comparison between the DTM-AV from 2000 and the historical DTM identifying the same area of subsidence.Comparing the two recent DTMs from 2016/2018 and 2000, as represented with the histogram in Fig. 9B, reveals a normally distributed curve with a small bias of − 0.04 m.Fitting a standard normal distribution ('fitdist' in Matlab, significant > 95% from One-sample Kolmogorov-Smirnov test) to the data results in µ = − 0.06 m and σ = 0.18 m.The two spikes at around + 0.25 cm are attributed to the changes in the lake water level of the two lakes Biel and Neuchâtel.Due to the small differences in height of the two recent DTMs in comparison with the stated accuracies of both DTMs, a detailed comparison is omitted.

Significance of changes
The original map from 1920 exhibits a high grade of topographical detail, and the point grid of the individual measurements follows a 20-30 m grid.The point grid is sufficiently dense to capture small-scale variations in surface height, and additional points were recorded as needed to capture specific surface features.This allows for a comparison of the 1920 surface height with a recent DTM from 2016/2018, in particular regarding the large-scale changes of the very flat agricultural areas.Within the last 100 years, the topographic characteristics across large areas have strongly changed.This can be well seen, for instance, with changes in the newly urbanized areas (Fig. 5B).Nevertheless, the historical structure of drainage channels and pathways was largely retained in large parts of the agricultural fields (Fig. 4), except in the northeast of the mapped region (Fig. 5C) or in the southern part near Galmiz.The comparison of the topographical structure also confirmed the reliable georeferencing of the historical map with a horizontal accuracy of ~ 5 m, allowing for a direct comparison of the height information by subtracting one DTM from the other.The differences in height between the historical and most recent DTMs are distributed between 0 m and − 2 m (histogram in Fig. 8A).A detailed error analysis was thus needed to confirm that the changes in surface height were larger than the error ranges of the individual methods.
The 2016/2018 DTM had a designated accuracy of 0.5 m in the product specification (Swisstopo, 2022a).This range of accuracy can be confirmed when comparing the height of the 10 selected reference points from the historical and recent DTM that showed a bias of − 0.4 m and a standard deviation of 0.41 m.This analysis is based on the assumption that main streets and roads stayed stable over this time period (except for instance changes in road cover).An accuracy of 0.41 m, with a slight overall tendency of 0.4 m subsidence even for roads, appears reasonable for this region.Considering the long timescale, it is plausible that roads could have undergone slight subsidence as well.This also confirms the very high-quality maps from 1920 as well as the DTMs from 2016/2018.Even though the most recent DTM data set originates two years apart (the northern part in 2016, the southern part in 2018), no systematic effect can be seen in the comparison of the merged DTM with the historical data (Fig. 8A).The subsidence rate is certainly too low to be detected within the 2-year difference in acquisition time, and swisstopo assured the temporal consistency.
The systematic double peak shown in the histogram in Fig. 8B can be clearly attributed to different regions of subsidence (see Fig. 8A) and is independent of the acquisition date.The subsidence areas, with a maximum of -0.31 m in the areas around the old Aare riverbed (blue in Fig. 8A), may not be significant regarding the accuracy of the specified accuracy of 0.5 m for the 2016/2018 DTM.However, the red areas with a subsidence between 0.84 m and 2.4 m are certainly above the measurement accuracy, and these areas can be identified as the most critical areas.
The additional comparison of the 2016/2018 DTM with the DTM from 2000 assists with the interpretation of this double peak and gives indications of short-and long-term effects.Due to the specified accuracies of the two data sets, the subsidence of the last 20 years was too small to be detected with the presented methods (i.e., less than a few decimeters).The large-scale subsidence in the range of meters was, however, not an effect of the last 20 years only, also discernible with the lacking double peak when comparing Fig. 8B (100 years difference) with Fig. 9B (less than 20 years difference).

Subsidence pattern
With our analysis, we can confirm a strong peat degradation as well as related subsidence of large parts of the investigation area within the last 100 years exhibiting a strongly heterogenous spatial pattern.Owing to the different soil properties, restoration efforts, surface cover, or agricultural exploitation history, different subsidence rates are expected and confirmed also in other studies (e.g., Egli et al., 2021;Leifeld et al., 2011).When looking at the height-difference map in Fig. 7, three main types of changes emerge.There is a large-scale pattern within the area-wide subsidence (red colors) with some increase in surface height (green/blue colors), and small-scale features related to changes in paths and creeks.When doing the interpretation, particularly for these small-scale features, one needs to consider that the individual measurement points from the 1920-DTM are interpolated with a TIN model without using break-lines.Thus, the paths and streets that are well discernible are partly attributed to the digitalization procedure.A more one-to-one comparison of the finest structures would be needed to see if they can be deduced to stable surfaces, differences in the surface material, changes in the course of the paths, or sampling resolution of the data sets.Additionally, areas remaining stable or experiencing an increase in surface height can be attributed to additional housing, unknown restoration efforts (e.g., by locally adding soil material), different soil properties (e.g., Dawson et al., 2010), different land use and management (e.g., different type of vegetables or crop) or other unknown effects.For example, the area around the Broye channel might have been restored or stabilized and the village of Sugiez at the end of the same channel is well detectable due to many new houses.
In the study area of the Three Lakes Region, we observed a subsidence rate of up to 2.4 cm yr −1 estimated over a duration of 100 years.The subsidence map indicates that the area surrounding the old Aare meanders underwent a lower subsidence rate, (i.e., < 0.84 cm yr −1 ) even though the old riverbed subsided more strongly (red in Fig. 8B).Areas within this old riverbed are still (or once again) covered by standing water (i.e., water ponds) and restored biotopes (e.g., regions of Bellechasse, Agriswilmoos, or Ziegelei).Thus, some of the changes in surface height can be attributed to restoration and changing water surfaces (e.g., Leiser et al., 2009).In general, large-scale heterogenous subsidence of the agricultural fields in the regions of Kerzersmoos (Kerzers, FR), Burgmoos (Galmiz, FR), or Birkenhof (Ins, BE) of 0.84 m to 2.4 m (red areas in Fig. 8B) was confirmed within our study.Our subsidence rates, as well as the strong heterogeneity in the spatial pattern, agree with other studies in the same region using different techniques such as Leifeld et al. (2011) who performed a study in the Lindenhof region by measuring surface level changes and soil profiles.Furthermore, the heterogenous distribution of surface lowering was also confirmed by Egli et al. (2021), who estimated in the same study region a surface lowering in the range of 76 cm to 249 cm within about 150 years.
A similar magnitude of degradation was found in other study areas in Europe with annual rates between 0.75-1.64cm yr −1 (Dawson et al., 2010;Eggelsmann, 1978;van Asselen et al., 2018).Dawson et al. (2010) suggested a decreasing rate over time that depends on the peat type (semi-fibrous and fibrous peats versus humified peat) with the driving factor for soil degradation being shrinkage or being controlled by biochemical processes.Van Asselen et al. (2018) also investigated the surface changes over time and found present-day values of up to 1.4 cm yr −1 .On a timescale of 1000 years (total rate of ~ 4 m) the current rates are strongly increased.In general, the rate of degradation thus depends on local conditions and changes as well as the timescale considered.Detailed knowledge of the local situation is thus needed to interpret measured subsidence rates and to predict a future change in surface height and degradation of the soils in a particular region.Within our study area, we can confirm that the rate of subsidence was below a few decimeters in the last two decades with a total change of up to 2.4 m within 100 years.Changes over the last 20 years fall below the error range of the techniques used, with a difference of only a few decimeters when comparing the two DTMs from 2016/2018 and 2000.Hence, we cannot confirm an increased rate of degradation over the last 20 years.
All the previously mentioned studies share the common observation that the subsidence rate is spatially heterogeneously distributed.However, the maximal subsidence rate is highest in our study.Moreover, our study is unique in terms of spatial and temporal resolution (i.e., 20-30 m spacing), enabling the detection of areas of degradation that might have been missed in other studies.Compared to other publications, our investigation presents the highest spatial resolution and precision over time (including clearly specified error ranges).

Conclusion
The analysis demonstrates high subsidence rates within the last 100 years in the Grosses Moos region in Switzerland with a high and unique spatial resolution and coverage.Due to the large-scale trends in addition to the measurement and processing methods, we can exclude the possibility that the height changes were caused by (a) systematic effects due to the evolution of the coordinate system, (b) the digitalization of the maps and georeferencing process, (c) different cultivation states (e.g., freshly plowed soil against high crop plants) and/or (d) measurement uncertainties.We can also confirm that subsidence rates are heterogeneously distributed.This data can thus serve as a basis for the future exploitation of the vegetation belt and restoration attempts.Identifying the regions having the highest subsidence rate helps to focus on certain regions where investigations and potential restoration efforts are needed so that the vegetable belt in Switzerland can be maintained in the future.

Fig. 1
Fig. 1 Overview of the study area Grosses Moos in the western part of Switzerland near Ins (inlet with location within Switzerland).The military coordinate system of Switzerland is shown with a north-oriented grid with 1 km spacing.The study area, surrounded by a red line, is located between the three lakes of Murten (south of the study area), Neuchâtel (west of the study area), and Biel (north of the study area) and with Aarberg located in the northeast corner of the map.(background maps: swisstopo)

Fig. 2
Fig. 2 Section of the historical map with each measurement point labeled with the relative difference to 430.0, measured in relation to the reference point labeled with a triangle and 431.62 m from a higher-ranking Swiss height network.Contour lines and other topographic details such as water channels are drawn in color.

Fig. 3
Fig. 3 Overview of available area covered by the historical maps.Tiles 4, 6, 7, 8, 9, and 11 were digitized and used for investigation of the area of Grosses Moos.(background map: swisstopo)

Fig. 5 AFig. 7
Fig. 5 A Overview of the reference points (red points; see also Table 3) with the historical map in the background and outline of the investigation area.B1 Validation point (red, Nr.5) located in Kerzers with the historical map in the background.B2 Orthoimage with validation point Kerzers (2018).C1 Validation point (red, Nr.8) Brühlteilen with historical map in the background.C2 Aerial image with validation point Brühlteilen (2016).(aerial image in B2, C2: swisstopo)

Fig. 8 AFig. 9 A
Fig. 8 A Difference between the two DTMs with color scale adopted to the two peaks in the histogram with red colors < − 0.84 and blue colors < − 0.84 with -0.84 as local minimum separating the two peaks visible on the histogram (B) of pixel values height difference between 2016/2018 and 1920 DTMs of the entire region of interest.Using this separation, the main area of soil degradation is visible with the red areas.(aerial image: swisstopo)

Table 1
Area covered by the different maps and number of digitized points.The last column includes the responsible person for drawing the individual map

Table 2
Residuals per map from the georeferencing procedure in QGIS using the plugin Georeferencer GDAL

Table 3
Result of pointwise comparison of the different DTMs with 10 points distributed over the area of investigation that are expected to have stayed stable over 100 years.Points are mainly selected on a street located near still-existing farms or buildings