Improved Arctic Sea Ice Freeboard Retrieval From Satellite Altimetry Using Optimized Sea Surface Decorrelation Scales

A growing number of studies are concluding that the resilience of the Arctic sea ice cover in a warming 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 30–220 (mean 105) km. 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 sea ice radar freeboard precision estimated at orbital crossovers by up to 20%. In regions


Introduction
Sea ice extent in the Northern Hemisphere has been declining at an increasingly alarming rate for more than two decades now (Parkinson & DiGirolamo, 2016). Recent studies have recognized that trends and interannual variations in ice extent are extremely sensitive to the pan-Arctic distribution of sea ice thickness (Castro-Morales et al., 2014;Rae et al., 2014). The transition from a sea ice cover dominated by thicker, multi-year ice in the 1980s to an ice cover dominated by thinner, first-year ice in the present day (Tschudi et al., 2016) has amplified interannual fluctuations in the sea ice extent . It has been demonstrated that seasonal forecasts for the sea ice area can be strikingly improved by initializing numerical models with ice thickness observations (Allard et al., 2018;Blockley & Peterson 2018;Fritzner et al., 2019;Msadek et al., 2014;Massonnet et al., 2015;Schröder et al., 2019). Abstract 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 30-220 (mean 105) km. 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 sea ice radar freeboard precision estimated at orbital crossovers by up to 20%. In regions of high SSH uncertainty, biases in CryoSat-2 radar 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.
One of the largest sources of uncertainty in sea ice thickness estimates from altimetry is introduced in the measurement of the sea surface height (SSH) (Armitage & Davidson, 2014;Tilling et al., 2018). The SSH is defined as the ocean free surface elevation with respect to a reference ellipsoid at a sea ice floe and is conventionally interpolated for all ice-covered locations from sea surface tie-points located at the closest leads, that is, openings in the sea ice pack, along the altimeter's orbital track (Laxon et al., 2003;Kwok et al., 2007). Uncertainty in the SSH can be estimated from height variations derived from altimeter returns at leads within a moving window applied along the track (Ricker et al., 2014) (further details in Section 2). However, distances between an ice-covered sample and its closest lead can exceed 200 km along the track, particularly in the compact pack ice of the Central Arctic Ocean (Wernecke & Kaleschke, 2015). In these cases, the SSH uncertainty is constrained only by the deviation of the interpolated sea surface from the local mean measured elevation and can reach 50 cm, varying considerably across the Arctic (Ricker et al., 2014). Averaged to a 25-km grid, the SSH error ranges from approximately 1 to 12 cm (Kwok & Cunningham, 2015;Landy et al., 2020;Ricker et al., 2014).
For several applications, including the reconciliation of sea ice mass balance, polar sea level, climate, and oceanography for scientific purposes and for commercial activities such as operational navigation, accurate determination of the SSH and its uncertainty in ice-covered waters are crucial. In regions with low lead density and high SSH uncertainty, derived freeboards can include long-distance spatially correlated biases (Xia & Xie, 2018). Such biases may either amplify or cancel each other out in different locations, for instance when estimating snow depth from centimeter-scale differences between radar and laser sea ice freeboards . Here, we present an approach for the accurate determination of the SSH which exploits all available lead observations from both the orbital track in focus and additional neighboring tracks. Although we use CryoSat-2 observations to demonstrate the method, the approach can be applied to any contemporary or historical satellite laser or radar altimetry mission (both pulse-limited and SAR). In Section 2 we give an overview of the conventional approaches for estimating the SSH and its uncertainty in sea ice-covered regions. In Section 3 we introduce the satellite and airborne data sets used within our study. In Section 4 we analyze multi-year mean patterns of the Arctic Ocean SSH's spatial and temporal decorrelation length scales from CryoSat-2 lead observations. We then compare these length scales to patterns derived from a high-resolution, nominally 1/12 deg. (∼4 km in the Arctic), simulation of a coupled sea ice-ocean model within the Nucleus for European Modeling of the Ocean (NEMO) framework, where the sea ice component is Louvain-le-Neuve sea ice model version 2 (LIM2) and the ocean component is Ocean Parallélisé (OPA). In Section 5 we apply these length scales in an objective mapping (OM) scheme to determine the optimal SSH at sea ice samples along the altimeter track, from all proximal lead observations on both the track in focus and neighboring tracks. We compare SSH and radar freeboard derived from the objective and conventional analyses in two case months: November 2012 and March 2013. Section 6 validates our results both

Estimating Sea Surface Height in Sea Ice-Covered Locations
Satellite altimeter returns from leads within the sea ice pack are identified by their reflectivity and roughness, in the case of laser altimetry (Kwok et al., 2007(Kwok et al., , 2019, or by their microwave scattering properties in the case of radar altimetry (Laxon et al., 2003(Laxon et al., , 2013. The classification algorithms for identifying leads are generally based on thresholds of parameters for the returning laser or radar echo (Quartly et al., 2019). For CryoSat-2, returns from leads typically make up 1%-15% of all the valid samples Wernecke & Kaleschke, 2015), with higher densities of returns in zones of first-year ice at the pack ice margins where the sea ice concentration is lower. The mean distance from a sea ice sample to the nearest lead along track is approximately 30 ± 60 km (based on our processing chain, see Section 3). However, interpolation distances for the SSH between lead samples actually depend on the 'strictness' of the waveform classifier, with a trade-off between the number of lead samples available and the precision/ accuracy of those height observations. For instance, in the performance analysis of (Wernecke & Kaleschke, 2015) the most liberal classifier produced a lead sample density of 26% but included 13% false positive leads and high variance between proximal SSH observations. With the same dataset, a conservative classifier produced a sample density of only 1% but included zero false positives and low variance between observations. Linear interpolation (Guerreiro et al., 2017;Landy et al., 2017;Lee et al., 2016;Ricker et al., 2014;Xia & Xie, 2018) or regression (Kwok et al., 2007;Lawrence et al., 2018;Tilling et al., 2018) is used to estimate the SSH between lead tie-points. For instance, lead observations from the yellow CryoSat-2 footprints are used to interpolate the SSH at the green 'estimator location', along the altimeter's orbital track, in Figure 1. Data may be discarded where insufficient lead returns are available to reliably interpolate the SSH at a sea ice location (Tilling et al., 2018). Uncertainty on the derived SSH is estimated from the root-mean square (RMS) height of lead returns within a moving window (25 km for instance) along track (Ricker et al., 2014), or by analyzing the RMS of SSH pairs at orbital crossovers (Tilling et al., 2018). For CryoSat-2, the uncertainty on a single SSH measurement has been estimated in the range of 2-50 cm (Ricker et al., 2014). However, this uncertainty is likely to be correlated over wavelengths ∼100 km owing to the length-scale of the SSH interpolation, to the typical distances between lead observations, and to errors in the satellite orbit determination or geophysical corrections (Wingham et al., 2006). Consequently, the uncertainty of the gridded SSH does not reduce as a function of the number of lead observations within a grid cell, it only reduces as a function of the number of individual tracks crossing the grid cell (Lawrence et al., 2018;Tilling et al., 2018).
In the current approach for estimating SSH from pulse-limited radar altimeters, significant positive biases can also be added to the radar range when leads are located outside the nadir point of the satellite 'snag' the radar (Armitage & Davidson, 2014). This sea surface elevation bias ranges from −1 to −4 cm, depending closely on the strictness of the lead classification algorithm, and results in a 10-40 cm overestimate in sea ice thickness if uncorrected (Armitage & Davidson, 2014). The bias can be reliably removed by using the information on the interferometric phase difference of the radar wave travel-time to an off-nadir lead scatterer (Di Bella et al., 2018Bella et al., , 2020; however, of all the radar altimeters only the CryoSat-2 mission has had this capability, and only operating over a small part of the Arctic Ocean. Taking advantage of the interferometric SARIn-mode, around 35% of the lead returns discarded in SAR-mode can be retained, leading to a ∼40% reduction in SSH uncertainty (Di Bella et al., 2018). In the conventional approach only the four leads along the central track (yellow footprints) are used to interpolate sea surface height (SSH) at the sea ice floe location (green footprint). In the proposed approach, all 14 leads (yellow plus blue footprints) acquired within ±2 days at neighboring tracks inside a prescribed sea surface height covariance limit (green circle) around the ice floe are used to compute the SSH. Background is a SAR image from Sentinel-1.

CryoSat-2 Level 2 Processing
We use Baseline-C Level 1B CryoSat-2 waveform observations, details in , for the winter periods between October 2010 and April 2019 obtained from the official ESA science server (accessed in June 2019 at https://science-pds.cryosat.esa.int). SAR-and SARIn-mode observations are retracked by fitting waveforms to echoes simulated from a numerical model for the delay-Doppler SAR altimeter waveform (Landy et al., 2019), using the Lognormal Altimeter Retracking Model (LARM) algorithm described in . A local interpolation of the mean sea surface (MSS) is then removed from the profile of surface heights. The MSS model is a 10-km field obtained from the linear interpolation of all CryoSat-2 lead observations between 2010 and 2019. We apply a three-parameter classification routine to separate CryoSat-2 returns from sea ice and leads, based on the calibrated backscattering coefficient ( 0 E ), the pulse peakiness (PP), and the waveform stack standard deviation (Laxon et al., 2013;Paul et al., 2018;Ricker et al., 2014), as described in . From the multi parameter analysis of (Wernecke & Kaleschke, 2015) the tradeoff for this classifier is around 70% 'leads detected as leads' versus <3% 'ice detected as leads'. Surface heights at leads referenced to the MSS, that is, sea level anomaly (SLA) observations, are retained at this point for further analysis.
To obtain estimates for the radar freeboard, the long-wavelength (>200 km) median profile (which we assume contains residual error from the satellite orbital determination and/or geophysical corrections (Kwok & Cunningham, 2015) or largest-scale features of the dynamic ocean topography) is removed from each CryoSat-2 elevation track. SSH is estimated at sea ice locations by linear interpolation between lead tie-points (Landy et al., 2017). We apply a 25-km low-pass filter to smooth sea surface topography at dense lead clusters and estimate the SSH uncertainty from the RMS height of leads within a 50 km window (Ricker et al., 2014). Radar freeboard is then estimated from the sea ice floe elevation minus the SSH, and the 'single-shot' uncertainty on a freeboard measurement is the RMS of the SSH uncertainty and the total instrument contribution to the elevation error (which is 11.6 cm for SAR-mode and 15.3 cm for SARIn-mode, including orbit and range errors (Wingham et al., 2006)).
Furthermore, we use the SSH observations at leads within the sea ice pack to estimate monthly fields of the Arctic Ocean mean geostrophic current, following the method of (Armitage et al., 2017). Dynamic ocean topography (DOT) within the sea ice-covered zone is estimated from the difference between CryoSat-2 SSH observations and the GOCO5S geoid (Kvas et al., 2019), referenced to the same WGS-84 ellipsoid. The DOT, therefore, contains the long-term offset of the SLA with respect to the geoid. Estimates of the DOT greater than ±2 m are removed before the remaining estimates are sampled onto a 25-km Northern Hemisphere EASE2 grid and smoothed with a 300-km width Gaussian convolution filter. We calculate gradients of the smoothed DOT grid along zonal and meridional axes and convert these to E u and E v vectors of the surface geostrophic current following (Armitage et al., 2017). For the purposes of this study, we calculate the average 'climatological' October-April Arctic Ocean surface current over the entire CryoSat-2 2010-2019 period and mask the region north of 87° latitude due to measurement noise ( Figure 2). The climatological field illustrates the major components of the long-term Arctic Ocean circulation in the winter, including the Beaufort Gyre, Transpolar Drift, East Greenland, and Baffin Island currents, along with the Atlantic and Pacific inflows to the Arctic.

Airborne OIB Ku-Band Data Level 2 Processing
To validate the CryoSat-2 sea ice freeboard observations derived from our new method, we use geolocated Level 1B echograms from the Center for Remote Sensing of Ice Sheets (CReSIS) ultrawideband (UWB) Ku-band airborne radar altimeter, operated on Arctic campaigns by NASA Operation IceBridge (OIB), to generate airborne estimates of radar freeboard coinciding with the satellite. The higher accuracy and resolution airborne freeboard estimates are used as a reference 'truth' from which to evaluate the coinciding satellite estimates of freeboard. The data were accessed from https://data.cresis.ku.edu/#KBRA in January 2020. We selected five airborne campaigns in  that were flown to coincide in space and time with CryoSat-2 overpasses and covered both first-year and multi-year sea ice in the Chukchi and Lincoln Seas, respectively. The CReSIS Kuband radar has a central frequency of 15 GHz (Rodriguez- Morales et al., 2013) and therefore should, in theory, produce a comparable estimate for the radar scattering horizon over snow-covered sea ice to the 13.6 GHz Cryo-Sat-2 radar (Willatt et al., 2011). The flat-surface range resolution of the UWB radar is approximately 4.9 cm in snow and the sensor has an along track sample spacing of approximately 5 m (Paden et al., 2017).

10.1029/2021JC017466
5 of 23 Our processing methodology for the CReSIS radar is built on the algorithm detailed in  to derive snow-ice interface elevation, with several additional steps required to determine the sea ice radar freeboard which we introduce here. We exclude all aircraft segments where the variability of the detrended aircraft altitude is > 0.6 m, or where the mean aircraft pitch or roll is >6°. The local CryoSat-2 MSS is removed from the retracked elevation profile. Radar returns from leads are classified by thresholding waveforms with  0 E and PP above dynamic thresholds. For instance, radar returns classified as leads are shown as red points on the example aircraft profile in Figure 3a. Each threshold is determined from the 99th percentile of  0 E or PP samples, but are no higher than 34 dB or 0.25, respectively, calculated recursively over groups of 12 radar segments (60 km total length). The SSH (red dashed line in Figure 3a) is estimated at ice floes using the method described in Section 3.1, including removing the long-wavelength (200 km) median elevation profile, and radar freeboard (blue bars in Figure 3b) is obtained from ice floe elevations minus the SSH. We exclude all samples located more than 5 km from their nearest lead to prevent the introduction of correlated freeboard biases away from leads. A single airborne freeboard estimate is calculated per ∼300 m coinciding CryoSat-2 footprint (orange boxes in Figure 3b) following (Di Bella et al., 2018).

Auxiliary Data
The EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI SAF) global sea ice concentration climate data record (OSI-450) daily EASE2 gridded observations (accessed from ftp://osisaf.met.no/reprocessed/ ice/conc/v2p0 in January 2020) (Lavergne et al., 2019), are used to filter valid CryoSat-2 observations from the sea ice zone. The OSI SAF global sea ice type record (OSI-403-c) daily polar stereographic gridded observations (accessed from ftp://osisaf.met.no/archive/ice/type in June 2019) (Breivik et al., 2012), are used to identify whether CryoSat-2 or OIB airborne observations are located over first-year or multi-year sea ice.

Decorrelation Length Scales for the Arctic Sea Level Anomaly From Satellite Altimetry
To determine which leads can be used for interpolating the local SLA at sea ice floe locations, we must first define the typical spatial and temporal length scales of the Arctic SLA. The CryoSat-2 lead observations present several challenges for accurately resolving characteristic wavelengths of the SLA signal. Generally, the observations are strongly clustered into groups of 1-10 consecutive valid specular lead returns along the track. The distances between clusters of valid lead returns can also be in excess of 100 km along the track, using our lead classification routine. Adjacent tracks are sampled every 1-2 hr and generally spaced hundreds of kilometers apart. Consequently, at small time and distance lags, we are limited by these sampling considerations and cannot accurately resolve the higher-frequency scales of the SLA. One might expect this to include SLA signatures of ocean circulation features (such as mesoscale eddies and meanders) caused by instability of ocean currents at the scale of the local, first-mode Rossby deformation radius, estimated to be around 5-15 km in the Arctic (Nurser & Bacon, 2014). These mesoscale features can alias the SLA signal, increasing the uncertainty in SLA predicted for nearby leads. However, through the present analysis, we will demonstrate that the Arctic Ocean SLA spatial decorrelation length scales are generally much larger than the local Rossby deformation radius.
Using CryoSat-2 observations for the Arctic SLA at leads within the sea ice pack, we map the winter decadal average spatial (in zonal and meridional directions) and temporal decorrelation length scales of the Arctic Ocean SLA signal. We map the decorrelation length scales onto a 50 km EASE2 grid (Brodzik et al., 2012) covering the Northern Hemisphere above 50° latitude. We select a minimum lag distance of 2.5 km and a lag time of 0.5 days based on the sampling limits of the CryoSat-2 data. The covariance  E of the SLA at lag distance or time E r is defined as: where E n is the number of paired observations at lag distance or time E r , E z is the instantaneous SLA, i E x is the location or time of observation i . A Gaussian function is then fit to the binned covariance to extract key parameters of the SLA decorrelation, following previous studies in the equatorial oceans, for example, (Jacobs et al., 2001). This function can account for non-zero covariance between observation pairs within the few shortest lag bins, which we expect due to the uncorrelated speckle noise properties of 20 Hz CryoSat-2 observations, and asymptotic limit at the covariance amplitude, that is, the random variance of the field. We fit the following Gaussian function to the binned zonal, meridional, and time-dependent covariance estimates obtained from Equation 1: where E a is the covariance amplitude, E s is the covariance at E r = 0, and L/ 3 is the e-folding scale of the SLA. E L is the 'effective range', which defines the lag where  E drops to 5% of the covariance amplitude. The function in Equation 2 is fit to the empirical covariance functions with a bounded nonlinear least-squares optimization algorithm. The lower bound for E s is zero and E L is bounded at the maximum lag distance or time. The quality of fit is determined from the coefficient of determination between empirical and model covariance functions.

Spatial and Temporal Length Scales
To obtain the SLA spatial length scales over the entire Arctic Ocean, we perform the following analysis at monthly intervals for the entire October-April 2010-2019 CryoSat-2 lead observation dataset. For every cell of the 50-km EASE2 grid, we identify all SLA observations in the month within 500 km. Lag distances are employed at 5 km intervals from 2.5 to 502.5 km, including pairs from the same and different tracks. A time limit of 3 days between observation pairs is imposed to maximize the likelihood that observations are correlated in time (Pujol et al., 2016) with sufficient observations remaining available for analysis. By doing this we are limiting the chances of decorrelation in time but searching for spatial correlations over a very wide range. The derived length scales are therefore representative of averaged conditions over a 3-day time window. We construct a matrix of the zonal and meridional distances of all valid observation pairs and sample the covariance for each lag bin along both directional axes using Equation 1. We then fit the Gaussian function in Equation 2 to each empirical function and determine the e-folding length, covariance amplitude, and minimum covariance for the grid cell. Only grid cells with a model r 2 fit >0.3 are retained.
To determine the SLA temporal length scales, we again perform the following analysis at monthly intervals of the CryoSat-2 data. For every cell of the grid, we identify all SLA observations within 100 km and ±30 days of the 15th of the month. The spatial limit of only 100 km is chosen to maximize the likelihood that observations are correlated in space (Pujol et al., 2016), with sufficient observations remaining available for analysis. By doing this we are limiting the chances of decorrelation in space but searching for temporal correlations over a very wide range. The derived time scales are therefore representative of averaged conditions over a 100 km radius. We have to apply a low-pass median filter to observations within 100-km window clusters along the orbital track to reduce the impact of small-scale signal noise along track on the time-dependent decorrelation of the SSH between tracks (Jacobs et al., 2001). Lag times are employed at 1-day intervals from 0.5 to 30.5 days. We construct a matrix of the time difference between all valid observation pairs and sample the covariance for each lag bin using Equation 1. Then fit the Gaussian function in Equation 2 to each empirical function and again determine the e-folding length, covariance amplitude and minimum covariance for the grid cell. Only grid cells with a model r 2 fit >0.3 are retained.

Mean Decorrelation Scales for the Sea Level Anomaly
From the analysis of monthly-mean SLA covariance fields, we find no clear seasonal or interannual patterns in the variability of both the spatial and temporal decorrelation scales. Therefore, as a first estimate here we calculate 2010-2019 'climatological' zonal, meridional, and temporal decorrelation length scales from the weighted mean of the e-folding lengths of all 62 fields (i.e., the total number of analyzed months), with the fit statistics from Equation 2 providing the weights (Figure 4). These climatological scales are used for all remaining analyses. Finally, we smooth the climatological fields with a 3 × 3 grid cell median filter to remove a few remaining anomalies.
The derived zonal, meridional, and temporal e-folding scales compare closely to the estimates of (Pujol et al., 2016) obtained from multiple altimeter missions for sub-polar oceans. For example (Pujol et al., 2016), estimated zonal length scales of 45-100 km for the latitude band between 50 and 70° north, which are comparable to our estimates of 40-120 km for the Barents, Kara, and Laptev Seas (Figures 4a and 4b). Temporal scales (for the same latitude band) of 3-7 days (Pujol et al., 2016) are marginally higher than our estimates from CryoSat-2 of 1-5 days (Figure 4c). Our estimates for the zonal and meridional decorrelation scales (Figures 4a and 4b) match patterns for the first-mode baroclinic Rossby radius obtained from hydrographic observations (Nurser & Bacon, 2014), with higher scales in the Western Arctic (Beaufort Sea region) than on the eastern side north of Svalbard (Figures 4a and 4b). However, the CryoSat-2 e-folding scales of 50-200 km are an order-of-magnitude higher than the baroclinic deformation radius (see Section 4.4) supporting the sub-polar observations of (Chelton et al., 2011). For instance (Chelton et al., 2011), show that eddies can be three times larger than the Rossby radius, suggesting that deformation radii cannot be directly associated with the size of eddies. Our CryoSat-2 data appear to characterize mesoscale anomalies at a scale between baroclinic instabilities and larger features of the geostrophic circulation field. However, the CryoSat-2 data do appear to resolve smaller 10s km features of the SLA signal over the shelf seas, for instance. The SLA decorrelation timescales (Figure 4c) match the typical 1-7 days synoptic period of passing weather systems (Hutchings et al., 2011).
Variations in the characteristic spatial and temporal length scales of the SLA are controlled by the Arctic Ocean's bathymetry, with shallower bathymetry on the shelves introducing additional tidal signals to the SSH that may be uncorrected and cause the signal to decorrelate more rapidly (Armitage et al., 2017). Generally, the patterns of the zonal and meridional length scales are quite similar, although it is particularly evident in the Siberian Seas and also in Hudson Bay and the Greenland Sea that the meridional scale is significantly shorter than the zonal scale (Figures 4a and 4b). This makes sense as the SSH will be less well correlated across the shelf-break than along it and emphasizes the need to apply these two scales independently in the analyses below. It is also interesting that the space and time decorrelation scales appear to be considerably longer in regions covered by ice for most of the year (i.e., the perennial ice zone) than areas of the marginal ice zone (MIZ) with lower sea ice concentrations (see Section 4.4).
The covariance amplitude (Figure 5a) characterizes the standard deviation of the SLA outside the decorrelation timescale shown in Figure 4c, that is, the variability present in the SSH signal over long time periods. It ranges from approximately 6 cm over the Central Arctic Ocean to 15+ cm on the shelf seas. If the SSH is estimated at a location from leads exclusively outside the correlated zone, the uncertainty on the SSH estimate can be no better than this value, which is a salient point because the conventional methods for interpolating SLA (Section 2) have often used length scales well above those shown in Figure 4. The covariance at zero lag ranges from around 2 cm over the central ocean to 6 cm on the shelf seas (Figure 5b). These values represent the characteristic uncertainty on an estimate for the SSH using only lead observations in the immediate vicinity of a location and close in time, from all available tracks. Generally, this includes only a small number of lead observations but with a low sample variance, and the Arctic Ocean mean of 3.8 cm is similar to the estimate of ∼4 cm SSH uncertainty derived from orbit crossover analysis (Tilling et al., 2018).

Interpreting the Decorrelation Scales
We can expect the ocean surface to be 'flat' over a length scale defined by the first mode baroclinic Rossby radius of vertical deformation, which characterizes the approximate scale of boundary currents, eddies, and fronts.
In the weakly-stratified Arctic Ocean and shallow shelf seas, the baroclinic Rossby radius has been determined as only 2-16 km from a climatology of hydrographic observations (Nurser & Bacon, 2014). This is around an order-of-magnitude smaller than the length scales over which SLA is conventionally interpolated along the altimeter's orbital track when deriving sea ice freeboard (see Section 2). Therefore, small-scale dynamic features of the ice-covered Arctic Ocean surface topography cannot always be resolved from dispersed lead observations in along-track altimeter data, only where adjacent tracks including leads are separated by less than a few 10s km. However, sea ice floes can transmit internal stress through the ice cover for several 10s km, interacting with and suppressing dynamic features such as eddies (Meneghello et al., 2017) and enhancing their diameter into the pack ice (Kozlov et al., 2019). Besides, ocean dynamics in ice-covered areas can be influenced by large scale coherent ocean structures for example, (Wilson et al., 2021), so the SLA in ice-covered waters may-in reality-covary over much longer distances than the baroclinic Rossby radius predicts (Chelton et al., 2011;Nurser & Bacon, 2014).
To examine whether this is likely to be the case, we have further analyzed the covariance of SSH fields from a 1/12° global simulation (the ORCA0083-N06 run) of the coupled ocean-sea ice model OPA-LIM2 (Madec et al., 1998;Fichefet & Morales Maqueda, 1997;Goosse & Fichefet, 1999), applying an identical method to the one we applied here for CryoSat-2 ( Figure S3 in Supporting Information S1), every 5 days between 2011 and 2015. The model uses the quasi-uniform, tri-polar ORCA grid (Madec & Imbard 1996) to avoid the singularity associated with the convergence of meridians at the north pole. The grid has 75 vertical levels and a lateral resolution of 2-5 km in the Arctic region ( Figure S1 in Supporting Information S1), which should be sufficient to capture decorrelation length scales of the SLA of order 10s km, indicative of dynamic features (such as eddies, e.g., see Figure S2 in Supporting Information S1) that we may be missing with CryoSat-2. The uniform model SSH fields also do not suffer from the same nonuniform clustered sampling limitations of the altimeter data. This configuration of NEMO has been widely used for Arctic Ocean studies, for example, (Bacon et al., 2015;Kelly et al., 2019;Tsubouchi et al., 2018). The SLA is calculated from the SSH fields with reference to a mean sea surface height model derived from all time slices between 2011 and 2015.
We find the smallest e-folding length scales from NEMO are 10-20 km in the North Atlantic ( Figure S4 in Supporting Information S1), which suggests the model can resolve small-scale dynamical features if they are present. Patterns of the zonal and meridional length scales are remarkably similar to CryoSat-2, with the largest scales in the Central Arctic and much smaller scales in the sub-polar seas. The range of length scales between NEMO and CryoSat-2, of around 20-200 km, are almost identical. There are relatively higher length scales in the East Siberian Sea, Central Arctic, and Hudson Bay, and relatively lower scales in the Southern Beaufort Sea and Baffin Bay, in NEMO. These model findings support previous idealized simulations of the Beaufort Gyre that resulted in eddies emerging with about 100 km scale (Manucharyan & Spall, 2016). Large-scale variability can still dominate the SLA due to basin and gyre scale mechanisms that exaggerate the decorrelation lengths (Jacobs et al., 2001). To examine whether our CryoSat-2 observations may be picking up only the largest gyre-scale features of the SSH, we try low-pass filtering the NEMO SLA to remove features greater than 250 and 125 km and recalculating the length scales ( Figure S5 in Supporting Information S1). Even after removing features >125 km, the derived scales remain 20-100 km within the Arctic and do not reduce to Rossby-like radii (despite these decorrelation scales appearing at other locations, such as the North Atlantic where the model grid is actually coarsest). This implies that length scales obtained from our analysis of the NEMO and CryoSat-2 data without filtering are the dominant length scales of the SLA. A final notable result from the NEMO analysis is that length scales are almost always higher when a grid cell is ice-covered than ice-free. The presence of sea ice reduces covariance amplitudes by 65% and increases decorrelation scales by 20% on average when we test the same locations with and without sea ice ( Figure S6 in Supporting Information S1). This may partly explain the enhanced decorrelation scales measured by CryoSat-2 in the more frequently sea ice-covered Central and Western Arctic.
The apparent decorrelation length and time scales observed by CryoSat-2 are also supported by previous observations of sea ice motion from ice-mounted buoy arrays. Multi-scale drifter arrays deployed in the Beaufort Sea as part of the 2007 SEDNA experiment showed little coherence in ice deformation patterns across spatial scales of 10-100 km, with coherence only appearing at scales exceeding 100 km (Hutchings et al., 2011). The observed coherence between buoys is also typically only lost over synoptic time periods longer than 3-8 days (Hutchings et al., 2011). These evident spatial scales of coherent sea ice motion are >10 times larger than the first mode baroclinic Rossby radius of deformation, reflecting more closely the apparent decorrelation scales of the SSH signal observed by CryoSat-2 ( Figure 4). For example, we find characteristic scales for the SSH signal of 100-150 km and 2-5 days in the Beaufort Sea.

Objectively Mapping the Arctic Ocean Sea Level Anomaly and Radar Freeboard
In a sea ice freeboard processing algorithm, the sea level anomaly at ice-covered locations is conventionally interpolated from adjacent lead observations along the altimeter track. Here we use an objective mapping methodology (Le Traon et al., 1998Pujol et al., 2016) as an alternative approach to estimate the SLA along the CryoSat-2 altimeter track. The method accounts for spatiotemporal decorrelation length scales (Section 4.2) and propagation velocities (Section 3.1) of the sea surface height signal, to determine the best linear SLA estimate and its error from CryoSat-2 observations at proximal leads (on the orbital track in focus and from adjacent tracks).
To limit the computational demand, we use up to 2001 leads to determine SLA at a single ice-covered location. Full details on the OM method and processing algorithm are given in Supporting Information S1. One month of CryoSat-2 Arctic Ocean observations takes approximately 4 days to process on a 56 core 256 GB RAM cluster with this algorithm.
We compare the SLA and radar freeboard obtained for November 2012 and March 2013 with the OM algorithm versus a conventional processing algorithm. The conventional method uses an external MSS model (DTU18) for deriving the SLA at leads, linear interpolation between leads along track smoothed with a low-pass filter, with the SSH uncertainty obtained from the RMS height of leads within a 25-km moving window along track (Landy et al., 2017. To demonstrate the advantages of the OM method in areas of low lead density, in Section 5.1 we do not set a distance limit on the along-track linear interpolation. However, for the analysis of monthly freeboard observations in Section 5.2, we remove all samples from the along-track interpolation >100 km from their nearest lead to follow convention.

Case Study Track on 3 March 2013
We first select a single ascending-orbit CryoSat-2 SAR-mode track at 03:46:51 on 3 March 2013 to illustrate the advantages of the OM method. This track crosses the Arctic Ocean from the Lincoln Sea to the East Siberian Sea. Although the Eastern Arctic sector of the track contains dense lead clusters, our waveform classification algorithm produces only five valid lead returns for the remaining 1,800 km ( Figure 6a). This track represents a case with particularly low lead density and requires interpolation of the SLA over distances of up to 500 km to ice floes from their nearest lead. Normally, most of these observations would be removed because the interpolation distances are too large. Owing to the low lead density, uncertainty on the derived SLA is > 6 cm for the majority of the track (Figure 6a), representing 20%-50% of the final derived radar freeboard (Figure 6d). In areas with sparse leads, the estimated SLA can be tied to single lead observations ( Figure 6a) despite each observation having a random uncertainty up to ∼10 cm (Wingham et al., 2006) and possible bias >4 cm (Armitage & Davidson, 2014).
By applying the OM approach, we sample up to 2001 local observations at leads for every sea ice floe along track and estimate the SLA from the optimal interpolation of them all. Figure 6b illustrates the covariance between the location of every 80th sea ice floe along track and its local sample of SSH observations. Generally, the SLA of the distribution of lead observations around a single ice floe ranges from approximately −0.2 to +0.2 m but is higher in the shallower East Siberian Sea sector (1,800-2,500 km along profile). For this track, 56% of all SLA observations used in the analysis are within half the distance of both x E L and y E L decorrelation length scales and 83% of observations are within the whole distance of x E L and y E L .
The final optimal interpolation (Figure 6c) predictably coincides with most of the lead observations on the focus track because they have a time lag close to zero. However, the covariance matrix between the neighboring lead observations in a local sample provides a weighting on the SLA estimate that reduces the influence of anomalous observations, that is, leads with high estimated measurement error with respect to their neighbors. For instance, the objective SLA estimate is 5 cm lower than the lead at (a) in Figure 6c, indicating this SLA observation may contain a significant error. Isolated leads or lead clusters do not over-influence the objective SLA estimate (Figure 6c) in the same way they can for the conventional method (Figure 6a). In such instances where the objective analysis indicates a lead is under-or over-estimated, as it does at (i), the derived radar freeboards between the OM and conventional approaches contain long wavelength correlated offsets, typically of between −20 and + 20 cm (Figure 6d). The uncertainty estimate for the objective analysis (from Figure S5 in Supporting Information S1) is generally <2 cm, representing <15% of the final derived radar freeboard (Figure 6d) because the SLA is estimated from tens-tohundreds of times more observations than in the conventional approach (Figure 6c). At many locations along this profile, the distance to the nearest lead along track is > 100 km; however, the objective analysis can still resolve a precise estimate for the SLA from leads on other neighboring tracks at short space and time lags.
The OM scheme for determining SLA appears to be particularly successful in resolving discontinuities in SLA (and its uncertainty) at the shelf break and other areas of complex bathymetry. For instance, the SLA does not become aliased when there are insufficient leads to resolve the detailed ocean surface topography, for example, at (iii) in Figures 6a and 6c because it utilizes lead observations from other tracks.

Impact on Radar Freeboard Results
We complete the same comparison between the radar freeboard estimated from a linear interpolation of the SLA along-track and from the objective analysis of all proximal leads from adjacent tracks, for every Cryo-Sat-2 SAR and SARIn mode track in November 2012 and March 2013. For all the analyses in this section, radar freeboards are removed from the along-track method if the interpolation distance from an ice floe to its nearest lead is >100 km. Figure 7 illustrates the differences in radar freeboard and SSH uncertainty between the two methods for March 2013. Pairwise differences in the radar freeboard obtained from the conventional and objective methods can diverge by 2-20 cm along large segments of individual tracks (Figure 7a), where the conventional SLA interpolation is constrained by a low density or biased lead observations along track, but are normally distributed (Figure 7b). On average, the conventional method underestimates the OM method by only 4 mm (Figure 7b; Table 1), constituting an estimated sea ice thickness difference of <10 cm. However, the mean absolute radar freeboard difference is 2.3 cm, which represents a 19% local uncertainty on the mean freeboard and constitutes ∼20 cm uncertainty in sea ice thickness estimated from these freeboards. The largest biases between methods, for instance in the Lincoln Sea, Laptev Sea, and Foxe Basin (Figure 7a), coincide with regions that have a low density of SSH observations at leads (Figure 8c). The SSH uncertainty at ice floes is highest for both methods in these locations (Figure 8d), where the along-track SLA interpolation can be tied to isolated single or clusters of leads (e.g, (i) in Figure 6c). Table 1 shows that the numbers of valid radar freeboard observations used to produce monthly 25 km freeboard grids are 13% higher in November 2012% and 2025% higher in March 2013 when using the OM method. Most of these additional observations have been removed from the along-track interpolation because ice floes are >100 km from their nearest lead. The larger difference in March can be explained by the higher sea ice concentration and lower lead density than in November (Figures 8a and 8c). The higher number of along-track observations (a) 25-km gridded distribution of CryoSat-2 radar freeboards from the objective mapping SSH method minus the conventional approach, including the limit of the multi-year sea ice area in black. (b) Radar freeboard differences between the two methods, from the raw along-track CryoSat-2 sea ice observations. (c) Ratio of the SSH uncertainty estimate from the objective mapping method to the conventional along-track uncertainty estimate, also from the along-track observations. (Note that all observations are removed where the along-track distance to a lead is >100 km). Note. For the conventional method, we remove all freeboard observations where the interpolation distance to the nearest lead is greater than 100 km. 'Number of valid observations' refers to the number of valid along-track freeboard/uncertainty observations used to generate monthly grids, whereas 'gridded data coverage' refers to the total area covered by the grids. ultimately leads to a small improvement in gridded data coverage for both months (Table 1), particularly in regions of low lead density (red markings on Figures 8a and 8c). The mean and standard deviation monthly gridded freeboards are both around 3-6 mm higher for the OM method than using along-track interpolation (Table 1). This would constitute a change in derived sea ice thickness of <10 cm. However, by utilizing many times (typically 1-2 orders-of-magnitude) more SSH observations to determine the optimal SSH at a sea ice floe, the OM method produces around a factor of three reductions in the mean SSH uncertainty estimated at sea ice floe locations (Figure 7c). Gridded on a monthly timescale the SSH uncertainty is estimated to be 47% lower in November 2012 and 40% lower in March 2013 through OM. The objective method accounts for uncorrelated random and correlated systematic errors in the observations, and their reduction to the error estimate at a single ice floe scales directly with the number and covariance of observations used to determine the SLA (see Supporting Information S1). This objective estimate for the uncertainty is therefore based entirely on the observations themselves and does not suffer from the assumptions or conditions of the conventional method, for instance, that the SSH uncertainty is uniform across the Arctic or depends only on the RMS of SSH observations along track (Section 2). The estimated SSH uncertainty is most improved (blue regions in Figures 8b  and 8d) in areas with a high lead density (Figures 8a and 8c) and high SLA decorrelation length scale (Figure 4), such as the Beaufort Sea. The uncertainty remains similar (white regions in Figures 8b and 8d) at locations with a high CryoSat-2 track density, for example, near to the pole hole, whereas the uncertainty is apparently degraded (red areas) in coastal or marginal areas of the sea ice pack. In these locations, the uncertainty obtained from the RMS of lead observations along track may underestimate the true uncertainty on the SSH.

Analysis at Orbital Crossovers
As a first assessment of the precision of the OM method for deriving SLA at ice floes, we identify crossovers of the CryoSat-2 orbit over the Arctic Ocean sea ice pack. Around 4,500 unique crossovers are identified in November 2012 and 13,000 in March 2013 where orbits intersect within 24 hr and CryoSat-2 measurements for each track are no more than 5 km apart. The crossover locations are clustered around the north pole because the polar-orbiting satellite disproportionally crosses itself within a small region north of ∼84° latitude; however, there are rings of crossovers at around 66 and 79° too (Figure 9c).
All pairwise differences in the radar freeboard estimated at crossover locations are normally distributed (Figure 9). The widths of the distributions represent a combination of random noise in the measurements, different orientations of the footprint illuminated by the nadir beam, aliased tidal signals, and additional errors relating to ice motion and possibly radar signal penetration for example, (Willatt et al., 2011). The RMS of the SLA estimated by along-track interpolation at crossover locations are 5.3 and 4.6 cm in November 2012 and March 2013, respectively. This is between the 4 cm reported by (Tilling et al., 2018) for SLA crossovers in the Arctic Ocean and the 6.7 cm reported by (Calafat et al., 2017) for the open ocean with CryoSat-2 in LRM mode. The OM scheme reduces the RMS of the radar freeboard measured at crossover locations by 12% (from 7.0 to 6.2 cm) and 19% (from 6.9 down to 5.5 cm) compared to the conventional approach, in November 2012 ( Figure 9a) and March 2013 (Figure 9b), respectively. The improved SLA estimation reduces a portion of the radar freeboard uncertainty. However, because the RMS of the radar freeboard is not as low as the RMS of the SLA at crossovers, this suggests a significant portion of the freeboard uncertainty is ice-related (i.e., including the effects of ice motion, signal penetration, speckle noise, and retracking uncertainties).

Independent Validation of Radar Freeboards
We use independent radar freeboards derived from the CReSIS airborne Ku-band radar flown on OIB Arctic campaigns (Section 3.2) to compare the accuracies of the conventional and objective SLA mapping techniques for deriving sea ice radar freeboards. The Ku-band radar freeboards are used here rather than the official OIB L4 total (snow plus sea ice) freeboard and thickness product, so that we do not have to correct freeboards for snow depth and delayed radar wave propagation through the snow layer . The OIB L4 snow depths contain known biases (Newman et al., 2014) and fixed snow densities may introduce further systematic uncertainties (Mallett et al., 2020). So, we use airborne radar freeboards to mimic the CryoSat-2 radar measurements as closely as Where there is a sufficient density of leads along track to resolve the SSH topography with low uncertainty, the conventional along-track method typically produces similar results to the objective analysis. However, for two of the campaigns, on 26th March 2012 (CryoSat-2 in SARIn mode, with a 1-hr time difference between aircraft and satellite passes) and on 26th March 2014 (SAR mode, with a 4.5-hr time difference), satellite radar freeboards from the objective analysis and conventional approaches diverged significantly. In these cases, the density of leads along track was lower. Here, we want to analyze which, if either, the method accurately reproduces the airborne radar freeboards.
The OIB campaign on 26 March 2012 measured mostly multi-year sea ice with some first-year ice in the "Wingham Box" (Figure 10). The SSH estimated for this track with objective analysis was between 4 and 10 cm lower (Figure 10a) than the SSH estimated with the conventional along-track approach but, owing to a low density of leads in the region, both methods included a relatively high uncertainty estimate (Figure 10b). (Note the uncertainty shown in Figures 10a and 10b characterizes the precision of the estimated SSH, whereas the remaining analyses here characterize its accuracy). Figure 10c illustrates the airborne and two satellite radar freeboard profiles after a moving average filter with 2 km width is applied. The distribution of radar freeboards from the OM method matches the airborne observations better than the distribution obtained from the conventional along-track method (Figures 10d and 10e). The conventional method appears to underestimate the airborne freeboards by 8.9 cm (mean difference, MD) because it overestimates the SSH (Figure 10b). In comparison, the MD between the objectively mapped CryoSat-2 freeboards and OIB is −3.4 cm. The accuracy of the new method (RMSE = 11.2 cm) is improved by around 25% versus OIB relative to the conventional method (RMSE = 14.9 cm), at the 2-km lengthscale of the freeboard observations, after the SSH bias is reduced.
The OIB campaign on 26 March 2014 measured predominantly multi-year ice in the Lincoln Sea ( Figure 11). The SSH estimated for this track with objective analysis was between 0 and 8 cm lower (Figures 11a) than the SSH estimated with the conventional along-track approach and both methods produced lower uncertainty estimates at one end of the section, owing to a cluster of leads to the north (Figures 11b). Figure 11c illustrates the airborne and two satellite radar freeboard profiles, after a moving average filter with 2 km width is applied. The CryoSat-2 freeboards from each method exhibit long-wavelength (>100 km) correlated differences (Figures 11b). Again, the distribution of radar freeboards from the OM method matches the airborne observations better than the distribution obtained from along-track interpolation between leads (Figures 11d and 11e). The conventional method underestimates radar freeboard by MD = 11.1 cm, in comparison to 4.8 cm for the new method. The accuracy of the OM method (RMSE = 13.8 cm) is improved by around 20% versus OIB relative to the conventional method (RMSE = 17.1 cm), at the 2-km length-scale of the freeboard observations, after the SSH bias is reduced.
Our independent evaluation of the CryoSat-2 radar freeboards for both OIB campaigns demonstrates that long-wavelength errors, caused for example, by a low density of valid lead returns along track, off-nadir lead errors, or errors in the L1B CryoSat-2 orbital/geophysical corrections, can introduce significant biases into derived radar freeboards using conventional SSH interpolation. In both cases where the conventional and objective SSH  mapping techniques diverged, the objective estimate more accurately reproduced the radar freeboards observed from OIB aircraft observations.

Prospects for Further Improvement
There are several avenues worth exploring to further improve the objective mapping of SLA and sea ice freeboard in ice-covered seas. It may be valuable to use leads at adjacent tracks for mapping SLA at ice floes with ICESat-2 because regions of sea ice further than 10 km from their nearest lead reference along track are currently discarded (Kwok et al., 2019). This leaves some areas such as the densely-concentrated multi-year ice of the Central Arctic occasionally missing valid observations for example, . However, the higher density of SSH observations from ICESat-2 may enable improved characterization of the spatiotemporal characteristics of the SLA signal, and possibly also its seasonal variation, for other altimetry missions. It may also enable discarded ICESat-2 segments in lead-sparse regions like the Canadian Arctic Archipelago or Lincoln Sea to be retained (Kwok et al., 2019). Now that Sentinel-3A and -3B are operating together with CryoSat-2 over a portion of the Arctic Ocean (Lawrence et al., 2019), there is strong potential for characterizing the SLA signal in more detail combining all three sensors.
Our current implementation of the OM scheme takes approximately five days to process one month of CryoSat-2 SLA estimates at ice floes over the Arctic Ocean. This is around 3-4 orders-of-magnitude longer than conventional along-track SSH interpolation methods. Since end users often use a pan-Arctic gridded sea ice freeboard or thickness product on monthly timescales, the gridded freeboard could be obtained from the difference of optimally interpolated SLA for example, (Rose et al., 2019) and ice floe height at grid locations, constrained by the spatiotemporal length scales derived here (Section 4). This would reduce the computational cost of deriving freeboards along the satellite track. However, there are several avenues of research requiring sea ice freeboard/ thickness observations closer to the altimeter 'sensor resolution' (i.e., ∼300 m along track) rather than a downsampled (10s km) monthly product. For instance, to obtain daily-weekly freeboards (Gregory et al., 2021) or snow depth (Guerreiro et al., 2016) from multiple altimeters and for data assimilation into short-term sea iceocean forecasting models (Fiedler et al., 2021).
It may then be possible to obtain results with similar accuracy but only processing one in E n sea ice floe samples along track, or using a reduced sample size of SSH observations, with considerable improvements in computation speed. Validation results in Section 7.2 indicate that the OM scheme provides a significant benefit over conventional along-track interpolation in regions of lower lead density (for instance, in zones of high sea ice concentration like the Central Arctic). However, the benefits are lower in regions of high lead density, such as the marginal ice zone, and simple ocean topography. So, it may be possible to apply the OM scheme only in the areas it provides the most benefit, reducing computing requirements.

Implications of the Objective Mapping Method for Deriving Inter-Mission Data Products and to the Reanalysis of Historic Altimeter Missions
The depth of snow on Arctic sea ice has been estimated from the offset between laser freeboards from ICESat-2  or radar freeboards from the Ka-band AltiKa (Lawrence et al., 2018) and radar freeboards from CryoSat-2. Whilst we do not expect errors in the determination of SLA to introduce pan-Arctic uniform biases between satellites (as we identify only small freeboard biases in Table 1), the along-track correlated errors from interpolation between leads (the non-zero regions in Figure 7a) could realistically introduce local biases to the derived inter-mission snow depths. These biases may either amplify or cancel each other out. Objective mapping, therefore, offers the prospect of combining SSH observations from multiple altimeter missions (Le Traon et al., 2003;Pujol et al., 2016): calculating constant inter-mission biases if present but, more importantly, preventing local biases where leads are sparse or have high uncertainty. Errors will be limited at mission crossover locations ( Figure 9) and systematic uncertainties should also be reduced on gridded freeboard differences.
The objective SLA mapping scheme offers the most improvement over conventional techniques where sea ice concentrations are highest and/or SLA observations at leads have the largest height uncertainty. For instance, the most obvious changes in a gridded freeboard in Figure 7a occur in the perennially-ice covered zone north of Greenland and Canada. Historic pulse-limited radar altimeter missions, such as Envisat or ERS-1 and -2 (and the ongoing mission AltiKa), have an effective footprint of 2-8 km and are therefore more sensitive to 'snagging' than the SAR-focused CryoSat-2 (Section 2). The instruments on ERS-1/-2 also had a larger bandwidth than recent missions, meaning their range resolution was lower with specular lead reflections more likely to be aliased in the recorded waveforms. By estimating the optimal local SLA from a greater number of proximal lead observations, accumulated from multiple tracks, the OM scheme should effectively reduce the random uncertainty from noise and waveform aliasing and the systematic uncertainty from snagging. Since a higher number of leads are used for each SLA interpolation, a more conservative lead classification can be employed for a smaller sample of higher accuracy SLA observations. Improvements on the conventional SLA interpolation schemes for these missions should, in theory, be larger than we have found for CryoSat-2.

Conclusions
The conventional method for estimating sea ice freeboard with altimetry uses only lead observations along the satellite orbital track to interpolate the local sea surface height at ice floes. The SSH uncertainty is typically estimated from the RMS of lead height observations within a moving window along track. Here we have introduced a method to determine the optimal interpolation of local SLA at sea ice floes using all valid proximal lead observations both on the orbital track in focus and other adjacent tracks. The principal novelty of our approach lies in the fact that, instead of 1D along-track SLA interpolation, the objective mapping is applied locally with the purpose to improve sea ice freeboard retrievals. The OM method assumes that spatial and temporal properties of the SLA in ice-covered waters can be predicted with a characteristic Gaussian covariance function. The decorrelation length scales and signal propagation velocities that constrain this function in the Arctic Ocean are obtained by analyzing historic SLA measurements acquired by the CryoSat-2 radar altimeter from lead locations between 2010 and 2019. The best linear least-squares solution for the SLA at each ice floe is determined from all valid SLA observations, weighted by their covariance with the floe location, their covariance with neighbors (i.e., to identify anomalies), and their random and systematic observation errors.
By exploiting a greater number of leads for interpolating the SLA, it is possible to use a stricter pulse peakiness classification threshold-discarding more ambiguous lead waveforms without compromising the height estimate and its uncertainty. For instance, the OM method can effectively reduce off-nadir lead biases on the derived Cry-oSat-2 SLA when corrections from the interferometric phase are not available. Applying our method to the Arctic Ocean in November 2012 and March 2013, the monthly gridded SSH uncertainty can be reduced by around two times in comparison to conventional uncertainty estimates. The RMS of the radar freeboard at orbital crossovers is reduced by 10%-20%. Where independent airborne observations are available and the coinciding OM and conventional SSH estimates from CryoSat-2 give different results, we find the objective method improves satellite-derived freeboard accuracy by 20%-25%. The OM method is also capable of resolving much finer-scale detail of the SSH signal in areas of complex ocean topography such as the circumpolar shelf break. However, inversion of the SSH observation matrix is computationally expensive, so our current software takes around five days on a cluster to process SSH at ice floes for one month of pan-Arctic CryoSat-2 data.
Objective SSH mapping produces the largest improvements at local scales and may therefore enable accurate sea ice freeboards to be estimated at kilometer-scale resolutions along the satellite track. With CryoSat-2 maneuvered to align with ICESat-2 from August 2020, it will be valuable to inter-compare the SSH between these two satellites and test whether OM can reduce local biases in the freeboard offsets. Furthermore, the scheme offers considerable potential for new missions such as CRISTAL and for reprocessing ice freeboards from historic pulse-limited radar altimetry missions, including AltiKa, Envisat, ERS-1 and -2, where SSH observations are more likely to have off-nadir lead biases, higher noise, and are regularly spaced >100 km along track.

Data Availability Statement
This work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk) for producing NEMO model simulations. We thank Andrew Coward for developing and running the NEMO model experiments and George Nurser for providing global fields for the baroclinic Rossby radius from (Nurser & Bacon, 2014). The along-track CryoSat-2 observations of SSH for 2011-2019 are available by request from the authors. Gridded