1. Introduction
Research on salt flats dynamics has gained much importance in recent years mainly due to environmental and mining issues. The worldwide increasing demand for lithium has put the focus on salt flats, which are the main source of this element, considered strategic and fundamental in the transition to an economy less dependent on fossil fuels [
1,
2,
3,
4,
5]. The construction of rechargeable, lightweight, lithium-based batteries with a high storage capacity and its demand in numerous green technologies, such as electric vehicles, turns lithium into a fundamental resource nowadays [
2,
6,
7,
8]. The Andean highlands shared by Argentina, Bolivia, and Chile concentrates about 60% of the lithium reserves of the planet within an area known as the “Lithium Triangle”. From a strategic and economic standpoint, surveys for lithium exploration in the Andean highlands are complex since they are remote, difficult-to-access areas, leading to time-consuming and costly field campaigns. The valuable capabilities of remote sensors help circumvent these issues, enabling information at different spatial and time scales through spaceborne platforms, thus avoiding a large number of field campaigns [
9,
10]. While optical sensors make use of the reflection properties of targets, radar sensors can obtain information with a certain degree of penetration under the Earth’s surface [
11,
12].
Among the different types of radars, the Synthetic Aperture Radar (SAR) highlights its spatial resolution of the order of a few meters and its all-weather capabilities. SARs actively transmit waves onto a linear basis with a single wavelength, being the most common free-space wavelengths for land applications between 3.8 and 7.5 cm (C-band) and between 15 and 30 cm (L-band). The larger the wavelength emitted by the radar, the greater the penetration on the land surface [
11,
12,
13,
14,
15], thus potentially allowing subsurface exploration. Besides observation geometry and acquisition configuration, the radar response depends on interface configuration such as roughness, bulk media properties such as composition, dielectric constant, and the condition of pore saturation. Overall dielectric properties of saline soils and related implications to radar remote sensing at C- and L-band were the subject of several studies [
16,
17,
18].
Dielectric properties of saline water at concentrations up to 188 g/L (grams of solid salt dissolved in one liter of solution) are successfully modeled by a single Debye relationship combined with an ionic conductivity term [
11,
19]. However, when dealing with higher salt concentrations such as those found in brines, which may exceed the aforementioned limit, a dedicated semiempirical model that extends the validity range is formulated in [
20]. The model relies on a Cole–Cole relationship for complex permittivity of sodium chloride (NaCl) solutions from pure water to salt concentrations up to 292 g/L at microwave frequencies up to 20 GHz and in the range 5 °C to 35 °C.
Most investigations with spaceborne radars focused on salt flats are two-fold. On one hand, “snapshot” studies exploit a few radar acquisitions combined with simultaneous field surveys to interpret radar responses. Pioneering research using satellite microwave active and passive acquisitions of the Utah Great Salt Lake Desert suggested the influence of subsurface layers of sediment saturated with brine on the observations [
21]. A study over Lop Nur Lake Basin (China) demonstrated the capabilities of SAR images to recognize different types of salt crusts by surface roughness parameters using Polarimetric Synthetic Aperture Radar (PolSAR) with ALOS-PALSAR images [
22]. Using field measurements and co-polar signatures derived from analytical simulations, the influence of soil salinity as a function of soil moisture on the dielectric constant was assessed for airborne and spaceborne acquisitions at C- and L-bands in a salt pan in Death Valley, California (USA) [
17]. In this same respect, radar responses of salt-affected soils were modeled for spaceborne imagery at the C-band over a salty depression located in the Egyptian desert [
23]. Existing methods that rely on evaluating scattering models require certain input parameters, which, in turn, imply in situ fieldwork. Extrapolation of published datasets might serve as an alternative to the fieldwork, although their availability on salt pans might be limited. Also, the validity range of the scattering model might narrow its applicability. However, scattering models provide a complete characterization of the salt pan by radar waves.
On the other hand, multi-temporal studies make use of long-term, dense sets of radar acquisitions to monitor salt flat dynamics, usually with the only ancillary information provided by weather stations and visual information collected on a field visit. As soil moisture largely affects dielectric properties of soils, multi-temporal studies become important to analyze seasonal variations of the SAR response on salt pans, as was demonstrated in Chott El Djerid playa deposit in Tunisia using multi-temporal RADARSAT-2 C-band full polarimetric imagery [
24,
25] and in Salar de Aguas Calientes Sur in Chile using Sentinel-1 and ALOS-2/PALSAR-2 [
26]. These investigations contributed to a better understanding of crust salt flat dynamics and related evaporitic processes. The methods using time series strongly depend on the availability of satellite platforms and their continuity on image provision. To date, few satellite missions fulfilled this constraint. On the other hand, insights gained by multi-temporal studies contribute to a comprehensive understanding of salt pans.
Penetration capabilities of radar waves were also exploited in a few research studies using fully polarimetric images. In the same Lop Nur region mentioned above, a study found some sensitivity in retrieving the depth of the subsurface brine layer with field measurements using co-polarized phase difference in ALOS-PALSAR L-band images [
27]. In 2004, using airborne data over Pyla Dune in Arcachon Basin close to Bordeaux, France, subsurface moisture information related to wet structures (paleosoils) was detected using phase signature of polarimetric L-band SAR data [
28]. Conclusions drawn from these studies and from [
21] highlight the relevance of modeling subsurface layers of the salt pan, where the pore saturation with high-salinity brines is a unique feature of these targets.
In understanding the radar response, electromagnetic models can make a difference by simulating the backscattered power, thus allowing a composite analysis of the different parameters involved in the characterization of the target. One of them is the Small Slope Approximation (SSA), which provides a solution for wave scattering both at small and large scales provided that surface roughness has small slopes (i.e., the ratio of vertical to horizontal scales is smaller than some function of the wavelength). It is shown that the obtained approximation offers a unified approach to wave scattering problems by combining perturbation theory with the tangent plane approximation, thus bridging the gap between the classical approaches small perturbation method and the Kirchhoff approximation [
29,
30]. A limitation of this model is that it will only be applicable to the crusts with the smaller slopes, usually the halite and the smoother earthy crusts. As a counterpart, the model allows a valuable, coherent, two-layer description of two half-spaces and a third in-between layer in which the top and below interfaces are rough. In this work, numerical simulations from second-order, two-layered SSA, as described in [
31,
32], were used to model the salt pan configuration.
Baseline information for a deeper understanding of crust dynamics is collected through a five-year multi-temporal dataset (January 2018–May 2023) with a dense time series of Sentinel-1 (C-band), ALOS-2/PALSAR-2, and SAOCOM-1 (L-band). Ancillary information from an in situ weather station and from a distributed network of bores aided the interpretation of crust formation and later development. In addition, in situ parameters, such as crust roughness, horizon configuration, water table depth, and brine salinity, were collected at sampling sites across the salt pan by means of a field campaign carried out on 28 and 29 May 2023.
This study is aimed at assessing the capabilities of a scattering model in relating the SAR backscattered signal with the surface and subsurface salt pan configuration of the Pozuelos salt flat located at Puna Salteña in northern Argentina.
The remainder of this paper is organized as follows. The study area and the fieldwork are presented in
Section 2.1 and
Section 2.2, respectively. Detailed information on SAR image processing, crust classification, backscattering models, and the Bayesian approach for estimating model parameters are presented in
Section 2.3 to
Section 2.7. The influence of weather on SAR backscattering over time is described in
Section 3.1 and
Section 3.2. Surface and subsurface modeling are analyzed in
Section 3.3 and
Section 3.4. Discussion is provided in
Section 4, followed by concluding remarks in
Section 5.
3. Results
3.1. Water Dynamics after Heavy Rainfalls
After heavy rainfall events, the salt pan flooded, except for the depocenter, and gradually dried up in the following days. Water remained longer over the earthy crusts than over the halite ones, the latter being covered with a waterbed, and displayed in bluish hues north-, east- and centerward (type I and V crusts) in
Figure 8. For those crusts with a large porous pattern, multiple depolarizations increase the backscattering coefficient at cross-polarization when they remain partially filled with water.
3.2. Time Series Analysis
The SAR backscatter dynamics are depicted in
Figure 9 and
Figure 10, where the average backscattering coefficients for each crust type and the 25th- and 75th-percentile were plotted. Rainfall was also plotted when available. For type II (orange), III (blue), and IV (green) earthy crusts, C-band VV-polarized backscattering coefficient temporal evolution depicted an annual periodic pattern triggered by the summertime rainfalls in January and then an overall decrease until the end of the year. The type I halite crust (red) exhibited the opposite dynamics, with a marked drop at the first heavy rainfall around January and then a steady increase until the end of the year. Among the earthy crusts, type IV had larger backscattering coefficients, followed by those of type III and type II, respectively. The aforementioned trend can be explained with the surface roughness slope summarized in
Table 1, where
s/
l~0.4–0.6 for type IV and
s/
l~0.2 for type III. Despite the type II crust having a large slope (
s/
l~0.45), the C-band scatter relies more on the underlying, moderately rough surface than on the protrusions, which are many wavelengths in size at the C-band.
For the cross-polarized case (
Figure 10), the annual pattern is similar to the co-polarized one, although type II and III crusts exhibited similar VH-polarized backscattering coefficients, opposite to the co-polarized case where at least a 2 dB difference between those crusts was observed. Also, the larger differences between these crusts occurred close to the rainfall events, when the voids below the bosses in type III crust filled with water, thus enhancing the volume scattering mechanism [
11]. In this same respect, the large porous pattern in type IV crust seemed to be responsible for the large backscattering coefficient at VH.
Backscatter dynamics at the L-band are depicted in
Figure 9 and
Figure 10 with colored triangle and diamond markers. Overall, type I crust is generally lower than the corresponding at the C-band for both the co- and cross-polarized modes. In contrast, type II (orange triangles) and type IV (green triangles) crusts exhibit the opposite, with the co-polarized L-band backscattering coefficients for ALOS-2/PALSAR-2 exceeding those at the C-band as well as most of the SAOCOM-1 ones. For type III, the distinction is less pronounced, at least in the case of ALOS-2/PALSAR-2 data.
In some months, namely August 2019, January 2020, and April 2021, ALOS-2/PALSAR-2 backscattering coefficient for type II (orange triangles) is larger than that of type IV (green triangles). The same is true in some months for SAOCOM-1 (diamond markers). This happened mostly for dried-up crusts and could be related to subsurface scattering mechanisms enhanced by a shallow water table at the depocenter in those months.
At cross-polarization, the backscattering coefficient for types II and IV at the L-band is consistently larger than that at the C-band as in the case of co-polarization. However, the ordering of crust types did not reverse and followed that at the C-band: σ0IV > σ0II > σ0I, and σ0II is usually larger than σ0III. For type III crust, the SAOCOM-1 backscattering coefficients are lower than those observed with Sentinel-1, with one exception on 5 January 2023 where the different orbit pass may explain that VH at the C-band is less than HV at the L-band.
Comparing the temporal curves obtained for the different SAR images, it can be observed that the classification of the crusts is better defined in the co-polarized data (
Figure 9) than in the cross-polarization data (
Figure 10), where temporal curves intersect for type II and type III crusts, clearly evident in the year 2018.
The backscatter dynamics of the different crusts are driven by the precipitation events, which are very scarce (10 to 50 mm/year in Puna Salteña and 50 to 80 mm/year in Pozuelos Salt Flat) but concentrated in the summertime from December to February [
33,
44,
45], leading to high rates of rainfall that affect salt crusts differently. Type I and V crusts experience a decrease in their backscattering coefficient, reaching a minimum. Water droplets cause a disruption in the crystalline structure of the halite. When exposed to rainfall, the crust smoothens as the crystalline structure dissolves, gradually recrystallizing and increasing its roughness over time.
Figure 11 displays backscattering coefficients for Sentinel-1 ascending (VV and VH) and ALOS-2/PALSAR-2 ascending (HH and HV), averaged over a circular area with a 50 m radius at each sampling site, and the corresponding depth of the water table. SAR images and in situ measurements corresponded to the period of the fieldwork. It can be observed that for type I crust, values for the C-band are greater than those of the L-band. However, in earthy crusts (II, III, and IV), values for the C-band are either similar or lower. A subtle correlation of the backscattering coefficient with the depth of the water table is observed, decreasing as the water table becomes shallower. This is consistent with the large attenuation of the propagating waves into a highly lossy media due to the high-salinity brine. When the water table is very shallow, large dielectric contrast onto a very rough upper boundary, as in the case of type II crust, results in a large backscattering coefficient. The dependence of SAOCOM-1 backscattering coefficients on water table depth was very similar to the corresponding ALOS-2/PALSAR-2 shown in
Figure 11.
3.3. Upper Layer Roughness from the Two-Layer SSA Model
Considering that the C-band has less penetration depth into the soil than the L-band, a first attempt to assess the SSA model in retrieving surface parameters s1 and l1 on a type I (ID-1) and two type III (ID-7 and ID-16) crusts was performed. The SSA model was used, considering a layer thickness d given by the water table. With the aid of pictures taken on the soil profile at the trenches, overall observations of the water table inclusions on the soil led to lower layer parameters s2 = 2 × s1 and l2 = l1/2. A number of model simulations showed that backscattering coefficients have very low sensitivity to s2 and l2 variations.
Contour levels in
Figure 12 show the model computations (thick contours) for Sentinel-1 (blue) at VV and for SAOCOM-1 (black) at HH. The co-polarized measured backscattering coefficients (dotted contours) are from
Table 3. Intersection of measured contours resulted in the (
s1,
l1)-pair combination compatible with the SSA model and the spaceborne observations. With that intersection close to the measured in situ roughness parameters at the red crosses, this first attempt at assessing the SSA model yielded satisfactory outcomes. Similar results were found for the intersection of the cross-polarized backscattering coefficients as well as the combination of Sentinel-1 and ALOS-2/PALSAR-2. Although type I crust is better described by an exponential power spectrum, only a Gaussian one led to intersecting contours as shown in
Figure 12a.
3.4. Subsurface Parameter Estimation
With the SSA enabling the modeling of the backscatter radar response at dual-polarization (HH and HV) and a Bayesian inference scheme at hand, marginal and posterior distributions of the subsurface parameters
Sb and
d constrained to the radar observations can be computed to assess their compatibility with the measured ones at fieldwork.
Figure 13 shows the joint probability generated by sampling out of the posterior as blue contours, whereas the marginal distributions for
Sb and
d are in the diagonal, with their corresponding Kernel Density Estimation (KDE) indicated in light gray. The quartiles Q1 to Q3, each representing a fourth of the distributed sampled population, are shown as vertical dotted lines. The red plus marks indicate the in situ measurements. Sampling locations are ID-1, ID-7, and ID-16, where the upper panel corresponds to ALOS-2/PALSAR-2 and the lower panel to SAOCOM-1. The lowest contour drawn is 0.05, such that the integral over the area within is 0.95.
The extent of the posterior is overall related to the variance of the likelihood given by the uncertainties of the measured HH and HV. Thus, the more precise the measured backscattering coefficient, the less extended the posterior and, therefore, the more precise the parameter estimation. Two contrasting cases are given by ID-1 and ID-16, both for SAOCOM-1. While the posterior computed at ID-1 spanned a small area, the one corresponding to ID-16 spanned an area roughly four times bigger, in accordance with the HH and HV relative errors for these sites, as readily computed from
Table 3. The variance of the priors for
Sb and
d are given in [
2] and by a two-year-long record of water table depths, respectively, as stated in
Section 2.7.
Overall, the measured parameters are close to the median (Q2 quartile, black cross marker). Additionally, in those posteriors that are multimodal, they correspond to one of the maximums of probability. The sharp decrease of the KDE for the parameter d is in accordance with the bounds of the uniform distribution used prior.
The estimation of the salt concentration of brine is very precise since the modes of the probability are vertically aligned with almost no spread. Accuracy computed as (Sb_Q2 − Sb_insitu)/Sb_insitu ranges from 2% to 5%, with that of ID-16 (SAOCOM-1) being the poorer agreement. When replacing Q2 with the mode, accuracy slightly improves, ranging from 2% to 4%. Similar results have been found for the remaining sampling sites, with accuracies ranging from 1% to 8% for both Q2 and mode estimators. On the other hand, the estimation of water table depth has one or two elongated areas corresponding to the larger contour levels of the posterior distribution, each one compatible to some extent with the measured backscattering coefficients.
4. Discussion
The multi-temporal analysis of backscatter variations in SAR imagery over the Pozuelos salt flat reveals a strong relationship with precipitation data recorded at the on-site weather station (
Figure 9 and
Figure 10). A plausible explanation for this observation is the alteration of salt crust properties due to the interaction between their constituent minerals and liquid water droplets. The physicochemical characteristics of the diverse salt crusts are likely influenced by the presence of water, leading to changes in their backscattering behavior.
The backscatter dynamics of the different crusts were driven by the precipitation events, which are very scarce (10 to 50 mm/year in Puna Salteña and 50 to 80 mm/year in Pozuelos salt flat) but concentrated in the summertime from December to February [
32,
43,
44], leading to high rates of rainfall that affect salt crusts differently. Type I and V crusts, mostly composed of halite, underwent a dissolution process due to rainfall, resulting in a smoother surface, which implied a sharp decrease in the backscattering coefficient. Remarkably, observed rates of decrease are very similar across the different summertime periods in 2019–2022, ranging from 0.6 dB to 0.9 dB per 10 mm of daily rainfall. Subsequent drying led to the growth of salt crystals, which, in turn, increased the surface roughness and consequently raised the backscatter response [
25,
46].
In addition to the growth of salt crystals, the formation and later development of halite polygons led to an increase in the centimeter-scale roughness over time. This contribution to the overall surface roughness, as seen by the SAR, was not accounted for in the crust profiles with the gridded board, therefore explaining the mismatch found in the intersection of the SSA-based contour levels and the fieldwork estimate.
Conversely, earthy crusts (types II to IV) exhibited an increase in their microwave response, reaching a maximum on the first rainy days after a long dry period (
Figure 9 and
Figure 10). The rates of increase of the type III and IV crusts are about 0.3 dB to 1.2 dB per 10 mm of daily rainfall. A saturation value slightly above −5 dB at the C-band VV is consistently reached by the type IV crust. Although naturally rougher than halite crust, these crusts exhibited a greater resistance to rain due to their hardness. However, they also contained halite crystals in their composition, and some were covered with a halite layer, which dissolved upon contact with rainfall, resulting in a layer with increased roughness due to the formation of pores. Over time, these pores became refilled through the halite crystallization process [
47,
48]. In addition, large pores remain filled with water, increasing the dielectric constant of the surface prior to water evaporation or vertical runoff. The size and extent of the pores are also related to the amount of backscattering power at cross-polarization.
The effects of rainfall on salt crusts are discussed in [
48], where it is indicated that thick salt crusts with significant surface relief (>10 cm) are primarily formed by rainfall, which is consistent with the increase in the backscatter coefficient due to a marked increase in surface roughness observed in crust types II, III, and IV. A process known as efflorescence, by which mineral salts crystallize on the surface of a material when the water containing them evaporates, modifies the surface roughness afterward. This occurs naturally when brine penetrates and then evaporates, leaving the salts deposited on the surface, thus covering the underlying relief. Therefore, the surface smoothens and decreases its backscatter. Efflorescence can appear as a white or crystalline layer on the surface (e.g.,
Figure 6c).
Regarding the backscatter response at the L-band, type II crust backscattering coefficients were above type III, indicating that wavelength–size protrusions at the depocenter largely contributed to the backscatter. The lack of pores in type II prevented this behavior from being observed at HV polarization.
Remarkably, co-polarized backscattered power at the L-band was larger than that of the C-band for the type II and IV earthy crusts, indicating some overall large-scale roughness superimposed to a small-scale Gaussian power spectrum, despite the different incidence angles (34° for ALOS-2/PALSAR-2 and 41° for Sentinel-1). Type III did not display a clear trend, since it varied below and above the corresponding dynamics at the C-band. This seemed to be indicative of multiscale roughness. On the other hand, surface roughness on halite crusts seemed to be of single scale, since L-band backscattering coefficients were systematically lower with respect to that of the C-band (the normalized roughness k × s at the L-band is lower than at the C-band for the same roughness s, i.e., for single scale surfaces). The wavenumber k is defined as 2π/λ. Exponential and Gaussian power spectrums related to type I and type II to IV crusts, respectively, were indicative of the different physicochemical processes related to both the halite and earthy crust growing.
Crust classification with a dense time series at the C-band was better suited with the co-polarized backscatter response (
Figure 8) than that of the cross-polarized one (
Figure 9), where differentiation between type II and type III crusts may be misleading. A potential radar characterization of crusts was first introduced by [
26] as the monthly rate of increase of co-polarized backscattering coefficient at the C-band over large no-rainfall periods. For the type I crust analyzed in this work, the mentioned rate is around 0.2 dB/month, very different from the rates of 1–2 dB/month found in [
26] corresponding to a soft pan crust with thrust polygons, possibly indicating different local features such as brine concentration, water availability, and/or solar radiation patterns. An additional overall feature of crust dynamics might be the rate of variation around the summertime heavy rainfalls observed in the backscattering coefficient at the C-band, which seemed very uniform across the seasons for the halite crust. More research is needed in this regard for the earthy crusts.
Enhanced contrast of the backscattering coefficients among the different crust types is observed by comparing images acquired around heavy rainfall events, which are concentrated in austral summertime. Thus, the optimal period for radar characterization of highland salt flats is in late December and early March, on a yearly basis.
Inference of the subsurface model parameters brine salinity and average layer thickness were conducted on the posterior. The overall extent of the posterior and, therefore, the precision of the model parameters that are estimated with, depends primarily on the variance of the likelihood, which is ultimately related to the speckle filtering technique and the spatial homogeneity of the backscatter. In effect, backscatter heterogeneity hinders the efficiency of the speckle filter by adding variance other than the expected from the coherent illumination of the SAR. Thus, the heterogeneity of the backscatter itself is statistically different from the one that speckle filters were designed to deal with. To circumvent this shortcoming, sampling sites over suitable homogeneous areas of the crust surfaces were taken into consideration. The final uncertainties in the backscattering coefficients were similar to the 0.5 dB precision (stability) of the radar sensors ALOS-2/PALSAR-2 [
49] and SAOCOM-1 [
50].
While the salt concentration of subsurface brine was estimated precisely and with an accuracy better than 8% when the median (Q2) was considered, a multimodal or elongated posterior anticipated water table depth estimation to be more of a struggle. Besides dielectric loss driven by the salt concentration of the brine, clay material in the soil composition also reduces the penetration capabilities of SAR backscatter at the L-band [
47]. Therefore, when considering the water table, a weak correlation between the backscattering coefficients and the water table depths was expected (
Figure 11). Polarimetric observables, such as the co-polarized phase difference used to gain insight into the feasibility of water table depth estimation [
27], required full-polarimetric SAR images, which are largely less available than dual polarimetric images, such as those used in this paper for subsurface parameter retrieval.
In two-layered media, resonances of the radar wave occur as it repeatedly bounce between the media boundaries, with a strong dependency on the layered geometry under study [
51]. However, simulation studies on low-loss media showed that attenuation played a key role in fading the power returned for layer distances beyond half wavelength [
39]. Yet, the high loss of the salt pan prevented this resonance effect from happening even in the case of type I crust with the shallowest water table. Furthermore, this was also related to the lack of sensitivity of the backscattering power in changing the lower layer roughness found in the simulation study carried out in this paper.
5. Conclusions
Supported by a dedicated field campaign and the processing of 316 SAR images, the dynamics of a highland salt flat have been characterized as a yearly cycle closely linked to summertime rainfalls, which alter roughness configuration depending on the growth process and chemical composition of the crusts.
C-band co-polarized (VV) long-term backscatter response has proven to be effective in differentiating crusts, although a quantitative classification in this respect requires more research. On the other hand, cross-polarization (VH) seems a good proxy for pore patterns when radar imagery close to rainfall events is available.
This paper also demonstrated the potential for subsurface estimation with L-band dual-polarization images, constrained to crusts compatible with the feasibility range of the layered model at hand. Compared with available studies on salt flats with spaceborne radars, this research provided a new perspective on using microwave scattering modeling on dual polarimetric SAR data over salt flats and allowed for better exploitation of radar imagery from existing and upcoming satellite radar missions.
Coupling available multi-sensor SAR data with scattering models allows for the mapping of surface and subsurface configuration of a highland salt pan in a cost-efficient way. This integration could be beneficial for lithium exploration in these environments.