Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
39,349 characters · 14 sections · 43 citation commands
Optimal Combination of Arctic Sea Ice Extent Measures: A Dynamic Factor Modeling Approach
\setcounter{page}{1} \thispagestyle{empty}
Climate change is among the most pressing issues of our time, with many severe economic, environmental, and geopolitical consequences. Recently, the application of time series analytical methods to this topic -- and, more broadly, a “climate econometrics" -- has emerged as a vibrant research literature, as highlighted, for example, in Hillebrandetal2020 and the references therein. One important issue that these methods can address is the loss of Arctic sea ice. The loss of Arctic sea ice is a vital focus point of climate study. It is both an ongoing conspicuous effect of climate change and a cause of additional climate change via feedback loops. In particular, reduced Arctic sea ice boosts solar energy absorption via decreased albedo due to darkening color (e.g., Stroeve2012, PistoneEtAl2019, DRice) and increased methane release due to melting permafrost (e.g., VaksEtAl2020).\footnote{For a broad and insightful overview of the evolution and causes of reduced Arctic sea ice cover, see SeaIceBookCh4.}
There are, however, several alternative measures of Arctic sea ice extent based on different processing methodologies of the underlying satellite-based microwave measurement data, and the choice among these measures is not clear-cut (BunzelEtAl16). In this paper, we study four such sea ice extent ($SIE$) measures, which we denote as Sea Ice Index ($SIE^S$), Goddard Bootstrap ($SIE^G$), JAXA ($SIE^J$), and Bremen ($SIE^B$). The top panel of Figure (ref) provides time series plots of these four measures of Arctic $SIE$ for the satellite measurement era, which started in 1978. The four measures appear almost identical, because their scale is dominated by large seasonal swings. However, the effects of seasonality can be removed by plotting each month separately for the four series, as done in the lower twelve panels of Figure (ref). Of course, the Arctic $SIE$ measures all trend down in every month, with steeper trends for the low-ice “summer" months (e.g., August, September, October). (Note the different axis scales for different months.) There are also systematic differences across indicators. $SIE^G$, for example, tends to be high, and $SIE^J$ tends to be low, while $SIE^S$ and $SIE^B$ are intermediate. But the deviations between various pairs of measures are not rigid; that is, they are not simply parallel translations of each other. Instead, there are sizable time-varying differences among the various measures.
All of this suggests treating the various measures as noisy indicators of latent true sea ice extent, which in turn suggests the possibility of blending them into a single combined indicator with less measurement error. Indeed some prominent studies have used simple equally-weighted averages of competing indicators, with precisely that goal. For example, a recent report on the state of the cryosphere (IPCC19) uses a simple average of three indicators.\footnote{See the notes for their Figure 3.3, page 3-13.} Simple averages, however, are often sub-optimal. Optimality generally requires use of weighted averages giving, for example, less weight to noisier indicators. Motivated by these considerations, in this paper, we propose and explore a dynamic factor state-space model that combines the various published indicators into an optimal measure of sea ice extent, which we extract using the Kalman smoother.
We proceed as follows. In section (ref), we describe the four leading Arctic sea ice extent indicators that we study, and the satellites, sensors, and algorithms used to produce them. In section (ref), we propose a basic dynamic-factor state-space model for sea ice extent and use it to obtain optimal extractions of latent extent. We conclude in section (ref).
Sea ice extent ($SIE$) indicators are constructed from satellite measurements of the earth's surface using passive microwave sensing, which is unaffected by cloud cover or a lack of sunlight. Several steps are necessary to convert raw reflectivity observations into final $SIE$ measurements. First, for a polar region divided into a grid of individual cells, various sensors record a brightness reading or “brightness temperature" for each cell. An algorithm then transforms these brightness readings into fractional surface coverage estimates -- sea ice concentration ($SIC$) values -- for each grid cell. Finally, $SIE$ is calculated by summing the area of all cells with at least 15 percent ice surface coverage.\footnote{PARKINSONETAL2008 discuss reasons for using a 15 percent cutoff.} This up-rounding in $SIE$ is effectively a bias correction, as determining the edge between ice and water can be especially difficult in the summer, when, for example, melting pools on summer ice surfaces can be mistaken for ice-free open water (MeierAndStewart2019).
Different algorithms for processing the raw measurements importantly shape the final $SIE$ estimates. In addition, the $SIE$ series are not based on identical raw data because they use somewhat different satellites and sensors (ComisoEtAl17, NSIDC_Boot_Explained). In this section, we review some aspects of the satellites, sensors, and algorithms that underlie the $SIE$ measures.
Table (ref) summarizes the operative dates of the various satellites and sensors relevant for Arctic sea ice measurement. The first multi-frequency sensor equipped on a satellite was the {Scanning Multichannel Microwave Radiometer} (SMMR) launched in 1978 (NSIDC_Caval). Starting in 1987, later sensors -- the Special Sensor Microwave Imager (SSM/I) and the {Special Sensor Microwave Imager/Sounder} (SSMIS) -- offered higher resolution images.\footnote{For detailed discussion of sensor characteristics see \url{https://nsidc.org/ancillary-pages/smmr-ssmi-ssmis-sensors}.} In 2002 and 2012, respectively, the {Advanced Microwave Scanner Radiometer for EOS} (AMSR-E) and {Advanced Microwave Scanner Radiometer 2} (AMSR2) sensors were launched and provided further improvements in resolution (ComisoEtAl17).\footnote{Early in the sample, operational problems prevented data delivery for several days during 1986 and between December 1987 and January 1988 (NSIDC_Boot). For more recent technical difficulties, see \url{https://www.nrl.navy.mil/WindSat/Description.php}.} Given the inclinations of the satellite orbits and the spherical shape of the earth, all of the satellites share an inability to observe the Arctic “pole hole" -- a circular region at the very top of the world. The size of the pole hole varies across sensors, but historically, there is full confidence that the area covered by the pole hole fulfills the 15 percent $SIC$ requirement (MeierAndStewart2019).
Table (ref) also describes the underlying source data for our four $SIE$ indicators. These measures use algorithms to transform the raw satellite brightness data into $SIC$ and $SIE$ values. We now turn to a more detailed discussion of these algorithms to illuminate the differences across $SIE$ indicators.
Once brightness data have been recorded by satellite sensors, an algorithm converts the measurements into estimates of $SIC$. Here, we discuss the algorithms and other details of the various $SIE$ indicators.
Updated on a daily basis and distributed by the {National Snow and Ice Data Center} (NSIDC), the Sea Ice Index (SII or $SIE^S$) combines two separate Sea Ice indicators: (1) the {Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data} (NASA Team) (NSIDC_Caval) -- produced at the {Goddard Space Flight Center} -- and (2) the {Near-Real-Time DMSP SSMIS Daily Polar Gridded Sea Ice Concentrations} (NRTSI) (NSIDC_NearReal) -- produced by the NSIDC itself.\footnote{See \url{https://doi.org/10.7265/N5K072F8}.} A time-lag of about one year between the $SIE$ estimates by NASA Team and its publication in the NSIDC database requires the NRTSI to complement the SII.
The NRTSI follows the NASA Team algorithm as closely as possible, but inconsistencies between the two series cannot be ruled out entirely (NSIDC_SII). In particular, the two sub-indicators use brightness temperatures from different providers.\footnote{NSIDC_NearReal takes the data from the National Oceanic and Atmospheric Administration Comprehensive Large Array-data Stewardship System (NOAA CLASS); (NSIDC_Caval) uses data processed at the NASA Goddard Space Flight Center.} These raw readings can be distorted by weather effects, making open water look like sea ice cover. Therefore, post-calculation quality checks apply land and ocean masks, to remove errorneous and implausible ice covers. However, NASA Team and NRTSI do not apply the exact same filters (NSIDC_SII). The former algorithm additionally screens the data manually for falsely detected ice formation (NSIDC_Caval), which can enhance accuracy but also reduce transparency of the final measurements.
As Table (ref) shows, the SII obtains the raw data from different generations of satellites and sensors. To make the data comparable, a linear least-squares model on the brightness temperatures, as reported by the two distinct sensors for an overlapping period of operation, is intended to adjust the reference points of 100 percent sea ice and 100 percent open water. These tie points then remain fixed over the lifetime of the new system (NSIDC_Caval_tie).
Another sea ice indicator, distributed by the NSIDC, relies on $SIC$ estimates from the Goddard Bootstrap algorithm ($SIE^G$).\footnote{For detailed algorithm description see ComisoEtAl17b. For the data, see \url{https://doi.org/10.5067/7Q8HCCWS4I0R}.} Despite the NASA Team and the Goddard Bootstrap algorithms having both been developed at the NASA Goddard Space Flight Center, there are some differences between the two approaches. These arise mostly from the calibration of tie points: While the NASA Team adjusts these reference points for 100 percent open water and 100 percent ice only when a new satellite or sensor becomes operational, the Goddard Bootstrap algorithm adjusts these reference points on a daily basis to account for varying weather conditions (ComisoEtAl17). Differing weather filters and sensitivities to varying physical temperature also lead to differences in the final measurements(ComisoEtAl97). In contrast to the NASA Team, the strength of the Goddard Bootstrap algorithm is the identification of melting sea ice. Therefore, the Goddard Bootstrap algorithm provides more accurate estimates of the edge of the ice cover (goldsteinetal2018).
Although differences between $SIE^S$ and $SIE^G$ are generally assessed to be small, they cannot necessarily be neglected (goldsteinetal2018).\footnote{See also \url{https://nsidc.org/support/faq/nasa-team-vs-bootstrap-algorithm}.} The differences between the NASA Team and Goddard Bootstrap algorithms occur especially during the melting season, when the former generally reports larger deviations from ship or radar observations. However, the relative accuracy of the two algorithms is not clear-cut, as the Goddard Bootstrap algorithm is highly sensitive to physical temperature and underestimates $SIC$ during winter times in the higher latitudes of the Arctic region.
As listed in Table (ref), both $SIE^S$ and $SIE^G$ rely on the same set of instruments, which have been criticized for their low spatial resolution (goldsteinetal2018). The Japan Aerospace Exploration Agency (JAXA) sea ice measure ($SIE^J$) uses an adapted version of the Goddard Bootstrap algorithm to derive $SIC$ measures from satellite readings with higher spatial resolution (ComisoEtAl17b).\footnote{For description and data see \url{https://kuroshio.eorc.jaxa.jp/JASMES/climate/index.html}.} However, readings from these high resolution satellites are only available since 2000, so their data must merged with observations from older sensors to extend the data coverage to 1978 (ComisoEtAl17). $SIE^J$ also distinguishes itself from $SIE^S$ and $SIE^G$ by using a 5-day moving average of observations to compensate for potentially missing data.
Using observations delivered by the high-resolution AMSR-E sensor, a group of researchers at the University of Bremen developed the ARTIST Sea Ice (ASI) algorithm (SpreenEtAl08) to estimate daily $SIC$.\footnote{Monthly data from \url{https://seaice.uni-bremen.de/data/amsr2/today/extent_n_19720101-20181231_amsr2.txt}.} The time-series ($SIE^B$) uses different algorithms for different sensors. Until the launch of the AMSR-E sensor in 2003, $SIE^B$ used the NASA Team algorithm to transform brightness readings into $SIC$ values. From then on, the NASA Team algorithm is replaced by the ASI.\footnote{See \url{https://seaice.uni-bremen.de/sea-ice-concentration-amsr-eamsr2/time-series/}.}
The four sea ice indicators discussed above differ in terms of the raw data sources and the algorithms used to process the raw data. They can be viewed as distinct indicators of an unobserved or latent “true” sea ice extent, $SIE^*$. Blending such a set of noisy indicators can produce a single series with less measurement error. Here we formalize this intuition in a state-space dynamic-factor model, from which we extract an optimal composite estimate of $SIE^*$ from the four component indicators.\footnote{Dynamic factor analysis is closely related to principal components analysis, but the dynamic factor model provides a fully-specified probabilistic modeling framework, in which estimation and factor extraction via the Kalman filter are statistically efficient. Under conditions the two approaches coincide in large samples, but those conditions include a very large number of indicators. Those conditions are violated in our case as we have only four indicators, so dynamic factor analysis is preferable. For a much more complete account, see SW2011.}
We work in a state-space environment, modeling each of the four indicators ($SIE^S$, $SIE^J$, $SIE^B$, and $SIE^G$) as driven by latent true sea ice extent, $SIE^*$, with an additive measurement error.\footnote{The approach parallels ADNSS2016, who extract latent “true" U.S. GDP from noisy expenditure-side and income-side estimates.} As discussed previously, some indicators present level shifts with respect to one another, most often resulting from how they respectively deal with tie points. It is thus preferable to constrain both the trend and dynamics to follow a common factor, while leaving level offsets unconstrained. The measurement equation is
where
with
Moreover, $c_i = 0$ when we normalize $\lambda_i = 1$ for $i \in \{S, J, B, G\}$.
The transition equation is
where $\eta_t \sim iid (0, \sigma^2_{\eta \eta})$ is orthogonal to $\varepsilon_t$ at all leads and lags. Various modeling approaches are distinguished by their treatment of $TREND_t$ and $SEASONAL_t$. We follow DRice and allow for 12 monthly deterministic seasonal effects, each of which is endowed with (possible) deterministic quadratic trend.\footnote{We emphasize that our model is meant to be a simple benchmark, and that many potentially-important variations and extensions are possible. For example, one could alternatively entertain stochastic as opposed to deterministic trend and seasonality. A simple approach would be separate month-by month modeling so that there is no seasonality, whether with one unit root, as in (for month $m$) $ TREND_{m,t} = d_m + TREND_{m,t-1} + u_{m,t}$, or two unit roots, as in $TREND_{m,t} = d_{mt} + TREND_{m,t-1} + u_{m,t}$, where $d_{m,t} = d_{m,t-1} + v_{m,t}$. } This results in a blended deterministic “trend/seasonal" given by
where $D_i$ indicates month $i$ and $TIME$ indicates time. Hence the full transition equation is
The model is already in state-space form, and one pass of the Kalman filter, initialized with the unconditional state mean and covariance matrix, provides the 1-step prediction errors necessary to construct the Gaussian pseudo-likelihood, which we maximize using the EM algorithm and calculate standard errors from the analytic Hessian matrix.\footnote{Note that we do not assume Gaussian shocks, and that it is not necessary to assume Gaussian shocks, as we can still maximize the Gaussian likelihood even if the shocks are not truly Gaussian, and the resulting (pseudo-)MLE still has the good properties of consistency, asymptotic normality, etc.} Following estimation, we use the Kalman smoother to obtain the best linear unbiased extraction of $SIE^*$ from the estimated model. The smoother averages across indicators, but it desirably produces optimally weighted averages rather than simple averages. The smoother also averages over time, using data both before and after time $t$ to estimate $SIE_t^*$, which is also necessary for optimal extraction, due to the serial correlation in $SIE^*$. For details see Harvey1989.
One or more restrictions are necessary for identification. The standard approach is to normalize a factor loading, which amounts to an unbiasedness assumption. Normalizing $\lambda_S {=} 1$, for example, amounts to an assumption that $SIE^S$ is unbiased for $SIE^*$. Whether there truly exists such an unbiased indicator (and if so, which) is of course an open question -- one can never know for sure. $SIE^S$ and $SIE^G$ are the most widely used indicators (MeierEtAl14, PengEtAl13), so it is natural to consider normalizing on $\lambda_S$ or $\lambda_G$. We explore both.
The estimated measurement equation ((ref)), normalized with $\lambda_S {=} 1$, is
where standard errors appear beneath each estimated loading. All indicators are estimated to load heavily on $SIE^*$, with all $\hat{\lambda}'s$ very close to 1. $SIE^J$ and $SIE^G$ load least heavily ($\hat{\lambda_J} {=} 0.950$, $\hat{\lambda_G} {=} 0.961$)), in accord with their generally less abrupt trend in Figure (ref). $SIE^B$ loads with an estimated coefficient marginally different from 1 but significant at the 5% level. Of course, the $SIE^S$ loading is 1 by construction. While Bremen's level offset (with respect to $SII$) is arguably negligible, those of Jaxa and especially Goddard are sizable. Hence, the estimation results accord with Figure (ref), with $SIE^S$ and $SIE^B$ more in the center of the range and $SIE^J$ and $SIE^G$ more extreme.
Alternatively, the estimated measurement equation normalized with $\lambda_G {=} 1$ is
The estimated loadings in equations ((ref)) and ((ref)), corresponding to $\lambda_S {=} 1$ and $\lambda_G {=} 1$ respectively, are effectively identical up to the normalization.
Now consider the associated measurement error covariance matrix ((ref)). The estimate for the $\lambda_S {=} 1$ normalization is
with implied estimated correlation matrix
Note that $\hat{\sigma}^2_{GG}$ is much higher than any of $\hat{\sigma}^2_{SS}$, $\hat{\sigma}^2_{JJ}$, and $\hat{\sigma}^2_{BB}$, potentially due to different indicators using different methods to determine tie points, i.e., reference points of brightness for 100% sea ice and 100% open water. The choice is crucial for accurate measurement of $SIC$ within grid cells. Tie points, moreover, need not be constant, as brightness readings are sensitive to weather effects and atmospheric forcings (Ivan15). Dynamic tie-point calibration is potentially desirable because it can decrease the bias of $SIC$ measurements (ComisoEtAl17). The latest version of the Goddard Bootstrap algorithm, in particular, calibrates tie points daily. One would expect, however, that the bias reduction from dynamic tie-point calibration may come at the cost of potential discontinuities that increase measurement error variance. Our results confirm that conjecture. Our estimate of $\hat{\sigma}^2_{GG}$ is about 13 times that of $\hat{\sigma}^2_{BB}$ (which uses constant tie points).
Alternatively, the estimated measurement error covariance matrix for the $\lambda_G {=} 1$ normalization is
with implied estimated correlation matrix
Now let us move to the transition equation ((ref)). Using the $\lambda_S {=} 1$ normalization we obtain $\hat{\rho} {=} 0.704 ~ [0.041]$ and trend/seasonal parameter estimates ($\hat{a}_i$, $\hat{b}_j$, and $\hat{c}_k$) as reported in the $\lambda_S {=} 1$ columns of Table (ref). Trends for all months are highly significant and downward sloping. The trends for summer months (August-November) display a notable negative, and generally statistically significant, {curvature}, whereas for non-summer months the quadratic trend terms are generally small and statistically insignificant. For the $\lambda_G {=} 1$ normalization we get $\hat{\rho} {=} 0.719~ [0.042]$ and trend/seasonal parameter estimates ($\hat{a}_i$, $\hat{b}_j$, and $\hat{c}_k$) as reported in the $\lambda_G {=} 1$ columns of Table (ref). The $\lambda_S {=} 1$ and $\lambda_G {=} 1$ results are very similar.
In Figures (ref) ($\lambda_S {=} 1$) and (ref) ($\lambda_G {=} 1$) we show optimal latent sea ice extent extractions ($\widehat{SIE^*}$) in black, together with the four raw indicators in color, by month. First consider Figure (ref). Of course $\widehat{SIE^*}(\lambda_S {=} 1)$ is centered on $SIE^S$ due to the $\lambda_S {=} 1$ normalization. Moreover, $\widehat{SIE^*}(\lambda_S {=} 1)$ is always very close -- almost identical -- to $SIE^S$ (and close to $SIE^B$, because $SIE^B$ tends to be very close to $SIE^S$).\footnote{Indeed when making Figure (ref) we added a tiny constant (0.1) to $SIE^S$ to make it easier to distinguish $SIE^S$ from $\widehat{SIE^*}(\lambda_S {=} 1)$.}
Now consider Figure (ref) ($\lambda_G {=} 1$). Due to the different normalization, $\widehat{SIE^*}(\lambda_G {=} 1)$ is centered not on $SIE^S$ but rather on $SIE^G$, so $\widehat{SIE^*}(\lambda_G {=} 1)$ is shifted upward relative to $\widehat{SIE^*}(\lambda_S {=} 1)$. The location of $\widehat{SIE^*}(\lambda_G {=} 1)$ relative to $SIE^G$, moreover, clearly varies by month. In winter months, it tends to be greater than $SIE^G$, whereas in summer months it tends to be less than $SIE^G$. Note in particular that the variation of $\widehat{SIE^*}(\lambda_G {=} 1)$ around $SIE^G$ is noticeably greater than the variation of $\widehat{SIE^*}(\lambda_S {=} 1)$ around $SIE^S$. Clearly $\widehat{SIE^*}(\lambda_G {=} 1)$ is influenced more by movements in other indicators ($SIE^S$, $SIE^J$, and $SIE^B$) than is $\widehat{SIE^*}(\lambda_S {=} 1)$.
The fact that $\widehat{SIE^*}(\lambda_S {=} 1)$ and $\widehat{SIE^*}(\lambda_G {=} 1)$ are different, both in terms of level and variation around the level, limits their usefulness for research focusing on level, as the level depends entirely on identifying assumptions. However, and crucially, in an important sense $\widehat{SIE^*}(\lambda_S {=} 1)$ and $\widehat{SIE^*}(\lambda_G {=} 1)$ are highly similar: The model and identification scheme makes $\widehat{SIE^*}(\lambda_S {=} 1)$ and $\widehat{SIE^*}(\lambda_G {=} 1)$ identical up to a linear transformation. This is clear in Figure (ref), which plots the two competing extracted factors, month-by-month. Regressions of $\widehat{SIE^*}(\lambda_G {=} 1)$ on $\widehat{SIE^*}(\lambda_S {=} 1)$ yield highly-significant intercepts not far from zero, highly-significant slopes near 1.0, and $R^2$ values above 0.999, for each month.
Because $\widehat{SIE^*}(\lambda_S {=} 1)$ and $\widehat{SIE^*}(\lambda_G {=} 1)$ are identical up to a linear transformation, it makes no difference which $\widehat{SIE^*}$ we use for research focused on linear relationships between $SIE^*$ and other aspects of climate (e.g., various radiative forcings). The obvious choice, then, is $\widehat{SIE^*}(\lambda_S {=} 1)$, which is just $SIE^S$ itself, dispensing with the need to estimate the factor model.
We propose a dynamic factor model for four leading Arctic sea ice extent indicators. We estimate the model and use it in conjunction with the Kalman smoother to produce a statistically-optimal combination of the individual indicators, effectively “averaging out" the individual measurement errors. We explore two identification strategies corresponding to two different factor loading normalizations. The corresponding two extracted combined measures (latent factors) are identical up to a linear transformation, so either one can be used to explore relationships of sea ice extent and other variables. Interestingly, however, the extracted factor for one of the normalizations puts all weight on the Sea Ice Index. Hence the Sea Ice Index alone is a statistically optimal “combination" and one can simply use it alone with no loss, dispensing with the need to estimate the factor model. That is, there is no gain from combining the Sea Ice Index with other indicators, confirming and enhancing confidence in the Sea Ice Index and the NASA Team algorithm on which it is based, and similarly lending credibility -- in a competition against very sophisticated opponents -- to the NSIDC's claim that the Sea Ice Index is the “final authoritative SMMR, SSM/I, and SSMIS passive microwave sea ice concentration record" (NSIDC_SII).
\addcontentsline{toc}{section}{References}