Mapping Arctic Sea Ice Thickness: A New Method for Improved Ice Freeboard Retrieval from Satellite Altimetry

A growing number of studies are concluding that the resilience of the Arctic sea ice cover in a warming climate is essentially controlled by its thickness. Satellite radar and laser altimeters have allowed us to routinely monitor sea ice thickness across most of the Arctic Ocean for several decades. However, a key uncertainty remaining in the sea ice thickness retrieval is the error on the sea surface height (SSH) which is conventionally interpolated at ice floes from a limited number of lead observations along the altimeter’s orbital track. Here, we use an objective mapping approach to determine sea surface height from all proximal lead samples located on the orbital track and from adjacent tracks within a neighborhood of 10s of kilometers. The patterns of the SSH signal’s zonal, meridional, and temporal decorrelation length scales are obtained by analyzing the covariance of historic CryoSat-2 Arctic lead observations, which match the scales obtained from an equivalent analysis of high-resolution sea ice-ocean model fields. We use these length scales to determine an optimal SSH and error estimate for each sea ice floe location. By exploiting leads from adjacent tracks, we can increase the SSH precision estimated at orbital crossovers by a factor of three. In regions of high SSH uncertainty, biases in CryoSat-2 sea ice freeboard can be reduced by 25% with respect to coincident airborne validation data. The new method is not restricted to a particular sensor or mode, so it can be generalized to all present and historic polar altimetry missions.


26
A growing number of studies are concluding that the resilience of the Arctic sea ice cover in a warming 27 climate is essentially controlled by its thickness. Satellite radar and laser altimeters have allowed us to 28 routinely monitor sea ice thickness across most of the Arctic Ocean for several decades. However, a 29 key uncertainty remaining in the sea ice thickness retrieval is the error on the sea surface height (SSH) 30 which is conventionally interpolated at ice floes from a limited number of lead observations along the 31 altimeter's orbital track. Here, we use an objective mapping approach to determine sea surface height 32 from all proximal lead samples located on the orbital track and from adjacent tracks within a 33 neighborhood of 10s of kilometers. The patterns of the SSH signal's zonal, meridional, and temporal 34 decorrelation length scales are obtained by analyzing the covariance of historic CryoSat-2 Arctic lead 35 observations, which match the scales obtained from an equivalent analysis of high-resolution sea ice-36 ocean model fields. We use these length scales to determine an optimal SSH and error estimate for each 37 sea ice floe location. By exploiting leads from adjacent tracks, we can increase the SSH precision 38 estimated at orbital crossovers by a factor of three. In regions of high SSH uncertainty, biases in 39 CryoSat-2 sea ice freeboard can be reduced by 25% with respect to coincident airborne validation data. 40 The new method is not restricted to a particular sensor or mode, so it can be generalized to all present 41 and historic polar altimetry missions. 42 43 1. Introduction leads typically make up 1 to 15% of all the valid samples (Wernecke & Kaleschke, 2015;Passaro, et 149 al., 2018), with higher densities of returns in zones of first-year ice at the pack ice margins where the 150 sea ice concentration is lower. The mean distance from a sea ice sample to the nearest lead along track 151 is approximately 30 ± 60 km (based on our processing chain, see Section 3). However, interpolation 152 distances for the SSH between lead samples actually depend on the 'strictness' of the waveform 153 classifier, with a trade-off between the number of lead samples available and the precision/accuracy of 154 those height observations. For instance, in the performance analysis of (Wernecke & Kaleschke, 2015) 155 the most liberal classifier produced a lead sample density of 26% but included 13% false positive leads 156 and high variance between proximal SSH observations. With the same dataset, a conservative classifier 157 produced a sample density of only 1% but included zero false positives and low variance between 158 observations. 159 Linear interpolation (Ricker, et al., 2014;Lee, et al., 2016;Guerreiro, et al., 2017;160 Xia & Xie, 2018) or regression (Kwok, et al., 2007;Tilling, et al., 2018;Lawrence, et al., 2018) is used 161 to estimate the SSH between lead tie-points (Fig. 1). A low-pass filter can be used to smooth the final 162 surface at clusters of leads, where noise may introduce artificially rough sea surface topography. Data 163 may be discarded where insufficient lead returns are available to reliably interpolate the SSH at a sea 164 ice location (Tilling, et al., 2018). Uncertainty on the derived SSH is estimated from the root-mean 165 square (RMS) height of lead returns within a moving window (25 km for instance) along track (Ricker, 166 et al., 2014), or by analyzing the RMS of SSH pairs at orbital crossovers (Tilling, et al., 2018). For 167 CryoSat-2, the uncertainty on a single SSH measurement has been estimated in the range of 2-50 cm. 168 However, this uncertainty is likely to be correlated over wavelengths >100 km owing to the length-169 scale of the SSH interpolation, to the typical distances between lead observations and to errors in the 170 satellite orbit determination or geophysical corrections (Wingham, et al., 2006). Consequently, random 171 errors in the SSH observations in sea ice zones cannot simply be reduced by accumulating observations 172 of the leads along the track. 173 In the current approach for estimating SSH from pulse-limited radar altimeters, significant positive 174 biases can also be added to the radar range when leads located outside the nadir point of the satellite 175 'snag' the radar (Armitage & Davidson, 2014). This sea surface elevation bias ranges from -1 to -4 cm, 176 depending closely on the strictness of the lead classification algorithm, and results in a 10-40 cm 177 overestimate in sea ice thickness if uncorrected (Armitage & Davidson, 2014). The bias can be reliably 178 removed by using information on the interferometric phase difference of the radar wave travel-time to 179 an off-nadir lead scatterer (Di Bella, et al., 2018;Di Bella, et al., 2020); however, of all the radar 180 altimeters only the CryoSat-2 mission has had this capability, and only operating over a small part of 181 the Arctic Ocean. Taking advantage of the interferometric SARIn-mode, around 35% of the lead 182 returns discarded in SAR-mode can be retained, leading to a ~40% reduction in SSH uncertainty (Di 183 Bella, et al., 2018). 184 algorithm described in (Landy, et al., 2020). A local interpolation of the mean sea surface (MSS) is 192 then removed from the profile of surface heights. The MSS model is a 10-km field obtained from the 193 linear interpolation of all CryoSat-2 lead observations between 2010 and 2019. We apply a three-194 parameter classification routine to separate CryoSat-2 returns from sea ice and leads, based on the 195 calibrated backscattering coefficient ( 0 ), the pulse peakiness (PP) and the waveform stack standard 196 deviation (Laxon, et al., 2013;Ricker, et al., 2014;Paul, et al., 2018), as described in (Landy, et al., 197 2020). Surface heights at leads referenced to the MSS, i.e. sea level anomaly (SLA) observations, are 198 retained at this point for further analysis. 199

Data & Preprocessing
To obtain estimates for the radar freeboard, the long-wavelength (>200 km) median profile (which we 200 assume contains residual error from the satellite orbital determination and/or geophysical corrections 201 (Kwok & Cunningham, 2015) or largest-scale features of the dynamic ocean topography) is removed 202 from each CryoSat-2 elevation track. SSH is estimated at sea ice locations by linear interpolation 203 between lead tie-points . We apply a 25-km low-pass filter to smooth sea surface 204 topography at dense lead clusters and estimate the SSH uncertainty from the RMS height of leads 205 within a 50 km window (Ricker, et al., 2014). Radar freeboard is then estimated from the sea ice floe 206 elevation minus the SSH, and the 'single-shot' uncertainty on a freeboard measurement is the root-207 sum-square of the SSH uncertainty and speckle noise (which is 11.6 cm for SAR-mode and 15.3 cm for 208 SARIn-mode (Wingham, et al., 2006)). 209 Furthermore, we use the SSH observations at leads within the sea ice pack to estimate monthly fields of 210 the Arctic Ocean mean geostrophic current, following the method of (Armitage, et al., 2017). Dynamic 211 ocean topography (DOT) within the sea ice-covered zone is estimated from the difference between 212 CryoSat-2 SSH observations and the GOCO5S geoid (Kvas, et al., 2019), referenced to the same WGS-213 84 ellipsoid. The DOT therefore contains the long-term offset of the SLA with respect to the geoid. 214 Estimates of the DOT greater than ±2 m are removed before the remaining estimates are sampled onto 215 a 25-km Northern Hemisphere EASE2 grid and smoothed with a 300-km width Gaussian convolution 216 filter. We calculate gradients of the smoothed DOT grid along zonal and meridional axes and convert 217 these to and vectors of the surface geostrophic current following (Armitage, et al., 2017). For the 218 purposes of this study, we calculate the average 'climatological' October-April Arctic Ocean surface 219 current over the entire CryoSat-2 2010-2019 period and mask the region north of 87° latitude due to 220 measurement noise (Fig. 2) Operation IceBridge (OIB), to generate airborne estimates of radar freeboard coinciding with the 228 satellite. The data were accessed from https://data.cresis.ku.edu/#KBRA in January 2020. We selected 229 five airborne campaigns in 2011, 2012 and 2014 (all in March) that were flown to coincide in space 230 and time with CryoSat-2 overpasses and covered both first-year and multi-year sea ice in the Chukchi 231 and Lincoln Seas, respectively. The CReSIS Ku-band radar has a central frequency of 15 GHz 232 (Rodriguez-Morales, et al., 2013) and therefore should, in theory, produce a comparable estimate for 233 the radar scattering horizon over snow-covered sea ice to the 13.6 GHz CryoSat-2 radar (Willatt, et al., 234 2011). The flat-surface range resolution of the UWB radar is approximately 4.9 cm in snow and the 235 sensor has an along track sample spacing of approximately 5 m (Paden, et al., 2017). 236 Our processing methodology for the CReSIS radar is built on the algorithm detailed in (Landy,et al.,237 2020) to derive snow-ice interface elevation, with several additional steps required to determine the sea 238 ice radar freeboard which we introduce here. We exclude all aircraft segments where the variability of 239 the detrended aircraft altitude is >0.6 m, or where the mean aircraft pitch or roll is >6°. The local 240 CryoSat-2 MSS is removed from the retracked elevation profile. Radar returns from leads are classified 241 by thresholding waveforms with 0 and PP above dynamic thresholds (Fig. 3a). Each threshold is 242 determined from the 99 th percentile of 0 or PP samples, but are no higher than 34 dB or 0.25, 243 respectively, calculated recursively over groups of twelve radar segments (60 km total length). The 244 SSH is estimated at ice floes using the method described in Section 3.1 and radar freeboard is obtained 245 from ice floe elevations minus the SSH (Fig. 3b). We exclude all samples located more than 5 km from 246 their nearest lead to prevent the introduction of correlated freeboard biases away from leads. A single 247 airborne freeboard estimate is calculated per ~300 m coinciding CryoSat-2 footprint (Fig. 3b) (Breivik, et al., 2012), are used to identify whether 256 CryoSat-2 or OIB airborne observations are located over first-year or multi-year sea ice. 257

258
To determine which leads can be used for interpolating the local SLA at sea ice floe locations, we must 259 first define the typical spatial and temporal length scales of the Arctic SLA. The CryoSat-2 lead 260 observations present several challenges for accurately resolving characteristic wavelengths of the SLA 261 signal. Generally, the observations are strongly clustered into groups of 1-10 consecutive valid specular 262 lead returns along the track. The distances between clusters of valid lead returns can also be in excess 263 of 100 km along track, using our lead classification routine. Adjacent tracks are sampled every 1-2 264 hours and generally spaced hundreds of kilometers apart. Consequently, at small time and distance 265 lags, we are limited by these sampling considerations and cannot accurately resolve the higher-266 frequency scales of the SLA. One might expect this to include SLA signatures of ocean circulation 267 features (such as mesoscale eddies and meanders) caused by instability of ocean currents at the scale of 268 the local, first-mode Rossby deformation radius, estimated to be around 5-15 km in the Arctic (Nurser 269 & Bacon, 2014). These mesoscale features can alias the SLA signal, increasing the uncertainty in SLA 270 predicted for nearby leads. However, through the present analysis we will demonstrate that the Arctic 271 Ocean SLA spatial decorrelation length scales are generally much larger than the local Rossby 272 deformation radius. 273 Using CryoSat-2 observations for the Arctic SLA at leads within the sea ice pack, we map the winter 274 decadal average spatial (in zonal and meridional directions) and temporal decorrelation length scales of 275 the Arctic Ocean SLA signal. We map the length scales onto a 50 km EASE2 grid (Brodzik,et al.,276 2012) covering the Northern Hemisphere above 50 degrees latitude. We select a minimum lag distance 277 of 2.5 km and lag time of 0.5 days based on the sampling limits of the CryoSat-2 data, although 278 following our analysis the smallest scales identified were several times larger than these values. The 279 covariance of the SLA at lag distance or time is defined as: 280 (1) Where is the number of paired observations at lag distance or time , is the instantaneous SLA, 281 is the location or time of observation . The SLA signal is modelled with a Gaussian function, 282 following previous studies in the equatorial oceans, e.g. (Jacobs, et al., 2001). This model can account 283 for non-zero covariance between observation pairs within the few shortest lag bins, which we expect 284 due to the uncorrelated speckle noise properties of 20 Hz CryoSat-2 observations, and asymptotic limit 285 at the covariance amplitude, i.e. the random variance of the field. We fit the following Gaussian model 286 to the empirical zonal, meridional, and time-dependent covariance functions obtained from Eq. (1): 287 Where is the covariance amplitude, is the covariance at = 0, and /√3 is the e-folding scale of the 288 SLA. is the 'effective range', which defines the lag where drops to 5% of the covariance amplitude 289 and is applied as the first zero-crossing of the imposed SLA signal decorrelation in Eq. 6 (Section 5). 290 The model in Eq. (2) is fit to the empirical covariance functions with a bounded nonlinear least-squares 291 optimization algorithm. The lower bound for is zero and is bounded at the maximum lag distance or 292 time. The quality of fit is determined from the optimized coefficient of determination between 293 empirical and model covariance functions.  To obtain spatial patterns for the characteristic SLA spatial length scale over the entire Arctic Ocean, 296 we perform the following analysis at monthly intervals for the entire Oct-Apr 2010-2019 CryoSat-2 297 lead observation dataset. For every cell of the 50-km EASE2 grid we identify all SLA observations in 298 the month within 500 km. Lag distances are employed at 5 km intervals from 2.5 to 502.5 km, 299 including pairs from the same and different tracks. A time limit of 3 days between observation pairs is 300 imposed to maximize the likelihood that observations are correlated in time (Pujol, et al., 2016) with 301 sufficient observations remaining available for analysis. By doing this we are limiting the chances of 302 decorrelation in time but searching for spatial correlations over a very wide range. The derived length 303 scales are therefore representative of averaged conditions over a 3-day time window. We construct a 304 matrix of the zonal and meridional distances of all valid observation pairs and sample the covariance 305 for each lag bin along both directional axes using Eq. (1). We then fit the Gaussian model in Eq.
(2) to 306 each empirical function and determine the e-folding length, covariance amplitude and minimum 307 covariance for the grid cell. Only grid cells with a model r 2 fit >0.3 are retained. 308 To determine the SLA temporal length scales, we again perform the following analysis at monthly 309 intervals of the CryoSat-2 data. For every cell of the grid, we identify all SLA observations within 100 310 km and ±30 days of the 15 th of the month. The spatial limit of only 100 km is chosen to maximize the 311 likelihood that observations are correlated in space (Pujol, et al., 2016), with sufficient observations 312 remaining available for analysis. By doing this we are limiting the chances of decorrelation in space but 313 searching for temporal correlations over a very wide range. The derived time scales are therefore 314 representative of averaged conditions over a 100 km radius. Since the orbit time for a single CryoSat-2 315 pass over the Arctic Ocean is a matter of minutes, we do not analyze the time-dependent correlation 316 between SLA observations along the track. In contrast, we apply a low-pass median filter to 317 observations within 100-km window clusters along the orbital track, to reduce the impact of small-scale 318 signal noise along track on the time-dependent decorrelation of the SSH between tracks. Lag times are 319 employed at 1-day intervals from 0.5 to 30.5 days. We construct a matrix of the time difference 320 between all valid observation pairs and sample the covariance for each lag bin using Eq. (1). Applying 321 Eq. (2) we then fit the Gaussian model to each empirical function and again determine the e-folding 322 length, covariance amplitude and minimum covariance for the grid cell. Only grid cells with a model r 2 323 fit >0.3 are retained. 324

Mean decorrelation scales for the sea level anomaly
325 From the analysis of monthly-mean SLA covariance fields, we find no clear seasonal or interannual 326 patterns in the variability of both the spatial and temporal correlation scales. Therefore, as a first 327 estimate we calculate 2010-2019 'climatological' zonal, meridional, and temporal decorrelation length 328 scales from the weighted mean of the e-folding lengths of all 62 fields (i.e., the total number of 329 analyzed months), with the optimized model fit statistics providing the weights. We use these 330 climatological scales for all remaining analysis. Another reason why climatological scales must be used 331 is because there are typically too few lead observations proximal to sea ice samples from which to 332 calculate local contemporary spatiotemporal correlation length scales, which can lead to high 333 uncertainty in the derived SLA. Finally, we smooth the climatological fields with a 3 x 3 grid cell 334 median filter to remove a few remaining anomalies. 335 The derived zonal, meridional, and temporal e-folding scales compare closely to the estimates of 336 (Pujol, et al., 2016) obtained from multiple altimeter missions for sub-polar oceans. For example, 337 (Pujol, et al., 2016) estimated zonal length scales of 45-100 km for the latitude band between 50 and 70 338 degrees north, which are comparable to our estimates of 40-120 km for the Arctic peripheral seas, the 339 Barents, Kara, and Laptev Seas ( Fig. 4a and b). Temporal scales (for the same latitude band) of 3-7 340 days (Pujol, et al., 2016) are marginally higher than our estimates from CryoSat-2 of 1-5 days for the 341 peripheral seas (Fig. 4c). Our estimates for the zonal and meridional decorrelation scales ( Fig. 4a and  342 b) match patterns for the first-mode baroclinic Rossby radius obtained from hydrographic observations 343 (Nurser & Bacon, 2014), with higher scales in the Western Arctic (Beaufort Sea region) than on the 344 eastern side north of Svalbard ( Fig. 4a and b). However, the CryoSat-2 e-folding scales of 50-200 km 345 are an order-of-magnitude higher than the baroclinic deformation radius (see Section 4.4) supporting 346 the sub-polar observations of (Chelton, et al., 2011). For instance, (Chelton, et al., 2011) show that 347 eddies can be three times larger than the Rossby radius, suggesting that deformation radii cannot be 348 directly associated with the size of eddies. Our CryoSat-2 data appear to characterize mesoscale 349 anomalies at a scale between baroclinic instabilities and larger features of the geostrophic circulation 350 field. However, the CryoSat-2 data do appear to resolve smaller 10s km features of the SLA signal over 351 the shelf seas, for instance. The SLA decorrelation timescales (Fig. 4c) match the typical 1-7 day 352 synoptic period of passing weather systems (Hutchings, et al., 2011). 353 Variations in the characteristic spatial and temporal length scales of the SLA are controlled by the 354 Arctic Ocean's bathymetry, with shallower bathymetry on the shelves introducing additional tidal 355 signals to the SSH that may be uncorrected and cause the signal to decorrelate more rapidly (Armitage, 356 et al., 2017). Generally, the patterns of the zonal and meridional length scales are quite similar, 357 although it is particularly evident in the Siberian Seas and also in Hudson Bay and the Greenland Sea 358 that the meridional scale is significantly shorter than the zonal scale ( Fig. 4a and b). This makes sense 359 as the SSH will be less well correlated across the shelf-break than along it and emphasizes the need to 360 apply these two scales independently in the analyses below. It is also interesting that the space and time 361 decorrelation scales appear to be considerably longer in regions covered by ice for most of the year (i.e. 362 the perennial ice zone) than areas of the marginal ice zone (MIZ) with lower sea ice concentrations. 363 The covariance amplitude (Fig. 5a) characterizes the standard deviation of the SLA outside the 364 correlation timescale shown in Figure 4c, i.e. the variability present in the SSH signal over long time 365 periods. It ranges from approximately 6 cm over the Central Arctic Ocean to 15+ cm on the shelf seas. 366 If the SSH is estimated at a location from leads exclusively outside the correlated zone, the uncertainty 367 on the SSH estimate can be no better than this value, which is a salient point because the conventional 368 methods for interpolating SLA (Section 2) have often used length scales well above those shown in 369 Figure 4. The covariance at zero lag ranges from around 2 cm over the central ocean to 6 cm on the 370 shelf seas (Fig. 5b). These values represent the characteristic uncertainty on an estimate for the SSH 371 using only lead observations in the immediate vicinity of a location and close in time, from all available 372 tracks. Generally, this includes only a small number of lead observations but with a low sample 373 variance, and the Arctic Ocean mean of 3.8 cm is similar to the estimate of ~4 cm SSH uncertainty 374 derived from orbit crossover analysis (Tilling, et al., 2018). 375 376 We can expect the ocean surface to be 'flat' over a length scale defined by the first mode baroclinic 377

Interpreting the decorrelation scales
Rossby radius of vertical deformation, which characterizes the approximate scale of boundary currents, 378 eddies, and fronts. In the weakly-stratified Arctic Ocean and shallow shelf seas, the baroclinic Rossby 379 radius has been determined as only 2-16 km from a climatology of hydrographic observations (Nurser 380 & Bacon, 2014). This is around an order-of-magnitude smaller than the length scales over which SLA 381 is conventionally interpolated along the altimeter's orbital track when deriving sea ice freeboard (see 382 Section 2). Therefore, small-scale dynamic features of the ice-covered Arctic Ocean surface 383 topography cannot reliably be resolved from dispersed lead observations in along-track altimeter data 384 (let alone in adjacent time-lagged tracks). However, sea ice floes can interact with and suppress 385 dynamic features such as eddies (Meneghello, et al., 2017), so the SLA in ice-covered waters mayin 386 realitycovary over much longer distances than the baroclinic Rossby radius predicts (Chelton, et al., 387 2011;Nurser & Bacon, 2014). 388 To examine whether this is likely to be the case, we have further analyzed the covariance of SSH fields 389 from a 1/12° global simulation (the ORCA0083-N06 run) of the coupled ocean-sea ice model OPA-390 LIM2 (Madec, et al., 1998 The covariance amplitudes of the NEMO SLA also has a similar pattern to those derived from CryoSat-419 2 ( Fig. S4) but are consistently ~5 cm lower reflecting the absence of measurement noise in model SSH 420 fields. The NEMO amplitudes also underestimate the high CryoSat-2 amplitudes measured in Hudson 421 Bay and the Canadian Arctic, for example. One final notable result from the NEMO analysis is that 422 length scales are almost always higher when a grid cell is ice-covered than ice-free. Presence of sea ice 423 reduces covariance amplitudes by 65% and increases decorrelation scales by 20% on average when we 424 test the same locations with and without sea ice (Fig. S6). This may partly explain the enhanced 425 decorrelation scales measured by CryoSat-2 in the perennially ice-covered Central and Western Arctic.
is an observation, i.e. the true SLA Φ and its observation error . is the covariance 450 matrix of all selected observations, and is the covariance vector between the observations and field 451 to be estimated: 452 The covariance matrix in Eq. 5 is constructed between all SLA observations. The 'single-shot' random 485 error associated with a 20 Hz CryoSat-2 observation is 11.6 cm for SAR mode and 15.3 cm for 486 SARIn mode (Wingham, et al., 2006). This is combined with the long-wavelength error estimated 487 as 25% of the signal variance Eqs. (4) and (7). Finally, after deriving individual SLA estimates for every CryoSat-2 footprint 492 classified as sea ice along a track, we smooth the resulting profile with a low-pass filter whose window 493 is limited to the mean of the local SSH e-folding scales our waveform classification algorithm produces only five valid lead returns for the remaining 1800 km 506 (Fig. 7a). This track represents a case with particularly low lead density and requires interpolation of 507 the SLA over distances of up to 500 km to ice floes from their nearest lead (if all floes are to be 508 included in the analysis). Owing to the low lead density, uncertainty on the derived SLA is >6 cm for 509 the majority of the track (Fig. 7a), representing 20-50% of the final derived radar freeboard (Fig. 7d). 510 In areas with sparse leads, the estimated SLA can be tied to single lead observations (Fig. 7a)  By applying the objective mapping approach, we sample up to 2001 local observations at leads for 514 every sea ice floe along track and estimate the SLA from the optimal interpolation of them all. Figure  515 7b illustrates the covariance between the location of every 80 th sea ice floe along track and its local 516 sample of SSH observations. Generally, the SLA of the distribution of lead observations around a 517 single ice floe ranges from approximately -0.2 to +0.2 m but is higher in the shallower East Siberian 518 Sea sector (1800-2500 km along profile). For this track, 56% of all SLA observations used in the 519 analysis are within half the distance of both and correlation length scales and 83% of 520 observations are within the whole distance of and . 521 The final optimal interpolation (Fig. 7c) predictably coincides with most of the lead observations on the 522 focus track because they have a time lag close to zero. However, the covariance matrix between the up 523 to 2001 neighboring lead observations in a local sample provides a weighting on the SLA estimate that 524 reduces the influence of anomalous observations, i.e., leads with high estimated measurement error 525 with respect to their neighbors. For instance, the objective SLA estimate is 5 cm lower than the lead at 526 (i) in Figure 7c, indicating this SLA observation may contain significant error. Single isolated leads or 527 lead clusters do not over-influence the objective SLA estimate (Fig. 7c) in the same way they can for 528 the conventional method (Fig. 7a). In such instances where the objective analysis indicates a lead is 529 under-or over-estimated, as it does at (i), the derived radar freeboards between the new and 530 conventional approaches contain long wavelength correlated offsets, typically of between -20 and +20 531 cm (Fig. 7d). The uncertainty estimate for the objective analysis (from Eq. 7) is generally <2 cm, 532 representing <15% of the final derived radar freeboard (Fig. 7d), because the SLA is estimated from 533 tens-to-hundreds of times more observations than in the conventional approach (Fig. 7c). The 534 uncertainty is notably higher at (ii) in Figure 7c because and are <150 km at this location, so the 535 number of available SLA observations is significantly lower than the maximum 2001 permitted and 536 their variance is larger (Fig. 7b). 537 The new scheme for determining SLA enables the radar freeboard to contain greater along-track 538 variability than the conventional scheme (Fig. 7d) because the estimated SLA is not fixed over long 539 (>100 km) distances along track by isolated single or clusters of leads. The new scheme appears to be 540 particularly successful resolving discontinuities in SLA (and its uncertainty) at the shelf break and 541 other areas of complex bathymetry. For instance, the SLA does not become aliased when there are 542 insufficient leads to resolve the detailed ocean surface topography, e.g., at (iii) in Fig. 7a and c. 543 544 We complete the same comparison between the conventional SLA estimated from a linear interpolation 545 along-track and from the objective analysis of all proximal leads from adjacent tracks, for every 546 CryoSat-2 SAR and SARIn mode track in March 2013. Pairwise differences in the radar freeboard 547 obtained from the conventional and objective methods are normally distributed ( Fig. 8a and b) but 548 comprise long-wavelength (10-500 km) correlated offsets between the methods in either direction. 549 (Note we do not discard any freeboard observations based on their distance to the nearest along-track 550 lead for this analysis). The radar freeboards can diverge by >5 cm along large segments of individual 551 tracks (Fig. 8a), where the conventional SLA estimate is poorly constrained through biased 552 observations or a low density of lead observations along track (Section 6.1). The conventional method 553 is essentially as likely to underestimate the objectively mapped SLA as overestimate it (Fig. 8a). On 554 average, the conventional method underestimates the objective mapping method by ~1 cm (Fig. 8b), 555

Analysis of entire month
constituting an estimated sea ice thickness difference of only ~10 cm. However, the mean absolute 556 radar freeboard difference is 3.3 cm, which represents a 27% local uncertainty on the mean freeboard 557 and constitutes >30 cm uncertainty in sea ice thickness estimated from these freeboards. The biases 558 between SSH mapping techniques appear to be independent of location, although there is a pattern of 559 positive freeboard differences (mean = +2.5 cm) in the multi-year ice zone north of Canada and the 560 largest differences are evident in coastal regions (Fig. 8a). These areas coincide with shallower tidal 561 zones that have high SLA variability over short temporal and spatial scales (e.g. Fig. 5) and/or have the 562 lowest available density of SSH observations at leads (Wernecke & Kaleschke, 2015). 563 By utilizing many times (typically 1-2 orders-of-magnitude) more SSH observations to determine the 564 optimal SSH at a sea ice floe, the objective mapping method produces a factor of three reduction in the 565 estimated SSH uncertainty (Fig. 8c) CryoSat-2 measurements for each track are no more than 5 km apart. The crossover locations are 580 clustered around the north pole, because the polar-orbiting satellite disproportionally crosses itself 581 within a small region north of ~84 degrees latitude; however, there are rings of crossovers at around 66 582 and 79 degrees too (Fig. 9c). 583 All pairwise differences in the SLA or radar freeboard estimated at crossover locations are normally 584 distributed (Fig. 9). The widths of the distributions represent a combination of random noise in the 585 measurements, orthogonal sensing footprints for crossing orbits, aliased tidal signals, andin the case 586 of the radar freeboardsadditional errors relating to ice motion and possibly radar signal penetration 587 e.g. (Willatt, et al., 2011). The new objective mapping scheme reduces the RMS of the SSH estimated 588 at crossover locations by 70% compared to the conventional approach, from 4.6 down to 1.4 cm (Fig.  589   9a). The RMS of crossovers for the conventional scheme is close to the 4 cm reported by (Tilling, et 590 al., 2018). It is not surprising that the RMS is reduced through objective analysis, as the SSH is 591 estimated at ice floes from all available leads on all proximal tracks. However, with our new scheme 592 the SSH compared at a crossover is still, in all cases, from an optimal interpolation of nearby 593 observations rather than actual lead observations, so the improvement remains impressive. The 594 objective mapping scheme reduces the RMS of the radar freeboard measured at crossover locations by 595 19% compared to the conventional approach, from 6.9 down to 5.5 cm (Fig. 9b). The improved SSH 596 estimation reduces a portion of the radar freeboard uncertainty. However, because the new scheme 597 reduces the RMS of radar freeboard at crossovers by significantly less than it reduces the RMS of SSH 598 at crossovers, this suggests around three quarters of the total uncertainty in freeboard measurements at 599 crossovers is ice-related (i.e. including the effects of ice motion, signal penetration, speckle noise and 600 retracking uncertainties). 601 602 We use independent radar freeboards derived from the CReSIS airborne Ku-band radar flown on OIB 603

Independent validation of radar freeboards
Arctic campaigns (Section 3.2) to compare the accuracies of the conventional and objective SSH 604 mapping techniques. The Ku-band radar freeboards are used here rather than the official OIB L4 total 605 (snow plus sea ice) freeboard and thickness product, so that we do not have to correct freeboards for 606 snow depth and delayed radar wave propagation through the snow layer (Landy, et al., 2020). The OIB 607 L4 snow depths contain known biases (Newman, et al., 2014) and fixed snow densities may introduce 608 further systematic uncertainties (Mallett, et al., 2020). So, we use airborne radar freeboards to mimic 609 the CryoSat-2 radar measurements as closely as possible and limit the chances of introducing further 610 systematic biases into our comparisons. Satellite radar freeboards are obtained with both the objective 611 and conventional SSH interpolation methods along CryoSat-2 tracks coinciding with five processed 612 OIB campaigns in 2011, 2012 and 2014. Of the five coinciding tracks, three produced similar radar 613 freeboard profiles between the objective and conventional methods suggesting that the conventional 614 along-track approach was sufficient to resolve the SSH in these cases. The more sophisticated but less 615 computationally efficient objective analysis is not always necessary. However, for two of the 616 campaigns, on 26 th March 2012 (CryoSat-2 in SARIn mode, with a 1-hour time difference between 617 aircraft and satellite passes) and on 26 th March 2014 (SAR mode, with a 4.5-hour time difference), 618 satellite radar freeboards from the objective analysis and conventional approaches diverged 619 significantly. Here, we want to analyze which, if either, method accurately reproduces the airborne 620 radar freeboards. 621 The OIB campaign on 26 th March 2012 measured mostly multi-year sea ice with some first-year ice in 622 the 'Wingham Box' (Fig. 10). The SSH estimated for this track with objective analysis was between 4 623 and 10 cm lower (Fig. 10a) than the SSH estimated with the conventional along-track approach but, 624 owing to a low density of leads in the region, both methods included a relatively high uncertainty 625 estimate (Fig. 10b). (Note the uncertainty in Figure 10a and b characterizes the precision of the 626 estimated SSH, whereas the remaining analysis here characterizes its accuracy). Figure 10c illustrates 627 the airborne and two satellite radar freeboard profiles, after a moving average filter with 2 km width is 628 applied. There is some correlation between the airborne and satellite observations in places, but it is 629 very clear that the distribution of radar freeboards from the objective mapping method matches the 630 airborne observations far better than the distribution obtained from the conventional along-track 631 method ( Fig. 10d and e). The conventional method appears to underestimate the airborne freeboards by 632 8.9 cm (mean difference, MD), because it overestimates the SSH (Fig. 10b). In comparison, the MD 633 between the objectively mapped CryoSat-2 freeboards and OIB is -3.4 cm. The accuracy of the new 634 method (RMSE = 11.2 cm) is improved by around 25% versus OIB relative to the conventional method 635 (RMSE = 14.9 cm), at the 2-km length-scale of our averaged freeboard observations. 636 The OIB campaign on 26 th March 2014 measured predominantly multi-year ice in the Lincoln Sea (Fig.  637   11). The SSH estimated for this track with objective analysis was between 0 and 8 cm lower (Fig. 11a) 638 than the SSH estimated with the conventional along-track approach and both methods produced lower 639 uncertainty estimates at one end of the section, owing to a cluster of leads to the north (Fig. 11b). 640 Figure 11c illustrates the airborne and two satellite radar freeboard profiles, after a moving average 641 filter with 2 km width is applied. Like in the 2012 comparison, the CryoSat-2 freeboards from each 642 method exhibit long-wavelength (>100 km) correlated differences (Fig. 11b). Again, it is clear that the 643 distribution of radar freeboards from the objective mapping method match the airborne observations 644 better than the distribution obtained from along-track interpolation between leads (Fig. 11d and e). The 645 conventional method underestimates radar freeboard by MD = 11.1 cm, in comparison to 4.8 cm for the 646 new method. The accuracy of the objective mapping method (RMSE = 13.8 cm) is improved by around 647 20% versus OIB relative to the conventional method (RMSE = 17.1 cm), at the 2-km length-scale of 648 our averaged freeboard observations. 649 Our independent evaluation of the CryoSat-2 radar freeboards for both OIB campaigns demonstrates 650 that long-wavelength errors, caused for example by a low density of valid lead returns along track, off-651 nadir lead errors, or errors in the L1B CryoSat-2 orbital/geophysical corrections, can introduce 652 significant biases into derived radar freeboards using the conventional SSH mapping approach. In both 653 cases where the conventional and objective SSH mapping techniques diverged, the objective estimate 654 more accurately reproduced the radar freeboards observed from OIB aircraft observations. 655 ICESat-2, because regions of sea ice further than 10 km from their nearest lead reference along track 660 are currently discarded (Kwok, et al., 2019). This leaves some areas such as the densely-concentrated 661

Discussion
multi-year ice of the Central Arctic occasionally missing valid observations e.g., . 662 However, the higher density of SSH observations from ICESat-2 may enable an improved 663 characterization of the spatiotemporal characteristics of the SLA signal, and possibly also its seasonal 664 variation, for other altimetry missions. It may also enable discarded ICESat-2 segments in lead-sparse 665 regions like the Canadian Arctic Archipelago or Lincoln Sea to be retained (Kwok, et al., 2019). Now 666 that Sentinel-3A and -3B are operating together with CryoSat-2 over a portion of the Arctic Ocean 667 (Lawrence, et al., 2019), there is strong potential for characterizing the SLA signal in more detail 668 combining all three sensors. Moreover, CryoSat-2 has been maneuvered to coincide more frequently 669 with the ground track of ICESat-2 (as part of the CRYO2ICE Project) which could enable the direct 670 intercomparison of SLA characteristics. 671 It is unlikely our assumption that systematic error between altimeter tracks is a maximum 25% of the 672 signal variance holds in all situations (Section 5). The systematic offset between tracks will be greater 673 when (i) orbital errors are higher (Wingham, et al., 2006) are therefore more sensitive to 'snagging' than the SAR-focused CryoSat-2 (Section 2). The 705 instruments on ERS-1/-2 also had a larger bandwidth than recent missions, meaning their range 706 resolution was lower with specular lead reflections more likely to be aliased in the recorded 707 waveforms. By estimating the optimal local SLA from a greater number of proximal lead observations, 708 accumulated from multiple tracks, our new scheme should effectively reduce the random uncertainty 709 from noise and waveform aliasing and the systematic uncertainty from snagging. Since a higher 710 number of leads are used for each SLA interpolation, a more conservative lead classification can be 711 employed for a smaller sample of higher accuracy SLA observations. Improvements on the 712 conventional SLA interpolation schemes for these missions should, in theory, be larger than we have 713 found for CryoSat-2. 714 demonstrate that the SSH uncertainty can be reduced by around three times in comparison to 734 conventional uncertainty estimates. The root-mean square of interpolated SSH pairs at orbital 735 crossovers is reduced by a factor of three and radar freeboard crossover RMS is reduced by 20%. 736

Conclusions
Where independent airborne observations are available and the coinciding new and conventional SSH 737 estimates from CryoSat-2 give different results, we find the objective method improves satellite-738 derived freeboard accuracy by 20-25%. The new method is also capable of resolving much finer-scale 739 detail of the SSH signal in areas of complex ocean topography such as the circumpolar shelf break. 740 However, inversion of the SSH observation matrix is computationally expensive, so our current 741 software takes around five days on a cluster to process SSH at ice floes for one month of pan-Arctic 742 CryoSat-2 data. 743 Objective SSH mapping produces the largest improvements at local scales and may therefore enable 744 accurate sea ice freeboards to be estimated at kilometer-scale resolutions along the satellite track. With 745 CryoSat-2 maneuvered to align with ICESat-2 from August 2020, it will be valuable to inter-compare 746 the SSH between these two satellites and test whether objective mapping can reduce local biases in the 747 freeboard offsets. Furthermore, the scheme offers considerable potential for new missions such as 748 CRISTAL and for reprocessing ice freeboards from historic pulse-limited radar altimetry missions, 749 including AltiKa, Envisat, ERS-1 and -2, where SSH observations are more likely to have off-nadir 750 lead biases, higher noise and are regularly spaced >100 km along track. 751 752 Figure 1. Schematic diagram illustrating the conventional and proposed new methods for interpolating 939 the sea surface height at sea ice locations. In the conventional approach only the four leads along the 940 central track (yellow footprints) are used to interpolate SSH at the sea ice floe location (green 941 footprint). In the proposed approach, all 14 leads (yellow plus blue footprints) acquired within ±2 days 942 at neighboring tracks inside a prescribed sea surface height covariance limit (green circle) around the 943 ice floe are used to compute the SSH. Background is a SAR image from Sentinel-1. 944    SSH method minus the conventional approach, including the limit of the multi-year sea ice area in 974 black. (b) Radar freeboard differences between the two methods, from the raw along-track CryoSat-2 975 sea ice observations. (c) Ratio of the SSH uncertainty estimate from the objective mapping method to 976 the conventional along-track uncertainty estimate, also from the along-track observations. 977