Sand Dune Dynamics Exploiting a Fully Automatic Method Using Satellite SAR Data

: This work presents an automatic procedure to quantify dune dynamics on isolated barchan dunes exploiting Synthetic Aperture RADAR satellite data. We use C-band datasets, allowing the multi-temporal analysis of dune dynamics in two study areas, one located between the Western Sahara and Mauritania and the second one located in the South Rayan dune ﬁeld in Egypt. Our method uses an adaptive parametric thresholding algorithm and common geospatial operations. A quantitative dune dynamics analysis is also performed. We have measured dune migration rates of 2–6 m / year in the NNW-SSE direction and 11–20 m / year NNE-SSW for the South Rayan and West-Sahara dune ﬁelds, respectively. To validate our results, we have manually tracked several dunes per study area using Google Earth imagery. Results from both automatic and manual approaches are consistent. Finally, we discuss the advantages and limitations of the approach presented.


Introduction
An understanding of dune dynamics becomes critical in cases where dunes represent a real hazard due to their proximity to villages or agricultural fields and when dunes move towards them, as is the case for [1] and [2][3][4] respectively.
In order to increase our understanding of this, the accumulation of windblown sand into sand dunes has been the subject of numerous studies using both traditional field surveys and remote sensing techniques. Field surveys may include the use of measuring tape roads [5,6], optical and electronic levelling [7], DGPS [8][9][10], RTK-GPS [11,12], total station [13], terrestrial laser scanning [14] and ground penetrating radars [8,15]. However, field surveys are expensive and have limited spatial coverage and revisits compared to satellite remote sensing techniques, which can cover very wide areas systematically, repetitively, and at a very low cost. Nowadays, satellite data is systematically acquired globally with a repetition frequency that can range from 1 to 16 days with spatial resolutions from below a meter up to tens of meters, such as PlanetScope, Sentinel-1/2 or Landsat 8 [16]. Furthermore, satellite RS has the advantage of being able to analyze less accessible dune regions and/or regions with extreme climatic conditions that may restrict field surveys. This even includes extraterrestrial dune systems, such as those on Mars [17] or in the equatorial regions of Titan [18][19][20][21][22][23].
Monitoring sand dunes on Earth is mostly performed by optical RS data, applying multi-spectral datasets [24], using multi-temporal false RGB color techniques [4], the manual delineation of dunes, temporal acquisitions, e.g. COSI-CORR [26]. The potential of microwave (SAR) sensors for monitoring dune dynamics has not been fully exploited, despite being able to acquire data independently from daylight, cloud coverage and weather conditions. Early works focused on deriving linear dune attributes [27], such as dune height using JERS-1 and ERS-1 data, studies on detecting dunes using SAR or SRTM data [28,29], the correlation of SAR images for the automatic detection of dune areas [30] and more recently, the use of interferometric SAR techniques [31,32]. In this work we proceed further by automatizing the full dune dynamics' analysis.
As demonstrated in [29], dunes can be identified easily in radar images, appearing darker than their surroundings when they are located in either rougher or vegetated environments. In optical images, sand dunes are more difficult to delineate as their composition can be the same as their surroundings (see Figure 1). How dark they appear within an SAR image can vary depending on the incidence angle, signal wavelength, dune type and orientation, and their being able to return a higher signal that becomes a white pixel (like a corner reflector) on the dunes' crests and edges (see Figure  1). However, the predominant dune characteristic is that they appear darker than their environment. Furthermore, even when interdune areas are also characterized by sand, with the same optical properties as the dune itself, the dune can still be detected on SAR images when the SAR signal attenuation produced in the dune returns a backscatter signal much lower than its surroundings. A recent study [33] modelled the radar backscatter response of sand-covered objects to radar signals, measuring the backscatter response of a surface covered with very dry sand, showing a predominant volume scatter mechanism with values lower than −15 dB for sigma0 VV channel. The potential and advantages of SAR imagery over in situ and optical methods is related to the geomorphology of dunes and the difficulties of monitoring their dynamics. Sand dunes may be classified by their (i) size and shape, (ii) location (coastal, desert, polar), (iii) growth stages and degree of complexity and (iv) the wind direction responsible for their formation [34]. Additionally, a geomorphological classification groups dunes by shape, number of slip faces and the wind directions that form them, resulting in six categories: barchan, transverse, barchanoid, longitudinal, parabolic and star [35].
In this work, we focus on the dune dynamics of barchan dunes (see Figure 2), which are characterized by a crescent shape with a concave slip face and "horns" or arms extending downwind [34]. Barchan dunes can form when the terrain is flat while winds blow from one dominant direction (15 degrees difference or less), vegetation cannot grow and sand is available but limited [34]. If the sand supply increases, barchan dunes begin to connect with others, forming barchanoid ridges and The potential and advantages of SAR imagery over in situ and optical methods is related to the geomorphology of dunes and the difficulties of monitoring their dynamics. Sand dunes may be classified by their (i) size and shape, (ii) location (coastal, desert, polar), (iii) growth stages and degree of complexity and (iv) the wind direction responsible for their formation [34]. Additionally, a geomorphological classification groups dunes by shape, number of slip faces and the wind directions that form them, resulting in six categories: barchan, transverse, barchanoid, longitudinal, parabolic and star [35].
In this work, we focus on the dune dynamics of barchan dunes (see Figure 2), which are characterized by a crescent shape with a concave slip face and "horns" or arms extending downwind [34]. Barchan dunes can form when the terrain is flat while winds blow from one dominant direction (15 degrees difference or less), vegetation cannot grow and sand is available but limited [34]. If the Remote Sens. 2020, 12, 3993 3 of 22 sand supply increases, barchan dunes begin to connect with others, forming barchanoid ridges and eventually transforming into transverse dunes. Barchan dunes have the ability to migrate long distances with only a minor change in form when the above conditions are in place [36].
Remote Sens. 2020, 11, x FOR PEER REVIEW 3 of 24 eventually transforming into transverse dunes. Barchan dunes have the ability to migrate long distances with only a minor change in form when the above conditions are in place [36].
In this study, we propose a fully automatic method to analyze the dynamics of isolated barchan dunes using SAR data. We have applied this method over two dune fields with multi-temporal Cband SAR satellites, allowing us to identify a consistent, efficient and reliable way to derive dune shapes, dunes migration rates, and directions, that is applicable globally. We focus our analysis only on this type of dune for the sake of simplicity and to test the proposed approach.

Study Area
To test the proposed methodology, we selected two sandy deserts in Western Sahara-Mauritania (WSM) and Egypt (see Figure 3). We will focus our work on the dynamics of isolated barchan dunes found at the two sites. These frequently propagate as a group, sometimes interacting with one another through collisions and indirect sand exchange [38], making an automatic analysis more difficult. Isolated barchans in equilibrium move without changing their shape [36], enabling the usage or the proposed automatic analysis approach.
Both sites have different topography and wind conditions. The yearly average wind in the WSM area is much higher than in Egypt, being, N-NE 4-6 and NW 1-2, corresponding to 20-50 km/h and 1-10 km/h, respectively [39]. The WSM site comprises a sand dune field with barchan dunes over an almost flat topography with an average slope of 7 m/km in a westerly direction [40] and its substrate layer is a desert pavement formed by gravel and coarse sand [41]. The South Rayan dune field (SRDF) in Egypt lays over a substrate similar to WSM desert pavement, also formed by gravel and coarse sand [42,43] and starts with linear dunes that develop into isolated barchan dunes, which encroach on the western part of the Nile valley. Due to the tilted topography with an average slope of 50 m/km in an easterly direction, the barchan dunes are not as symmetric as those from WSM, but have an elongated eastern horn, as illustrated in Figure 3G [44].

Materials
The selected datasets are from previous and current European SAR satellites, i.e., the European Remote Sensing satellite (ERS), the Environmental Satellite (Envisat) and the Copernicus Sentinel-1 (S-1). Due to the spatial resolutions (about 20 m for all three sensors) and the expected dune migration rate on the study area, we selected time intervals that were distant enough to measure the dune movement. Detailed information on the dataset employed is provided in Table 1. In this study, we propose a fully automatic method to analyze the dynamics of isolated barchan dunes using SAR data. We have applied this method over two dune fields with multi-temporal C-band SAR satellites, allowing us to identify a consistent, efficient and reliable way to derive dune shapes, dunes migration rates, and directions, that is applicable globally. We focus our analysis only on this type of dune for the sake of simplicity and to test the proposed approach.

Study Area
To test the proposed methodology, we selected two sandy deserts in Western Sahara-Mauritania (WSM) and Egypt (see Figure 3). We will focus our work on the dynamics of isolated barchan dunes found at the two sites. These frequently propagate as a group, sometimes interacting with one another through collisions and indirect sand exchange [38], making an automatic analysis more difficult. Isolated barchans in equilibrium move without changing their shape [36], enabling the usage or the proposed automatic analysis approach.
Both sites have different topography and wind conditions. The yearly average wind in the WSM area is much higher than in Egypt, being, N-NE 4-6 and NW 1-2, corresponding to 20-50 km/h and 1-10 km/h, respectively [39]. The WSM site comprises a sand dune field with barchan dunes over an almost flat topography with an average slope of 7 m/km in a westerly direction [40] and its substrate layer is a desert pavement formed by gravel and coarse sand [41]. The South Rayan dune field (SRDF) in Egypt lays over a substrate similar to WSM desert pavement, also formed by gravel and coarse sand [42,43] and starts with linear dunes that develop into isolated barchan dunes, which encroach on the western part of the Nile valley. Due to the tilted topography with an average slope of 50 m/km in an easterly direction, the barchan dunes are not as symmetric as those from WSM, but have an elongated eastern horn, as illustrated in Figure 3G

Materials
The selected datasets are from previous and current European SAR satellites, i.e., the European Remote Sensing satellite (ERS), the Environmental Satellite (Envisat) and the Copernicus Sentinel-1 (S-1). Due to the spatial resolutions (about 20 m for all three sensors) and the expected dune migration rate on the study area, we selected time intervals that were distant enough to measure the dune movement. Detailed information on the dataset employed is provided in Table 1.

Method
The key element of our automatic method consists of identifying the dunes using SAR images, assuming that in an SAR image, dunes look like dark areas surrounded by brighter pixels (Figure 4).
The method is divided into three steps (see Figure 5): (i) SAR data preprocessing; (ii) dune identification using the Hierarchical Split Based Approach (HSBA) algorithm [28] and; (iii) dune movement computation using Geographic Information System command line tools and validation. The validation is performed using existing data and Google Earth and validated data over the South Rayan dune made available from [43].

Method
The key element of our automatic method consists of identifying the dunes using SAR images, assuming that in an SAR image, dunes look like dark areas surrounded by brighter pixels (Figure 4). The method is divided into three steps (see Figure 5): i) SAR data preprocessing; ii) dune identification using the Hierarchical Split Based Approach (HSBA) algorithm [28] and; iii) dune movement computation using Geographic Information System command line tools and validation. The validation is performed using existing data and Google Earth and validated data over the South Rayan dune made available from [43].

SAR Data Preprocessing
SAR data has been calibrated and geocoded using the ESA calibration method incorporated in the Sentinel-1 Toolbox [45] integrated into the SeNtinel Application Platform (SNAP), using the Precise Orbits and the SRTM DEM for accurate geocoding determination. Output SAR data has been converted to a decibel scale, enhancing the contrast of the image, and making the backscattering distributions more symmetrical and Gaussian, suitable for the subsequent identification step.

Dune Identification Using the HSBA Algorithm
This step uses the calibrated and geocoded SAR data in the decibel scale as input and, using the HSBA algorithm, produces a raster binary image where dunes and groups of dunes, from now on called "sand bodies", are identified. The HSBA algorithm assumes that dune backscattering is very low, positioning its probability density function (PDF) in the lowest part of the image histogram (see Figure 6B,D).
To separate the target class from the background, a thresholding approach is frequently used [46]. The drawback of this method is that the classification relies heavily on the adequacy of the selected threshold. To this end, parametric thresholding algorithms are preferable because they estimate PDFs of the target class and its background, and based on these, the threshold is selected [47]. In this context, one of the main difficulties in parameterizing these functions originates from the fact that the target class often represents only a small fraction of the image. Under such circumstances, the histogram of the image values is often not obviously bimodal, and it becomes difficult to parameterize PDFs. In the case of water detection, HSBA searches for tiles of variable sizes allowing the parameterization of PDFs of two classes [48]. This is applied here for dune detection (see Figure 6). The approach has also been successfully applied to detect buildings on a global scale [48,49] and volcanic lava flow-induced changes [50] from SAR intensity images.
To fit the PDF of the dunes class (PDF D ), we use HSBA, a statistically based algorithm which makes use of a hierarchical tiling of the image. Based on PDF D , we combine histogram thresholding and region-growing processes to identify the dunes. The parameters of the region-growing and thresholding processes are automatically derived from PDF D. The definition of PDF D starting from the entire SAR acquisition is possible if the class is sufficiently represented, i.e., identifiable, and this generally depends on the shape of the histogram. The PDF D may not be easily fit from the histogram when dunes represent only a small percentage of the entire image. Therefore, it is necessary to focus on those areas of the image that are composed of a similar number of pixels belonging to dunes and background (PDF B ) classes, respectively. To this end, HSBA has been used to automatically identify regions in any given SAR image where the two PDFs of dunes and background are well-separated and give rise to a bimodal histogram of the region. The main scope of this is to obtain a robust parameterization of the PDF D and PDF B , making the classification of the dune class more reliable. Here, the PDF D and PDF B are assumed to be Gaussian, motivated by the fact that the input data are multi-looked and log-transformed SAR intensity images to increase the equivalent number of looks (ENL) and to have more Gaussian distributions [51].

Method
The key element of our automatic method consists of identifying the dunes using SAR images, assuming that in an SAR image, dunes look like dark areas surrounded by brighter pixels (Figure 4). The method is divided into three steps (see Figure 5): i) SAR data preprocessing; ii) dune identification using the Hierarchical Split Based Approach (HSBA) algorithm [28] and; iii) dune movement computation using Geographic Information System command line tools and validation. The validation is performed using existing data and Google Earth and validated data over the South Rayan dune made available from [43].   The proposed approach is described in detail in Chini et al. [48] In HSBA, a hierarchical tiling of the scene is initiated starting with 4 0 tiles (i.e., the entire image) on the first level and then continuing by iteratively subdividing the image into 4 L sub images, with L being the hierarchical level of splitting. At L = 1, the image is split into quarters; with L = 2, the image is subdivided into sixteenths; and so on. Depending on L, the tiles will thus be characterized by their different sizes. At each level, descending from the upper level to the lower one, only tiles fulfilling fixed criteria are selected, while the others will be further split.
The criteria to select a tile are: (1) The histogram of the tile has to be bimodal and composed of two Gaussian PDFs ( Figure 6D); (2) The number of pixels belonging to a dunes class must be composed at least of 20% of the considered tile; (3) The mode of PDF D has to be lower than a predefined value, which is set to −17 dB, because the backscattering of the dune class is expected to be low.
To fit PDF D and PDF B from the histogram of a given tile, we use the Levenberg-Marquardt algorithm, a technique to solve nonlinear least square problems [52]. To check the bimodality of the histogram, we use the Ashman D coefficient [53], which quantifies how well two Gaussian distributions are separated, e.g., PDF D and PDF B , by considering the distance between their main modes and their dispersions, i.e., standard deviations. To consider the PDFs as well-separated, the Ashman D coefficient has to be higher than 2.
All tiles selected based on the three previously defined criteria are merged together and the resulting histogram is used to fit the final PDF D and PDF B to be used in the next steps for thresholding the image. Although the resulting PDFs are well-separated, some overlap is still present, consequently setting the threshold where PDF D is equal, PDF B can produce some over-and under-detection. To reduce these latter drawbacks, contextual information is also used via a region-growing approach ( [48,54]). The region-growing algorithm starts from seed pixels and adds connected pixels to them that lie within a predefined tolerance value. The choice of the threshold values for the seeds and for stopping the growth is a crucial point. The strategy we follow to select these two parameters is based on PDF D . We select seed pixels with a high likelihood of belonging to the dunes class, e.g., pixels with backscattering values lower than the PDF D mode, while many different thresholds to stop the growth are tested, and we select the one that minimizes the RMSE between the theoretical PDF D and the empirical histogram resulting from the region-growing.

Dune Movement Using GIS Functions
This last step uses a pair of binary images obtained after the dune identification step as input and computes the movement of individual dunes based on centroid differences of the intersecting dunes (see Figure 7) as done in [43]. The algorithm assumes that a dune shifts by less than its width in the direction of movement, making an adequate SAR time sampling necessary. This latter will be guaranteed by new satellite constellations, which are able to systematically image the surface with a very small repeat cycle, e.g., the six days of Sentinel-1.
Remote Sens. 2020, 11, x FOR PEER REVIEW 8 of 24 modes and their dispersions, i.e., standard deviations. To consider the PDFs as well-separated, the Ashman D coefficient has to be higher than 2. All tiles selected based on the three previously defined criteria are merged together and the resulting histogram is used to fit the final PDFD and PDFB to be used in the next steps for thresholding the image. Although the resulting PDFs are well-separated, some overlap is still present, consequently setting the threshold where PDFD is equal, PDFB can produce some over-and underdetection. To reduce these latter drawbacks, contextual information is also used via a region-growing approach ( [48], [54]). The region-growing algorithm starts from seed pixels and adds connected pixels to them that lie within a predefined tolerance value. The choice of the threshold values for the seeds and for stopping the growth is a crucial point. The strategy we follow to select these two parameters is based on PDFD. We select seed pixels with a high likelihood of belonging to the dunes class, e.g., pixels with backscattering values lower than the PDFD mode, while many different thresholds to stop the growth are tested, and we select the one that minimizes the RMSE between the theoretical PDFD and the empirical histogram resulting from the region-growing.

Dune Movement Using GIS Functions
This last step uses a pair of binary images obtained after the dune identification step as input and computes the movement of individual dunes based on centroid differences of the intersecting dunes (see Figure 7) as done in [43]. The algorithm assumes that a dune shifts by less than its width in the direction of movement, making an adequate SAR time sampling necessary. This latter will be guaranteed by new satellite constellations, which are able to systematically image the surface with a very small repeat cycle, e.g., the six days of Sentinel-1. The functions employed in this step are: i) fill the holes-of the binary images; ii) polygonize them; iii) compute x and y centroid coordinates for the different polygons; iv) find the intersecting polygons at two dates; and v) estimate movement and direction. Directions of dune movement range, clockwise, from 0 to 360 degrees, where 0 degree corresponds to the north direction. The QGIS processing modeler [55] was employed to define and create the automatic processing employed in this step (see Figure 8). The functions employed in this step are: (i) fill the holes-of the binary images; (ii) polygonize them; (iii) compute x and y centroid coordinates for the different polygons; (iv) find the intersecting polygons at two dates; and (v) estimate movement and direction. Directions of dune movement range, clockwise, from 0 to 360 degrees, where 0 degree corresponds to the north direction. The QGIS processing modeler [55] was employed to define and create the automatic processing employed in this step (see Figure 8).
Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 24 The results obtained are finally filtered based on an area consistency criterion that we defined in line with the area consistency principle of the isolated barchan, which moves conserving its shape and hence also its area. This area consistency criterion can be translated as the dune area loss between the observed dates, as illustrated in Equation (1). Based on that criterion, variations in the area measured using our method can be values ranging in the interval [0,1], which can help us discriminate individual dunes from other cases such as dune merging and dune splitting (see Figure  9), which could produce high variability of the observed dune area. The area consistency criterion measures the area loss, given by the following equation: (1) The results obtained are finally filtered based on an area consistency criterion that we defined in line with the area consistency principle of the isolated barchan, which moves conserving its shape and hence also its area. This area consistency criterion can be translated as the dune area loss between the observed dates, as illustrated in Equation (1). Based on that criterion, variations in the area measured using our method can be values ranging in the interval [0, 1], which can help us discriminate individual dunes from other cases such as dune merging and dune splitting (see Figure 9), which could produce high variability of the observed dune area. The area consistency criterion measures the area loss, given by the following equation: where the area is the dune area for each time observed t1 and t2 and measures the area loss between two observations compared to the maximum area observed. The results obtained are finally filtered based on an area consistency criterion that we defined in line with the area consistency principle of the isolated barchan, which moves conserving its shape and hence also its area. This area consistency criterion can be translated as the dune area loss between the observed dates, as illustrated in Equation (1). Based on that criterion, variations in the area measured using our method can be values ranging in the interval [0,1], which can help us discriminate individual dunes from other cases such as dune merging and dune splitting (see Figure  9), which could produce high variability of the observed dune area. The area consistency criterion measures the area loss, given by the following equation: where the area is the dune area for each time observed t1 and t2 and measures the area loss between two observations compared to the maximum area observed.

Validation
We have manually delineated dunes from both the WSM and SRDF areas for two different time periods using very high resolution optical data available in Google Earth, enabling the comparison of the results obtained from our automatic approach with the manual delineation and analysis, in terms of dune locations, and movement characteristics such as distance, average velocity, direction and area loss for both delineated and satellite-derived data.
Additionally, and for the specific case of the South Rayan dune field, we also used the existing dataset obtained from [43].

Validation
We have manually delineated dunes from both the WSM and SRDF areas for two different time periods using very high resolution optical data available in Google Earth, enabling the comparison of the results obtained from our automatic approach with the manual delineation and analysis, in terms of dune locations, and movement characteristics such as distance, average velocity, direction and area loss for both delineated and satellite-derived data.
Additionally, and for the specific case of the South Rayan dune field, we also used the existing dataset obtained from [43].

West Sahara-Mauritania
We applied the aforementioned methodology to the West Sahara dataset, obtaining after our dune identification step the detected sand bodies for the different satellite acquisitions. Figure 10 shows one example of the sand body identification obtained using Envisat data acquired in 2003. Figure 10B illustrates an Envisat ASAR image overlaid with the detected sand bodies in brown. Figure 10C shows a small subset of the SAR image overlaying it with the contour of such detections.

West Sahara-Mauritania
We applied the aforementioned methodology to the West Sahara dataset, obtaining after our dune identification step the detected sand bodies for the different satellite acquisitions. Figure 10 shows one example of the sand body identification obtained using Envisat data acquired in 2003. Figure 10B illustrates an Envisat ASAR image overlaid with the detected sand bodies in brown. Figure 10C shows a small subset of the SAR image overlaying it with the contour of such detections. In addition, we manually delineated six dunes over VHR optical data on 2 dates, one in 2003, matching one of our datasets, and the other in 2013, as this was the latest data available in Google Earth for that area (see Figure 11). By doing this, we obtained a realistic and true threshold for the area loss computation (Equation (1)) to use in our study. Moreover, we were able to compute their long-term movement characteristics. Furthermore, Figure 11G shows a small area of the WSM site where the six dunes employed as ground truth are located, showing the crescent shape of the barchan dunes, which is symmetric in most cases. In addition, we manually delineated six dunes over VHR optical data on 2 dates, one in 2003, matching one of our datasets, and the other in 2013, as this was the latest data available in Google Earth for that area (see Figure 11). By doing this, we obtained a realistic and true threshold for the area loss computation (Equation (1)) to use in our study. Moreover, we were able to compute their long-term movement characteristics. Furthermore, Figure 11G shows a small area of the WSM site where the six dunes employed as ground truth are located, showing the crescent shape of the barchan dunes, which is symmetric in most cases.
Besides this, we computed the distance, velocity, area loss and direction of the delineated dunes, which are found in Table 2, confirming the area consistency principle of isolated barchan dunes (with area loss < 0.2), and also confirming the average direction obtained using the automatic method (direction SSE, from 180 to 190 degrees). This threshold (area loss < 0.2) was used to identify the isolated barchan dunes in equilibrium among all detected sand bodies with overlap in Table 3.
Using our automatic method, we have detected a total amount of sand bodies (not only isolated dunes) which varies from one image to another ( Table 3). The biggest differences among the detected sand bodies correspond to images acquired with different sensors and without the exact same geographic coverage (Table 3).
We applied the area consistency criterion to Equation (1), ensuring that the dune measurements reported are obtained from isolated barchan dunes, whose characteristic is its invariant shape in an ideal situation, the results of which are reported in Table 3. We obtained average dune velocities ranging from 15-26 m/yr for the analyzed area. Note that the resulting number of isolated dunes with an overlap between the dates is much lower than the total amount of detected sand bodies. Figure 12 illustrates the histogram of the heading direction of movement of these isolated barchans in WSM.     Looking at the values reported in Table 3 and Figure 12, it seems that there could be an error in the average heading obtained, which could also interfere with the velocity calculated. We will analyze this in the discussion section.
The values obtained using the proposed method for the delineated dunes are shown in Table 4. There is a high level of consistency between the velocities obtained with the automatic method and the ground truth, as well as the heading (max difference of 7 degrees) and distance. Despite this consistency, there are differences in the mean area detected, and these are higher for the smallest dunes (C) detecting almost 36% less compared to 12-20% for the others. These differences in the detected area are linked to: (i) the spatial resolution of the sensor (20 m × 20 m); and (ii) the non-detected areas being located in the borders of the dunes, as a certain volume of sand accumulated is needed in order to appear dark in the SAR data. The proposed method could not provide information for dune B in an automatic manner due to its lack of overlap in the analyzed period. However, with the scope of the comparison, its dune migration information (normally performed in step 3) was measured manually and is also presented in Table 4.

South Rayan Dune Field in Egypt
We also applied our methodology to the SRDF dataset. A general overview of the results obtained in the dune detection step is illustrated in Figure 13C, and Figure 13A shows the Google Earth overview of the SRDF (within the red polygon). Figure 13B illustrates how dunes can be easily identified using SAR data, compared to optical data ( Figure 13A), as they appear darker than their surroundings, as explained in the introduction section.

South Rayan Dune Field in Egypt
We also applied our methodology to the SRDF dataset. A general overview of the results obtained in the dune detection step is illustrated in Figure 13C, and Figure 13A shows the Google Earth overview of the SRDF (within the red polygon). Figure 13B illustrates how dunes can be easily identified using SAR data, compared to optical data ( Figure 13A), as they appear darker than their surroundings, as explained in the introduction section.
Furthermore, we delineated three dunes using Google Earth, selecting the ones which were specially monitored and validated with GPS and VHR optical data in [43]. These three manually delineated dunes in the South Rayan dune field are illustrated in Figure 14 and their values obtained are shown in Table 5. Table 6 reports the information over the same three dunes observed in [43].
Applying our automatic method, we have detected a total amount of sand bodies reported in Table 7. Again, the biggest differences among the detected sand bodies with overlap corresponded to images acquired with different sensors, and not always to the same exact geographical extent.
According to [43], the average velocity of the dune movement for the SRDF is expected to be around 2-6 m/yr . Similar values were found for some of the analyzed periods reported in Table 7. These values can be affected by the merging or splitting of several dunes [56]. This was observed when smaller dunes moving faster reached larger dunes moving at a slower pace [56]. Additionally, dunes in the SRDF could have been affected by new fields that had appeared in the proximities of the dune field, even surrounding dunes, in the Western Desert.   Furthermore, we delineated three dunes using Google Earth, selecting the ones which were specially monitored and validated with GPS and VHR optical data in [43]. These three manually delineated dunes in the South Rayan dune field are illustrated in Figure 14 and their values obtained are shown in Table 5. Table 6 reports the information over the same three dunes observed in [43].  Applying our automatic method, we have detected a total amount of sand bodies reported in Table 7. Again, the biggest differences among the detected sand bodies with overlap corresponded to images acquired with different sensors, and not always to the same exact geographical extent. According to [43], the average velocity of the dune movement for the SRDF is expected to be around 2-6 m/yr. Similar values were found for some of the analyzed periods reported in Table 7. These values can be affected by the merging or splitting of several dunes [56]. This was observed when smaller dunes moving faster reached larger dunes moving at a slower pace [56]. Additionally, dunes in the SRDF could have been affected by new fields that had appeared in the proximities of the dune field, even surrounding dunes, in the Western Desert.
In Table 6, information extracted from [43] on the dunes A, B and C are reported and it is possible to appreciate that both velocity and area are consistent with the dunes extracted from Google Earth. Figure 14 shows the dunes manually delineated using Google Earth from 2007 and 2010, in green and red, respectively. In the middle image (D), we show the delineated dunes over the Envisat ASAR data acquired on 2010/08/24, highlighting the high consistency between the outline polygons of the dunes on the left (A, B and C) and the black dunes in the SAR image. In the right-hand diagram (E), the Envisat image is overlapped with the outline of the dunes detected using our automatic method. We can observe some differences in the detected shapes in D and E, and this will be discussed in Section 5.
From Figure 14E, we can see that detected dunes do not exactly match the delineated ones; some dunes do not appear totally black in the SAR images, as they have white pixels due to a high backscatter produced over the sharp edges created with the help of the topography (average slope of 50 m/km). These bright spots on the dunes were not encountered in the uniform and regular barchan shapes on the WSM dune field (see Figure 11).
In Table 7, the average values for distance, velocity and direction obtained for isolated dunes are reported after applying the area consistency criterion. For the periods from 2000 to 2004 and from 2014 to 2019, values are higher than in the other periods observed. In the South Rayan dune field, dunes move slower, and in some cases, the movement observed is near the limit of the technique sensitivity. A good area detection becomes critical when the expected movement is only 1 or 2 pixels (with 20m pixel spacing). In this case, the topography influences the dune morphology, creating sharp edges on dunes. These sharp edges appear as white pixels in the SAR image, and since they are not detected as dunes, they affect the dune area measured and finally its movement. However, long-term migration is well detected, both in speed and direction.
Additionally, in Table 7, we analyzed the dune migration within consecutive observations, as well as the observations with more time difference for both sensor changes, i.e., ERS2-Envisat (2000-2010) and ERS2-S1 (2000-2019). It is also important to highlight that different sensor data had different extents, for which the maximum spatial overlap is obtained within the same sensor, i.e., Envisat-Envisat (2004-2010) and S1-S1 (2014-2019) where we detected a higher number of individual smaller dunes.
For the dunes remaining after applying the area consistency criterion, in order to further validate our results, we compared our dune map with the data from [43], obtaining 30 overlapping dunes. From these 30 dunes, we analyzed the variation in their area detection for the different dates, obtaining the results illustrated in Figure 15. Finally, our method provided dune migration information also for the dunes A, B and C, whose values are reported in Table 8. These values are not totally consistent with those found in the literature, probably due to the detected dune area being 85%, 80% and 65% of the area detected in the SAR images for dunes A, B and C, respectively. This underestimation of the dune area directly affects their centroid coordinates and hence the distance and velocity. Despite this, the values obtained are still in the range of the expected dune migration and close to the ground-truth values. Table 8. Dune information provided by the automatic method on SAR data acquired in 2004 and 2010 for the same dunes which were also manually delineated, located in the South Rayan dune field.

Discussion
The total number of detected sand bodies reported in Table 3 and Table 7 showed high variability among the different periods. One of the possible causes for this could be explained by the differences in area detected which are minimized after applying the proposed criterion.
An intrinsic characteristic of free dunes belonging to the same dune field is that the distance that the dunes rove has an inverse relationship with their size, i.e., small dunes move for longer and are faster than bigger ones (see Figure 16), with the characteristic of keeping their shape invariant in the event of barchans in equilibrium [36]. Dunes can fuse when smaller and faster dunes reach bigger and slower ones, and this fusion can be temporary or permanent [56]. Figure 16 illustrates the curve fitting of the aforementioned characteristic (the distance that the dunes rove has an inverse relationship with their size) of our manually delineated dunes on the WSM site (in blue) and the dunes from SRDF delineated in [43] (in orange), and clearly shows the differences in area and migration rate of the dunes of both analyzed sites. Finally, our method provided dune migration information also for the dunes A, B and C, whose values are reported in Table 8. These values are not totally consistent with those found in the literature, probably due to the detected dune area being 85%, 80% and 65% of the area detected in the SAR images for dunes A, B and C, respectively. This underestimation of the dune area directly affects their centroid coordinates and hence the distance and velocity. Despite this, the values obtained are still in the range of the expected dune migration and close to the ground-truth values.

Discussion
The total number of detected sand bodies reported in Tables 3 and 7 showed high variability among the different periods. One of the possible causes for this could be explained by the differences in area detected which are minimized after applying the proposed criterion.
An intrinsic characteristic of free dunes belonging to the same dune field is that the distance that the dunes rove has an inverse relationship with their size, i.e., small dunes move for longer and are faster than bigger ones (see Figure 16), with the characteristic of keeping their shape invariant in the event of barchans in equilibrium [36]. Dunes can fuse when smaller and faster dunes reach bigger and slower ones, and this fusion can be temporary or permanent [56]. Figure 16 illustrates the curve fitting of the aforementioned characteristic (the distance that the dunes rove has an inverse relationship with their size) of our manually delineated dunes on the WSM site (in blue) and the dunes from SRDF delineated in [43] (in orange), and clearly shows the differences in area and migration rate of the dunes of both analyzed sites. Small dunes rove faster than larger dunes, as illustrated in Figure 16. This may produce their non-intersection between two time periods, which excluded them from our analysis due to the non-intersection condition. In order to provide dune migration information, our method needs each dune shape to intersect in space for the analyzed time periods. Nowadays, and due to the revisit times of the new satellite constellations such as Sentinel-1, we could also analyze faster dunes, just by decreasing the time interval between images. Barchan dunes located in the almost flat WSM area have a more regular shape, helping the proper detection of the dune area. Additionally, the average velocity measured on the intersected dunes is over 18 m/yr, and almost as fast as the barchan dunes in Chad, considered some of the fastest dunes on Earth, moving at an average velocity of about 20 m/yr [57]. Additionally, we measured an average ratio of the dune migration rates (automatically computed versus ground truth) of 0.96, being very consistent with the estimated migration rates. Moreover, we also computed the differences in the dune area measured, which varied from 11% to 36% of area difference compared to the ground truth, being more critical for smaller dunes. This is directly linked to the spatial resolution of the data employed. However, this misdetection of the dune area seems systematic and is located on the dune edges, where thinner sand layers are located, allowing us to compute reliable migration rates despite the dune area difference.
We also computed the dune migration information for the different periods analyzed for the intersecting dunes in the SRDF dune field, obtaining an average movement of 4-5 m/yr in the SSE direction in agreement with [43], which reported an average migration rate of 4.4 m/yr also in the SSE direction. Additionally, we observed that the dune area on the SRDF site is slightly underestimated, probably due to the bright pixels that appear, preventing better detection (see Figure 1). However, long-term migration information is properly detected.
We also noticed a misalignment of the ERS-2 SAR data compared to the Envisat ASAR and Sentinel-1 dataset. This is visible when computing the dunes' displacement direction for the period 1999 to 2003 in Figures 12B and 17C. This shift, we found, could not be avoided and could be attributed to the precision of the orbital information available and employed during the pre-processing of the Hence, despite the fact that dunes may appear slightly smaller in SAR images compared to optical images, we have demonstrated that SAR data can be used to analyze dune dynamics at the dune field scale, as shown in the consistency of the dune migration rates between our results and the ground-truth data.

Overall Comments for the Proper Exploitation of the Proposed Approach
For all the sandy areas detected, we discarded larger volumes of sand corresponding to dune groups in our analysis, analyzing only individual dunes for which the detected area has varied less than 20%. By doing this, additional dunes are discarded, for example, when such dunes join and separate from other dunes through the analyzed period.
Considerations identified for a proper usage of our approach can be grouped into two categories: i) dune detection; and ii) dune migration analysis. For the former, obtaining an accurate detection is linked directly to the spatial resolution of the employed sensor and to bright pixels that sometimes appear on dunes, while for the latter, dune migration information can only be provided for isolated dunes that have spatial overlap on the analyzed images.
To maximize the number of dunes with overlap, it is suggested that the number of images to be processed is increased, thus reducing the time interval between two consecutive ones, in order to ensure a geographical overlap of the moving dunes. This can be foreseen by knowing the wind speed of the area of interest and by selecting the proper dataset with spatial-temporal resolution to ensure the success of the method presented here.
As seen in the case of the South Rayan dune field, dunes located there are more likely to show SAR bright pixels on their large tail, which our approach sometimes cuts off from the dune, affecting their detected area, and hence their movement. However, this has not been found in dunes located in the West Sahara-Mauritania dune field.
It is important to note that on the SRDF site, for some cases we have tried to measure movements near the limit of the detection error, as one pixel counts for 20 meters, and this is, on average, the distance an average dune in that field will move in four4 years, at the speed of 5 m/yr. Hence, we obtained better results for longer periods (i.e., 2000 to 2019 or 2004 to 2014).
Moreover, special attention must be paid to the proper geolocation of the satellite data, to avoid misalignments among data coming from different sensors, as is the case for the older dataset employed, i.e., ERS-2 SAR, and as we have shown happens on the WSM test site, but not on the SRDF. Hence, despite the fact that dunes may appear slightly smaller in SAR images compared to optical images, we have demonstrated that SAR data can be used to analyze dune dynamics at the dune field scale, as shown in the consistency of the dune migration rates between our results and the ground-truth data.

Overall Comments for the Proper Exploitation of the Proposed Approach
For all the sandy areas detected, we discarded larger volumes of sand corresponding to dune groups in our analysis, analyzing only individual dunes for which the detected area has varied less than 20%. By doing this, additional dunes are discarded, for example, when such dunes join and separate from other dunes through the analyzed period.
Considerations identified for a proper usage of our approach can be grouped into two categories: (i) dune detection; and (ii) dune migration analysis. For the former, obtaining an accurate detection is linked directly to the spatial resolution of the employed sensor and to bright pixels that sometimes appear on dunes, while for the latter, dune migration information can only be provided for isolated dunes that have spatial overlap on the analyzed images.
To maximize the number of dunes with overlap, it is suggested that the number of images to be processed is increased, thus reducing the time interval between two consecutive ones, in order to ensure a geographical overlap of the moving dunes. This can be foreseen by knowing the wind speed of the area of interest and by selecting the proper dataset with spatial-temporal resolution to ensure the success of the method presented here.
As seen in the case of the South Rayan dune field, dunes located there are more likely to show SAR bright pixels on their large tail, which our approach sometimes cuts off from the dune, affecting their detected area, and hence their movement. However, this has not been found in dunes located in the West Sahara-Mauritania dune field.
It is important to note that on the SRDF site, for some cases we have tried to measure movements near the limit of the detection error, as one pixel counts for 20 m, and this is, on average, the distance an average dune in that field will move in four4 years, at the speed of 5 m/yr. Hence, we obtained better results for longer periods (i.e., 2000 to 2019 or 2004 to 2014).
Moreover, special attention must be paid to the proper geolocation of the satellite data, to avoid misalignments among data coming from different sensors, as is the case for the older dataset employed, i.e., ERS-2 SAR, and as we have shown happens on the WSM test site, but not on the SRDF. Note that this misalignment does not happen with the newer satellite data, as seen with Envisat and Sentinel-1.
Despite our both dune fields moving mainly in the N-S direction, we consider that our method should work as well in dune fields where the main direction is E-W as there is no constraint in our method that limits the migration direction.
Anywhere where no criterion is applied based on dune area consistency (area loss < 0.2), individual dune migration information might be perturbed.
Future research directions will lead to an enhancement of the method we propose to enable dune matching even for dunes that do not overlap using computer vision algorithms and analyzing dunes in both ascending and descending orbits when possible.

Implications for Other Domains
Climate change studies point towards the increase in wind speed on Earth, which will directly affect the dune dynamics of entire dune fields. Our automatic approach could be used to measure such wind variations of dune fields from all the world, providing results more rapidly and in a systematic manner. This systematic and global monitoring is not possible using other non-automatic techniques, such as traditional field campaigns or using optical data, the limitations of which were already mentioned in the introduction section.
Deriving the dune migration information from multiple dates, and obtaining the dune speed and movement direction could help to derive wind parameters, such as average direction, speed or the possible variations of them over time. This could help us to understand whether these conditions have changed over time during the analyzed periods.
Moreover, extra-planetary research, as already mentioned in the introduction section, has observed dunes on Mars and Titan. Hence, applying the proposed approach when sufficient data is available could help measure wind speeds on their surface. These wind speeds could become crucial pieces of information for the planning of space exploration missions, especially for those planning to land on them.

Conclusions
We have developed a totally automatic approach to measure dune dynamics, and have applied it to different test sites where barchan dunes move freely, but at different speeds and in different directions, showing the robustness of our method.
We have demonstrated the possibility of extending the applicability of techniques that were originally developed for other purposes, such as flood mapping, to dune dynamic analysis, due to the similar behavior of dunes and water to the SAR co-polarized signal, and hence we have also demonstrated the usability of SAR data for dune migration analysis at scale.
The proposed methodology has been demonstrated to work on isolated barchan dunes and provides dune migration information about dunes with an overlap between images. This methodology can be exported to other barchan dune fields. Additionally, as a future work, the applicability of this method to other dune types could be explored.
We have measured, for our analyzed time period, dune migration rates of 11-20 m/yr moving in a south-south westerly direction for the West Sahara/Mauritania study area and 2-6 m/yr moving in a south-easterly direction for the South Rayan dune field, being very consistent with the available ground truth being 13-28 m/yr and 3-5 m/yr respectively.