SHARAD Observations of Temporal Variations of CO2 Ice Deposits at the South Pole of Mars

Mars’s polar regions are covered by kilometers-thick layered deposits which carry a record of the planet’s climate history. The deposition and volatilization of the shallow CO2 deposits in the south pole have a large impact on the planet’s atmosphere and environment. This research focuses on the timing variation of the thickness of the shallow deposits based on the SHARAD data collected from the past 11 terrestrial years, and analysis of the contributing factors based on the volatilization and deposition mechanisms of surface and subsurface materials. In this work, we selected more than four thousand data points, covering several seasons and Martian years, to extract radar echoes and calculate the thickness changes in the subsurface layer over time. We found that the thickness of the CO2 layer becomes thinner in the summer, with seasonal variation in the range of ~16–45 m. The thickness variations have a Gaussian-like distribution and do not increase with the distance between the compared node pair, implying that the phenomenon is not caused by regional differences. The overall thickness within the 11 terrestrial years does not show a clear trend of thickening or thinning, indicating a moderate vertical change of the southern deposits.


Introduction
The Martian polar regions are covered by a large volume of sediments, mainly composed of water ice and dust, which hold a record of past climate variations on Mars [1]. The sediments are in direct contact with the atmosphere, so their volume variation affects the pressure, atmospheric composition, and temperature of the planet. The layered structure and the temporal variation of the polar sediments impose important constraints on the atmosphere modeling and understanding of Mars's climate history.
The surface of the south pole is covered by a layer of seasonal CO 2 frost. The sublimation and deposition of the seasonal cap results in periodic changes in the Martian atmosphere. Underlying the seasonal frost is the South Pole Residual Cap (SPRC). The depression features on the surface of the SPRC erode at a constant rate in the summer, releasing CO 2 into the atmosphere. The South Polar Layered Deposits (SPLD) are beneath the SPRC and dominate Mars's south polar ice cap [2]. The Shallow Radar (SHARAD) [3] instrument from the Mars Reconnaissance Orbiter (MRO) can penetrate the southern polar deposits and reveal the internal layered structure ( Figure S1). It has unveiled a substantially large amount of CO 2 -ice layered deposits in the AA 3 unit which could double the atmospheric pressure on Mars. The CO 2 -ice deposits are preserved in some of the scarps and troughs within the SPLD formed during periods of strong ablation [4]. Although CO 2 ice is sheltered by the thin water-ice layer, the depression and pit features on the surface of AA 3 indicate sublimation of the internal CO 2 .
The deposition or volatilization of the internal CO 2 -ice layers could be controlled by the periodic changes of obliquity, perihelion, and eccentricity of Martian orbits, which

Temporal Variations of South Polar Cap
The Martian atmospheric pressure has seasonal and periodic changes due to the removal and deposition of the seasonal ice cap [12]. The CO2 deposits on the surface become thicker and thinner due to frost in winter and sublimation in summer, respectively [13]. CO2 frost reaches maximum accumulation at solar longitude, Ls = 155°. The solar  [8,9]. Black dots in AA 4a and BL 3  The AA 2 unit underlying AA 3 consists of evenly bedded planar sequences, and unconformably overlies most surfaces of unit AA1, which is the main component of SPLD with maximum thicknesses of >3.5 km. The layering structure revealed by SHARAD is suggestive of the erosion and deposition of Martian polar ice caps that interacted with midlatitude sediments [10,11], which are accompanied by unconformable contacts between different units.

Temporal Variations of South Polar Cap
The Martian atmospheric pressure has seasonal and periodic changes due to the removal and deposition of the seasonal ice cap [12]. The CO 2 deposits on the surface become thicker and thinner due to frost in winter and sublimation in summer, respectively [13]. CO 2 frost reaches maximum accumulation at solar longitude, Ls = 155 • . The solar longitude (Ls) is the Mars-Sun angle, and for the southern hemisphere, the autumn equinox is Ls = 0 • , followed by the winter solstice at Ls = 90 • , spring equinox at Ls = 180 • and southern summer equinox at Ls = 270 • . MOLA data indicate that the seasonal variation of surface elevation predominantly causes the seasonal cycle of CO 2 condensation and sublimation to be less than 1 m [14,15]. Generally, the SPRC has been balanced with the atmosphere, so it is sensitive to changes in the Martian climate [12,16]. Its temporal variations can be observed from MRO Context Camera (CTX) images, which provide context with~6 m/pixel resolution. The quasi-circular depressions known as "Swiss Cheese" [1] features appear on the surface of SPRC (Figure 2a). In the summer when CO 2 frost is removed, the walls of the depressions composed by an 8 m thick CO 2 -ice layer expand outward with a constant erosion rate [17], as shown in Figure 2b,c. Based on the observed erosion rates and the size of the erosion pits, studies estimated that SPRC has existed on Mars for a century [17,18]. The area of SPRC did not fluctuate significantly during the long period of erosions of pits and slopes [19,20].
Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 23 longitude (Ls) is the Mars-Sun angle, and for the southern hemisphere, the autumn equinox is Ls = 0°, followed by the winter solstice at Ls = 90°, spring equinox at Ls = 180° and southern summer equinox at Ls = 270°. MOLA data indicate that the seasonal variation of surface elevation predominantly causes the seasonal cycle of CO2 condensation and sublimation to be less than 1 m [14,15]. Generally, the SPRC has been balanced with the atmosphere, so it is sensitive to changes in the Martian climate [12,16]. Its temporal variations can be observed from MRO Context Camera (CTX) images, which provide context with ~6 m/pixel resolution. The quasi-circular depressions known as "Swiss Cheese" [1] features appear on the surface of SPRC (Figure 2a). In the summer when CO2 frost is removed, the walls of the depressions composed by an 8 m thick CO2-ice layer expand outward with a constant erosion rate [17], as shown in Figure 2b,c. Based on the observed erosion rates and the size of the erosion pits, studies estimated that SPRC has existed on Mars for a century [17,18]. The area of SPRC did not fluctuate significantly during the long period of erosions of pits and slopes [19,20]. Sublimation and erosion also occur in the AA3 unit. The polygonal patterns formed by highly fractured water-ice layer BL3 are unique features found at the surface of AA3, which might allow the buried CO2 ice to volatilize and form depressions. The troughs, Sublimation and erosion also occur in the AA 3 unit. The polygonal patterns formed by highly fractured water-ice layer BL 3 are unique features found at the surface of AA 3 , Remote Sens. 2022, 14, 435 5 of 21 which might allow the buried CO 2 ice to volatilize and form depressions. The troughs, depressions, and pits observed in the AA 3 unit have a depth of~75 m, indicating that some areas have internal CO 2 volatilization and release [9]. SHARAD observations at different times can help us understand the temporal variation of the AA 4 and AA 3 units.

Study Regions
Three study regions indicated by the white squares in Figure 3a were determined based on two criteria. The first criterion was to select the region containing subsurface CO 2 deposits within unit AA 3 . Its layering structure can be clearly observed in radargram. The AA 3 unit overlaps the AA 4 horizontally [9] in several instances. Therefore, we selected S1-S3 in Figure 3b as our study area, which is mainly located at 93 • -101 • W in longitude, 86 • 46 30 -87 • 02 30 S in latitude and is~108 km 2 in total area. The side length of the three regions is~6 km, respectively. Each region is within the horizontal resolution of SHARAD to obtain the time difference of the layer thickness at the same location. The parameters on SHARAD will be described in 2.1.2. The second criterion was to select SHARAD data with clear horizontal reflectors, so an area with a smooth surface was selected to minimize the measuring uncertainties of layer thickness caused by roughness. Figure 4 shows the CTX images of three study areas mainly located on the geologic units B6 and B7, where eroded landforms such as 2 m sized pits are commonly seen. The erosion features in B7 are mostly isolated mesas and low scarps, with few remaining pits, so surface roughness has little impact on radar signal in this area. In this region, surface clutter is barely seen, providing a clear view of surface and subsurface echoes used for data analysis. S2 belongs to B6 [21], the surface of which is coarser than those of S1 and S3 ( Figure 4), with numerous pits of a depth of~2 m [21]. The smooth surface of S1 and S3 has an advantage for improving the accuracy of the layer thickness, as measurement uncertainties are sensitive to roughness, and temporal variations are close to SHARAD resolution. depressions, and pits observed in the AA3 unit have a depth of ~75 m, indicating that some areas have internal CO2 volatilization and release [9]. SHARAD observations at different times can help us understand the temporal variation of the AA4 and AA3 units.

Study Regions
Three study regions indicated by the white squares in Figure 3a were determined based on two criteria. The first criterion was to select the region containing subsurface CO2 deposits within unit AA3. Its layering structure can be clearly observed in radargram. The AA3 unit overlaps the AA4 horizontally [9] in several instances. Therefore, we selected S1-S3 in Figure 3b as our study area, which is mainly located at 93°-101°W in longitude, 86°46′30′'-87°02′30′'S in latitude and is ~108 km 2 in total area. The side length of the three regions is ~6 km, respectively. Each region is within the horizontal resolution of SHARAD to obtain the time difference of the layer thickness at the same location. The parameters on SHARAD will be described in 2.1.2. The second criterion was to select SHARAD data with clear horizontal reflectors, so an area with a smooth surface was selected to minimize the measuring uncertainties of layer thickness caused by roughness. Figure 4 shows the CTX images of three study areas mainly located on the geologic units B6 and B7, where eroded landforms such as 2 m sized pits are commonly seen. The erosion features in B7 are mostly isolated mesas and low scarps, with few remaining pits, so surface roughness has little impact on radar signal in this area. In this region, surface clutter is barely seen, providing a clear view of surface and subsurface echoes used for data analysis. S2 belongs to B6 [21], the surface of which is coarser than those of S1 and S3 ( Figure 4), with numerous pits of a depth of ~2 m [21]. The smooth surface of S1 and S3 has an advantage for improving the accuracy of the layer thickness, as measurement uncertainties are sensitive to roughness, and temporal variations are close to SHARAD resolution. From top to bottom are S1, S2, S3, respectively. The red and black points are the trajectories of the radar data. The red points are the data used for analysis, and the black points are the data that were discarded. From top to bottom are S1, S2, S3, respectively. The red and black points are the trajectories of the radar data. The red points are the data used for analysis, and the black points are the data that were discarded. Figure 4. CTX images of three study areas: (a) S1; (b) S2; (c) S3. CTX image B06_012061_0930_XI_87S093W. Red dots represent SHARAD data used for analysis, and the black dots represent SHARAD data that were viewed but not used due to rough surface.

Shallow Radar Data
SHARAD instrument carried by the MRO began to work in October 2006, collecting radar echoes reflected from the Martian surface and subsurface layers. Its frequency is 15-25 MHz, and its bandwidth is 10 MHz. SHARAD can penetrate about 1 km beneath the Martian surface [8]. The vertical resolution of SHARAD in the subsurface layer can reach 10 m to 20 m [9] and has a cross-track resolution of 3-6 km and an along-track resolution of 0.3-1 km [3].
More than 6000 data points from SHARAD have been viewed for this study area. These cover Martian years 28-34, or an observation period from 2007 to 2018. The red dots in Figure 4 show the 4862 points with clear surface and subsurface echoes, which were selected for analysis. A rough surface could cause the delay offset and increase the halfheight width of SHARAD echoes. To reduce the influence of roughness on the data uncertainties, we set the upper threshold value of FWHM (full width at half maximum) to filter radar data. In this step, the radar echo trace of each subsatellite point was extracted and the corresponding pulse width measured. Data points, both surface and subsurface echoes of which have a pulse width of less than 150 ns, were used for subsequent statistical analysis. Figure 5 illustrates the depth of subsurface echo from the bottom of the AA3c and the pulse width before and after filtering. CTX images of three study areas: (a) S1; (b) S2; (c) S3. CTX image B06_012061_0930_XI_87S093W. Red dots represent SHARAD data used for analysis, and the black dots represent SHARAD data that were viewed but not used due to rough surface.

Shallow Radar Data
SHARAD instrument carried by the MRO began to work in October 2006, collecting radar echoes reflected from the Martian surface and subsurface layers. Its frequency is 15-25 MHz, and its bandwidth is 10 MHz. SHARAD can penetrate about 1 km beneath the Martian surface [8]. The vertical resolution of SHARAD in the subsurface layer can reach 10 m to 20 m [9] and has a cross-track resolution of 3-6 km and an along-track resolution of 0.3-1 km [3].
More than 6000 data points from SHARAD have been viewed for this study area. These cover Martian years 28-34, or an observation period from 2007 to 2018. The red dots in Figure 4 show the 4862 points with clear surface and subsurface echoes, which were selected for analysis. A rough surface could cause the delay offset and increase the half-height width of SHARAD echoes. To reduce the influence of roughness on the data uncertainties, we set the upper threshold value of FWHM (full width at half maximum) to filter radar data. In this step, the radar echo trace of each subsatellite point was extracted and the corresponding pulse width measured. Data points, both surface and subsurface echoes of which have a pulse width of less than 150 ns, were used for subsequent statistical analysis. Figure 5 illustrates the depth of subsurface echo from the bottom of the AA 3c and the pulse width before and after filtering. the width of the echo at 3dB before and after filtering the data, respectively. The red line represents the width of the surface echo, and the black line represents the width of the subsurface echo. Corresponding data in S1 and S3 can be found in the TABLE( Figure S2).

Measurement of CO2 Deposits' Thickness
The thickness measurement is illustrated in Figure 6. When electromagnetic waves emitted by the radar enter the Martian crust, the speed of propagation in the first layer v is reduced relative to the vacuum, here = / . n is the refractive index of the medium, associated with the real part of the dielectric constant , n = √ . If a subsurface layer with a discontinuous dielectric constant is present at a depth ∆d2 below the surface layer, the reflection and transmission of electromagnetic waves will occur again at the interface.
The secondary reflection echo penetrates layer one and should be received by the radar, although the signal is much weaker than the surface echo. The time delay ∆t2 between the subsurface echo and surface echo can be used to calculate the thickness of layer one ∆d2= ∆ 2 /2 √ 1 ( Figure 6).
Suppose that the surface of Mars is smooth and flat, only echoes of sub-satellite nadir point A can be received. However, in a real case, radar might also receive the reflection from C on the surface ( Figure 6). If its time delay is the same as the echo of the subsatellite point on the subsurface B, then it is likely to submerge the signal of B, which brings the problem of clutter. Clutter simulation, which uses a digital elevation model (DEM) to model the reflections from the non-nadir surface, is necessary to identify the subsurface echo. Figure 1a shows the stratigraphic map in the study areas. Since the thickness of surface layer AA4 is about 10 m, less than the vertical resolution of SHARAD, the lower interface of AA4 cannot be clearly seen in radargram. Therefore, the first layer recognized in the radargram includes AA4 if present, Bl3, and AA3c, meaning the measurement results given in this work are the overall thickness of three sublayers. For convenience, we call the first layers AA3c. We took the dielectric constant of 2.2, a typical value for the mixture of water ice and dry ice. The temporal change of thickness within seven Martian years in each small area was calculated and is shown in the results section. (c,d) the width of the echo at 3dB before and after filtering the data, respectively. The red line represents the width of the surface echo, and the black line represents the width of the subsurface echo. Corresponding data in S1 and S3 can be found in the TABLE ( Figure S2).

Measurement of CO 2 Deposits' Thickness
The thickness measurement is illustrated in Figure 6. When electromagnetic waves emitted by the radar enter the Martian crust, the speed of propagation in the first layer v is reduced relative to the vacuum, here v = c/n. n is the refractive index of the medium, associated with the real part of the dielectric constant ε r , n = √ ε r . If a subsurface layer with a discontinuous dielectric constant is present at a depth ∆d 2 below the surface layer, the reflection and transmission of electromagnetic waves will occur again at the interface.
The secondary reflection echo penetrates layer one and should be received by the radar, although the signal is much weaker than the surface echo. The time delay ∆t 2 between the subsurface echo and surface echo can be used to calculate the thickness of layer one Figure 6). Suppose that the surface of Mars is smooth and flat, only echoes of sub-satellite nadir point A can be received. However, in a real case, radar might also receive the reflection from C on the surface ( Figure 6). If its time delay is the same as the echo of the subsatellite point on the subsurface B, then it is likely to submerge the signal of B, which brings the problem of clutter. Clutter simulation, which uses a digital elevation model (DEM) to model the reflections from the non-nadir surface, is necessary to identify the subsurface echo. Figure 1a shows the stratigraphic map in the study areas. Since the thickness of surface layer AA 4 is about 10 m, less than the vertical resolution of SHARAD, the lower interface of AA 4 cannot be clearly seen in radargram. Therefore, the first layer recognized in the radargram includes AA 4 if present, Bl 3 , and AA 3c, meaning the measurement results given in this work are the overall thickness of three sublayers. For convenience, we call the first layers AA 3c . We took the dielectric constant of 2.2, a typical value for the mixture of water ice and dry ice. The temporal change of thickness within seven Martian years in each small area was calculated and is shown in the results section.

Surface Roughness
To analyze the influence of surface roughness on uncertainties of measurements, we performed the simulation of radar signals with two models. GprMax, a software based on the finite-difference time-domain (FDTD) method [22] to model ground-penetrating radar (GPR) responses from arbitrarily complex targets was used in this study. Two-layer models were established, as illustrated by Figure 7. The first layer represents the model of AA4b, AA4a, Bl3, and AA3c layers. The dielectric constant of AA3c is in the range of 2.0 to 2.2 [9] and in this paper the average of these layers 1 = 2.2 is used. The second layer models Bl2, nearly pure water ice, 2 = 3.15 [8,9]. In the smooth surface model, the thickness of the first layer is 325 m, and the thickness of the second layer of water ice is 50 m. According to the observation results of SHARAD, the thickness difference between winter and summer is in the 0-45 m range, and the specific description is in 3.1.1. Therefore, the elevation of the rough surface is lower than that of the smooth model surface. The simulation results of a smooth surface ( Figure 1b) and a rough surface (Figure 1c) are compared to show the influence of surface roughness on radar signals.
The surface roughness is calculated with Equation (1). This parameter can reflect the elevation change of each vertex of the surface cell grid in the analysis area. The value from elevation data is highly correlated with the clutter of the SHARAD profile [23]. Roughness = Root mean square elevation/Average elevation (1) The elevation data in the study region were generated with CTX image pairs and the NASA Ames Stereo Pipeline (ASP) tool [24], and the resolution is 18 m/grid point [25]. Although a DEM model derived from HiRISE data is more accurate and trusted in estimating the influence of seasonal variations of surface roughness from SHARAD radar signal [26], the best available DEM model in the study area we obtained was developed from CTX images, which also shows clear changes of geomorphology in different seasons.

Surface Roughness
To analyze the influence of surface roughness on uncertainties of measurements, we performed the simulation of radar signals with two models. GprMax, a software based on the finite-difference time-domain (FDTD) method [22] to model ground-penetrating radar (GPR) responses from arbitrarily complex targets was used in this study. Two-layer models were established, as illustrated by Figure 7. The first layer represents the model of AA 4b , AA 4a , Bl 3 , and AA 3c layers. The dielectric constant of AA 3c is in the range of 2.0 to 2.2 [9] and in this paper the average of these layers ε 1 = 2.2 is used. The second layer models Bl 2 , nearly pure water ice, ε 2 = 3.15 [8,9]. In the smooth surface model, the thickness of the first layer is 325 m, and the thickness of the second layer of water ice is 50 m. According to the observation results of SHARAD, the thickness difference between winter and summer is in the 0-45 m range, and the specific description is in 3.1.1. Therefore, the elevation of the rough surface is lower than that of the smooth model surface. The simulation results of a smooth surface ( Figure 1b) and a rough surface (Figure 1c) are compared to show the influence of surface roughness on radar signals.   Table 1. The parameters of the second set of subsatellite points shown in (e-f) are given in Table 2. Figure (a,e) are the  Table 1. The parameters of the second set of subsatellite points shown in (e,f) are given in Table 2. Figure   The surface roughness is calculated with Equation (1). This parameter can reflect the elevation change of each vertex of the surface cell grid in the analysis area. The value from elevation data is highly correlated with the clutter of the SHARAD profile [23]. Roughness = Root mean square elevation/Average elevation (1) The elevation data in the study region were generated with CTX image pairs and the NASA Ames Stereo Pipeline (ASP) tool [24], and the resolution is 18 m/grid point [25]. Although a DEM model derived from HiRISE data is more accurate and trusted in estimating the influence of seasonal variations of surface roughness from SHARAD radar signal [26], the best available DEM model in the study area we obtained was developed from CTX images, which also shows clear changes of geomorphology in different seasons.

Case Studies of Seasonal Variations
Previous studies show that seasonal variations and yearly erosions occur on the AA 4 layer, and the CO 2 -ice layer of AA 3c might sublimate through the cracks of AA 4a in summer and lead to a collapse in some areas [9]. To obtain the temporal variation of the same layer, we select one point pair with a subsatellite distance of less than 3 km. The radar profile of each point provides the positions of the surface and lower interface of AA 3c , and then the thickness can be obtained. The seasonal variation is estimated by comparing the thickness at different seasons, e.g., thickness in summer and winter. Figure 7 presents two comparison cases of radar profiles in different seasons.
The frosts on the surface sublimate in summer and are deposited in winter at the south pole. The surface roughness is usually higher in summer than in winter. In addition, AA 4b is eroded at a constant rate by expanding outwards from the steeply sided walls of the depressions, and the AA 3c unit might sublimate through the cracked BL 3 layer in the spring and summer. In the radargram from Figure 7, the surface radar reflections appear wider, and clutter is more abundant in summer than in winter, causing uncertainties in measuring the thickness. Another prominent phenomenon caused by seasonal change is that the subsurface echo also becomes wider and fuzzier in the summer, leading to measurement errors. The "foggy" signals appear at depth in the summer radargram of Figure 7a, resulting in difficulties in identifying other AA 3 subunits. The causes of "foggy" echoes, whether they are formed due to surface roughness or are also related to subsurface roughness, e.g., sublimation of the AA 3 layer, are analyzed with the simulation results (details in Section 3.2.2). In Figure 7d, although echoes of the surface and the upper interface of AA 3c are sharper than those in the summer season, the lower interface of AA 3b in the winter disappears in both the radargram and radar single trace profile. Table 3 lists six data pairs of SHARAD for observations of seasonal thickness variations in S1 and S2. Table 1 shows that the thickness of the AA 3 layer does have seasonal variations ranging from~16-45 m, which could reach 3 times the vertical resolution of SHARAD data. The results given in Tables 1 and 2 indicate that the thickness of shallow deposits in winter is higher than that in summer, and the FWHM of surface layers in summer increases by~20% compared to winter.

Statistical Study of Seasonal Variations
Through the SHARAD observations of 11 terrestrial years from May 2007 to September 2018, we can see some variations in the thickness of AA 3c during different seasons based on the statistical analysis of these track intersections in S1, S2, and S3.
From M28 in summer to M31 or M32 in winter, most AA 3 layers have a positive change in thickness,~16 m to 45 m in S1, as shown in Figure 8a,b. In the studied S1 region, the majority of the thickness variation is larger than the SHARAD vertical resolution, and a few of them reach three times the resolution. As we expected, the pulse width of the surface echo in winter is narrower than that in summer, meaning their differences are all below 0. We speculate that a layer variation larger than the MOLA observation might be caused by surface roughness. However, no clear correlation between the width variation and thickness variation is observed from Figure 9, which means that the large change of surface roughness does not correspond to a large thickness variation. Figure 8c shows the average thickness of shallow deposits in each season in S1, S2, and S3. We found that the thickness of each region in the winter is greater than that in the summer, which is consistent with previous case studies. The range of the average thickness variation is from~7 m to 20 m, higher than prior studies [15] which reported~1 m of seasonal variation on the surface AA 4 layer. The error is calculated from the standard deviation, whose range is~9-21 m. The error in winter is relatively smaller than in others, and may be due to the stable data quality on the smooth surface in winter. The source of the measured thickness is difficult to model as it represents a comprehensive change of multiple deposition layers: AA 4b , Bl 3 , or AA 3c . In addition, the uncertainties of the measurements, especially the influences of surface and subsurface roughness, should be considered to obtain the quantity of the temporal variations of the shallow deposits in the south polar region.  tical axis on the right shows the width difference about FWHM, or the width of winter subtracted from that in summer. M28summer: Martian year 28 in summer; M31winter: Martian year 31 in winter; M32winter: Martian year 32 in winter. (c) the average thickness of AA3c in different seasons. S1, S2 and S3 regions are indicated by the colors blue, red, and yellow, respectively. The horizontal axis is the solar longitude, and the vertical axis is the depth of AA3c. Please note that data points with FWHM of the surface and subsurface echoes >150 ns are not included. For example, the vacancy of the summer data in S3, or the solar longitude is 270-360 degrees, is because no data meet the standard for the season.

Seasonal Variations of Surface Roughness
The temporal variations of the thickness measured by SHARAD could be affected by the seasonal changes of surface roughness. Table 4 lists the maximum value of the surface roughness obtained with Equation (1) and elevation data using CTX image pairs. Figure  9 shows the average maximum surface roughness and corresponding season in S1, S2, and S3. Due to the lack of winter/fall CTX stereo images in the study region, only the roughness data in the summer/spring seasons were analyzed. Based on Table 4 and Figure 9, we can conclude that the surface roughness in summer is greater than that in the spring, especially in S1, where the difference could reach 5%. This is because the surface dry-ice layer (AA4b) begins to sublimate in spring, and then later in the summer, more surface CO2 sublimates to the atmosphere and the rough surface is exposed.
To find the relationship between the surface roughness and the radar signal, we performed a statistical analysis on the 3 dB pulse width of surface echoes. Figure 10 shows the statistic study of all the data including those with FWHM greater than 150 ns. For the data points with FWHM larger than 150 ns, the numbers of autumn points and winter points account for more than 60% of the total data, while for the data points with FWHM smaller than 150 ns, data of spring and summer seasons have a larger proportion. This seems due to the fact that the increased width of the radar signal is caused by the rough surface in the summer. The pulse width of the subsurface shows a similar pattern, with

Seasonal Variations of Surface Roughness
The temporal variations of the thickness measured by SHARAD could be affected by the seasonal changes of surface roughness. Table 4 lists the maximum value of the surface roughness obtained with Equation (1) and elevation data using CTX image pairs. Figure 9 shows the average maximum surface roughness and corresponding season in S1, S2, and S3. Due to the lack of winter/fall CTX stereo images in the study region, only the roughness data in the summer/spring seasons were analyzed. Based on Table 4 and Figure 9, we can conclude that the surface roughness in summer is greater than that in the spring, especially in S1, where the difference could reach 5%. This is because the surface dry-ice layer (AA 4b ) begins to sublimate in spring, and then later in the summer, more surface CO 2 sublimates to the atmosphere and the rough surface is exposed.
To find the relationship between the surface roughness and the radar signal, we performed a statistical analysis on the 3 dB pulse width of surface echoes. Figure 10 shows the statistic study of all the data including those with FWHM greater than 150 ns. For the data points with FWHM larger than 150 ns, the numbers of autumn points and winter points account for more than 60% of the total data, while for the data points with FWHM smaller than 150 ns, data of spring and summer seasons have a larger proportion. This seems due to the fact that the increased width of the radar signal is caused by the rough surface in the summer. The pulse width of the subsurface shows a similar pattern, with the exception of S1. From these results, the width of the surface/subsurface echo of the radar has seasonal variations and is related to the seasonal change of the surface roughness. Remote Sens. 2022, 14, x FOR PEER REVIEW 15 of 23 Figure 10. The proportions of echoes with FWHM >150 ns (orange color) and ≤150 ns (grey color) in different seasons in S1-S3. Data include the surface, subsurface, and all the echoes in the study area. The horizontal axis represents seasons, and the vertical axis shows proportional distribution within each season.   Figure 11 shows the thickness variations of shallow deposits at all subsatellite points with the observation period in the three study areas. The thickness values span multiple levels with the fixed gap determined by radar data-sampling rate. Although the thickness of AA 3c has periodic variations, it did not show an obvious increasing or decreasing trend during these 11 terrestrial years from the moving average. The reason could be that the slight yearly change of the thickness is not detectable by SHARAD.  Figure 11 shows the thickness variations of shallow deposits at all subsatellite points with the observation period in the three study areas. The thickness values span multiple levels with the fixed gap determined by radar data-sampling rate. Although the thickness of AA3c has periodic variations, it did not show an obvious increasing or decreasing trend during these 11 terrestrial years from the moving average. The reason could be that the slight yearly change of the thickness is not detectable by SHARAD.

Variations of CO2 Layer Thickness in 11 Terrestrial Years
Another observation from Figure 11 is that in the same area, the thickness of AA3c has been constantly fluctuating in value. It raises a reasonable suspicion that thickness variation is mainly due to regional differences of two compared data points. To find out the possible cause, we show the distribution of the thickness value and the corresponding average subsatellite distance between the data pair in Figure 12. The statistical results display a normal distribution of temporal variations of CO2 ice in S1, S2, and S3. We expect that the average distance increases with the thickness variations, as it is likely the two points further apart have large variations in the AA3c layer thickness. However, such correlation between thickness variations and the average distance is not very clear in Figure  12. The systematic variations of thickness may not be the main reason and we could not exclude the possibilities that observed variations are caused by measurement errors such as low radar resolution and the influence of surface roughness. To evaluate the latter, a simulation of the radar signal was performed, and the results are given in the following subsection.  Another observation from Figure 11 is that in the same area, the thickness of AA 3c has been constantly fluctuating in value. It raises a reasonable suspicion that thickness variation is mainly due to regional differences of two compared data points. To find out the possible cause, we show the distribution of the thickness value and the corresponding average subsatellite distance between the data pair in Figure 12. The statistical results display a normal distribution of temporal variations of CO 2 ice in S1, S2, and S3. We expect that the average distance increases with the thickness variations, as it is likely the two points further apart have large variations in the AA 3c layer thickness. However, such correlation between thickness variations and the average distance is not very clear in Figure 12. The systematic variations of thickness may not be the main reason and we could not exclude the possibilities that observed variations are caused by measurement errors such as low radar resolution and the influence of surface roughness. To evaluate the latter, a simulation of the radar signal was performed, and the results are given in the following subsection.
Remote Sens. 2022, 14, x FOR PEER REVIEW 17 of 23 Figure 12. Statistical result of the variations of AA3c thickness and the corresponding subsatellite distance between the compared node pair in study areas. The horizontal axis presents the thickness change of AA3c of each pair of subsatellite points, whose distance is less than 3 km. The vertical axis on the left shows the percentage. The vertical axis on the right shows the average distance for each difference.

Simulation Results of Rough Surface
From the results of the surface roughness simulation on the radar signal, the surface echo of the smooth surface (Figure 13a) is clear and sharp, and the position of the surface echo can be accurately determined. When the surface roughness increases (Figure 13b), the width of the surface echo also increases significantly, and the width can reach ~30 m at the widest point. We believe that the increase in surface roughness results in a rich surface echo signal and the appearance of signal superimposition, which causes errors regarding the position of the surface echo, thereby affecting the accuracy of measurement on the thickness of AA3c. The simulation results show that the surface roughness does affect the measurement of layer thickness by increasing the width of surface and subsurface echoes. In addition, not only is the surface echo affected by the roughness, but the subsurface interface in the radargram becomes wider and rougher, which contributes further to the uncertainties of the measurement. It also indicates that the blurred subsurface boundary in the summer may not be caused by the rough subsurface but by the sublimation of the surface CO2 layer. Figure 12. Statistical result of the variations of AA 3c thickness and the corresponding subsatellite distance between the compared node pair in study areas. The horizontal axis presents the thickness change of AA 3c of each pair of subsatellite points, whose distance is less than 3 km. The vertical axis on the left shows the percentage. The vertical axis on the right shows the average distance for each difference.

Simulation Results of Rough Surface
From the results of the surface roughness simulation on the radar signal, the surface echo of the smooth surface (Figure 13a) is clear and sharp, and the position of the surface echo can be accurately determined. When the surface roughness increases (Figure 13b), the width of the surface echo also increases significantly, and the width can reach~30 m at the widest point. We believe that the increase in surface roughness results in a rich surface echo signal and the appearance of signal superimposition, which causes errors regarding the position of the surface echo, thereby affecting the accuracy of measurement on the thickness of AA 3c . The simulation results show that the surface roughness does affect the measurement of layer thickness by increasing the width of surface and subsurface echoes. In addition, not only is the surface echo affected by the roughness, but the subsurface interface in the radargram becomes wider and rougher, which contributes further to the uncertainties of the measurement. It also indicates that the blurred subsurface boundary in the summer may not be caused by the rough subsurface but by the sublimation of the surface CO 2 layer.  Figure 13. Comparisons of simulated radargram on a smooth (a) and a rough (b) surface. The horizontal axis represents the points number, and we set the distance between every two points as 10 m. The vertical axis shows the depth in meters.

Discussion
With 11 terrestrial years of observation data, we found that the thickness of the AA3c layer (which also includes surface AA4 layer if present) is essentially between 300 m to 400 m, and the layer of the thin water ice between AA3c and AA3b, or Bl2, is mainly 30 m to 60 m thick. There are two-or even three-layers of CO2 ice in some regions, and each CO2ice layer is separated by a thin layer of water ice. This observation was consistent with previous research results. The diagram of the stratigraphic structure of AA3 is shown in Figure 1, which is drawn based on the SHARAD radargram (Figure 7). The underground layered sediment measured in this work includes AA4b, AA4a, Bl3, and AA3c layers.
Through measuring the thickness variations of AA3c, the results in Section 3.1 show that there is indeed a seasonal thickness difference, meaning that the thickness of the AA3c layer in the same location may vary in different Martian seasons. From the results, we found that the thickness of the shallow deposits thickens in winter; consequently, the dryice layers decrease in the summer. The maximum variation value is about three times of SHARAD resolution; hence, it is unlikely that the observed seasonal variations are random noises due to detection limitations. In addition, this phenomenon is not caused by regional factors as shown in Figure 12.

Discussion
With 11 terrestrial years of observation data, we found that the thickness of the AA 3c layer (which also includes surface AA 4 layer if present) is essentially between 300 m to 400 m, and the layer of the thin water ice between AA 3c and AA 3b , or Bl 2 , is mainly 30 m to 60 m thick. There are two-or even three-layers of CO 2 ice in some regions, and each CO 2 -ice layer is separated by a thin layer of water ice. This observation was consistent with previous research results. The diagram of the stratigraphic structure of AA 3 is shown in Figure 1, which is drawn based on the SHARAD radargram (Figure 7). The underground layered sediment measured in this work includes AA 4b , AA 4a , Bl 3 , and AA 3c layers.
Through measuring the thickness variations of AA 3c , the results in Section 3.1 show that there is indeed a seasonal thickness difference, meaning that the thickness of the AA 3c layer in the same location may vary in different Martian seasons. From the results, we found that the thickness of the shallow deposits thickens in winter; consequently, the dryice layers decrease in the summer. The maximum variation value is about three times of SHARAD resolution; hence, it is unlikely that the observed seasonal variations are random noises due to detection limitations. In addition, this phenomenon is not caused by regional factors as shown in Figure 12.
The previous study shows that seasonal sublimation and deposition of dry ice occurs on the surface near the Martian southern polar ice sheet. However, the change in the surface dry ice measured with MOLA data is within 10 m, and most is concentrated in 1 m [15], much lower than the upper bound of results given by SHARAD observations. The internal dry-ice layer can be sealed by a~15-45 m thick water-ice layer [8] and may sublimate in the spring or summer if there are cracks on the surface [8,9]. Unfortunately, we did not observe any cracks on AA 4 in our study area at present. A conceptual model was proposed that the subsurface dry ice sublimates in the event of no gaps being on the surface in the northern polar cap. The sunlight shines through the water ice into the lower layer of dry ice, and the dry ice sublimates until the pressure is large enough to break the thin upper layer of water ice [27]. According to this model, if the water ice of the upper layer is brittle, then the AA 3c can escape. In late spring, the dry-ice layer of the seasonal ice cap was evaporated, and then the lower layer of water ice or Bl 3 will be exposed. When the summer comes, sunlight will penetrate the water-ice layer, and it is possible that the AA 3c sublimates in the same way to form cracks on the overlying layer. Eventually, the underlying dry ice will reduce, and the surface water-ice layer will collapse. However, even if the internal dry ice in the study area escapes, can the total thickness variation of the shallow deposits reach tens of meters?
To investigate whether the SHARAD observation data reflect the real volume change or are caused by uncertainties of measurement, we conducted roughness analysis on the elevation data in different seasons and found that roughness in summer is greater than that in spring. At the same time, we observed that the FWHM is wider in the summer than that in other seasons. From the consistency of the variation for signal width on the surface and thickness of dry ice, we suspect that the surface roughness will lead to the broadening of the width of radar signals and an increase in errors of measurement. Specifically, it is the surface roughness which affects the width of the radar's echo signal, which in turn affects the position of the echo, and finally results in a difference in layer thickness. The simulation results also demonstrate that the surface roughness can also affect the width of subsurface echoes, which further contributes to the deviation of the measurements. Therefore, the resolution of SHARAD is not high enough to quantify the vertical changes of the thickness of shallow deposits within the geologically short period.
The variation in mass for the dry ice buried in the south pole of Mars is closely related to the history of the Martian climate's evolution. The intersection of every two tracks was used to extract the subsatellite points of different years in the same area for our comparative analysis. In the 11 terrestrial years from 2007 to 2018, the thickness of the shallow dryice layer, or AA 3c , has had no obvious trend of thickening or thinning. It is currently believed that although the layer of dry ice may experience erosion, the variations during a geologically short period are not detectable by SHARAD. In the future, an instrument with high vertical resolution can provide more data on the elevation changes of SPLD for the modeling of the evolution of the Martian atmosphere.

Conclusions
The sediments in the polar regions of Mars carry a record of its past climate. Imagery and spectrum data can be used to assess annual changes in the lateral coverage of volatile deposits. However, vertical variations, in particular the changes of the buried CO 2 -ice layer, cannot be assessed this way.
In this work, we used data from the SHARAD, which has been measuring the Martian stratigraphy since 2007, to obtain temporal variations of the shallow deposits in the south polar region. Analyzing the data produced through different seasons during seven successive Martian years in the region of interest, we have observed strong seasonal variations showing that shallow deposits are thicker in the winter than in the summer. However, the maximum variation values are much higher than found in previous estimations based on MOLA data. The reason could be the measurement uncertainties introduced by variations in the surface roughness, as this work shows that there is a correlation between the seasonal changes of surface roughness, calculated with the elevation data generated with CTX image pairs, and the width of the surface radar echoes. In addition, the simulation results suggest that the subsurface radar echoes are also affected by the surface roughness, leading to further measurement deviations. SHARAD was thus successful in detecting clear seasonal variation of the southern polar deposits, but its limited resolution and influences of surface roughness make it difficult to derive quantitative vertical variations within a short geologic time. In the 11 terrestrial years from 2007 to 2018, the thickness of the shallow dry-ice layer, or AA 3c , shows no evident trend in thickening or thinning.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14030435/s1. Figure S1: This group of pictures shows the thickness of AA3c and the echo's width of the subsatellite point before and after filtering data in S1 ( Figure S1(a-d)) and S3 ( Figure S1(e-h)). Figure S2: The white box shows the subsurface reflectors in SPLD, the track number is 0564201. Table S1: This is the supplementary material of Table 4. The index corresponds to that in Table 4. The third column CTX_number indicates the source of each DEM. Ls is the solar longitude corresponding to each CTX.