A method for sea ice thickness and concentration analysis based on SAR data and a thermodynamic model

An analysis of ice thickness distribution is a challenge, particularly in a seasonal sea ice zone with a strongly dynamic ice motion field, such as the Gulf of St. Lawrence off Canada. We present a novel automated method for ice concentration and thickness analysis combining modeling of sea ice thermodynamics and detection of ice motion on the basis of space-borne Synthetic Aperture Radar (SAR) data. Thermodynamic evolution of sea ice thickness in the Gulf of St. Lawrence was simulated for two winters, 2002–2003 and 2008–2009. The basin-scale ice thickness was controlled by atmospheric forcing, but the spatial distribution of ice thickness and concentration could not be explained by thermodynamics only. SAR data were applied to detect ice motion and ice surface structure during these two winters. The SAR analysis is based on estimation of ice motion between SAR image pairs and analysis of the local SAR texture statistics. Including SAR data analysis brought a significant added value to the results based on thermodynamics only. Our novel method combining the thermodynamic modeling and SAR yielded results that well match with the distribution of observations based on airborne Electromagnetic Induction (EM) method. Compared to the present operational method of producing ice charts for the Gulf of St. Lawrence, which is based on visual interpretation of SAR data, the new method reveals much more detailed and physically based information on spatial distribution of ice thickness. The algorithms can be run automatically, and the final products can then be used by ice analysts for operational ice service. The method is globally applicable to all seas where SAR data are available.


Introduction
Detailed information on sea ice conditions is crucial for navigation in ice-covered waters. The most important sea ice variables are the ice concentration, ice thickness, and the amount and height of ridges. Hence, national operational ice services try to collect all available information on present ice conditions (ice analysis) and, with the help of prognostic dynamic-thermodynamic models, to forecast future ice conditions. The ice analysis is usually based on one or more of the following sources of information: coastal and shipbased in-situ observations on ice thickness, concentration, and ridges; satellite or airborne remote sensing data on ice thickness and concentration; and numerical model products. So far, however, the observational and model-based information has been treated more or less in parallel, without many activities in interactive combination of the approaches in production of high-resolution ice charts. Assimilation of remote sensing ice measurements into models has been performed widely, but not much using SAR data. Some attempts to assimilate SAR data into sea ice analysis have been reported, e.g. in Kasapoglu (2011). Typically the data assimilated to sea ice analysis have been from microwave radiometers. In our approach we, however, use the SAR data ice concentration as an input for the thermodynamic ice model. The ice analysis is typically illustrated via operational ice charts. The ice information illustrated in the ice charts differs between various national ice services, depending on the observations available and the practices developed. They may also vary depending on the ice analyst; each ice analyst has an individual practice of making ice charts. Some analysts may for example include more details than some others. As an example, the Finnish ice charts show the ice concentration, the mean, maximum and minimum ice thickness, as well as the spatial distribution of different ice types. In addition, the level of ridging is characterized by a number from 1 to 5. In the Canadian Ice Service, the ice chart is in the form of the WMO egg code (Sandven and Johannesen, 2006) presenting the concentration, stage of development, and predominant forms of up to 3-5 ice types for each selected geographical area. The Chinese ice charts for the Bohai Sea illustrate the isolines of ice thickness and ice edge as well as the ice drift velocity field (Luo et al., 2004).
In large sea areas with little navigation, remote sensing is the most powerful tool to collect information needed for an ice analysis, but challenges still remain. These depend on the type of remote sensing data available. In the case of passive microwave data, the detection of ice concentration is a challenge when it exceeds approximately 85 % (Andersen et al., 2007). Retrieval of high ice concentrations by microwave radiometers is degraded by surface effects, sensor smearing and sensor noise. The development of Synthetic Aperture Radar (SAR) methods have yielded a lot of advance, but problems still exist, among others, in the detection of ridges and measurement of ice thickness and above all in the discrimination between thin and thick ice when both are deformed (Mäkynen and Hallikainen, 2004;Shokr, 2009). Also, discrimination of sea ice from open water based on SAR signature is complicated by the ambiguities in the SAR signature of sea ice in comparison to the SAR signature of open water. The SAR signature of open water is a function of the size of the openings and the prevailing wind conditions. However, methods utilizing texture features have yielded better results in SAR-based open water detection, these methods have been described, for example, in Karvonen et al. (2005Karvonen et al. ( , 2010. The challenges in ice analysis also depend on the local conditions. The analysis is more difficult if the ice thickness and concentration are strongly affected by both thermodynamics and dynamics. In seasonally ice-covered seas, the production of ice (basin-scale temporal evolution of ice thickness) is typically dominated by thermodynamic processes; the air-ice-ocean thermal interaction is crucial to determine the columnar and granular ice growth (Bitz et al., 2001) (columnar growth is due to freezing of sea water, whereas granular growth is due to refreezing of melted snow). The spatial distributions of ice thickness and concentration are, however, often dominated by ice dynamics (Fichefet and Maqueda, 1997). A preliminary investigation of ice thickness analysis based on SAR data and a thermodynamic ice model has been made for the Baltic Sea (Karvonen et al., 2007(Karvonen et al., , 2008. Also, the RADARSAT geophysical processor system (RGPS) (Kwok, 1998) has utilized the thermal growth to estimate the ice thickness distributions of Lagrangian cell ice age classes. The RGPS ice age classes are computed based on the ice drift between to adjacent SAR images over the same area. The method used in RGPS re-quires reliable tracking of each Lagrangian cell from image to image in time. In practice, this kind of tracking using available SAR data, e.g. for Gulf of St. Lawrence, is not possible because of too sparse temporal coverage. In addition, the local SAR backscattering dynamics is attenuated by occasional above-zero temperatures, and for this reason the SAR images typically contain less textured areas for which meaningful correlations between pairs of SAR images can be obtained and reliable ice drift estimates be given. Due to these restrictions, our approach for ice parameter estimations is different from the RGPS.
In this study, a new method is presented for sea ice thickness and concentration analysis based on SAR data and a thermodynamic model. Here, we address the seasonal sea ice cover in the Gulf of St. Lawrence (GSL) off Canada, but the methodology is applicable for all ice-covered seas with weather forcing and SAR data available. In the GSL the ice cover is more dynamic than in the Baltic Sea. Hence, the SAR-based method for the Baltic Sea (Karvonen et al., 2007(Karvonen et al., , 2008 cannot be directly applied for the GSL, because ice motion detection from SAR data is also needed. First we present a new method to derive sea ice concentration and drift vectors on the basis of SAR data. In locations where the SAR data indicate presence of sea ice, a thermodynamic sea ice and snow model HIGHTSI is used to calculate the ice growth, forced by operational model products of the European Centre for Medium-Range Weather Forecasts (ECMWF). The thermodynamically based ice thickness is used as background information for the remote sensing algorithm. For the first time, we combine the SAR-based information and thermodynamic modeling of sea ice and its snow cover. The resulting sea ice thickness and concentration analyses are finally compared against airborne Electromagnetic Induction (EM) measurements and the present operational products of the Canadian Ice Service (CIS).

Study area, ice seasons, and data
The GSL is a seasonally ice-covered gulf at latitudes as low as 45.6-52.0 • N (Fig. 1). Winter navigation in the GSL is active, as it is the pathway to St. Lawrence River providing access to Quebec City, Montreal, and the North American Great Lakes. Due to the high economic importance of navigation, there is a strong need to provide accurate sea ice analysis for the GSL. The dynamic sea ice cover in the GSL represents a challenging environment for ice analysis.
Two winter seasons, 2002-2003 and 2008-2009 were selected for this study. The first one was a severe and the second one a milder winter. In this study, we defined the winter season starting from the day when the operational ice chart service began and ending when the ice chart production ended.   The WMO egg code given for each polygon (area). C t is the total concentration; C a , C b , and C c are the fractional concentrations; S a , S b , and S c are the fractional stages of development corresponding to the fractional concentrations; and F a , F b , and F c are the fractional predominant forms of ice. S o and S d are the stages of development of the possible remaining ice types within the the area. All three fractional values are not necessarily assigned for each polygon if there are less ice types present according to the ice analyst interpretation. agery and ice reconnaissance flights. The ice conditions are qualitatively illustrated by the standard WMO egg code over polygons defined by the CIS ice analysts in various parts of the GSL. A diagram of the WMO egg code is given in Fig. 2. The CIS ice charts are drawn daily and they rely on EO data and ship and aircraft observations. RADARSAT SAR data with a wide coverage are their major data source. For making ice charts the SAR data is first analyzed and then a forecast (nowcast), based on other sources of information (including other EO data, ice reports from ships), is made. The CIS ice charts have not been systematically verified objectively, but they represent the best and most accurate ice information available in GSL. The daily CIS ice charts are issued at 18:00 UTC. They at least give a qualitative description of the prevailing ice conditions. Additionally, we have also tried to use the ice charts for quantitative evaluation by using the total ice concentration and the ice thickness ranges given by the stage of development for each CIS ice chart segment (defined by the bounding polygons drawn by the ice analysts).
The quantitative mapping of the egg-code stage of development to ice thickness was made following the criteria given in Table 1, and the conversion of the egg-code total concentration (C t ) was made according to Table 2. The egg code gives for each polygon the total concentration, fractional concentrations, fractional stages of development (related to the ice age), and fractional predominant forms of 1st-year ice 50 7 Thin 1st-year ice 70 8 1st-year thin -1st stage 50 9 1st-year thin -2nd stage 70 1.
Thick 1st-year ice 140 Fast ice 120 ice. The segment total concentration value C t was used as the concentration value in the comparisons. The ice thickness was derived from the stage of development, and we have here used the lower boundary of the ice thickness ranges given by the stage of development. For each ice chart segment, the stage of development corresponding to the most frequent fractional concentration is used in the mapping of ice thickness. In the case of multiple ice types with different ice thickness ranges present within one ice chart segment, the used ice thickness mapping may to some extent deviate from the mean of the ice thickness lower boundaries computed based on the ice chart fractions of the different ice types. The quantified daily ice charts were interpolated to the HIGHTSI model domain (3111 grid cells, Fig. 1) and the ice thickness and concentration from each grid cell were extracted on each day. Due to the dynamic sea ice cover in the GSL, the spatial distribution of ice thickness at a given time step is very complex and difficult to assess. The seasonal ice evolution, however, can be illustrated by the integrated basinscale mean daily ice thickness H a and daily total ice cover area H A : where H i is the ice thickness in the grid cell i, N is the total number of grid points in the model domain, t is the time (day), C i is the ice concentration, and G i is the grid cell i area. Figure 3 shows the quantified CIS-based daily time series of ice area as well as maximum and minimum basin-scale ice thicknesses for winter seasons -2009. Ice winter 2002 was severe. The maximum values of ice area and mean basin-scale ice thickness occurred in late March. In the early freezing season, the ice cover developed faster with respect to area than thickness; from 1 January to 1 February, the ice area increased by 11 times whereas the ice thickness only increased by 5 times. During the melting season, however, ice thickness and area decayed with the same rate. Ice winter 2008-2009 was milder; in early winter the ice area increased faster than ice thickness. After the maximum ice cover was reached, the ice area and thickness varied with the same rate until the final ice melt. Compared with winter 2008-2009, the maximum ice area in winter 2002-2003 was almost 3.5 times as large and the basin-scale sea ice mean thickness was 3.1 times as thick. The maximum ice area and mean ice thickness appeared in mid-March to be in line with the ice climatology in the GSL.
The electromagnetic induction (EM) (Prinsenberg et al., 2002) ice thickness measurements have been carried out over the GSL during several seasons. The measurements for the two particular seasons studied in this experiment were used as a reference ice thickness data set. EM measurements are considered as one of the most accurate methods to measure sea ice thickness over a large area that cannot be covered by drilling Haas et al., 2006). In both the winters 2002-2003 and 2008-2009, the EM measurements were made within a relatively short period of 1-2 weeks from late February to early March. The measurements were made close to the Prince Edward Island, and were mainly headed to the north of the island (Fig. 4). The sampling rate of the EM measurements was 3-4 m and the resolution (footprint) was around 20-30 m. The EM flight lines used in this study with their acquisition dates shown by color codes are shown in Fig. 4.

HIGHTSI model
A high-resolution thermodynamic snow and ice model (HIGHTSI) is used in this study to simulate thermodynamic ice growth and melt. HIGHTSI was initially developed to study snow and ice thermodynamics in seasonally ice-covered seas (Launiainen and Cheng, 1998;Vihma et al., 2002;Cheng et al., 2003Cheng et al., , 2006. The model has been further developed to investigate snow and sea ice thermodynamics in the Arctic (Cheng et al., 2008a,b) and also lake ice Semmler et al., 2012).
Ice and snow thickness, heat conductivity and temperature are simulated, solving the heat conduction equations for multiple ice and snow layers. In addition, special attention is paid to the parametrization of the air-ice fluxes and the solar radiation penetrating into the snow and ice. The turbulent surface fluxes are parametrized taking the thermal stratification into account. The penetration of solar radiation into The Cryosphere, 6, 1507-1526, 2012 www.the-cryosphere.net/6/1507/2012/     the snow and ice depends on the cloud cover, albedo, color of ice (blue or white), and optical properties of snow and ice. The shortwave radiation penetrating through the surface layer is parametrized, making the model capable of quantitatively calculating sub-surface melting. Short-and long-wave radiative fluxes can either be parametrized or prescribed based on results of numerical weather prediction (NWP) models. The surface temperature is then solved from a detailed surface heat/mass balance equation, which is defined as the upper boundary condition and also used to determine whether surface melting occurs. At the lower boundary, the ice growth/melt is calculated on the basis of the difference between the ice-water heat flux and the conductive heat flux in the ice. The ice-water heat flux (refers to oceanic heat flux in this study) is assumed as a constant. Model experiments in the Arctic Ocean have suggested that albedo schemes with sufficient complexity can more realistically reproduce basin-scale ice distributions (Liu et al., 2007). We applied a sophisticated albedo parametrization according to temperature, snow and ice thickness, solar zenith angle and atmospheric properties (Briegleb et al., 2004). The thermal conductivity of sea ice is parametrized according to Pringle et al. (2007), and snow density and heat conductivity are calculated according to Semmler et al. (2012).
In seasonally ice-covered seas, snow plays an important role in sea ice growth. On one hand, snow generates a strong insulation and prevents sea ice growth. On the other hand, snow may transform to granular ice due to refreezing of ocean flooding and snow melt, enhancing ice growth. Snow-to-ice transformation through re-freezing of flooded or melted snow are considered in the model. A reliable calculation of snow depth heavily relies on precipitation, which is given as external forcing.

Forcing data and studied winter seasons
The atmospheric forcing for HIGHTSI was based on the operational analysis and short-term forecasts of the ECMWF. The reason for applying ECMWF products is their global availability. This is essential as we want to develop an ice analysis method applicable in all ice-covered seas. The analyses were available with 6-h intervals. The 2 m air temperature and humidity as well as the 10 m wind speed were applied to parametrize the turbulent fluxes of sensible and latent heat. The downward components of the solar shortwave and thermal long-wave radiation were based on 12-h operational forecasts. To avoid spin-up problems in precipitation (Tietäväinen and Vihma, 2008), snow fall was based on 24-h forecasts. We selected the ECMWF products as forcing for HIGHTSI because previous studies have demonstrated that the accuracy of ECMWF products is high compared with other re-analysis data, such as NCEP/NCAR, at high latitudes (Bromwich et al., 2009;Cheng et al., 2008a;Jakobson et al., 2012). To get an idea of the accuracy of the ECMWF products in the GSL, the wind speed (V a ) and air temperature (T a ) at selected locations were compared with results of the Environment Canada regional operational forecast model (in 15 km resolution). The results of the ECMWF and Environment Canada models agreed very well. The correlation coefficient between ECMWF and CIS modeled basin-scale wind speed and air temperature were 0.97 and 0.98, respectively, for both winter seasons. The temporal standard deviation of ECMWF and CIS wind speed was the same within 0. . The wind direction is not investigated as HIGHTSI does not include ice dynamics. As a summary, either of the atmospheric models could have been used to provide forcing for HIGHTSI. We selected the ECMWF model because of our objective to develop a globally applicable method for sea ice analysis.
To carry out HIGHTSI simulations over the entire GSL, atmospheric forcing is needed over a high-resolution grid. The ECMWF analysis and forecasts were available with a horizontal resolution of 0.358 • for winter 2002-2003 and 0.225 • for winter 2008-2009, whereas HIGHTSI was applied in a resolution of 0.225 • in longitude and 0.1125 • in latitude. Hence, the ECMWF forcing fields were linearly interpolated to the HIGHTSI grid.

Model simulations
HIGHTSI simulations were carried out at each of the 3111 grid cells. Each model run starts 1 December and lasts until 31 May. The simulations are initialized with a thin ice layer (0.01 m) at each grid cell. The ice thickness is calculated as a mean over the ice-covered area in the grid cell. During the simulation, if the external weather data do not favor ice growth, HIGHTSI maintains the ice thickness of the previous time step. Off the coastal land-fast ice zone, ice drifts significantly in the GSL. The impact of ice drift is taken into account by incorporating sea ice concentration (A) derived from SAR data into the HIGHTSI model (see Sect. 4). When a grid cell is at least partly covered by ice (A > 10 %), the ice growth is calculated by applying the atmospheric forcing at that particular grid cell. If the ice concentration at a certain grid cell is reduced below 10 %, the ice thickness will remain at the value of the previous time step. The calculation of ice growth is resumed once the ice concentration is again larger than 10 %. This approach implicitly excludes new ice formation in open leads. In winter season, however, open leads or polynyas cover a small fraction of the GSL. In the GSL, the fresh water runoff drives circulation in the Gulf, and the upper layer ocean tends to be less saline and cold in early winter (approximately as in the estuary). As a first order estimation, we assume the oceanic heat flux is 2 W m −2 in early winter. Far from the cost, the oceanic heat flux is parametrized as a function of the SAR-based ice concentration. It is assumed to increase from 2 W m −2 to 6 W m −2 when ice concentration decreases from 100 to 10 %. This assumption is related to the low latitudes of 46-52 • N, where increasing amounts of open water are usually associated with increasing solar energy absorbed by the ocean, consequently increasing the oceanic heat flux. Our model run excluded the new ice formation in open leads. Figure 5 shows the modeled basin-scale average ice thickness compared with the quantified mean ice thickness from CIS ice charts. From the operational point of view, the maximum ice thickness is more important for navigation. Hence, we consider it as a good result, if the overall calculated and observed maximum ice thickness are in good agreement. For the cold winter season 2002-2003, the modeled ice thickness showed a reasonable agreement with CIS ice charts, except in March when the ice thickness was underestimated by around 35 %. For the milder winter 2008-2009, the HIGHTSI model run showed a better agreement with CIS ice charts. For a cold winter, the large ice growth may be associated with frequent high pressure conditions, which are associated with less precipitation. On the contrary, a milder winter is usually associated with low pressure conditions, which means frequent cases of advection of warm, moist marine air masses, bringing more precipitation (snowfall). Overall, the amounts of snow and winter air temperatures are the most important factors for the thermodynamic ice growth. The ECMWF accu-  equivalent, which was 45 % more than in winter 2002-2003 (231 mm). The role of snow may generate uncertainties due to potential errors in precipitation based on a NWP model. We emphasize that the HIGHTSI-based ice thickness only accounts for the thermodynamic growth and melt. The dynamic ice growth such as rafting and ridging are taken into account by ice motion detection using SAR data. The HIGH-TSI model runs are used as background information for the remote sensing algorithm described in the following section.

SAR remote sensing
In the experiment we used HH-polarized (SAR antenna transmitting horizontally polarized and receiving horizontally polarized signals) SAR data from RADARSAT-1 (season 2002RADARSAT-1 (season -2003 and RADARSAT-2 (2008RADARSAT-2 ( -2009). RADARSAT-2 has an option to acquire dual-polarized images in the wide swath imaging mode, but we only had dualpolarized data (HH/HV polarization combination, HV denotes transmitting horizontally polarized and receiving vertically polarized signals) over a short period (a few weeks in late February and early March 2009) at our disposal, and in this experiment we only utilized HH-polarized data, also including the HH-channel data of the 2009 HH/HV data. Three SAR algorithms were used in this experiment. The first SAR algorithm estimates the ice concentration, based on SAR texture features. This concentration is used as an input for the HIGHTSI model and also for the ice thickness estimation. The SAR ice drift algorithm estimates the ice drift between two co-registered SAR images, in this case two adjacent daily SAR mosaics, acquired over the same area at different time instants. The ice thickness estimation method presented here is a novel method combining the remote sensing and model data. The SAR ice thickness algorithm estimates the ice thickness using the modeled ice thickness distribution, kinematic features based on the SAR ice drift, SAR features, and SAR ice concentration as inputs.

SAR data preprocessing
The SAR images are first calibrated to get the logarithmic backscattering coefficient values, which are presented in decibels. The backscattering coefficient, σ 0 , describes the properties of the target area producing the backscatter of each SAR pixel. The σ 0 values are then linearly sampled to eight bits per pixel (8 bpp) images; the scaling interval is from −35 dB to 0 dB; backscattering values higher than 0 dB are saturated to pixel value 255; and values lower than −35 dB are saturated to pixel value one. The general calibration equations are K is a SAR calibration coefficient, and typically given in the SAR metadata; A is the SAR amplitude value; and I = A 2 is the SAR intensity value. α is the SAR incidence angle, which increases from the near range to far range, and typically varies about in the range of 20-50 • . The 8 bpp SAR images are then rectified to Lambert Conformal Conic (LCC) projection. After rectification, a land masking to mask off all the land areas is applied. The land masking was based on the GSHHS (A Global Self-consistent Hierarchical Highresolution Shoreline Database) shoreline (Wessel and Smith, 1996). After calibration, rectification and land masking, we also apply an incidence angle correction (Mäkynen et al., 2002), which is based on an empirical linear relationship between mean sea ice backscattering and incidence angle value. The same incidence angle correction defined for the Baltic Sea has performed well for Arctic SAR data also, and according to visual inspection of some test mosaics over GSL, this linear correction also works fine for the GSL SAR data. The incidence angle correction was made to a reference incidence angle of 30 degrees, which is approximately in the middle of the RADARSAT ScanSAR mode incidence angle range. The incidence angle correction is important in building mosaics www.the-cryosphere.net/6/1507/2012/ The Cryosphere, 6, 1507-1526, 2012 used in classification, this correction smooths the SAR image boundaries within the SAR mosaic. After the incidence angle correction, we build daily mosaics of the SAR data in a constant grid with a 500 m resolution in the LCC projection. The daily mosaics are built such that they always contain the most recent SAR values, i.e. the newer SAR images are laid over the older data. In the areas without new data, the older data remains, and these still represent the most recent SAR measurements over these areas. The SAR mosaics are used as inputs for the ice concentration and ice thickness estimation algorithms. All the SARrelated computations are made using the mosaics instead of single SAR images. In practice, the entire study area is not covered daily, and the SAR value at a location can be 0-3 days old -in some cases for the 2002-2003 data set (with less SAR data) even five days old in some areas of the mosaic (worst case). In general, the SAR mosaic update rate for a certain location was 0-3 days for the 2002-2003 data set, and 0-2 days for the 2008-2009 data set, except for a few worse cases. Zero days indicates that there exist SAR data acquired twice over the same area within one day, e.g. SAR data acquired in the morning and evening of a same day, or in the evening of a day and in the next morning. An example of a SAR mosaic (5 March 2009) is shown in Fig. 6. After the incidence angle correction, the SAR image boundaries in the ice areas do not affect the segmentation and classification significantly, as they would without an incidence angle correction. In the open water areas, the incidence angle correction defined for sea ice does not work, but these areas can rather reliably be located by our ice concentration algorithm.
In the next phase the SAR mosaics are segmented. Here we use a simple k-means clustering (MacQueen, 1967) for the segmentation with the parameter value k = 10. The small segments are then removed iteratively, starting from the smallest ones, merging them to the larger adjacent segments with the closest mean (incidence angle corrected) backscatter value. The size (in pixels) for the segments to be merged is increased iteratively starting from one and ending up to the given upper boundary. We have used the upper boundary of 100 pixels in the segment merging, i.e. the segments smaller than 100 pixels are merged to the neighboring segments. After the merging the lower boundary of the segment size is 100 pixels, and there is no defined upper boundary for the segment size. The segmentation is an important stage of the SAR processing, because the ice concentration and thickness are estimated for each SAR segment, instead of single pixels. We use a lower bound for the segment size because it is necessary to have a large enough number of points to compute reliable statistics. Computing features over segments also maintains the natural ice segment boundaries, unlike computing over regular data windows of a fixed size.

Ice motion detection
The ice motion estimation algorithm is based on phase correlation computed in two resolutions. First the larger-scale ice motion is detected in the coarse resolution, and then it is refined in the fine resolution. A multi-resolution phasecorrelation algorithm for sea ice drift estimation was also presented in Thomas et al. (2004). In our approach, the phase correlation computation is performed only for areas containing edges in the SAR images, because these are the areas where reasonable cross-correlations can be found. In practice, these are locations with more edge points, detected by an edge detecting algorithm, in their neighborhood within a given radius, relative to the area of the neighborhood, than a given threshold value. In smooth or random areas, the crosscorrelations are too low for a reliable ice drift detection.
The edges in SAR images are here detected using the Canny edge detection algorithm (Canny, 1986), but in practice most of the edge detection algorithms produce similar results. Typically the edge detection is based on pixel value gradient magnitude in the images. After edge detection, we perform a filtering of the edges, filtering out small edge segments. An edge segment is here defined as a set of connected edge pixels, in the sense of the pixel's 8-pixel neighborhood. All the edge segments smaller than a given size (an integer number of pixels) are removed from the edge image. A suitable size threshold for a full-resolution SAR edge image is five. The small edge segments are typically due to speckle and do not describe any actual SAR features. An example of the edge detection and filtering is shown in Fig. 7.
The ice motion is then determined for sampled data windows from the two images using the phase correlation in both resolutions. The window size is W × W pixels; we use W = 16. To compute the phase correlation, 2-D FFT is applied to the data windows sampled from the two adjacent daily SAR mosaic images at the same location. Then FFT coefficients of the two image windows are normalized by their magnitudes, and then the FFT coefficients of the two image windows are multiplied and the inverse 2-D FFT is applied, i.e. the phase correlation (PC) array is computed from the the normalized cross power spectrum. The best matching displacement in a Cartesian (x, y) coordinate system is defined by the maximum PC: (4) A and B are the Fourier transforms of the two sampled data windows from the two temporally adjacent daily SAR mosaic images; k and l are the coordinates in the Cartesian coordinates in the FFT transform domain; and FFT −1 denotes the inverse Fourier transform, using the FFT algorithm. To get the low-resolution image, the two SAR mosaic images are first down-sampled to the given low resolution. The downsampling ratio R S is a power of two (because the FFT algorithm requires this); we use R S = 8. The two low resolution images are generated by successively applying an optimal half-band low-pass FIR filter designed for multi-resolution image processing (Biazzi et al., 1998) and re-sampling the image (down-sampling by two in both x-and y-directions). This filtering preserves the edges in the coarser resolution. Multiple coarse-resolution ice drift candidates, corresponding to several detected coarse-resolution phase correlation maxima, are typically found for each data window pair. The coarse-resolution motion candidates, corresponding to the PC maxima, are refined in the fine resolution and one of them is selected to represent the final drift estimate in the fine resolution. The details of the ice drift estimation algorithm are presented in Karvonen (2012a).
The ice drift is computed using a window size of W × W pixels in the 500 m mosaic resolution (see Sect. 4.1), and the motion is given in a grid size of W/2 pixels in both x-and y-directions. We have used the value W = 16, i.e. the grid resolution is 4 km and the window size at the fine resolution level is 8 km × 8 km. The coarse resolution is 500 m down-sampled by the factor R s , in this case 4 km, and the corresponding coarse resolution window size is thus 64 km × 64 km.
Finally, a vector median filtering (Astola et al., 1990) is performed in the high-resolution ice motion grid to obtain the final motion estimate.
Divergence (D) and shear (S) of ice motion are calculated from the ice drift vector field approximating the spatial derivatives by finite differences in the grid. The total deformation (Lindsay et al., 2002) can then be computed as Another feature derived from the ice drift is the ice drift ratio, denoted by R M at each grid position (r, c): In the equation, f i (r, c) is the ice motion vector at location (r, c) at the moment i, i.e. f i describes the drift between two adjacent SAR images. The sums at each location (r, c) are computed over the motion vectors from SAR image pairs during a certain time period, e.g. one week or two weeks. The indexes t 0 and t n refer to start and end times of the summing, i.e. the values are accumulated over a certain period starting at t 0 and ending at t n (for daily mosaics the time is counted in days). One is subtracted just to change the lower bound of the values of R M to zero instead of one (the numerator is always larger than or equal to the denominator). R M is typically higher in the areas of thicker ice and lower in the areas of thinner ice. In the areas where this ratio is higher, the motion has not been in a certain direction, but the motion direction has been oscillating; and in the areas of low ratio the motion has been more directional, either through or out of the area. We here call D T and R M kinematic features. The kinematic features are used in the ice thickness estimation algorithm.

SAR-based ice concentration
For the ice concentration detection, we have used two SAR texture features. The first is the relative amount of edge and corner points with respect to the segment size for each SAR segment. The edges and corners are detected based on a variant of the Local Binary Patterns (LBP) texture measure (Ojala et al., 1996). In this version the LBP patterns are formed such that a threshold (T ) for the absolute difference between the middle point and the reference points is applied; www.the-cryosphere.net/6/1507/2012/ The Cryosphere, 6, 1507-1526, 2012 here we use T = 10. The reference points are the eight points with the angular step of /4 around the middle point within a given distance R; we use R = 2. After computing LBP for each pixel location, the minimum of all the cyclic shifts of the eight-bit LBP number is selected to be the final LBP for the pixel; this operation makes the values rotationally invariant. These LBP values can then be interpreted as edges (LBP value 15), blunt corners (31) and sharp corners (63). The relative number of edge points within a segment is then the number of the points with the above-mentioned three LBP values (denoted by N e ) divided by the area of the segment A, computed in pixels: The other texture feature is based on the three gradient features, computed in fixed-size, round-shaped (here radius R = 5) data windows over the SAR image at each image pixel location. The three gradient features basically describe the magnitude and directional coherence of the local image pixel gradient. The gradient texture features are described in detail in (Karvonen, 2012b). Here, we used a combination M g of the three gradient features: This combined gradient feature M g is high in the areas where the SAR signal gradient is high and oriented, i.e. in the areas where the SAR signal is not smooth and there is some kind of non-random structure present. In the concentration algorithm we use segment means of the combined gradient texture feature. An example of the two features used is shown in Fig. 8. As an input for the concentration estimation, we use the maximum of the two texture features, i.e. F = max(F a , M g ). The concentration C is estimated with an experimental relationship: Here, we have used T 1 = 10 and T 2 = 30. The thresholds do not have a physical unit; they are related to the SAR texture features. The factor a c is defined such that C is continuous at T 1 and T 2 . The threshold T 1 and T 2 were experimentally defined based on visual comparison of six 2003 SAR mosaics and the corresponding CIS ice charts. The fraction of open water (concentration 0 %) can be increased by increasing the value of T 1 , and the amount of 100 % concentration ice can be increased by decreasing T 2 . T 1 is smaller than or equal to T 2 . In the case of equal values for the thresholds T 1 and T 2 , the algorithm produces only either open water or ice with 100 % concentration.
Open water areas can then simply be found by thresholding the concentration. Here, we use 50 % as the threshold. We use this low threshold to make sure all the ice-covered fields are classified as ice, because this discrimination between open water and ice is used as an input for the HIGH-TSI model, and we want to guarantee inclusion of all the ice areas in the HIGHTSI computation.
On the basis of comparison with a visual interpretation of CIS ice charts, the ice concentration estimation seems to work rather well. In general, compared to the CIS ice charts, the concentration was overestimated. We compared the algorithm concentrations to CIS ice chart concentrations, extracted as described in Sect. 2, during the period from the beginning of January to the end of March for both the test seasons. The signed mean estimation error, i.e. mean of the CIS ice chart concentration, subtracted from the estimated concentration for this period, was 9.7 percent units for the 2003 data, and 4.4. percent units for the 2009 data set over the same period of the year. During the freeze-up and melting seasons, the overestimation was smaller because the apparent larger open water areas are correctly mapped to the concentration of 0 % by the algorithm. During the melting season the algorithm performance is degraded by the wet ice cover and melt ponds, causing underestimation of concentration in these areas. Wet snow cover attenuates the backscattering, and areas of smooth wet snow cover can have SAR signatures similar to calm water. Also, large melt ponds have a similar SAR signature as open water. However, according to our experience, the results were not significantly degraded by warm weather conditions, and we were still able to give reasonable ice concentration estimates. The main problem for the algorithm seemed to be the overestimation of non-zero concentrations during the mid-winter period.
The reliability of the concentration can be improved by using multiple SAR images over the same area with a relatively short temporal difference. Then, for example, the temporal medians of the features can be used instead of the The Cryosphere, 6, 1507-1526, 2012 www.the-cryosphere.net/6/1507/2012/ features from a single SAR measurement, making the estimation more robust, especially over open water areas with occasional short (wavelengths in the magnitude of the SAR wavelength) sea wave patterns.

SAR-and model-based ice thickness
The correlation between SAR backscattering and different ice types is not always very high (Mäkynen and Hallikainen, 2004;Shokr, 2009); for this reason it is necessary to utilize ice motion and SAR texture features in the ice thickness estimation. Our ice thickness algorithm utilizes the SAR features, ice thickness modeled by HIGHTSI, kinematic ice features (D T , R M ) derived from the SAR-based ice drift, and ice concentration from the SAR ice concentration algorithm. A general flow chart of the ice thickness algorithm is presented in Fig. 9. The initial steps of the algorithm are forming and updating the image mosaic. In this process each new SAR image is drawn over the existing mosaic (which in the beginning of each season is empty). The mosaic is used as a basis for the daily ice concentration and thickness estimates. The inputs in the diagram are shown in cyan color, and the final output is indicated by red color. The NWP data is received from the ECMWF model; the SAR data is the RADARSAT data from CIS; and the linear mapping coefficients (a and b), used in the final ice thickness estimation phase, are defined based on a training data set, consisting of the 2003 EM measurements and SAR mosaic data of the same days with the EM measurements.
To get reasonable estimates of the ice thickness, training data over a representative ice season is required to define the algorithm parameters. In the next phase, the HIGHTSI ice thickness is remapped based on several segment-wise image features, computed over the SAR mosaics of a period of two weeks before the mosaic of each day. The ice thickness is mapped such that the cumulative feature distribution is mapped to correspond to the HIGHTSI ice thickness distribution computed over the whole study domain. After this mapping computation the SAR feature values in the image grid are mapped to the corresponding HIGHTSI ice thickness values according to the cumulative distribution matching, i.e. each SAR feature value has a corresponding unique ice thickness value based on the distribution mapping. These statistics represent the cumulative ice dynamics during each two-week period. The features computed from the ice motion are the drift ratio R M over the period, total deformation D T , mean SAR pixel value µ (with incidence angle correction applied) over the period, and Fig. 9. A general block diagram of the ice thickness estimation scheme one step to another are also shown (with the same names as in the equat 46 Fig. 9. A general block diagram of the ice thickness estimation scheme. The parameters delivered from one step to another are also shown (with the same names as in the equations). mean relative amount of edges E, with respect to the segment area, over the period. The geometric mean of these features is then used to form a new feature F 1 for each pixel with row and column coordinates (r, c): The feature F 1 is used as a segment-wise means. It describes the ice accumulation during the preceding two-week period updated by a SAR texture feature describing ice deformation within the segments. The cumulative distribution of the HIGHTSI ice thickness, computed over the whole study area, is mapped to correspond to the distribution of the feature F 1 over the study area such that the Kolmogorov-Smirnov distance (DeGroot, 1991) between the HIGHTSI ice thickness distribution and the mapped feature distribution, both computed over the whole study domain, is minimized. After this mapping of distributions, the feature F 1 values at each pixel are mapped to ice thickness values according to the distribution mapping; each feature F 1 value is uniquely mapped to an ice thickness value. The remapped ice thickness grid is then used as an input for the daily ice thickness estimation. This step represents remapping of the thermodynamic growth based on medium-scale temporal dynamics, i.e. based on currents and wind-driven drift during the two-week period, combined with the SAR-based information on the local deformation. This gives a segment-scale input for the www.the-cryosphere.net/6/1507/2012/ The Cryosphere, 6, 1507-1526, 2012 final local ice thickness estimation. This is our background ice thickness grid. An example of this mapping is shown in Fig. 10. This remapped HIGHTSI ice thickness at (r, c) is here denoted by H h (r, c). Because winds, e.g. winter storms, can cause rather rapid changes in the GSL ice cover (Prinsenberg et al., 2006), the segment-scale background ice thickness field still needs to be refined based on the most recent available local SAR backscattering and texture information in the fine scale (here in 500m mosaic pixel scale). In this final phase of the ice thickness estimation, the segment-scale remapped HIGHTSI ice thickness is used as an input, and it is refined by using the pixel-wise SAR features to get the final local ice thickness estimates. The basic assumption here is, like for the Baltic Sea ice thickness charts (Karvonen et al., 2003), that the local SAR backscattering and its texture are functions of the ice thickness. In general, it is assumed here that the backscattering and the texture at a local scale are stronger for older and thicker ice, which has gone through more deformation events, conditioned by the background ice thickness based on the remapped feature and HIGHTSI distributions. In this final phase the coefficient of variation C v , i.e. the ratio of the local (here the segment-wise) standard deviation and the local mean (given in percents), is used as an additional texture feature: where µ(r, c) and σ (r, c) are the local mean and standard deviation, respectively, computed for the segment in which the pixel at (r, c) belongs. This feature is also quite robust against variations due to incidence angle. The incidence angle corrected SAR pixel value mosaic P (r, c) is modulated by C v (r, c) (geometric mean), and this value F 2 (r, c) is used in the fine-scale ice thickness estimation: F 2 (r, c) = (C v (r, c)P (r, c)) 1/2 .
The remapped HIGHTSI ice thickness values H h (r, c) are first modulated by the value of F 2 , and the ice thickness H (r, c) at (r, c) is estimated as a linear function of the geometric mean of this modulated term: We use the values a = 5.52 and b−112.85 for the linear mapping parameters. These parameters were defined based on the 2003 EM data (training data set) and the (H h F 2 ) 1/2 distributions, using data derived from the same day SAR data with the EM measurements. The distributions were first estimated by smoothing the distributions by a Gaussian kernel with a standard deviation of 10.0. Then the means and standard deviations for the two distributions (denoted by µ 1 , µ 2 , σ 1 , and σ 2 ) were computed. Assuming, for simplicity, that the distributions are Gaussian, we get the mapping coefficients (a and b) between the two distributions: 6 Comparisons against CIS ice charts and EM data

CIS ice charts
The total ice concentrations C t extracted from CIS ice chart and SAR algorithm for the fifth day of the winter months (January-April) for both ice seasons used in this study are shown in Figs. 11 and 12. The overall agreements are reasonably good, although the CIS charts give smaller ice concentrations in both winter seasons, particularly in February and April. In each February and April, the difference for the cold winter (2002)(2003) is slightly higher compared with results for the milder winter (2008)(2009). However, because the distinction between ice and open water is based on the ice concentration algorithm output, and this result is used as an input for HIGHTSI, we rather slightly overestimate the concentration to include all the ice-covered areas in the HIGH-TSI computation than leave ice-covered areas out of the computation. The difference in the mid-winter (January-March) was from 4.4 (2009) to 9.6 percentage points, the overall L1 (absolute mean error) was 15.9 (2003) and 13.9 percentage points. Figures 13-14 show the daily ice thickness estimation from CIS ice charts (lower bounds of the ice stage of development), the HIGHTSI model alone and the SAR algorithm results for the fifth day of each winter months (January-April). Visual inspection reveals that differences between CIS ice chart and SAR-based ice thickness do exist when making point to point comparisons, and they may also be quite large in certain areas. The CIS ice chart attributes are given for large areas (polygons in the ice charts), and they can merely be used to indicate the ice thickness within the area relative to the ice thickness of other ice chart polygons, because the ice thickness ranges of the degree of development categories are wide and usually overlap each other. This crude analysis also contributed to differences between CISand SAR-based ice thickness. The ice drift was estimated in a resolution of 8 km (corresponding to the window size of the ice motion estimation algorithm, multiplied by the SAR mosaic resolution, i.e. 16 × 500 m); in this scale the motion in narrow straits and river outflows cannot be detected, and the HIGHTSI remapping in these areas fails. This could probably be improved by estimating the motion in a higher resolution (we have used the 500 m mosaic resolution). On the other hand, there are similarities between CIS-and SARbased ice thickness, particularly in terms of ice thickness distributions. The distributions for the area where the annual EM measurements have been made match each other rather well. For the areas with larger deviations between the CIS ice charts and the thickness estimates, the differences can be due to failure of our algorithm to capture all the ice motion, misinterpretation of the SAR features by our algorithm, or by the coarse presentation in the CIS ice charts. The ice thickness distributions from HIGHTSI represent the thermodynamic ice growth, and considerably differ from distributions based on the CIS ice charts. Because HIGHTSI literally produces the basin-scale thermodynamic ice growth, the effect of ice dynamics such as ice deformation, rafting, and ridging are not included. On the other hand, the basinscale sea ice production agreed reasonable well compared with CIS data (Fig. 5).
We want to emphasize that there is still room for further development of our approach. Currently, we cannot identify whether the differences between CIS-and SAR-based ice thicknesses are due to inaccuracies of the ice thickness extracted from the CIS ice chart or due to estimation inaccuracy of our algorithm. At least, one limitation of our algorithm is In general, we can say that the CIS ice charts give a useful reference for comparing the ice concentration; the ice thickness is only a rough first order estimate, but it can at least be used in qualitative visual comparison. Ice concentration has a different nature than ice thickness, because it is a quantity relative to the segment area, and reasonable values to large segments can always be given. However, because no other ice thickness measurements covering large spatial and temporal scales exist, we have also made quantitative comparisons between CIS and SAR algorithm thicknesses. The differences in ice thickness between the SAR products and CIS ice charts can be explained by multiple factors: (a) the ice thickness in CIS ice charts is only implicitly expressed by the typical ice thickness range for the specific stage of development, (b) the resolution of the ice charts is significantly lower than that of the SAR products, and (c) the ice properties in ice charts are only described for large polygons whose sizes are at least tens of square kilometers or even more. We can see that in some areas the ice thickness given by our algorithm is quite different from the CIS ice charts. These cases are partly due to the above-mentioned reasons, and also partly because our algorithm fails in the estimation in some areas (see above for the limited training of the algorithm; the CIS ice charts could not be used as a training data set because of the coarse resolution and the ambiguities in the ice thickness extracted). A systemic evaluation of CIS-and SAR-based ice thickness charts reveals that the correspondence between the SAR ice thickness products and the CIS ice charts for the winter 2002-2003 is better than for 2008-2009. This also applies to the other shown examples (not only to March 2003, which includes some of the training data points). For example, in February 2009 the thick ice area visible in the ice chart north of Prince Edward Island was not detected by the algorithm, and the ice thickness was underestimated. On the other hand, there are smaller, thicker ice segments in the thick ice area based on our algorithm, and also small very thin ice areas. It is impossible to say which interpretation actually corresponds the truth. Further, some areas of thicker and thinner ice are found from different locations, for example in March 2003 ( Fig. 13g and h) the thin and thick ice are on different sides of the strait of Belle Isle (separating Labrador peninsula from Newfoundland) in the northern part of GSL, in the ice chart and ice thickness estimate by our algorithm. Such estimation errors can be due to some rapid strong ice motions which have not been detected by the ice drift estimation algorithm. Such errors can only be corrected by having a higher temporal resolution to enable continuous tracking of the ice drift.

EM measurements
The EM measurements were acquired during the relatively short time periods in late February and early March for both ice seasons studied, as described in Sect. 2.
EM ice thickness measurements essentially include the snow thickness. Unfortunately we did not have any statistics on the snow cover over GSL sea ice. Typically, snow cover over sea ice is not very thick in Northern Hemisphere. HIGH-TSI modeled basin-scale mean snow thickness was roughly 8 % and 20 % of the total mean ice thickness over the winter seasons 2002-2003 and 2008-2009, respectively. Statistical analysis suggests that for ice thicker than 20 cm, the snow thickness is roughly 9-10 % of the ice thickness over Arctic sea ice (Doronin, 1971;Mäkynen et al., 2012). We, therefore, decided not to make any snow corrections for the EM data, i.e. the EM thicknesses used in this study are thus slightly biased by the snow thickness. On the other hand, as we have used the EM data with the snow cover included in the training of the algorithm, the snow cover is also included in the thickness estimates.
The ice thickness values were compared to EM thickness measurements made in the GSL in 2003 and 2009. Applying a sliding median window along the EM measurement line, the EM measurements were smoothed to a resolution similar to that of the SAR data mosaics (500 m). The combined comparisons along all the EM measurement lines are shown in Fig. 15. In the figure the algorithm thickness estimates at The EM data from 2009 are independent, but the 2003 data were used for defining the mapping coefficients from the SAR features and HIGHTSI results to the ice thickness estimates. In general, the SAR algorithm overestimates the ice thickness; for the test data of 2009 the mean overestimation was 11.4 cm, and for the 2003 data 2.7 cm. It also seems that there are some shifts between the EM measurements and SAR thickness estimates. This is probably due to the ice motion and temporal differences (the mosaics contain data from multiple SAR images with different acquisition times) between the EM measurement and the SAR acquisition. Registration inaccuracies between the EM and SAR mosaic data sets are also possible. Assuming the EM measurements are precisely geocoded, the possible registration errors are made in making the mosaic based on the SAR geocoding, which may contain inaccuracies.
Because it seems that direct pixel-wise comparison does not give very convincing error statistics (due to the abovementioned shifts), we also compared the ice thickness distributions based on the SAR algorithm, CIS ice charts, HIGH-TSI experiments, and EM measurements. The comparison was limited to the EM measurement lines. We computed the Kolmogorov-Smirnov statistics between the EM ice thickness distribution, the CIS ice chart mean ice thickness, and the HIGHTSI ice thickness. The Kolmogorov-Smirnov statistics D n between two cumulative distribution functions F n (x) and F (x) is where the cumulative distributions are also functions of the ice thickness. The Kolmogorov-Smirnov statistics are tabulated in Table 3, the mean ice thickness along the EM flight lines are tabulated in Table 4, and the comparison of the ice thickness distributions is shown in Fig. 16. Compared to the CIS ice charts and the solely thermodynamic results, the agreement between the EM data and the new method is excellent. The results of the new method for the test season 2008-2009 ice thickness distribution even have a slightly better Kolmogorov-Smirnov statistic than that for the the season 2002-2003, used in the estimation of the algorithm parameters. The mean ice thickness for the season 2008-2009 data is, however, greater than for the season Fig. 15. The ice SAR algorithm thickness compared to HIGHTSI ice thickness, mean ice thickness derived from digitized CIS ice charts, and EM ice thickness measurements along the EM measurement lines for the two test seasons. The x-axis value is the number of the EM measurement re-sampled to the SAR mosaic resolution of 500 m, in the figure all the EM flight lines for each year have been concatenated in temporal order. 52 Fig. 15. The ice SAR algorithm thickness compared to HIGHTSI ice thickness, mean ice thickness derived from digitized CIS ice charts, and EM ice thickness measurements along the EM measurement lines for the two test seasons (2002-2003 in left panel, 2008-2009 in right panel). The x-axis value is the number of the EM measurements resampled to the SAR mosaic resolution of 500 m; in the figure all the EM flight lines for each year have been concatenated in temporal order.

Conclusions
We developed a novel methodology for deriving sea ice parameters based on a thermodynamic ice model and SAR data in a fully automated manner over GSL. Our algorithms produce ice concentration and thickness estimations utilizing input data from a thermodynamic ice model, SAR-based ice motion, and SAR texture features. The method has been applied for GSL with promising results. The HIGHTSI-based ice thickness is only affected by thermodynamic growth and melt. Ice dynamics was considered by incorporating ice concentration information from SAR data analysis. In the basin scale, the modeled evolution of ice thickness was comparable with that based on the CIS ice charts. The results suggested that the seasonal evolution of ice thickness and area are largely dominated by the weather forcing conditions. The main sources of uncertainty in the model experiments were the amount of snow fall and the oceanic heat flux.
The results of the SAR-based ice concentration estimation were promising. If multiple SAR images are available over The Cryosphere, 6, 1507-1526, 2012 www.the-cryosphere.net/6/1507/2012/ the same area with a reasonable temporal difference, the ice concentration estimates can be further improved by using a temporal minimum (or some other temporal filtering) over a short time period. Using a temporal minimum would reduce the overestimation compared to CIS ice charts, but also, for example, temporal mean or median could be used instead. The reasonable temporal difference depends on the prevailing ice conditions; in the freezing and melting periods with rapid changes, it could be about one day, and in the midseason even several days. Recently, posterior to this study, a more advanced method for ice concentration estimation from C-band HH-channel SAR data has been developed (Karvonen, 2012c), and applying this method would probably improve the ice concentrations estimates. The ice thickness estimation also produced promising results: validation against the local EM measurements showed that the ice thickness based on the new method was much more realistic than that based on the CIS ice charts. A problem in SAR-based ice thickness estimation is, however, to locate areas of thin ice. SAR backscattering and texture for deformed thin ice is sometimes very similar to the values for deformed thick ice areas. We have tried to address this problem by utilizing the ice drift from SAR data together with the thermodynamic ice model. However, the coarse resolution of the ice drift used here cannot yield proper ice drift information over narrow sea areas (straits) present in GSL. Also, the temporal resolution of SAR data at a given location varies from half a day to several days, and thus the SAR-based ice drift can only yield reliable statistics over periods of a few weeks. Rapid dynamic events can currently only be detected on the basis of SAR texture interpretation. It has been shown that radiometer data can be used for estimation of the thin ice thickness (Yu and Rothrock, 1996;Kaleschke et al., 2010;Mäkynen et al., 2010), and combining this information with the SAR data would be useful in better estimation of the thin ice thickness.
The ice thickness in the CIS ice charts is based on the stage of development, and it does not provide very good quantitative reference values. The CIS ice chart concentration seems to be a more reasonable reference for our new products. However, the CIS ice chart thicknesses at least give some idea of large-scale spatial variations in the ice thickness. Because of the coarse nature of the ice thickness description of CIS ice charts (large polygons, and large ice thickness ranges within each polygon) local ice thickness at the resolution of the SAR algorithm is, however, not well described in the ice charts. This fundamental difference also partly explains differences between the SAR algorithm results and the CIS ice charts. Comparisons to ice thickness distributions from resampled EM data and SAR algorithm produced relatively good results, but the results for single point measurements were not good because there seemed to be shifts between the ice thickness estimates and EM values along the flight lines. These shifts are mainly due to temporal differences between the SAR acquisition time at each SAR mosaic grid location and the corresponding EM measurements.
An alternative approach would be to make the ice analysis on the basis of observations and a dynamic-thermodynamic model. Modeling of ice dynamics in a complex coastal or archipelago region is, however, very challenging (Gao et al., 2011). The coarse resolution of the models also involves the problem that ice drift in narrow straits cannot be properly modeled. Our method is, however, not affected by the error sources of modeling sea ice dynamics.
Our estimates are directly affected by the temporal resolution of the SAR data. Because in the SAR mosaics we always use the most recent SAR data, poor temporal resolution may sometimes cause the problem that data over some areas is several days old, and not well representative for the current ice conditions. Sparse temporal data also affect the ice drift estimation, and with too sparse temporal coverage some ice dynamics may not be captured by the algorithm. In the near future satellite constellations, like the X-band Italian COSMO-SkyMed (already in operation), ESA's C-band Sentinel-1 (estimated launch 2013), and Canadian C-band RADARSAT constellation (estimated launch 2015), will provide better temporal coverage, and result in an improved SAR-derived product quality via improved ice drift detection and faster updating of mosaics.
In our opinion the SAR ice thickness estimation algorithm has potential for operational use. The main challenge towards operational sea ice thickness estimation is probably the lack of in-situ measurements for proper algorithm calibration and validation. To build a reliable operational system, almost simultaneous in-situ data with the SAR acquisition (delay no more than a few hours) throughout the ice season would be required for calibration and validation of an algorithm. In this study we only had EM measurements over a short temporal period for both study seasons, and also the spatial coverage of the EM data was limited. For this reason we have also tried to make some quantitative comparison to the ice information extracted from the CIS ice charts. Typically, in-situ measurements are collected in campaigns with a limited spatial and temporal coverage. For development of improved EO algorithms and numerical models, systematic collection of reference sea ice data with a representative temporal and spatial coverage would be desirable.
In the future the algorithm can be improved by computing the ice motion in a finer spatial and temporal scale, taking into account smaller details (narrow straits and river inflows). Also, incorporating data from other EO sources for improved thin ice identification will be considered. The algorithms will also be applied for other ice-covered seas, with validation against available observations.
With the acquisition frequencies of SAR data over GSL during the studied seasons, the proposed algorithms could probably not be used as a stand-alone operational ice monitoring system, but some additional manual corrections would be required. However, we have here demonstrated a methodology that could be used for automated sea ice information retrieval, presuming a high enough SAR data acquisition frequency to produce continuous ice drift estimates and more frequent SAR mosaic updates. With the current SAR acquisition frequency, the mosaic data in some areas are older than one day, even several days. In practice, the SAR data should cover the whole area daily to get reliable ice parameter estimates.
Some current (RADARSAT-2) and forthcoming SAR missions (Sentinel-1 and RADARSAT constellations) also have the possibility to acquire dual-polarized data, i.e. both copolarized (either HH or VV polarization combination) and cross-polarized channels (HV or VH) can be acquired simultaneously. Using cross-polarized SAR data, containing complementary information in addition to the co-polarized SAR data, can improve the sea ice classification, especially in locating open water, thin ice, and deformed ice areas (Scheuchl et al., 2004;Geldsetzer and Yackel, 2009). Utilization of dual-polarized SAR data in our algorithms would produce improved ice concentration and ice thickness estimates. Development and evaluation of an ice thickness estimation algorithm based on dual-polarized SAR data would require dualpolarized SAR data from at least one whole ice season over the area of interest. The major advance of the forthcoming SAR constellations will, however, be that the temporal coverage will be considerably improved as the planned multiple SAR constellations will be realized. Their realization will accomplish the objective of covering the Arctic waters daily at a high resolution.
Acknowledgements. The Gulf of Saint Lawrence experiment was funded and the Canadian NWP data were provided by Environment Canada.
The EM ice thickness data was provided by Fisheries and Oceans Canada. ECMWF is acknowledged for providing us the NWP model results.

RADARSAT
The work of Timo Vihma was supported by the Academy of Finland (contract 259537).
Edited by: D. Feltham