A New Method for Automated Measurement of Sand Dune Migration Based on Multi-Temporal LiDAR-Derived Digital Elevation Models

While remote sensing methods have long been used for coastal and desert sand dune studies, few methods have been developed for the automated measurement of dune migration in large dune fields. To overcome a major limitation of an existing method named “pairs of source and target points (PSTP)”, this paper proposes a toe line tracking (TLT) method for the automated measurement of dune migration rate and direction using multi-temporal digital elevation models (DEM) derived from light detection and ranging (LiDAR) data. Based on a few simple parameters, the TLT method automatically extracts the base level of a dune field and toe lines of individual dunes. The toe line polygons derived from two DEMs are processed using logical operators and other spatial analysis methods implemented in the Python programming language in a geographic information system. By generating thousands of random sampling points along source toe lines, dune migration distances and directions are calculated and saved with the sampling point feature class. The application of the TLT method was demonstrated using multi-temporal LiDAR-derived DEMs for a 9 km by 2.4 km area in the White Sands Dune Field in New Mexico (USA). Dune migration distances and directions for three periods (24 January 2009–26 September 2009, 26 September 2009–6 June 2010, and 24 January 2009–6 January 2010) were calculated. Sensitivity analyses were carried out using different window sizes and toe heights. The results suggest that both PSTP and TLT produce similar sand dune migration rates and directions, but TLT is a more generic method that works for dunes with or without slipfaces that reach the angle of repose.

While remote sensing methods have been used for studies of coastal and desert sand dunes for decades, the objectivity of visual image interpretation and manual measurement methods has been questioned [36]. In addition, visual image interpretation and manual measurement can be very labor-intensive and time-consuming for large dune fields. In In comparison with traditional medium-and low-resolution remote sensing m light detection and ranging (LiDAR) has provided unprecedented high-resolut high-accuracy digital elevation models (DEM) for sand dune studies, from coast [37][38][39] to desert dunes [40][41][42][43][44]. However, it appears that traditional manual del and measurement of sand dune features is still a common practice, and that the capability of LiDAR in providing three-dimensional information on topographi is not fully utilized. As reviewed by Hugenholtz et al. [45], the progress in spatial methods for characterizing dune field patterns has been limited. A review of LiD for change detection in Earth sciences can be found in Okyay et al. [46].
To facilitate the investigation of large dune fields using LiDAR data, Dong veloped a new approach, pairs of source and target points (PSTP), for the au measurement of sand dune migration rates and directions in a geographic info system (GIS). Using LiDAR data collected for a 2.4 km by 9 km area in the Whi Dune Field (WSDF) in New Mexico, USA on 24 January 2009 and 6 June 2010, D demonstrated that sand dune migration rates and directions in large dune field automatically detected and measured using multi-temporal LiDAR-derived D thousands of locations in minutes. A PSTP software tool was also developed by Dong [47] for the automated measurement of sand dune migration rates and dire A major limitation of the PSTP method developed by Dong [44] is that it obvious dune slipfaces that reach the angle of repose (30° to 34°) in both DEMs. quirement can be difficult to meet if one or both DEMs are created from data c when sand dune slipfaces do not reach the angle of repose due to seasonal winds few slipfaces reached the angle of repose in the DEM derived from LiDAR coll WSDF on 26 September 2009 due to the influence of reverse winds or crosswinds 1 shows three west-east profiles of a dune in WSDF extracted from LiDAR-derive of 24    To overcome the above major limitation of the PSTP method, this paper proposes a new method called toe line tracking (TLT) for the automated measurement of sand dune migration distance and direction using multi-temporal DEMs derived from LiDAR data in a study area in the White Sands Dune Field (WSDF) in New Mexico, USA. Compared with the PSTP method which relies on the angle of repose for sand dune, the new method aims to extract toe lines of sand dunes based on several simple parameters, and conduct automated measurement of dune migration rates and directions through random sampling and spatial analysis in geographic information systems (GIS). The study area, data, methods, results, discussion, and conclusions are presented in the following sections.

Study Area
A study area of 9 × 2.4 km in the White Sands Dune Field (WSDF) in New Mexico, USA is selected. WSDF is a gypsum dune field with an area of about 500 km 2 in the Tularosa Basin between the San Andres Mountains and the Sacramento Mountains ( Figure 2). The deflation of evaporite beds of Lake Otero and other playa lakes in the basin is the major sediment source of WSDF [48]. The selected study area has three types of dunes-barchans, transverse, and parabolic-with up to 15 m in height. The general trend of the crestlines is NNW (345 • ), and the average spacing of the crestlines is 136 m (Ewing et al. 2006). The dominant winds from the SW-W are strongest during the winter-spring, while winds from the N-NW occur during the fall-winter, and winds from the S-SE occur during the spring-summer [48]. WSDF has been well documented in numerous studies [41,42,[48][49][50][51]. To overcome the above major limitation of the PSTP method, this paper proposes a new method called toe line tracking (TLT) for the automated measurement of sand dune migration distance and direction using multi-temporal DEMs derived from LiDAR data in a study area in the White Sands Dune Field (WSDF) in New Mexico, USA. Compared with the PSTP method which relies on the angle of repose for sand dune, the new method aims to extract toe lines of sand dunes based on several simple parameters, and conduct automated measurement of dune migration rates and directions through random sampling and spatial analysis in geographic information systems (GIS). The study area, data, methods, results, discussion, and conclusions are presented in the following sections.

Study Area
A study area of 9 × 2.4 km in the White Sands Dune Field (WSDF) in New Mexico, USA is selected. WSDF is a gypsum dune field with an area of about 500 km 2 in the Tularosa Basin between the San Andres Mountains and the Sacramento Mountains ( Figure 2). The deflation of evaporite beds of Lake Otero and other playa lakes in the basin is the major sediment source of WSDF [48]. The selected study area has three types of dunesbarchans, transverse, and parabolic-with up to 15 m in height. The general trend of the crestlines is NNW (345°), and the average spacing of the crestlines is 136 m (Ewing et al. 2006). The dominant winds from the SW-W are strongest during the winter-spring, while winds from the N-NW occur during the fall-winter, and winds from the S-SE occur during the spring-summer [48]. WSDF has been well documented in numerous studies [41,42,[48][49][50][51].

Data
LiDAR point clouds acquired on 24 January 2009; 26 September 2009; and 6 June 2010 in the study area were downloaded from the OpenTopography Facility (www.Open-Topography.org, accessed on 2 March 2015) at the San Diego Supercomputer Center. Li-DAR data acquisition and processing was completed by the National Center for Airborne Laser Mapping (NCALM-http://www.ncalm.org (accessed date 2 August 2021)). The point densities of the three datasets are 4.19 points/m 2 , 5.63 points/m 2 , and 4.62 points/m 2 , respectively. The horizontal coordinate system is NAD83 UTM Zone 13N, and the vertical coordinate system is NAVD88. The vertical accuracy of the datasets are less than 0.35 m. Figure 3 shows the digital elevation models (DEM) created from the multi-temporal Li-DAR point clouds using Inverse Distance Weighted (IDW) interpolation with a cell size of 1 × 1 m. The IDW interpolation method was selected because (1) the LiDAR point densities are high enough for a cell size of 1 × 1 m; and (2) the results in this paper will be

Data
LiDAR point clouds acquired on 24 January 2009; 26 September 2009; and 6 June 2010 in the study area were downloaded from the OpenTopography Facility (www. OpenTopography.org, accessed on 2 March 2015) at the San Diego Supercomputer Center. LiDAR data acquisition and processing was completed by the National Center for Airborne Laser Mapping (NCALM-http://www.ncalm.org (accessed date 2 August 2021)). The point densities of the three datasets are 4.19 points/m 2 , 5.63 points/m 2 , and 4.62 points/m 2 , respectively. The horizontal coordinate system is NAD83 UTM Zone 13N, and the vertical coordinate system is NAVD88. The vertical accuracy of the datasets are less than 0.35 m. Figure 3 shows the digital elevation models (DEM) created from the multi-temporal Li-DAR point clouds using Inverse Distance Weighted (IDW) interpolation with a cell size of 1 × 1 m. The IDW interpolation method was selected because (1) the LiDAR point densities are high enough for a cell size of 1 × 1 m; and (2) the results in this paper will be compared with the results from Dong [44] which was based on IDW interpolation of LiDAR points. compared with the results from Dong [44] which was based on IDW interpolation of Li-DAR points.

Methodology Flowchart
The methodology flowchart for toe line tracking (TLT) is shown in Figure 4. The process is implemented using the Python programming language and ArcPy for ArcGIS 10.7 with ArcGIS Desktop Advanced License. The user input parameters include workspace (folder for output files), area of interest (AOI, optional), first DEM (D1), second DEM (D2), window size (w), toe height (h), area threshold (a), sample rate (d), and search radius (r). With these initial parameters, the whole process in the flowchart is executed automatically. The flowchart can be implemented in other GIS packages. More details of the major steps in the flowchart are explained in the following sections.

Methodology Flowchart
The methodology flowchart for toe line tracking (TLT) is shown in Figure 4. The process is implemented using the Python programming language and ArcPy for ArcGIS 10.7 with ArcGIS Desktop Advanced License. The user input parameters include workspace (folder for output files), area of interest (AOI, optional), first DEM (D 1 ), second DEM (D 2 ), window size (w), toe height (h), area threshold (a), sample rate (d), and search radius (r). With these initial parameters, the whole process in the flowchart is executed automatically. The flowchart can be implemented in other GIS packages. More details of the major steps in the flowchart are explained in the following sections.

Extraction of Base Level and Toe Lines
As shown in Figure 1, reverse winds or crosswinds can modify the dune shape, and slipfaces may not be detected in DEMs using the angle of repose. The reverse winds or crosswinds can also cause sand to accumulate at the base of the dune, creating a concaveup apron ( Figure 5). The point between the apron surface and the slipface (which may be modified by wind) is the toe, the area between dunes is the interdune, and the base level is a surface created from local minimum elevations of interdunes ( Figure 5). For multitemporal sand dune DEMs with or without slipfaces, the distance between two toes should provide a better approximation for the migration distance, because the dune crestline and the upper part of the leeward slope are more likely to change under reverse winds or crosswinds, while defining an identifiable point on the concave-up apron surface can be difficult. Therefore, a method is developed to extract the base level and toe lines based on the elevation of interdunes. For an input DEM D(x, y) for sand dunes, the base level raster M(x, y) can be calculated using the minimum elevation in a w by w square window that moves across D(x, y). The user-specified window size w is determined by measuring the approximate distance

Extraction of Base Level and Toe Lines
As shown in Figure 1, reverse winds or crosswinds can modify the dune shape, and slipfaces may not be detected in DEMs using the angle of repose. The reverse winds or crosswinds can also cause sand to accumulate at the base of the dune, creating a concave-up apron ( Figure 5). The point between the apron surface and the slipface (which may be modified by wind) is the toe, the area between dunes is the interdune, and the base level is a surface created from local minimum elevations of interdunes ( Figure 5). For multitemporal sand dune DEMs with or without slipfaces, the distance between two toes should provide a better approximation for the migration distance, because the dune crestline and the upper part of the leeward slope are more likely to change under reverse winds or crosswinds, while defining an identifiable point on the concave-up apron surface can be difficult. Therefore, a method is developed to extract the base level and toe lines based on the elevation of interdunes.

Extraction of Base Level and Toe Lines
As shown in Figure 1, reverse winds or crosswinds can modify the dune shape, and slipfaces may not be detected in DEMs using the angle of repose. The reverse winds or crosswinds can also cause sand to accumulate at the base of the dune, creating a concaveup apron ( Figure 5). The point between the apron surface and the slipface (which may be modified by wind) is the toe, the area between dunes is the interdune, and the base level is a surface created from local minimum elevations of interdunes ( Figure 5). For multitemporal sand dune DEMs with or without slipfaces, the distance between two toes should provide a better approximation for the migration distance, because the dune crestline and the upper part of the leeward slope are more likely to change under reverse winds or crosswinds, while defining an identifiable point on the concave-up apron surface can be difficult. Therefore, a method is developed to extract the base level and toe lines based on the elevation of interdunes. For an input DEM D(x, y) for sand dunes, the base level raster M(x, y) can be calculated using the minimum elevation in a w by w square window that moves across D(x, y). The user-specified window size w is determined by measuring the approximate distance For an input DEM D(x, y) for sand dunes, the base level raster M(x, y) can be calculated using the minimum elevation in a w by w square window that moves across D(x, y). The user-specified window size w is determined by measuring the approximate distance between two dune crestlines ( Figure 6). This is to ensure that, for any location in D(x, y), a minimum elevation in a neighboring interdune area can be used as output in M(x, y). An elevated base level raster E(x, y) is obtained by adding a constant toe height (h) to M(x, y) between two dune crestlines ( Figure 6). This is to ensure that, for any location in D(x, y), a minimum elevation in a neighboring interdune area can be used as output in M(x, y). An elevated base level raster E(x, y) is obtained by adding a constant toe height (h) to M(x, y) The toe height h is a user-specified constant which may vary with different dune areas. However, it is not necessary to know the exact toe height in a specific dune field. It is estimated that a default h value of 1 m should work for many sand dunes. Discussion on the sensitivity of the toe height h can be found in the discussion section. Like the window size w, h can be determined by examining the 3D dune profiles. The elevated base level raster E(x, y) can be converted to a binary raster B(x, y) showing dunes (0's) and interdunes (1's) by comparison with D(x, y) As shown in the flowchart (Figure 4), the binary raster B(x, y) is then converted to polygon features, and the polygons are smoothed and converted from multipart to single part to represent the toe lines. Sampling points can be generated on toe lines for the automated measurement of dune migration distance and direction ( Figure 6), which is explained in the next section.

Calculation of Migration Distance and Direction
Once the toe line polygon feature classes G1 and G2 are created from two DEMs (Figure 7A), two logical operators, "AND" and "XOR", and spatial relationships between the polygons can be used to derive source toe lines from the first DEM, and target toe lines The toe height h is a user-specified constant which may vary with different dune areas. However, it is not necessary to know the exact toe height in a specific dune field. It is estimated that a default h value of 1 m should work for many sand dunes. Discussion on the sensitivity of the toe height h can be found in the discussion section. Like the window size w, h can be determined by examining the 3D dune profiles.
The elevated base level raster E(x, y) can be converted to a binary raster B(x, y) showing dunes (0's) and interdunes (1's) by comparison with D(x, y) As shown in the flowchart (Figure 4), the binary raster B(x, y) is then converted to polygon features, and the polygons are smoothed and converted from multipart to single part to represent the toe lines. Sampling points can be generated on toe lines for the automated measurement of dune migration distance and direction ( Figure 6), which is explained in the next section.

Calculation of Migration Distance and Direction
Once the toe line polygon feature classes G 1 and G 2 are created from two DEMs ( Figure 7A), two logical operators, "AND" and "XOR", and spatial relationships between the polygons can be used to derive source toe lines from the first DEM, and target toe lines from the second DEM. "XOR" is an exclusive "or" (or and only or) operator to produce the difference between the two polygons, while "AND" produces the intersection between the two polygons. The difference polygons K (white areas in Figure 7B) are created using x FOR PEER REVIEW 8 of 18

Sensitivity Analysis
Among the five user-defined parameters in the proposed method (namely window size, w; toe height, h; area threshold, a; sample rate, d; and search radius, r), the last three parameters are based on extracted toe lines, and therefore do not have major impact on the output. Three different window sizes (w) and three different toe heights (h) are tested for sensitivity analysis, and the details are presented in 5.1. The multipart polygon feature class K is then converted to a single part polygon feature class L. The following three types of polygons are then removed from L: (1) small polygons with an area less than the area threshold a; (2) polygons that are not within G 1 ; and (3) polygons that do not intersect G 2 . It is necessary to remove polygons that do not intersect G 2 , because Equation (3) may produce polygons that are part of G 1 but do not intersect any polygon in G 2 , which may result in incorrect measurements of dune migration distance and direction. The white areas in Figure 7C show polygons in L.
The final source toe lines S (dashed lines in Figure 7D) and target toe lines T (solid lines in Figure 7D) are created by converting multipart lines in U and V to single part lines, respectively. Sampling points are generated on polylines in S based on the sample rate d. The distance from each sampling point to the nearest point on the target line is found using the search radius r. In the study by Yao et al. (2007), migration directions were perpendicular to the dune ridges. Similarly, migration directions in this study were calculated from the sampling points to the nearest points on the target lines. The directions are converted to the range of 0 • (north) to 360 • . If no nearest point is found, the distance value will be-1. Figure 7E,F show sample dune migration distances (meters) and directions (degrees) automatically calculated for sampling points on source lines. Since the number of days between the two DEMs is known, the migration distances can be further converted to migration rates (meters per year).

Sensitivity Analysis
Among the five user-defined parameters in the proposed method (namely window size, w; toe height, h; area threshold, a; sample rate, d; and search radius, r), the last three parameters are based on extracted toe lines, and therefore do not have major impact on the output. Three different window sizes (w) and three different toe heights (h) are tested for sensitivity analysis, and the details are presented in 5.1.

Migration Distance/Rate
The window size (w) can be estimated by measuring distances between adjacent dune ridges. Usually, the mean value of several measurements can be used as window size. The window size is a rough estimate, and does not have to be very accurate (see 5.1 for more discussion on this). Using the basic measurement tools such as distance measurement and profile generation in ArcGIS 10.7, the following parameters were estimated and employed for calculating dune migration distance and direction: w = 121 m, h = 1 m, a = 20 m 2 , d = 20 m, and r = 30 m. Figure 11 Figure 11C). Figure 11D shows the toe line change from 24 January 2009 to 26 September 2009 to 6 June 2010. Similar to Figures 9 and 10, it can be seen that the toe lines match the dune landforms very well. m/year in the study area, very similar to the results calculated from the slipfaces of 24 January 2009 and 6 June 2010 using the PSTP method [44]. The migration distances in Figure 11 can be further converted to continuous rasters for migration rates using spatial interpolation. Figure 13 is a raster map showing the dune migration rates in three periods: 24 January 2009 to 26 September 2009, 26 September 2009 to 6 June 2010, and 24 January 2009 to 6 June 2010. The rasters were created using Kriging interpolation of the migration rates for the sampling points in the study area, with a search radius of 300 m and an output cell size of 120 m. The results are very similar to the raster for migration rates created by Dong [44]: high migration rates mainly concentrate on the upwind, western edge of the dune field, whereas low migrations rates occur in the vegetated parabolic dune area which is relative stable in the east part of the study area. It seems that there are no significant changes in the spatial patterns of dune migration rates in the above three periods, although some locations may show slightly higher migration rates in the period of 26 September 2009 to 6 June 2010 compared with the other two periods. Of particular interest is the two areas with relatively low migration rates in the central part of the study area, which needs further investigation. The migration distances in Figure 11 can be further converted to continuous rasters for migration rates using spatial interpolation. Figure 13 is a raster map showing the dune migration rates in three periods: 24 January 2009 to 26 September 2009, 26 September 2009 to 6 June 2010, and 24 January 2009 to 6 June 2010. The rasters were created using Kriging interpolation of the migration rates for the sampling points in the study area, with a search radius of 300 m and an output cell size of 120 m. The results are very similar to the raster for migration rates created by Dong [44]: high migration rates mainly concentrate on the upwind, western edge of the dune field, whereas low migrations rates occur in the vegetated parabolic dune area which is relative stable in the east part of the study area. It seems that there are no significant changes in the spatial patterns of dune migration rates in the above three periods, although some locations may show slightly higher migration rates in the period of 26 September 2009 to 6 June 2010 compared with the other two periods. Of particular interest is the two areas with relatively low migration rates in the central part of the study area, which needs further investigation. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 18  Figure 14 shows the histograms of the source directions (from the source lines to the target lines), where 0° is the north. It can be seen that the migration of the toe lines is mainly from SW-W to SE-E.

Migration Direction
The mean migration directions were calculated for three time periods: (1)  ). However, it should be noted these mean migration directions are calculated from the migration directions of the toe lines, and do not represent the mean wind direction. In fact, reverse winds can change the dune crestlines but do not change the toe lines in the reverse direction in most cases (also see Figure 1). The relationship between the dune migration direction and the wind directions can be more complex. For example, the narrow intervals between the toe lines at ①, ②, and ③ in Figure 11D could be caused by wind flows not perpendicular to the migration direction, as discussed in Tsoar et al. [52].  Figure 14 shows the histograms of the source directions (from the source lines to the target lines), where 0 • is the north. It can be seen that the migration of the toe lines is mainly from SW-W to SE-E.

Migration Direction
The mean migration directions were calculated for three time periods: (1)  ). However, it should be noted these mean migration directions are calculated from the migration directions of the toe lines, and do not represent the mean wind direction. In fact, reverse winds can change the dune crestlines but do not change the toe lines in the reverse direction in most cases (also see Figure 1). The relationship between the dune migration direction and the wind directions can be more complex. For example, the narrow intervals between the toe lines at 1 , 2 , and 3 in Figure 11D could be caused by wind flows not perpendicular to the migration direction, as discussed in Tsoar et al. [52].

Sensitivity Analysis
As explained in 3.1, there are five user-defined parameters in the proposed window size (w), toe height (h), area threshold (a), sample rate (d), and search ra The last three parameters do not have major impact on the output because they a on the extracted toe lines. However, window size (w) and toe height (h) may a binary rasters (B1 and B2 in Figure 4 Figure 15 shows interdune areas derived from the 24 January 2009 and 6 J DEMs using the TLT method, with h = 1.0 m and w = 106 m, 121 m, and 136 m. Th toe lines (dashed lines) and target toe lines (solid lines) obtained using w = 121 m shown as reference. When the window size w changes from 106 m to 121 m to can be seen that identical toe lines can be obtained using the TLT method. There minor changes in the raster cells along the borders of interdune areas, which are s during raster to polygon conversion. Such results are not surprising, because

Sensitivity Analysis
As explained in 3.1, there are five user-defined parameters in the proposed method: window size (w), toe height (h), area threshold (a), sample rate (d), and search radius (r). The last three parameters do not have major impact on the output because they are based on the extracted toe lines. However, window size (w) and toe height (h) may affect the binary rasters (B 1 and B 2 in Figure 4 Figure 15 shows interdune areas derived from the 24 January 2009 and 6 June 2010 DEMs using the TLT method, with h = 1.0 m and w = 106 m, 121 m, and 136 m. The source toe lines (dashed lines) and target toe lines (solid lines) obtained using w = 121 m are also shown as reference. When the window size w changes from 106 m to 121 m to 136 m. It can be seen that identical toe lines can be obtained using the TLT method. There are only minor changes in the raster cells along the borders of interdune areas, which are smoothed during raster to polygon conversion. Such results are not surprising, because the base level rasters are obtained by calculating the minimum elevations in the w × w windows, following by calculating the mean values in the w × w windows. As a result, different window sizes, such as 106 × 106, 121 × 121, and 136 × 136, may produce very close base level rasters. To show the similarities between the base levels generated using different window sizes, 384 random points were generated in the study area to extract elevations from the three base level rasters. Figure 16 shows the scatter diagram between the base level rasters obtained using w = 106 m and 136 m, respectively. It can be seen that the two base level rasters are almost identical. The results suggest that an estimated average distance between sand dune ridges can be used as the window size, and that average distance does not have to be very accurate. from the three base level rasters. Figure 16 shows the scatter diagram between the base level rasters obtained using w = 106 m and 136 m, respectively. It can be seen that the two base level rasters are almost identical. The results suggest that an estimated average distance between sand dune ridges can be used as the window size, and that average distance does not have to be very accurate.    from the three base level rasters. Figure 16 shows the scatter diagram between the base level rasters obtained using w = 106 m and 136 m, respectively. It can be seen that the two base level rasters are almost identical. The results suggest that an estimated average distance between sand dune ridges can be used as the window size, and that average distance does not have to be very accurate.     To get a better picture of the toe line variations with different toe heights in the whole study area, the following two tests were conducted: (1) A total of 2670 random sampling points were generated on source toe lines (h = 1.0 m), and a radius of 3 m around each sampling point was used to search the nearest point on the source toe lines (h = 1.5 m) and measure the distance (DIST1). (2) Similarly, a total of 2297 random sampling points were generated on target toe lines (h = 1.0 m) to measure the distances (DIST2) from each point the nearest target toe lines (h = 1. 5 m). In other words, random points were generated to measure distances between the two dashed lines (DIST1) and between the two solid lines (DIST2) in Figure 16 over the whole study area. The results are summarized in Table 1.  Table 1 shows that almost 90% of toe line points are within 2 m when toe height changes from 1.0 m to 1.5 m, while a majority of them (> 60%) are within 1 m. It should be noted that such displacements in toe lines caused by toe height change do not necessarily mean changes in migration distances, because both the source toe line and the target toe line usually change in the same direction as shown in Figure 17, maintaining a relatively stable distance. The results suggest that the toe height (h) change does not significantly affect the dune migration measurements, as long as h is not very low (close to the vertical accuracy of LiDAR data) and not very high (for example, less than 1.5 m so that small dunes can also be measured). A toe height of 1 m should work for measurement of many dune fields.

Limitations
Like the PSTP method, the application of the TLT method can be affected by dune To get a better picture of the toe line variations with different toe heights in the whole study area, the following two tests were conducted: (1) A total of 2670 random sampling points were generated on source toe lines (h = 1.0 m), and a radius of 3 m around each sampling point was used to search the nearest point on the source toe lines (h = 1.5 m) and measure the distance (DIST1). (2) Similarly, a total of 2297 random sampling points were generated on target toe lines (h = 1.0 m) to measure the distances (DIST2) from each point the nearest target toe lines (h = 1.5 m). In other words, random points were generated to measure distances between the two dashed lines (DIST1) and between the two solid lines (DIST2) in Figure 16 over the whole study area. The results are summarized in Table 1.  Table 1 shows that almost 90% of toe line points are within 2 m when toe height changes from 1.0 m to 1.5 m, while a majority of them (>60%) are within 1 m. It should be noted that such displacements in toe lines caused by toe height change do not necessarily mean changes in migration distances, because both the source toe line and the target toe line usually change in the same direction as shown in Figure 17, maintaining a relatively stable distance. The results suggest that the toe height (h) change does not significantly affect the dune migration measurements, as long as h is not very low (close to the vertical accuracy of LiDAR data) and not very high (for example, less than 1.5 m so that small dunes can also be measured). A toe height of 1 m should work for measurement of many dune fields.

Limitations
Like the PSTP method, the application of the TLT method can be affected by dune form, dune interaction, dune migration rate, and data acquisition interval. The TLT method can be used for barchan, barchanoid ridge, transverse dunes, and even linear dunes, but cannot be directly applied to star dunes where the complex wind regime creates multiple slopes. Depending on the migration rates of a dune field, the TLT method may not be suitable for measuring dune migration distance and direction if two datasets are acquired with a substantial time interval (e.g., several years) for a fast-moving dune field, because dune interaction during migration may involve complex changes in dune shape. It is hoped that terrestrial (stationary or mobile) and unmanned aerial vehicle (UAV) LiDAR systems can provide new options for monitoring sand dune changes.

Conclusions
Automated or semi-automated change detection of large dune fields is a challenging problem. After reviewing the limitation of the existing PSTP method for measuring sand dune migration rate and migration direction, this paper introduced a toe line tracking (TLT) method and implemented the new method using the Python programming language for automated measurement of dune migration distance and direction based on multi-temporal LiDAR-derived DEMs. Using a few simple parameters, the TLT method automatically extracts the base level of a dune field and toe lines of individual dunes. The toe line polygons derived from two DEMs are processed using logical operators and other spatial analysis methods in a geographic information system, and source toe lines and target toe lines are extracted from the toe line polygons.
The application of the TLT method was demonstrated using multi-temporal LiDARderived DEMs for a 9 × 2.4 km area in the White Sands Dune Field in New Mexico (USA). By generating over 7000 random sampling points along source toe lines, dune migration distances and directions are calculated and analyzed. Dune migration distances and directions for three periods ( Sensitivity analyses were carried out using different window sizes (106, 121, and 136 m) and toe heights (0.5, 1.0, and 1.5 m). The results suggest that the TLT method is relatively insensitive to changes in window size. Although a toe height of 0.5 m (which is close to the vertical accuracy of LiDAR data, 0.35 m, in this study) does not work well in some areas due to elevation variations near the interdune surfaces and presence of shrubs, a toe height between 1.0 to 1.5 m can produce very close outputs, as shown in statistics from thousands of sampling points. The results show that TLT is a more generic method than PSTP in that it works for dunes with or without slipfaces that reach the angle of repose, and can be used for effective measurement of dune migration distance/rate and direction in large dune fields.