Mapping of sea ice concentration using the NASA NIMBUS 5 Electrically Scanning Microwave Radiometer data from 1972–1977

The Electrically Scanning Microwave Radiometer (ESMR) instrument onboard the NIMBUS 5 satellite was a one-channel microwave radiometer that measured the 19.35 GHz horizontally polarized brightness temperature (T_B) from 11 December 1972 to 16 May 1977. The original tape archive data in swath projection have recently been made available online by the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC). Even though the ESMR was a predecessor of modern multi-frequency radiometers, there are still parts of modern processing methodologies which can be applied to the data to derive the sea ice extent globally.

Here, we have reprocessed the entire dataset using a modern processing methodology that includes the implementation of pre-processing filtering, dynamical tie points, and a radiative transfer model (RTM) together with numerical weather prediction (NWP) for atmospheric correction. We present the one-channel sea ice concentration (SIC) algorithm and the model for computing temporally and spatially varying SIC uncertainty estimates. Post-processing steps include resampling to daily grids, land-spillover correction, the application of climatological masks, the setting of processing flags, and the estimation of sea ice extent, monthly means, and trends. This sea ice dataset derived from the NIMBUS 5 ESMR extends the sea ice record with an important reference from the mid-1970s. To make it easier to perform a consistent analysis of sea ice development over time, the same grid and land mask as used for EUMETSAT's OSI-SAF SMMR-based sea-ice climate data record (CDR) were used for our ESMR dataset. SIC uncertainties were included to further ease comparison to other datasets and time periods.

We find that our sea ice extent in the Arctic and Antarctic in the 1970s is generally higher than those available from the National Snow and Ice Data Center (NSIDC) Distributed Active Archive Center (DAAC), which were derived from the same ESMR dataset, with mean differences of 240 000 and 590 000 km², respectively. When comparing monthly sea ice extents, the largest differences reach up to 2 million km². Such large differences cannot be explained by the different grids and land masks of the datasets alone and must therefore also result from the differences in data filtering and algorithms, such as the dynamical tie points and atmospheric correction.

The new ESMR SIC dataset has been released as part of the ESA Climate Change Initiative (ESA CCI) program and is publicly available at https://doi.org/10.5285/34a15b96f1134d9e95b9e486d74e49cf (Tonboe et al., 2023).

1 Introduction

Several sea ice concentration (SIC) algorithms have been developed for passive microwave data (PMW). They differ in the usage of, e.g., frequencies and polarizations of the PMW data (Comiso et al., 1997) or the usage of static or dynamic tie points (Parkinson et al., 2004; Tonboe et al., 2016; Lavergne et al., 2019). From the resulting datasets, it is apparent that the Arctic sea ice extent (SIE) in September has been decreasing at a rate of about 12 % per decade since the launch of modern satellite multi-frequency microwave radiometers in 1978 (Onarheim et al., 2018; Stroeve and Notz, 2018). This negative sea ice trend in the Northern Hemisphere started in the 1970s (Rayner et al., 2003; Walsh et al., 2017), though regional trends can differ, as seen for example within the Barents Sea (Chapman and Walsh, 1991). In the Antarctic there are large regional differences in SIE trends, but, until recently, the overall trend was positive due to sea ice dynamics (Turner et al., 2009; Sun and Eisenman, 2021). This has changed, however, in the last decade as a result of several record lows, and, as such, overall trends have shifted to a more homogeneous pattern (Schroeter et al., 2023) and are now slightly negative in summer, November to February (NDJF). Until recently, the slightly positive trend was believed to be part of the long-term natural variability that overshadowed the effects of global warming starting in the 1960s (Wang et al., 2019; Stammerjohn et al., 2008; Thompson and Solomon, 2002; Ferreira et al., 2015; Singh et al., 2019; Fogt et al., 2022). In order to fully understand the drivers of sea ice variability, extending the sea ice data record backwards in time is essential.

Globally, SIE information prior to the satellite data record came largely from ice charts and ship observations. While there have been efforts to include these data in long-term assessments of sea ice change, the data are typically provided at relatively coarse spatial and temporal resolutions (1° grid) (Walsh et al., 2019), interpolating in both time and space (Titchner and Rayner, 2014). Only satellite-based datasets offer the ability to cover both hemispheres at improved spatial and temporal resolutions and generally show consistency in processing methods (Lavergne et al., 2019). The SIC derived from NIMBUS 5 Electrically Scanning Microwave Radiometer (ESMR) data was previously processed by Parkinson et al., 2004. Here, we apply a new processing method that is comparable to the EUMETSAT/ESA CCI SIC record from 1978 and onwards (see Andersen et al., 2006; Tonboe et al., 2016). Compared to Parkinson et al. (2004), this method reduces atmospheric noise regionally over both ice and water surfaces and uses the pre-processed data to develop a SIC algorithm calibration that is effective in removing both instrument drift and offsets. Seasonal sea ice signature variations are removed by using dynamical tie points. Lastly, the algorithm calculates temporally and spatially varying uncertainty estimates. The ESMR SIC data are presented on the same grid and with the same masking as the EUMETSAT/ESA CCI record, which makes these two records directly comparable. This and the modern processing chain mentioned above warrant the reprocessing presented in this article.

In Sect. 2, the satellite and reanalysis data are described, including the formatting and initial filtering of the data. Section 3 introduces the radiative transfer model (RTM) used for the atmospheric correction, while Sect. 4 describes the dynamical tie points, the SIC algorithm with uncertainty estimations, the land-spillover method, and data flags assigned during post-processing. In Sect. 5, the resulting SIC dataset is presented and compared to other datasets. Finally, Sect. 6 consists of a discussion, and Sect. 7 provides the conclusions of this work.

2 The NIMBUS 5 ESMR instrument and data

The NIMBUS 5 ESMR instrument was a cross-track scanner that measured at 78 scan positions perpendicular to the flight track with a maximum incidence angle of about 64° on both sides (NASA GSFC, 2016). No direct observations at nadir were made; the closest positions were at ±0.7°. Its near-circular orbit had a height of about 1112 km and an inclination of 81°. The phased-array antenna dimensions were 85.5×83.3 cm, and the spatial resolution was about 25 km near nadir, increasing to about 160×45 km at the edges of the swath (NASA GSFC, 2016). The full swath was about 3100 km, with the varying incidence angle and the spatial resolution giving very good (unprecedented) daily coverage in polar regions with no gaps, i.e., no pole holes. The ESMR onboard the NIMBUS 5 satellite was a one-channel 19.35 GHz horizontally polarized microwave radiometer operating from 11 December 1972 until 16 May 1977 (1617 d) with some interruptions (see the list of days with missing files in Appendix A2). Due to a hot-load anomaly, there are major data gaps between March to May and again in August 1973. Another major data gap occurs from 3 June until 14 September 1975 because the ground segment was used for receiving NIMBUS 6 data instead (Parkinson et al., 1999). When operation was resumed in September 1975, the instrument was only operated approximately every other day. From late 1976 to the end of the mission, operation was highly irregular. The last file in the dataset is from 16 May 1977. The data have recently been made available online by NASA in the original tape archive format (TAP files).

2.1 Formatting and co-location of brightness temperatures and ECMWF ERA5 data

The ESMR data were retrieved from the online data archive of the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC) (NASA GSFC, 2016). This dataset contains, along with a number of instrument and geographical parameters, calibrated 19.35 GHz brightness temperatures expressed in units of kelvin. The raw data were recovered by NASA from the magnetic tapes, called calibrated brightness temperature tapes (CBTTs), where they were stored in the original binary TAP file format, each file corresponding to a particular orbit (NASA GSFC, 2016).

All variables in the TAP files were read using online NASA software and converted to NetCDF format without changing the original data structure, creating raw ESMR NetCDF files. Each data point in the TAP file was matched with the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 re-analysis data (Bell et al., 2020; Hersbach et al., 2020) that were nearest in time and space and appended to the raw ESMR NetCDF file serving as input to the processing chain. The resulting data are structured in arrays line by line (across track). Appendix A1 summarizes the variables included in the NetCDF files.

2.2 Initial filtering and correction of brightness temperatures

NASA provides a correction of the brightness temperature data to account for the lobe structure, antenna loss, and angular brightness temperature (T_B) variation (NASA CR, 1974). According to NASA, the correction was needed because “The cause of the gross variations in antenna properties which were observed soon after launch has been determined to be a cross-polarized grating lobe […] The problem does not exist for the near-nadir beam positions so these positions are unaffected. […] An empirical calibration has been developed which removes the effect of the lobe structure and antenna loss, which vary with position, and roughly corrects for angular variations in viewing geometry.” (NASA CR, 1974, p. 400).

Originally, it was planned to use only lobe-corrected T_B values with their natural angular dependency, but we did not find a way to extract this from the NASA-provided dataset. Essentially, only the combined lobe and angular correction, which is a function of brightness temperature, can be removed altogether from the data NASA provides. Thus, the T_B values do not vary as a function of incidence angle, as would be expected for T_B values from the sea surface and sea ice.

Despite the corrections done by NASA, the ESMR data still contain erroneous T_B values, scan lines, sudden jumps in the calibration, and other obvious artifacts.

Since the ESMR data contain corrupted data and erroneous scan lines, filtering is needed before the data can be used for sea ice mapping. The filters that we apply are described in Eqs. (1)–(4). They are applied in the same order as described here. If only a single data point or a couple of scan lines are affected, only those data points and scan lines are removed from the file. If the whole file is corrupted then it is deleted.

An initial analog filter is used for filtering erroneous T_B values and scan lines. The filter is based on the 16 analog voltage entries in the data. The NIMBUS 5 ESMR user's guide (Sabatini, 1972) does not explain very well what the 16 entries really are, or what unit the voltages are stored in, but jumps in these analog signals correspond to anomalous T_B values. Our analog filter computes the absolute gradient in the analog signals, and anything over a value of 10 is removed. This threshold was estimated experimentally.

Next, several other filters are employed using the processed T_B (in K) from the previous step. The filters are applied in the following order.

Data that are non-physical and outside the expected range for sea and ice surfaces are removed. Only data points that lie inside the range specified in Eq. (1) are kept:

$$\begin{array}{}\text{(1)}& \mathrm{90}\mathrm{K}<T_{\mathrm{B}}\le \mathrm{273.15}\mathrm{K}.\end{array}$$

The next filter removes erroneous scan lines (across track rows). Consecutive scan lines should not differ by more than 50 K, as shown in Eq. (2). The threshold of 50 K was estimated experimentally.

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{{\sum }_j^n|T_{{\mathrm{B}}_{j,i+\mathrm{1}}}-T_{{\mathrm{B}}_{j,i}}|}{n}}\le \mathrm{50}\mathrm{K},\end{array}$$

where $T_{{\mathrm{B}}_i}$ in Eq. (2) is an across-track row of T_B values and i is an index along the track, while j is an index across the track. n is the maximum across-track index of a row, i.e., for a complete row with valid data points for all 78 incidence angles, n=78.

Afterwards, single T_B outliers are removed if they do not satisfy Eq. (3):

$$\begin{array}{}\text{(3)}& |p_i-p_{i-\mathrm{1}}|+|p_{i+\mathrm{1}}-p_i|\le \mathrm{150}\mathrm{K},\end{array}$$

where p is a single-pixel T_B and i is an along-track index. The threshold of 150 K was selected manually after identifying erroneous single-pixel outliers in the data. The last filter in Eq. (4) removes neighboring T_B values which are locked to the same value, i.e., T_B values for which the following equation equals zero:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}|p_{i+\mathrm{6}}& -p_{i+\mathrm{5}}|+|p_{i+\mathrm{5}}-p_{i+\mathrm{4}}|+|p_{i+\mathrm{4}}-p_{i+\mathrm{3}}|+|p_{i+\mathrm{3}}\\ & -p_{i+\mathrm{2}}|+|p_{i+\mathrm{2}}-p_{i+\mathrm{1}}|+|p_{i+\mathrm{1}}-p_i|\ne \mathrm{0},\end{array}\end{array}$$

where p is again a single-pixel T_B with an along-track index i. The decision to compare seven consecutive T_B values was made based on qualitative experiments. Since the filter is used universally for all incidence angles, the search window varies but covers a minimum distance of 175 km.

The outer data points of the swath edges show significantly higher noise levels and a coarser resolution than the near-nadir data points (Veng, 2021). Therefore, after the filtering, we additionally remove the four outermost data points of the swath, corresponding to incidence angles between 57° and 64° on both sides. The new outer edge of the swath data is then at ∼56°, similar to that of modern microwave radiometers used for sea ice retrieval: 50–55° (AMSR, Meisner and Wentz, 2000; SMMR, Wentz, 1983; SSM/I, Wentz, 1997).

Before the filtering, the dataset contained 13 496 orbital data files; after filtering, there are 10 649 (∼79 %) good files left. A complete list of days where data are missing after filtering and no SIC could be retrieved is given in Appendix A3.

3 The radiative transfer model

The radiative transfer model (RTM) was developed specifically for atmospheric noise reduction and is comparable to the RTMs used in Andersen et al. (2006) and Tonboe et al. (2016) but with the addition that this ESMR RTM (Eq. 5) can be applied for different incidence angles over both ocean and ice. The RTM takes as input the atmospheric columnar water vapor V [mm or kg m⁻²], the 10 m wind speed W [m s⁻¹], the atmospheric columnar cloud liquid water L [mm or kg m⁻²], the sea surface temperature T_s [K], the ice emitting layer temperature T_i [K], the sea ice concentration c_ice [0–1], and the incidence angle [°]. In return, it simulates the top-of-the-atmosphere 19.35 GHz T_B at horizontal polarization:

$$\begin{array}{}\text{(5)}& T_{\mathrm{B}}=\mathrm{RTM}(V,W,L,T_{\mathrm{s}},T_{\mathrm{i}},c_{\mathrm{ice}},\mathit{\theta }).\end{array}$$

The RTM uses the atmospheric part of the model described in Wentz (1997) to compute the atmospheric emission, transmissivity, and reflectivity at the sea surface (open water) together with added modules for simulating the sea ice emissivity (Fig. 1) and open water reflectivity as a function of incidence angle.

Figure 1 The first-year ice emissivity used in the RTM.

For the sea ice emissivity, a look-up table is produced from a simulation using a combined sea ice thermodynamic and emission model during the Arctic winter on first-year ice. The thermodynamical and emission model setup and the simulations are described in Tonboe (2010). The emissivities as a function of incidence angle are shown in Fig. 1, and the look-up table is given in Table 1.

Table 1 Sea ice emissivity look-up table (θ is the incidence angle and E_ice is the first-year ice emissivity).

Sea water permittivity, which is used to estimate the calm sea reflectivity (Eq. 6) as a function of temperature, is computed using equation E64 (p. 2046) in Ulaby et al. (1986). The permittivity is almost invariant with respect to the water salinity at 19 GHz, and a constant value of 34 ppt is used here.

The calm sea (Fresnel) power reflection coefficient, r_h, as a function of the relative permittivity, ϵ, and the incidence angle, θ, for a lossy medium is computed using Eq. (1.52) in Schanda (1986), i.e.,

$$\begin{array}{}\text{(6)}& r_h\left(\mathit{\theta }\right)={\displaystyle \frac{(p-\cos (\mathit{\theta }){)}^{\mathrm{2}}+q^{\mathrm{2}}}{(p+\cos (\mathit{\theta }){)}^{\mathrm{2}}+q^{\mathrm{2}}}},\end{array}$$

where p and q are abbreviations for:

$$\begin{array}{}\text{(7)}& p={\displaystyle \frac{\mathrm{1}}{\sqrt{\left(\mathrm{2}\right)}}}{\left((\mathit{ϵ}{{}^{\prime }}_r-{\sin }^{\mathrm{2}}\mathit{\theta }{)}^{\mathrm{2}}+{{\mathit{ϵ}}^{\prime \prime }}_r^{\mathrm{2}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}+(\mathit{ϵ}{{}^{\prime }}_r-{\sin }^{\mathrm{2}}\mathit{\theta })\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\end{array}$$

and

$$\begin{array}{}\text{(8)}& q={\displaystyle \frac{\mathrm{1}}{\sqrt{\left(\mathrm{2}\right)}}}{\left((\mathit{ϵ}{{}^{\prime }}_r-{\sin }^{\mathrm{2}}\mathit{\theta }{)}^{\mathrm{2}}+{{\mathit{ϵ}}^{\prime \prime }}_r^{\mathrm{2}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}-(\mathit{ϵ}{{}^{\prime }}_r-{\sin }^{\mathrm{2}}\mathit{\theta })\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

The relative permittivity ${\mathit{ϵ}}_r=\mathit{ϵ}{{}^{\prime }}_r+\mathit{ϵ}{{}^{\prime \prime }}_{r}j$ of the water surface is a complex number. The calm sea emissivity, E₀, is then

$$\begin{array}{}\text{(9)}& E_{\mathrm{0}}=\mathrm{1}-r_h.\end{array}$$

The rough water surface emissivity component, E_W, is added to the calm sea emissivity, E₀, to produce the total sea water emissivity, E_water. Between ESMR incidence angles of 0 and 63°, the sensitivity of the 19.35 GHz horizontal polarization, E_W, to wind speed is an almost linear function $\left(\frac{\mathrm{\Delta }\left(E_{\mathrm{W}}{T}_{\mathrm{s}}\right)}{\mathrm{\Delta }W}=\mathrm{0.0094}\mathit{\theta }+\mathrm{0.3}\right)$ of the incidence angle θ (Meisner and Wentz, 2012), i.e.,

$$\begin{array}{}\text{(10)}& E_{\mathrm{W}}={\displaystyle \frac{W(\mathrm{0.0094}\mathit{\theta }+\mathrm{0.3})}{T_{\mathrm{s}}}},\end{array}$$

and

$$\begin{array}{}\text{(11)}& E_{\mathrm{water}}=E_{\mathrm{0}}+E_{\mathrm{W}},\end{array}$$

where θ is the incidence angle in degrees, W is the wind speed, and T_s is the sea surface temperature [K].

This combination of E₀ and E_W follows the procedure described in Wentz (1997).

The resulting brightness temperature is a linear combination of the sea water and sea ice emission weighted by the SIC following Andersen et al. (2006):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{T}_{\mathrm{B}}& =T_{\mathrm{BU}}+\mathit{\tau }\left(\right(\mathrm{1}-c_{\mathrm{ice}})E_{\mathrm{water}}\times T_{\mathrm{s}}+(\mathrm{1}-c_{\mathrm{ice}})\\ & (\mathrm{1}-E_{\mathrm{water}})(\mathrm{\Omega }{T}_{\mathrm{BD}}+\mathit{\tau }{T}_{\mathrm{BC}})+c_{\mathrm{ice}}{E}_{\mathrm{ice}}{T}_i+c_{\mathrm{ice}}\\ & (\mathrm{1}-E_{\mathrm{ice}})(T_{\mathrm{BD}}+\mathit{\tau }{T}_{\mathrm{BC}})),\end{array}\end{array}$$

where T_BU is the upwelling brightness temperature from the atmosphere, τ is the atmospheric transmissivity, E_water is the water surface emissivity, E_ice is the sea ice emissivity, Ω is the reflection reduction factor due to water surface roughness, T_BD is the downwelling brightness temperature, and T_BC is the cosmic background brightness temperature (2.7 K).

EMSR-simulated T_B values and emissivities have been compared with other T_B values simulated using other RTMs for AMSR (Meisner and Wentz, 2000), SMMR (Wentz, 1983), and SSM/I (Wentz, 1997) for a constant incidence angle of 55°, which is close to the incidence angles of other instruments and RTMs (53–55°). The comparison showed that the T_B values of the ESMR RTM are within the range (within approx. 2 K) of values obtained with the other models and therefore seem to be reasonable given the differences in instrument center frequencies and measurement geometries. Note that in the correction procedure, the difference between two simulated T_B values is used to minimize model biases. Even if the absolute values of the RTM-simulated T_B values were biased, this bias would be removed by taking the difference between two simulated T_B values, which is the only part used in the correction.

4 Derivation of the sea ice concentration

The RTM is an essential part of the SIC derivation for applying an atmospheric noise reduction. The following section presents the calculations of dynamical tie points and the SIC algorithm along with uncertainty estimations. Lastly, the post-processing, including the land-spillover method and data flag assignment, is described. A flow chart illustrating this processing chain can be found in the ESA CCI ESMR product user guide (PUG) (ESA CCI, 2022).

4.1 Tie points and geophysical noise reduction

Tie points are typical signatures of 100 % ice and open water (0 % ice) and are used in SIC algorithms as a reference for estimating the total ice fraction per satellite pixel c_ice. The aim of using dynamical tie points is to reduce SIC biases that may result from seasonal and interannual variations in T_B (Kongoli et al., 2011), instrument drift, and RTM and ERA5 biases. For example, Comiso and Zwally (1980) argue that the variations in the average open-water T_B values near the ice edge are affected primarily by variations in instrument calibration, and they describe the drop followed by a sharp increase seen in 1975 (see Fig. 2) as an instrument drift issue.

Our ESMR tie points are derived on a daily basis from the swath files. Regions of open water and high SICs are selected for each hemisphere, resulting in two regions of sea ice and two of open water for both hemispheres. The ERA5 SIC prior to 1979 is based on the HadISST2.0.0.0 dataset (Bell et al., 2021), which mainly utilizes digitized sea ice charts for this period (Rayner et al., 2003). The two main data sources are the Walsh dataset (Walsh, 1978; Walsh and Johnson, 1979; Walsh and Chapman, 2001) and National Ice Center (NIC) charts (Knight, 1984). The datasets also consist of several data types besides ice charts, e.g., ship observations and satellite data (both infrared and microwave observations, including data from ESMR). The selection of the four tie points is based on a criteria set for SIC from ERA5 – the distance from the ice edge, the observed brightness temperature, the latitude, and the sea surface temperature – as shown in Table 2. The distance from the ice edge criterion is imposed by putting a threshold on the mean SIC of a 5×5 grid point box; this should be larger than 80 % for ice tie points or less than 1 % for open water points. While computed daily, these are subsequently combined into 15 d running-mean tie points that are 7 d ahead and 7 d behind the processed date shown in Fig. 2. The 15 d averaging period is maintained even at the beginning and end of the dataset and when there are data gaps; i.e., if there are gaps, the averaging is done over the days available in the 15 d period (this is also the case in the beginning, where the first 7 d are missing). Figure 2 depicts the 15 d averaged tie points through time. It shows that the ice tie points follow a seasonal pattern, while the water tie points are relatively constant. The tie-point criteria from Table 2 ensure that each daily tie point is based on many observations, which results in stable tie points. The 15 d interval was chosen experimentally, so the T_B variations seem reasonable and one is still able to identify calibration issues such as jumps (e.g., that seen in 1976).

Table 2 Criteria for the four different tie points.

Figure 2 The 2-weekly tie points for Arctic and Antarctic ice and water after T_B correction. The boxes show the period May–July 1976, when there were obvious instrument calibration issues.

Figure 3 Mean atmospheric water vapor for all grid points included in the four tie points. Water vapor data are from ERA5 (Hersbach et al., 2020).

The per-grid-point T_B correction term ΔT_B,simulated is the difference between a simulated reference T_B obtained using mean values of the total column water vapor [kg m⁻²] in the atmosphere ($\stackrel{\mathrm{‾}}{V}$), 10 m wind speed [m s⁻¹] ($\stackrel{\mathrm{‾}}{W}$), total column cloud liquid water [kg m⁻²] in the atmosphere ($\stackrel{\mathrm{‾}}{L}$), sea surface temperature ($\stackrel{\mathrm{‾}}{T_{\mathrm{s}}}$), and ice emitting layer temperature ($\stackrel{\mathrm{‾}}{T_{\mathrm{i}}}$) as input to the RTM and a simulated T_B obtained using the actual ERA5 values (V, W, $\stackrel{\mathrm{‾}}{L}$, T_s, and T_i) for the grid point. T_B is not corrected for the cloud liquid water, L, so the mean L is input to both the reference and the actual simulation. ΔT_B,simulated can both be negative and positive, and the T_B values have a reduced sensitivity to the geophysical noise sources (V, W, T_s, T_i) after correction. The fact that the correction term is the difference between two RTM simulations minimizes the impact of biases in the model and the ERA5 data.

The correction term is added to the measured T_B, i.e.,

$$\begin{array}{}\text{(13)}& T_{\mathrm{B},\mathrm{corrected}}=T_{\mathrm{B},\mathrm{measured}}+\mathrm{\Delta }{T}_{\mathrm{B},\mathrm{simulated}}.\end{array}$$

Here,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{B},\mathrm{simulated}}& =\mathrm{RTM}(\stackrel{\mathrm{‾}}{V},\stackrel{\mathrm{‾}}{W},\stackrel{\mathrm{‾}}{L},\stackrel{\mathrm{‾}}{T_{\mathrm{s}}},\stackrel{\mathrm{‾}}{T_{\mathrm{i}}},c_{\mathrm{ice}},\mathit{\theta })\\ & -\mathrm{RTM}(V,W,\stackrel{\mathrm{‾}}{L},T_{\mathrm{s}},T_{\mathrm{i}},c_{\mathrm{ice}},\mathit{\theta }),\end{array}\end{array}$$

where c_ice=0 is the open water tie point and c_ice=1 the ice tie point. Following Svendsen et al. (1983), T_i is computed as

$$\begin{array}{}\text{(15)}& T_{\mathrm{i}}=\mathrm{0.4}\cdot T_{\mathrm{2}\mathrm{m}}+\mathrm{0.6}\cdot \mathrm{272},\end{array}$$

where T_2 m is the 2 m air temperature, which is taken from the ERA5 data. Horizontal bars above the variables indicate that they are daily mean values for the cluster of points selected for the tie point. The mean water vapor, $\stackrel{\mathrm{‾}}{V}$, at the tie point is shown in Fig. 3.

Figure 4 shows the correction term, ΔT_B,simulated, for 1 January 1974 over open water in the Northern and Southern hemispheres, respectively. The path length through the atmosphere is longest at high incidence angles and shortest near nadir, and thus the absolute value of the correction is largest at high incidence angles. For example, when the atmosphere is driest in the reference compared to the actual simulation, the ends of the corrected T_B turn negative, while they turn strongly positive when the reverse is true.

Figure 4 Difference between the T_B values before and after correction with a mean reference for open water tie points in the (a) Northern Hemisphere and (b) Southern Hemisphere.

The correction works best over open water areas, where it acts as only an atmospheric correction. The RTM appears to better simulate the relevant emission processes in the atmosphere, and the ERA5 data more accurately quantify the atmospheric noise sources. Over sea ice, geophysical noise sources are related to processes in the snow and ice profile (Tonboe et al., 2021) which are not characterized by the RTM, except for the emitting layer temperature T_i. T_i, which is used as input to the RTM, is estimated from the 2 m air temperature in the ERA5 data using Eq. (15). This is important because the ESMR is a single-channel instrument and thus the T_B and also the derived c_ice are sensitive to the emitting layer temperature.

The standard deviations of the brightness temperature for water points in both hemispheres before and after the correction are shown in Fig. 5.

Figure 5 Standard deviations of T_B values before and after correction in January 1974 for the (a) Northern Hemisphere and (b) Southern Hemisphere. The plots only show ocean points which fulfill the ERA5 SIC and sea surface temperature (SST) criteria described in Table 2.

4.2 The sea ice concentration (SIC) and its uncertainty

SIC (c_ice) is estimated using the measured brightness temperature (T_B,measured) and the open water (T_p,water) and ice (T_p,ice) tie points, i.e.,

$$\begin{array}{}\text{(16)}& c_{\mathrm{ice}}={\displaystyle \frac{T_{\mathrm{B},\mathrm{measured}}-T_{\mathrm{p},\mathrm{water}}}{T_{\mathrm{p},\mathrm{ice}}-T_{\mathrm{p},\mathrm{water}}}}.\end{array}$$

Because the RTM requires c_ice as input, c_ice is processed iteratively in two steps:

1. The c_ice is first estimated using uncorrected T_B values and tie points derived from uncorrected data. The c_ice estimate is truncated to the interval between 0 and 1 and an open water filter is applied, forcing all c_ice values of less than 0.15 to 0.

2. The c_ice estimate from step (1) is used in the RTM calculation (Eq. 5) together with ERA5 data for the geophysical noise reduction of the T_B values, and c_ice is then estimated again in a second iteration, this time using the corrected T_B values and corrected tie points. The mean values of $\stackrel{\mathrm{‾}}{V},\stackrel{\mathrm{‾}}{W},\stackrel{\mathrm{‾}}{L},\mathrm{\dots }$ used in the reference simulation are averages weighted by the c_ice obtained using the mean water and ice tie-point values respectively, i.e., c_ice is used as a ratio to mix the two tie-point values to create mean values of the numerical weather prediction (NWP) data for any sea ice concentration.

Iterations to update c_ice could continue in principle. However, tests show that updates are small after one iteration, and so we only iterate once (e.g., Lavergne et al., 2019).

The total SIC uncertainty is the combination of two components: (1) the algorithm uncertainty, which includes instrument noise and tie-point variability (geophysical noise), and (2) the resampling uncertainty, which is the uncertainty due to data resampling.

The algorithm uncertainty is the squared sum of three independent components following Parkinson et al. (1987):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}\mathit{\delta }{c}_{\mathrm{ice},\mathrm{algorithm}}& =\left({\left({\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{B}}}{T_{\mathrm{p},\mathrm{ice}}-T_{\mathrm{p},\mathrm{water}}}}\right)}^{\mathrm{2}}\right.\\ & +{\left({\displaystyle \frac{-\left(\mathrm{1}-c_{\mathrm{ice}}\right)\mathit{\delta }{T}_{\mathrm{p},\mathrm{water}}}{T_{\mathrm{p},\mathrm{ice}}-T_{\mathrm{p},\mathrm{water}}}}\right)}^{\mathrm{2}}\\ & {\left.+{\left({\displaystyle \frac{-c_{\mathrm{ice}}\mathit{\delta }{T}_{\mathrm{p},\mathrm{ice}}}{T_{\mathrm{p},\mathrm{ice}}-T_{\mathrm{p},\mathrm{water}}}}\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where the first term in Eq. (17) represents variations due to instrument noise, estimated to a brightness temperature error δT_B of 3 K (Parkinson et al., 1987).

Without the instrument noise term, which is already included in the two tie-point uncertainties, the second and third terms in Eq. (17) are used to compute the algorithm uncertainty, δc_{ice,algorithm}:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathit{\delta }{c}_{\mathrm{ice},\mathrm{algorithm}}& =\left({\left({\displaystyle \frac{-\left(\mathrm{1}-c_{\mathrm{ice}}\right)\mathit{\delta }{T}_{\mathrm{p},\mathrm{water}}}{T_{\mathrm{p},\mathrm{ice}}-T_{\mathrm{p},\mathrm{water}}}}\right)}^{\mathrm{2}}\right.\\ & +{\left.{\left({\displaystyle \frac{-c_{\mathrm{ice}}\mathit{\delta }{T}_{\mathrm{p},\mathrm{ice}}}{T_{\mathrm{p},\mathrm{ice}}-T_{\mathrm{p},\mathrm{water}}}}\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where δT_p,water is the water tie-point error (1 standard deviation of the daily tie point here) and δT_p,ice is the ice tie-point error (e.g., 1 standard deviation of the daily tie point). The water and ice tie-point errors are weighted by the SIC, and all three errors are normalized with the ice–water brightness temperature contrast and the 2-weekly tie points. The algorithm uncertainty is computed based on swath data.

The resampling uncertainty δc_{ice,re-sampling} is the maximum difference c_ice- minimum c_ice in a 3×3 pixel window. The resampling uncertainty is computed based on resampled data (e.g., Lavergne et al., 2019).

The total uncertainty is the squared sum of the algorithm and resampling uncertainties, i.e.,

$$\begin{array}{}\text{(19)}& \mathit{\delta }{c}_{\mathrm{ice},\mathrm{total}}={\left(\mathit{\delta }{c}_{\mathrm{ice},\mathrm{algorithm}}^{\mathrm{2}}+\mathit{\delta }{c}_{\text{ice,re-sampling}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

The two uncertainty components and the total uncertainty are included in the data file.

4.3 Land-spillover correction and post-processing

The land-spillover correction follows the procedure described in Markus and Cavalieri (2009). A 5 by 5 pixel neighborhood of the land mask (EASE2 version 2, by OSI-SAF) is analyzed to determine which coastal points should be corrected. The land mask is divided into two classes: land points, which are given a value of 90 % SIC, and open ocean points. If the difference between the original land mask and the mean mask calculated for the 5 by 5 window is larger than the previously estimated SIC (the RTM-corrected and resampled SIC), i.e., the SIC is smaller than the theoretical value of the land spillover, the SIC value is set to 0 % and the status flag variable of the dataset is raised to 8.

Additionally, a monthly climatology (also from OSI-SAF; the same version as for the land mask) is used to set SIC to 0 % and to mark open water points by a climatology boundary, which is indicated by a status flag value of 64. Afterwards, the land mask is also used to mark lakes and coastal areas with status flags 2 and 32, respectively. An overview of all status flag values is shown in Table 3.

Table 3 Description of values of the status flag variable of the dataset.

The results of the post-processing are included in the daily NetCDF files for the Northern and Southern hemispheres, respectively.

5 Results

A list of all output variables in the daily SIC files and a short description of them can be seen in Appendix A2. Examples of monthly means of the SIC and the mean uncertainty can be seen in Figs. 6 and 7.

It is worth noting that the coverage in Fig. 6 is complete and, because of the ESMR's wide swath width of 3100 km and its inclination, the North Pole is covered, in contrast to the satellite microwave radiometers that have followed the NIMBUS 5 ESMR, which have a “pole hole”. The area covered by multi-year ice in the central Arctic has a lower SIC than the first-year ice regions. This is a consequence of the one-channel SIC algorithm, which has inherent ambiguity in SIC, ice type, and emitting layer temperature variations. In Figs. 6b and 7b, it can be seen that the uncertainties are largest near the ice edges, as expected. This is due to the resampling uncertainty, which dominates near the ice edge, where the spatial variability of T_B is high. Coastal regions also show higher uncertainties for this reason, since the land-spillover correction is first applied to the SIC after the uncertainty estimations.

Figure 6 Monthly mean (a) SIC and (b) uncertainty for February 1973 in the Northern Hemisphere. Water areas with no ice and uncertainty due to the ocean climatology are displayed in dark blue.

Figure 7 Monthly mean (a) SIC and (b) uncertainty for October 1974 in the Southern Hemisphere. Water areas with no ice and uncertainty due to the ocean climatology are displayed in dark blue.

The SIC shows interesting sea ice features in the years 1972–1977. One such feature is the Odden ice tongue (Comiso et al., 2001) extending eastward from the East Greenland Current, which is visible in Fig. 6 (around 73° N, 0° E), while another feature is the Maud Rise Polynya (Jena et al., 2019), an open water area encircled by sea ice in the Southern Hemisphere, which can be seen in Fig. 7 (around 65° S, 0° E). Both examples were much larger in extent in the 1970s and occurred more frequently then than they do today (Comiso et al., 2001; Cheon and Gordon, 2019; Jena et al., 2019).

The daily coverage of valid data points that passed all filtering varied a lot through the ESMR's operating period. While there was nearly full coverage for the first few months, it got much worse after the summer of 1975, when the instrument only recorded data every second day. An example of the poor coverage is shown for May 1976 in Fig. 8.

Figure 8 Example of poor monthly coverage in May 1976 in the Northern Hemisphere for (a) the mean SIC and (b) the number of days with valid values. Water areas that showed full daily coverage due to the ocean climatology are displayed in dark blue to avoid misleading comparisons. The number of days with valid data is indicated by the colorbar.

Monthly averaged SICs are derived to compare our results against other datasets. Only months with 99 % coverage were used in the comparison, i.e., 99 % of all grid points are covered at least once per month. Monthly mean SIEs are calculated from the monthly SICs using the threshold c_ice>30 %. In Fig. 9, the ESMR dataset (orange line) is shown together with the OSI-SAF climate data record (CDR; blue line) for 1979–2022 (EUMETSAT, 2017a, b) and the sea ice extent derived from the NSIDC's NIMBUS 5 ESMR ice concentration product (green line) (Parkinson et al., 2004), using the same threshold for all products (c_ice>30 %).

The comparison shows comparable SIE levels around 1980. In general, our ESMR dataset has slightly higher monthly SIE values than the NSIDC's ESMR product, even though the seasonal pattern is the same.

The mean difference between our sea ice extent and that from NSIDC is 0.24 million km² in the Arctic and 0.59 million km² in the Antarctic for the whole dataset.

For the Northern Hemisphere, the SIE during the operational period of NIMBUS 5 ESMR (1972 to 1977) seems to have been similar in magnitude to the SIE during the operational period of NIMBUS 7 SMMR (from 1978 to 1987), with the ESMR's minimum extents being slightly higher than the SMMR's ones. In the Southern Hemisphere, the values for the second half of the 1970s seem to have been around the same magnitude as the largest SIE during the 2014/2015 season.

Figure 9 Monthly sea ice extent time series for the (a) Arctic and (b) Antarctic based on the c_ice>30 % threshold. The orange curve shows values in the ESMR dataset which have 99 % coverage of the hemisphere monthly, while the blue curve is based on the SIC products from the OSI-SAF (EUMETSAT, 2017a, b), where the c_ice>30 % threshold was applied as in EUMETSAT (2017c). The green line represents the NSIDC's ESMR SIC product (Parkinson et al., 2004).

6 Discussion

Comparisons between different sea ice products and the new ESMR dataset proved to be more difficult than initially expected, since not only the processing algorithms but also the land masks, map projections, and dataset grids differ. We were not able to find two independent SIC datasets for 1978 onwards and 1972–1977 which share exactly the same land mask.

Thus, it was decided at the beginning of the processing to use the same land mask as the OSI-420 product (1978 onwards) (EUMETSAT, 2020) for our ESMR dataset, i.e., a 25 km equal-area grid (EASE-2 version 2) land mask, to at least ensure a fair comparison between these two datasets. The NSIDC ESMR dataset (green line in Fig. 9) used a different land mask with a polar stereographic projection (Parkinson et al., 2004), but it was still compared to our ESMR dataset.

The difference in SIE is therefore also influenced by the different projections of the data; however, the area difference between the projections is relatively small (only a few thousand km²), so even a re-projection is expected to yield minimal differences compared to the differences caused by the use of different land masks. The comparison of different land masks is complicated by the varying sea ice extent, which exposes more or less land throughout its annual cycle and thus changes the number of grid points affected by the land mask.

The land mask land area differs between the OSI-SAF and NSIDC ESMR land masks. A comparison of land and ocean points between the NSIDC land mask and the OSI-SAF land mask showed a difference for the Northern Hemisphere of 460 000 km² (north of 60° N), where the NSIDC land mask has more land, while the difference is the opposite and much smaller in the Southern Hemisphere – only 79 000 km² (south of 60° S) – with the OSI-SAF's land mask having slightly more land. More land points in the land mask result in fewer available grid points for potential sea ice. The difference in Southern Hemisphere sea ice extent is significantly larger and opposite to the expected contrast from the land mask difference. In the Northern Hemisphere, it is not very clear how much of the SIE difference can be accounted for by the land mask or by algorithm differences. However, since the SIE differences vary a lot, and the differences in the Southern Hemisphere clearly cannot be explained by the land mask difference alone, it is likely that most of the SIE differences come from the algorithms and processing methods, such as the atmospheric correction and tie point calculation.

The blue curve in Fig. 9 is based on OSI-SAF's SIC products OSI-450 and OSI-430-b (EUMETSAT, 2017b, a). The shown extend corresponds to the SIE of the OSI-420 product (EUMETSAT, 2020) but with a 30 % instead of a 15 % sea ice threshold applied, matching the OSI-402-d sea ice extent product (EUMETSAT, 2017c).

Compared to the more common 15 %, the 30 % threshold was better suited to a comparison between different ESMR SIE datasets due to the relatively high noise level, which can be seen from the total uncertainty in Fig. 5. The uncertainty algorithm has been applied for easier data assessment and comparison to other datasets.

A large number of ESMR data are currently filtered out, and the 99 % threshold for inclusion in the monthly timeline in particular filters out the second half of the ESMR data, where large data gaps occurred, as seen in Fig. 8.

The filters worked as expected and removed erroneous T_B values from the raw data. A reprocessing of the data to rescue more of the currently filtered out data points from the 20 % of ESMR data files is planned. Two approaches will be used to improve the data selection during the next processing. One approach consists of making changes to the set of filters, i.e., performing an adjustment of the filter thresholds, testing different filters, and the possible incorporation of an incidence angle dependency into the data selection. The other approach is to recalibrate some of the erroneous data files that have shown some systematic offsets to rescue whole swaths.

To reduce the uncertainty caused by atmospheric noise, the brightness temperatures were corrected with an RTM that used several atmospheric parameters from NWP (ERA5) data, such as water vapor and wind, as input. This correction showed a consistent reduction of the standard deviation of the brightness temperatures for water points in both hemispheres, as can be seen from Fig. 5. The correction was less steady over ice surfaces, since the RTM did not describe all relevant processes related to the snow and ice, which are the main noise source over sea ice. Such processes include sea ice deformation, the creation of leads or ridges, and changes in the snow layer, e.g., changes in snow depth, snow density, and grain size but also melting and refreezing, which influences scattering and emission processes inside the layer. Atmospheric noise caused by, e.g., water vapor and cloud liquid water influences the T_B over sea ice, but its influence is smaller: only around one-third of the total noise (Tonboe et al., 2021). By correcting for atmospheric effects with ERA5 data, we might have introduced some noise into the angular dependency of the SIC due to the use of an incident-angle-dependent emissivity in the RTM (Fig. 1 and Table 1).

To avoid biases from the RTM and the NWP data, dynamical tie points have been used, which also calibrate the algorithm for seasonal variations and instrument drift. However, we currently use mean tie points that are independent of the incident angle. Therefore, a possible improvement to a future version of the dataset might be accomplished by using angle-dependent tie points instead.

Even after filtering the data for obvious errors, it is clear that there are still issues with the absolute calibration of the instrument (Comiso and Zwally, 1980). For example, after the hot-load anomaly in 1973, the ocean T_B in the Southern Hemisphere is several kelvin below the T_B level before the anomaly, and in 1976 there is a dip in May and June followed by a sharp increase in T_B (Zwally et al., 1983). Low-frequency (with timescales ≥ days) T_B variations and regional variations at hemispherical scales are compensated for by the dynamical SIC algorithm tie points (Tonboe et al., 2016).

In spite of data gaps and calibration issues, the experimental NIMBUS satellite program was very successful. Applying modern processing methodologies, including dynamical tie points and atmospheric noise reduction of the T_B values, reduces the noise over both ice and open water consistently. This newly processed ESMR sea ice dataset extends the existing sea ice climate data record (CDR) with an important period from the 1970s. This extension of the SIE record contributes to the United Nations Sustainable Development Goals (SDGs) related to climate change by providing more observations for longer-term assessments of Arctic and Antarctic sea ice changes.

7 Data availability

The newly processed ESMR data have been released through the ESA CCI Open Data Portal: https://climate.esa.int/en/odp/#/project/sea-ice (last access: 8 March 2024) and at DOI: https://doi.org/10.5285/34a15b96f1134d9e95b9e486d74e49cf (Tonboe et al., 2023).

8 Conclusions

In this paper, we have presented a new SIC dataset covering the period 1972–1977, obtained using the ESMR data from the NIMBUS 5 satellite. The dataset consists of resampled daily NetCDF files for the Northern and Southern hemispheres. The SIC, associated uncertainties, and processing flags are included in the dataset. The uncertainties follow the same principles as those of the EUMETSAT SIC CDR, including both algorithm and resampling uncertainties. The choice to use the same land mask, spatial grid, and projection as for EUMETSAT's SIC CDR make comparisons between the time periods easier.

A comparison to NSIDC's ESMR SIC product and the OSI-SAF CDR showed that the seasonal pattern is very similar to NSIDC's ESMR SIC product, but our product shows systematically larger SIE values which can not be explained by differences between land masks alone. For the Northern Hemisphere, our SIE values match the levels in the 1980s from the OSI-SAF CDR with the same land mask, while values in the Southern Hemisphere in the 1970s are larger than those in the 1980s.

Compatibility with the EUMETSAT's SIC CDR was achieved by using a similar processing chain. The processing included atmospheric noise reduction with the use of an RTM and the ERA5 atmospheric data, which lowered the standard deviation of the T_Bs consistently. Additionally, dynamical tie points were used to avoid biases from the RTM and NWP data as well as to adjust for seasonal variability and instrument biases. To ensure better data assessment in analysis and models and easier comparison to other datasets, temporally and regionally varying uncertainty estimates have been included in our ESMR dataset.

Appendix A

Table A1 The data variables in the NetCDF processing input file and a description of each variable.

Table A2 The output variables stored in the daily NetCDF files and a description of each variable.

Table A3 Missing dates.