PIOMAS-20C, an Arctic sea ice reconstruction for 1901–2010, is produced by forcing the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) with ERA-20C atmospheric data. ERA-20C performance over Arctic sea ice is assessed by comparisons with measurements and data from other reanalyses. ERA-20C performs similarly with respect to the annual cycle of downwelling radiation, air temperature, and wind speed compared to reanalyses with more extensive data assimilation such as ERA-Interim and MERRA. PIOMAS-20C sea ice thickness and volume are then compared with in situ and aircraft remote sensing observations for the period of ~1950–2010. Error statistics are similar to those for PIOMAS. We compare the magnitude and patterns of sea ice variability between the first half of the twentieth century (1901–40) and the more recent period (1980–2010), both marked by sea ice decline in the Arctic. The first period contains the so-called early-twentieth-century warming (ETCW; ~1920–40) during which the Atlantic sector saw a significant decline in sea ice volume, but the Pacific sector did not. The sea ice decline over the 1979–2010 period is pan-Arctic and 6 times larger than the net decline during the 1901–40 period. Sea ice volume trends reconstructed solely from surface temperature anomalies are smaller than PIOMAS-20C, suggesting that mechanisms other than warming, such as changes in ice motion and deformation, played a significant role in determining sea ice volume trends during both periods.
Changes in Arctic sea ice are an important fingerprint of natural and anthropogenic climate change. The dominant signal in sea ice variability over the satellite era (1979–present) is the reduction of sea ice extent, area, and thickness. While the first two characteristics are well measured from satellites, a basinwide record of sea ice thickness and volume is not available from direct measurements over the same period. Instead, this record is either pieced together from a variety of in situ measurements and remote observations from an array of platforms (e.g., upward-looking sonar, aircraft-based sensors) (Kwok 2018; Kwok and Cunningham 2015; Kwok et al. 2009; Laxon et al. 2013; Lindsay and Schweiger 2015; Rothrock et al. 1999) or reconstructed by driving an ice–ocean model with atmospheric reanalysis data while assimilating any sea ice data that are suitable. This constitutes a model-based sea ice reconstruction, which provides sea ice thickness and volume, as well as other ice and ocean variables, at daily or monthly time scales, typically from 1979 to the present. (e.g., Chevallier et al. 2017; Fučkar et al. 2015; Johnson et al. 2007; Kauker et al. 2008; Kauker et al. 2009; Schweiger et al. 2011; Tietsche et al. 2014). Such reconstructions are also sometimes referred to as sea ice reanalyses.
A widely used sea ice reconstruction of this type is the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) (Schweiger et al. 2011; Zhang and Rothrock 2003). Estimated sea ice thickness and volume uncertainties from PIOMAS are of similar magnitude to those currently observable from satellite (Labe et al. 2018; Laxon et al. 2013; Schweiger et al. 2011; Stroeve et al. 2014; Tilling et al. 2015; Wang et al. 2016). Sea ice volume variability since 1979 is dominated by a continuous decline that is in large part attributable to anthropogenic global warming (Min et al. 2008; Notz and Marotzke 2012), although contributions from internal variability are not negligible (Ding et al. 2019; Jahn et al. 2016; Notz and Stroeve 2016; Screen and Deser 2019; Swart et al. 2015; Winton 2011) and may account for as much as 50% of the variability in fall sea ice extent (Ding et al. 2017; Kay et al. 2011). The impact from a substantial contribution from internal variability at decadal time scales is difficult to assess, both due to the relative scarcity of historical sea ice data and because it is difficult to determine whether a mismatch between observed and modeled sea ice trends is due to systematic problems with the model or to internal variability (Ding et al. 2019; Winton 2011). Therefore, an extended record of sea ice parameters, including volume and thickness, is desirable because it provides the opportunity to assess both forced and internal variability at longer time scales.
One important and longstanding question is how the present-day loss of Arctic sea ice compares to a notable historical event now known as the early-twentieth-century warming (ETCW) episode that occurred roughly between 1920 and 1940 (Bengtsson et al. 2004; Brooks 1938; Hegerl et al. 2018; Scherhag 1937; Wood and Overland 2010). Sea ice datasets have been pieced together from ship and shore station reports reaching back to the 1850s (Chapman and Walsh 1991; Mahoney et al. 2011; Rayner et al. 2003; Titchner and Rayner 2014; Walsh et al. 2016; Zakharov 1997), but they provide no information about thickness and total volume changes needed to fully assess the impact of the ETCW on Arctic sea ice. Kauker et al. (2008) developed a model-based reconstruction of sea ice for the Arctic spanning from 1900 to 1997. Their approach relied on statistically reconstructed atmospheric forcing fields for the North Atlantic–Arctic Ocean–Sea Ice Model (NAOSIM) from gridded and in situ observations, which they then used to model sea ice extent for comparison with historical observations, particularly for the ETCW period in the Arctic. Their analysis focused on sea ice extent and concentration, not on thickness and volume. Other reconstructions have used the correlation between surface air temperatures and sea ice extent to reconstruct sea ice extent for summer (Alekseev et al. 2016) or for all seasons (Connolly et al. 2017). However, reconstructions based on temperature alone cannot fully represent the role that both thermodynamics and dynamics play in modulating sea ice thickness and volume (e.g., Koberle and Gerdes 2003; Rothrock and Zhang 2005).
The advent of extended atmospheric reanalyses (Compo et al. 2011; Hersbach et al. 2015; Poli et al. 2016) offers an alternative approach to statistical reconstructions by providing physically constrained atmospheric reanalysis fields that can be used to directly drive an ice–ocean modeling and assimilation system. Here we present first results from such a sea ice reconstruction, utilizing the European Centre for Medium-Range Forecasts (ECMWF) twentieth-century reanalysis (ERA-20C) to force a PIOMAS integration from 1901 to 2010. We refer to the output as the PIOMAS-20C sea ice reconstruction, which provides a record of Arctic sea ice thickness and volume variation over the entire twentieth century and into the early twenty-first century.
2. Approach and organization
We first introduce the ice–ocean model and the atmospheric forcing data used for driving the model (section 3). We then evaluate individual parameters in the atmospheric forcing fields of ERA-20C relative to measurements and other reanalysis datasets (section 4). Although ERA-20C has been used to examine atmospheric variability in the Arctic (Belleflamme et al. 2015; Wegmann et al. 2017; Wegmann et al. 2018), to our knowledge an examination of ERA-20C for the purpose of providing atmospheric forcing for sea ice models has not been conducted. Since reconstruction results depend on this forcing, this is a critical step. This is particularly important since direct model output validation data (e.g., thickness) for the period before the satellite record are sparse. The following data sources are employed. For radiation measurements we use data from the Surface Heat Budget of the Arctic Experiment (SHEBA), which to date still provide the most comprehensive dataset of the ice-covered Arctic Ocean covering the annual cycle. Temperature and wind data are compared with former Soviet Union North Polar drifting stations (hereafter NP), which provide a data record from 1957 through 1990 (Lindsay 1998) and begins before satellite information substantially influences reanalyses.
We then demonstrate the predictive capabilities (hindcast) for sea ice thickness and other sea ice parameters using ERA-20C for a calibration and validation period for which extensive validation data exists (section 5). PIOMAS-20C is then constructed for the entire period and compared to existing ice thickness and sea ice drift records with particular attention on earlier periods when data are usually limited and sparse in time and space (section 6). For example we use data from the Norwegian North Polar Expedition with the Maud for the 1922–24 period, including air temperature, wind, sea ice drift, and sea ice thickness (Sverdrup 1927). We then examine the resulting PIOMAS-20C reconstruction with a focus on sea ice thickness and volume and compare them to direct observations as well as the standard PIOMAS sea ice reconstruction (section 7a and 7b). We then compare ice thickness variability during the ETCW period with the more recent sea ice thickness decline (section 7c). We also utilize an entirely new dataset of sea ice information extracted from the logbooks of U.S. government vessels that operated in the Arctic in the early twentieth century (section 7d). This information is compared with model-derived sea ice information for the early twentieth century in the Pacific sector. Finally we conduct an uncertainty analysis of our sea ice volume trend estimates and to test the robustness of our conclusions (section 7e).
3. Methods and data
a. Ice ocean model
The numerical modeling system underlying PIOMAS and PIOMAS-20C consists of coupled sea ice and ocean modeling components. The sea ice model is a multicategory thickness and enthalpy distribution sea ice model that employs a teardrop viscous plastic rheology (Zhang and Rothrock 2005), a mechanical redistribution function for ice ridging (Hibler 1980; Thorndike et al. 1975), and a LSR (line successive relaxation) dynamics solver (Zhang and Hibler 1997). The model features 12 ice thickness categories covering ice up to 28 m thick. The sea ice model is coupled with the Parallel Ocean Program (POP) model developed at the Los Alamos National Laboratory. The PIOMAS model domain is based on a curvilinear grid with the north pole of the grid displaced into Greenland. It covers the area north of 49°N with an average grid-cell size of ~40 km and is one-way nested into a similar, but global, ice–ocean model (Zhang 2005). Note that the domain configuration excludes some areas in the Sea of Okhotsk, the Gulf of St. Lawrence, and the Labrador Sea that are covered by sea ice during some winter months. Therefore, comparisons with total volume from other sources need to take account the different domain sizes.
PIOMAS is capable of assimilating ice concentration data using an optimal interpolation approach (Lindsay and Zhang 2006). The assimilation procedure in PIOMAS-20C was slightly changed from PIOMAS so that satellite ice concentrations are assimilated only near the ice edge (defined as 0.15 ice concentration). This means that the assimilation is allowed only in the areas where either model or satellite ice concentration is at or below 0.15. In other words, no assimilation is conducted in the areas where both model and satellite ice concentrations are above 0.15. If the observed ice edge exceeds the model ice edge, then sea ice is added to the thinnest sea ice thickness category. If the model ice edge exceeds observations, excess ice is removed in all thickness categories proportionally. In the areas where observed ice concentration is above 0.15 and model ice concentration is below 0.15, the surface ocean temperature is set to the freezing point. The ice-edge assimilation approach forces the simulated ice edge close to observations, while preventing satellite-derived ice concentrations, which can be biased low during the summer (e.g., Ivanova et al. 2015), from inaccurately correcting model ice concentrations in the interior of the ice pack. In addition, sea ice concentration data prior to the routine satellite observations typically are based on observed ice-edge information with concentrations derived from climatological gradients based on satellite data. Assimilating the ice edge rather than concentrations therefore prioritizes the actual observations. PIOMAS is also capable of assimilating observations of sea surface temperature (SST) following Manda et al. (2005); however, SST assimilation is not used for the PIOMAS-20C project since we found it to add no additional skill.
The PIOMAS framework has undergone substantial validation (Labe et al. 2018; Laliberté et al. 2018; Laxon et al. 2013; Schweiger et al. 2011; Stroeve et al. 2014; Wang et al. 2016) and has been shown to simulate sea ice thickness with error statistics similar to the uncertainty of the observations. Although the default PIOMAS reconstruction is driven with NCEP–NCAR reanalysis (version 1; hereafter NCEP-R1) data, PIOMAS has been successfully integrated using atmospheric forcing data from different atmospheric reanalysis projects (e.g., CFSR, CSFv2, MERRA-1, ERA-Interim) (Lindsay et al. 2014).
b. Atmospheric forcing data, comparison with observations, and other reanalysis data
PIOMAS and PIOMAS-20C are driven with atmospheric forcing data consisting of downwelling longwave, shortwave fluxes, 10-m wind speed, 2-m surface air temperature and humidity at daily time resolution. Precipitation minus evaporation (P − E) is calculated from precipitation and latent heat fluxes provided by the reanalysis model and specified at monthly time resolution in order to allow the calculation of snow depth over sea ice and input of freshwater into the ocean. River inflow into the model domain is specified from climatology (Hibler and Bryan 1987). Turbulent and momentum transfer is calculated using a surface layer model (Briegleb et al. 2004) that is part of the PIOMAS framework. For PIOMAS-20C, atmospheric forcing data come from ERA-20C (Hersbach et al. 2015; Poli et al. 2016). ERA-20C provides a global atmospheric reanalysis for the period from 1900 to 2010. It relies on the Integrated Forecast System (IFS cy38r1). The reanalysis model runs at T159 with an approximate spatial resolution of 125 km. Prior to running the model, atmospheric forcing data are interpolated to the model grid using bilinear interpolation. Sea surface temperature and sea ice concentrations are prescribed using HadISST, version 184.108.40.206 (Titchner and Rayner 2014). Other external forcing data (aerosols, ozone, greenhouse gases) are specified according to the CMIP5 protocol. ERA-20C only assimilates in situ observations of surface level pressure and marine 10-m wind speeds. For consistency PIOMAS-20C assimilates the same ice concentration data that were used in the ERA-20C project.
For comparison and uncertainty assessment we also use the NOAA-CIRES Twentieth Century Reanalysis, version 20CRv2c (Compo et al. 2011), and the new CERA-20C (Laloyaux et al. 2016) from ECMWF. CERA-20C is a 110-yr reanalysis based on a coupled model and assimilation scheme using only marine winds and surface pressures as in ERA-20C. A 10-member ensemble is available for CERA-20 based on perturbed initial SST and sea ice conditions. For the uncertainty assessment in both cases only ensemble means are used. An examination of surface air temperature variability among CERA-20C members over the Arctic north of 60°N shows only small differences and no trend differences between ensemble members, suggesting that the assimilated data significantly constrain atmospheric temperature variability.
4. Validation of atmospheric forcing data
When embarking on this project, we had two choices for atmospheric reanalysis fields covering the twentieth century: the NOAA-CIRES Twentieth Century Reanalysis (version 20CRv2c) (Compo et al. 2011) and ERA-20C. The new CERA-20C only became available recently and we therefore only used it to assess uncertainties. To help decide which reanalysis was more suitable, or whether it would be useful to generate parallel sea ice simulations in order to assess uncertainty, we examined how both products represented atmospheric variables over the sea ice covered area of the Arctic Ocean. We used data from the SHEBA experiment, which includes surface measurement of winds, temperature, downwelling longwave and shortwave radiation (Persson et al. 2002) for nearly a complete annual cycle from 1997 through 1998. In addition, we used surface air temperature and wind data from the former Soviet Union NP drifting ice-station data (Lindsay 1998). This dataset provides nearly continuous measurements from one to two simultaneous drifting ice stations from 1968 to 1991. Data from the Maud expedition for the 1922–24 period provided information about sea ice drift and thickness, as well sea level pressure and surface temperature.
During our investigation, we discovered several issues with 20CRv2c in the Arctic. Although some variables of the 20CRv2c show excellent fidelity in comparisons with validation data and significant improvement in its Arctic performance over the previous version (Lindsay et al. 2014), surface air temperatures in winter are still strongly biased warm. Moreover, we identified spurious surface air temperature variability due to the way the thermodynamic sea ice model in 20CRv2c is initialized as part of the processing streams that break the processing into 5-yr segments. Because the sea ice thickness in the 20CRv2c framework is interactive, changes in the surface energy balance lead to a change in sea ice thickness until the next initialization (beginning of the 5-yr stream). Differences between initialized and resulting sea ice thickness therefore introduce a spurious 5-yr cycle into some of the 20CRv2c variables (Fig. S1 in the online supplemental material). Experiments to force our ice ocean model with 20CRv2c data did not generate credible sea ice fields largely due to excessively warm winter temperatures (Fig. 1a), more than 5°C in November, which yielded unrealistically thin ice that did not survive the summer. Model tuning strategies that have yielded realistic simulations with other forcing datasets (Lindsay et al. 2014) were not successful with 20CRv2c. We therefore abandoned this effort and focused on ERA-20C. 20CRv2c is currently undergoing reprocessing which will hopefully mitigate the above listed issues (G. Compo 2019, personal communication). Atmospheric forcing validation focuses on ERA-20C, but where appropriate 20CRv2c meteorological variables are shown for comparison. For all intercomparisons, reanalysis data were reduced to a daily mean and interpolated to a regular pole-centered equal area grid with 40-km resolution. Data for comparison with station observations were then interpolated linearly to the corresponding daily mean position of the surface station. To remove annual cycles that strongly affect correlations, we employed anomaly correlations and a daily differencing scheme.
a. Downwelling radiative fluxes
Figures 1c and 1e show the annual cycle for downwelling longwave and shortwave fluxes from five reanalysis datasets and from the measurements made at the SHEBA Atmospheric Surface Flux Group (ASFG) tower (Persson et al. 2002). ERA-20C differences with measurements are on the order of 20 W m−2 and on par with other reanalysis products except for NCEP-R1, which shows previously documented (e.g., Makshtas et al. 2007; Schweiger et al. 2008; Serreze et al. 1998) large compensating discrepancies with measurements in both long and shortwave fluxes are due to underestimated summer cloud cover in NCEP-R1. Daily differences (Figs. 1d and 1f), which remove the annual cycle, show that considerable amounts of synoptic variability are captured in ERA-20C. Seven-day smoothed daily differences are correlated at r = 0.60 for both long and shortwave fluxes. In terms of the annual cycle, downwelling radiative fluxes, both long-term reanalysis datasets (20CRv2 and ERA-20C) show equivalent performance to the reanalysis projects that assimilate all available data. However, correlations of daily differences are slightly higher for ERA-Interim (hereafter ERA-INT), NCEP-R1, and MERRA (not shown).
b. Surface air temperature
The annual cycle of surface air temperatures from ERA-20C is very close to observations (Fig. 1a). Notably, winter temperatures tend to be lower and closer to observations than the other reanalysis datasets except NCEP-R1. ERA-20C, using a similar atmospheric model and sea ice representation as ERA-INT, is cooler by several degrees during the fall and winter months and closer to observations than the ERA-INT. NCEP-R1, which does not assimilate observed surface air temperatures over Arctic sea ice, similarly has colder surface air temperatures than observations. The 20CRv2c temperatures are substantially too warm during the winter and fall months with biases greater than 5°C in November. The ERA-20C surface temperatures are highly correlated (0.97) with SHEBA observations with an RMS error of 3°C and very little overall bias (Fig. S2). Removing the annual cycle by taking daily differences (differences from one day to the next) yields correlations of r = 0.40 and r = 0.46 (7-day smoothing) (Fig. 1b). For the NCEP-R1 the RMS error is slightly smaller than for the ERA-20C. For comparison, 20CRv2c surface air temperatures have an RMS error of 4°C, correlation of 0.98, but substantial positive biases in winter. Comparison with NP station data shows a similar relationship between observations and reanalysis datasets as found at the SHEBA site, although ERA-20C values are warmer than observations by about 2°C and warmer than NCEP-R1 from January through March (Fig. S3). ERA-20C surface air temperature validation statistics for a much more data-sparse period can be computed using measurements from the Maud expedition for the period 1922–24. RMS errors (5.2°C) during that time increase markedly from later periods, but ERA-20C temperatures remain highly correlated with observations (0.92) and show relatively little bias (Fig. 2). Daily differences and monthly anomalies clearly indicate that ERA-20C has skill beyond capturing the annual cycle with monthly anomaly correlations at r = 0.72.
c. Wind speed
Sea ice dynamics and thermodynamics are directly influenced by surface winds, and thus an accurate representation is critical to achieve realistic sea ice thickness variability. To assess ERA-20C we compare 10-m wind speeds with equivalent observations from NP stations. Figure 3 shows the comparison of the annual cycle of mean monthly wind speed for ERA-20C and other reanalysis datasets and observations from the NP stations. The ERA-20C wind speeds are close to observed wind speeds with a slight positive bias and slightly higher than the ERA-INT wind speeds. Daily ERA-20C wind speed RMS error is 2.5 m s−1 with a correlation of 0.50. This compares to RMS errors of 1.3 m s−1 and correlation of 0.85 for ERA-INT, which performs best relative to the NP stations. The NCEP-R1 wind speeds are consistently lower than observed and typically lower than other reanalysis datasets. Using only NP stations prior to 1979 (when satellite data begin to be assimilated in modern reanalyses) yields very similar errors. We therefore consider those values to be representative for ERA-20C prior to the satellite period. NP station results are similar to comparisons with measurements at the ASFG tower during SHEBA (not shown), although less noisy due to the greater number of observations. Interestingly, 20CRv2c shows substantially better error statistics for SHEBA than ERA-20C with RMS errors of 1.8 m s−1 and correlations of 0.76 versus RMS and correlation values of 2.3 m s−1 and 0.54 for ERA-20C. Wind speeds are also available for the Maud expedition. They were originally recorded in wind force and subsequently converted to wind speed and are noisy in part due quantization errors. Applying a 7-day running mean yields correlations of r = 0.61 and an RMS error of 2.0 m s−1 (Fig. 4). This indicates that ERA-20C is capable of capturing a considerable part of wind speed variability at weekly time scales during the early twentieth century when very little data was available to constrain the reanalysis.
5. PIOMAS-20C calibration
Following the assessment and selection of ERA-20C to supply atmospheric forcing for PIOMAS-20C integration, we perform a two-step model calibration. The calibration and validation approach follows the strategy developed for the standard PIOMAS integration (Schweiger et al. 2011) and ice–ocean model experiments using alternate atmospheric forcing fields (Lindsay et al. 2014). During the development of PIOMAS we have previously tuned parameters such as water drag turning angle, ice strength coefficient, and magnitude of tensile stress. We have found that the two parameters, melting ice albedo and surface roughness length, are most effective to reduce PIOMAS bias and RMS errors relative to observations when adjusting to different atmospheric forcing datasets. Therefore only those two were adjusted for the calibration of PIOMAS-20C.
a. Ice drift calibration
PIOMAS-20C ice speed is tuned to minimize differences with a set of drifting buoys from the International Arctic Buoy Program (IABP) following Zhang et al. (2012). This drift validation dataset consists of daily averaged drift speeds from all available IABP buoys from 1979 to 2010 that have been screened for consistency. Differences in mean daily drift speed between observations and model are minimized by adjusting the aerodynamic surface roughness in the surface layer model that couples the atmospheric forcing fields to the sea ice surface. A roughness length value of 0.0008 m is selected for PIOMAS-20C after calibration using buoy drift data, which results in a mean model ice speed bias of −0.005 m s−1 and model–data correlation of 0.71 over the period of 1979–2010.
b. Mean ice thickness calibration
Surface albedo for melting sea ice is adjusted to minimize mean model ice thickness bias against in situ measurements available from the Unified Sea Ice Thickness Climate Data Record (ThickCDR) (Lindsay 2010; Lindsay and Schweiger 2015). ThickCDR contains ice thickness observations from in situ, submarine, airborne, and satellite remote sensing platforms. Only a subset of measurements (1975–2009) and only upward-looking sonar (ULS) and airborne electromagnetic induction (EM) measurements are used for model calibration. Mean ice thickness differences for N = 3101 model–observation pairs is minimized using a manual process involving sequential test integrations with the model. The mean monthly difference in (effective) ice thickness for PIOMAS-20C for the calibration dataset is 0.11 m. This compares favorably with a mean ice thickness difference of −0.06 m for PIOMAS. The albedo for melting ice is set to 0.65 for PIOMAS-20C.
6. PIOMAS-20C assessment
a. Sea ice thickness
PIOMAS-20C is then compared to all available observations from the ThickCDR (version 20170601), which contains sea ice thickness, sea ice draft, and sea ice + snow thickness records from 1948 through the present. Ice thickness measurements in the ThickCDR usually reflects the “effective” ice thickness or draft for a 50-km spatial average including the open water category. PIOMAS-20C snow water equivalent is used to convert ThickCDR observations of draft to sea ice thickness. For EM-based measurements, ThickCDR provides a combined sea ice plus snow thickness. Model snow thickness is therefore subtracted from EM measurements to compute sea ice thickness (Schweiger et al. 2011). ThickCDR is augmented with ice thickness observations made during the Maud expedition from 1922 to 1924 (Sverdrup 1927).
To adjust measurements from the Maud expedition and those from the stations in the ThickCDR that are based on individual “ice-only” measurements (data points in ThickCDR identified as “Canadian coastal”), model ice concentration is used to convert the observed ice thickness to effective ice thickness. One of the Canadian coastal stations located on Isachsen Island provides two separate, partially overlapping, records with large differences in mean ice thickness. Both station records are highly correlated, allowing an adjustment via a linear regression (Fig. S5). The adjustment does not affect overall validation statistics significantly but decreases the mean difference between PIOMAS-20C and observations between 1948 and 1954 from 0.6 to 0.4 m.
Figure 5 shows the time series of a comparison with PIOMAS-20C and corresponding ThickCDR observations. The largest differences in the time series occur from 1948 to 1954 when the observations come exclusively from a few Canadian coastal landfast ice stations. Some of these stations show a very good agreement with PIOMAS-20C, while others show poor correlation and biases; however, overall PIOMAS-20C clearly captures a portion of the ice thickness variability at these stations at monthly and interannual time scales with annual anomalies correlated at r = 0.52 (Fig. 6). Rather than selecting stations based on their agreement with PIOMAS-20C or some other criteria, we elected to include all available stations. Although PIOMAS-20C does not explicitly simulate landfast ice through grounding, the teardrop plastic rheology used in PIOMAS-20C allows biaxial tensile stress (Zhang and Rothrock 2005), which makes it easier for ice to “stick” to the coast under wind and current forcing, thus behaving like landfast ice (Lemieux et al. 2016). A recent evaluation of landfast ice representation in sea ice models found good agreement between PIOMAS and landfast ice observations (Laliberté et al. 2018).
Mean annual ice thickness from PIOMAS-20C is within the range of, or close to, the standard errors of the annual mean thickness of observations. Figure 7 shows observed and modeled mean annual ice thickness with model data sampled at observation times and locations. PIOMAS-20C ice thickness on average is thicker than observations by 0.15 m with an RMS difference of 0.26 m and a correlation of 0.85. Using monthly statistics, the mean ice thickness error is 0.09 m, the correlation drops to 0.67, and the RMS error increases to 0.88 m (excluding satellite data). Monthly statistics are substantially influenced by the Canadian coastal station record, which provides landfast ice thickness. Excluding these measurements improves the correlation to 0.72, reduces the mean error to 0.04 m, and slightly increases the RMS error to 0.91 m. However, without the Canadian coastal stations the observational record does not start until 1960 when ice thickness measurements from the first submarine cruises become available in ThickCDR. PIOMAS-20C validation statistics for the 1979–2010 period that overlaps with the standard PIOMAS product are nearly identical with correlations of 0.76, a mean error of 0.06 m, and RMS error of 0.81 m (not shown). The previously documented tendency of PIOMAS to overestimate areas of thin ice and underestimate areas of thick ice (Schweiger et al. 2011) is slightly more pronounced in PIOMAS-20C. This tendency is apparent in the spatial patterns of ice thickness with sea ice being too thick in the Beaufort Sea area and too thin elsewhere (Fig. 8). This spatial bias is common among current sea ice models (Chevallier et al. 2017; Johnson et al. 2012; Uotila et al. 2018). There is considerable uncertainty regarding the magnitude and exact spatial pattern of this bias since sea ice thickness retrievals from ICESat are affected by biases in the used snow depth, which has substantial uncertainties (Blanchard-Wrigglesworth et al. 2018). Other isolated sea ice thickness measurements not yet integrated into the ThickCDR also support that PIOMAS-20C sea ice thickness simulations capture mean sea ice conditions and variability (Fig. S6).
b. Sea ice drift
Direct observations of meteorological and sea ice variables are sparse, particularly prior to the 1950s, when regular measurements from NP and other international drifting stations began. An indirect validation is possible by using the model simulation to compare daily drift speed and direction at ship/station locations with daily positions logged at the time. We examine the records of two multiyear stations from the presatellite period, the drift of the United States Arctic Research Laboratory Ice Station 2 (ARLIS-2; 1961–65) and the drift of the Maud expedition between 1922 and 1924. Daily drift speeds from the ARLIS-2 station for 1961–65 (Fig. 9) show an excellent match during some periods (1961–62) but are a bit noisier during other times. Overall, daily drift speeds are correlated at 0.4 with and RMS error of 2.4 km day−1. Applying a 7-day smoothing improves correlations to 0.7. Drift agreement for the earlier Maud expedition is lower, with daily correlations of 0.37, or 0.53 after a 7-day smoothing is applied, likely reflecting the less well-constrained atmospheric forcing and possible errors in the ship’s navigation during that time (Fig. 10). Even though the timing and intensity of individual drift events appear to be less well simulated during this time, the distribution of drift speed and direction of simulated and observed drift show strong similarities. Figure 11 shows simulated and observed daily drift speed and direction distributions in the form of drift rosettes. Both show very similar characteristics in the distribution of drift speed and direction. This comparison provides a measure of confidence that using ERA-20C winds to drive PIOMAS-20C generates a realistic representation of sea ice drift distribution and therefore capture sea ice variability associated with advective processes, even during a time when data coverage was sparse. However, it is difficult to assess whether this example is representative for other locations within the Arctic during the early part of the twentieth century before more frequent observations became available for assimilation. Examination of the feedback records for ERA-20C indicates that SLP measurements from the Maud were assimilated in ERA-20C, thereby likely positively affecting the resulting wind field.
7. PIOMAS 20C sea ice volume
a. Total ice volume: Comparison of PIOMAS-20C with PIOMAS
Figure 12 shows the total ice volume time series from PIOMAS-20C for April and September from 1901 through 2010. Standard PIOMAS products from 1979 to 2016 are shown for comparison. April sea ice volume for PIOMAS-20C is larger by 4.3 × 103 km3, a 13% difference in total volume for that time of the year. PIOMAS-20C September volume for the 1979–2010 period is 2.6 × 103 km3 or 17% larger than PIOMAS. These volume differences are larger than the previously established uncertainty estimates of 2.8 × 103 and 1.2 × 103 km3 for April and September, respectively (Schweiger et al. 2011).
The differences between the PIOMAS and PIOMAS-20C ice volume results can be attributed to three sources that influence the resulting sea ice volume. Those are the atmospheric forcing data, the input data used for the assimilation of sea ice information, and the method by which sea ice information is assimilated. To investigate the role of each of these mechanisms on total volume we conduct two separate experiments (EXP1 and EXP2). The potential influence of model tuning is included in the sensitivity to forcing data and is considered minor since the same tuning procedure is followed. The details for each experiment are provided in Table 1.
Table 1 shows a comparison of volume differences and trends relative to the PIOMAS-20C integration. The annual volume difference between PIOMAS and PIOMAS-20C is 15% of the total volume. EXP1 uses forcing from NCEP-R1 but leaves everything else as in the PIOMAS-20C integration. This experiment generates sea ice volume that is 6% lower on annual average than PIOMAS-20C, indicating that 40% of the total difference between PIOMAS-20C and PIOMAS is due to different forcing data. EXP2 uses the same forcing and input data as PIOMAS but assimilates sea ice observations near the ice edge. It has similar results as EXP1, with annual ice volume lower than PIOMAS-20C by 6%. This suggests that for annual sea ice volume, forcing data account for 40% of the difference between PIOMAS-20C and the assimilation method accounts for 60% of the difference. The effect of the assimilated sea ice dataset is relatively small, even though differences exist (Titchner and Rayner 2014). However, it is important to note that these factors are not entirely independent since the effect of the assimilation method will depend on the ice extent and the ice concentration gradients at the ice edge in the dataset being used. The overall effect of the assimilation procedure is to remove sea ice relative to an integration without assimilation. Sea ice is removed in the North Atlantic where the ice-edge position tends to be too far south without assimilation. The removal of sea ice at the edge also has an impact on the sea ice thickness and total volume, because the effect of corrections that occurred at the ice edge persists and propagates through the system. However, it does not affect the ice volume trends during the ETCW while it increases the simulated ice losses during the recent period (see the supplemental material for a further discussion).
b. Total ice volume trends
Atmospheric forcing, sea ice information data source, and assimilation method also affect total volume trends. The PIOMAS-20C volume loss from 1979 to 2010 is −0.37 × 103 km3 yr−1 while PIOMAS volume loss for the same periods is only −0.28 × 103 km3 yr−1, a relative change of 25%. For annual volume trends from 1979 to 2010, the atmospheric forcing has the largest impact, accounting for a relative differences of 16%. The assimilation method accounts for 9% of the volume trend difference and the input dataset for 7%. The differences in trends are similar in size as previously found using different forcing datasets or assimilation methods (Schweiger et al. 2011; Lindsay et al. 2014). They are also consistent with our previous assessment that PIOMAS trends have relative uncertainties of about 30% (Schweiger et al. 2011).
Using the trends from the respective runs to estimate relative sea ice volume loss, we have lost 65% of the total Arctic sea ice volume in September from 1979 to 2010 using PIOMAS-20C. This number is reduced to 55% for PIOMAS. Ice loss for the ice volume maximum in April is smaller, with 41% for PIOMAS-20C and 35% for PIOMAS. The greater volume loss trends in PIOMAS-20C relative to PIOMAS arise from the fact that ice thickness in the Beaufort and Chukchi Seas is larger in PIOMAS-20C during the early 1980s than for PIOMAS. Extending the period through 2016 shows a total volume loss of 72% for September and 35% for April from PIOMAS.
Trend differences arise mostly from sea ice volume differences during the early 1980s when the two reconstructions show the largest differences. Examining ice thickness differences from ThickCDR provides additional clues. ThickCDR data from 1979 to 1984 from U.S. submarines shows a small bias (−0.03 m) for PIOMAS-20C and a larger one for PIOMAS (−0.61 m). Therefore, PIOMAS-20C ice thickness may indeed reflect ice thickness during the early 1980s more accurately than PIOMAS. Ice-ocean model integrations with atmospheric forcing fields from different reanalysis datasets (Lindsay et al. 2014; see their Fig. 13) but without data assimilation, also showed that the integration using NCEP-R1 forcing had substantially lower volume during the early 1980s. While this is supportive of our previous result that PIOMAS provides a conservative estimate of sea ice volume loss (Schweiger et al. 2011), we do not believe it provides sufficient evidence to decide which of the simulations provides the more accurate estimate of sea ice volume and trend. Instead the differences should be viewed as a measure of the uncertainty in the reconstructed sea ice volume in either dataset.
c. Early-twentieth-century warming
In the early twentieth century the Arctic experienced a well-documented warming period (Beitsch et al. 2014; Bengtsson et al. 2004; Brooks 1923; Kincer 1933; Suo et al. 2013; Wood and Overland 2010; Zubov 1948). The start and end points of the ETCW period in the Arctic vary somewhat between authors, although starting dates around 1920 are commonly given, and ending dates in the 1940s. Recognizing that the start and end points for this period are not clearly defined, we consider a longer period from the start of the record in 1901 through 1940. Although trends over the period are sensitive to the selection of the start and end points of the period, the overall conclusions we draw here are not. While the warming is well established, less well documented is the impact of this warming on sea ice, particularly on sea ice thickness and total volume. PIOMAS-20C sea ice volume anomalies (Fig. 13) show a downward trend from 1901 through about 1940. Sea ice volume over this period decreased by 600 km3 decade−1. This compares with a decrease of 3810 km3 decade−1 over the 1979–2010 period. Ice volume increased from 1940 through the mid-1950s, providing some justification for selecting 1940 as the end point of the ETCW.
Information from ice charts, vessel logs, and reports published at the time (e.g., CIEM/ICES 1948; Koch 1945; Zubov 1948) suggests that sea ice loss during the ETCW may not have been Arctic wide but featured a strong imprint on sea ice in the Atlantic sector. Little evidence is shown for a similar imprint in the Pacific sector (Wood and Overland 2010). PIOMAS-20C simulations provide further support for these results. Figure 14a shows the trend of sea ice thickness from 1901 through 1940 for September. Ice thickness in the Atlantic sector decreased considerably, by as much about 0.3 m decade−1 in an area reaching from Fram Strait to the Barents Sea. A simultaneous thickening occurred in the East Siberian, Chukchi, and Beaufort Seas. The decrease in the Atlantic sector qualitatively matches the decline in ice extent over the same period reported by others (Divine and Dick 2006) who provide ice extent from 1750 to 2002 for the sector covering 30°–70°W. The sea ice thickness trend pattern largely matches the surface air temperature pattern in the ERA-20C data (Fig. 15), with strong warming reaching from the Atlantic side of the Arctic deep into the Arctic Ocean. A simultaneous cooling in the Chukchi and Beaufort Seas is consistent with sea ice thickening in this area.
Further support for this temperature and ice thickness trend pattern comes from the NOAA 20RCv2c reanalysis, which also shows warming in the Atlantic sector and cooling on the Pacific side (Fig. S7). Although details and magnitude of the temperature anomaly are different, the general pattern with warming in the Atlantic sector and cooling in the Pacific is similar. Regional sea ice thickness fluctuation associated with the ETCW in the Arctic can be contrasted with the more recent decline in sea ice thickness from 1979 to 2010 (Fig. 14b), which shows a decrease in sea ice thickness throughout the Arctic, with a maximum thickness decrease along the Pacific side of the Arctic (see Fig. S9 for comparison the surface air temperature trends for 1979–2010). Note that the standard PIOMAS thickness trend pattern over the 1979–2010 period is very similar to the PIOMAS-20C reconstruction, but it shows a weaker decline in sea ice thickness on the Siberian side. This is because PIOMAS-20C has thicker ice in this area in the early 1980s, as discussed above. The spatial pattern of the sea ice thickness differences between the model and ICESat observations discussed earlier has the potential to interact with the spatial pattern in the thermal forcing due to the ice growth feedback (Bitz and Roe 2004), which dictates that thinner ice grows more rapidly than thick ice. However, since this ice growth feedback is not linear it is possible that the negative model bias in the area of thinning during the ETCW may have reduced the rate of modeled ice volume loss.
d. Examining sea ice variability in the Pacific sector: Logs from U.S. Revenue Cutters
While the sea ice variability in the Atlantic sector of the Arctic during the ETCW is well characterized on the basis of historical records, relatively few observations have been available for the Pacific sector of the Arctic. Many records that exist have so far not been incorporated into the Hadley ISST2 V220.127.116.11 data, which in turn provide the ice concentrations assimilated into PIOMAS-20C.
Using a new dataset compiled as part of this project, we examine Pacific sector sea ice variability from PIOMAS-20C during the first four decades of the twentieth century, which includes the ETCW period. Historical ship logbooks from U.S. government ships, including the U.S. Navy, U.S. Revenue Cutter Service/U.S. Coast Guard, U.S. Coast Survey, and U.S. Fisheries Service located at the National Archives were digitally imaged as part of a joint NOAA/National Archives effort, and then transcribed by citizen scientists participating in the Old Weather (www.oldweather.org) project. The period captured spans 1844–1970, with near-comprehensive coverage of logbooks relating to the Arctic to 1955.
From these, Pacific sector logbooks where sea ice was reported were selected for enhanced analysis, including reconstruction of the ship tracks from navigational data recorded in the logs (e.g., course/distance run and bearing/range information), and then coded for sea ice conditions at hourly resolution. Sea ice information is coded as ice present or absent, along with descriptive terms for both sea ice and vessel operating condition that together provide a qualitative assessment of sea ice conditions. For the purpose of this project we only use the ice/no ice information in this dataset. Forty ship logs from 1901 through 1938 are included here. Note that the U.S. Coast Guard cutter Bear logbook for 1920 is missing, and there were no Arctic cruises in 1927 and 1939. The U.S. Coast Guard shifted assets to the Atlantic sector in 1940 in response to the outbreak of World War 2 in Europe; regular operations in the Pacific sector of the Arctic resumed in 1946.
Figure 16 shows a comparison of observed ice versus ice-free conditions in the Bering/Chukchi Sea from these logs, overlaid on sea ice concentrations from the PIOMAS-20C reconstruction for June and August for four periods, 1901–10, 1911–20, 1921–30, and 1930–38.
We draw the following conclusions from this comparison: ship observations of ice are generally located north of the ice edge, while those reporting ice-free conditions are located south of the ice edge. A quantitative comparison is shown in Table 2. Of 8091 ship observations, 74% were correctly represented in PIOMAS-20C and 26% were misclassified. The vast majority of misclassified observations are those where ship observations indicate open water while PIOMAS-20C reports sea ice. To an extent this is an expected bias, given that observations of ice presence are definitive while in many cases no ice reported is not. For example, persistent low visibility may have hampered the sighting of sea ice that was actually present at some distance from the ship. This type of open-water mischaracterization bias near the ice edge is also fostered by the navigating officer’s explicit selection of leads and open-water routes wherever possible. No apparent trend in the sea ice edge or the quality of the comparison statistics is apparent over the four decades, suggesting no discernable multidecadal trend. This provides further evidence that the sea ice decline during the ETCW period was largely confined to the Atlantic sector, while the Pacific sector showed no strong trend (Wood and Overland 2010). Since the ice edge in PIOMAS-20C is assimilated using the HadISST2 dataset, this comparison indirectly provides some validation of the HadISST2 dataset, which during that period entirely relies on Walsh and Chapman (2001), who utilized sea ice charts produced by the Danish Meteorological Institute (DMI; https://nsidc.org/data/G02203) for this time period.
e. Uncertainty in ice volume trends for the ETCW
Given the lack of sufficient validation data for the ETCW period it is difficult to directly quantify the uncertainty of the sea ice volume trend. Since ERA-20C only provides a single-member simulation, an uncertainty estimate based on ensemble simulations using this dataset is not possible, and uncertainties due to biases in the atmospheric reanalysis could not be addressed in this fashion in any case. To provide some measure of the uncertainty in the ETCW volume trends we use a different approach. Previous research has utilized the relationship between sea ice extent and surface air temperatures as a way to reconstruct sea ice conditions over centuries (Alekseev et al. 2016; Connolly et al. 2017). We also utilize this approach to establish the relationship between 1901 and 2010 sea ice volume anomalies (all months and September) from PIOMAS-20C and surface air temperatures temperature anomalies (all months and September) from ERA-20C over the PIOMAS-20C model domain (ocean areas north of ~49°N). These linear relationships have correlations of −0.49 for all months and −0.64 for September (coefficients given in Table 3). The regression equation is then used to reconstruct sea ice volume based on air temperature anomalies from ERA-20C, CERA-20C, and 20CRv2c. For the latter two, ensemble means are used. Figure 17 shows PIOMAS-20C sea ice volume, temperature anomalies for all months in the year over the PIOMAS-20C model domain from ERA-20C, CERA-20C, and 20CRv2c. Although air temperature anomalies are clearly correlated with sea ice volume anomalies, there is also considerable variability that is not related to average temperature anomalies. This variability is generated by dynamic processes including ocean and sea ice dynamics, by responses to thermodynamic forces not captured by the domain averages, or by interactions between them. In fact, thermodynamic and dynamic (wind) forcing contribute in equal parts to Arctic sea ice volume variability based on dedicated experiments (Koberle and Gerdes 2003; Rothrock and Zhang 2005). While temperature-based sea ice reconstruction therefore can only account for part of the sea ice variability, they can provide some measure of the uncertainty for our derived ETCW sea ice volume trend due to differences in trends in the forcing data. Using the range of reconstructed sea ice volume anomalies as a measure of uncertainty, we can characterize the PIOMAS-20C ice volume anomaly time series (Fig. 18). An examination of temperature and reconstructed sea ice volume anomalies during the ETCW (Fig. 17a) shows that ERA-20C has the weakest temperature increase and smallest sea ice volume decrease during the ETCW. Table 4 provides PIOMAS-20C and trends reconstructed from air temperature anomalies for 1901–40 and 1980–2010. PIOMAS-20C sea ice volume trends from 1980 to 2010 period are about 6 times larger than during the 1901–40 ETCW period. ERA-20C temperature-based volume trends for the ETCW period are about half (−0.27 × 103 km3 decade−1) than when sea ice dynamics (−0.56 × 103 km3 decade−1) are included. Using CERA-20C and 20CRv2c temperature anomalies to reconstruct sea ice volume during the ETCW yields significantly larger sea ice volume losses (−0.43 × 103 km3 decade−1 and −0.45 × 103 km3 decade−1, respectively) because of the relatively stronger warming during that period in both of these reanalysis products. Applying the same approach to reconstructing sea ice volume for September based on air temperature anomalies yields similar results (Fig. 18b). While these temperature-based reconstructions indicate that PIOMAS-20C ice volume losses during the ETCW are likely conservative estimates, the strong contrast with the more recent warming remains robust. Conservatively estimating the uncertainty of the ETCW sea ice loss based on temperature sensitivities as 100% would yield a maximal 1200 km3 decade−1 sea ice volume loss during the ETCW. Comparing this number to the more reliably known sea ice losses during the more recent period, shows that the 1979–2010 losses are still larger by a factor of 3.
Our conclusion that sea ice loss during the ETCW period was drastically smaller than during the 1980–2010 period is further supported by comparing PIOMAS-20C ice concentration anomalies with a new ice concentration dataset compiled by Walsh et al. (2016) (hereafter SIBT-1850) (Fig. 19). PIOMAS-20C ice concentration anomalies closely correspond to the SIBT-1850 during the ETCW and show perhaps a weak downward trend in both datasets from 1920 through 1940. This comparison indicates that ice concentration anomalies for the relevant periods compared are relatively unaffected by the differences between the more data-rich SIBT-1850 dataset and the HadISST v18.104.22.168 that was assimilated into PIOMAS-20C. The largest differences occur during the World War II period (1939–45) when SIBT-1850 include additional datasets that had not been incorporated into the version of the HadISST v22.214.171.124 dataset, which provides the ice concentration information for PIOMAS-20C. This discrepancy in ice concentration during this period needs to be taken into consideration when using PIOMAS-20C.
8. Summary and conclusions
We have reconstructed a 110-yr-long record of Arctic sea ice thickness and volume using the PIOMAS ice-ocean assimilation system forced with ERA-20C atmospheric data. ERA-20C, which only assimilates surface level pressure and marine winds, performs similar to other, more comprehensive reanalyses with respect to the annual cycles of downwelling radiation, surface air temperature and wind speed. For surface air temperature, RMS errors are similar to RMS errors from NCEP-R1 but nearly double to those for ERA-Interim. Wind speed validation statistics are not as good for ERA-20C than for other reanalysis projects, including 20CR. Validation statistics for temperature, wind, and sea ice drift for the early data sparse presatellite period, show the considerable skill of ERA-20C in capturing daily to monthly variability. Sea ice drift from PIOMAS-20C shows that significant fractions of the variance of the historical observations from the Maud expedition (1922–24) and the ARLIS-2 drifting station (1960–65) are captured.
Comparisons of PIOMAS-20C ice thickness with historical sea ice thickness observations going back to the Maud expedition (1922–24) provide similar error statistics as for the standard PIOMAS reconstruction over 1979–2010. Mean annual ice thickness differences between observations and reconstruction have a small bias of 0.15 m, an RMS error of 0.26 m, and they are correlated with an r value of 0.85. Total sea ice volume from PIOMAS-20C is generally larger than PIOMAS over 1979–2010.
Substantial differences between PIOMAS and PIOMAS-20C in total sea ice volume occur in spring and can be attributed to roughly equal parts to atmospheric forcing and differences in the assimilation method.
Ice thickness patterns for the ETCW periods show a decline in sea ice thickness in the Atlantic sector of the Arctic but increases in the Pacific sector. This pattern of thickness variability is consistent with the pattern of surface temperature anomalies that occurred over this time period. Comparison of PIOMAS-20C with sea ice information from newly transcribed logs from U.S. Revenue Service/Coast Guard cutters operating in the Bering and Chukchi Seas between 1901 and 1938 support the notion that Pacific sea ice did not show a strong trend during the ETCW. Total Arctic sea ice volume loss during the ETCW period from 1901 to 1940 is only 600 km3 decade−1, 6 times smaller when compared with the loss of 3800 km3 decade−1 for the more recent 1979–2010 period. Using temperature-based reconstructions of sea ice volume using 20CRv2c and CERA-20C we show that considerable uncertainty remains for sea ice volume trends during the ETCW. PIOMAS-20C sea ice losses during the ETCW are likely conservative estimates because of the subdued warming in ERA-20C during that period relative to other reanalyses. Nevertheless, a much stronger decline in sea ice volume during the 1980–2010 period compared the ETCW period is robust feature of these reconstructions.
PIOMAS-20C provides a first step toward century-scale sea ice reanalysis. An important limitation of the approach using a coupled ice–ocean modeling and assimilation system in combination with an atmospheric reanalysis is the lack of coupling between the atmosphere and the ice–ocean system. The sea ice dataset used to provide boundary conditions in the atmospheric reanalysis will damp to some extent the ice–ocean model solution forced with the atmospheric reanalysis. The resulting PIOMAS-20C ice thickness is therefore not entirely independent of the sea ice information utilized in the atmospheric reanalysis. Coupled air–ice–ocean reanalyses may remove this limitation in the future but other difficulties, such as the uncertainties in the coupling of various model components, will have to be overcome first to assure realistic sea ice thickness simulations.
The data assimilation approach used in this study is a very simple one. It blends model-estimated ice concentration with satellite observations at the ice edge. Because of the complete coverage of the available satellite ice concentration data in space and time, this simple approach allows the model to adjust to the observations at each time step and at each ice-edge grid cell in an efficient manner. As a result, the simulated ice edge after assimilation resembles the observed one, while ice concentration and thickness in the interior of the ice pack are not directly constrained by observations. However, like many other data assimilation techniques, this approach does not conserve mass or energy. In addition, this approach does not directly address the source of the biases that may exist in the system including forcing, model physics, and uncertainties in the assimilated data. Because of these potential biases/uncertainties, we rely on extensive model calibration and validation using a variety of data over a significant period of time as mentioned in section 5. Other data assimilation methods that attempt to conserve mass and energy (Fenty et al. 2017; Kauker et al. 2009; Koldunov et al. 2017) may have the potential to create better ice thickness and volume products in the future. Data assimilation approaches that yield direct estimates of uncertainty of the integrated parameter (e.g., ice volume) by generating ensembles might be explored.
Another path forward is the direct assimilation of historical sea ice information such as qualitative sea ice observations from shipping logs. Sea ice information from in situ datasets is typically sparse in both time and space. In the present approach, gridded sea ice datasets that are space and time complete are generated from ship logs, charts, and satellite data, using various techniques to assemble the data into a gridded product. This typically involves cross calibration of datasets, derivation of ice concentration from extent based on historical gradients across the ice edge, and ample application of extrapolation and default to climatology (Mahoney et al. 2008) or the use of gap-filling (Walsh et al. 2016). Observed ice thickness information is not currently used at all in the development of this sea ice reconstruction other than for calibration and validation. This may be improved in the future by assimilating available thickness data. As part of this study we have assembled a host of historical sea ice information from ship logs, historical charts, and field programs. We expect future version of PIOMAS-20C to utilize this information through direct assimilation of these data types as well.
We thank Gil Compo for early access to 20CRv2c data and for useful discussions regarding this data product. We like to thank the four anonymous reviewers for their constructive and detailed reviews. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC02-05CH11231. Financial support for this work came from NSF Grants ARC-1203425, PLR-1416920, PLR-1603259, and OPP-1744587 and NASA Grants NNX15AG68G and NNX17AD27G, and the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under NOAA Cooperative Agreement NA15OAR4320063. The joint NOAA/National Archives digital imaging program, which provided online access to primary source logbooks, was supported by North Pacific Research Board Grants 1124, 1226, and 1322. This work uses data provided by Old Weather (https://www.oldweather.org), a Zooniverse project. Data access for several reanalysis data sets was through the Data Support Section of the Computational and Information Systems Laboratory at the National Center for Atmospheric Research. NCAR is supported by grants from the National Science Foundation. Sea ice thickness observations are available from http://psc.apl.uw.edu/sea_ice_cdr/. PIOMAS-20C sea ice thickness fields are available at http://psc.apl.uw.edu/research/projects/piomas-20c/. CERA-20C data were obtained from http://apps.ecmwf.int/datasets/data/cera20c-edmm/levtype=sfc/type=an/.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-19-0008.s1.