Lower Bounds on the Thickness and Dust Content of Layers within the North Polar Layered Deposits of Mars from Radar Forward Modeling

The Mars Reconnaissance Orbiter's Shallow Radar (SHARAD) emits radar signals and records their reflections from layer boundaries within the Martian north polar ice cap. Previous studies have suggested that the ice cap is composed of thin dust-rich layers between thicker layers of nearly pure water ice. The prevailing hypotheses suggest that each dust-rich layer represents either a period of ice sublimation at the poles or a period of reduced ice deposition relative to dust deposition. To test whether thin dust beds are a plausible hypothesis for the observed SHARAD reflectors, we use RadSPy (radar sounding simulator in Python, https://github.com/scourvil/RadSPy.git), an open-source N-layer radar sounder forward-modeling software that we have developed and describe herein. We forward model radar data from thin dust-rich beds interspersing pure ice, and compare them to observed radar reflection data over Gemina Lingula in the north polar layered deposits (NPLD). We consider two end-member cases: (1) thin beds composed entirely of dust, but with thickness varying from 0.05 m to 0.4 m; and (2) dust beds all with the same thickness, but with varying dust content. We find that the observed reflections can be explained by either scenario, i.e., varying thickness or varying dust content, and we conclude that a combination of both is likely. More importantly, our results provide lower bounds on the layer thickness and dust fraction for the flat-lying reflectors of Gemina Lingula in the NPLD. Our findings support the thin dust layer hypothesis, providing new constraints on layer composition and geometry for Mars climate researchers.


Introduction
Similar to how tree rings track wet and dry seasons, Mars' north polar layered deposits (NPLD) and the pattern of dust and ice within them may track periods of warm and cool climates. However, the relationship between the layers and the Martian climate is unclear. The NPLD are kilometer-thick structures of nearly pure water ice interspersed with either impure water ice or dust layers (Howard et al. 1982;Nunes & Phillips 2006;Phillips et al. 2008;Grima et al. 2009;Putzig et al. 2009;Lalich et al. 2019). As with large-scale ice sheets on Earth, these polar layered deposits are thought to contain millions of years of climatological data, the deciphering and understanding of which is paramount to Mars polar science (see Smith et al. 2020 and references therein). The deposition and ablation of the NPLD is likely a function of the orbital obliquity of Mars, with periods of low obliquity associated with large-scale deposition and periods of high obliquity associated with ablation (e.g., Levrard et al. 2007). However, ablation and deposition could be taking place concurrently in different areas of the polar layered deposits, depending on the local climate and other complicating factors (e.g., Perron & Huybers 2009). Furthermore, the observed layering in the NPLD may be more dependent on dust availability and temperature gradients than on the climate of the polar region as a whole (e.g., Newman et al. 2005;Toigo et al. 2020). Regardless, to uncover the relation between the NPLD and the past climate of Mars, one needs estimates of the NPLD composition.
Orbital sounding radars provide the ability to explore the subsurface of these features in great detail since water ice is nearly transparent to the part of the electromagnetic spectrum utilized by the sounders. In fact, a large body of science has amassed concerning the interpretation of radar data at the Martian poles (e.g., Putzig et al. 2018 and references therein). As part of that effort, attempts have been made to correlate the layers observed in High Resolution Imaging Science Experiment and Context Camera imagery with the layers observed in the radar data with limited success (Christian et al. 2013). A better correspondence may be obtained by comparing radar layers to protrusion profiles obtained from topographic data (Becerra et al. 2019(Becerra et al. , 2020. However, visible layers may not correlate with layering observed in radar data because radar reflections originate from stark contrasts in dielectric properties, while visible layering could be present due to a number of factors independent of dielectric properties. Therefore, correctly interpreting how radar reflections relate to physical layering processes is challenging. Two prevailing hypotheses exist for the formation of the observed radar layers. The first suggests accumulation of ice and dust occurs simultaneously, however the relative rate of ice and dust accumulation varies such that some periods are dominated by dust deposition resulting in less icy layers that produce radar reflections (Putzig et al. 2009;Hvidberg et al. 2012). The second supposes that occasional periods of ice instability at the poles forces ice to sublimate/ablate leaving behind a lag deposit composed of dust that reflects radar signals (Levrard et al. 2007). In reality, the layers may be constructed through a combination of these two scenarios. However, if the first hypothesis is dominant, then the variation of observed radar reflection power should be driven primarily by variations in the dielectric properties caused by varying dust to ice ratio in the layers. In such a case, the thicknesses of the layers may range between >1 m and 8 m, which is below the 8 m resolution of SHARAD in ice (Seu et al. 2004). Alternatively, if the second hypothesis is dominant, we should expect the reflections to be created from the thinnest possible layers of dust lag because modeling suggests that a thin layer of dust is enough to prevent further sublimation of ice (Skorov et al. 2001). In this case, the layers would all have approximately the same dielectric permittivity, but have variable thicknesses of approximately <1 m which control the reflection power. However, as Lalich et al. (2019) have shown by modeling radar reflections from dust layers thinner than the resolution of SHARAD that are bound above and below by water ice, one cannot distinguish one hypothesis from the other using solely radar sounding data because the reflectivity can be controlled by either the thickness or the permittivity of the dust layer.
One way to combat nonuniqueness of radar sounding data is to employ numerical models to test various geological interpretations. While a variety of simulators exist for orbital radar sounders in the scientific literature (e.g., Leuschen et al. 2003;Nunes & Phillips 2006;Xu et al. 2006;Mastrogiuseppe et al. 2016;Gerekos et al. 2018;Thakur & Bruzzone 2019), most are instrument specific and are not publicly available, which greatly restricts their use and reduces the ability of others to cross-check the quality of radar data analysis and interpretation. To advance the understanding of the formation mechanisms and composition of polar layered deposits and to allow testing of more complicated subsurface hypotheses, we developed RadSPy, an open-source 1D electromagnetic wave forward-modeling software in Python 3.6 and based on the reflectivity method described in Wait (1970). The simulator is applicable to any user-defined 1D N-layer model, where each layer is characterized by a relative dielectric permittivity, loss tangent or conductivity value, and thickness. Additionally, the simulator is compatible with any user input radar source (e.g., pulse, chirp, continuous wave, etc.) making it applicable to any orbital or ground-based sounder. The simulator assumes that (1) the incident radar source is a plane wave, (2) the dielectric permittivity and loss tangent are constant with respect to frequency (i.e., no dispersion), and (3) the surface and subsurface interfaces are completely flat-lying (i.e., no roughness). The program is freely available under the Apache 2.0 license on GitHub: https://github.com/scourvil/RadSPy.git. Version 1 has been preserved on Zenodo (Courville & Perry 2021).
Using the simulator, we test the two aforementioned polar layered deposits formation scenarios over the Gemina Lingula region of the NPLD. Since we cannot directly distinguish one hypothesis from the other, we generate plausible models for both hypotheses. By assuming that the reflectivity is driven by variation in dielectric properties, we can investigate the minimum possible dielectric permittivity of the layers. By assuming the reflectivity is driven by thickness variations, we can investigate the minimum thickness required for each layer. In Section 2, we describe the mathematical theory we use to simulate radar sounding data. In Section 3, we describe our procedure for generating plausible physical models for the observed radar data using our radar simulator. In Section 4, we illustrate two model results, one representing the minimum permittivity possible, and the other representing the minimum thickness possible. Finally, in Section 5, we discuss how our results can be used to help future researchers understand the past climate of Mars.

1D Electromagnetic Wave Propagation Theory
Radar sounders provide a means to probe the subsurface of planetary bodies. Whereas visible light cannot typically pass through the surface of planetary bodies, radio waves (3 kHz-300 GHz) can transmit up to several kilometers into the subsurface, depending on the properties of the transmitted radio waves and of the ground. The speed that the electromagnetic wave travels through the subsurface is dependent on the real part of the dielectric permittivity, ¢  (also known as the dielectric constant), and the decay rate, or attenuation, is dependent on the imaginary part of the dielectric permittivity, ò″ (Campbell 2002). In orbital radar sounding, attenuation is typically expressed in terms of the loss tangent, d tan , or the medium's conductivity, σ, rather than ò″ (Campbell 2002). The three quantities are related by where ω is the angular frequency of the radar wave. In this article, we express attenuation in terms of the mediumʼs conductivity, such that all further ò values refer strictly to the real part of the dielectric permittivity. Additionally, we refer to the relative dielectric permittivity, ò r , such that the real dielectric permittivity is given by ò = ò r ò 0 , where ò 0 = 8.85 × 10 −12 F m −1 . In a homogeneous medium with a given dielectric permittivity, a radar plane wave propagates downward until it attenuates away. However, if the wave encounters a layer with a differing dielectric permittivity and/or conductivity, a portion of the wave's energy is reflected upward while the remaining energy is transmitted downward (Campbell 2002). Thus, if a radar wave from a sounder travels through the ground and encounters subsurface geologic layering with a significant dielectric contrast, the portion of the radar energy that is reflected is subsequently recorded by the instrument.

Electromagnetic Wave Reflection from Horizontal Layers
A monochromatic electromagnetic plane wave with angular frequency ω that is normally incident upon a layered medium has an electric field that obeys a 1D wave equation, where E is the electric field amplitude in a direction orthogonal to the propagation direction z, σ is the conductivity of the medium, μ is the magnetic susceptibility, and ò is the real dielectric permittivity of the medium. Figure 1 illustrates the 1D N-layer reflectance problem. We assume that the solution to the electric field in Equation (3) within each layer n is given by the summation of a down-going wave and an up-going wave, where a n is the amplitude of the down-going wave, and b n is the amplitude of the up-going wave. Since the source is known, the amplitude of the down-going wave in the first layer, a 0 , is also known. Additionally, since there is no interface beyond the final layer N, there is no up-going wave in the Nth layer, thus b N = 0. For a monostatic radar system (i.e., transmitter and receiver are colocated) we wish to solve for the electric field that returns back to the source point after reflecting from the layered medium. This is determined by solving for b 0 . Each layer interface contributes to the amplitude b 0 of the reflected up-going wave. Thus, we must solve for every a and b coefficient. Since there are N layers, there are N total a coefficients and N total b coefficients. However, since a 0 and b N are known, this leaves a total of 2(N − 1) unknown coefficients.
For an N-layer model, there are N − 1 interfaces. Since the electromagnetic wave must be continuous, at each interface between a layer n − 1 and a layer n, we know that the electromagnetic wave must satisfy two boundary conditions, where B n and E n are the magnetic and electric fields in the nth layer, respectively. Given that n n the first boundary condition then becomes Since there are N − 1 interfaces, there are 2(N − 1) total boundary conditions. Substituting Equation (4) into each boundary condition then results in a system of 2(N − 1) equations as a linear function of the 2(N − 1) unknown a and b coefficients. Thus, we are able to solve for the coefficient b 0 , the amplitude of the reflected wave at the location of the source. The solution to this system of equations, as derived by Wait (1970), is In this solution, h n is the thickness of the nth layer, and the term γ n is given by

Full Radar Pulse Simulation
The previous section displays the solution for the N-layer reflectance problem for a monochromatic wave. To solve for the reflectance from a radar pulse spread over a range of frequencies, we repeat the calculation in Section 2.1 for all frequencies in the pulse. Approximating the pulse as an array of discrete frequency components, the reflection for an angular frequency component, ω, is is the electric field amplitude of the reflected pulse for the angular frequency component ω. This means that the total reflected pulse in the time domain is given by the inverse Fourier transform of the coefficients b 0 . We assume that the layer's conductivity and permittivity are constant over all frequencies. We also assume that the material has no magnetic susceptibility such that all the material has a magnetic permeability of μ = 4π10 −7 H m −1 . These are acceptable approximations for most planetary surfaces (Campbell 2002).

Source Pulses and Pulse Compression
It is common for orbital radars to use a chirped pulse, also known as a linear frequency modulated (FM) pulse. From Cumming & Wong (2005), a standard linear FM pulse when discretized and expressed in the frequency domain is given by where f 0 is the center frequency of the pulse, Δf is the bandwidth of the pulse, and T is the pulse length in time. The values of a 0 [ω] are the downward-going source pulse amplitude coefficients to be put into Equation (14). To provide high resolution, reflections of chirped pulses are compressed using convolution of a matched filter (Cumming & Wong 2005). The matched filter is the complex conjugate of the original pulse. Therefore, the final signal in the frequency domain after pulse compression is given by 1 6 out 0 0 * In accordance with standard SHARAD processing developed by the US SHARAD Science team (Campbell et al. 2011), we apply a Hann window to b 0 [ω] before compression. For the purpose of this NPLD study, we use the ideal SHARAD source pulse and matching filter. The ideal SHARAD source pulse is a linear FM downswept pulse from 25 MHz to 15 MHz that has a duration of T = 85.05 μs (Seu et al. 2007). So for this study, f 0 = 20 MHz and Δf = 10 MHz.

Model Generation and Data Matching Methodology
We create models of the aforementioned two NPLD composition hypotheses that match radar reflection data from the NPLD using the forward-modeling theory described in Section 2. To do this, we choose SHARAD reflection data from a nearly flat location (83°. 942N, 36°.947E; Figure 2) within Gemina Lingula, which obviates the need to consider sloping reflectors or off-nadir returns. In order to use our forwardmodeling theory, we assume that the surface and subsurface interfaces are smooth such that scattering from surface roughness can be ignored. If this assumption is incorrect, then our models may overestimate the reflected power from each layer and subsequently cause us to underestimate the layer thickness and/or permittivity values. However, this is acceptable as long as we interpret our results as lower bounds.
We identify 20 SHARAD observations that intersect near the location of interest, and process them through the Colorado SHARAD Processing System (CO-SHARPS, see Putzig et al. 2016) using default processing parameters for the US SHARAD Science Team processor (US Focused Processor B, see Campbell et al. 2011). For each radargram (e.g., Figure 3), we interpret reflection events based on their continuity across Gemina Lingula and their consistent power above the SHARAD noise floor (≈−35 dB), while avoiding likely side lobes (pulse-compression artifacts). From each radargram, we extract six frames centered on the investigation location ( Figure 3; left panel, darkened-gray region), and average them together into one frame in order to improve the signal-to-noise ratio. Averaging over more than six frames may distort the shape of reflection events because reflection events in separate frames may not be aligned in time the further apart the frames are in space due to gently sloping terrain. Our goal is to generate layering models that can be used to recreate the data in the averaged frame. Because each radargram has different ionospheric interference, is not centered exactly on the same location, and has different surface reflectivity, the shape and location of reflection events are not perfectly consistent across every radargram (Foss et al. 2017). Therefore, we choose to generate models for each radargram separately rather than average all radargrams together. To test the first hypothesis that the reflections are made by dust layers with variable dusty content surrounded by pure ice, we generate models with dusty layers of constant thickness but varying permittivity. Because layer thickness also controls the reflected power, we choose the thickness that would generate the strongest possible reflection for a layer with a given permittivity, thereby ensuring that our layer permittivity results for this test are low bounds. Specifically, we use a brute-force approach to find the model that visibly fits the data best under the following constraints: 1. Each dusty layer has a thickness of 1.7 m, which produces the maximum possible reflection power for a given permittivity (Lalich et al. 2019). 2. The relative permittivity value of each layer can be between 3.1 and 8.8 at intervals of 0.05. 3. There is one dusty layer for each reflection event in the data, and the best location of each dusty layer is checked on intervals of 1 m.
Under these constraints, there is one clear best fitting model. Additionally, assuming that the dusty layers are ice with rock inclusions, where pure water ice has a relative dielectric permittivity of 3.1 and the rock has a relative dielectric permittivity of 8.8, we can convert the model dielectric permittivity results into a dust volume fraction using the formula, where ò eff is the dielectric permittivity of the effective medium in the model result, ò m is the dielectric permittivity of the surrounding ice medium, and ò i is the dielectric permittivity of the rock inclusions (Choy 1999).
To test the hypothesis that the reflections are caused by thin dust lags, we generate models consisting of layers with constant permittivity, but variable thickness. Again, we use a brute-force approach to find the model that visibly fits the data best, but under a different set of constraints: 1. Each dust lag layer has a relative permittivity value of 8.8, which corresponds to no ice (Lalich et al. 2019). 2. The thickness of each lag layer can be any value between 0 m and 1.7 m on intervals of 0.005 m. 3. There is one dust lag layer for each reflection event in the data, and the best location of each dust layer is checked on intervals of 1 m.
We choose the range between 0 m and 1.7 m thickness because in this range, reflection strength strictly increases with increasing thickness (Lalich et al. 2019), and therefore our model thickness results represents a lower bound. Under these constraints, there is one clear best fitting model. For all models, we assume the material between the dust layers is pure water ice with a relative permittivity of 3.1 and a loss tangent of 0.001, which we base on Grima et al. (2009). However, we use a lower value of loss tangent than their reported 0.0026 because that value is an absolute high end that does not consider losses due to reflections or losses from attenuation within the dust layers. Our simulations assume that all dust layers have a significantly higher loss tangent of 0.01, which matches values observed by Campbell & Morgan (2018). Lastly, our models assume the first layer is 300 km of free space to represent the average orbital altitude of MRO. By repeating this exercise for each radargram over the same location, we estimate the uncertainty of our derived model results by displaying the mean and standard deviation of the properties for each model layer.  Figure 4(b) illustrates a model with varying permittivity that closely fits the observed data in Figure 4(a). This model represents the lowest dielectric values possible to fit the observed data because this model assumes that all layers are the same thickness, 1.7 m, which produces the maximum reflection possible for a given layer permittivity. Under this assumption, we find that all of the reflectors can be explained with dielectric permittivities ranging from 3.5 to 6.0. Due to increasing radar energy losses with depth, the dielectric contrasts must increase with depth to explain deep bright reflectors. Figure 4(d) shows a model with varying thickness that closely fits the same observed data. This model represents the minimum possible thickness that the layers could be. This model assumes that all layers are composed of the same compacted dust lag material with a relative dielectric permittivity of 8.8 alternating with layers of pure water ice (dielectric permittivity of 3.1). Under this assumption, we find that all of the reflectors can be explained with dust layer thicknesses of less than 0.5 m. Since there is loss as the wave propagates with depth, bright deep reflectors must be caused by thicker dust layers than those causing bright shallow reflectors. Figure 5 summarizes the modeling results from all 20 modeled radargrams. Additionally, Figure 5 converts the model dielectric permittivities into dust fraction values via Equation (17). Due to noise and the varying relative strength of the first surface return, to which all subsequent reflectors are normalized, the modeled results for each radargram are different. Thus, Figure 5 depicts the mean model result of all 20 radargrams and the corresponding standard deviation. Because the reflection events for each radargram vary slightly in timing due to ionospheric interference with SHARAD radar signals, the layers also have an uncertainty in their depth placement. Additionally, we find that as a layer increases in reflectivity, its model parameter result is more uncertain. This is because, for example, a small increase in thicknesses from 0.1 m to 0.2 m creates a large increase in reflectivity, but the same small increase from 0.4 m to 0.5 m creates less of an increase in reflectivity. The same argument holds for permittivity values. In general, we observe that there are three packets of reflectors. The middle packet of reflectors requires fewer reflective layers, and therefore can be composed of thinner material, or material with less permittivity.

Discussion
We find that thin dust-rich layers are a valid model for SHARAD reflections in the NPLD. Figure 5 demonstrates that the SHARAD reflectivity can be explained either by layers of varying thickness, or layers with varying permittivity. It is likely that the true model is a combination of both thickness and dielectric permittivity variation. However, by separating the two limiting cases, we arrive at two important results, the minimum dust concentration possible for the layers (Figure 5(a)), and the minimum possible thickness for the layers (Figure 5(b)). These provide excellent constraints for future studies wishing to relate the NPLD to the ancient Martian climate. If it is true that the layers within the NPLD were created by less ice deposition relative to dust deposition, and if a relationship is discovered that relates the dust fraction percent to the warmth and duration of a Martian obliquity cycle, then the dust fraction results in Figure 5(a) can characterize the warmth and variability of the past Martian climate. Alternatively, if it is true that warm Martian climates cause ice ablation on the poles that leave behind dust lag deposits, and if a suitable relationship is discovered between dust lag thickness and deposition time, then the thickness results in Figure 5(b) can provide an estimate of the period of time that the climate was warm.
Despite differing methodology, our minimum layer permittivity results in Figure 5(a) agree well with the results from Lalich et al. (2019), and suggest that some layers must be at least ≈ 50% dust. Our minimum layer thickness results in Figure 5(b) indicate that dust lag deposits formed from sublimation must be on the order of 30 cm thick. Although plausible, some studies (e.g., Skorov et al. 2001) suggest that lag deposits this thick could not form because even a lag of a few millimeters would be sufficient to insulate further ice from sublimating. Further research is necessary to determine how thick dust lag deposits may become. Also, although we choose a maximum dust lag layer relative permittivity of 8.8, which is consistent with the rock permittivity value chosen in previous modeling studies (e.g., Nunes & Phillips 2006;Lalich et al. 2019), there is not yet conclusive data to support this value. If it was found that a dust lag could have a higher permittivity, then the layers could be thinner than those we present. The more likely scenario, however, is that a dust lag deposit would have a lower permittivity, in which case, the layers would need to be thicker than those we present, which is consistent with our thickness values representing a lower bound.
Along with radar reflection data, high-resolution images of outcrops on the edge of the north polar ice cap also reveal layering. However, attempts to correlate radar reflections with images of outcrops have been inconclusive (Christian et al. 2013). Such studies have found that the NPLD is likely composed of layers that are anywhere from less than a meter to more than 10 m thick. If the scale of layers observed in imagery are indicative of that of the layers creating the radar reflections, then the thickness of the layers generating radar reflections could be much greater than those considered in our modeling effort. However, the thicker the beds are, the more lossy the NPLD as a whole would be due to a larger portion of the total volume affected by more lossy materials. If all of the layers are at a thickness close to the resolution limit of SHARAD, then the effective loss of the NPLD may be at odds with the effective loss values determined in Grima et al. (2009).
Our study underlines the importance of considering an Nlayer full-radar-pulse simulation rather than simply interpreting the reflectivity of layers individually. Even when assuming that simple thin layers are the cause of the reflection events, there is nuance in the correct model layer parameter. In general, a brighter reflector requires either a thicker bed, or a bed with a higher dielectric contrast. However, there are some locations where a dimmer reflector requires a greater contrast. For example, in Figure 4(a), the reflectors at 1211 m, 1255 m, and 1389 m all produce approximately the same reflection power, but require significantly different dielectric contrasts, a nonintuitive result. These scenarios typically occur in layers that are closely spaced (<50 m), and are a result of reflected pulses interfering with each other. This demonstrates that the interference between reflections from adjacent layer interfaces affects the observed reflected power. Indeed, in some cases, faint reflectors can occur that are entirely the result of interference or multiple reflections, such as the labeled event ii in Figure 4(a). Also, even though the models in Figures 4(a) and (b) produce approximately the same data, they do not have the same multiples or interference patterns because the models' layers have different locations and thicknesses. This result emphasizes the need to consider the full physics of a radar pulse reflecting from layered media when interpreting the composition of radar reflectors. Although we assume that each reflector is the result of one layering interface, this may not be the case. There is evidence that the geometry that produces the radar reflections is more complicated than single layers. Radar reflections could be caused by two or more closely spaced thin layers. Additionally, a layer whose thickness is on the cusp of SHARAD's resolution may produce one or two distinct reflections. One such ambiguous reflection is demonstrated by event i in Figure 4(a), where the denoted reflector appears to have two peaks. Campbell & Morgan (2018) have shown promising results toward distinguishing layer thickness and complexity by analyzing the lowand high-frequency portions of SHARAD's pulse separately. Future study combining their method with the layering modeling techniques in this paper may prove it possible to distinguish between thick reflectors and multiple closely spaced reflectors.
Our results are also influenced by the reflectivity of the surface. We assume that the polar cap surface is a clean boundary between Martian atmosphere and ice. If the surface has significant dust overlying the ice that raises the reflectivity of the surface, then our subsurface layer property results are underestimated. If the surface has significant carbon dioxide deposits that lower the reflectivity of the surface, then our subsurface layer property results are overestimated. Future investigation into the surface properties on the polar cap could help refine the treatment of the surface reflection in radar models.
We present our modeling software as a tool for other radar researchers. Forward modeling is valuable for testing subsurface layering hypotheses. However, there are limitations to the applicability of this modeling code. Geophysical data is inherently nonunique, and any given data set likely has multiple models that can satisfactorily explain the data. We present our methodology as a tool to test composition hypotheses. The modeling approach presented in this paper requires extensive human input and analysis through a trial-and-error like approach. Future work could allow more objective results using inverse theory. Nonlinear inversion algorithms could allow impartial recovery of model parameters, but also would require additional inputs and assumptions (e.g., Oldenburg & Li 2005). Finding the optimal model fit requires searching the parameter space for a global minimum. Unfortunately, the space has many local minima, as demonstrated by being able to find two models with good fits that are significantly different. Future work should focus on finding input constraints and nonlinear global-minimum search functions.

Conclusion
Previous studies have suggested that the Martian polar layered deposits are composed of thin dust-rich layers between nearly pure water ice. Using our open-source forward-modeling program, we demonstrate that this is a plausible model for the Gemina Lingula region of Marsʼs north polar layered deposits. We model two cases, thin beds with varying thickness but all with the same relative permittivity of 8.8, and 1.7 m thick beds with varying dielectric permittivity. We find that the observed reflections can be explained by either scenario, varying thickness or varying dielectric permittivity, and is likely a combination of both. If reflection power is driven by the varying thickness of thin beds with high permittivity, we find that all layer reflection events within Gemina Lingula can be explained by layers with thicknesses of 0.5 m or less. If reflection power is driven by the permittivity of the layers, we find that all reflections can be explained by relative permittivities of less than 6.0. More importantly, our results provide lower bounds on the thickness and dust fraction of layers within Gemina Lingula in the NPLD. For each reflection event in the >1 km ice cap, we provide a lower bound on the layer thickness and dust volume fraction. Our forward-modeling software adds a tool for the radar sounding community to test layering models, and our findings support the thin dust layer hypothesis, providing new layer composition constraints for Mars climate researchers.