Monitoring Lakes Surface Water Velocity with SAR: A Feasibility Study on Lake Garda, Italy

Abstract

The SAR Doppler frequencies are directly related to the motion of the scatterers in the illuminated area and have already been used in marine applications to monitor moving water surfaces. Here we investigate the possibility of retrieving surface water velocity from SAR Doppler analysis in medium-size lakes. ENVISAT images of the test site (Lake Garda) are processed and the Doppler Centroid Anomaly technique is adopted. The resulting surface velocity maps are compared with the outputs of a hydrodynamic model specifically validated for the case study. Thermal images from MODIS Terra are used in support of the modeling results. The surface velocity retrieved from SAR is found to overestimate the numerical results and the existence of a bias is investigated. In marine applications, such bias is traditionally removed through Geophysical Model Functions (GMFs) by ascribing it to a fully developed wind waves spectrum. We found that such an assumption is not supported in our case study, due to the small-scale variations of topography and wind. The role of wind intensity and duration on the results from SAR is evaluated, and the inclusion of lake bathymetry and the SAR backscatter gradient is recommended for the future development of GMFs suitable for lake environments.

1. Introduction

Microwave methodologies based on the use of Synthetic Aperture Radar (SAR) sensors have been widely developed for many Earth surface monitoring applications. In the last decade, images acquired by SAR have been increasingly exploited by the scientific community for deformation surveys [1] even at single facilities scale, that is, buildings and infrastructures [2]. More recently, increasing attention has been paid to the observation of sea surface, complementing the traditional use of optical or multispectral images.

Differently from optical sensors, SAR allows the obtaining of all weather and day–night 2D images of the illuminated area of the Earth’s surface [3]. SAR remote sensing for ocean, seas and coastal applications mostly exploits the amplitude of the backscattered signal for, for example, monitoring oil-spills [4] and sea-ice [5], ship detection [6] and high-resolution wind fields retrieval [7]. However, SAR uses coherent radiation and the complementary information carried in the phase of the received complex signal can be also exploited. An example from land applications is the use of the complex information from a single SAR image to reconstruct infrastructures micro-motion [8]. By analysing the complex backscattered (received) signal, it is possible to measure the Doppler properties of the scatterers. The latter are directly related to the motion of the scatterers along the radar line of sight (LOS) in both land [9] and marine [10] contexts. In ocean applications, these properties are currently used to retrieve near surface wind speed [11] and surface current [12].

Several factors contribute to the measured Doppler frequency, such as the bulk surface current, the wind waves and the wave–current interactions at different scales [13]. The drift current induced by the wind on the sea surface (wind drift) and the sea state (more precisely the high-frequency waves) depend on the near-surface wind field. In this regard, [14] performed a global ocean analysis and and found that an increase of the wind speed (projected along the LOS) corresponds to a first-order increase of the Doppler anomaly. Such a wind-wave bias is empirically estimated and removed to retrieve the radial component of the bulk surface current. To this end, empirical Geophysical Model Functions (GMF) have been determined [12].

This analysis is going toward a consolidated state for oceanic applications, such that surface velocity maps in open waters are distributed by the European Space Agency (ESA) among the second level products of Sentinel 1 [15]. On the contrary, the extension of a similar approach to coastal zones is more delicate. At least two main assumptions of the analysis of the Doppler Centroid Anomaly (DCA) in the open ocean do not hold for coastal areas: (i) uniform surface flow field and (ii) fully developed wind-waves. In the near-shore areas, small- scale variations of topography and wind can cause higher variability in the currents and the waves to develop in fetch-limited condition.Thus, the application to other contexts of the above mentioned methods for the removal of the wind waves bias may not be straightforward as in the open ocean case, although examples along this line are given in the recent literature [10,13,14,16,17,18,19]. At the smaller spatial scale of enclosed basins (e.g., lakes), no studies have yet explored the possibility of inferring surface currents from SAR Doppler frequencies.

More generally, SAR images of lakes have been seldom exploited for several reasons, among which the minor intensity of hydrodynamic processes (e.g., smaller wave heights and currents than in open seas and oceans) and the effect of surrounding topography (e.g., foreshortening and layover effect [20]).

Currently, SAR applications for lakes are limited to the analysis of the backscattered signal amplitude. SAR amplitude images have been used in support of optical or thermal imagery for identifying and mapping inland waters surface features, for example, ice [21] and cyanobacteria blooms [22,23]. A few tests also showed the potentialities of using SAR images for retrieving the size of surface eddies in large lakes [24] and spatially distributed information on wind intensity [25].

The SAR Doppler analysis could represent a useful way to retrieve surface velocity in lakes and would compensate for the lack of synoptic measurements of the key physical quantities determining lake hydrodynamics. Indeed, while the use of remote sensing products is well consolidated to infer lake surface water temperature (LSWT) [26], water velocity, as well as all wind velocity, is traditionally monitored with in-situ point measurements. In most cases, traditional instrumentation requires direct access to the lake, with strong operational limitations, and provide data at single locations. Remote sensing techniques provide instead synoptic coverage, fine spatial detail and repeated regular sampling.

Periodic snapshots of lake surface currents, at the present day, cannot be achieved alternatively than with remote sensing instruments, for example, airborne or space, the latter being traditionally exploited in marine applications [27,28,29]. Few attempts towards the reconstruction of the surface transport from thermal infrared imagery [30] showed encouraging results in this direction, but such approach requires high temporal resolution, which is not always guaranteed. Moreover, the use of thermal infrared remote sensing is hindered by cloud cover. In this regard, the use of SAR images combined with the DCA technique could overcome the technical issues typically encountered for reconstructing the surface flow field in lakes, asides from being rather innovative.

This study aims to investigate the feasibility of extracting the surface velocity from DCA in lakes by considering as test case Lake Garda, a large and deep lake in northern Italy. The test case is a clear oligotrophic lake, thus allowing for neglecting turbidity and algal blooms which might affect the SAR signal.

We analyze and process images acquired by the ENVISAT C-Band sensor. Results are compared with the outputs of a three-dimensional (3D) numerical model validated for Lake Garda [31]. Water surface temperature products from the MODIS (Moderate Resolution Imaging Spectroradiometer) Terra sensor are also considered to further validate the simulation results. The numerical results are used to interpret the SAR retrieved signal and to evaluate the different factors influencing the Doppler shift analysis in the case study. Based on this, we draw some guidelines on the variables to be considered for a future definition of GMFs for lakes.

2. Materials and Methods

In this section we describe the test site (Section 2.1), the sensors and datasets for SAR and thermal infrared images (Section 2.2), the procedure for the analysis of the DCA (Section 2.3) and the adopted numerical models (Section 2.4).

2.1. Case Study

Lake Garda (Figure 1a) is one of the large European perialpine lakes and is the largest lake in Italy by surface extension and volume. It covers an area of 370 km${}^2$, with a water volume of 50 km${}^3$ and a maximum depth of 346 m. The lake represents an essential water supply for many sectors of the local economy (e.g., agriculture, industry, fishing and drinking [32]). Thanks to its attractive landscape, mild climate and water quality, Lake Garda is also an important resource for recreation and tourism.

The lake has a heterogeneous bathymetry (Figure 1a). A narrow (average width 4 km) and deep (maximum depth 346 m) northern trunk, enclosed between steep mountains, is connected to a southern large (maximum width 18 km) and shallower (average depth 65 m) basin, which is characterized by more gentle slopes and lies in a flat plain.

The typical winds in Lake Garda are thermal breezes [33]. In the northern region, the most frequent provenance direction of high speed winds is along the longitudinal axis of the lake, where winds are channeled by the steep lateral mountains (see the orography in Figure 1a). During summer, northerly and southerly breezes alternate in the night/morning and in the afternoon/evening, respectively [34]. Such breezes are characterized by moderate to high velocities (≈5–10 m/s), alternating directions and not necessarily uniform spatial distribution. During winter and mid-seasons, winds are usually weaker. However, it is not uncommon for this lake to be subject to long-lasting storm winds at a synoptic scale, for example, Föhn winds. These winds frequently come from North-East, have almost uniform spatial distribution on the lake surface, reach speeds of more than 10 m/s and often last for more than one day [35].

2.2. Sensor and Dataset

In the last years, there has been a consistent proliferation of SAR systems with different technical characteristics, that is, resolution, operational frequencies and revisiting time [36]. Among them we selected the ENVISAT, following the experience gained in marine environment with C-Band sensors for the estimation of the surface current velocity. In this work, we adopt a methodology validated in the coastal area of the Gulf of Naples [16,27,37] and the Gulf of Trieste [38]. ASAR ENVISAT operated in C Band (wavelength $\lambda$ equal to 5.6 cm), in multiple polarization mode (HH/VV,HH/HV,VV/VH). In stripmap mode it covered up to 100 km with a single swath. The spatial resolution of the sensor’s Single Look Complex (SLC) product is 5 m in azimuth and 20 m ground range (see Figure 1b for SAR acquisition geometry). A considerable archive of data exists from the ENVISAT operative period, that is, from 2002 to 2012, on ascending and descending orbits. Within the broad landscape of SAR sensors currently available, data acquired by Sentinel 1 (S1) could also be exploited. S1 operates in C-Band and regularly provides very wide swaths of about 250 km by using the TOPSAR (Terrain Observation with Progressive Scans SAR) operation mode. However, due to the steering of the antenna along the azimuth, this operation mode complicates the estimation of the Doppler Centroid. The presence of discontinuities in the estimated Doppler Anomaly, which in turn affects the retrieved radial velocity, has been already observed [39]. Based on these considerations and on the experience we gained in marine environments, ENVISAT data were chosen for this exploratory study on Lake Garda. The ENVISAT archive was screened to eliminate images characterized by very low backscattering coefficients. We selected six images showing amplitude patterns, that is, sufficiently intense and spatially varying backscattering, potentially associated with interesting hydrodynamic features and fulfilling the wind criteria specified in detail in Section 3.2. The selected images are acquired over ascending orbit, corresponding to evening acquisitions (around 21:00 UTC), and span from 2008 to 2010 in summer, autumn and winter seasons. The choice of ascending acquisitions is due to the favourable acquisition geometry related to the orientation of the lake. In fact, the ground projection of the ascending sensor LOS is more aligned with the longitudinal axis of the lake than the descending one. We note that the longitudinal axis of the lake is expected to be the main direction of lake currents and waves propagation, as it also corresponds with the typical direction of winds (see Section 2.1). The images have been classified based on the wind velocity registered at the acquisition time. In

Table 1. Mean, standard deviation and maximum of the simulated key geophysical quantities for each date considered. Module of wind speed from WRF ($W^{\mathrm{SIM}}$); module of the surface water velocity from Delft3D + SWAN ($v^{\mathrm{SIM}}$); significant wave height ($H_{\mathrm{s}}$) and peak wavelength (${\lambda }_{\mathrm{peak}}$) from SWAN.

Class	Date	${\mathit{W}}^{\mathbf{SIM}}$ [m/s]		${\mathit{v}}^{\mathbf{SIM}}$[cm/s]		${\mathit{H}}_{\mathit{s}}$ [cm]		${\mathit{\lambda }}_{\mathbf{peak}}$ [m]
		Mean ± Std	Max	Mean ± Std	Max	Mean ± Std	Max	Mean ± Std	Max
low	29 October 2009 ${}^{*}$	2.3 ± 1.5	5.5	8.3 ± 6.8	30.3	8.6 ± 5.6	19.6	2.04 ± 1.23	4.40
low	5 August 2010	3.1 ± 0.9	5.1	10.6 ± 5.4	33.5	11.4 ± 4.2	17.4	2.34 ± 0.94	3.74
low	14 October 2010 ${}^{*}$	1.3 ± 0.5	2.1	3.2 ± 2.0	13.6	3.3 ± 0.9	5.3	0.81 ± 0.18	1.14
high	7 February 2008 ${}^{*}$	8.2 ± 3.2	16.4	16.4 ± 8.1	61.0	49.7 ± 24.4	92.4	9.71 ± 5.04	19.24
high	13 November 2008	10.2 ± 5.2	22.2	33.1 ± 21.2	117.5	68.1 ± 42.7	145.1	12.79 ± 8.31	27.8
high	11 June 2009 ${}^{*}$	7.9 ± 4.6	17.2	17.3 ± 8.4	53.2	48.9 ± 35.0	109.7	9.34 ± 7.00	29.67

${}^{*}$ MODIS product available.

the selected dates are listed and their simulated geophysical characteristics are summarized. For a more detailed description of the investigated dates, we refer to Section 3.2.

For the selected dates, thermal images from MODIS were also analyzed to validate the simulations. In particular, we preferred nighttime Level 2 MODIS Terra Sea Surface Temperature (SST) products, because they are closer in time to SAR acquisitions. MODIS images are cloud free for four out of six investigated dates. The four dates are indicated with an asterisk in Table 1.

2.3. SAR Doppler Centroid Anomaly Estimation

For the retrieval of information about the ocean surface velocity from SAR data, two methodologies have been developed: the Along-Track Interferometry (ATI) [18] and the DCA method [14]. Here we exploit the technique of the DCA, which was successfully applied in open oceans [14,40,41] and ocean coastal areas [14,18,42]. The same approach was also adopted for the estimation of water currents in the Mediterranean sea, a semi-enclosed sea, where they show a smaller range of intensity, [16,27,38]. The DCA technique estimates the surface velocity of the illuminated area by measuring the variation that such a velocity produces in the spectrum of the SAR image [14]. The principle is rather well known: let us consider a SAR system acquiring data on a scene where a target is in relative motion with respect to the sensor. When the transmitted radiation hits the target, it is backscattered toward the radar along the antenna pointing direction, also called Line of Sight (LOS) or radial direction (see Figure 1b). If the velocity associated with the motion has a non-zero component along the radial direction ($v_{\mathrm{r}}$), this component determines a shift in the azimuth spectrum of the received signal. Such a shift is called Doppler Centroid frequency $f_{DC}$ and is proportional to the radial component of the velocity associated to the relative motion between target and sensor.

Some remarks on the formation of the Doppler shift, which is influenced by the target and sensor motion along the LOS direction, are in order. First of all, the radial direction is influenced by the platform attitude as well as by possible electronic steering introduced to point the beam in a direction different from the broadside one (squinted acquisition). Secondly, the measurement of $f_{DC}$ from the power spectral density allows the estimation of the radial component of the (relative) target motion [43,44]. For each pixel, say with azimuth x and slant range r, $f_{DC}(x,r)$ can be considered the sum of two contributions: the Doppler Centroid $f_{DC0}(x,r)$ and Doppler Centroid Anomaly $f_{DCA}(x,r)$. The first one ($f_{DC0}$) corresponds to a stationary scene and is referred to as “background DC”, due to a scene in which all the scatterers move at the same velocity (background motion). The second one ($f_{DCA}$) is the Doppler component associated with the target motion with respect to the background one.

From what is stated above, it follows that for the satellite case the background DC, $f_{DC0}(x,r)$, is contributed by both the platform movement and the Earth rotation, in a way that depends on the platform attitude and beam orientation. This term is generally a slow varying function along the space, typically well approximated by a polynomial in r. It is worth noting that, strictly speaking, ionospheric variations could also impact $f_{DC0}(x,r)$: this phenomenon is, however, more relevant for low frequency systems, for example, L- and P-Band SAR sensors [45]. The latter are usually not exploited for water surface velocity retrieval, due to the reduced sensitivity associated with the increase of the wavelength.

As for ASAR-Envisat, the term $f_{DC0}(x,r)$ can be retrieved from the ancillary data through annotated polynomial coefficients, which are computed knowing the platform attitude, the antenna pointing and the scene latitude. Previous experiences carried out with other platforms showed that this information is frequently not sufficiently accurate, especially for the range variation. This fact might be due to inaccurate attitude information, to unexpected beam bending as well as to other spurious effects. Following the strategy adopted in [16] for marine applications, we estimated $f_{DC0}(x,r)$ directly from the data via a polynomial Least Square fitting along r. More details on the estimation of the background Doppler component associated with the stationary scene are given in [16]. In the same work it is also shown how to handle unwanted azimuth oscillating components, frequently observed for ASAR-ENVISAT acquisitions and evidently not associated with any geophysical process.

The velocity at which the target proceeds toward or away from the radar, with respect to the background scene motion, can be thus obtained by removing from the measured Doppler Centroid its component associated to the stationary scene, as follows:

$$\begin{array}{cc}\hfill v_{\mathrm{r}}(x,r)& =\frac{\lambda }{2}[f_{DC}(x,r)-f_{DC0}(x,r)]\hfill \\ \hfill & =\frac{\lambda }{2}{f}_{DCA}(x,r),\hfill \end{array}$$

(1)

where $\lambda$ is the sensor wavelength. The scatters’ surface velocity along the ground range direction $v_{\mathrm{gr}}$ can be finally obtained from $v_{\mathrm{r}}$ by considering that the latter is its projection along the LOS (see also Figure 1b):

$$v_{\mathrm{gr}}=\frac{v_{\mathrm{r}}}{sin\theta },$$

(2)

where $\theta$ is the incidence angle.

In order to retrieve the actual lake ground range surface velocity, any possible spurious effect, such as the wind-wave bias [12], must be precisely identified, estimated and removed from $v_{\mathrm{gr}}$. Thus, $v_{\mathrm{gr}}$ represents a raw variable to be processed with an appropriate GMF to obtain the ground range surface velocity, hereinafter surface velocity. The identification of the spurious effects on $v_{\mathrm{gr}}$ is the goal of this study, where no GMFs are used, and it is performed with the help of the numerical model.

In Figure 2a, we summarize the procedure for the generation of $v_{\mathrm{gr}}$ maps. After the SAR focusing operation on the raw data, the measurement of $f_{DC}$ is performed pixel by pixel by estimating the azimuth self-correlation function along the azimuth direction [43]. The Doppler spectrum is evaluated on patches, hereafter referred to as Doppler Resolution Cell (DRC), sliding along the whole SAR image. The DRC has been chosen to guarantee a good balance between spatial and spectral resolutions. In particular, we considered patches of width 512 azimuth by 128 range pixels, corresponding to a DRC of 2.5 km × 2.5 km. In this condition the spectral resolution is 2.85 Hz, which corresponds to a velocity resolution of 7.98 cm/s for $v_{\mathrm{r}}$.The Doppler Centroid estimation accuracy can be estimated through the calibration of the radial velocity over the land area (i.e masking out Lake Garda). For all the acquisitions here analyzed (see Table 1) $v_{\mathrm{r}}$ has zero mean and a standard deviation ranging from 27 cm/s to 30 cm/s. Note that these values account for the variability of the residual signal all over the scene with the lake only covering a small portion.

Subsequently, the compensation of the stationary component of the Centroid is performed and the final DCA is estimated (Equation (1) and then $v_{\mathrm{gr}}$ from Equation (2)). The DCA map is converted to surface velocity and geocoded through a SRTM DEM on a geographical grid with a spacing of about 92.7 m N and 64.9 m E [46]. At this stage, the geocoded surface velocity map is comparable with geo-referred simulations results and MODIS products.

The output product is supplemented by the basckscattering information. A standard calibration procedure is applied to the focused amplitude data to get the normalised measure of the radar cross-section per unit area on the ground, that is, the backscatter coefficient ${\sigma }^0$. A multi-look algorithm is then applied to reduce the noise (or “speckle”) of SAR images by averaging adjacent pixels. The resulting ${\sigma }^0$ with resolution 100 m × 100 m is finally geo-referred (Figure 1c).

2.4. Numerical Modeling

For simulating the flow field, we use a modeling chain (Figure 2b) composed by three numerical models: an atmospheric, a hydrodynamic and a wave model. In this work, we start from the outputs of an existing model for the case study [35,47] in the configuration presented in [31], which consists of the hydrodynamic model Delft3D [48] fed by the results of the Weather Research and Forecasting (WRF) atmospheric model [48].

The WRF model provides the atmospheric variables as time and space varying fields of wind velocity, air pressure, temperature and relative humidity, incoming shortwave radiation and cloudiness at the lake surface. The resolution of the simulated weather forcing is 2 km in space and 1 h in time.

From the WRF outputs, the hydrodynamic model Delft3D computes the heat fluxes between air and water and wind stress at the lake surface, which are necessary for simulating the flow field, turbulence, heat and mass transport within the lake. The nominal horizontal resolution of the hydrodynamic model is ≈100 m, with smaller cells in the northern and larger cells in the southern part. The vertical layers are 1 m thick near the surface (over the upper 40 m) and smoothly increase their thickness up to 25 m at the deepest point. The outputs of the described model are available over the period 2004–2018 and will be henceforth referred to as “long-run data”. The long-run data are saved on a daily basis at 10:00 UTC.

For the present work, the available long-run data were used to restart ad hoc simulations of the selected dates (Table 1) with the same setup as in [31]. The restart was necessary to produce outputs synchronous with the ENVISAT acquisitions (21:00 UTC) and to couple the Delft3D model with the wave model. We adopted a third-generation wind-wave model known as Simulating Waves Nearshore (SWAN) [49]. The wave spectrum employed in the simulations was composed of 101 logarithmically spaced frequencies in the range of 0.15–3 Hz. The frequency range was specifically chosen to be representative of wave field expected in Lake Garda and therefore included higher frequencies than those typically employed in coastal applications, but neglected very low frequencies that hardly occur in medium-size lakes. Since no data were available for the calibration of the model, the parameters were chosen in consistency with previous experiences in other large and deep lakes (e.g., Lake Constance [50] and Michigan [51]).

The results of the six Delft3D + SWAN simulations were saved at hourly intervals, over the whole computational domain, and covered a period of two days (the SAR image acquisition day plus one day before).

3. Results

In this work we use the model outputs with the aim of interpreting the results obtained from the SAR images. Thus, we rely here on the assumption that the model is capable of reproducing wind fields and surface transport. The reliability of the wind fields simulated by the WRF model has been extensively demonstrated by [31,52]. As for the surface transport, we base our work on [31], who showed that the Delft3D model results are coherent with in-situ and remotely sensed water temperature observations. We stress that [31] provided an indirect verification of the flow field, since in-situ measurements of water velocity are not available. In Section 3.1 we perform a similar verification by comparing the spatial patterns of surface temperature simulated by the model and reconstructed from the MODIS images. In Section 3.2, we introduce the main features of the investigated dates and in Section 3.3 and Section 3.4, we evaluate the outputs of our procedure for the extraction of lake surface currents from SAR.

3.1. Model Verification Against MODIS

We provide here the comparison between the Level 2 MODIS Terra SST products ($T^{\mathrm{MODIS}}$) in Celsius degrees and the simulated lake surface water temperature ($T^{\mathrm{SIM}}$). Figure 3 shows that the model correctly reproduces the remotely sensed patterns of lake surface temperature and the spatial variability of this quantity in the different seasons. This confirms that the heat fluxes between atmosphere and lake are well simulated, and that the numerical flow field is coherent with the transport patterns responsible for horizontal anomalies in lake surface temperature. The difference between the spatial mean values of $T^{\mathrm{SIM}}$ and $T^{\mathrm{MODIS}}$ (reported in each panel of Figure 3) can be due to model errors as well as to the skin-bulk bias [53]. However, in all investigated dates the mean error is below 1 ${}^{\circ }$C, well in line with the performances of the long-run model in [31].

3.2. Main Features of Investigated Dates

The wind field has been recognized as playing a major role for the use of SAR in the open ocean [54]. In particular, in order generate capillary waves detectable from C-band sensors, the wind forces have to overcome the viscous ones. A wind speed threshold value for this process to occur is estimated to be about 3.25 m/s (at 10 m above the surface, [55]).

Therefore, we classify the six cases considered in this study (and listed in Table 1) based on the wind intensity. In Figure 4 the daily cycle of wind speed and direction simulated in a mid-lake point (point P in Figure 1) is displayed for the selected dates. The SAR acquisition time is indicated as a red vertical line in the figure panels.

We define “low wind” dates (Figure 4a–c) the three acquisitions in which the spatial mean of the simulated wind speed is equal or smaller than this threshold value (29 October 2009, 5 August 2010, 14 October 2010). In these dates mild evening breezes blew at 21:00 UTC over the lake surface (≈1–3 m/s). This condition is frequent in sunny and warm days, when thermal breezes develop with alternating directions during daylight hours and calm down after the sunset. Under low wind conditions, simulations predict mean surface velocities below 10 cm/s, significant wave heights $H_{\mathrm{s}}$ of the order of few cm, and peak wavelengths ${\lambda }_{\mathrm{peak}}$ of the order of 1 m (Table 1). Such values suggest that gravity effects do not completely dominate the dynamics as surface tension still plays a significant role [56].

The remaining three acquisitions are classified as “high wind” dates (7 February 2008, 13 November 2008, 11 June 2009, Figure 4d–f). On those days, the mean value of wind speed is above 7 m/s and its maximum value exceeds 15 m/s in large regions of the lake. Both the mean and the maximum wind speed values are significantly higher than the above-mentioned threshold. Such high speeds are often due to large scale synoptic winds, blowing almost uniformly over the lake’s surface for more than 12 h. Figure 4d,e shows that the latter is indeed the case for the dates 7 February 2008 and 13 November 2008, when a northerly Föhn wind persistently blew for one day before the SAR acquisition time. On the contrary, the date 11 June 2009 (Figure 4f) coincides with the onset of a similar wind, which started one hour before the acquisition time and lasted all night long (not shown).

The simulated gravity waves on high wind speed days have $H_{\mathrm{s}}$ of the order of 50 cm on average (Table 1), overcoming 1 m in some areas of the lake, and ${\lambda }_{\mathrm{peak}}$ of the order of 10 m. These values are fully consistent with the gravity waves field that can be theoretically expected [57] and observed [58] in a medium-size lake.

In the next two sections, we will provide for all investigated dates the amplitude SAR backscatter image (${\sigma }_0$), the ground range surface velocity retrieved from SAR DCA (Figure 2a, $v_{\mathrm{gr}}$) and simulated by the numerical model (Figure 2b, $v_{\mathrm{gr}}^{\mathrm{SIM}}$), and the simulated wind field (direction and ground range component magnitude $W_{\mathrm{gr}}^{\mathrm{SIM}}$). The represented quantities are obtained as follows:

• ${\sigma }_0$ is processed as described at the end of Section 2.3 expressed in dB to enhance the visibility of its variation over lake;
• $v_{\mathrm{gr}}$ is obtained as in Equation (2).
• $v_{\mathrm{gr}}^{\mathrm{SIM}}$ is obtained by computing the component of the simulated surface velocity along the ground range;
• $W_{\mathrm{gr}}^{\mathrm{SIM}}$ is obtained by computing the component of the simulated wind velocity vector (at 10 m a.w.s.) along the ground range direction.

3.3. Low Wind Dates

The results for the “low wind” dates are reported in Figure 5. In marine/ocean observations, a low-wind regime corresponds to a low mean amplitude and a limited spatial variation of the backscatter. In our case study, Figure 5 clearly shows a different behaviour. The SAR backscatter retrieved on these dates (first column of the Figure 5, panels a, e, i) shows recurring patterns in all examined images. Higher values are detected in the northwestern basin, lower values in the south-eastern region. A high intensity feature located at the NE edge of the lake is present in all the ascending SAR images analysed. It consists in a cross-track radiometric compression, known as foreshortening effect, which occurs when the radar beam impinges on the foreslope of a mountainous area [20]. The component of surface motion retrieved from SAR DCA $v_{\mathrm{gr}}$ is represented in the second column of Figure 5, panels b, f, j. Red and blue colors (redshift and blueshift) indicate a motion towards and away from the radar antenna, respectively. We recall here that the surface velocity has been estimated with a DRC of 2.5 km × 2.5 km (see Section 2.3) and that any variation of the resulting $v_{\mathrm{gr}}$ at a smaller spatial scale has to be interpreted with care. While this is not a critical issue along the longitudinal axis of the lake and in the southeastern basin, where both the length and the width of the lake are significantly larger than the DRC size, difficulties arise in the northern trunk of the lake. Here the DRC size is comparable to the width of the lake and the estimation might be affected by the presence of portions of land inside the DRC. For such a reason we recommend caution in the interpretation of the variations of $v_{\mathrm{gr}}$ along the crosswise direction.

The spatial patterns of $v_{\mathrm{gr}}$ resemble those from ${\sigma }^0$. In the northern trunk of the lake, $v_{\mathrm{gr}}$ shows in all dates a recurring pattern: it overcomes 2 m/s in the central and western region, and is almost null along the eastern shore. This behaviour suggests a south-westward surface transport which is not confirmed by the numerical model (third column of Figure 5, panels c, g, k). In fact, the simulated field of $v_{\mathrm{gr}}^{\mathrm{SIM}}$ shows significant differences from one date to another, mostly due to the differences in the spatial distribution of the wind field (fourth column of Figure 5, panels d, h, l). Moreover, the simulations show significantly lower surface velocities, below 0.1 m/s, as it is reasonably expected in a lake of the dimensions of Lake Garda under low wind conditions [59].

In the southern sub-basin, similar patterns are found in $v_{\mathrm{gr}}$, retrieved on 29 October 2009 and 14 October 2010 (Figure 5b and j respectively), with high and positive velocities at the center of the sub-basin and almost null or negative velocities close to the shores. The spatial distribution of $v_{\mathrm{gr}}$ in this region is similar to that of $v_{\mathrm{gr}}^{\mathrm{SIM}}$ on the same dates (Figure 5c,k). On 5 August 2010 $v_{\mathrm{gr}}$ (Figure 5f) is strongly positive over the entire southern region, while $v_{\mathrm{gr}}^{\mathrm{SIM}}$ (Figure 5g) is strongly negative. The spatial distribution of $v_{\mathrm{gr}}^{\mathrm{SIM}}$ on this date reveals a predominant eastward surface transport in the southern basin and an opposite tendency in the northern basin, which is the direct consequence of two divergent wind fields with respect to the ground range direction (Figure 5h).

3.4. High Wind Dates

The results for the “high wind” dates are reported in Figure 6. The SAR backscatter images (Figure 6, panels a,e,i) and the surface velocity $v_{\mathrm{gr}}$ (Figure 6b,f,j) present similar patterns, which are in good agreement with the simulated surface velocity $v_{\mathrm{gr}}^{\mathrm{SIM}}$ (c,g,k), especially in the elongated region.

On all dates ${\sigma }^0$ shows a pattern characterized by increasingly high values as moving south-west in the northwestern basin (Figure 6a,e,i). Such pattern can be observed also in the $v_{\mathrm{gr}}$ images (Figure 6b,f,j), with a redshift located mostly along the western shore indicating a south-westerly surface current moving towards the sensor. The value of $v_{\mathrm{gr}}$ in this area is higher than 2 m/s and slightly decreases in the southeastern part of the region. This behavior finds qualitative correspondence with the patterns of $v_{\mathrm{gr}}^{\mathrm{SIM}}$ (Figure 6c,g,k), where the signature of wind waves can be detected from localized peaks of velocity up to 0.5 m/s. The black arrows displayed in Figure 6d,h,l indicate the presence of a wind blowing from north-east to south-west, which is responsible for the surface transport observed both from SAR and the model.

4. Discussion

In all cases examined, the surface ground range velocity retrieved from SAR ($v_{\mathrm{gr}}$) is much higher than the simulated $v_{\mathrm{gr}}^{\mathrm{SIM}}$. Regardless of the wind intensity, $v_{\mathrm{gr}}$ always overcomes ±2 m/s, which is too high for a medium-size lake environment. Despite the empirical relations linking lake surface velocity and wind speed are often site-specific, the order of magnitude of the first quantity is generally assumed to depend on the second. A wind factor can be computed as the ratio between water surface velocity and wind speed. In open waters the wind factor is typically of the order of 3%, but it can also be lower for the nearshore areas [60]. In the investigated dates, the simulated water velocity was found to be, on average, between 2% and 3.5% of the wind speed. Even assuming that the observed $v_{\mathrm{gr}}$ is the only component of the surface current (i.e., assuming that the water is moving precisely along the ground range direction), a value of 2 m/s would require a wind speed of more than 60 m/s, which is an absolutely unrealistic wind speed. Such wind speed would be even higher if we take into account also the other component of surface motion, which is likely to exist and contributes to the module of surface current, despite it can not be retrieved with the DCA technique. Hence, the velocity retrieved from SAR is a large overestimation of the actual lake surface velocity.

Such overestimation is also found in open ocean applications [61], where $v_{\mathrm{gr}}$ is affected by the surface bulk velocity and by a term of the order of the 30% the ground range component of the wind velocity $W_{\mathrm{gr}}^{\mathrm{SIM}}$ [14]. The latter term represents a bias to be removed for the effective estimate of the surface ground range velocity. Such bias is typically predicted and removed using an ad hoc geophysical model function (GMF), for example, CDOP [12], which exploits the correlation between DCA and the ground range component of the wind velocity $W_{\mathrm{gr}}^{\mathrm{SIM}}$. The key geophysical process underlying this function is represented by wind waves, due to their orbital velocities.

Following this approach, the existence of a correlation between $v_{\mathrm{gr}}$ and $W_{\mathrm{gr}}^{\mathrm{SIM}}$ is investigated. Scatter plots of $v_{\mathrm{gr}}$ vs $W_{\mathrm{gr}}^{\mathrm{SIM}}$ are shown in the first row of Figure 7 (panels a, b from the low and high wind dates respectively). The colour scale reflects the probability density estimation from data computed via a Kernel Density Estimation (KDE) algorithm. The two variables appear not to be correlated in all cases examined here. This was reasonably expected, as in small to medium-size lakes the wind–waves field is typically underdeveloped, with complex dynamics often subject to local effects (such as topography) [57]. Thus, the validated approaches operationally used in the open ocean [10,15] cannot be applied to our case study and other quantities possibly affect $v_{\mathrm{gr}}$.

Patterns observed in Figure 5 and Figure 6 for ${\sigma }^0$ and $v_{\mathrm{gr}}$ suggest a possible contribution from bathymetry to the SAR products. Figure 7c,d shows that a correlation between $v_{\mathrm{gr}}$ and bathymetry exists in all the investigated dates: higher correlation between $v_{\mathrm{gr}}$ and water depth is found the shallower areas (0 to 100 m deep), and in the deepest area (below 300 m). It is recognized that underwater bottom topography can be sensed indirectly via surface effects [62], although the electromagnetic waves emitted by the radar penetrate water only to a depth that is small in comparison to the radar wavelength. In fact, the underwater bottom topography generates short-scale surface roughness variations, which in turn modulate the radar reflectivity. As a complement to the analysis of the correlation between $v_{\mathrm{gr}}$ and water depth, we note that the deepest region of the lake (depth > 300 m, where highly correlated values are visible in Figure 7c,d) coincide with the northern narrow area of the lake. In this area, the width of the lake ranges between 3 and 5 km, that is, between one and two times the DRC size (2.5 km). Hence, we assume that an effect of the mixed signals from water and land shall also be considered. To what extent such an effect combines with the bathymetry and the wind waves field and alters the measured $v_{\mathrm{gr}}$ is beyond the aims of this work and we recommend further investigation in this sense.

Given the specific geometry of Lake Garda, the observed significance of bathymetry could, in principle, also be ascribed to the presence of longer (meters scale) wind waves. As shown in Figure 7c,d there is a peak of correlated values around ≈170 m depth and $v_{\mathrm{gr}}\approx 1.5$ m/s. Such a feature, in the case of high wind, is likely related to the wave field development. In this regard, in Figure 8a–c we report the peak wavelength map in the high wind dates. As it is clear from the first two days (panels a,b), ${\lambda }_{\mathrm{peak}}$ reaches its maximum at the end of the northern trunk, where water depth is below the two 150 m isobaths. This area also corresponds to maximum effective fetch for a wind blowing from North-East. The relative weight of the short and long waves also affects $v_{\mathrm{gr}}$ and is at the base of the applications of the DCA method to oceans [11,61]. Nevertheless, it cannot be easily quantified in lake applications, despite the fact that it can be inferred from ${\sigma }^0$ variations. According to the Bragg scattering model, the radar backscatter is proportional to the surface wave spectral density at the Bragg wavelength. As these centimeters-scale waves are tilted by longer (meters-scale) waves, ${\sigma }^0$ varies along with such longer wave profiles. This leads to a correlation between ${\sigma }^0$ and the orbital velocities of the peak waves [63]. Theoretical findings predict an average increase in the mean spectral density of the short waves due to the interaction between short “Bragg” waves and meters-scale waves [64]. Such a relationship also holds in the high wind cases analyzed here, as shown in Figure 7f, where a clear correlation exists between ${\sigma }^0$ and ${\lambda }_{\mathrm{peak}}$. However, no correlation is found in the case of low wind dates (Figure 7e), when the wave field is not developed.

Figure 7g,h clearly shows that a tight relationship exists between $v_{\mathrm{gr}}$ and ${\sigma }^0$. A remark about the compensation of the effects of the azimuth variation of the backscattering on the DCA is now in order. The Doppler shift is estimated by exploiting the spectral tapering of the azimuth antenna beam pattern. The algorithms for the DC estimation assume a backscattering spectral flatness, that is, a flat backscattering spectral power density [43,44]. Evidently, the presence of azimuth variations of the backscattering within the Doppler Resolution Cell (DRC) can negatively impact this assumption leading to a spatially varying Doppler shift bias, that is, a disturbance in the measured DCA. In [12], this phenomenon was physically explained by referring to the unbalancing, appearing in the presence of a backscattering azimuth gradient, between the contributions to the Doppler Centroid associated with the scatterers ahead and behind the zero Doppler scatter. In [12], the authors proposed a solution to this problem, which is based on the subtraction of the contribution correlated to the azimuth backscattering gradient from the estimated DCA. A detailed analysis of the effect of azimuth backscattering variations on the estimated DCA, as well as its mitigation, is beyond the scope of the present work and is certainly worth future research, especially for the specific application to lakes. The compensation procedure proposed in [12] was not adopted here due to the possible consequences on the cancellation of the useful signal.

5. Conclusions

In recent years, the literature has reflected increasing awareness of the opportunities of satellite-based monitoring of lake surfaces. In this study, we investigate the possibility of inferring lake surface water velocity from existing methods, exploiting the SAR Doppler Centroid Anomaly. To our knowledge, this is the first application of such methodologies to inland waters. The analysis of ENVISAT SAR images and the related DCA maps for Lake Garda returns evident structures with intense and spatially varying values of backscattering and Doppler frequency. These are compared with the results of a numerical model. In the three dates when a high wind blows uniformly over the lake, the spatial patterns retrieved from SAR resemble those from the numerical model, while this is not the case for the low wind dates. In all cases, the surface velocity retrieved from DCA appears to always be significantly overestimated.

The factors potentially influencing the extraction of surface velocity on the investigated dates are identified as are the intensity of the wind forcing, the bathymetry, and the wind waves’ wavelength. The simulated wind waves are also found to be well correlated with the backscatter amplitude. Some hypotheses on the potential biases to be removed are addressed, to preliminarily evaluate the applicability of existing GMFs.

Lake dynamics can be significantly different to those in open seas, especially in medium-size lakes where the response to the wind forcing varies locally. In our case study, wind space and time variability is mainly responsible for the surface flow field [59]. The cases analyzed so far suggest that the GMFs derived for inferring the surface velocity from DCA in open seas [12] cannot be applied to our case study, as they account for the wind vector solely and assume a fully developed wind waves field. Indeed, none of the cases analyzed show a clear correlation between ground range velocity measured from SAR and the correspondent wind component. On the other hand, the data processed suggest that a specific GMF could be defined for our case study, and potentially for the more general case of medium-size lakes. This function should include parameterizations handling the small scale variations associated with wind forcing, bathymetry and waves. The possibility of including the information carried by the SAR backscatter should be evaluated, and the interpretation of the environmental variables signature on the SAR backscatter should be further improved.

The application to lakes of Doppler methods commonly exploited in the open sea is not only challenging for the site related features, for example, specific environmental conditions, hydrodynamic and topographic aspects, but also from the data processing perspective. Further developments shall be carried out in the processing chain to account for possible size limitations of the lake compared to the Doppler resolution cell. Adaptive strategies based on multi-resolution Doppler estimation could be investigated, as well as approaches for the correction of the bias induced by a non-uniform backscattering within the Doppler resolution cell. Such data processing developments are considered priorities for improving the application of SAR-based current estimation methods to lakes.

The results and the interpretations provided in this work must be deeply validated and massive research activity is required, as well as dedicated field campaigns. In this sense, we hope that our contribution will pave the way for the development of a procedure for the extraction of the current signature from SAR that is adequately re-adapted for lakes.