Retrieval of the thickness of undeformed sea ice from simulated C-band compact polarimetric SAR images

In this paper we introduce a parameter for the retrieval of the thickness of undeformed first-year sea ice that is specifically adapted to compact polarimetric (CP) synthetic aperture radar (SAR) images. The parameter is denoted as the “CP ratio”. In model simulations we investigated the sensitivity of the CP ratio to the dielectric constant, ice thickness, ice surface roughness, and radar incidence angle. From the results of the simulations we deduced optimal sea ice conditions and radar incidence angles for the ice thickness retrieval. C-band SAR data acquired over the Labrador Sea in circular transmit and linear receive (CTLR) mode were generated from RADARSAT-2 quad-polarization images. In comparison with results from helicopter-borne measurements, we tested different empirical equations for the retrieval of ice thickness. An exponential fit between the CP ratio and ice thickness provides the most reliable results. Based on a validation using other compact polarimetric SAR images from the same region, we found a root mean square (rms) error of 8 cm and a maximum correlation coefficient of 0.94 for the retrieval procedure when applying it to level ice between 0.1 and 0.8 m thick.


Introduction
Sea ice covers about one-tenth of the world ocean surface and significantly affects the exchanges of momentum, heat, and mass between the sea and the atmosphere.Not only is sea ice extent a significant indicator and effective modulator of regional and global climate change, but sea ice thickness is also an important parameter from a thermodynamic and kine-matic perspective (Soulis et al., 1989;Kwok, 2010).The decline of sea ice extent recently observed in the Arctic, e.g., is linked with a decrease of ice thickness and increasing fractions of seasonal ice areas (e.g., Kwok et al., 2009).Measurements of sea ice thickness are compared with model results to control and validate the model capabilities of reproducing recent and predicting future trends of sea ice conditions in the Arctic (e.g., Laxon et al., 2013).Although sea ice thickness is only several meters at most, it forms an effective thermal insulation layer due to its high albedo and low thermal conductivity, leading to a significant reduction in the heat flux from the ocean to the atmosphere, especially in winter ( Vancoppenolle et al., 2005).Besides investigations focusing on the entire Arctic or Antarctic region, other studies analyze ice thickness variations on local scales to improve regional ice thickness retrievals (e.g., Haapala et al., 2013).Operational services charged with providing sea ice maps and forecasting ice conditions for marine transportation and offshore operations need near-real-time regular information about local and regional ice thickness distributions.The use of sensors providing high spatial resolutions on the order of 100 m or better for ice thickness retrieval, such as synthetic aperture radar (SAR), is an important topic of recent research (Dierking, 2013).
Unfortunately, the sea ice thickness distribution is also one of the most difficult parameters to measure.The most direct and accurate measurement technique is in situ drilling with an ice auger.Although it provides data with sufficient accuracy (in the range of centimeters), it is time-consuming and spatially limited.Therefore, this method is mainly used for the calibration of other sensors or methods.To obtain ice thickness distributions at larger spatial scales, remote sensing methods are requisite tools.There are generally different strategies.
1. Measurements of ice draft using upward-looking sonar on ocean moorings or submarines (Wadhams, 1980;Behrendt et al., 2013) from which thickness is estimated based on assumptions about buoyancy, ice density, and snow load (e.g., Rothrock et al., 1999) are used.Such data provide information about detailed temporal thickness variations (daily or even hourly) at a fixed location.An example for using in situ measurements of ice thickness from the New Arctic Program initiated by the Canadian Ice Service (CIS) starting in 2002, and sea ice draft measurements from moored upward-looking sonar (ULS) instruments in the Beaufort Gyre Observing System for testing a method of ice thickness retrieval from optical methods is provided by Wang et al. (2010).
2. Measurements of sea ice freeboard (i.e., the part of the ice above the water level) plus snow layer thickness with laser altimetry (e.g., Wadhams et al., 1992;Dierking, 1995) are used.From such data, the average ice thickness can be estimated, or the probability density function (PDF) of ice freeboard can be converted to a PDF of ice thickness.However, the estimation of ice thickness from freeboard data is less reliable than from ice draft because of a relatively stronger impact of errors in the freeboard measurements (Goebell, 2011).
Although ULSs and EMSs have all contributed greatly to our knowledge about ice thickness distributions on local and regional scales, such data can only be obtained at specific locations over a limited time period.Satellite remote sensing, on the other hand, is useful to monitor ice thickness variations regularly over much larger areas.On a still experimental basis, data of L-band passive microwave sensors, such as for example the Soil Moisture and Ocean Salinity mission (SMOS) radiometer, have been employed to retrieve thickness of sea ice thinner than about 0.5 m.The limitation of this approach is that it is only possible for very high (almost 100 %) sea ice concentration and in cold freezing conditions (Tian-Kunze et al., 2014;Huntemann et al., 2014).A space-borne altimeter has been used primarily to map ice thickness, and to monitor and study trends in thickness changes.The capabilities of laser and radar altimeter systems (such as CryoSat-2 and ICESAT) of measuring ice freeboard have been extensively investigated during the last decade (e.g., Kwok and Cunningham, 2008;Kwok et al., 2009;Laxon et al., 2013).Compared with radiometers, which only collect data at a coarse spatial resolution of a few to tens of kilometers (e.g., 25 km for the Special Sensor Microwave Imager (SSM/I) 37 GHz data), the spatial resolutions of radar altimeter systems are about 250 m alongtrack and 1.5 km across-track for CryoSat-2, and the footprint for ICESAT is about 70 m diameter.The sea ice products derived from altimeters usually focus on large-scale spatial and temporal variations.While the large-scale ice thickness product is important for climate research, the support of marine navigation and offshore operations in polar areas are crucially dependent on precise and reliable sea ice thickness maps with spatial resolutions better than 1 km.
Space-borne synthetic aperture radar (SAR), which operates in the microwave frequency band, provides all-weather and day-night high-resolution imagery (within a range of 1-100 m) with 1-3 days' temporal coverage.Hence, SAR is in general very useful for operational mapping tasks on regional and local spatial scales (Dierking, 2013).The disadvantage of SAR systems is that higher spatial resolutions are linked with a limited coverage between 10 and 500 km, compared, for example, to more than 1000 km for passive microwave radiometers.SAR measures the intensity of the radar signal backscattered from the ice surface and volume at different polarizations.The backscattered intensity depends on the dielectric constant of the ice and small-scale (mmdm range) ice properties such as ice surface roughness and air bubble fractions and sizes.If at least two polarizations are measured simultaneously, the SAR, which is a coherent device, can also provide the phase difference between the differently polarized channels.The most recent SAR sensors have polarimetric capabilities.A fully polarimetric radar transmits and receives both linear horizontal (H) and vertical (V) polarized electromagnetic waves.Amplitude and phase information of the backscattered signal are recorded for four transmit/receive polarizations (HH, HV, VH, and VV).This mode is commonly referred to as "quad-pol".Quad-pol scenes are usually acquired at very high spatial resolution.A RADARSAT-2 quad-pol scene has a spatial resolution of 4.7 m (slant range) × 5.0 m (azimuth) at a swath width of 25/50 km.Dual-pol scenes contain two polarimetric channels (e.g., HH and HV or VV and VH).In operational ice-charting services, dual-pol scenes are preferred because of their wider areal coverage (Geldsetzer et al., 2015).The RADARSAT-2 ScanSAR wide mode, e.g., can have a swath width of 500 km with 160-72 m (ground range) × 100 m (azimuth) resolution.Despite their currently very limited coverage, the quad-pol images are important information sources to understand the scattering mechanisms of sea ice.
Recently a number of investigators noted correlations between ice thickness and the co-polarization ratio, which is the ratio of measured intensities at VV and HH polarization (here we use VV / HH).The sensitivity between co-polarization ratio and thin ice thickness was first demonstrated by Onstott (1992), based on C-band radar data from the eastern Arctic region.Kwok et al. (1995) estimated the thin ice thickness (0 to 0.1 m) from L-and C-band fully polarimetric airborne SAR data acquired over the Beaufort Sea.Their approach included the training of a neural network.L-band polarimetric characteristics of ice in the Sea of Okhotsk were investigated by Wakabayashi et al. (2004), and the L-band co-polarized ratio was used to estimate ice thicknesses between 0 and 2 m (their Fig. 13).The investigation was further extended to other sensors, e.g., to the airborne Pi-SAR (X-and L-band data from the Sea of Okhotsk; Nakamura et al., 2009a;Toyota et al., 2009) and to Envisat ASAR, using radar intensity and ice thickness data from 0.2 to 2.5 m, the latter acquired from a research vessel in Lützow-Holm Bay, Antarctica (Nakamura et al., 2009b).The good correlations were attributed to the fact that the co-polarized ratio values are sensitive to the dielectric constants of the ice surface layer which change due to the process of desalination during ice growth.The relationship between relatively thick multiyear ice (thickness between 2 and 5 m), on the one hand, and co-polarized correlation and cross-polarized ratio HV/HH or VH/VV, on the other hand, was also investigated in the Arctic Ocean employing RADARSAT-2 and TerraSAR-X data (Kim et al., 2012).They found that the degree of depolarization is linked to the thickness of the multi-year ice as ice surface roughness increases and salinity decreases.
Although the above-mentioned parameters derived from polarimetric SAR imagery have shown the potential for estimating sea ice thickness under certain conditions, polarimetric SAR data can presently only be acquired at limited swathwidths.The quad-pol mode on RADARSAT-2 has a swath width of only 25-50 km, as mentioned above.The swath width of the VV/HH dual-polarization StripMap mode on TerraSAR-X is 15 km.Therefore, they are insufficient for operational use which requires a large-scale coverage (Scheuchl et al., 2004).The limited swath width also restricts scientific investigations to local domains.An alternative is to use compact polarimetry.
The methods of generating compact polarimetric (CP) information (explained below) are based on receiving data at two different polarizations (Souyris et al., 2005;Raney, 2007).Compared with the "traditional" dual-polarization modes described above, CP data include a greater amount of polarization information (but less than quad-polarization data).They can cover much greater swath widths compared to quad-polarization modes due to reduced power consumption and data storage requirements.
The term "CP system" refers to a unique polarization in transmission and coherent dual-orthogonal polarizations in reception.There are three different CP configurations (Nord et al., 2009).The first one is the π /4 mode, with a slant linear transmission and horizontal (H) and vertical (V) receptions (Souyris et al., 2005).The second is the dual circu-lar (DC) mode, i.e., transmitting at a single circular polarization and receiving two orthogonal circular polarizations.The last approach is circular transmit and linear (H and V) receive (called CTLR mode).Among these three compact polarization modes, the latter has been ranked to be the most promising in terms of performance and receiver complexity.The current Indian RISAT-1, the Japanese ALOS-2, and the planned Canadian RADARSAT Constellation Mission (RCM) also support the CTLR mode.According to the description in Geldsetzer et al. (2015), the coming CTLR mode of RCM will be particularly tailored to sea ice applications by offering a medium-resolution mode with a swath width of 350 km and a resolution of 50 m, a low-noise mode with the same swath width and a resolution of 100 m, or a low-resolution mode with a swath width of 500 km and a resolution of 100 m.Hence, the CTLR modes of RCM are well suited for operational sea ice monitoring.
However, one apparent disadvantage of the CP mode as compared to dual-or quad-polarization mode is the fact that the HH, VV, and HV signal combinations are not directly measured.This means that the co-polarized ratio (Wakabayashi et al., 2004;Nakamura et al., 2009a;Toyota et al., 2009) and the cross-polarized ratio (Kim et al., 2012) which are often used as an ice thickness proxy cannot be directly calculated from CP-mode SAR data.Although CP SAR images have been used to distinguish sea ice types (Dabboor and Geldsetzer, 2014;Charbonneau et al., 2010;Geldsetzer et al., 2015), to our knowledge there have been no published studies on its use for ice thickness detection in the open literature until now.Therefore, in this study, we considered the CTLR mode and developed an approach to directly retrieve the thickness from CP SAR data (hereafter we assume that the CP SAR is operated in CTLR mode).The paper is organized as follows: in Sect. 2 we introduce a new parameter to estimate ice thickness and demonstrate its sensitivity to different ice parameters by numerical modeling in Sect.3. In Sect.4, an empirical relationship based on a comparison of CP-SAR signatures with ice thickness data obtained from electromagnetic induction sounding is presented, and the retrieval performance of this algorithm is described.Further discussions and conclusions are presented in Sect. 5.

Full polarimetry and compact polarimetry
The full polarimetric radar scattering return can be represented by the scattering matrix S: where S pq denotes the p transmit and q received linear polarization.In the monostatic case and considering that reciprocity can be assumed for sea ice and snow, S HV = S VH .
X. Zhang et al.: Retrieval of the thickness of undeformed sea ice We use the coherency matrix T to evaluate the second-order statistics of the scattering matrix S. The coherency matrix T formed from the elements of the scattering matrix S is where * denotes the complex conjugate and denotes the ensemble average.
We consider the CTLR mode for which the scattering vectors are given by (e.g., Nord et al., 2009) (3) As usual, R denotes that the transmitted polarization is right circular, while H and V stand for the linear reception.We set From Eq. (3) it then follows that The terms | H | 2 and | V | 2 can be expressed as Under the assumption of reflection symmetry, the cross-and co-polarized scattering coefficients are uncorrelated.This assumption is reasonable for snow and sea ice surfaces at various frequencies and for different spatial scales (Souyris et al., 2005).Hence and Eq. ( 6) can be rewritten by the elements of coherency matrix T:

X-Bragg model and X-SPM model
According to the results obtained by the Cold Region Research and Engineering Laboratory (CRREL'88), the typical ranges of root mean square (rms) height and correlation lengths for smooth level sea ice are 0.02-0.143and 0.669-1.77cm respectively (Fung, 1994).For C-band SAR, the small perturbation method (SPM) can be applied for explaining the surface scattering characteristics from smooth level sea ice.By doing so, the underlying assumption is that the received radar signatures are typical for Bragg scattering.However, the SPM fails to describe cross-polarization and de-polarization effects that are observed in real SAR data.In order to overcome these limitations and to widen the SPM range of validity, an extended Bragg model (termed X-Bragg model) was presented by Hajnsek et al. (2003).In the X-Bragg model, the scattering surface is composed of rough randomly tilted facets that are large with respect to the wavelength but small with respect to the spatial resolution of the sensor (for RADARSAT-2 fine-quad mode, the wavelength is 5.6 cm and the resolution is 8 or 25 m).Scattering from each rough facet is evaluated by employing the SPM, whereby for the facets, a random tilt is assumed which causes both a random variation θ of the incidence angle θ and a random rotation β of the local incidence plane around the line of sight.In the X-Bragg model, the random incidence angle variation θ is ignored, and the incidence plane angle of rotation β is assumed to be uniformly distributed in an interval (−β 1 , β 1 ), where the parameter β 1 is used to characterize the large-scale roughness (del Monaco et al., 2009).
In order to improve the range of validity of the X-Bragg model, different approaches (termed X-SPM model) were proposed by del Monaco et al. (2009) and Iodice et al. (2011).In those studies, more realistic distributions of β and θ were derived by assuming that the range and azimuth facet slopes are Gaussian random variables.The coherency matrix of the X-SPM model (T X-SPM ) after ensemble averaging over the local incidence angle θ l and rotation angle β can be expressed as follows (del Monaco et al., 2009): Here, < > θ l means averaging over the local incidence angle; θ l which is used to characterize the random slope varia- tions of the facets; < > | β means averaging over the rotation angle β; ρ includes small-scale roughness effects, and σ is the standard deviation of the surface slope which is a Gaussian random variable.Erfc{ q } is the complementary Gauss error function.R P and R S are the Bragg scattering coefficients perpendicular and parallel to the incident plane, respectively.Both are functions of the complex permittivity ε and the incidence angle θ (Iodice et al., 2011): In the paper by del Monaco et al. (2009) it is demonstrated that the X-SPM model coincides with the X-Bragg model when the standard deviation of the surface slope is zero, and that the X-Bragg model can only be applied for standard deviations of the surface slope σ < 0.1.When σ > 0.1, the effects of incidence angle fluctuations, which are ignored in the X-Bragg model, are significant (del Monaco et al., 2009).
Because of its wider range of validity, we used the X-SPM model in our study.

Inversion model
For ice thickness retrievals we propose to exploit the ratio between | V | 2 and | H | 2 (here denoted as the CP ratio).
The CP ratio can be written as (see Eq. 4) By relating the CP ratio to the elements of the coherency matrix given for the X-SPM we obtain CP ratio= Equation ( 12) shows that the CP ratio is controlled by ensemble averages of the difference and sum of the Bragg coefficients with respect to the incidence angle.From del Monaco et al. (2009), the probability density function for cos θ l is a normal distribution with mean cos θ and standard deviation equal to σ sin θ .After averaging over variations of the local incidence angle θ l , the CP ratio is dependent on the dielectric constant of the surface ε, the incidence angle θ , and the standard deviation of the surface slope σ .By using the model of del Monaco et al. (2009), the results of SAR measurements can be better explained than with the SPM.We calculated the CP ratio as a function of the standard deviation of surface slope σ , assuming ε = 3.9 + j 0.15 which is suggested in Fung and Eom (1982) for first-year sea ice.The results show that the CP ratio increases with increasing standard deviation of the surface slope at fixed incidence angles and with increasing incidence angle at fixed σ (Fig. 1).The relationship between the CP ratio and the dielectric constant is presented in Fig. 2. When the incidence angle is constant, the CP ratio reveals monotonically increasing values with increasing dielectric constant.A similar trend can also be found in the co-polarization ratio (Iodice et al., 2011).With respect to our simulated results shown in Figs. 1 and 2, it is important to note that the proposed parameter CP ratio is sensitive to the variation of the dielectric constant and almost insensitive to surface slope variations if σ < 0.15.
For our analysis we use the fact that a dry snow layer is transparent at C and L frequencies (meaning that our method is only applicable under freezing conditions), and we do not consider metamorphosis of the basal snow layer due to brine wicking effects or due to melt-freeze cycles.We focus on undeformed Arctic young and first-year ice for which volume scattering is low because of the relatively high ice salinity, which means that the ice surface is the dominant scattering source.Then the backscattering coefficients depend on the small-scale surface roughness and the dielectric constant of the ice surface.Desalination of the ice occurs parallel to its growth due to brine drainage (Kovacs, 1996).The desalination process causes a decrease of the dielectric constant.Hence the basic idea of our method for retrieving ice thickness is to relate changes of the dielectric constant to ice thickness growth.Because the CP ratio is sensitive to the variation of the dielectric constant, it is well-suited for detecting changes of the ice thickness of smooth first-year level ice.

Forward scattering model
In this section, we describe the combined use of an ice growth model and an electromagnetic scattering model for level sea ice to study sensitivities of the CP ratio to different ice and radar properties.We applied the scattering model proposed by Nghiem et al. (1995) to simulate the sea ice volume scattering and absorption by brine inclusions.The surface contribution was calculated with the polarimetric two-scale model (Iodice et al., 2011(Iodice et al., , 2013) ) and incoherently added to the volume term.
The sea ice scattering model configuration is presented in Fig. 3.Note that we do not explicitly include a snow layer (see also Sect.2.3).The effects of a dry snow layer are that (1) the dielectric contrast between ice and snow is lower than between ice and air, hence the reflectivity of the ice surface is lower; (2) the radar wavelength in the snow is shorter than in air, hence the ice surface appears rougher to the radar; (3) the incidence angle gets steeper (depending on the dielectric constant of the snow), which (relatively) causes a stronger backscattering.Since we carry out simulations with different dielectric constants (by varying temperature and brine volume fraction), surface roughness parameters, and radar incidence angles, the results obtained without snow can be transferred to cases with dry snow layers.
In our model, the uppermost layer is air with permittivity ε 0 ; the lowermost medium is seawater with complex permittivity ε 2 , both enclosing the ice layer.The sea ice background is assumed to be pure ice with complex permittivity ε i .The complex permittivity of brine inclusions is ε b , and their fractional volume is f v .The relative permittivity of the sea ice ε eff is a function of the volume fraction of brine inclusions (Arcone et al., 1986;Vant et al., 1978).The ice surface roughness is described by the correlation length l, rms height s, and the standard deviation of surface slope σ .
The thickness and surface temperature of the sea ice layer are H and T 0 , respectively.Lastly, the magnetic permeability of free space is µ 0 .Thickness and permittivity of sea ice are subject to dynamic changes during the ice growth process.The small-scale surface roughness (on a centimeter scale) may also vary temporally and spatially.This, however, can hardly be measured in the field with sufficient spatial density over larger areas.Here we do not consider deformation processes causing surface roughness components on the order of meters.Furthermore, we assume that the scattering contribution of the ice-water interface can be neglected because of the relatively high salinity of Arctic young and first-year ice.Very thin ice, for which reflections of the radar waves between surface and bottom have to be considered, is excluded from this study.In our simulations, we do not take snow cover into account.We restrict our analysis to temperatures well below freezing point, which means that a dry snow layer would change the incidence angle and the dielectric contrast at the ice surface.In the case of the ice growth simulations described below, the snow has an insulating effect that changes the rate of ice thickness growth.Hence, various scenarios can be constructed, which is beyond the scope of this paper, which we regard as a first step towards developing a methodology for ice thickness retrieval using CP SAR.
For ice growth simulations we use a 1-D thermodynamic model developed by Maykut (1978Maykut ( , 1982) ) based on the energy balance equations at the atmosphere-ocean boundary.The balance of the heat fluxes at the upper surface of the ice can be expressed as where F r is the incident short wave radiation, αF r is the short wave radiation reflected by ice, and α is the albedo.I 0 is the amount of shortwave radiation absorbed in the interior of the ice layer, F L is the incoming long wave radiation, F E is the long wave radiation emitted by the ice, F s is the sensible heat flux, and F E is the latent heat flux.The last term F C is the upward conductive heat flux that is the heat from the bottom interface conducted through the ice to the upper surface.We assume that the temperature at the ice-water interface is at −1.8 • C. The equations and parameters used in this study are listed in Table 1.
Substituting the equations and parameters listed in Table 1 into Eq.( 13) and using the Newton-Raphson iteration method, the sea ice surface temperature T 0 is obtained.Once T 0 is known, F E , F s , F E , and F C can be easily calculated.A linear temperature profile within the sea ice layer is assumed.For volume scattering and absorption calculations we use a mean ice temperature (T ) calculated from the melting temperature at the ice-water interface temperature (T b = −1.8• C) and the ice surface temperature (T 0 ).Furthermore, the thickness H (cm), density ρ (kg m −3 ), brine volume fraction f vb , and permittivity ε eff of sea ice, which are directly related to the volume scattering and absorption in the ice, are obtained by the equations given in Table 2.The ice The long wave radiation ζ is the Stefan-Boltzman constant (unit: w/(m 2 K 4 )); T 0 is the surface temperature of sea ice (unit: K); T a is the air temperature (unit: K); ε i is the emissivity of sea ice; ε a is the emissivity of atmosphere; e is the water vapor pressure at T a (unit: HPa).(Maykut, 1978) k=0.0017 The sensible heat flux (Cox and Weeks, 1988) ρ a is the air density (unit: kg/m 3 ); C p is the specific heat at constant pressure (unit: J/(kg•K)); C s is the sensible heat transfer coefficient; u is the wind speed; The latent heat flux F e =ρ a LC e u(q a -q 0 ) (Cox and Weeks, 1988) L is the latent heat of vaporization (unit: J/kg); The albedo of sea ice Cox and Weeks, 1988) H is the sea ice thickness (unit: cm).
The upward conductive heat flux and Weeks, 1988)   Fukusako,1990) is the sea ice growth rate when ice thickness is H (unit: m/s); Ice thickness is the sum of ice growth rate.ΔTime is the time lag (unit: hour).
The sea ice density and brine volume fraction (Cox and Weeks, 1983) ρ is sea ice density (unit: kg/m 3 ); f vb is the relative brine volume fraction ρ i is pure ice density (unit: kg/m 3 ); T i is the temperature of sea ice (unit: °C); T a is the air temperature (unit: K); S i is ice salinity.The functional forms of F 1 and F 2 can be found from the work of Cox and Weeks (1983).(Arcone et al., 1986;Vant et al., 1978) f vb is the relative brine volume fraction.

Simulation results
To assess theoretical possibilities and limitations of icethickness measurements by CP ratio, we simulated the evolution of ice growth for given temperature and wind conditions based on the growth model described in Sect.3.1.The air temperature and wind speed were set to −12 • C and 10.5 m s −1 , respectively, throughout this simulation, based on reports from the field measurements that are described in Sect. 4 below.The simulation started at an initial ice thickness of 1.0 cm.A finite difference scheme was used to calculate the increase of ice thickness at every 1 h step.After executing about 25 days' simulation, the following parameters were extracted as a function of time to drive the sea ice scattering model: ice permittivity ε eff , thickness of ice layer H , and volume fraction of brine inclusions f vb .For evaluating the rough surface scattering contribution, we took roughness data reported in Onstott (1992, Table 5-3) who listed them for different stages of ice growth: (1) s = 0.031 cm and l = 1.26 cm (kl = 0.035, ks = 1.4 for the radar frequency of 5.4 GHz, k -wavenumber) for dark nilas, (2) s = 0.12 cm and l = 1.45 cm (kl = 0.14, ks = 1.6) for light nilas, and (3) s = 0.11 cm and l = 0.54 cm (kl = 0.12 and ks = 0.6) for smooth first-year ice.We note that we will use these roughness values for first-year ice in general, considering the large variability of small-scale ice surface roughness.The values are in the validity range of the original Bragg scattering theory and should hence be fully covered by the X-SPM model presented in Iodice et al. (2011).The standard deviation of the large-scale slope σ ranges according to the validity range of the X-SPM model (Iodice et al., 2011).
At this point we note that a systematic relationship between small-scale surface roughness and ice thickness has never been reported.Weathering effects, melt events, and snow metamorphism influence the millimeter-to-centimeter ice surface roughness to a highly variable extent, independent of ice thickness.As we will show below, the influence of the small-scale roughness on the CP ratio is moderate to low; hence the issue of varying small-scale surface roughness is not very critical.
Figure 4 illustrates the simulated sea ice thickness as a function of time and ice temperature, and the volume fraction of brine inclusions as functions of ice thickness.Figure 4 clearly shows that the volume fraction of brine inclusions reduces due to desalination processes as the ice thickness increases.
To investigate the dependence of the CP ratio on the radar incidence angle and ice thickness, the complex scattering coefficients (S HH , S VV , and S HV ) were computed for the Cband (5.4 GHz) at incidence angles of 20-60 • .Then the CP ratio was calculated from Eq. ( 12).The relationship between the CP ratio and sea ice thickness in case 3 (first-year ice

CP ratio
Figure 5.The relationship between the CP ratio and ice thickness at different incidence angles for C-band radar (x axis on a log scale).The incidence angle varies from 20 to 60 • .The small-scale roughness parameters are set to s = 0.11 cm and l = 0.54 cm (case 3), the standard deviation of the surface slope σ = 0.1.roughness conditions given above) and σ = 0.1 is shown in Fig. 5.It reveals that the CP ratio exhibits a monotonically decreasing trend with growing ice thickness at constant incidence angles.It should be noted that the sensitivity of the CP ratio to vertical ice growth is much higher at smaller ice thickness values up to approximately 0.4 m.This can be explained by fact that the ice salinity is calculated according to the relationship proposed by Cox and Weeks (1983).Their parameterization of salinity as a function of ice thickness reveals a discontinuity at a thickness of 0.4 m.
Figures 6 and 7 indicate the roughness dependencies of the CP ratio.In Fig. 6 the standard deviation of surface slope σ is varied from 0.05 to 0.4 and the small-scale roughness is fixed at the case 2 roughness condition (s = 0.12 cm, l = 1.45 cm).When σ is smaller, the effect of the variability of the icesurface slope on the sensitivity of the CP ratio to ice thickness is small; however, at larger values of σ , this effect becomes significant and weakens the capability of the CP ratio to estimate thickness.Given the same σ values, the magnitude of the CP ratio is higher at larger than at smaller incidence angles, while the sensitivity (given by the local slope of the curves) hardly changes, as depicted in Fig. 7, where we show examples for case 2 (moderate) ice roughness.The sensitivity as a function of ice thickness remains basically the same for all incidence angles.A larger magnitude of the CP ratio means that it is less affected by noise (see Eq. 11).From the results of these simulations, we expect that the proposed new parameter for thickness retrieval has a strong correlation with the thickness of smooth undeformed sea ice over all incidence angles, and the sensitivity is larger for thinner Figure 6.Sensitivity of the CP ratio to the standard deviation of the surface slope σ (x axis in log scale).The standard deviation of the surface slope σ varies from 0.05 to 0.4, while the small-scale roughness is fixed at s = 0.12 cm and l = 1.45 cm (case 2).The top panel is for the 20 • incidence angle and the bottom panel is for the 40 • incidence angle.
(< 0.4 m) than for thicker sea ice.At larger incidence angles, the reduction of the radar wavelength in a snow layer on top of the ice is not a critical issue, since the effect of the smallscale roughness on the CP ratio is low in this case.However, the snow layer also changes the incidence angle of the radar beam on the ice surface, which can have a considerable impact on the thickness retrieval, in particular at thickness values larger than 0.3 to 0.4 m where the slope of the curves theoretically decreases to a low value (Fig. 7).In practice, this limitation is less critical as we show below.On first-year CP ratio Figure 7. Sensitivity of the CP ratio to the small-scale roughness (x axis in log scale).The standard deviation of the surface slope σ is fixed at 0.1.Black, blue, red, green and cyan colors are for 20, 30, 40, 50, and 60 • incidence angles, respectively.In the legend, C1, C2, and C3 denote the three cases of small-scale surface roughness respectively (C1: s = 0.031 cm, l = 1.26 cm; C2: s = 0.12 cm, l = 1.45 cm; C3: s = 0.11 cm, l = 0.54 cm).sea ice, the bottom part of the snow layer can be saline due to brine wicking, possibly creating a dielectric interface within the snow, or resulting in brine volumes large enough to influence the radar backscatter (Barber and Nghiem, 1999;Galley et al., 2009).This may also affect the accuracy of the thickness retrieval using the CP ratio.Finally, we note that the model simulations include interactions between the ice surface and the ice-water interface, which result in oscillations of the CP ratio for an ice thickness < 0.16 m.In the field measurements discussed below, this effect was not observed.We assume that the actual ice thickness is rarely exactly constant over larger areas.

Field study
On 19-20 March 2011, a field program was conducted by the Department of Fisheries and Oceans Canada (DFO) along the mid-Labrador coast (Fig. 8) (Prinsenberg et al., 2012a).As part of the field survey, snow thickness and ice thickness were measured with a helicopter-borne sensor package which consists of a laser altimeter, an electromagnetic induction sounder (EMS), and a ground-penetrating radar (GPR).The laser altimeter provides the distance to the snow or ice surface, whereas the induction sounder measures the distance from the sensor to the ice-water interface.Hence the snow plus ice thickness can be obtained (Prinsenberg et al., 2012a, b).Comparisons with drill hole data showed that  3.
the ice thickness values derived from such soundings agree well within ±0.1 m over flat homogeneous ice (Haas et al., 2006;Prinsenberg et al., 2012b).The accuracy decreases over ridges and deformed ice, where the maximum thickness can be underestimated by as much as 50 % (Haas et al., 2006;Prinsenberg et al., 2012b).Snow thickness profiles were collected concurrently with a ground-penetrating radar (GPR) and the laser altimeter measurements.The ground-penetrating radar, which was operated at a frequency of 1 GHz, receives returns from the ice-snow and air-snow interfaces, though the return from air-snow surface is very weak.The laser altimetry is superior for defining the airsnow interface.Therefore, the combination of the GPR and laser altimetry allows the snow depth on sea ice to be retrieved.For a 1 GHz GPR system, the minimum detectable snow layer thickness is 0.12 m and the measurement error is 0.08 m in light dry snow (Lalumiere, 2006).By subtracting the GPR snow thickness measurements from the EMS snow plus ice thickness measurements, sea ice thickness can be estimated.

Data sets and data processing
All data are available on the website of DFO including pictures, notes, and reports of the survey (http://www.bio.gc.ca/science/research-recherche/ocean/ ice-glace/data-donnees-eng.php).
During the field survey, four C-band RADARSAT-2 quadpolarization images were acquired nearly coincident with the DFO airborne survey flight lines (Fig. 8).The RADARSAT-2 data were provided by the MacDonald, Dettwiler and Associates Ltd (MDA).Important SAR parameters are listed in Table 3.For our processing we used the RADARSAT-2 single-look slant range complex format as starting point.A speckle reduction filter (13 × 13 Lee filter) and radiometric calibration procedures were applied for the calculation of the scattering matrix.With the quad-polarization data, the CTLR compact polarimetry mode can be generated via Eq.(3).Subsequently the CP ratio was extracted by Eq. ( 11).Lastly, the geometric registration of the simulated CP SAR images (i.e., their representation in geographical coordinates) was performed based on longitude and latitude data provided in SAR metadata.
Figure 8 presents the ice condition at the study site, flight paths and four nearly coincident RADARSAT-2 fine quadpolarization images.Eight EMS profiles were measured within the coverage of the four SAR images, and the time differences between the SAR acquisitions and EMS flights are summarized in Table 4.The images in Fig. 8 show the RADARSAT-2 data overlain with the EMS flight tracks over the fast ice and drifting pack ice.According to the ice charts, the total ice concentrations in fast ice and pack ice regions are 10/10 and 9/10, respectively.The main ice type in land fast is first-year ice of 70-120 cm in thickness, and the drift ice region contains gray ice (10-15 cm thick), gray-white ice (15-30 cm), thin first-year ice (30-70 cm), and again firstyear ice, 70-120 cm thick.In the drifting ice region several openings can be seen in the SAR images.The extent of land-fast ice evolves in the offshore direction and can be visually separated from the pack ice.Most of the rougher land-fast ice is brighter in the SAR images than the thinner undeformed land-fast ice.According to the meteorological  (Prinsenberg et al., 2012a).Hence the snow can be regarded as dry.We also note that thermodynamically driven effects on the bottom snow layer such as brine wicking take place at temperatures higher than −7 • C ( Barber and Nghiem, 1999) which means that we can ignore them here for the freshly fallen snow.However, we do not have any information about elder snow layers changed by metamorphosis processes, which may have an influence on the effective backscattering signature; nor can we exclude the fact that sea ice flooding took place in some smaller areas.Figure 9 shows the ice thickness and snow depth profiles of the land-fast and drift ice, indicating that the ice freeboard was mostly above the water level.The histograms shown in Fig. 9 confirm that the land-fast mean ice thickness is smaller than the one of the drifting pack ice.The percentages of areas with snow thickness above 0.2 m for land-fast and drift ice are 26.4 and 18.2 % respectively.The flight profiles also show that there are deformed ice or ridges (ice thickness exceeded 2.0 m) in the survey field.
A direct comparison between SAR imagery and flight profiles' data may cause errors due to the time differences of the data acquisitions (the time difference between SAR and flight data is shown in Table 4).In addition, spatial differences may be caused by the different sampling and spatial resolutions of the measurement instruments.The sampling rate for the EMS and the laser is 10 Hz, which, given a typical helicopter survey speed of 80 mph, corresponds to a spatial sampling interval of about 3-4 m.While the footprint size of the laser is very small (several centimeters), the footprint of the EMS is around 20 m at a typical operation height of 5-6 m.For this experiment, the GPR was configured to a scan rate of approx- imately 30 scans per second.When flying at 60-80 knots, the ground sample spacing is approximately one sample per 1.0-1.5 m.Moreover, according to the DFO survey report, the floating ice drifted 1.4-1.8knots towards the southeast, as measured by ice beacons (Prinsenberg et al., 2012a).In order to mitigate the errors caused by time and spatial resolution differences, we used the following processing chain for linking SAR and airborne data.
1.The correction of the time difference was only implemented for the drifting ice region.The boundary between fast ice and drifting pack ice was taken from ice charts of the Canadian Ice Service (Fig. 8).Of the eight EMS profiles, P1, P2, P5, and P7 are in or near the landfast ice region, whereas P3, P4, P6, and P8 are from the drift ice zone.With an ice drift speed of 1.5 knots, and drift direction southeast taken from the DFO survey report and considering the respective time differences, the profiles P3, P4, P6, and P8 are shifted to their approximate positions at the acquisition time of the SAR images.The shifted profiles are presented in Fig. 8 (dotted line).It should be noted that 28 h passed between the acquisition times of the P8 and SAR data, and the corrected location of P8 is beyond the coverage of the SAR image.Hence P8 was discarded from further analysis.
2. The EMS (ice plus snow) thickness values below 0.1 cm were removed to consider the measurement accuracy of the EMS.Regions for which only EMS data but no GPR data are available were also removed.
3. Regions with GPR snow thickness values higher than 0.20 m were removed, because snow layers thinner than 0.20 m are nearly transparent to C-band radar waves, and the backscatter from the snow surface and volume can be neglected (Hall et al., 2006).
4. By combining the field survey data (ice charts and aerial photos), a visual interpretation of RADARSAT-2 SAR was made, and regions of open water, land, and deformed ice were masked in the SAR images.Land was identified using the coastal line; open water areas were interpreted via backscattering and texture.Deformed ice was brighter than level ice in single-polarization SAR images, and revealed a higher entropy, which was extracted using H /A/α decomposition (Scheuchl et al., 2002).We emphasize that in step 2, most open water areas are already excluded from further analysis.
5. For ice zones of 50 m in length, averages of different parameters were evaluated.Firstly, we used the H /A/α unsupervised Wishart classifier to segment the SAR images, and each patch was regarded a homogeneous ice area with respect to its radar signature.Then the snow thickness, snow plus ice thickness profiles were cut into flight track segments 50 m long.The CP ratio values were evaluated from the co-located, drift-corrected, segmented SAR images, provided that the 50 m segment contained a homogeneous piece of ice.The segment length of 50 m was chosen according to the spatial resolution of the SAR image.Range and azimuth spacing of a RADARSAT-2 fine quad-polarization product are 4.7 m × 4.9 m respectively.Since we applied a 13 × 13 window for speckle reduction (see above), the effective spatial resolution is about 50 m.For the averages along transects, 13 SAR pixels, 15 EMS samples, and 45 GPR samples were used.
6.The sea ice thickness was extracted from the averaged GPR snow depth and EMS snow plus ice thickness values.
7. Finally we calculated the CP ratio from Eq. ( 11) using the averaged complex backscattering coefficients.
This processing chain ensures that only level ice is considered for which the EMS system delivers reliable thickness data with an acceptable accuracy.The total length of the profile segments that we used in this study amounts to about 16 km (320 samples).Compared with the original data, almost 60 % of the data were discarded in this processing chain (step 1: 17 %, step 2: 10 %, step 3: 23 %, step 4: 10 %).

Ice thickness retrieval
To investigate the possibility of using the proposed polarimetric parameter CP ratio to estimate sea ice thickness from SAR images, we plotted ice thickness values obtained during the field campaign against the corresponding values of the CP ratio derived from the RADARSAT-2 images in Fig. 10 (using all 320 samples).It can be seen that at C-band, the CP ratio shows a negative trend relative to the ice thickness as the simulated results given in Sect.3.2 predicted.Figure 10 reveals that the highest sensitivity occurs between 0 and 0.5 m and saturates with thickness values exceeding 1.5 m.As shown in Figs. 5 to 7, the sensitivity should be smaller for ice thickness exceeding 0.4 m.However, the slope change of the curves at 0.4 m is not as abrupt as in the theoretical curves predicted in Sect.3.2.This can be presumably explained by the fact that we average over segments with different values of ice roughness parameters s, l, and σ .We also need to consider that the salinity-thickness parameterization proposed by Cox and Weeks (1983) includes a discontinuity in the slope of the salinity curve at a thickness of 0.4 m, which may not exist in reality.Since our data comprise different incidence angles (29, 42, and 49 • at the survey positions, Table 3), we can construct the relationships between ice-thickness and the CP ratio dependent on the incidence angle.We applied two different fits, a linear and a logarithmic one, to obtain an empirical relationship between the ice thickness and CP ratio.The best regression was obtained using a logarithmic function (Fig. 10).
www.the-cryosphere.net/10/1529/2016/The Cryosphere, 10, 1529-1545, 2016 where all data points (320 samples) are used to derive the empirical regressions in the thickness range from 0.1 to 1.8 m.The reason to include larger ice thickness values is that they can be measured with a larger accuracy, hence leading to a more robust relationship at least for the moderate thickness values between 0.4 and 0.8 m.However, to our knowledge the distribution of the CP ratio due to speckle has not been derived yet which makes it difficult to judge its variation.
The smallest values of the CP ratio observed are about 0.03, which may indicate the noise level of the CP ratio.The measured values of the CP ratio for ice thickness values > 0.2 m shown in Fig. 10 are lower than the theoretical computations.This can presumably be explained by the fact that underlying theoretical models are an oversimplification of the actual situation.We note that due to the limitation of sample points, the fit for 49 • incident angle is mainly determined by ice thickness values > 0.5 m.
We found that the level of the CP ratio increases as the incidence angle increases at a given value of the sea ice thickness.This observation compares well with the forward simulation studies as shown in Fig. 5.These high correlations enable us to derive reliable thickness information for smooth level ice from radar images, assuming winter conditions (dry snow, no brine wicking).The ice thickness can be estimated using an exponential function, which can be described as follows: where a and b are the coefficients of the exponential fit.
At the next stage, we focused on the RADARSAT-2 images no. 2 and no. 3 (which have the same incidence angle of 42 • ) to validate our method.Out of a total of 320 samples, 159 samples belong to images no. 2 and no. 3.According to the principle of independent sample tests, we divided these 159 samples into two data sets in an arbitrary way.The first set includes 79 samples that are used to fit the model for estimating ice thickness, and the second one comprises 80 samples that serve to retrieve ice thickness and compare the results with the data from the field campaign.The coefficients a and b of the empirical fit generated from the first data set are 0.068 and 0.077 respectively.Note that these coefficients are different from those derived in Eq. ( 14) from the same two SAR images because now fewer points could be used to derive the fit.The fitted curve and validation results are presented in Fig. 11a and b, respectively.The correlation coefficient for the fit shown in Fig. 11a is 0.93 for the thickness range from 0.1 to 1.8 m and 0.94 for the thickness range from 0.1 to 0.8.The rms error and the relative error between the observed and the estimated ice thickness, shown in Fig. 11b, are 12 cm and 20 % in the thickness range from 0.1 to 1.8 m, and 8 cm and 17 % for 0.1 to 0.8 m.The relative rms error implies, e.g., that the absolute rms error is 0.2 m at an ice thickness of 1.0 m (for the range 0.1 to 1.8 m). Figure 11b also demonstrates that the error of the retrieved ice thickness is very large at values > 0.8 m which is to be expected from the theoretical curves, considering the significantly decreased sensitivity of the CP ratio to larger ice thickness.

Discussion and conclusion
This paper provides a first analysis of sea ice thickness retrieval using compact polarimetric SAR.We developed a new parameter that we call the CP ratio to estimate the thickness of undeformed first-year level ice from C-band radar images, under dry snow conditions (snow depth < 20 cm).Numerical model simulations showed that this parameter is sensitive to changes of the dielectric constant that are linked to the growth of sea ice.We developed empirical relationships for the retrieval of level ice thickness from CP ratios.For the validation of our results we also employed RADARSAT-2 images for which thickness values were available.The optimal regression between the CP ratio and ice thickness was achieved with an exponential fit.The rms error was 12 cm, and the relative error amounted to 20 % for a thickness range between 0.1 and 1.8 m, and 8 cm and 17 % for the range between 0.1 and 0.8 m.This indicates that the proposed parameter is very useful for the retrieval of first-year level ice thickness between 0.1 and 0.8 m.
Since the thickness of deformed ice can be underestimated by the EMS measurements by as much as 50 or 60 % in the worst cases, we could only study the case of level ice.The capability of CP SAR to retrieve the thickness of deformed ice, which reveals a larger variation of large-scale roughness with respect to the sensor resolution, needs to be further discussed and studied.
Although our tests are performed on a limited sample of images, our findings demonstrate that the C-band compact polarimetric SAR has a potential for sea ice thickness retrievals over level first-year ice covered by a thin dry snowpack.The issue of environmental factors affecting the retrieval accuracy, e.g., brine wicking in the snow, or snow layers with different dielectric properties, has to be investigated further in more detail.The several planned Earth-observing satellite missions supporting compact polarimetry (e.g., the RCM operated at C-band) will provide the wide swath cover-age necessary for operational sea ice monitoring.Hence our approach potentially provides a new operational tool for sea ice thickness measurements with a large areal coverage.In this case, the resulting thickness products are also of interest for the development, improvement, and validation of forecast models for the prediction of ice conditions, or of interest for seasonal and climate simulations that consider Arctic and Antarctic ice conditions.

Figure 1 .
Figure1.Variations of the CP ratio as a function of the standard deviation of surface slope σ for different incidence angles and ε = 3.9 + j 0.15.The red line marks the maximum threshold of σ for the validity of our approach.

Figure 2 .
Figure 2. CP ratio as a function of dielectric constants for different σ and incidence angle = 30 • .The results for other incidence angles follow the similar trends.

Figure 3 .
Figure 3. Structure and geometric model of the configuration of sea ice.

Figure 4 .
Figure 4.The simulated sea ice growth process.Blue: sea ice thickness; red: sea ice surface temperature; green: the volume fraction of brine inclusions.

Figure 8 .
Figure 8. Location of the study site in the Labrador Sea, with Pauli RGB (HH + VV for blue, HH − VV for red, and HV for green) decompositions of the RADARSAT-2 images © MDA.The specifications of the SAR data used are given in Table3.

Figure 9 .
Figure 9. Histogram of ice and snow thickness in the Labrador Sea.(a, b) Ice plus snow thickness collected with EMSs in the pack ice (a) and in the fast ice area (b).(c, d) Snow thickness collected with GPR in the pack ice (c) and in the fast ice area (d).The bin widths of ice and snow thickness and snow thickness are 0.1 and 0.02 m respectively.The histograms of the fast ice area are generated from flight tracks of P1, P2, P5, and P7.The histograms of the pack ice area are generated from flight tracks of P3, P4, P6, and P8.These histograms include both level and deformed ice.

Figure 10 .
Figure 10.Regressions relating ice thickness to CP ratio at different incidence angles.The solid lines represent the fits, dashed lines the 90 % confidence intervals.The black, green and red colors are used for the incidence angles of 29, 42 and 49 • , respectively.

Figure 11 .
Figure 11.(a) Relationship between the CP ratio and the observed EM sea thickness.(b) Comparison between the observed and estimated ice thicknesses, and the error bars show the standard deviation with respect to the observation data for every 0.05 m segment of ice thickness.

Table 1 .
Equations and parameters used for the sea ice thermodynamic model.

Table 2 .
Equations and parameters used for the sea ice properties.
* Resolution is nominal.Ground range resolution varies with incidence angle.

Table 4 .
Specifications of helicopter-borne EMS ice thickness data sets.), the air temperature was around −9 to −17 • C on 15-16 March 2011, and snowfall was registered during 2 days in the period 17-19 March with average air temperature around −15 • C. Therefore, a large fraction of the sea ice was covered with snow, which can be clearly seen in aerial photos (not shown).On 19-20 March 2011, the average air temperature was around −8 to −12 • C and the wind speed around 11-15 m s −1