The physical properties of Mars polar deposits have been studied for several decades. Some studies show that North Polar Layered Deposit (NPLD) and South Polar Layered Deposit (SPLD) might be rich in water ice   . The stratified NPLD and SPLD media were formed due to varying amounts of dust impurity mixed with the water ice   . The varying impurity ratio is likely related to historical climate change . Techniques to study the dielectric properties of the regolith media in NPLD and SPLD are important for the study of Mars climate. One such technique is the inversion of dielectric constants using HF radar to penetrate through the regolith.
HF radar waves can penetrate through the Mars regolith media several kilometers. The MARSIS (Mars Advance Radar for Subsurface and Ionospheric Sounding) onboard the Mars Express operates in 4 bands centered at 1.8, 3, 4, 5 MHz with a bandwidth of 1 MHz, and can penetrate through the media as deep as 4 km . The HF radar sounder of the SHAllow RADar (SHARAD) onboard Mars Reconnaissance Orbiter (MRO) operates with a 20 MHz central frequency with a bandwidth of 10 MHz. Its vertical resolution is higher than MARSIS, and its penetration depth is usually less than about 1 - 2 km   (As water ice with a small loss tangent, the SHARAD signals may reach depths more than 2 km  ). China is planning to explore Mars multi-layer structure using HF and VHF radar as well  .
To study the radar sounder data, Mouginot et al.  adopted the MARSIS data (1 - 5 MHz) to retrieve global surface reflectivity with kilometer-scale surface roughness to estimate the Mars surface dielectric constant. Nouvel et al. , Lauro et al. , Mouginot et al.  particularly studied the surface dielectric constants of NPLD and SPLD. Grima et al. , Zhang et al.  , Alberti et al.  evaluated the average dielectric constant within 2 km depth of Mars polar deposits, using the radar echo time delays through the surface/subsurface of the regolith media. However, the study of the Mars cratered rough surface with multi-stratified interfaces, the development of the physical parameters inversion, and the validations using HF radar data are remained to be further studied .
Based on numerical simulations of the radar sounder echoes from the one-layer model with rough surface/subsurface media, Ye and Jin  found that under the Kirchhoff approximation with a mean zero slope, the received echo at nadir direction preserves the functional dependence of the surface reflectivity. It leads to the inversion of the surface dielectric permittivity derived from the ratios of the received echo powers and the medium reflectivity. Furthermore, Liu and Jin  proposed a numerical approach of radar echoes from rough surface/multi-subsurface and inversions.
Based on  , this paper presents a model of parallel-stratified media to describe the multi-layer structure of Mars polar region. A relationship between the received radar sounder echoes and reflectivities of rough surface/multi-subsurface is presented. The inversions of the thickness and dielectric permittivity of each layer where the radar wave can reach are designed. As data validation, inversions are applied to SHARAD data on Promethei Lingula of Mars SPLD.
2. Model of the Stratified Media and Inversion of Dielectric Constants
1) A Model of Parallel Stratified Media
The SHARAD radar sounder data has shown that there are multi-layer structures in Mars Polar regions . When the sub-layer thickness varies little within one radar footprint, the multi-layer structure from the top to hundreds meters below can be modeled as a parallel-stratified media, as shown in Figure 1. As the electromagnetic (EM) wave of radar sounder is vertically incident upon the top surface through ionosphere, there would be multi-reflection and transmission through the media. Based on the difference of time delay of each echo from media interfaces, and separating the rough surface clutter and the echoes of the interfaces, the echoes from each interface can be identified.
Suppose that multiple reflection and transmission between interfaces are neglected. It means that the echo from the n-th interface experienced one-reflection, round trip of 2(n − 1) transmissions (i.e. including round trip attenuation) through previous (n − 1) layer media. This assumption is based on small difference on final surface reflectivity caused by underlying multi-layer structures. Thus, as the incident radar power through ionosphere is directly on the top surface , the echo from the n-th interface is written as
where is denoted as the transmitted power of dipole antenna (with notation (a)), and similarly, denotes the n-th reflected peak power received by the radar antenna, i.e. observation. In Equation (1), is the transmittivity between the (m − 1)-th and the m-th media, is the reflectivity of the n-th layer, is the thickness of the m-th layer, is the imaginary part of the wave number of the m-th layer, is the round-trip attenuation in the m-th layer , is the one-way attenuation through ionosphere layer. The wave number of the m-th layer is written as
Figure 1. A model of rough surface and stratified media.
Here is the wave number of free space. and are the real and imaginary parts of the m-th layer dielectric constant, respectively. is written as
The transmittivity satisfies
and the layer thickness can be expressed as
where c is the light speed in free space, and is the time delay as EM wave propagates through the m-th medium twice, i.e. the time delay between two echoes from two successive interfaces.
Substituting Equations (3)-(5) into Equation (1), the reflected power from the n-th interface is derived as
where is the (attenuation) loss-tangent of the m-th layer.
2) Calculation of Loss Tangent
Since roughness of underlying interfaces is totally unknown, the model makes all underlying interfaces as plane-stratified, as shown in Figure 1. Equation (6) presents a set of total n equations, where and can be found from the radar data , but the incidence power directly upon the top surface , the attenuation , the reflectivities , and the loss tangents are to be solved. Considering the ionosphere attenuation is wrapped into , there are totally 2n independent unknowns in Equation (6).
It has been known from Mars studies that the loss tangents of NPLD and SPLD are actually very small, as usually 0.001 - 0.005  . Some extreme cases of the media with high conductivity (dipolar and conductive) are excluded. In radar sounder technology for exploring multi-layering structure, whole media should be at a very low loss to make the wave penetration to reach the interfaces.
Certainly, the reflectivity is a main factor to affect , comparing with and . To reduce the number of unknowns and reach final inversion, all small loss tangents of the n-layers are trivial and seen as the same within one illuminated area (e.g. 1.81 km in the next example). The low value of the loss tangent makes such approximation reasonable. Thus, Equation (6) is simplified as
Taking natural log of both sides of Equation (7), it gives
It means that the echo from each interface is a linear function of the time delay
where is the time delay from the n-th interface echo to the top
surface echo, b is an unknown constant, is a random variable to take account of different reflectivities, , of the interfaces.
Equation (9) can be seen as a linear regression model with the regress or . Using the radar range echoes from all interfaces and their respective ranges, the linear fitting is obtained with the least square method. Then the loss tangent can be calculated by the slope of linear function in Equation (9).
3) Solution of Dielectric Constant of Each Layer
Since the loss tangent is obtained, the number of unknowns now becomes n + 1. The set of Equation (7) can be directly solved, and the reflectivity is written as
Equation (10) can be solved, iteratively.
Reflectivity is a function of the dielectric constants of layering media. Because the roughness of all sub-interfaces is totally unknown, it is practicable to model all sub-interfaces between layers as flatly stratified within a limited area. This is a good and workable assumption within a limited area illuminated by the radar waves, even if a small error might be caused due to small scale roughness of the interfaces.
It is also noted that in derivation of Equation (1), multiple-reflection and transmission are neglected. Thus, the reflectivity from the (n − 1)-th layer to the n-th layer is derived based on a half-space model, i.e. the reflectivity of this interface is written as
Since the loss tangent is very small, it yields , Equation (11) becomes
But, there would be two solutions from Equation (12) due to the term . Using the phase change of the echoes, an unique solution may be obtained.
The echo phase from the n-th interface is written as
where is the phase of EM incidence, denotes the phase from each transmission, and the phase of each reflection.
As EM wave is vertically incident from the (n − 1)-th layer to the n-th layer, the reflection coefficient and transmission coefficient are, respectively, written as
Since the loss tangent is very small, it can be seen that if , it makes and ; otherwise, if , it makes and . In transmission, the phase keeps unchanged, i.e. .
As incident upon the top surface, . Equation (13) gives . Thus, all phases due to reflections from all interfaces can be calculated from the data of radar range echoes as
Based on these approximations, it yields the dielectric constant of each layer as
4) Calculation of Ionospheric Attenuation and Dielectric Constant of the Surface Medium
Propagation through the ionosphere causes phase distortions and attenuation. The SHARAD Reduced Data Record (RDR) data has already corrected the phase distortion using Phase Gradient Autofocus (PGA) method . The ionosphere is excited by solar radiation, and its effect during daytime is much stronger than during night   . To avoid the additional error caused by phase distortion, only data acquired during the night are specifically used in the following example. Moreover, ionospheric attenuation is taken into account using an uniform constant when the SHARAD data acquired with the similar solar zenith angle (SZA).
Since the SHARAD data, as available, have not been absolutely calibrated , the simulated echo power from the interfaces used in the inversion  , , should be adjusted to match the observations on the radar receiver. It gives
where calibration constant can be seen to take account the one-way ionospheric attenuation.
It has been studied  that there is purely CO2 ice covering an area of the South Pole of Mars, and its dielectric constant is known as about 2.2. Thus, based on the observation data , as available, and our simulated data at this South Pole location, C of Equation (18) is obtained and applied to the whole inversion region.
The top surface is modeled as a rough surface, described by the known DEM data. From the radar equation , the echo from the top surface (with ) under radar EM wave incidence (nadir incidence ) is written as
where is the backscattering coefficient of rough surface.
The Kirchhoff approximation (KA) of rough surface scattering requires the curvature radius of the surface much larger than the radar wavelength , it is a gentle rough surface on most areas on Mars. It has been discussed  that based on derivations of rough surface scattering with surface mean zero slope and the KA, the received echoes power at nadir direction is simply proportional to the surface reflectivity. It presented the inversion that the ratio of the received echo powers from one unknown and another tested surface permittivity can present inversions of the surface reflectivity, and dielectric unknown.
Thus, it is derived as 
where is the incident wave vector, is the distance vector from the pixel center of integral to the nadir point.
Suppose that the top surface has a test value , the simulation gives . The ratio of the observation with an unknown over the simulation with assumed gives 
where C was defined in Equation (18), and actually can be evaluated in the next approach. Substituting Equation (11) into Equation (21), it gives the inverted , as follows
It is noted that this inversion is applicable for surfaces with gentle roughness and zero mean slope. Indeed, there are large areas on the SPLD to fit this description as shown from DEM data. Our approach focuses on those cases and ignores those cases with steep slopes or highly varying roughness.
Using the inverted and Mars Orbiter Laser Altimeter (MOLA) elevation data, the backscattering coefficient can be calculated in our numerical simulation  . It yields the incident power on the top surface, of Equation (10) as follows,
Substituting the inverted and observation into Equations (10), it gives . Then, Equation (17) gives . Sequentially, it yields , etc. The thickness of each layer can be then calculated by Equation (5).
3. Inversion of Dielectric Constants at the Promethei Lingula of Mars South Polar Region
1) SHARAD Radar Echoes Data from the Promethei Lingula
As shown in Figure 2, the Promethei Lingula around Mars South Pole is a part of SPLD. Figure 3 presents one SHARAD observation over this region. The radargram in Figure 3 is composed of about 2500 echoes after pre-summing on board the spacecraft. The radar flight distance is about 90.5 km during 2500 times vertical sounding. A stratified structure can be well identified in SHARAD radargram. Figure 4 is the optical photo of High Resolution Image Science Experiment (HiRISE) around the cliff of Chasma Australe canyon on the west edge of the Promethei Lingula, which clearly shows the existence of stratifications of layering media .
Figure 2. Location of SHARAD track 17485_01 data in Promethei Lingula on a MOLA elevation map.
Figure 3. A portion of SHARAD track 17485_01 data from stratified media in Promethei Lingula.
Figure 4. The stratified structure on HiRISE photo of Chasma Australe canyon eastern cliff. (ESP_023590_0975_RED.abrowse.jpg on HiRISE website).
We choose the track 17485_01 of SHARADRDR data on PDS Geosciences node (filename: r_1748501_001_ss11_700_a.dat on website http://pds-geosciences.wustl.edu/), which passes by Promethei Lingulanear the Chasma Australecanyon during the night (SZA is about 112˚), for inverting the parameters of multi-layer media. The vertical resolution of the data is 15m in vacuum and about 8.5 m in pure ice. The latitude and longitude of each frame are indicated in the SHARAD RDR data, and the along-track distance between each frame can be calculated. Especially, the distance between each frame in Figure 3 is about 36.2 m (0.00056˚ in latitude and 0.002˚ in longitude). According to the radargram and optical image, the topography and interfaces below the top surface look almost flat for one radar footprint. It validates our flatly stratified media model for Promethei Lingula.
2) Echoes from the Surface/Sub-Surfaces and Interface Locations
In the radargram, the nadir surface echo, off-nadir clutters of rough surface, and echoes from layering interfaces must be identified and treated, separately. While the strongest peak is often from the nadir echo, the surface roughness and DEM geometry can make bright off-nadir returns . Here, we specifically choose Promethei Lingula, where whole surface slopes are rather small and make the brightest echoes from nadir direction. However, it is always difficult to separate the echoes from underlying interfaces and clutter of rough top surface, because both amplitudes might be in the same order. As the radar is moving, approaching or leaving, the clutters of surface roughness might show variation as a point or short arc in radargram. In contrast, the time delay of the echoes from flat interfaces has little change as radar is moving, and usually shows straight line in radargram. Based on this intuition, we define a threshold S to decide if the echo is from the subsurface or not:
For example, let n = 25 and m = 1, and hence S indicates the ratio of local maximum from totally nearby 50 frames with the similar time delay (not exceeding 1 sampling interval). If , it means that more than 70% of adjacent frames have reflector with the same time delay, and the reflector is judged to be from the interface. Otherwise, the isolated reflector is seen as the surface clutter. In this way, the surface echoes and interface echoes are distinguished from frame 69,301 - 69,350 of the track 17485_01, which extends about 1.81 km.
Finally, the stratification of multi-interfaces is shown in Figure 5. Not all the reflectors in Figure 5(a) are kept in Figure 5(b). Some short lines and bright points are erased, and only are kept those long lines. If a few adjacent pixels on the long line happen to be weak and break the continuity of the extracted interface, the weak pixels are artificially fixed and judged as interface to avoid layering discontinuity. Linear interpolation of the interface line is used to decide the interface ranges of weak pixels. And the power of the nearest pixel is judged as the echo power of the interface.
Thus, the surface echoes, the echoes from underlying interfaces, i.e. , can be obtained. Sometimes, the echoes from different locations of the same interface might be quite different. It might be caused by different interface- topography, or happens to be mixed by the surface clutters. Dimmer or brighter radar echoes may also be caused by the change of interface reflectivity or the change of the interface time delay, for the time delay changes lead to different interference in the radar signal. To avoid such fluctuations of the interface echoes to affect final inversion, the echoes from the same interface is taken as an averaged value.
Figure 5. Location of multi-interfaces extracted from the SHARAD data. (a) SHARAD radargram from track17485_01; (b) stratification Location of multi-layers.
3) Dielectric Constant of the Surface Medium
Figure 6(a) shows the surface radar echoes power acquired from SHARAD in South Polar Region, where the central white circle is no-data region. The surface echo (the strongest peaks in each frame) from each area is taken from SHARAD data are used to generate Figure 6(a), which covers most part of SPLD. Similar figures had been done by Grima et al.  (Mars surface reflectivity map of SHARAD) and Mouginot et al.  (Mars surface reflectivity map of MARSIS).
Using the MOLA (Mars Orbiter Laser Altimeter) elevation data, the surface echoes from this rough surface can be numerically simulated  . At the beginning, we take a proposed dielectric constant over whole area, which is similar to the water ice. And is used for the echoes simulation around Mars South Polar region. The simulated is obtained as shown in Figure 6(b).
shown in Figure 6(c). It has been studied  that there is thick purely CO2 ice cap covering the South Pole, which can be also seen from blue colors of Figure 6(a) and Figure 6(c). The dielectric constant of CO2 is much lower than water ice and rocks. So the echo power of CO2 cap is much less than other places. The constant of C, Equation (18), is actually obtained from these figures.
Figure 6. (a) The SHARAD surface echoes data; (b) simulated echoes with dielectric constant 3; (c) ratio of SHARAD echoes/simulated echoes; (d) dielectric constant of Mars South Polar Region.
Using the algorithm aforementioned, the dielectric constant of Mars surface over the South polar region, , can be inverted, as shown in Figure 6(d). The inverted results are influenced by the seasonal CO2 ice, because the data used in the inversion are all acquired during autumn and winter, and the Mars polar regions are covered by a thin CO2 ice layer less than 1 m thick at that time . Some studies show that the radar echo power decreases when a thin CO2 ice covers on the top player . Since the covering of seasonal CO2 ice is variable and uncertain for data acquired in different time, which is much thinner in terms of SHARAD resolution, the inverted dielectric property of the first layer can be seen as the effective dielectric constant including thin CO2 ice layer as available. The inverted results are the equivalent dielectric constants of a cluster of thin layers comprised of CO2 ice and water ice with dust. The of most areas of SPLD is between 3 and 4. These results show that SPLD is mainly comprised of water ice, which is consistent with previous studies  . However, there are still some areas of SPLD in Figure 6(d) having larger than 4, including the Promethei Lingula. It can be seen that the areas with large are all located in the outer part of SPLD, which does mean that these areas may contain more impurity than the center of SPLD. Certainly, low resolution of MOLA data (except the data near the pole, which has a high resolution) as available, to describe complex rough surface might cause inversion error. The typical Root Mean Square heights over the SPLD at SHARAD scales is ~0.30 m , which is smaller than the MOLA resolution of 1 m, which might cause an error in inversion.
4) Calculation of Loss Tangent
To reduce the fluctuation of different reflectivities of the interfaces ( in Equation (9)), the echoes powers of the interfaces with the same time delay are averaged. Using the least square method to make a linear fitting for the average echoes power with different time delay, a linear equation of is obtained, as shown in Figure 7. From Equation (9), it yields . The confidence interval of loss tangent with the confidence level of 95% is [0.0004, 0.0014].
Figure 7. Relationship between each echo power and its time delay.
Figure 8. Dielectric constant of each layer (along track).
Figure 9. Dielectric constant of each layer (across track).
Figure 10. Dielectric constant range of each layer.
The hypothesis of linear regression model to well fit the data is tested using F distribution. It takes for total 42 data points in the linear fitting. Setting the significant level , it gives and . Thus, the linear regression is good to fit the data. The loss tangent is much smaller than previous result  , which contained both absorption of martials and signal loss at each interface reflection , while in our inversion, the signal losses of the interface reflections are totally removed, as shown in Equation (8).
5) Echoes Phase Estimation and Unique Solution Determination
The sampling interval of the SHARAD data is 0.075 μs , which corresponds to the depth of 0.75 wavelength in free space according to Equation (5). It would be difficult to obtain enough accurate phase information as the radar wave propagates through the layers. Firstly, we interpolate each frame data using sinc function . Then, the local maximum of the interpolated data is judged as the new interface instead of the interface extracted before. Thus, the time delay of each interface can be estimated accurately. Then, the phases are calculated by Equation (16). Due to surface clutter interference, there are noises in the phases which make errors during the phase evaluation. We take an average of the phases of echoes from the same interface to reduce the errors.
where is the phase of frame f and n-th interface.
Substituting Equation (25) into Equation (16), the phase of reflection is calculated. Considering that usually is not exactly equal to 0 or π due to the phase accuracy, Equation (17) is changed as
6) Inversions and Validation
The dielectric constant of surface medium inverted from the frame 69,301 - 69,350 (83˚S 102˚E in Figure 6(d)) of the track17485_01 is about 5 (see the location of Layer 1 at elevation 2100-2200 in Figure 8 and Figure 9. Note that the layer 0 on the top of Figure 8 and Figure 9 is Mars atmosphere). This result may be overestimated as discussed in part C of section II. However, it is possible that Layer 1 of Promethei Lingulacontain more impurity than the center of SPLD. Thus we choose .
Substituting #Math_130# into Equation (25), and changing the loss tangent from 0.0004 to 0.0014, the varying range of the layered dielectric constants are iteratively calculated, as shown in Figure 10. Since the loss tangent is quite small, the error caused by using the same loss tangent is not larger than 25% for the inverted dielectric constants in Figure 10. The dielectric constant of Layer 7 is 2.5, which is much lower than ice. As the error is taken into account, this layer might be seen as pure ice.
Using Equation (5), the thickness of each layer is calculated, as shown in Figure 8. There won’t be echo from the place below Layer 9, since EM wave therehas become very weak due to transmission and reflection. Thus, the dielectric constant below 1500 m in Figure 8 is actually unknown, and the thickness of layer 9 cannot be inverted.
The inversion accuracy depends on the evaluation of interface locations and echo powers. The phase accuracy might be also interfered by surface clutter.
Suppose the depth of each interface do not change across the track, as shown in Figure 9. The red point in Figure 9 is the nadir of track 17485_01 and the dashed line is the stratification in Figure 8. Then the position of each interface exposed on the cliff can be calculated using MOLA elevation. To validate our inversions, the calculated interface exposed on the cliff is projected to a map-projected optical image, as shown in Figure 11. Note that Figure 8 gives the vertically layering structure, while Figure 11 shows a slant section of the cliff photographed vertically downwards
It can be seen that the inverted multi-layering structure is well described on the cliff of Chasma Australe canyon near track 17485_01, which was indicated in a HiRISE optical image, even not exactly the same matching. The optical layers are much finer than the radar resolution, and they cannot be a 1:1 correlation. The layers in radargram correspond to the packets of thinner layers in optical image. Many layers can be seen on optical image but cannot be detected by radar, because this HF radar technology is capable only to detect the interfaces with a significant change of dielectric media. The inversion model presented is based on the implicit assumption that there is no more than one such interface within a SHARAD vertical resolution (10 - 15 m). If SHARAD reflections are caused by merged reflections of packets of thinner layers with different dielectric constants, the inversed layer-depths and dielectric properties are understood on the average or effective sense for radar sounder echoes. Moreover, details of the technology parameters and measurements, such as SNR for each echo and high resolution elevation data might further improve the inversions.
Figure 11. Cliff stratification of layering media overlapped on a HiRISE photo (a portion of PSP_006264_0970_RED.JP2).
The inversion results show that the dielectric deposits in Promethei Lingula are not uniformly stratified, seen along the dashed line indicated in Figure 9. Most layers seem nearly pure ice (e.g. blue colors), but some thin layers including the surface layer are with a lot of dust (e.g. yellow colors). As the deposit is seen as a component of ice and dust, its dielectric constant is written as 
where is the dust fraction, and can be calculated as
Taking the dust basalt ( ), the dust fraction is calculated using Equation (28), as shown in Figure 11. Layer 3 and Layer 5 have a dust load of 60% - 70%, which means these layers may be a compact frozen ground (permafrost) rather than dirty ice. However, SHARAD reflections might be caused by thinner layers, on the order of a meter thick or so  . Packets of low dust thin layers with a few thin permafrost layers might be another possible existence of Layers 3 and 5. Thus, it is possible to overestimate the dust fraction in the inversion.
Making a dielectric average of top 8 layers media, it yields
It gives the dielectric constant of regolith impurity, , which corresponds to the dust fraction about 12%.
This paper presents a model of rough surface and stratified interfaces to describe the multi-layer structure of Mars Polar Layered Deposit. The range echoes of HF radar sounder from rough surface/subsurface is numerically calculated. And under the Kirchhoff approximation with a mean zero slope, the received echo at nadir direction preserves the functional dependence of the surface reflectivity. In radargram to show the radar range echoes, the nadir surface echo and the echoes from layering interfaces are separated from off-nadir surface clutters. As the surface dielectric constant is derived from the ratio of the received echo peak power, the inversion approach is designed to obtain the dielectric constant and layer thickness of next layers, sequentially.
As a validation example, the SHARAD data are adopted to inversions of the Promethei Lingula stratified media of Mars South Polar region. The vertical profile of the dielectric constants of layering media is obtained. Correspondingly, the layered structure is obtained, which is visually similar to the optical image on the cliff of Chasma Australe canyon. The inversion results show that the surface media of Promethei Lingula has larger dielectric constant with impurity, while as the media below the surface might contain much less impurity and more pure water ice. As more details of the technology parameters and measurements can be taken into account, the inversion accuracy can be further improved.
The inversion model is based on the assumption that there is no more than one interface within a SHARAD vertical resolution. As SHARAD reflections might be caused by packets of thinner layers, multi-layer model under this resolution won’t be able to see them, individually, which are merged with other stronger adjacent reflections. The inversed layer-depths and dielectric properties are understood on the average or effective sense for radar sounder echoes.
This work was supported by the National Key Research and Development Program of China 2017YFB0502703.
The MOLA elevation data and SHARAD Reduced Data Record (RDR) data are all from PDS Geosciences Node (website: http://pds-geosciences.wustl.edu/). The HiRISE optical data are from the HiRISE website (http://hirise.lpl.arizona.edu/).