Assessment of new spectral indices and multi-seasonal ASTER data for gypsum mapping

Advanced spaceborne thermal emission and reflection radiometer (ASTER) data and different spectral indices were employed to map gypsum minerals. However, most proposed gypsum indices are designed based on the 2.21 μm gypsum absorption, overlapping with most hydroxy-bearing minerals. Moreover, the seasonal mutual transformation between gypsum, bassanite, and anhydrite may lead to seasonal reflectance variability of gypsum formation pixels, affecting the classification accuracy of gypsum indices. In this research, the feasibility of 2.26 μm (ASTER band7) reflectance absorption for gypsum mapping was assessed, using lab and ASTER reflectance. On the basis of this, two new ASTER gypsum spectral indices (GI1: B4*B8/B6*B7; GI2: B4*B8/B7*B7) were proposed and applied to exclude the interference of hydroxyl-bearing minerals effectively. Seasonal reflectance variability of gypsum formation pixels was confirmed, and it causes the accuracy difference of gypsum indices for multi-seasonal ASTER data. The GI1 achieves the most robust accuracy for multi-seasonal ASTER data with average areas under receiver operator characteristic (ROC) curve of 98.5% and 98.7% for summer and winter ASTER data. Therefore, the GI1 can be used for gypsum mineral mapping, especially in the areas where clay minerals and other hydroxyl-containing minerals are widely distributed.


Introduction
Gypsum (CaSO 4 ⋅2H 2 0), Bassanite (CaSO 4 ⋅0.5H 2 0), and anhydrite (CaSO 4 ) are common Ca-sulfates minerals on the Earth (Bishop et al. 2014). Ca-sulfates (gypsum, bassanite, and anhydrite) typically formed under different temperature conditions in arid evaporation environments (Billo 1987), and they are very important not only for paleo-environmental and paleo-climatic studies but also for industrial applications (components of gypsum boards and cement (Singh and Garg 1995)). Ca-sulfates rocks also play an important role in oil and natural gas exploration in Iran, as they may form the trap where hydrocarbons are accumulated (Gürbüz 2019). In addition, gypsum and bassanite have even been identified in several regions on Mars by both the Observatoire pour la Minéralogie, l´Eau, les Glaces et l´Activité (OMEGA) onboard Mars Express and the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) onboard Mars Reconnaissance Orbiter (MRO) (Wray et al. 2010(Wray et al. , 2011Sefton-Nash et al. 2012;Weitz et al. 2015). Thus, remote detection of Ca-sulfates is essential to mineral exploration, paleo-climatic studies, and minerals mapping on Mars' surface.
Different data-driven methods, such as data fusion (Kavak 2005), Decorrelation Stretching (Chapman et al. 1989;Kavak 2005), Linear Mixture Modeling (Bryant 1996), Support Vector Machines (SVM), Spectral Angle Mapper (SAM) (Soltaninejad et al. 2018), Mixture Tuned Match Filtering (MTMF) (Pour et al. 2018), Spectral Feature Fitting (SFF), and Spectral Mixture Analysis (SMA) (Robertson et al. 2016), were applied to multi-spectral and hyper-spectral images for gypsum identification. However, the spectral properties of gypsum may not be adequately considered in these studies. No targeted enhancement based on gypsum spectral properties was tried in these data-driven methods. And the accuracy of these methods relied on the correct selection of training samples, reference spectrum, or spectral endmember, which was always complicated and challenging. In addition, knowledge-based approaches were applied for gypsum mapping, as well. Feature-oriented principle component selection (FPCS) was applied to ASTER data for gypsum mapping in Turkey (Öztan and Süzen 2011;Özyavaş 2016;Gürbüz 2019) and Pakistan (Khan et al. 2020). Same bands (4,6,8,9) were chosen based on the spectral features of gypsum in these studies, but different principle components were utilized for gypsum detection. Least-squares spectral band-fitting algorithm was applied to Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data for playa evaporite mapping in Death Valley, California, and multiple spectral features of gypsum were chosen for band-fitting analysis (Crowley 1993). Gypsum spectral properties were considered in the FPCS and least-squares spectral band-fitting algorithm; however, the problem of these methods is time-consuming and unstable.
All above image processing methods for gypsum mapping have ignored the characteristic 2.26 μm (ASTER band 7) gypsum absorption, not overlapping with most hydroxylbearing minerals and other major rock-forming minerals. Although the ankerite also shows a 2.26 μm absorption feature (Mars and Rowan 2011), it is easy to distinguish ankerite and gypsum according to the significant absorption feature difference at 2.21 μm or 2.34 μm between the two minerals.
Previous studies held different opinions about the seasonal variation of image reflectance of rocks. The reflectance of Granite outcrops, covered with lichen and algae, was reported to keep consistent all year across the southwest Australian floristic region (Alibegovic et al. 2015). On the contrary, the image reflectance of rocks was considered variable with seasons of the year, because the physical and chemical properties of rocks may change under heavy precipitation and heat depending on their crystal structure (Yao et al. 2017). Namibian seasonal salt pan crust dynamics were monitored by a dense time-series consisting of 77 Landsat 8 cloud-free scenes, and climate controls on crust formation were discussed (Milewski et al. 2020). Moreover, the performance of ASTER mineral spectral indices was proven to be varying according to the data acquisition season (Hewson and Cudahy 2011). In addition, the mutual transformation between Ca-sulfates ( Fig. 1) under laboratory conditions and natural conditions has been presented. Experiments of Harrison (2012) showed that completely dehydrated gypsum (anhydrite) rehydrated to form bassanite at ambient conditions within 1 day. Meanwhile, particulate anhydrite may also be slowly hydrated to form gypsum in a wet environment at low temperatures (Sievert et al. 2005). Under natural conditions, gypsum can be dehydrated to form bassanite and anhydrite under intense solar heating in the arid regions (Warren 2006;Rezaee et al. 2016). Moreover, it was reported that bassanite is present in the dry summer season in Kuwait and hydrates to gypsum in the wet winter months (Gunatilaka et al. 1985). Furthermore, the SWIR absorption feature of Ca-sulfates is closely related to structural water molecules contents (Crowley 1991) (Fig. 2). Consequently, the reflectance of Ca-sulfates may change according to the seasonal change of Ca-sulfates facies. However, to our knowledge, no research has attempted to explore the seasonal variation in image reflectance of Ca-sulfates formations in arid regions.
The objectives of this paper are (1) to assess the SWIR absorption feature for gypsum mapping based on published spectral library and ASTER image reflectance, and to design new gypsum indices (the GI 1 and GI 2 ), eliminating the interference of hydroxyl-containing minerals for gypsum mapping; (2) to analyze seasonal reflectance variability of gypsum formation pixels, and to access the accuracy of spectral indices for multi-seasonal ASTER data.

Regional setting
The study area is located in the Zagros fold-thrust belt, south of Iran (Salati et al. 2019). Climate appears to be a critical factor affecting the distribution of evaporite mineral facies (Gunatilaka et al. 1985). The average rainfall in the study area is about 190 mm/y, falling mostly between November and April. Summer average temperature is about 37 ℃, the ground surface temperature is up to 60 ℃.
The gypsum formations or gypsum-bearing formations in the study area are Gachsaran Formation (Miocene) and the Hormoz Complex (Precambrian). The Gachsaran Formation contains mainly evaporites, and it also contains marls, limestones, and shales (Bahroudi and Koyi 2004). Traditionally, the Kazerun-Qatar fault divides the Gachsaran basin into two different segments: (1) Fars province in the SE and (2) a Northwestern basin (Stocklin 1968). In Fars Province, Gachsaran Formation consists of three members from bottom to top, namely the Chehel Member (anhydrite or gypsum), the Champeh Member (gypseous marl, chalky-gypseous limestone or dolomite, and nodular to massive gypsum), and the Mol Member (gypseous marl and gypseous limestone) (James and Wynd 1965). The Hormoz Complex appears as Salt plugs on the ground (Bosák et al. 1998). The Bam plug locates in the west of the study area. Petrographic composition of the plug is variable with the occurrence of clastic sedimentary, magmatogenic rock types, and chemogenic, including shales, siltstones, limestones, dolostones, gabbro, diorite, quartz diorite, and propyllitized andesite (Bosák et al. 1998). Gypsum is common as layers and interbeds in sediments, and white gypsum in a tremendous amount covers terrace sediments at the eastern plug margin (Bosák et al. 1998). Figure 3b shows a simplified geological map of the study area, which was interpreted from a high-resolution Google image with reference to the LAMAZAN and ROYDAR-E-SOFLA geological maps (1/100000) (Moeini et al. 2004;Aghababaei and Motamedi 2006). Based on previous studies (James and Wynd 1965;Stocklin 1968;Bosák et al. 1998;Bahroudi and Koyi 2004;Mortazavi et al. 2017), the strata of the study area were divided into three categories: gypsum and anhydrite formations (GF), gypsum-bearing formations (BGF), and non-gypsum formations (NGF) ( Table 1). Fig. 1 The mutual transformation between gypsum and bassanite (Gunatilaka et al. 1985) Fig. 2 Lab reflectance of gypsum, bassanite, and anhydrite in SWIR region (1.5-2.9 µm) from USGS mineral spectral library. Gypsum sample is HS333.3B, bassanite is GDS145, and anhydrite is GDS42

Methodology
Six ASTER level 1B (L1B) scenes were acquired in dry summer and wet winter over the study area (Table 2). Digital number (DN) values of ASTER L1B were converted to calibrated spectral radiance at-sensor. Atmospheric correction was carried out, respectively, to VNIR bands and Fig. 3 The ASTER image and simplified geological map of the study area. a The ASTER false color composite image of the study area (bands 3, 2, and 1 assigned in red, green, and blue, respectively); b simplified geological map of the study area modified from Moeini et al. (2004) and Aghababaei and Motamedi (2006), horizontal lines in the polygons represent gypsum and anhydrite formations, black dot pattern represents gypsum-bearing formations 3) was used to approximate and remove the atmospheric contributions from TIR bands radiance data. The datasets and methods used in this study are described in a flowchart (Fig. 4).

Acquiring training samples
Spectral characteristics of "pure" pixels were usually used to identify land-cover types that affect target objects extraction accuracy (Feyisa et al. 2014;Fisher et al. 2016) while attempting to design an index to discriminate target objects and other objects accurately. However, different from water or plants, most of the rocks show spectral variability (Murphy et al. 2012; Zhang and Li 2014) (the difference in whole spectral shape, the variations in spectral absorption position, and depth) due to the difference of minerals proportions, chemistry, structure, grain size of same rock types. The spectral variability makes it difficult to obtain "pure" pixels for rocks and minerals. Here, we used approximate "pure" pixels as training samples to evaluate the capacity of ASTER bands to distinguish between gypsum formations and other formations.
The methods used to acquire approximate "pure" pixels (training samples) of geologic formations include Minimum Noise Fraction (MNF) Transform and the Sequential Maximum Angle Convex Cone (SMACC). MNF transform was used to determine the inherent dimensionality of image data, segregate noise in the data, and reduce the computational requirements for subsequent processing (Luo et al. 2016). SMACC is a tool to find spectral endmembers and their abundances. In this study, we applied MNF transform to ASTER NIR and SWIR reflectance (ASTER 1-9 bands, 2004-11-22). Then, the first 6 output channels of the MNF transform were applied to SMACC tool, and 20 spectral endmembers and their abundance images were acquired. A simplified geological map and Google image were used as a reference to identify spectral endmembers and their abundance images for each geologic formation. Thresholds of 0.60-0.75 were used to extract training samples from the abundance images. The thresholds were decided in a process of trial and error with the reference of a simplified geological map and Google image. As shown in Fig. 5, training samples of alluvial terraces and recent deposit are divided into three groups, based on their internal spectral variability and the difference of sediment source. The sediment sources

Spectral analysis
To assess the 2.26 µm absorption for gypsum mapping, and to develop an index that maximizes the differences of image reflectance values between gypsum and other minerals, lab and ASTER SWIR reflectance of gypsum and other minerals were compared.

Lab reflectance comparison between gypsum and other minerals
We compared lab reflectance between gypsum and several minerals with diagnostic absorption features in SWIR region. As shown in Fig. 6a, absorption positions of gypsum in SWIR region (1.4-2.8 μm) are at 1.49 µm, 1.75 µm, 1.9 µm, 2.21 µm, and 2.26 µm. And it can be identified that the 2.21 µm absorption position of gypsum overlaps with hydroxy-bearing minerals (illite, kaolinite, montmorillonite, chlorite, and muscovite) and weathered feldspar minerals, of which 2200 nm absorption features were considered to be a starting alteration to hydroxy-bearing clays (Graham et al. 1970;Hecker et al. 2010). The overlapping of absorption position is the main resource of commission errors for the gypsum indices, designed based on 2.21 μm (ASTER band6) absorption of gypsum. Moreover, the 1.9 µm absorption position of gypsum also overlaps with other hydrated minerals (illite, kaolinite, and montmorillonite), and 1.49 µm absorption position is outside the atmospheric window. To our knowledge, there is no multi-spectral satellites sensor covering 1.75 μm. Thus, it should be noted that the 2.26 µm gypsum absorption does not overlap with other minerals, and probably can be detected by ASTER data (center wavelength of ASTER band7 is exactly at 2.26 µm). However, ASTER band7 has not been used to calculate a gypsum spectral index.

ASTER reflectance comparison between GF and other formations
Spectral endmembers reflectance of GF, BGF, and NGF was compared in Fig. 6b  show absorption features at 2.33 μm for their lithology of carbonates. It can be noted from Fig. 6b that spectral member reflectance of Qh al1 and Qh al3 also has absorption features at 2.21 μm, attributing to clay minerals.

Calculating new gypsum indices
The primary purpose of designing new gypsum indices is to exclude the potential interference of hydroxyl-bearing minerals under the premise of mapping accuracy assurance. The ASTER SWIR reflectance (4-9 bands) distributions of training samples of GF, BGF, and NGF are shown in Fig. 7. It can be identified that the most significant difference of reflectance distributions between GF (Chm) and other geologic formations is observed in band6 (2.21 μm) and band7 (2.26 μm). However, the reflectance distributions of Chm severely overlap with HS 1 and Qhal 3 in both two bands. The lithology of HS 1 is complicated, except for gypsum, weathered gabbro, diorite, and quartz diorite (illite, montmorillonite, and kaolinite are the main secondary mineral (Pédro 1997;Younis et al. 1997;Deepthy and Balakrishnan 2005)), and propyllitized andesite (altered minerals are mainly chlorite and epidote) also shows absorption feature at 2.21 μm. Alluvium most probably contains a large number of clay minerals (Notesco et al. 2015), showing 2.21 μm absorption. Meanwhile, we found that band4 performed well at discriminating the three geologic formations, and Chm's reflectance in band4 is greater than Qh al3 and HS 1 .

Receiver operator characteristic (ROC) curves
The performance of the different gypsum spectral indices was evaluated using receiver operator characteristic (ROC) curves with ground truth information. ROC curves plot the true-positive rate (TPR) against the false-positive rate (FPR) for a range of threshold values (Fawcett 2006;Fisher et al. 2016). When an ROC curve of a spectral index is the closest to the top-left corner, it represents that the spectral index achieves more accurate mapping results than other indices with a limited range of thresholds. Besides, the proportion of area under the curve also indicates the accuracy of spectral indices. A proportion of 100% means a perfect spectral index; a proportion of 50% means a random spectral index (Fisher et al. 2016). When the area under ROC curve of a spectral index is the largest, it represents that the spectral index achieves more accurate mapping results than other indices with conservative or liberal thresholds. Ground truth regions (pixels) of GF, BGF, and NGF were determined from the simplified geological map of the study area (Fig. 1b).
The optimum threshold of each spectral index, which achieved the lowest combined error (100 − TPR + FPR) for GF range (pixels) (Fisher et al. 2016), was also applied to map GF for multi-seasonal ASTER data. And, the results were compared with the reference of the simplified geological map (Fig. 1b).

Analyzing seasonal reflectance variability of GF
To explore whether Ca-sulfates facies of GF change between dry summer and wet winter in the study area, 2.21 μm (ASTER band6) and 2.26 μm (ASTER band7) reflectance of GF spectral member and training samples were compared for winter and summer multi-seasonal ASTER data.
Seasonal reflectance variability of Chm spectral member and training samples are presented in Figs. 8 and 9. As shown in Fig. 8, Chm spectral member reflectance curves shapes, absorption positions, and absorption depth of summer ASTER data (2001-7-25, 2002-7-28, 2003-7-15) are similar. And the same phenomenon is observed for winter ASTER data (2000-11-18, 2004-3-11, and 2004-11-22). Nevertheless, the differences between winter and summer reflectance of Chm spectral member are as follows: (1) the summer ASTER reflectance is higher than winter overall, attributing to larger solar elevation and stronger solar radiation in summer; (2) 2.26 μm (band7) reflectance shows the most significant change, leading to different slopes between 2.21 μm (band6) and 2.26 μm (band7) for summer and winter ASTER data (Fig. 8). These results were further confirmed by comparing reflectance distributions of Chm training samples between summer and winter ASTER data. As shown in Fig. 9, the summer ASTER band6 (2.21 μm) and band7 (2.26 μm) reflectance of Chm training samples are higher than winter overall; furthermore, the reflectance change magnitude of band7 (2.26 μm) is greater than band6 (2.21 μm).
Seasonal ASTER reflectance variability of GF pixels was confirmed, attributing to the seasonal mutual transformation between Ca-sulfates, possibly. Furthermore, seasonal variability in absorption depth of band6 (2.21 μm) and band7 (2.26 μm), and the slope of reflectance between them is bound to affect the mapping accuracy of gypsum indices, calculated by band6 and band7. Thus, we evaluated the stability of gypsum indices for winter, summer, and multiseasonal ASTER data.
Omission errors were calculated as the number of GF pixels incorrectly classified as a percentage of all GF pixels in each GF category. Chm pixels had fewer omission Fig. 8 Multi-seasonal ASTER reflectance of Chehel (Chm) spectra member. The study area is dry summer between May and October, and wet winter between November and April Fig. 9 Reflectance distributions of training samples of Chehel (Chm) for multi-seasonal ASTER data. Each box plot shows the location of the 10th, 25th, 50th, 75th, and 90th percentiles using horizontal lines, and the circles are 5th and 95th percentiles. The boxes of dry summer ASTER data are marked as red, and blue for wet winter ASTER data errors than HS 2 pixels across most of the gypsum indices. The GI 1 had fewer omission errors for Chm pixels than other indices, while the B8/B6 had fewer omission errors for HS 2 pixels. The GI 1 and GI 2 had 4-6% omission error, the (B10*B12)/(B11*B11), (B4 + B8)/B6 and B4 + B8/B6 had 6-9% omission error, and the (B4 + B8)/B6, B8/B6 and B4/ B9 had 10-18% omission error for Chm pixels. While the B8/B6 had an omission error of less than 1%, the B4 + B8/ B6 had an omission error of less than 10%, the (B4 + B8)/ B6, GI 1 , and (B10*B12)/(B11*B11) had an omission error of 10-21%, the GI 2 , B4/(B6 + B9) and B4/B9 had an omission error of more than 50% for HS 2 pixels.
Commission errors were calculated as the number of GF pixels incorrectly classified as a percentage of all BGF or NGF pixels in each BGF or NGF category. The B4 + B8/B6 achieved greater commission errors for each BGF or NGF (22-94%). Most of the indices had fewer commission errors of Bk (1-9%), except the B8/B6 (25.21%) and B4 + B8/B6 (71.68%). Due to gypsum layers and interbeds contained in BGF, commission errors were greater for BGF (HS 1 , Cpm-Mln) across nearly all indices. And, the GI 1 and GI 2 showed less commission than (B4 + B8)/B6 for HS 1 , while the (B4 + B8)/B6 had less commission than the GI 1 and GI 2 for Cpm-Mln. Grm and Aj-Mn commission errors were the fewest for the (B4 + B8)/B6 (5.58% and 2.37%), and the GI 1 also had fewer commission errors for Grm (6.00%) and Aj-Mn (2.55%). Commission errors of Qh al were fewer for the GI 1 and GI 2 .
The GI 1 and (B4 + B8)/B6 showed better accuracy for GF mapping, using winter ASTER data (2004.11.22). Thus, they were applied to map GF in the study area, using the optimum index thresholds. The results were qualitatively compared, superimposing with the geologic formations boundary. As shown in Fig. 11, both the gypsum indices described the distributions of GF in the study area successfully, and the distribution areas of omission pixels for the two indices were similar, attributing to terrain shadows in ASTER data. We focused on comparing the distribution of commission pixels in the results. In the result of (B4 + B8)/B6, more pixels from alluvail terraces (Qh al ) were misclassified into GF, while fewer Qh al commission errors for the GI 1 , indicating the vital role of band7 (2.26 μm) in excluding the interference of clay minerals for GF mapping. The commission of HS 1 was a problem for both GI 1 and (B4 + B8)/B6, and less commission for the GI 1 . Moreover, the HS 1 commission pixel distribution of the GI 1 matched better with the reported distribution characteristics of gypsum layers and interbeds in HS 1 . Besides, more pixels of Cpm-Mln were classified into GF for the GI 1 , but these commission pixels were distributed as east-west stripes, consistent with distribution characteristics of gypsum interlayer in Cpm-Mln (Fig. 11).

Evaluation of gypsum indices for summer ASTER data
A scene of ASTER image, acquired on July 15, 2003, was selected as representative of summer ASTER data to access the accuracy of gypsum indices for GF mapping. Among the indices, the (B4 + B8)/B6 showed the best performance with the steepest ROC curve (Fig. 10b) and the greatest area under ROC curve. The GI 1 was second only to the (B4 + B8)/B6 in the accuracy of GF mapping, the ROC curve of GI 1 was the closest to (B4 + B8)/B6, and the area under ROC curve of GI 1 was second only to (B4 + B8)/B6. The (B4 + B8)/B6, GI 1 , and GI 2 showed better accuracy than other indices with an area under ROC curve of more than 95% (Table 4). While, the B8/B6, B4/B9, B4 + B8/B6 and B4/(B6 + B9) exhibited poor performance in GF mapping Fig. 10 The accuracy of mapping gypsum, anhydrite formations for a range of index thresholds, using Receiver Operator Characteristic (ROC) curves. The circles represent the optimum index thresholds. a ROC curves of GF mapping for winter ASTER data (2004.11.22). b ROC curves of GF mapping for summer ASTER data (2003.7.15) with areas under ROC curve of less than 90%. When optimum index thresholds were applied, the B4/B9 showed the greatest overall accuracy (96.8%) and TPR (96.8%) but the least GF users' accuracy (13.7%), indicating a large number of commission errors. The (B4 + B8)/B6 achieved the second greatest overall accuracy (95.9%), TPR (95.9%), and GF users' accuracy (72.5%), meanwhile, the least combined error (7.7%) and second least FPR (3.5%). Though the GI 1 had the greatest GF users' accuracy (76.3%) and the least FPR (2.7%), but lower overall accuracy (90.6%) and TPR (90.5%). Besides, the overall accuracy and TPR of the GI 2 were more than 90%. The B8/B6, B4/(B6 + B9), (B10*B12)/ (B11*B11) and GI 2 had 30-40% GF users' accuracy.
Mapping results of the GI 1 and (B4 + B8)/B6, applying to summer ASTER data of July 15, 2003, were also compared qualitatively. Both indices extracted the spatial distribution of Chm member in the study area successfully; however, more omission errors were observed in the marginal zone of Chm regions of the north study area for the GI 1 . And the distribution areas of other omission pixels in Chm were similar, for the two indices, attributing to terrain shadows or anhydrite distribution. While mapping results of HS 2 were quite different (Fig. 12 region A). The (B4 + B8)/B6 showed significantly better performance on recognition of HS 2 pixels than GI 1 , with fewer omission errors. Besides, more commission errors were observed in the results of (B4 + B8)/B6) for Quaternary deposits (Qh al ) and HS 1 , attributing to the overlapping absorption feature on 2.21 μm between gypsum and hydroxyl-containing minerals. While more As-Ja pixels were misclassified into GF for the GI 1 .

Evaluation of gypsum indices for multi-seasonal ASTER data
Three scenes of summer ASTER images (2001.7.15, 2002.7.28, 2003.7.15) and three scenes of winter ASTER images (2000.11.18, 2004.3.11, 2004.11.22) were selected to evaluate the accuracy of gypsum indices for multi-seasonal ASTER data. ROC curves of (B4 + B8)/B6, GI 1 and (B10*B12)/(B11*B11) were plotted (Fig. 13), areas under ROC curves and the lowest combined errors were calculated (Table 5). Overall, the (B4 + B8)/B6, B4/(B6 + B9), (B10*B12)/(B11*B11), GI 1 and GI 2 showed better accuracy than other indices for multi-seasonal ASTER data in the study area, with an average area under ROC curves of more than 95%. Furthermore, the (B4 + B8)/B6 and GI 1 exhibited the most robust mapping accuracy and with areas under ROC curves of more than 98% for six ASTER scenes. For summer ASTER data, the (B4 + B8)/B6 showed the best accuracy with the largest average areas under ROC curves of Fig. 12 The GF mapping performance comparison between the (B4 + B8)/B6 and the GI 1 for summer ASTER data (2003.7.15), applying the optimum index thresholds. a Mapping results of the (B4 + B8)/B6. b Mapping results of the GI 1 99.2%, the lowest combined error of 8.6%. The performance of GI 1 and GI 2 was secondary to (B4 + B8)/B6, with an average area under ROC curves of 98.5% and 95.6%, while, for winter ASTER data, the GI 1 acquired the best accuracy, with the largest average area under ROC curves of 98.7%, the lowest combined error of 10.1%. Besides, the areas under ROC curves of (B4 + B8)/B6, B4/(B6 + B9), (B10*B12)/ (B11*B11) and GI 2 were all more than 95%, exhibiting good mapping performance for winter ASTER data.
In addition, different seasonal variation characteristics of GF mapping accuracy were observed for different gypsum indices. The (B4 + B8)/B6, B8/B6 and B4 + B8/B6 had greater areas under ROC curves for summer ASTER data. While the B4/(B6 + B9), B4/B9, (B10*B12)/(B11*B11), GI 1 , and GI 2 exhibited greater areas under ROC curves for winter ASTER data. It can be observed in Fig. 13a that ROC curves of (B4 + B8)/B6, applied to summer ASTER data are closer to the left top than winter ASTER data. And the (B4 + B8)/B6 had a greater average area under ROC curves and lower combined error for summer ASTER data (Table 5), indicating better performance for summer ASTER data.
While the GI 1 showed better accuracy for winter ASTER data with larger areas under ROC curves and lower combined error. As shown in Fig. 13b, ROC curves of winter ASTER data are steeper. Although ROC curves of 2004.11.22 and2003.7.15 are close, GI 1 achieved better performance on the identification of GF (HS 2 ) pixels in the north of the study area for 2004.11.22 ASTER data ( Fig. 11 and Fig. 12). In addition, HS 2 pixels had more omission errors for 2003.7.15 ASTER data, and the overall accuracy and TPR of the GI 1 was greater for 2004.11.22 ASTER data (Tables 3, 4). For the (B10*B12)/(B11*B11), there was a huge accuracy difference between summer and winter ASTER data. As shown in Fig. 13c, ROC curves for winter ASTER data are much closer to the left top than summer ASTER data. The average area under ROC curve for winter ASTER data was 97.5%, while 93.5% for summer ASTER data. The average lowest combined error for winter ASTER data was 15.4%, while 28.9% for summer ASTER data. Moreover, TPR for winter ASTER data were all more than 90%, while less than 90% for summer ASTER data. Besides, the average area under ROC curve of B4/B9 for winter ASTER data (94.7%) was much greater than summer ASTER data (80.4%), indicating the considerable accuracy difference between summer and winter ASTER data.

Origin of gypsum formation reflectance seasonal variability
ASTER band7 (2.26 µm) appeared obvious seasonal variability for Chm spectral member and GF training samples in this study. It has been proven that the 2.26 µm absorption features of Ca-sulfates were attributed to the proportions of structure molecular water in Ca-sulfates (Hunt 1971;Crowley 1990Crowley , 1991Drake 1995;Cloutis et al. 2006;Harrison 2012;Bishop et al. 2014). Besides, experiments of Bishop Table 3 Gypsum and anhydrite formations classification accuracy statistics for the different index methods, using winter ASTER data (2004.11.22) TPR true-positive rate for the lowest combined error (identical to producers' accuracy). FPR false-positive rate for the lowest combined error (2014) and Harrison (2012) showed that the 2.21 µm and 2.26 µm absorptions decrease in depth, and the absorption positions shift to shorter wavelengths (Harrison 2012), as the gypsum samples dehydrate to bassanite or anhydrite. The climate of the study area is similar to Kuwait, and Casulfate mineral facies were reported to change seasonally in there (Gunatilaka et al. 1985). It can be inferred from these results that seasonal reflectance variability of GF pixels at 2.26 µm may be a response to alternating transformation between bassanite and gypsum with the season. High temperature, extreme drought, and very little vegetation cover provided possible conditions for transformation from gypsum to bassanite in summer (Rezaee and Salari1 2016;Warren 2006). Bassanite was generally considered to be a metastable Ca-sulfates during the conversion between gypsum and anhydrite (Sievert et al. 2005), and only existed steadily in extremely arid climate conditions (Bishop et al. 2014). The summer climate of study area is such arid conditions with few rainfalls, providing probably suitable environment for bassanite. Meanwhile, bassanite was also thought to be a more stable Ca-sulfate phase rather than continuing to rehydrate to form gypsum, without being directly exposed to liquid water (Harrison 2012). Thus, the concentrated precipitation in winter may be the key factor in converting bassanite to gypsum in the study area. Due to the limitation of the research conditions, we do not verify the mutual transformation of gypsum and bassanite between summer and winter in the field survey. And, it should be noticed that GF reflectance seasonal variability is also under the influence of the different solar elevations and radiation intensity between summer and winter ASTER data.
It has been proven that Landsat digital numbers for rocks are proportional to the cosine of the solar zenith angle (Kowalik et al. 1982), and reflectance values are proportional to the slope of the digital numbers versus cosine of the solar zenith angle (Kowalik et al. 1982).

The performance of the proposed new gypsum indices
In this study, two ASTER gypsum indices (the GI 1 and GI 2 ) were developed. Accuracy assessment showed that distributions of GF are identified successfully by the GI 1 for both summer ASTER data and winter ASTER data. And for winter ASTER data, the GI 1 improved the GF mapping accuracy of the study area, achieving the steepest ROC curve (Fig. 10a), the largest average area under ROC curve (98.7%, Table 5), and the lowest average combined error (10.1%, Table 5). While, for summer ASTER data, GF mapping accuracy of the GI 1 is only second to (B4 + B8)/B6 with an average area under ROC curve of 98.5% and average lowest combined error of 14.3%. Besides, the (B4 + B8)/B6 also showed better performance for gypsum mapping in the study area, consistent with the results of Özyavaş (2016).
Similar to the results of previous studies, ASTER band6based gypsum indices are susceptible to the interference of hydroxyl-bearing minerals. Not only pure gypsum rocks but also other authigenic clays and disseminated gypsum crystals were mapped by (B4 + B8)/B6 in the Kohat Plateau, north Pakistan (Khan et al. 2020). "False kaolinite" identification was reported for SWIR hyper-spectral data in Israel, attributing to the proximity of the absorption at 2.20 µm of  (Notesco et al. 2016). In this study, commission errors were more severe for Qh al and HS 1 , geologic formations with more hydroxyl-bearing minerals, across the (B4 + B8)/B6, B8/B6, and B4 + B8/B6 (Tables 3  and 4). While, the GI 1 and GI 2 showed less commission than the three gypsum indices for Qh al and HS 1 (Tables 3 and 4).
There are usually large amounts of clay minerals in alluvail terraces and recent deposite (Qh al ). As shown in Fig. 11, commission pixels of Qh al , distributed on river channels and alluvial fans, for the GI 1 are much fewer than the (B4 + B8)/ B6. Besides, fewer HS 1 pixels were misclassified into GF for the GI 1 . Previous studies (Bosák et al. 1998) showed that gypsum-containing layers or interlayers are distributed in HS 1 in strips, and there is also a large amount of chlorite, epidote, illite, and montmorillonite minerals, formed as secondary weathering minerals on the surface of HS 1 regions. These hydroxyl-bearing minerals have an overlapping 2.21 μm absorption position with gypsum. The mapping results of (B4 + B8)/B6 and GI 1 in the HS 1 region are compared in Fig. 11, it was found that commission pixels of HS 1 distributed in several strips for the GI 1 , consistent with the distribution of gypsum layers and interlayers in HS 1 . While, in the results of (B4 + B8)/B6, a large number of commission pixels scattered in HS 1 disorderly, attributing to distributions of both the gypsum layers, interlayers, and secondary weathering minerals on the surface of rocks. These indicated that the proposed GI 1 successfully identifies gypsum layers or interlayers, and effectively eliminates the interference of the hydroxyl-containing minerals in the background.
We also found that the GI 2 , established using the 2.26 µm absorption feature of gypsum alone, did not achieve the highest accuracy (the mapping accuracy of GI 2 is inferior to (B4 + B8)/B6 and GI 1 for multi-seasonal ASTER data with an average area under ROC curve of 96.2%). While the GI 1 , combined using the 2.21 µm and 2.26 µm absorption features of gypsum, showed the best mapping accuracy for winter ASTER data, and the second highest mapping accuracy for summer ASTER data. Furthermore, the GI 1 is the most robust gypsum index for multi-seasonal ASTER data with average areas under ROC curve of 98.5% and 98.7% for summer and winter ASTER data. In summary, the 2.21 µm strong absorption feature is critical for gypsum mapping, while the 2.26 µm absorption helps improve the interference of the hydroxyl-containing minerals for gypsum mapping. However, we also found that commission error of carbonates is a problem for the GI 1 and GI 2 , a few pixels of carbonate rocks (As-Ja) were misclassified into GF. Thus, background lithology also should be considered for the selection of gypsum indices.

The performance of gypsum indices for multi-seasonal ASTER data
Previous studies reported that MgOH/carbonate RBD products (B6 + B9)/B8 show better accuracy for summer ASTER data than winter ASTER data, attributing to the difference of terrain shadow and sun altitude (Hewson and Cudahy 2011). Nevertheless, the accuracy of different gypsum indices showed markedly different patterns in terms of the effect of ASTER acquisition season in this study. The B4/(B6 + B9), B4/B9, (B10*B12)/(B11*B11), GI 1 , and GI 2 showed better accuracy for winter ASTER data, and the GI 1 achieved the best accuracy. While the (B4 + B8)/ B6, B8/B6, B4 + B8/B6 showed better accuracy for summer ASTER data, and the (B4 + B8)/B6 achieved the best accuracy. We also found that gypsum indices, calculated mainly by dividing by the reflectance of band7 (2.26 µm) or band9 (2.40 µm), exhibited better accuracy for winter ASTER data. While gypsum indices, designed based on the band6 (2.21 µm) gypsum absorption, show better accuracy for summer ASTER data. Experimental results showed that seasonal variability of accuracy for gypsum indices is caused by the seasonal change of reflectance curve shape in SWIR Fig. 13 The accuracy of mapping GF for a range of index thresholds, using Receiver Operator Characteristic (ROC) curves. The circles represent the optimum index thresholds. a ROC curves of the (B4 + B8)/B6 for multi-seasonal ASTER data. b ROC curves of the GI 1 for multi-seasonal ASTER data. c ROC curves of the (B10*B12)/ (B11*B11) for multi-seasonal ASTER data region. The weakening of ASTER band7 (2.26 µm) reflectance absorption feature of Chm pixels for summer ASTER data is the reason for the better accuracy of the GI 1 and GI 2 , applied to winter ASTER data.

Conclusions
Three major conclusions are drawn from this study.
(1) Seasonal reflectance variability of gypsum pixels was confirmed, and the seasonal variability of 2.26 µm reflectance may be a response to mutual transformation between bassanite and gypsum with the season.
(2) The effectiveness of 2.26 µm (ASTER band7) reflectance in gypsum mapping was validated. On the basis of 2.21 µm and 2.26 µm absorption features of gypsum, two new ASTER gypsum indices (GI 1 and GI 2 ) were proposed, improving the interference of hydroxyl-containing minerals for gypsum mapping. Furthermore, the GI 1 achieved the most robust mapping accuracy for multiseasonal ASTER data.
The next step will be to verify the reason for the seasonal reflectance variability of gypsum pixels in the arid region and explore the possibility of subdividing Ca-sulfates types (anhydrite, bassanite, and gypsum) by ASTER data. The proposed GI 1 can be applied in similar arid regions for gypsum mapping.