Evaluation of Macrophyte Community Dynamics (2015–2020) in Southern Lake Garda (Italy) from Sentinel-2 Data

Featured Application: The improvement of effective remote sensing-based approaches to map macrophyte features can provide a baseline of adequate spatiotemporal resolution for 21st century monitoring applications equipped to play a prominent role in the context of medium–large-scale management programs of ecological conservation and scientiﬁc research. Abstract: Macrophytes are of fundamental importance to the functioning of lake ecosystems. They provide structure, habitat, and a food source and are a required component in monitoring programs of lake ecological quality. The key aim of this study is to document the variation in spatial extent and density of macrophytes seasonally between 2015 and 2020 of the Sirmione Peninsula (Lake Garda, Italy), using Sentinel-2 imagery. In addition to this, our results were compared to previous data from imaging spectrometry; individual parameters affecting macrophyte communities were tested, and the possible effect of the COVID-19 lockdown on macrophyte colonization was evaluated. Satellite images allowed the mapping of the spatiotemporal dynamics of submerged rooted macrophytes in order to support monitoring of the shallow water ecosystem under study. Substantial changes were found in both spatial extent and density over the period from 2015 to 2020, particularly in 2019 when there was almost a complete absence of dense macrophytes. Variables found to inﬂuence the amount of macrophytes included transparency, chlorophyll–a, water level, winter wave height, and grazing by herbivores. A separate analysis focusing on areas associated with boat transit found a recovery in macrophyte coverage during the period of COVID-19 lockdown. The outcome of the study highlights a decline in the density of the macrophytes and a shift towards deeper areas compared to the situation in 1997. The area examined is part of an internationally important site containing the highest abundance and diversity of overwintering water birds in Italy. Exploiting satellite data at high frequency provided an insight to understand the dynamic changes and interactions with herbivorous birds, environmental factors, and anthropogenic pressures, revealing a delicately balanced and threatened ecosystem. to assess the spatiotemporal dynamics of the macrophyte community’s colonization.


Introduction
The issue of the dynamics of primary producers of lake ecosystems is acquiring more and more relevance because of increasingly marked effects of climate change on the functioning of these systems [1,2]. Indeed, the progressive reduction in the availability of water resources sharpens conflicting interests on the water uses, imposing increasingly robust assessments of macroclimatic, biological, and human drivers of lacustrine biota and spectral properties of reflectance spectra [29,30]. The major technical issues are presented by submerged or flooded vegetation, as water can alter the spectral characteristics of the images [31].
Lake Garda is the largest Italian Lake, one of the most important lacustrine ecosystems in Europe [23]. It is located in one of the most densely populated and industrialized regions of Italy representing an essential resource due to the multiple uses of its waters, as well as being a strategic area due to its ecological, conservational, and socioeconomic value [23,32]. In particular, this study is focused on the Sirmione Peninsula located in the southern portion of the lake. This area includes key habitats, such as common reed stands and extensive submerged macrophyte meadows, in addition to hosting the most diverse populations of overwintering water birds in Italy.
The colonization patterns of the aquatic vegetation of this portion of the Lake Garda were previously studied by Giardino et al. (2007) [33] and Bresciani et al. (2012) [6], and described the macrophyte coverage evolution around the Sirmione Peninsula. Giardino et al. (2007) [33] reported, for the 8-year interval analyzed by airborne imaging spectrometry, a decline in the submerged macrophyte meadows from 1997 to 2005 (from 72% to 52%) with a concomitant increase in sandy substrates. Bresciani et al. (2012) [6] used in situ and airborne hyperspectral data to analyze the variation of colonized substrates in relation to hydrological and physicochemical variables. They found a substantial modification of the structural complexity of macrophyte communities and of the colonized areas: 98% of macrophyte meadows were lost, mainly the dense beds (density over 70%), which were replaced by moderate to rare communities with lower density (range 10-40%). Transparency, water level fluctuation, trophic status, grazing by herbivorous aquatic birds, and number of boats were identified as the main factors to explain the spatiotemporal variability of macrophyte communities in the period investigated (13 years, from 1997 to 2010).
The present case study fits into a broader discussion about the global trends of macrophytes in lakes [34] which are currently undergoing an accelerating rate of decrease, with losses especially acute for submerged species in lakes larger than 50 km 2 . Previous evidence indicated that in some eutrophic lakes the current phytoplankton dominance can be a result of submerged macrophyte decline rather than the cause [35,36], further suggesting that these phenomena are a conclusion of a process of years (over perhaps 10-100s), driven by many pressures rather than a rapid change. In fact, several studies have reported the occurrence of long-term trends of macrophyte diversity and abundance loss in many shallow lakes affected by eutrophication [37,38] and spatial and seasonal variation of macrophytes over the long term [7,39,40]. At the same time, the reduction of human-mediated pressures exerted on lake systems can stimulate the recovery of macrophytes, as in the case of Lake Constance [41]. Here, the reduction of phosphorus inputs progressively stimulated the recovery of macrophytes, with a significative increase in species richness and charophytes dominance.
To add new knowledge on the topic, the present work aims to update and expand the analysis and assessment performed by previous studies, adding the use of the recent Sentinel-2 for the assessment of the evolution of submerged aquatic vegetation cover and abundance in Lake Garda. These patterns will be discussed in relation to the main factors that can influence seasonal and interannual variation in the macrophyte community, such as hydrological and meteo-climatic factors, water quality, biological interaction, and tourism pressure.
The specific objectives of the study are (i) to produce seasonal maps of macrophyte density and cover for the period 2015-2020 with Sentinel-2 imagery; (ii) to analyze the spatiotemporal variation in macrophyte colonization; (iii) to test individual parameters in the context of temporal and seasonal macrophyte change; (iv) to analyze interannual changes by comparison with the previous results from imaging spectrometry; and (v) to evaluate the possible lockdown effect on macrophyte colonization.

Study Area
Lake Garda is the largest Italian lake, with a surface of 368 km 2 , and is located in the continental bioregion of Europe. It is characterized by a maximum depth of 350 m and a volume of 49 km 3 . Lake Garda contains 30% of Italy's whole natural and artificial lake water resources and represents an essential strategic water supply for agriculture, industry, energy, fishing, and drinking [42]. This oligotrophic lake [43] of glacial origin is located at an altitude of 65 m a.s.l. and counts a total of 25 tributaries, including the main inflow River Sarca located at the northern edge of the lake. The only emissary of Lake Garda is the River Mincio in the southeast of the lake, with average discharge in the winter season of 25 m 3 s −1 , while removal for irrigation ranges between 40 and 80 m 3 s −1 [44]. The hydrometric level of the lake is influenced by the regulation of its flow at Salionze dam.
The catchment of Lake Garda is composed of 60% sedimentary rocks, while the remaining 40% is equally distributed between crystalline rocks and secondarily deposited glacial and fluvial sediments [45]. In addition, the catchment is relatively small compared to the surface area of the lake (6:1). This characteristic, combined with low annual precipitation, has led to a long theoretical water renewal time (27 years) [46]. Moreover, Lake Garda is characterized by long periods of incomplete vertical winter water circulation, which are interrupted by full mixing of the water column after the occurrence of harsh winters [47]; consequently, Lake Garda can be defined as warm monomictic and oligomictic.
Regarding the physical-chemical characteristics, the average value of Secchi disk changes considerably on a seasonal basis: in summer it can vary between 4-5 m, while in winter reaches up to 16 m [48]. Furthermore, the detectable concentration of total phosphorus in the epilimnion is below or around 10 µg L −1 , while the average chlorophyll-a (Chl-a) concentration is about 3 mg m −3 [48]. Lastly, according to Premazzi et al. (2003) [49], the suspended particulate matter (SPM) concentration spans a winter-summer range between 0.1 and 5.5 g m −3 (average value 2.5 g m −3 ).
With respect to morphology, Lake Garda can be divided into two different portions, separated by an underwater ridge connecting Punta San Vigilio with the Sirmione Peninsula, our study area ( Figure 1): the largest part is located to the west and is characterized by a greater maximum depth (350 m), while the east part is much less extensive (representing less than 7% of the total volume of the lake) and is characterized by shallower waters (maximum depth of 80 m).
Sirmione Peninsula is one of the most popular tourist destinations in the Lake Garda region (with more than a million tourists per year), with attractions such as the ruins of Roman villas and the Scaligeri castle, as well as hot springs and spa facilities. This area extends for about 4 km into the lake and its shores are characterized by moderate slopes inhabited by diverse species of macrophytes [6]. Less urbanized shores are characterized by reedbeds dominated by Phragmites australis (Cav.) Trin. ex Steud., while submerged aquatic vegetation are dominated by several Chara species (Chara globularis Thuillier 1799, C. polyacantha A. Braun in Braun, Rabenhorst & Stizenberg 1859, C. tomentosa Linnaeus 1753), Vallisneria spiralis L., Lagarosiphon major (Ridl.) Moss, and a few taxa of the Potamogeton genus [6,50].

Sentinel-2 Images Processing and Validation
The twin satellites of Sentinel-2 (S2A and S2B) carry on board the MultiSpectral Instrument (MSI) optical sensor, which acquires images in the visible-near-infrared (VIS-NIR: 13 bands from 442 to 2202 nm) spectrum, with a spatial resolution of 10, 20, and 60 m (depending on the band), at a frequency of about 5 days were used in this study to assess the spatiotemporal dynamics of the macrophyte community's colonization. . 2022, 12, x FOR PEER REVIEW 5 of 22 Figure 1. Sentinel-2 image of Lake Garda in true color (23/03/2019). Above on the left is the location of Lake Garda in Italy. The blue and green rectangles represent the two areas analyzed to extract turbidity values (NTU). The red box represents the study area (Sirmione Peninsula). Below on the right shows an enlargement of the Sirmione Peninsula (grey): the yellow area represents the mask used to create the bottom coverage maps (3.97 km 2 ), while the lines with black squares represent the transects used for assessing the macrophytes gradients from the coast.

Sentinel-2 Images Processing and Validation
The twin satellites of Sentinel-2 (S2A and S2B) carry on board the MultiSpectral Instrument (MSI) optical sensor, which acquires images in the visible-near-infrared (VIS-NIR: 13 bands from 442 to 2202 nm) spectrum, with a spatial resolution of 10, 20, and 60 m (depending on the band), at a frequency of about 5 days were used in this study to assess the spatiotemporal dynamics of the macrophyte community's colonization.
A total of 34 images available between 2015 and 2020 were downloaded from the Copernicus Open Access Hub [51] or from the ONDA catalogue [52] (Table 1). In detail, for each year except 2015 (the S2 mission began in July 2015), six images were chosen: two for each season characterized by macrophyte presence (from spring to autumn). For this work, the first seven bands were used: from 442.7 to 782.8 nm for S2A and from 442.3 to 779.7 nm for S2B. The bands of each image were then resampled to the spatial resolution Above on the left is the location of Lake Garda in Italy. The blue and green rectangles represent the two areas analyzed to extract turbidity values (NTU). The red box represents the study area (Sirmione Peninsula). Below on the right shows an enlargement of the Sirmione Peninsula (grey): the yellow area represents the mask used to create the bottom coverage maps (3.97 km 2 ), while the lines with black squares represent the transects used for assessing the macrophytes gradients from the coast.
A total of 34 images available between 2015 and 2020 were downloaded from the Copernicus Open Access Hub [51] or from the ONDA catalogue [52] (Table 1). In detail, for each year except 2015 (the S2 mission began in July 2015), six images were chosen: two for each season characterized by macrophyte presence (from spring to autumn). For this work, the first seven bands were used: from 442.7 to 782.8 nm for S2A and from 442.3 to 779.7 nm for S2B. The bands of each image were then resampled to the spatial resolution of 10 m. All images were chosen without clouds and cloud-shadows, without sun-glint and other radiometric noise. All images were downloaded as Level 1 (L1); consequently, the first processing step was the conversion of top of atmosphere (TOA) reflectances into TOA radiances, through the SNAP tool "ReflectanceToRadianceOp" for allowing the atmospheric correction to be performed. To this aim, the Second Simulation of the Satellite Signal in the Solar Spectrum code vector version (6SV ver. 2.1) [53] was adopted in agreement with previous validating studies [54,55]. The 6SV was parametrized with aerosol values (AOD 550 ) referring to the acquisition dates and collected from a nearby AERONET station [56] or from NASA data [57] when AERONET data were not available. The 6SV-derived atmospherically corrected products were converted into remote sensing reflectances (Rrs) dividing by π.
The robustness of the atmospheric correction performed by 6SV on S2 images, already demonstrated in Bresciani et al. (2018) and Cazzaniga et al. (2019) [54,55], was further tested over land. To this aim, an evergreen wooded area, a sandy beach, and some rooftops were selected as invariant reflectance targets. The 6SV-derived reflectances of the targets were compared for the 34 atmospherically corrected images. Three regions of interests (ROIs; of a 2 × 2 pixel size) were created, one for each invariant surface ( Figure S1), for computing the average reflectances for each of the 34 S2 images.
Rrs products were then used as input for the BOMBER bio-optical model (Bio-Optical Model-Based tool for Estimating water quality and bottom properties from Remote sensing images) [58] for estimating the substrate coverage. BOMBER parameterization was based on the in situ inherent optical proprieties (IOP) of the water column and reflectances of bottom substrate collected in Lake Garda and previously used and validated in different work [54,[59][60][61]. The BOMBER code was applied in shallow water mode for each pixel included in the predefined mask of 3.97 km 2 (Figure 1), defining the shallow waters of Sirmione Peninsula reaching a maximum depth of 10 m. This mask was created to avoid areas too close to the coast that might contain mixed pixels (water-land interface).
The BOMBER outputs provide the fractional cover of bottom substrates through three classes: sand (b0), rocks (b1), and rooted macrophytes (b2), and in addition produce the optical closure error maps. For each image, the least represented class is rocks (b1) (not shown here), so we decided to combine the two classes b0 and b1 to obtain products characterized by two main classes: bare sediment (BS, b0 + b1) and rooted macrophytes (RM, b2).
The maps related to each season were merged in average maps representative of the season, as suggested in Roelfsema et al. (2018) [62]. In detail, starting from the 34 maps, 17 seasonal maps were obtained (three for each year, except for 2015) to represent spring, summer, and autumn. The maps were composed by averaging each pixel value of the individual maps, considering both classes: BS and RM. Finally, in order to analyze the spatiotemporal variations of macrophyte communities, the main class RM was partitioned into two subclasses: dense macrophytes (D RM ) and sparse macrophytes (S RM ). Assuming that each pixel has a value between 0 and 1 (for both BS and RM), thresholds were estimated as follows: D RM (RM > 0.66), S RM (0.33 < RM < 0.66), and BS (RM < 0.33).
The validation of the BOMBER outputs was performed analyzing optical closure error maps in two predefined ROIs: the first in an area (35 pixels, located in the west sub-basin) classified as colonized by macrophytes throughout the dataset, the second characterized by a constant absence of submerged vegetation (90 pixels, located in the east sub-basin) ( Figure S1). To identify these two ROIs, all maps depicting macrophyte meadows (RM > 0.33) and all maps relating to the bare sediment (RM < 0.33 or BS > 0.66) were intersected. Once identified, the mean value of the ROIs and standard deviation were calculated for the entire dataset.
In order to study the pattern of bottom cover classes in relation to distance from the coast, six transects were created, three on each side of the coast (Figure 1). Along each transect, every 100 m, ROIs (3 × 3 pixels) were created for a total of 38 ROIs. The ROIs were then applied to the seasonal maps to obtain the dominant bottom cover class. In detail: transect T1 (500 m) represents a port area characterized by large reed beds; transect T2 (400 m) is in the proximity of a reed bed area; transect T3 (500 m) crosses the route of the ferry that brings tourists to Sirmione (in Figure 2, a semicircle area lacking macrophytes can be observed in the northwest part of the mask of almost all seasonal maps); transect T4 (500 m) is located near a port and an urban area (car parking); finally, transects T5 (700 m) and T6 (1200 m) are located near port areas with a high presence of discharges from drains.  Finally, we investigated the possible impact of the lockdown period (March to May 2020) on the density of macrophytes stands. In Italy, the first case of COVID-19 was detected on 21/02/2020, in Codogno (Lombardy). As a result, a number of restrictions were implemented to contain the spread of the virus. These limitations were finally extended to the entire national territory from 10 March [63][64][65], leading to a drastic reduction of anthropogenic impact on the environment (e.g., travel restrictions, traffic reduction, and closure of many businesses, including the entire tourism sector). In our study, the areal ranges of the D RM and S RM classes were compared between the two spring images of each year (except for 2015): i.e., late March vs. late April-early May.  [6]. The slight mismatch of months between some of the dates compared should not affect the extent of macrophyte meadows too much, since many authors estimate the vegetative peak between July and August [66,67], while plant senescence begins towards the end of September [67,68]. Accordingly, the MIVIS-derived maps (1997-2010) were compared to the summer seasonal maps (2015-2020) of this study, which represent the period of maximum macrophyte abundance and density.

Ancillary Data
Data on chlorophyll-a, Secchi disk depth, and total phosphorus and nitrate were obtained from national monitoring programs implemented and reported by ARPA Veneto [69]. In addition, satellite estimates of turbidity (NTU) were extracted from two areas (of eight pixels each) east and west of Sirmione from 2016 to 2019 using data from the Lakes CCI project (Figure 1) [70]. Data on bird numbers were taken from the annual census of aquatic birds for Lake Garda carried out every January [71][72][73][74][75][76]. Herbivorous birds were dominated in 2019 by the Eurasian coot (Fulica atra) (13,295), mallard duck (Anas platyrhynchos) (1442), and red-crested pochard (Netta rufina) (123) [75]. Total number of piscivores and herbivores were calculated following Özgencil et al. (2020) [77]. Data on daily lake level from the hydrometric station located at Peschiera was obtained from AIPo [78]. Tourist numbers were obtained from the province of Brescia website [79]. Climatic data were obtained from ERA5, the fifth-generation ECMWF reanalysis for the global climate and weather (hourly data on single levels) [80]. Data used for analysis included 10 m u-component of wind, 10 m v-component of wind, and surface solar radiation downwards. The wind components were used to calculate wind speed and direction which, together with fetch, allowed the prediction of wave height (see equation 15 in Carter, 1982 [81]). All ERA5 data documentation is available online: [82].

Statistical Analysis
There are a large number of parameters that can potentially affect macrophytes, making it difficult to extract their influence over a short time period. For this reason, one of the approaches to analysis has been more descriptive, with a testing of individual parameters in the context of temporal and seasonal change. Nonparametric multiplicative regression (NPMR) models [83] for total macrophytes (RM), dense macrophytes (D RM ), and sparse macrophytes (S RM ) were used to examine the response to meteo-climatic and environmental parameters. Variables that were included alongside time and season included chlorophyll-a, Secchi disk depth, total phosphorus, nitrate, the sum of herbivorous and piscivorous birds, mean lake level, the winter (DJF) mean wave height, wind vectors u and v, tourist numbers, and surface solar radiation downwards. Statistical analysis of wind data was carried out in the package "openair" of R [84,85]. NPMR defines response surfaces using predictors in a multiplicative rather than in an additive way. This method is better at defining unimodal responses of ecological data than other methods such as multiple regression [83]. NPMR was applied using the software HyperNiche version 2.3 [86]. NPMR models were produced by adding predictors stepwise with fit expressed as a cross-validated R 2 (xR 2 ) which can be interpreted in a similar way as a measure of fit as a traditional R 2 .

Evaluation of S2-Derived Products
The results of the comparison of the spectral signatures (first seven bands: 442-780 nm) as produced by the 6SV for the three invariant surfaces attested to the quality of the atmospheric correction: for all the S2 images the coefficient of variation is in fact low: 0.26 for sandy beach, 0.23 for evergreen vegetation, and 0.14 for rooftops. In particular, the bands that present the highest coefficients of variation are the b6 and b7 bands, i.e., those that are closest to the near-infrared. As a result, only the first five bands of Rrs product (442-705 nm) were used as input for the BOMBER code. Furthermore, by removing the last two bands, the coefficients of variation of the invariant surfaces become 0.25 for sandy beach, 0.15 for evergreen vegetation, and 0.1 for rooftops. By adding the previous validating studies with S2 over Lake Garda [54,55] we might hence conclude that the atmospheric correction carried out by 6SV is robust and coherent in order to obtain temporal Rrs that are comparable.
The mean value of the error maps of the BOMBER products estimation for the two ROIs (bottom and macrophytes) was calculated for each year and for each season ( Figure S2) and were considered sufficient for the purposes of this study. In particular, on an annual basis, the average value of error is lower than 20%, while on a seasonal basis it has always been lower than 15%. More generally, considering all S2 images available, the mean error value is 7% and 10% for bottom ROI and macrophyte ROI, respectively. In addition, the mean error value of the bottom ROI is always lower than the macrophyte one. for the DRM class; while for the remaining transects (T2-T5-T6), no particular interannual differences were found.   The western portion of the Sirmione Peninsula was more favorable to macrophyte growth, while the southern portion of the eastern sub-basin was mainly characterized by bare sediment (BS), resulting in a spatial difference between the two sub-basins ( Figure 2). This is supported by transect analysis where the frequency of the two macrophyte subclasses (S RM and D RM ) in the seasonal maps is higher in the western transects than in the eastern transects. Specifically, considering all ROIs and the entire study period (2015-2020), in the western transects the average frequency of the RM class is 74.71%, while in the eastern transects it is 50.20%. T1 is the most favorable transect for macrophyte growth (92.94%), while T5 and T6 are the most unfavorable (42.86% and 34.80%, respectively). The eastern side, especially areas that are shallower-closer to the shore and also towards the main shoreline to the southeast of the peninsula, where the slope is more gentle, tend to experience the largest areas of bare lake substrate over the time series (see bathymetry map in Giardino et al., 2014 [59]). In contrast, the western side is somewhat in the lee of the peninsula head and can be seen to consistently retain a greater proportion of macrophytes throughout the seasons compared to the eastern side. In addition, areas with more dense macrophytes are typically located in deeper water off the peninsula (Figure 2). In fact, in all transects, an increase in the frequency of the D RM class can be detected in relation to distance from the coast. With the exception of T5 and T6, the D RM class becomes dominant at the end of each transect (surpassing even the S RM class). In particular, in T1 at 200 m, the D RM class overtakes the BS class, in T2-T3-T4, this occurs at 300 m, while in T5 and T6 the dominant class is BS for the entire length of the transects (with the only exception of the second ROI of T5, i.e., 200 m, where the S RM class is slightly more frequent).

Spatiotemporal Dynamics of Macrophyte Communities
Lastly, comparing 2020 (lockdown) with previous years, the following results were obtained: in T1, 2020 is the only year in which the D RM class is present both in spring and summer at 100 m; in T3, the improvement is found at 200 m with the presence of S RM both in spring and summer; in the first 300 m of T4, 2020 is the most favorable year, especially for the D RM class; while for the remaining transects (T2-T5-T6), no particular interannual differences were found.
The possible effect of the lockdown period (March to May 2020) on macrophyte meadow growth was estimated by comparing spring images from each year. Comparison of the early spring image of each year with the late spring image shows that from 2016 to 2019 there was always a decline in total macrophyte (RM) coverage: once for both subclasses (2016), sometimes only for D RM (2017 and 2018), and finally only once for S RM (2019). The only year characterized by a slight increase in both subclasses and, consequently, in total macrophytes was 2020 (Table 3). Table 3. Comparison of area (km 2 ) occupied by macrophyte classes (S RM , D RM , and RM) for the two spring images (early spring vs. late spring) of each year (2016-2020). In green are positive values (increase of area) and in red are negative values (decrease in area).

Year-(Day)
Area ( The interannual macrophyte cover trends are reported in Figure 4. Bare sediment bottom cover was greater than 50% in both 2005 and 2016, followed by 42% in 2019, while in the rest of the years investigated it was, on average, 37%. Macrophyte cover reached the highest values in 1997 and 2017 (72% for both), followed by coverage in the range from 65% to 69% in 2010, 2015, 2018, and 2020 ( Figure 4). Higher variation was found in the density of the macrophytes in the years analyzed, with dense macrophytes cover being higher compared to sparse density in 1997, 2015, 2017, and 2018 ( Figure 4).   Table 4 lists significant models where the environmental variables increased the xR 2 when added to NPMR. Variables that were included alongside time and season included Chl-a, Secchi disk depth, herbivorous birds, mean lake level, and the winter (DJF) mean wave height. All of the models were for total macrophyte area (RM), while only one model was included for dense macrophytes (DRM) with Chl-a, albeit with a p value of 0.06. No significant models were produced for sparse macrophytes (SRM). The xR 2 ranged from 0.51 for Chl-a to 0.39 for Secchi disk depth and winter mean wave height. The sensitivity val-  Table 4 lists significant models where the environmental variables increased the xR 2 when added to NPMR. Variables that were included alongside time and season included Chl-a, Secchi disk depth, herbivorous birds, mean lake level, and the winter (DJF) mean wave height. All of the models were for total macrophyte area (RM), while only one model was included for dense macrophytes (D RM ) with Chl-a, albeit with a p value of 0.06. No significant models were produced for sparse macrophytes (S RM ). The xR 2 ranged from 0.51 for Chl-a to 0.39 for Secchi disk depth and winter mean wave height. The sensitivity values, indicating the importance of the parameters, ranged from 0.05 for winter mean wave height (for total macrophytes, RM) to 0.23 for Chl-a (for dense macrophytes, D RM ). Several environmental parameters were not found to increase the xR 2 : total phosphorus, nitrate, wind vectors u and v, tourist numbers, surface solar radiation downwards, and the sum of piscivorous birds. Table 4. Results of NPMR (nonparametric multiplicative regression) models for total macrophytes (RM) and dense macrophytes (D RM ). No successful regressions were found for sparse macrophytes (S RM ). xR 2 = cross-validated R 2 ; Ave. size = average neighborhood size; Tol. = tolerance; Sen. = sensitivity; Secchi = Secchi disk depth (m); Level_mean = mean lake level (m); Sum Herb = sum of herbivorous bird species; DJF_Wave H = winter average wave height; Chl-a = chlorophyll-a (mg m −3 ). Successful models (additional variation explained by an increase in xR 2 ) were not found for total phosphorus, nitrate, wind vectors u and v, tourist numbers, surface solar radiation downwards, and sum piscivorous birds.  Table 5 interprets the changes that have occurred over the period for macrophytes, and the key parameters noted above. Of note is the fact that an increase in herbivorous birds (estimated from a winter (January) census) tends to be followed by a decline in macrophytes, and if the decline is substantial, this in turn leads to a reduction in bird numbers the following winter. For example, in 2016 the lowest total cover (RM) was recorded for all macrophyte growth seasons considered (spring, summer, and autumn) (Figure 3a). This decline followed directly after an increase in herbivorous birds the preceding winter ( Figure 3b); in addition, high northeasterly winter wind speeds were recorded ( Figure S3). Winter herbivorous birds declined in 2017, and macrophyte recovery followed in 2017 and 2018, after which the macrophytes declined again in 2019, with a notable reduction of the area covered by dense macrophytes-reaching the lowest values recorded for the dataset in summer and being almost absent in autumn (0.025 km 2 ) (Figure 3a). This was preceded by the highest numbers of winter herbivorous birds (Figure 3b). However, the Chl-a concentrations were also high for this period, as was the water level (Figure 3c,d). The minimum seasonal Secchi disk depth of 5.4 m for the time series was also recorded in the summer of 2019. In addition, a storm event occurred on 5 May, with maximum wind speeds recorded for the time series up to 7 m s −1 . Macrophyte recovery, notably the area of dense macrophytes, occurred in 2020 with lower herbivorous birds but also a decline of 55.7% of tourist numbers, owing to the coronavirus epidemic. Two regions of interest east and west of the peninsula were examined for differences in turbidity using the lakes CCI dataset to see if localized differences in water quality might be influencing the macrophytes (Figure 1). No large changes over the years available (2016-2019) were evident, although the western side tended to have slightly higher turbidity, despite it having typically more macrophytes (Figures S4 and 2). Table 5. Potential drivers, responses, and implications of changes in macrophyte cover in Lake Garda from 2015 to 2020 (MAM = March, April, May; JJA = June, July, August; SON = September, October, November). Upward/downward/sideways arrows indicate increase/decrease/stable trend in herbivores (blue) or macrophytes (green).

Year
Potential Drivers Response Implication for Following Year
Spring lake level was 1.06 m.
Macrophyte cover normal.
Lower herbivore bird grazing earlier in the year may have allowed an increase in macrophyte growth in subsequent months and lead to more foraging the following year.

2016
Herbivore population increases to 13,908. High winter wind speed recorded (red in Figure S3). Summer Chl-a increases to 4.3 mg m −3 . Spring lake level was 1.03 m.
Lowest total macrophyte cover in timeseries for MAM, JJA, SON recorded.
Lower macrophyte cover and lower density in SON will mean less foraging for birds.

2017
Herbivore population declines to 10,327. Winter wind mostly calm. Summer Chl-a was 1.9 mg m −3 .
Spring lake level was 1.08 m.
An increase in Piscivores recorded, Macrophytes show recovery in subsequent months.
Recovery of macrophytes will mean more foraging for birds.

2018
Herbivore population increases to 12,512. Winter wind mostly calm. Summer Chl-a was 3.0 mg m −3 .
Spring lake level was 1.07 m.
Macrophytes remain at similar or slightly lower levels.

2019
Herbivore population continues recovery and reaches highest level of 15,074. Winter wind mostly calm but storm in May. Spring had higher water levels (1.29 m). Summer Chl-a increases to highest for the timeseries of 4.8 mg m −3 .
Area of dense macrophytes is reduced to minimum of the timeseries. Total macrophytes in summer lower than previous two years.
Reduction in density of macrophytes may have implications for foraging.

2020
Herbivore population declines to 11,326. Winter wind mostly calm. 55.7% reduction in tourist numbers likely to reduce erosion/turbidity from boat wakes. Summer Chl-a was 2.0 mg m −3 .
Spring lake level was 1.2 m.
Area of dense macrophytes recovers. Total macrophytes in summer return to higher levels. Piscivores decline.

Discussion
The results of the present study confirmed the capability of Sentinel-2-derived products in tracking the bottom coverage of optically shallow waters and demonstrated the existence of striking inter-and intra-annual variations in the spatial cover patterns of macrophytes. In this context, remote sensing has proven to be a valid tool for investigating seasonal variation in macrophytes over a wide study area (3.97 km 2 ); processes that would be impossible to adequately assess using classical limnological approaches (e.g., by macrophyte survey along transects, as implemented under the Water Framework Directive) [9,87].
The comparison between our seasonal maps and those of the two previous works [6,33] shows that 1997 was the most favorable year for macrophyte growth (the RM class covered 72% of the entire mask), especially for the D RM class (69%) that was located mainly near the coast. These percentages of vegetation coverage were also reached in other years (e.g., 2010, 2015, 2017, 2018, and 2020), but the substantial difference concerns the D RM class, which after 1997 never reached 50%, even in the most favorable years. In addition, when present, it is mostly located far from the coast. As a result, a gradual change in macrophyte quality/density is observable from 1997 to the present day, and a shift of dense macrophytes to deeper areas. This is in line with previous works highlighting high intrinsic dynamic rates of macrophytes in deep lakes [4,39].
As indicated in the methods section, the approach to the interpretation of the drivers behind the changes in composition of the macrophytes over the time series is challenging and was focused on a more descriptive listing of potential drivers coupled with testing of individual parameters in the context of temporal and seasonal change.
The areas that tended to be more dynamic in terms of macrophyte loss and regrowth tended to be shallower areas, especially close to the main shore. This indicated that the influence of exposure to waves was one of the principal drivers (favoring the damage, breakage, or uprooting of macrophytes) as shallower areas of 3 m or less are often exposed to wave action [15,16]. In addition, the area to the west, partially in the lee of the head of the peninsula, had more macrophytes and experienced less loss, and recent models have indicated that this area may be more protected given the direction of wind vectors [88].
Transect analysis confirmed that the west sub-basin is more conducive to macrophyte growth than the east sub-basin, which has a large area to the southeast characterized by ports and many discharges from drains. In spring and autumn this area is dominated by the BS class. Specifically, transects T5 and T6 (located close to port areas with a high presence of discharges from drains) show how this class dominates the others at any distance from the shore (with some exceptions in T5, where the S RM class is also present). The presence of macrophytes is only evident in the summer during peak vegetation.
In addition, no particular differences emerged between 2020 (lockdown period) and other years, so it is plausible that in this area drains play a more important role than the presence/absence of boats. In contrast to these transects, in T1-T3-T4 (characterized by the presence of ports and ferry docking), improvements emerged in 2020 in ROIs located closer to the coast. In detail, looking at the values of 2020, T1 is characterized by the presence of D RM in the first 100 m (both in spring and summer), in T3 emerges an improvement especially at 200 m (both in spring and summer), while in T4 the improvement is visible in the first 300 m (both in spring and summer). Further evidence relating to the impact of boats on the growth of macrophytes can be observed in the first 100 m of T1 (one of the largest ports in Sirmione) where the BS class is present only in summer, i.e., when there is the greatest transit of boats.
Wind and wave erosion were always also likely to influence this zone given its location at the south of the lake and the prevailing northerly winds, coupled with a potential fetch of 50 km. The wind event in May 2019 was estimated using ERA5 data as no in situ data was available but estimated wave heights of 0.4 m were close to previous estimates on typical high-wind-speed days of 0.5 m [88]. The wind event is likely to have played a role in the loss of dense macrophytes in 2019, and local newspapers recorded the storm as causing substantial crop damage and felling trees. However, other factors are also likely to be important, as examining the spatial distribution map prior to this event in spring indicated that the cover was already low. Seasonal wind components were not significant when included in the NPMR models but estimated wave height in winter was, indicating the importance of including wind speed, direction, and fetch to calculate wave height.
The cyclical behavior of herbivorous bird increase and macrophyte decline, in turn followed by a decline in herbivores with subsequent macrophyte recovery, suggests an ecological balance at play. The minimum of macrophytes in 2016 and 2019 corresponded to the years where winter herbivore numbers were highest. Lake Garda is frequently recorded as one of the sites in Italy with the highest abundance and diversity of overwintering birds [76]. The influence of waterfowl in grazing macrophytes varies from lake to lake and has been implicated in hindering macrophyte recolonization in some restoration projects [19]. Years with high macrophyte abundance have been found to increase total bird numbers as well as herbivores, whereas piscivores were lower during such times [77]. This study also found a similar trend, including that of piscivores, which could suggest a benefit in predation in the absence of refuge which is also likely to be a cyclical benefit as macrophytes generally provide habitat and increase the abundance and diversity of fish and macroinvertebrates [89][90][91]. Also worthy of further study are the relative roles of invasive crayfish and fish predation on the macrophytes that could also contribute to eliciting a cyclical response in the macrophytes similar to that forwarded for herbivorous birds [19,92,93]. It is noteworthy that the levels of coverage and especially density of macrophytes observed in 1997 (particularly in shallower areas) have not been observed since, and that this coincides with the arrival of the invasive crayfish Orconectes limosus in 1998 and Procambarus clarkii in 2003 [93]. Their numbers are also likely to decline in response to a decrease in macrophytes and increasing piscivorous birds observed here.
It is likely that no relationships were found between sparse macrophytes (S RM ) and environmental variables because their variation was dominated by the transition to and from dense macrophytes (D RM ).
Of the indicators of trophic state, both chlorophyll-a and Secchi disk depth were significant in the models for predicting total macrophytes but not total phosphorus. The loss of dense macrophytes in 2019 corresponded to a summer where the highest Chl-a and lowest transparency were recorded. While the seasonal mean of 5.4 m still appears high, a value of 3 m was recorded on 7 May 2019; in comparison, excluding 2019, the average Secchi disk depth in May from 2015-2020 was 10.6 m. Using these extreme values to estimate the depth of colonization of Charophytes would estimate depths of 4.7 m (for 3 m Secchi) to 17 m (for 10.6 m Secchi) [94]. While not an attempt to accurately characterize the light climate over the growing season, it indicates that 2019 was a year with significant pressure that could have contributed to the loss of dense macrophytes. These communities appeared to be lost from areas at a distance off the peninsula in deeper water, indicating that the changing light climate may have been responsible given the exponential decline in light with depth ( Figure 2). It is possible that the storm event on 5 May increased both turbidity and released nutrients from the sediment layer via direct resuspension [95]. Such storms may have a double effect of reducing macrophytes in shallower water through wave action, while those in deeper water may suffer from poorer light climate through increased turbidity and nutrient addition causing phytoplankton and epiphytic growth with resultant shading [96].
The higher water level recorded in spring and summer in 2019 will also have had a direct influence on macrophytes through increased attenuation of light alone. For example, in a clearwater oligotrophic lake, if the water above a substrate is increased from 2 to 3 m depth then this would decrease the percentage of light reaching it from 51 to 36%, whereas for a substrate where water is increased from 7 to 8 m depth, then this would decrease the percentage of light reaching it from 9 to 7% [97]. In 2016, similar Chl-a concentrations in the summer (JJA) and comparable water level in spring and summer periods (MAM and JJA) to 2019 were recorded and were reflected in the prevalence of sparse macrophyte meadows also in the summer 2016 ( Figure 3).
While there are many potential drivers of the interannual variation in macrophytes, the influence of climate change is likely to be prominent. The principal change in Lake Garda is that warming winters are reducing the frequency of winter overturn, with no substantial mixing event since 2005, resulting in altered nutrient cycling and higher nutrients in the hypolimnion [47,98]. There is some evidence that this is leading to a change in the phenology of the phytoplankton with a shift from spring to summer maxima [61]. Summer is the period with maximum growth and abundance for macrophytes, so such changes to the light climate during this period can be expected to negatively affect the macrophytes.
There is also an increase in intensity of storms predicted with climate change [99,100], and this will likely alter nutrient cycling as well as increase physical exposure pressure and, if accompanied by significant rainfall, increase the lake level increasing light attenuation. The higher level of nutrients in the hypolimnion also represents a threat in the event of an overturn event if a cold windy winter occurs. The amount of nutrients released could lead to bloom events, less macrophytes, and less foraging for herbivorous birds.
While the number of tourists did not appear to be a significant factor in the analysis of the whole study area, it is likely to be a hidden significant driver, perhaps less through the addition of nutrient load within the catchment [47] than through increased boating pressure around the peninsula as revealed by examining specific transects subject to boating traffic, as discussed above. In this context, the observed drop in tourist numbers by 56% in 2020 owing to the coronavirus epidemic may have reduced such pressure and contributed to the recovery in macrophytes observed. Similar results emerged in Braga et al., 2020 [101], where the impact of lockdown on the Venice lagoon was studied. In detail, the drastic decrease in boats led to a reduction in erosion and particle resuspension, consequently promoting an increase in water transparency (a key factor for macrophyte growth).
Adequately understanding and describing the contribution macrophytes make to lake food webs has long remained a challenge to limnology, not least because of the large diversity of macrophytes of various growth forms and stages as well as grazing being carried out by multiple taxonomic groups [92]. This study demonstrated the dynamic seasonal and interannual variation in macrophytes, attributable also to physical and chemical factors as well as grazing, but more importantly demonstrated that such data can be adequately extracted from freely available satellite archives. Further research could extend this approach to cover larger areas to obtain holistic estimates for lakes, and, if coupled with additional data on fish and macroinvertebrates, it could significantly advance knowledge on structural and functional relationships essential for conservation of such internationally important sites.

Conclusions
The dynamic changes in the macrophyte community in the Sirmione Peninsula located in the south of Lake Garda indicate that this important ecosystem is delicately balanced and dependent on many factors. The changes in macrophyte total cover appear to be driven by a complex of wave erosion with eutrophication, turbidity, and lake level influencing the light climate, as well as grazing by herbivores. Satellite observation reveals the detail of seasonal and interannual variation in macrophyte coverage and density crucial to informing ecological assessments and understanding the functioning of dynamic, highvalue ecosystems. Further investigations involving the functional features of macrophytes (e.g., growth forms, life history strategies (e.g., Temmink et al., 2021 [102])) should be the next step needed to improve current evidence and ensure more robust insights to address the effects of global change on such fundamental ecosystems.
Specifically, S2 allowed us to accurately map the spatiotemporal dynamics of submerged rooted macrophytes in order to support monitoring of the shallow water ecosystem of the Sirmione Peninsula. Our results, in addition to capturing a marked fluctuation of macrophyte spatial patterns during the analyzed period (2015-2020), allowed us to highlight the gradual decline in the density of these macrophytes compared to the situation in 1997 and the shift of these organisms towards deeper areas. In addition, comparing the spring period of 2020 (lockdown) with the other years showed a clear improvement of the vegetation community, especially in the nearshore areas where ports are located, thus demonstrating the crucial role of boats in the growth dynamics of macrophytes.
In conclusion, the use of remote sensing has offered a new perspective regarding the comprehension of lacustrine macrophytes dynamics, allowing us to overcome the technical and economic limitations of traditional limnological techniques. In the near future, it could be interesting to expand the study, both exploiting the characteristics of the new hyperspectral satellite sensors (e.g., PRISMA and DESIS) and enriching ancillary data (e.g., lake water levels from the CCI Lake ECV project), for further improving the analysis and monitoring of macrophytes in the region.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/app12052693/s1, Figure S1: Sentinel-2 image of Sirmione Peninsula in true color (23/03/2019). The circled images represent the three invariant surfaces: roofs (red), sandy beach (yellow), and evergreen vegetation (dark green). The light green area located in the west sub-basin (35 pixels) and the brown area in the east sub-basin (90 pixel) represent the two ROIs used for the analysis of the error maps. Figure S2: (a) Mean annual error value for the two ROIs (2015-2020); (b) mean seasonal error value considering all years of the dataset (from spring to autumn). Figure S3: Winter wind rose for southern Lake Garda showing frequency (%), strength (m s −1 ), and direction of wind (2015-2020). DJF = December, January, February. Figure S4: Turbidity estimated from Lakes CCI dataset for the eastern and western portion of the Sirmione Peninsula for the period from 2016 to 2020. NTU = Nephelometric turbidity unit. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.