The Application of ALOS/PALSAR InSAR to Measure Subsurface Penetration Depths in Deserts

Spaceborne Synthetic Aperture Radar (SAR) interferometry has been utilised to acquire high-resolution Digital Elevation Models (DEMs) with wide coverage, particularly for persistently cloud-covered regions where stereophotogrammetry is hard to apply. Since the discovery of sand buried drainage systems by the Shuttle Imaging Radar-A (SIR-A) L-band mission in 1982, radar images have been exploited to map subsurface features beneath a sandy cover of extremely low loss and low bulk humidity in some hyper-arid regions such as from the Japanese Earth Resources Satellite 1 (JERS-1) and Advanced Land Observing Satellite/Phased Array type L-band Synthetic Aperture Radar (ALOS/PALSAR). Therefore, we hypothesise that a Digital Elevation Model (DEM) derived by InSAR in hyper-arid regions is likely to represent a subsurface elevation model, especially for lower frequency radar systems, such as the L-band system (1.25 GHz). In this paper, we compare the surface appearance of radar images (L-band and C-band) with that of optical images to demonstrate their different abilities to show subsurface features. Moreover, we present an application of L-band InSAR to measure penetration depths in the eastern Sahara Desert. We demonstrate how the retrieved L-band InSAR DEM appears to be of a consistently 1–2 m lower elevation than the C-band Shuttle Radar Topography Mission (SRTM) DEM over sandy covered areas, which indicates the occurrence of penetration and confirms previous studies.


Introduction
Hyper-arid regions remain one of the most desolate and inhospitable places on Earth. However, some of them host many subsurface natural resources, such as oil, gas, groundwater, etc. The exploitation of these subsurface resources needs substantial reconnaissance work, such as field surveys and in-situ measurements, which can be very time and labour consuming and hence costly. Nowadays, with the development of various spaceborne and airborne remote sensors at microwave frequencies, it may be feasible to study the subsurface part of hyper-arid regions. The successful acquisition of L-band SAR images by SEASAT in 1978 [1] led to several follow-up missions: SIR-A (1981), SIR-B (1984), and two SIR-C/X-SAR missions (1994). These missions revealed that at L-band wavelengths, radar has a limited penetration ability over the Eastern Sahara [2], although coverage with penetration was very limited.
Regarding the ability to identify subsurface features, SAR images are superior to optical images because their longer wavelength leads to subsurface penetration [1,2]. The penetration of waves occurs under the condition of a small scattering loss within the surface cover layer, which has grain sizes smaller than 10% of the incident wavelength and has a moisture content of less than 1% [3].

Study Site
Kufra is in southeastern Libya and is within the Eastern Sahara region. By using Radarsat-1 C-band (5.6 GHz/5.5 cm) images, Robinson et al. interpreted subsurface paleo-drainage, which passes through this region from south to east until Sarir Dalmah in Libya [5,28]. This paleo-drainage does not coincide with that recorded in the United States Geological Survey (USGS) HYDRO1k database in this area, which suggests that it may be of subsurface origin [5]. Following this, a more complete L-band radar coverage of the eastern Sahara by the Japanese JERS-1 satellite was used to create the first regional-scale radar mosaic covering Egypt, northern Sudan, eastern Libya, and northern Chad [29]. Subsequently, ALOS/PALSAR L-band data allowed, for the first time, an accurate mapping of a continuous 900 km-long paleo-drainage system, which is termed the Kufra River [10].
The sand cover layer in this area overlies the Nubian Sandstone Series, which are composed of sands, sandstones, clays, and shales [30,31]. The Nubian Sandstone Aquifer System (NSAS) is the largest known fossil water aquifer system on Earth. It is located underground in the Eastern Sahara Desert and spans just over 2,000,000 km 2 , including north-western Sudan, north-eastern Chad, south-eastern Libya, and most of Egypt [32]. It contains an estimated 150,000 km 3 of groundwater [33]. The aquifer is mainly composed of hard ferruginous sandstone with shale and clay intercalation, having a thickness that ranges between 140-230 m [33].
Our study site is selected to be situated over one of the eastern tributaries of the Kufra River, which arises in northern Uweinat close to the Sudanese border and is named the Uweinat tributary [5]. The coverage (centre latitude and longitude are 23.265 • N, 23.603 • E) of the InSAR image pair is shown as a red rectangle in Figure 1a. Figure 1b shows the ALOS/PALSAR HH amplitude image, along with three close-up pictures from Google Earth to show three bright spots in the ALOS/PALSAR HH amplitude image. Figure 1(c1,c2) display two areas which appear to be plants covered by sand and Figure 1(c3) shows an area which seems to be outcrop rocks. The background image is from Google Earth and is georeferenced in QGIS. The blue lines denote the paleo-drainage system which is visible within the ALOS/PALSAR image [12]. The red rectangle is the coverage of one ALOS/PALSAR frame, which shows the study site; (b) the ALOS/PALSAR HH amplitude image (70 × 110 km); (c1-c3) show three areas: (c1) appears to be the plants covered by sand (white diamond denotes the location of a reference point used in InSAR processing), (c2) showing a bright spot in the channel area (b), which also appears to be plants covered by sand and (c3) showing a bright area in the ALOS/PALSAR image, which may be outcrop rocks.

Comparison of Image Surface Morphology
The occurrence of radar penetration in the Selima Sand Sheet in the Eastern Sahara, where there is a fine-grained sand cover, was first revealed by comparing radar images with visible images [2]. The buried subsurface features appear in radar images, while remain unobserved in corresponding visible images. The radar waves can penetrate the sandy cover layer and are less likely to be reflected or scattered at the surface, mainly attributed to its longer wavelength. The occurrence and brightness of subsurface features in radar images depends on many factors, including the wavelength, incidence angle, aspect angle, polarisation of the incident wave, surface roughness, grain size, moisture, etc. Dark tones in the radar images indicate low backscatters, which can be explained by two scenarios. One scenario suggests a smooth surface of the radar wavelength , which causes specular reflection of the incident radar waves. Alternatively, the incident radar wave penetrates into the sand layer and is absorbed or scattered in this sandy cover material. In this condition, the volume scattering V is usually high enough to cause a large volume scattering loss within the covering material. On the other hand, the bright tone in the radar image, indicates higher backscatter, which may be caused by outcropping rocks with a large radar cross section. The medium tone can be explained by a rough surface if no penetration occurs or by the scenario in which the radar wave transmits into the subsurface and is reflected by a shallow layer, while part of the wave energy is absorbed or scattered during the transition between the surface and subsurface layer.
In this study, we compare the surface morphology of an L-band ALOS/PALSAR image against an optical image from Landsat, an SAR image, and DEM from C-band SRTM. The SRTM [21] was an 11-day mission, which aimed to produce a near global DEM with a resolution of up to 1 arc-second The background image is from Google Earth and is georeferenced in QGIS. The blue lines denote the paleo-drainage system which is visible within the ALOS/PALSAR image [12]. The red rectangle is the coverage of one ALOS/PALSAR frame, which shows the study site; (b) the ALOS/PALSAR HH amplitude image (70 × 110 km); (c1-c3) show three areas: (c1) appears to be the plants covered by sand (white diamond denotes the location of a reference point used in InSAR processing), (c2) showing a bright spot in the channel area (b), which also appears to be plants covered by sand and (c3) showing a bright area in the ALOS/PALSAR image, which may be outcrop rocks.

Comparison of Image Surface Morphology
The occurrence of radar penetration in the Selima Sand Sheet in the Eastern Sahara, where there is a fine-grained sand cover, was first revealed by comparing radar images with visible images [2]. The buried subsurface features appear in radar images, while remain unobserved in corresponding visible images. The radar waves can penetrate the sandy cover layer and are less likely to be reflected or scattered at the surface, mainly attributed to its longer wavelength. The occurrence and brightness of subsurface features in radar images depends on many factors, including the wavelength, incidence angle, aspect angle, polarisation of the incident wave, surface roughness, grain size, moisture, etc. Dark tones in the radar images indicate low backscatters, which can be explained by two scenarios. One scenario suggests a smooth surface of the radar wavelength λ, which causes specular reflection of the incident radar waves. Alternatively, the incident radar wave penetrates into the sand layer and is absorbed or scattered in this sandy cover material. In this condition, the volume scattering V is usually high enough to cause a large volume scattering loss within the covering material. On the other hand, the bright tone in the radar image, indicates higher backscatter, which may be caused by outcropping rocks with a large radar cross section. The medium tone can be explained by a rough surface if no penetration occurs or by the scenario in which the radar wave transmits into the subsurface and is reflected by a shallow layer, while part of the wave energy is absorbed or scattered during the transition between the surface and subsurface layer.
In this study, we compare the surface morphology of an L-band ALOS/PALSAR image against an optical image from Landsat, an SAR image, and DEM from C-band SRTM. The SRTM [21] was an 11-day mission, which aimed to produce a near global DEM with a resolution of up to 1 arc-second (≈30 m). Apart from the global DEM, the mission also provides a C-band combined image (sub-swath images of HH and VV are combined [34]).

Interferometic Phase Formation
SAR interferometry is used to calculate the interference pattern caused by the difference in phases between the two images acquired by a space-borne SAR, which is usually processed from two repeat passes of the satellite at two different times. Figure 2 shows the geometry of SAR interferometry. The resulting difference of phases (wrapped in the range of [-π, π]) forms a new kind of image called an interferogram. Equation (1) shows the formation of the interferogram. M and S are the master and slave images, respectively, and I m , I s , ϕ m , ϕ s denote the intensities and phases of the master and slave images. The interferometric phase, ∆ϕ, a pattern of fringes containing all of the information on the relative geometry, including contributions of the topography ϕ top , the orbital trajectories ϕ orb , atmospheric effects including ϕ tro from the troposphere and ϕ ion from the ionosphere, surface displacement ϕ dis , and phase resulting from speckle noise ϕ noi . These phase components can be expressed as Equation (2) [35].
Remote Sens. 2017, 9, 638 4 of 19 (≈30 m). Apart from the global DEM, the mission also provides a C-band combined image (sub-swath images of HH and VV are combined [34]).

Interferometic Phase Formation
SAR interferometry is used to calculate the interference pattern caused by the difference in phases between the two images acquired by a space-borne SAR, which is usually processed from two repeat passes of the satellite at two different times. Figure 2 shows the geometry of SAR interferometry. The resulting difference of phases (wrapped in the range of -π, π ) forms a new kind of image called an interferogram. Equation (1) shows the formation of the interferogram. and are the master and slave images, respectively, and , , , denote the intensities and phases of the master and slave images. The interferometric phase, , a pattern of fringes containing all of the information on the relative geometry, including contributions of the topography , the orbital trajectories , atmospheric effects including from the troposphere and from the ionosphere, surface displacement , and phase resulting from speckle noise . These phase components can be expressed as Equation (2)  In conventional differential SAR interferometry, the remaining phase after removing the orbital and topographic phase, which is termed the differential phase, should include , , , and . However, if an L-band SAR image pair is utilised to form an interferogram and an external Cband SRTM DEM is employed to simulate the topographic phase, an additional phase component, _ , which results from the difference between L-and C-band, should be included. Therefore, the differential phase in Equation (2) changes into Equation (3). Equations (4) and (5) express the phases caused by topographic variation and surface displacement , where , , are the perpendicular baseline, slant range, and incident angle, respectively. A flowchart of the InSAR processing and analysis of phase components is shown in Figure 3.  In conventional differential SAR interferometry, the remaining phase after removing the orbital and topographic phase, which is termed the differential phase, should include ϕ tro , ϕ ion , ϕ dis , and ϕ noi . However, if an L-band SAR image pair is utilised to form an interferogram and an external C-band SRTM DEM is employed to simulate the topographic phase, an additional phase component, ϕ ele_di f , which results from the difference between L-and C-band, should be included. Therefore, the differential phase in Equation (2) changes into Equation (3). Equations (4) and (5) express the phases caused by topographic variation p and surface displacement d, where B ⊥ , R, θ are the perpendicular baseline, slant range, and incident angle, respectively. A flowchart of the InSAR processing and analysis of phase components is shown in Figure 3.
Remote Sens. 2017, 9, 638 5 of 19 = 4 (5) Figure 3. Flowchart of InSAR processing and analyses of phase components used in this study. The interferograms contains phase contributions from topography, orbital error, tropospheric and ionospheric delay, and the remaining phase. The remaining phase can be caused by surface displacement, topographic phase residual, which may result from the different penetration depths detected by L-and C-band SAR.
SAR data from ALOS/PALSAR is processed using InSAR in this study. ALOS/PALSAR [24] is a radar system launched by the Japan Aerospace Exploration Agency (JAXA) that can also be used to derive high resolution DEMs at 30 m. Four acquisitions of ALOS/PALSAR data are collected and their key features are listed in Table 1  InSAR utilises the phase differences of the waves, so it requires signals from two different acquisitions to be coherent. Multi-look processing is a traditional method to reduce phase noise in SAR processing, while it reduces the resolution of interferograms. The phase noise can be estimated from the interferometric SAR pair using coherence. The coherence is the cross-correlation coefficient of the SAR image pair, which is estimated over a small window (a few pixels in range and azimuth). The theoretical elevation dispersion (standard deviation) can be calculated from the perpendicular baseline and coherence as Equation (6) [37]. The interferograms contains phase contributions from topography, orbital error, tropospheric and ionospheric delay, and the remaining phase. The remaining phase can be caused by surface displacement, topographic phase residual, which may result from the different penetration depths detected by L-and C-band SAR. SAR data from ALOS/PALSAR is processed using InSAR in this study. ALOS/PALSAR [24] is a radar system launched by the Japan Aerospace Exploration Agency (JAXA) that can also be used to derive high resolution DEMs at 30 m. Four acquisitions of ALOS/PALSAR data are collected and their key features are listed in Table 1  InSAR utilises the phase differences of the waves, so it requires signals from two different acquisitions to be coherent. Multi-look processing is a traditional method to reduce phase noise in SAR processing, while it reduces the resolution of interferograms. The phase noise can be estimated from the interferometric SAR pair using coherence. The coherence is the cross-correlation coefficient of the SAR image pair, which is estimated over a small window (a few pixels in range and azimuth). The theoretical elevation dispersion (standard deviation) can be calculated from the perpendicular baseline and coherence as Equation (6) [37].
where λ is the wavelength, R is the range, θ is the incidence angle, B ⊥ is the perpendicular baseline, and NL, r are the number of looks and coherence, respectively. Figure 4a shows the relationship between coherence and elevation dispersion for different perpendicular baselines when NL = 8, while Figure 4b presents the relationship between coherence and elevation dispersion for different numbers of looks when the perpendicular baseline is set to 2000 m. We set NL = 8 to improve the coherence and maintain the spatial resolution (range direction is~37 m and azimuth direction is~56 m) in this study.
We can see from Figure where is the wavelength, is the range, is the incidence angle, is the perpendicular baseline, and , are the number of looks and coherence, respectively. Figure 4a shows the relationship between coherence and elevation dispersion for different perpendicular baselines when NL = 8, while Figure 4b presents the relationship between coherence and elevation dispersion for different numbers of looks when the perpendicular baseline is set to 2000 m. We set NL = 8 to improve the coherence and maintain the spatial resolution (range direction is ~37 m and azimuth direction is ~56 m) in this study. We can see from Figure

Removal of Topographic and Orbital Phase
Topography is one of the main contributions to the interferometric phase . Topographic phase, , yields fringes that hug the topography that looks like contour lines. Topographic phase can be simulated by using an available external DEM. The external DEM used to simulate topographic phase is SRTM C-band DEM (Version 3, resolution of 1-arc second ~30 m). Table 2 lists the characteristics of these two systems. The SRTM DEM heights are referenced to the Earth Gravitational 1996 (EGM96) geoid, while the World Geodetic System (WGS84) ellipsoid is used in InSAR processing. We converted the datum of SRTM C-band DEM from EGM96 to WGS84 by using an EGM96 15' geoid height model [38].

Removal of Topographic and Orbital Phase
Topography is one of the main contributions to the interferometric phase ∆ϕ. Topographic phase, ϕ top , yields fringes that hug the topography that looks like contour lines. Topographic phase can be simulated by using an available external DEM. The external DEM used to simulate topographic phase is SRTM C-band DEM (Version 3, resolution of 1-arc second~30 m). Table 2 lists the characteristics of these two systems. The SRTM DEM heights are referenced to the Earth Gravitational 1996 (EGM96) geoid, while the World Geodetic System (WGS84) ellipsoid is used in InSAR processing. We converted the datum of SRTM C-band DEM from EGM96 to WGS84 by using an EGM96 15' geoid height model [38]. The orbital phase ϕ orb , which is another contribution to the ∆ϕ, is caused by the orbital trajectories. It can be removed given the knowledge of the satellite trajectories. For some missions, this knowledge is not accurate to the scale of the wavelength, leading to errors in the baseline estimation and leaving regular fringes in the interferogram. These periodic fringes can be used in turn to refine the relative separation between the two trajectories, which is termed as baseline refinement. In ROI_PAC, the baseline refinement is realised by a quadratic model fitted between the total phase and surface elevation from an external topographic data source [37]. Baseline refinement is applied to the interferograms when periodic fringes are observed after removing topographic orbital phase. The remaining phase is unwrapped by using the phase unwrapping algorithm, the Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU) [39]. Radar waves of very low frequency propagate through the ionosphere undergoing Faraday Rotation (FR), i.e., rotation of the polarisation vector. This is mainly caused by the anisotropy of Total Electron Content (TEC) in the ionosphere, the values of which were currently very low in the 11-year solar cycle [40]. In addition, most PALSAR acquisitions occur at a local time of 23:00, when TEC activity is also typically low [41].
The radar wave is delayed when it propagates through the troposphere and ionosphere, resulting in ϕ tro and ϕ ion respectively. The sum of these two components is called the atmospheric effects. Phase delay through the troposphere is primarily caused by Precipitable Water Vapour (PWV). In hyper-arid regions, such as the Sahara Desert, the temporal and spatial variation of the PWV within an SAR frame is likely to be less significant than over other areas because of the low humidity in this region. To confirm this, the MODerate-resolution Imaging Spectro-radiometer (MODIS) Precipitable Water product (MYD05 in collection 06) and the ERA-Interim (a global atmospheric reanalysis) data from the European Centre for Medium-Range Weather Forecasts (ECMWF) are collected to compute the phase delay differences between two acquisitions, which is called the Atmospheric Phase Screen (APS). The methods used to generate an APS from the MODIS Precipitable Water product and ECMWF data are discussed in [42,43], respectively. The dates and times when these datasets were acquired are listed in Table 3. The phase component ϕ dis is caused by surface displacement during the time period between the SAR image pair. Surface displacement might be due to crustal movement, groundwater withdrawal, or sand movements in this study site. However, there is no nearby fault [44] and no earthquake was reported from 13 August 2007 to 31 December 2008 [45] over this study site. In a wider region which covers this study area, the land deformation that is caused by groundwater withdrawal is on the level of a few millimetres per year [46]. This will cause an elevation error on the order of 0.5 m when it is wrongly interpreted as ϕ ele_di f . Sand movements are composed of a mixture of deposition and erosion. However, if the L-band ALOS/PALSAR were to observe sand movements, it would probably result in the de-correlation between two acquisitions.

Estimation of Penetration Depth
If the differential phase analysed in Section 3.2.3 is mainly caused by a difference between Land C-band SAR, then we can invert the L-band and C-band Elevation Difference (LCED, L-band minus C-band) by using Equation (4). The differential phase is converted to height by using the Delft Object-oriented Radar Interferometric Software (DORIS) developed by the Technical University of Delft (TUDelft) [47]. The resulting LCED map would indicate the L-band C-band penetration difference if the LCED has a negative shift in the channel areas compared to that of the non-channel areas. Theoretically, LCED on rock outcrops should be zero and be a negative value in the area where there is a fine-grained sand cover.
To confirm this idea, we classified the study area into three classes, i.e., hilly areas, channel areas, and plain areas, by applying a supervised classification method, the Support Vector Machine (SVM) [48], to an ALOS/PALSAR amplitude image, SRTM C-band DEM, and InSAR coherence. These three datasets are chosen for classification because the channel area is most pronounced in these images. The amplitude decreases when the elevation difference increases according to a study of elevation difference between SRTM DEM and ICESat elevations in a wider region in Kufra [49]. Samples fed to the SVM are chosen manually and separated in training and test parts to evaluate the accuracy of the classification.
Finally, histograms of LCED within each group are analysed. For each of these histograms, the mean value of the histogram will be regarded as the penetration depth, which represents all of the measurements in the surface type. The standard error of mean, σ SM = σ SD / √ N, can be used to estimate the accuracy of this mean value. The σ SD is the standard deviation of the histogram and N is the number of samples, which is the pixel numbers within the surface type area in this case. Due to the spatial correlation of the LCED map, a variogram can be used calculate the spatial correlation distance, R [50,51]. Therefore, the independent sample numbers, N, can be calculated by dividing the area A of each surface type into areas with a radius of correlation range, which is N = A/R 2 .

Validation of the Penetration Depth
The LCED indicates the differing penetration depths between the L-band PALSAR and C-band SRTM. Since SRTM C-band DEM may also detect elevations lower than the surface, it is necessary to validate the SRTM C-band DEM to surface elevations. In the literature, the SRTM C-band DEM shows an average penetration depth of 1.7 m when compared with SRTM X-band DEM and around 7-9 m when compared with ICESat data over glaciers in the Pamir-Karakoram-Himalaya [52]. The SRTM X-band DEM shows systematic bias from SRTM C-band DEM, which is pronounced in Africa and South America [53]. ICESat elevations are considered to represent the surface elevations over glaciers and peatland areas [52,54]. Therefore, we employ the land altimetry product (GLA14) from ICESat as a reference. ICESat is an earth observation mission primarily designed for measuring ice sheet mass balance from 2003 to 2009. ICESat/GLA14 elevation data is acquired by altimetry LiDAR and has a vertical accuracy up to 0.04 m [26,55,56]. The land altimetry of ICESat is carried out over footprints of a 70 m diameter and spaced every 150 m. The ICESat altimetry data is referenced to TOPEX/Poseidon. We transformed the vertical datum of all datasets to be the same WGS84 ellipsoid. Both the SRTM C-band amplitude image ( Figure 5b) and the ALOS/PALSAR image (Figure 5e) not only clearly display the main channel, but also show some branches. There is a very high contrast between the edges of the channels. As seen from the blue arrow in Figure 5b,e, the ALOS/PALSAR image shows a subtle branch, while in the SRTM amplitude image, no such feature is visible. This observation implies that the ALOS/PALSAR detects deeper subsurface features than the SRTM DEM.

Interferometic Phase Formation
The temporal and spatial baselines of the four ALOS/PALSAR images were calculated and are shown in Figure 6. Six radar image pairs can, in theory, be formed by using the four radar acquisitions. However, two of the pairs (  Both the SRTM C-band amplitude image (Figure 5b) and the ALOS/PALSAR image (Figure 5e) not only clearly display the main channel, but also show some branches. There is a very high contrast between the edges of the channels. As seen from the blue arrow in Figure 5b,e, the ALOS/PALSAR image shows a subtle branch, while in the SRTM amplitude image, no such feature is visible. This observation implies that the ALOS/PALSAR detects deeper subsurface features than the SRTM DEM.

Interferometic Phase Formation
The temporal and spatial baselines of the four ALOS/PALSAR images were calculated and are shown in Figure 6. Six radar image pairs can, in theory, be formed by using the four radar acquisitions.

Removal of Topographic and Orbital Phase
Topographic phase was simulated by using a C-band SRTM DEM and then subtracting this from the initial interferograms to form differential interferograms. These are shown in Figure 7(a1-a4). These interferograms include all the phase components in Equation (3). The phase ramps in these differential interferograms are along both the range and azimuth directions. These phase ramps may be caused by imprecise baselines which are estimated from inaccurate orbit information of the satellites or can be caused by both orbital inaccuracy and the different atmospheric delay between two acquisitions. In other words, orbital phases and atmospheric phase may be mingled together, and therefore behave as a phase ramp along both the azimuth and range directions [57]. Figure 7(b1-b4) shows the interferograms after baseline refinement. The phase ramps appear to have been successfully removed. Phase distributions in Figure 7(b2,b3) resemble each other, while Figure 7(b1) has more spatial variation and Figure 7(d4) is more uniform in space.

Analysis of Remaining Phase Components
For the analysis of ionospheric effects, we investigated the Faraday rotation angles from the PALSAR catalogue metadata located at the Alaska Satellite Facility (ASF). For SAR acquisitions of 13 February 2008 and 31 December 2008, the estimated Faraday rotation angles are 1.3° and 1.4°, respectively. Besides, the interferograms in Figure 7(b1-b4) do not show noticeable artifacts along the azimuth direction. Therefore, we think that the ionospheric effects have no significant impact on the interferograms and can thus be neglected in this study.
We analysed the tropospheric contributions to the interferograms by using an APS, which can be either calculated from thus ECMWF temperature and pressure data, or from the MODIS Precipitable Water product (MYD05 in collection 06). The results of these two methods are shown in Figure 7(c1-c4),(d1-d4), respectively. It was cloudy on 13 August 2007 in this region, so Figure 7(d1) shows almost no result from the MODIS PWV. We can see from Figure 7(c1-c4) that the APS modelled from the ECMWF temperature and pressure data are uniform in the range direction, but follow a linear trend along the azimuth direction. This may be the azimuth component to the phase ramp in Figure 7(a1-a4), and if it is, it may have been removed with orbital phase during the baseline refinement. The spatial variations of these four APS are much less than 1 radian, which corresponds to an elevation variation of 4.46 m when the perpendicular baseline is 2000 m. However, compared to the APS modelled by ECMWF data, the APS shown in Figure 7(d1-d4) have more spatial variation, which are within 1 radian. According to [58], errors in the derived water vapour of MODIS nearinfrared retrieval are estimated to be 5-10%, and errors can be up to 14% under hazy conditions. Considering that new systematic errors may be imported if these APS are not accurate, in this study, we do not subtract the APS from interferograms. We suspect that the interferometric phase in Figure  7(b1) suffers from the worst atmospheric effects among these four pairs. The spatial similarity

Removal of Topographic and Orbital Phase
Topographic phase was simulated by using a C-band SRTM DEM and then subtracting this from the initial interferograms to form differential interferograms. These are shown in Figure 7(a1-a4). These interferograms include all the phase components in Equation (3). The phase ramps in these differential interferograms are along both the range and azimuth directions. These phase ramps may be caused by imprecise baselines which are estimated from inaccurate orbit information of the satellites or can be caused by both orbital inaccuracy and the different atmospheric delay between two acquisitions. In other words, orbital phases and atmospheric phase may be mingled together, and therefore behave as a phase ramp along both the azimuth and range directions [57]. Figure 7(b1-b4) shows the interferograms after baseline refinement. The phase ramps appear to have been successfully removed. Phase distributions in Figure 7(b2,b3) resemble each other, while Figure 7(b1) has more spatial variation and Figure 7(d4) is more uniform in space.

Analysis of Remaining Phase Components
For the analysis of ionospheric effects, we investigated the Faraday rotation angles from the PALSAR catalogue metadata located at the Alaska Satellite Facility (ASF). For SAR acquisitions of 13 February 2008 and 31 December 2008, the estimated Faraday rotation angles are 1.3 • and 1.4 • , respectively. Besides, the interferograms in Figure 7(b1-b4) do not show noticeable artifacts along the azimuth direction. Therefore, we think that the ionospheric effects have no significant impact on the interferograms and can thus be neglected in this study.
We analysed the tropospheric contributions to the interferograms by using an APS, which can be either calculated from thus ECMWF temperature and pressure data, or from the MODIS Precipitable Water product (MYD05 in collection 06). The results of these two methods are shown in Figure 7(c1-c4),(d1-d4), respectively. It was cloudy on 13 August 2007 in this region, so Figure 7(d1) shows almost no result from the MODIS PWV. We can see from Figure 7(c1-c4) that the APS modelled from the ECMWF temperature and pressure data are uniform in the range direction, but follow a linear trend along the azimuth direction. This may be the azimuth component to the phase ramp in Figure 7(a1-a4), and if it is, it may have been removed with orbital phase during the baseline refinement. The spatial variations of these four APS are much less than 1 radian, which corresponds to an elevation variation of 4.46 m when the perpendicular baseline is 2000 m. However, compared to the APS modelled by ECMWF data, the APS shown in Figure 7(d1-d4) have more spatial variation, which are within 1 radian. According to [58], errors in the derived water vapour of MODIS near-infrared retrieval are estimated to be 5-10%, and errors can be up to 14% under hazy conditions. Considering that new systematic errors may be imported if these APS are not accurate, in this study, we do not subtract the APS from interferograms. We suspect that the interferometric phase in Figure 7(b1) suffers from the worst atmospheric effects among these four pairs. The spatial similarity between Figure 7(b2,b3) provides indirect evidence that these two pairs do not suffer from serious atmospheric effects.

Estimation of Penetration Depth
We converted the remaining phase shown in Figure 7(b1-b4) to the LCED maps which are shown in Figure 8(a1-a4), along with the InSAR coherence maps (Figure 8(b1-b4)). During the conversion, we utilised a reference point (denoted by the white diamond in Figure 8(a1-a4)), which is located at 23.52 • N, 23.604 • E and is of high coherence, to be the zero-phase reference. We suspect that the first InSAR pair of 13 August 2007-13 February 2008 suffers from severe atmospheric effects because it involves the acquisition on 13 February 2008 when it was cloudy, according to the MODIS observation on this day. Regarding 30 September 2008-31 December 2008, the perpendicular baseline is too short to derive the topographic variation of the interferometric phase. For instance, a topographic variation of 5 m only contributes to a phase change of 0.14 radians within the InSAR frame according to Equation (6). Therefore, we discard two of the InSAR pairs from 13 August 2007-13 February 2008 and 30 September 2008-31 December 2008. The elevation difference caused by the L-and C-band SAR system is a spatial and temporal low-frequency signal in the LCED maps. Two LCED maps in Figure 8(a2,a3) are averaged into a final LCED map to reduce the impact from high-frequency components.

Estimation of Penetration Depth
We converted the remaining phase shown in Figure 7(b1-b4) to the LCED maps which are shown in Figure 8(a1-a4), along with the InSAR coherence maps (Figure 8(b1-b4)). During the conversion, we utilised a reference point (denoted by the white diamond in Figure 8(a1-a4)), which is located at 23.52°N, 23.604°E and is of high coherence, to be the zero-phase reference. We suspect that the first InSAR pair of 13 August 2007-13 February 2008 suffers from severe atmospheric effects because it involves the acquisition on 13 February 2008 when it was cloudy, according to the MODIS observation on this day. Regarding 30 September 2008-31 December 2008, the perpendicular baseline is too short to derive the topographic variation of the interferometric phase. For instance, a topographic variation of 5 m only contributes to a phase change of 0.14 radians within the InSAR frame according to Equation (6). Therefore, we discard two of the InSAR pairs from 13 August 2007-13 February 2008 and 30 September 2008-31 December 2008. The elevation difference caused by the L-and C-band SAR system is a spatial and temporal low-frequency signal in the LCED maps. Two LCED maps in Figure 8(a2,a3) are averaged into a final LCED map to reduce the impact from highfrequency components.  Figure 7(b1-b4). These represent elevation differences between ALOS/PALSAR L-band DEMs and SRTM C-band DEMs; (b1-b4) Coherence maps of the four InSAR pairs. The ALOS/PALSAR amplitude, SRTM C-band DEM, and the InSAR coherence map are shown in Figure 9a-c. The study site is classified into three groups and the classification result is shown in Figure 9d. Three groups of samples are chosen manually. Samples are divided in a proportion of 50% to 50% for training and testing the SVM classification. The overall accuracy of the classification is 97.2803%.  Figure 7(b1-b4). These represent elevation differences between ALOS/PALSAR L-band DEMs and SRTM C-band DEMs; (b1-b4) Coherence maps of the four InSAR pairs. The ALOS/PALSAR amplitude, SRTM C-band DEM, and the InSAR coherence map are shown in Figure 9a-c. The study site is classified into three groups and the classification result is shown in Figure 9d. Three groups of samples are chosen manually. Samples are divided in a proportion of 50% to 50% for training and testing the SVM classification. The overall accuracy of the classification is 97.2803%.  Figure 10a-c show the histograms of the LCED within each group along with normal distributions fitting to the histograms. From these histograms and normal distributions, we can observe increasing negative elevation differences in the plain and channel areas compared to those in the hilly area. Although some locations in the hilly areas show positive elevation differences, the mean value is still negative, as we can see in Figure 10a. The elevation difference in the hilly area is −0.983 ± 2.284 m, whilst in the plain area and channel area, it is −1.348 ± 2.463 m and −2.744 ± 3.244 m, respectively. The standard deviation is relatively large, which may result from the inclusion of many measurements of low coherence. However, we cannot totally rule these out because the penetration process is more likely to occur where the coherence is low. Thus, we only rule out measurements of coherence lower than 0.218. The larger skewness of the hilly area and channel areas indicate that more measurements tend to show negative shift values of LCED. To evaluate the mean value obtained by the histogram, independent sample numbers need to be considered. We calculated the correlation range of the LCED map, which is about 2.2 km in this study area. The calculation and statistical information is listed in Table 4. According to the calculation in Table 4, the standard errors of the mean are all within 0.1 m, which indicates that the shift we derived is reliable.  Although some locations in the hilly areas show positive elevation differences, the mean value is still negative, as we can see in Figure 10a. The elevation difference in the hilly area is −0.983 ± 2.284 m, whilst in the plain area and channel area, it is −1.348 ± 2.463 m and −2.744 ± 3.244 m, respectively. The standard deviation is relatively large, which may result from the inclusion of many measurements of low coherence. However, we cannot totally rule these out because the penetration process is more likely to occur where the coherence is low. Thus, we only rule out measurements of coherence lower than 0.218. The larger skewness of the hilly area and channel areas indicate that more measurements tend to show negative shift values of LCED. Figure 10a-c show the histograms of the LCED within each group along with normal distributions fitting to the histograms. From these histograms and normal distributions, we can observe increasing negative elevation differences in the plain and channel areas compared to those in the hilly area. Although some locations in the hilly areas show positive elevation differences, the mean value is still negative, as we can see in Figure 10a. The elevation difference in the hilly area is −0.983 ± 2.284 m, whilst in the plain area and channel area, it is −1.348 ± 2.463 m and −2.744 ± 3.244 m, respectively. The standard deviation is relatively large, which may result from the inclusion of many measurements of low coherence. However, we cannot totally rule these out because the penetration process is more likely to occur where the coherence is low. Thus, we only rule out measurements of coherence lower than 0.218. The larger skewness of the hilly area and channel areas indicate that more measurements tend to show negative shift values of LCED. To evaluate the mean value obtained by the histogram, independent sample numbers need to be considered. We calculated the correlation range of the LCED map, which is about 2.2 km in this study area. The calculation and statistical information is listed in Table 4. According to the calculation in Table 4, the standard errors of the mean are all within 0.1 m, which indicates that the shift we derived is reliable. To evaluate the mean value obtained by the histogram, independent sample numbers need to be considered. We calculated the correlation range of the LCED map, which is about 2.2 km in this study area. The calculation and statistical information is listed in Table 4. According to the calculation in Table 4, the standard errors of the mean are all within 0.1 m, which indicates that the shift we derived is reliable.

Validation of the Penetration Depth by ICESat/GLA14
ICESat has three transects in this study area, which are labelled as ICESat-T1, ICESat-T2, and ICESat-T3 in Figure 11a,b. The base map of Figure 11a is the classification result shown in Figure 9d. The base map of Figure 11b shows the elevation difference, which is averaged from Figure 8(b2,b3). Along these transects, we compared the LCED and elevation difference between SRTM C-band DEM and ICESat elevations. Figure 11(c1-c3) shows the comparisons. From Figure 11(c1-c3), we can see that in the hilly area, both the SRTM DEM and PALSAR DEM agree closely with the ICESat elevations. However, when it comes to the plain and channel areas, both the SRTM DEM and the PALSAR DEM show lower elevations than the ICESat elevations and the PALSAR elevations deviate even more than the SRTM DEM from the ICESat elevation. The red rectangle highlights the area where a shift of elevation difference can be noticed between SRTM DEM and PALSAR DEM. Over the channel area, the PALSAR DEM is noisy due to low coherence; however, it clearly shows that there are lower elevations than either the SRTM DEM or ICESat elevations.

Validation of the Penetration Depth by ICESat/GLA14
ICESat has three transects in this study area, which are labelled as ICESat-T1, ICESat-T2, and ICESat-T3 in Figure 11a,b. The base map of Figure 11a is the classification result shown in Figure 9d. The base map of Figure 11b shows the elevation difference, which is averaged from Figure 8(b2,b3). Along these transects, we compared the LCED and elevation difference between SRTM C-band DEM and ICESat elevations. Figure 11(c1-c3) shows the comparisons. From Figure 11(c1-c3), we can see that in the hilly area, both the SRTM DEM and PALSAR DEM agree closely with the ICESat elevations. However, when it comes to the plain and channel areas, both the SRTM DEM and the PALSAR DEM show lower elevations than the ICESat elevations and the PALSAR elevations deviate even more than the SRTM DEM from the ICESat elevation. The red rectangle highlights the area where a shift of elevation difference can be noticed between SRTM DEM and PALSAR DEM. Over the channel area, the PALSAR DEM is noisy due to low coherence; however, it clearly shows that there are lower elevations than either the SRTM DEM or ICESat elevations.

Discussion
This study is aimed at estimating the different penetration depths detected by L-band and C-band SAR. InSAR can measure the terrain variation by using phase differences between two or more acquisitions. When it is applied to L-band and C-band SAR datasets, it can measure different elevations. The accuracy of the elevation difference measured by the InSAR technique depends on internal errors, such as unwrapping errors and baseline errors, and external interruptions, such as atmospheric effects and surface displacement.
The unwrapping error causes a phase ambiguity which leads to a height ambiguity. For the ALOS/PALSAR datasets used in this study, the height ambiguity is 28 m. It indeed exists during InSAR processing on some pixels. However, from the unwrapping results shown in Figure 7(a1-a4) and Figure 7(b1-b4), we can see that most of the pixels in this dataset have been correctly unwrapped. Besides, we analysed the LCED by partitioning the regions and showing all measurements as histograms. There are outliers whose LCED is larger than 28 m in the LCED maps. Their impact on analysing the regional difference of LCED is minimised by using the histograms. The baseline error, which can lead to an error in the height measurement, can be expressed as ∆h = ∆B ⊥ ·h/B ⊥ . The height variation in this study site is about 400-600 m, thus, ∆B ⊥ of 10 m can cause a height error of 1 m for the ALOS/PALSAR dataset when the perpendicular baseline is 2000 m.
According to Section 4.2.3, the atmospheric phase screen simulated by using ECMWF data features primarily as signals of low spatial frequency. This low-frequency signal can be removed during the baseline refinement. Even though a residual of the atmospheric phase may be preserved in the final unwrapped phase which is converted to the LCED maps, its impact on the LCED is depressed by using long perpendicular baselines, according to the Equation (5). We have discussed the surface displacement in Section 3.2.3. There has not been any earthquake in this region, which could be a source of surface displacement. According to Section 4.1, areas appearing to be rippled sand dunes are observed in SRTM C-band image whilst unobserved in ALOS/PALSAR image. This indicates that sand movements are very likely to be transparent to a L-band SAR system. As for any groundwater withdrawal, it is highly unlikely to cause the elevation shift between the channel and non-channel areas, which is described in Section 4.3.
After considering the internal and external errors of InSAR processing and the observed elevation shifts observed between the channel and non-channel areas, we believe that the remaining phase in Equation (3) is mainly caused by the differences in penetration depths between the L-and C-bands. This difference results from their different penetration depths or different incident angles. When the SRTM DEM has an incidence angle of 60 • , the index of refraction needs to be larger than 80 to meet an LCED of 2 m, which is impossible. Otherwise, it may be explained by sand erosion from 2000 to 2008. However, this implies further studies on the locations and changes in the sand deposits and erosions, which is beyond the scope of this study. It may be discussed in the future if sand movement data become available for this study site.
In this study, we regard that the offset of the LCED between the channel and non-channel areas is mainly caused by different penetration depths as it has been reported in several locations in the eastern Sahara [2,4]. Moreover, the increasing negative elevation differences from hilly areas to plain areas, and to channel areas, which are shown in Section 4.3, support this argument. Since the possibility of penetration occurring in the hilly area cannot be excluded, the negative shift of the LCED is about 1.76 m between the hilly and channel areas, which should be the penetration depth if the elevation difference in the hilly area is zero. During InSAR processing, the SRTM DEM was taken as a reference. The comparison between SRTM DEM and ICESat/GLA14 elevations in Section 4.4 demonstrated that the SRTM DEM also detected a subsurface elevation rather than a surface elevation. Therefore, the actual penetration depth of the L-band PALSAR DEM would be the previous estimation of 1-2 m plus SRTM DEM's penetration depth. In conclusion, the ALOS/PALSAR InSAR derived DEM shows an elevation that is at least 1-2 m lower compared with the SRTM and ICESat elevations.
The penetration depth is key to studies related to terrain measurements in regions where there is ice or sand cover. It is especially crucial for the longer wavelength SAR mission, such as L-band ALOS/PALSAR and future P-band BIOMASS. In the future, the methods proposed in this study can be repeated in other deserts, such as the Mojave Desert in South California, USA, Gobi Desert, etc., to estimate the penetration depth. Given more SAR datasets, the phase components from sources like atmosphere and surface displacement can also be analysed in more detail. Currently, limited to a single wavelength (L-band) used by space-borne remote sensing data, a penetration depth of only a few metres can be estimated. It should be noted that the aquifers in the Eastern Sahara have been reported to be of about a 25-100 m depth below ground level in research conducted by the British Geological Survey [59]. Therefore, it is unfeasible to map any aquifer below the Sahara for the present time. However, in the future, low-frequency radar systems including the ESA BIOMASS P-band, as well as even lower frequency radars, may be able to be applied in this field [60], particularly if good corrections can be found for ionospheric distortions at the P-band.

Conclusions
This study presents a comparison between ALOS/PALSAR L-band SAR image to Landsat optical image and SRTM C-band amplitude image regarding different surface morphologies that are caused by deeper penetration depth detected by L-band SAR. An L-band and C-band elevation difference map was produced using InSAR and used to estimate the penetration depth in Kufra region. The results show that over one tributary of the paleochannel in Kufra region, the ALOS/PALSAR detects surfaces at least 1-2 m lower than the corresponding C-band SRTM.
The subsurface mapping and the estimation of penetration depth are limited by currently available SAR datasets to L-band frequency, so it is hard to detect some important subsurface features, such as aquifers in hyper-arid regions. However, this study indicates potential for future low frequency SAR missions such as ESA BIOMASS P-band SAR system, which are under development and scheduled to be launched in the early 2020s.