Soluble phosphorus content of Lake Balaton sediments

ABSTRACT Lake Balaton has undergone rapid eutrophication in the last decades causing dramatic deterioration in water quality. Although water quality has been successfully improved, significant algal blooms have been experienced over the past years again. Since the high level of phosphorus is primarily responsible for algal blooms, it is necessary to explore the status of the nutrient content of the lake sediment. In this study, a 100 m resolution map was created by the reanalysis of an archive database (with measured ammonium-lactate extracted phosphorus concentration) using the regression kriging method with topography, hydrodynamics and vegetation auxiliary variables. The regression model with 0.23 R2 does not indicate a strong correlation with the auxiliary variables. However, based on the result of the validation, the accuracy of the final map after kriging was high (R2=0.81). The reanalysis of soluble phosphorus should represent a significant input to Lake Balaton's water quality management planning.


Introduction
Lake Balaton is the largest freshwater lake in Central Europe. Its unique natural values make it one of the most significant national treasures of Hungary, and after the capital city, it is the second most visited tourist destination in the country.
Eutrophication accelerated by human activities causes a global problem in the water quality of natural waters. From the nineteenth century, significant changes started in the area of Lake Balaton (railway construction, urban sprawl, excessive agricultural use of the surrounding fields), which led to a decrease in the water-covered areas of Kis-Balaton wetland. Kis-Balaton had played a key role in the natural filtering of the inflow waters. As a result of the drastic changes, Lake Balaton underwent a rapid eutrophication episode during the 1970s-1980s causing dramatic deterioration in water quality (Herodek, 1986). In the 1990s, eutrophication was halted by targeted environmental actions, which led to a noticeable improvement in the environmental condition of the lake (Honti et al., 2020). By the end of the 2000s, a water quality protection system was built in the area of Kis-Balaton to restore its natural filtering function (Istvánovics et al., 2007).
However, in the last two years, significant algal blooms have been experienced again, therefore the need for a detailed analysis of the environmental condition of Lake Balaton and the factors influencing the eutrophication has arisen.
During the 1970s and 80s, detailed research (taking more than 5.000 samples) was carried out on the physicochemical properties of Lake Balaton sediments including phosphorus concentration (Máté, 1987). In the early 2000s, digital sediment maps were generated from the point database using distance-based buffering or the Thyssen polygon formation method by the spatial extension of the measured chemical parameters (Csermák & Máté, 2004). The phosphorus map of this former map series can be seen in Figures 1-2. However, as it is clearly shown in the case of phosphorus map, formerly applied interpolation techniques did not take account of the spatial correlation existing between the data, which can be considered as a serious shortcoming of these techniques. It is also worth mentioning that detailed environmental  auxiliary variables, e.g. digital topography model for the Balaton bed as well as hydrodynamics data had not been available when this map series was edited.
This point database indeed reflects the state of the lake in the late 1970s and mid-1980s but it is the only available dataset that contains chemical parameter information on Lake Balaton sediments. The authors (Csermák & Máté, 2004;Máté, 1987;Sisák & Máté, 1993) who analysed the dataset earlier, could not use up-to-date geostatistical methods and GIS programs to create a more detailed synthesis and evaluation of the results.
GIS-based mapping methods have undergone tremendous development in recent decades. The most widely used and, at the same time, the most successful new method is regression kriging (RK). As auxiliary variables, which might have significant estimation power, are now available, it is worth re-evaluating the archive data with a more effective GIS method, and including these variables as predictors in order to create a more detailed sediment map about Lake Balaton.
The recurring serious algal blooms in the last two summers have drawn attention to the importance of further investigation of the nutrient dynamics of Lake Balaton. Therefore, in our study, we aimed to re-analyse this dataset with up-to-date methods of GIS and thus provide detailed spatial information on the phosphorus content of the bed sediment in the lake, which explicitly accounts not just for the environmental correlation existing between the auxiliary variables and the phosphorus content but for the spatial auto-correlation between the phosphorus data as well.

Study site
Lake Balaton is situated in the western part of the country. The surface area is 593 km 2 . Its elongated shape is 77 km long and about 7 km wide. Despite its large surface area, Lake Balaton is a shallow lake with an average depth of 3.2 metres. In the summer, the mean temperature of the water is 23°C. Out of ions, Ca 2+, Mg 2+ and HCO 3+ are present in the largest proportion in the water (Cholnoky, 1936;Herodek, 1992).
Zala River, which provides the largest inflow of water, discharges into the westernmost basin flowing through a wetland area called Kis-Balaton. The only outlet of the lake is the Sió Canal on the southern shore of Balaton, where the excess water is discharged through the sluice in Siófok to the Sió River (a canalized waterway connected to the Danube) in a controlled manner. The water volume of the lake is around 1.8 billion m 3 and the water spends about 4.7 years in the lake (lake retention time).
The trophic status and the degree of eutrophication give a good indication of the environmental condition of stagnant waters Boström & Petterson, 1982;Hutchinson, 1957;Wetzel, 1975;Williams, Syers, Harris, et al., 1971;. As a limiting factor, phosphorus determines the suspended and the settled biomass content, the amount of a-chlorophyll and algae, and the spreading of aquatic plants in the lakes. Organic carbon content also correlates with the phosphorus level. From the early 1970s, several studies focused on the investigation of organic carbon content in the sediments of Lake Balaton (Frankó & Ponyi, 1973;Herodek & Oláh, 1973;Herodek & Tamás, 1973;Kárpáti, 1970;1972;Pónyi et al., 1972). Oláh et al. (1977) later established that the high level of P was mainly responsible for the harmful algae bloom in the lake. It was also found that, under certain conditions, nitrogen and light can also be limiting factors (Herodek et al., 1988), but phosphorus is responsible for the negative change of the water quality (Auerswald, 1989;Jolánkai, 1986). In the history of Lake Balaton, the most serious state developed in the Keszthely basin, the smallest basin of the lake. The water of the Keszthely basin was classified as eutrophic in the first half of the 1970s and hypertrophic in 1982 (Somlyódy, 1983).
The nutrient input transported by inflow waters plays a key role in the eutrophication of Lake Balaton (Sisák & Máté, 1993). Joó and Lotz (1980) established that about two-thirds of the nutrients get into the lake due to floods in the catchment area. The authors considered the source of the nutrient to be of agricultural origin, however, they did not analyse what proportion derived from the natural nutrient content of the soils and what proportion from fertilizers.
According to Somlyódy (1983), the majority of the nutrients came from leaching of fertilizers but, based on time series fertilization data, Joó (1986) stated that accumulated phosphorus stock in the soil was responsible for 50% of the total nutrient load of the lake and the leaching from periodic fertilization made up only 12%. Máté (1984) established that the phosphorus from fertilizers could reach Lake Balaton only in combination with the natural phosphorus content of the soil in proportion to the amount of transported suspended sediment and it did not exceed 2-3% of the total load.
Several studies (Miller et al., 1982;Pusztai, 1978;Sharpley et al., 1987) described that nutrients suspended in the water moved more easily than coarser solid particles of the soil. According to Stokes' law, the colloid-sized clay fraction (with phosphorus adsorbed to the clay surface) from eroded and water-suspended soil, and the light organic matter particles could reach surface waters. The coarser and heavier particles of the soil, which usually did not have a significant phosphorus nutrient content, settled down earlier and did not enter the lake. According to the calculations of Csermák and Máté (2004), about one-twentieth of the phosphorus entering the lake left through the Sió canal and approximately 95% of it stayed in the lake.
Due to the fact that Lake Balaton is considered a shallow lake, the frequent stirring of the sediment caused by wind-induced waves plays a very significant role in the development of both the dissolved phosphorus in the water and the total phosphorus content in the sludge. The bed conditions of the lake result in a specific phosphorus metabolism, different from the deep, stratified lakes. Oláh et al. (1977) found that the phosphorus content was generally higher in sediments of large bays lined with reeds and lower in sandy sediments. Within this, a significant part of the nutrient was bound to calcium due to the high CaCO 3 content of the sediment, while a smaller proportion was bound to aluminium and iron.
Based on the previous analyses of the environmental condition of Lake Balaton, we assume that environmental auxiliary variables (topographic, hydrodynamics, vegetation) could help to develop a more detailed map of the lake sediment.

Description of the sediment dataset of Lake Balaton
Between 1978 and 1984, Máté (1987;2004) sampled the upper 10 cm layer of the recent bed sediments of Lake Balaton in order to investigate the environmental status of the lake. Upon designating the locations of sediment samples, the aim was to obtain the greatest possible representativeness of the lake bed. Therefore, the location of points was predetermined on a 1:10.000 scale topographic map by triangulation method, then the exact locations of the sampling points were designated with geodetic methods under the EOV (Unified National Projection of Hungary). A total of 5.560 sediment samples were taken with an average density of 7-10 points/km 2 . The characteristics of the data set containing the results of sediment mapping are summarized in Table 1.
The traditional Hungarian method for determining phosphorus status in soils is ammonium-lactate acetic acid (AL) extraction. AL is an acidic solution (buffered at pH 3.75), which is also able to dissolve P reserves but it is does not suitable to determine total phosphorus content. AL-P is expressed in P 2 O 5 but it can be converted to atomic P by division it with 2.29.
In our analysis, we filtered the vectorized phosphorus data set in two steps. First, we identified the data records that were due to digitization errors (e.g. mistyping), and removed them from the data set. Secondly, we identified data records affected by measurement error using a spatial filtering approach, which was done as follows: (i) we computed the mean of the data points that are no farther than 1 km from each other, (ii) generated a raster layer showing the computed mean values, and (iii) identified and then removed the data points where the difference between the measured phosphorus content and the computed mean value was larger than 20 mg kg −1 . In total, 3.340 data points of phosphorus content remained after filtering, i.e. 60 per cent of the original.

Auxiliary variables used for spatial estimation of phosphorus content
For mapping the phosphorus content of lake sediment, we used three types of auxiliary variables: morphology, hydrodynamics and vegetation datasets. The short descriptions of auxiliary variables are summarized below ( Table 2).
The auxiliary variables used for the RK were produced from different archival data sources. Several GIS data processing methods were used to make these variables usable for mapping purposes.
i. Morphology: To create the digital bathymetric model of Lake Balaton, Hydrological Survey maps from 1975 with a scale of 1:25,000 have been used (Sass, 1979). The 22 map sheets were scanned and georeferenced in QGIS software. The original datum and projection of the survey were the Hungarian Stereographic Projection that had to be defined in QGIS using parameters provided by  Timár et al. (2003). After geo-referencing, all map sheets were converted to UTM/WGS84 Zone 33 N, and the bathymetric contour lines were digitized. Avoiding false raster values, breakpoints were densified to 25 m. Raster creation tool and the natural neighbour interpolation method were used with 100 m horizontal resolution. This output resolution was considered adequate for the objectives of the research, furthermore, it would have been unreasonable to produce an even more detailed bathymetric model considering the 1.000 m spacing of the original survey transects. Contour lines represented water depth and had to be recalculated to elevation above sea level. According to the notes presented on the maps, 0 depth regarded to 105.59 m above sea level of the Adriatic reference point, which was equal to 104.915 m above the recently used sea level (Baltic). In order to get elevation, interpolated raster values were subtracted from this reference altitude. From the digital elevation model, we derived some geomorphometric parameters commonly used in the environmental mapping (Laborczi et al., 2020;Pásztor et al., 2016;Tóth et al., 2016). These were as follows: analytical hill shading, slope, aspect, slope length and steepness factor (LS), channel network base level, vertical distance to channel network, valley depth, and relative slope position. ii. Hydrodynamic data: Hydrodynamic data related to sediment motion were also used as auxiliary variables, produced by Monte Carlo simulation. Indeed, the phosphorus adsorbing capacity of bed sediment can be closely related to its grain size composition. The bed sediment has been naturally sorted horizontally over a 100-year time scale for two main reasons. First, the force required to stir up finer and coarser fractions is different. Second, the settling velocity significantly varies as a function of grain size (Rákóczi, 1987). Both stirring and settling can be directly related to the shear stress acting on the lake bed; thus, we chose bed shear stress as the proxy variable of bed composition. We aimed to reconstruct its spatial and probability distribution for the 1970s by hydrodynamic modelling.
Wind-induced waves and currents shape the bed shear stress in Lake Balaton simultaneously; therefore, model calculations had three main components. The time-and space-varying wind field above the lake was determined using a fetch dependent algebraic atmospheric model developed in-house (Torma & Krámer, 2017). This unsteady wind field was imposed to drive the depth-integrated flow model ADH and the spectral, period-averaged wave model SWAN (SWAN Team, 2010). ADH is a finite-element solver of the 2D shallow water equations (Berger et al., 2011), used here on a fixed irregular triangular mesh. Two strips of 180 m cells resolved the nearshore which gradually transitioned to a cell size of 360 m offshore. The mesh of the SWAN model was also composed of irregular triangles with a spacing of 100 m nearshore and 400 m offshore. The flow and the wave models were run independently, without coupling them at timestep level. The interaction and combined effect of flow and wave motions near the bottom were considered by post-processing the shear stresses obtained from the two models using an algebraic procedure (for details, see Soulsby & Clarke, 2005) at a 30-minute time interval, from which statistics were computed. In both models geometry was set up based on the bathymetry model described in the previous subsection. The key meteorological and hydrodynamic model parameters were assumed to be stationary, and thus they were adopted from our modelling studies of Lake Balaton and Lake Fertő performed in the last decade (Homoródi et al., 2012;Lükő et al., 2020;Torma & Krámer, 2017). Based on the ECMWF ERA-5 Land Reanalysis from 1979 to the present (Muñoz Sabater, 2019), we found that the lake's wind climate has not changed significantly over this period. The probability distribution of the combined bottom shear stress was evaluated with a Monte Carlo simulation of a representative six-month period. The water level was set to 90 cm (relative to the Siófok gauge), approximating its mean in the 1970s. Following studies of  Rákóczi (1987) and Mehta et al. (2015), the upper, dilute sediments can be stirred by critical shear stress in a range estimated at 0.05-0.2 Pascal. Therefore, we finally calculated the spatial distribution of the frequency of exceeding these two threshold shear stress values to be used as auxiliary variables in the regression analysis.
(iii) Vegetation map: The vectorised vegetation map of Lake Balaton was also available, which was derived from an aerial photographic survey with a resolution of 8 cm. It shows the coverage of water plants including the reed in the lake. The categories of the vegetation map were determined according to the NBmR ÁNÉR and C-NÉR classification systems (Dömötörfy et al., 2005). Littoral vegetation including submerged and emergent plants are known to exert an influence on deposited sediment and water quality. For this reason, the presence or absence of these plants was used as a binary auxiliary variable.
Due to the various data sources, we resampled the above-listed auxiliary variables into a common reference system with a raster resolution of 100 m. In the next step, we standardized the auxiliary variables. To avoid multicollinearity in spatial modelling, we carried out principal component analysis (PCA) on the continuous auxiliary variables. In further spatial modelling, we used the principal components that explained about 99% of the total variation.

Spatial modelling of phosphorus content
We randomly partitioned the sample points of phosphorus content into a calibration dataset (90%) and a validation dataset (10%). We used both the calibration dataset and the validation dataset for modelling the spatial distribution of phosphorus content in the lake sediments and for checking the performance of spatial predictions respectively. We used regression kriging for modelling the spatial distribution of phosphorus content in the lake sediments. To evaluate model performance, we calculated residual standard error and multiple Rsquared of the model. To evaluate estimation accuracy we computed the following measures: root means square error (RMSE) and determination coefficient (R-square).

Results
The result of the PCA can be seen in Table 3. In summary, the variance of phosphorus values can be explained rather by the bathymetry variables, and the hydrodynamic variables may have less estimation power. As a result of this analysis, the exceedance probability of critical shear stress of 0.2 Pa was excluded from the auxiliary variables due to its low performance.
The multiple R-squared of the regression model was 0.23 and the residual standard error was 37.86 mg/kg. As a result, we concluded that the auxiliary variables alone are not strong enough to estimate phosphorus values between the measured points. However, regression combining with kriging proved to be effective in mapping the spatial distribution of phosphorus in the bed sediment. This was demonstrated by the R-squared value calculated on the validation dataset because it was significantly higher (0.81), indicating a good estimation performance. The most important parameters of our model and the model performance can be seen in Table 4.
The predicted phosphorus map (Main Map) produced with regression kriging has a resolution of 100 m and provides phosphorus concentration data as a continuous variable covering the whole area of the lake. The map shows that the phosphorus content follows a characteristic spatial pattern, which is discussed in more detail in the next section.

Discussion
Compared to the formerly compiled phosphorus map (Figure 1), the new one (Main Map) provides more reliable spatial information about the spatial variability and distribution of phosphorus concentration in  the lake bed. The previously compiled map was created by using the Thyssen polygon interpolation technique, which assigns interpolated value equal to the value found at the nearest data point. This yields artefacts in the resulting map, which spatially harms its interpretation. The new map (Main Map) is based on a novel digital mapping approach, which used regression kriging to predict the spatial distribution of phosphorus. Indeed, regression kriging is able to consider not only the available, spatially exhaustive auxiliary information but also the spatial correlation existing between the phosphorus data in spatial modelling and mapping. Therefore, the resulting map could give more reliable spatial information on phosphorus concentration, which explicitly reflects the spatially continuous property of the phosphorus concentration of bed sediments.
The up-to-date mapping, which was based on geostatistics, along with the archival phosphorus dataset from the 1970s and 1980s confirms our previous knowledge about the lake sediment and the state of eutrophication of the lake in many aspects. Eutrophication, in the case of lakes in temperate climate zones, is mainly induced by phosphorus loading. Parallel to the longitudinal axis of Lake Balaton, the phosphorus concentration of the sediment shows a decreasing trend from west to east. The spatial pattern of the map corresponds to the experienced distribution and, at the same time, reflects the trophic (oligotrophic, mesotrophic, eutrophic and hypertrophic) levels and water quality differences of the lake.
The hydrodynamic variables and the morphological variables derived from the digital elevation model provided important information for the mapping of the phosphorus content in Lake Balaton's sediment. The effect of these variables can be recognized in the lower phosphorus content of the sandy shoals all along the southern and eastern shores and in the case of the Tihany Strait. Similar shoals are absent along the north shore because it is not frequently exposed to strong winds, thus bed shear stresses are moderated. In the Tihany Strait, shear stresses are high due to the currents, which promote the coarsening of the bed material, resulting in lower phosphorus content. The analysis has confirmed that high phosphorus concentrations tend to be related to a finer sediment composition, which, in turn, is linked to a lower exceedance probability of the critical shear stress.
The phosphorus map (Main Map) demonstrates that the western basin of the lake including the Keszthely Bay is the most vulnerable to eutrophication. The effect of the River Zala, which is responsible for the most significant phosphorus load on Lake Balaton, can be recognized. The plume of the river is carried to the east with the prevailing longshore currents, in accordance with the mean flow field obtained with the hydrodynamic simulations. It can also be observed that there is significant phosphorus enrichment in the sediment along the northern shore of the Keszthely Bay, despite the absence of major surface water inflow. This can be explained by the fact that most of the phosphorus load enters the lake in water-suspended form through the River Zala, and windinduced motionthe combined effect of circulation and seichetransports parts of phosphorus bound to clay fraction to the northern shore of the bay, where it is slowly accumulated under the wind-sheltered calmer water conditions.
In the middle basins of the lake (around Balatonakali and Balatonszárszó), but especially in the eastern basins (Balatonfüred and Csopak area), there are areas with a phosphorus content of less than 25 mg kg −1 . Since sediments here are also characterized by a high CaCO 3 content, it can be concluded that, similarly to soils in Lake Balaton, soluble phosphorus may be fixed in the calcareous sediments and transformed into less soluble mineral forms.
Our phosphorus content map (Main Map) does not only show the conditions of the mid-1980s but can also provide important information for professionals who research the current state of the lake and plan future water quality management interventions. As a result of the comprehensive environmental protection program of the second half of the 1990s, the external phosphorus load of Lake Balaton decreased significantly. Over the last 40 years, the amount of phosphorus in the lake has been more or less constant, and the majority of the phosphorus has become fixed in the sediment. However, phosphorus fixed in the sediment might be partially released, thus the internal phosphorus load can significantly increase and can exceed the external load. As no similar survey was conducted in the last decades, very limited information is available about the present phosphorus level in the lake. In an unpublished report from 2020 (KEHOP-1.1.0-15-2017-00011' of the General Water Directorate) there are ∼150 sample points containing total phosphorus test results mainly from the western part of the lake. With the conversion the total P values to AL-P of this datasetbased on the conversion equations of Sárdi et al. (2009) and Trinchera and Baratella (2019)the mean level of phosphorus in the western basin of the lake shows about ∼230 mg/kg which is very similar to the level of 1980s.
The re-occurrence of summer algal blooms in 2019 and 2020 indicates an increase in the internal phosphorus load, which requires further analysis to halt possible water quality deterioration of Lake Balaton in the future. The detailed phosphorus map presented here may help identify the areas where internal phosphorus loading should be reduced by dredging or any other water quality management methods. Software QGIS 3.10 was used to produce the digital elevation model, ESRI Desktop ArcGIS 10.2.2 was used to processing the other input variables and compiling the final map. SAGA GIS 6.3.0 was used to calculate the variogram model and to carry out the regression kriging.