Understanding and forecasting the climatic observational record (COR) is critical for a variety of tasks -notably, the planning of future mitigation and adaptation infrastructure and investment. The COR has many diverse elements, such as temperature (means, extremes, variability, etc. ), extreme events (storms, floods, droughts, wildfires, etc. ), and other physical attributes (sea level, cloud coverage, glacial extent, etc.). These elements differ importantly at a spatial level, and the Arctic region has emerged as a salient focal point for investigations of ongoing climate change. In particular, Arctic melting is a conspicuous effect of climate change, which has warmed the Arctic two to three times faster than the global average. The melting Arctic also promotes additional future climate change, in part due to an "albedo" feedback amplification loop as more reflective Arctic sea ice is replaced by darker open water. Along with global effects, a melting Arctic also has pervasive regional consequences. For example, in the Arctic region, less sea ice will reduce the cost of shipping, promote the extraction of natural resources, and expand tourism. Together, at a global and regional level, a melting Arctic will result in widespread and transformative impacts, opportunities, and risks (Alvarez et al., 2020;Lynch et al., 2022;Brock and Miller, 2023). Large-scale climate models have been the source of much of the forward-looking global and Arctic climate analysis. These models attempt to capture the fundamental physical drivers of the earth's climate with much spatial and temporal detail. Such structural models are invaluable for understanding climate variation; however, from a forecasting perspective, climate models have been less adept and, in particular, have largely underestimated the amount of lost sea ice in recent decades (e.g., Stroeve et al., 2012;Rosenblum and Eisenman, 2017). As a practical complementary approach, statistical projections from parsimonious dynamic time series representations often produce forecasts that are at least as accurate as detailed structural models (Diebold and Rudebusch, 2022). COR forecasting can be facilitated by using reduced-form econometric/statistical methods, which are designed to summarize the relevant data patterns (correlations) regardless of the deep underlying causal climatic mechanisms. This reduced-form agnosticism is valuable for COR forecasting because the granular climate dynamics are incompletely understood -witness the fifty or so different CMIP6 climate models and their widely-differing projections (e.g., Notz and SIMIP, 2020). Specifically, for forecasting Arctic sea ice, there is already some evidence that smallscale econometric/statistical models can have some success (Wang et al., 2016;Diebold and Rudebusch, 2022). In this paper, we therefore attempt to push forward the econometric/statistical analysis of the long-run future evolution of Arctic sea ice. Of particular interest are probability assessments of the timing of an ice-free Arctic, an outcome with vital economic and climate consequences (Jahn et al., 2016). A key contribution of our analysis is that we take a multivariate approach to the various measures of Arctic sea ice. Specifically, we consider three aspects of pan-Arctic sea ice: surface coverage (measured in two ways -area and extent, to be defined shortly), thickness, and volume. In theory, these three measures of sea ice should be perfectly related. Volume would equal the product of surface coverage and thickness, so modeling any two would imply a specification for the third. In practice, however, each sea-ice indicator is based on a blend of pure observations and model interpretations and interpolations. As a result, measurement differences and errors introduce discrepancies between the various indicators and drive wedges between them, even if the amount of measurement error has generally diminished over time as observational and modeling methods have improved. Thus, the various measures of sea ice are not perfectly related, and it may be beneficial to consider all of them. We characterize the various measures of September Arctic sea ice as time series processes and examine the likely future path toward seasonal ice-free conditions with continued global warming. September is the month with the least seasonal ice, and we focus on projecting the timing of the first ice-free September, exploring both a literally ice-free Arctic (IFA) and an effectively or nearly ice-free Arctic (NIFA). 1 We model sea-ice area, extent, thickness, and volume, and we implement multivariate models that constrain these measures to reach zero simultaneously. We anchor our analysis on area, sequentially considering pairwise blends of area with each of the other indicators. The area-extent pair is essentially based on only observed data; the area-thickness pair integrates the two key indicators characterizing sea ice; and the area-volume pair is closest to an all-indicator analysis that incorporates area, thickness, and extent. We implement our constrained joint forecasting procedure in "carbon-trend" models that relate sea ice to atmospheric carbon dioxide (CO 2 ) concentration (with this time series denoted as CO2C ). This allows us to condition our forecasts on various carbon scenarios. We also check our carbon-trend results against similarly-constrained "time-trend" models that relate sea ice simply to time. Finally, we verify the consistency of the two approaches, and we use them to make probabilistic forecasts of the arrival years of NIFA and IFA. There is of course a large literature on measurement and modeling of declining Arctic sea-ice extent and area, and insightful overviews include Stroeve et al. (2012), Shalina et al. (2020a), andFox-Kemper et al. (2022). For analyses of the measurement and modeling of declining thickness and volume, see, e.g., Stroeve and Notz (2018) and Shalina et al. (2020b). We will have more to say throughout the paper about the literature on Arctic sea ice and how our results relate to it. For now, however, we emphasize again that our multivariate modeling of area, extent, thickness, and volume -respecting in particular the fact that if and when they vanish, they must vanish simultaneously -sets our analysis apart. In particular, the present paper is quite different from our earlier research, Diebold and Rudebusch (2022), which compares structural climate models to reduced-form statistical models based on polynomial trends. The present paper avoids any statistical vs structural forecast model comparisons.
Instead, it dives more deeply into the statistical analysis and modeling, incorporating much more information to form statistical forecasts: (1) more sea-ice indicators, (2) bivariate models with equal-IFA constraints imposed, and (3) conditioning on a key geophysical covariate (CO 2 concentration), instead of just a generic polynomial trend on a single sea-ice indicator. We proceed as follows. In section 2, we introduce our four sea-ice indicators. In section 3, we model and forecast the indicators using regressions on atmospheric CO 2 concentration ("carbon trends"), producing point and interval forecasts. In section 4, we provide a major robustness check by modeling and forecasting using direct time trend regressions ("time trends"), assessing the consistency of the carbon-trend and time-trend forecasting approaches. In section 5, we use both approaches to provide full probabilistic forecasts of the two objects of ultimate interest: the arrival years of IFA and NIFA. In section 6, we provide a second major robustness check by redoing most of the carbon-trend analysis using an alternative carbon measure, cumulative anthropogenic CO 2 emissions (CO2E ), assessing the consistency of the two carbon-trend versions. We conclude in section 7. We consider four key aspects of Arctic sea ice: surface coverage (of which there are two measures, area and extent), thickness, and volume. Our analysis focuses on the seasonal minima for these indicators, specifically, the September monthly average at a pan-Arctic scale.2

Arctic sea-ice surface coverage has been well measured using satellite-based passive microwave sensing since late 1978. The Arctic's surface is divided into a grid of cells, and for each cell, satellite sensors measure the amount of reflected solar radiation. This reflected energy varies from cell to cell depending on the relative cell coverage of sea ice and open water, as the latter absorbs more radiation. The resulting measured "brightness" is converted into a fractional measure of sea-ice "concentration" for cell i at time t, denoted c i,t . The individual c i,t are then aggregated into two standard, somewhat different measures of overall sea-ice coverage. The first of these, "sea-ice area" (SIA), sets measured concentration in each cell to zero if it is below .15 and retains the satellite-recorded value otherwise:

All grid cell areas are then multiplied by c SIA i,t , and their weighted sum over the entire Arctic is denoted SIA t . As a result, SIA is the measured coverage of all Arctic grid cells with at least 15 percent sea-ice concentration. The second measure of coverage, "sea-ice extent" (SIE), also sets measured concentration below .15 to zero but otherwise adjusts it upward to full coverage, or 1:

In this case, all grid cell areas are multiplied by c SIE i,t , and their weighted sum is total Arctic SIE t . Thus, SIE is the total area of all of the ocean grid cells that are measured to have at least 15% sea-ice coverage. The up-rounding in the construction of SIE serves as a bias correction to offset satellite sensors' tendency to mistake shallow pools of sea-ice surface melt for open sea.3 There are costs and benefits of such a bias correction, and both SIA and SIE are widely employed in research and forecasting. We will examine both, using September monthly average SIA and SIE data from 1979 to 2021 from the National Snow and Ice Data Center (NSIDC). The up-rounding in the construction of SIE also implies that by definition SIE ≥ SIA, and indeed, in the historical data, SIE has been substantially higher than SIA. For our analysis, the crucial definitional constraint is that these two measures of sea-ice coverage can only equal zero under the same condition, that is, when all Arctic grid cells have no more than fifteen percent coverage (c i,t ≤ .15 for all i ). In that case, the two measures will be equivalent, so SIA = SIE = 0. This constraint that both indicators must reach zero at the same time will be the crucial timing restriction that allows us to construct a joint projection of pan-Arctic September SIA and SIE. Our next sea-ice indicator is thickness, denoted SIT . Sea-ice thickness -the distance between the ocean underneath and any snow coverage or air above the ice -is challenging to measure but allows a fuller understanding of sea-ice conditions. 4 In the past, data on sea-ice thickness have been obtained via a variety of methods including submarine records, ocean floor buoys, bore holes, helicopter and sled surveys, and satellite measures, but these sources have varying accuracy and availability by time and location. A popular alternative to these disparate thickness measures is produced by the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS) from the University of Washington's Polar Science Center (Zhang and Rothrock, 2003). The PIOMAS thickness indicator incorporates (or assimilates) data such as the NSIDC sea-ice concentration measure as well as atmospheric information such as wind speed and direction, surface air temperature, and cloud cover. However, the PIOMAS data are importantly driven by an ocean and sea-ice physical model that characterizes how sea-ice thickness evolves over location and time in response to dynamic and thermodynamic forcing and mechanical redistribution such as ridging. That is, the PIOMAS thickness data are effectively a model construct, although they have been largely validated against the available observational data. For example, Selyuzhenok et al. (2020), building on Schweiger et al. (2011), report that the spatial pattern of PIOMAS ice thickness agrees well with observations derived from in situ and satellite readings, and Labe et al. (2018) note that decadal trends in thickness appear realistically reproduced. Surveys that reach a similar conclusion include Leppäranta et al. (2020, section 8.4) and Shalina et al. (2020b, section 5.6). Our final sea-ice indicator is volume, SIV , which is also supplied via the model-based PI-OMAS measurement. Specifically, across a grid of cells, this series combines the PIOMAS thickness and the sea-ice coverage data to produce sea-ice volume cell by cell.
Of course, in theory, volume, which represents the entire mass of sea ice, should be the most comprehensive indicator of Arctic sea ice. On the other hand, SIV , which is produced as the product of surface coverage and thickness, inherits measurement error from both of those series and appears to be subject to the most uncertainty in measurement. As noted above, our joint estimation approach relies on a simultaneous zero-ice constraint on more than one measure of sea ice. The two sea-ice coverage series, SIA and SIE, share data sources and definitions that effectively hard-wire or force such a simultaneous ice-free Arctic. A similar exact relationship is less assured for thickness and volume as these depend on model-based imputations. Furthermore, SIT and SIV do not impose the 15% cutoff used in SIA and SIE, which potentially could lead to later ice-free dates for thickness and volume (though our empirical results below suggest the opposite timing). However, given that all four indicators depend on the amount of sea-ice coverage, it is likely that SIA, SIE, SIT , and SIV will all reach zero very closely in time if not simultaneously, and we will assess the imposition of that constraint in the estimation of sea-ice dynamics. 3 Linear Carbon-Trend Models, Fits, and Forecasts

Here we explore the empirical relationships between the four Arctic sea-ice indicators and atmospheric CO 2 concentration. We then use these representations to forecast the likely progression of a melting Arctic and, specifically, the years of first IFA and NIFA.5

We measure CO 2 using readings on atmospheric concentration rather than cumulative emissions, which is a widely used alternative (e.g., Notz and Stroeve, 2016). Atmospheric CO 2 concentration has notable advantages in this regard. It is a direct measure of the amount of heat-trapping gasses in the atmosphere -the source of the melting Arctic. In particular, unlike cumulative emissions, atmospheric concentration can account for time-varying CO 2 absorption rates, so as the efficacy of natural terrestrial and oceanic carbon sinks change over time, concentration can better capture the extent to which greenhouse gases are driving climate change. In addition, atmospheric CO 2 concentration, which is measured by direct air sampling, is subject to little measurement error. By contrast, fossil CO 2 emissions are indirectly calculated from energy, fuel use, and cement production data, while emissions from land-use changes are based on deforestation estimates. Our atmospheric CO 2 concentration series (measured in parts per million, ppm) is shown in Figure 1. The historical data  are from the NOAA Global Monitoring Laboratory, as measured at the Mauna Loa Observatory in Hawaii. The four projected scenarios (2022-2100, shaded) are from the SSP Public Database. 6 These are the standard SSP sce-narios denoted SSP1 2.6 (which has very low emissions or radiative forcing), SSP2 4.5 (low emissions), SSP3 7.0 (medium), and SSP5 8.5 (high).7 They are distinguished by the amount of progress made in reducing further concentration increases, as shown by the variation in their upward curvature. We focus on the SSP3 7.0 scenario as a baseline but also show that our results are qualitatively robust to use of the low and high emissions scenarios. Various researchers -including Johannessen (2008), Notz and Stroeve (2016), and Stroeve and Notz (2018) -have identified a linear empirical relationship between observed Arctic sea-ice area and CO 2 .8 This linear carbon-trend relationship, which fits remarkably well in the observed data, can be expressed as:

where SIA t is sea-ice area, CO2C t is atmospheric CO 2 concentration, and ε t represents deviations from the linear fit. The regression intercept, a, calibrates the level of sea-ice coverage, and the slope, b, provides a measure of the Arctic-carbon response. A negative value for b captures the diminishing coverage of Arctic sea ice as greenhouse gases accumulate in the atmosphere.9 Equation (3) also can be used to characterize the occurrence of first IFA in terms of CO2C. In particular, the level of atmospheric CO 2 concentration consistent with the first occurrence of an ice-free Arctic (SIA t = 0) occurs when CO2C t = -a/b. The "univariate" representation (3) has typically been used to fit and forecast a single sea-ice indicator -either sea-ice area or extent. We generalize this analysis to consider thickness and volume as well. More importantly, we also introduce a joint modeling strategy that allows two or more sea-ice indicators to be modeled together. Such a representation combines indicators using a simultaneous zero-ice constraint. For example, for two indicators, the "bivariate" linear carbon trend model with a "simultaneous first IFA" constraint is given by

where x t and y t are two sea-ice indicators. This regression model has jointly constrained slopes and intercepts so that x t = y t = 0 at the same level of CO2C t and can be estimated via non-linear least squares. Like the univariate equation ( 3), the bivariate system (4) can be used to characterize a first (simultaneous) IFA in terms of CO2C -that is, a concentration level for a seasonal ice-free Arctic. For our bivariate empirical analysis with (4), we always set x t to SIA t and y t to one of the other sea-ice indicators (SIE t , SIT t , or SIV t ).10 That is, our focus here is on the three bivariate combinations: SIA+SIE, SIA+SIT , and SIA+SIV . Each indicator pair brings additional information to the estimation of the first IFA year by incorporating additional indicators with SIA, which provides a common benchmark for comparing results. In addition, each pair has idiosyncratic characteristics that provide useful differentiation. The SIA+SIE pair is the only one that is based on data that are essentially directly observed, and these two indicators also share a clear-cut zero-ice constraint by definition. However, this pair is limited to only sea-ice coverage, while the other combinations account for thickness as well. The SIA+SIT pair integrates the two key elements necessary for measuring sea ice, but again, SIT is a model-processed measure that has greater measurement error.
The SIA+SIV pair combines area with the most comprehensive measure of sea ice, but again, at a likely cost of less measurement accuracy. It is worth highlighting three aspects of our approach. The first is our use of pure-trend bivariate models with no lags. In our view, such models are well suited to our dataset. Our use of annual (September) sea-ice indicators renders greater dynamic generality unnecessary. For example, the first autocorrelation of detrended September SIE is only 0.06. Also, degrees of freedom are scarce in our sample of roughly 40 annual observations. An alternative vector autoregression (VAR) specification would gain little from adding dynamics at the cost of a profligate parameterization; a 4-variable VAR(3), for example, would have more than 50 parameters. The second aspect of our approach is use of fixed deterministic carbon trends with no breaks. We intentionally skirt the unit-root vs. trend-break minefield. That debate remains unresolved in macroeconomics, for example, after more than half a century of research, as standard data simply are not very informative about the hypotheses of interest. 11 Nevertheless, one can interpret our carbon trend regressions as cointegrating regressions under a unit-root worldview, and our focus on fixed carbon trends as a simple and accurate empirical approximation, as per the fit described below in Figure 2. The third aspect of our approach is use of linear carbon trends, which might be violated by episodes of prolonged natural variation or by nonlinear shifts or tipping points in the ice-carbon linkage. Our stochastic modeling methodology accounts for natural variation, although its specification may be imperfect. Natural variation is emphasized in Miller and Nam (2020) as a potential source of a future slowdown in Arctic ice melting. However, the "hiatus" in global warming of 1998-2012 had a minor effect on the pace of Arctic melting, which may reflect its potential status as measurement error. Alternatively, as sea-ice melt depends importantly on polar conditions, Arctic amplification of temperature may mask such temporary global fluctuations. Regarding nonlinearities, following the earlier literature, we are comfortable with linear carbon trends as a suitable empirical approximation, especially in light of the excellent fit evident in Figure 2 and discussed below. Although many potential factors could produce a nonlinear ice-carbon relationship, linearity turns out to be a remarkably robust approximation for the entire satellite-era observational record, regardless of whether CO 2 is measured as cumulative emissions or atmospheric concentration (e.g., Eisenman and Wettlaufer, 2009;Tietsche et al., 2011;Winton, 2011;Notz and Stroeve, 2016;Notz and SIMIP, 2020). Moreover, the empirical climate-science evidence on linearity of the ice-carbon relationship continues to accumulate; see for example the recent work of Bennedsen et al. (2023), who work in a very general nonlinear state space environment, allowing for nonlinear carbon trend, but find no need for nonlinearity. The authoritative recent review of Wang et al. (2023) accurately summarizes the state of current knowledge:

The literature points toward a linear, predictable response of summer sea ice 11 Classic macroeconomic work on stochastic vs. deterministic trend, with or without trend breaks, includes Nelson and Plosser (1982) and Perron (1989), among many others; Stock (1994) provides a thorough survey. Climate economics work grappling with the same issues includes Kaufmann et al. (2006) and Estrada and Perron (2019), among many others. extent in response to greenhouse gas emissions, despite large inter-annual variability in weather patterns, rather than an abrupt transition to seasonally ice-free conditions (monthly average sea ice area of less than 1 million km 2) consistent with tipping point behavior (p. 44). That is, tipping points are unlikely to affect Arctic sea ice carbon-trend linearity, even if they affect other aspects of climate change. Of course, moving forward, the effects of diminished albedo, greater storm-generated wave action, or permafrost melt may notably accelerate sea-ice loss, or, alternatively, greater cloud cover in a warming planet may retard that loss. There are several subtleties regarding a potential melting permafrost tipping point (e.g., Vaks et al., 2020). Permafrost degradation can occur as gradual top down thaw or as abrupt collapse of thawing soil. Such melting permafrost releases both methane and CO 2 . The exact timing and proportion of these releases is uncertain, but to the extent that CO 2 and non-CO 2 greenhouse gases (GHGs) remain in a similar proportion, our linear carbon trends will still hold approximately. The CO 2 emissions time trajectory and the sea-ice time trends, however, would change. A permafrost melt that quickens the pace of warming strengthens our thesis that Arctic sea ice is melting more quickly than widely projected. In closing this subsection, we note that the theoretical climate science literature has also gravitated toward linearity, as summarized by Hillebrand (2023), precisely because of the compelling empirical evidence discussed above.12 Important linear modeling of carbon cycle land and ocean sinks, for example, includes not only the classic work of Bacastow and Keeling (1973) but also the more-recent work of Raupach (2013). Similarly, important linear modeling of atmosphere/ocean energy balance ranges from the classic work of Budyko (1969) and Sellers (1969), to more recent work that introduces two-layer energy-balance models with linear heat transfer between atmosphere/ocean surface and deep-ocean layers, such as Gregory (2000) and Held et al. (2010). The top panel shows univariate models for each indicator, and the bottom panel shows constrained bivariate models for SIA paired with SIE, SIT , or SIV that impose a simultaneous first IFA. We report parameter estimates, R 2 's, and projected concentrations and years at first IFA and first NIFA (the latter for SIA only). The first IFA and NIFA years are based on the time path under the SSP3 7.0 scenario. Throughout, standard errors appear in parentheses. See text for details. Each panel graphs a pair of sea-ice measures against the atmospheric CO 2 concentration. Unconstrained univariate results appear in the left panels, and constrained common-IFA bivariate results for various pairs of indicators appear in the right panels.
The black irregular lines show SIA (measured in 10 6 km 2 ) versus CO2C, and the blue, red, and green irregular lines are for SIE, SIT , and SIV (measured in 10 6 km 2 , m, and 10 3 km 3 , respectively). These various standard units of measurement are such that the pairs of SIA, SIE, SIT , and SIV can be plotted conveniently together in the various panels of Figure 2, but this convenience is irrelevant for our results. Trend regression coefficient estimates depend on the choice of measurement units, but, crucially, the implied estimates of our key objects of interest, namely, the extrapolated years of the first ice-free Arctic, do not depend on these units. In Figure 2, linear carbon trends, again colored black, blue, red, and green for SIA, SIE, SIT , and SIV , are fitted to the data in the unshaded sample. Their fit illustrates how well the linear regressions capture the relationships between the sea-ice indicators and atmospheric CO 2 concentration. For each measure, the historical data cluster quite tightly around the fitted carbon trends. This remarkably robust linearity has been noted for sea-ice area in the literature (e.g., Notz and Stroeve, 2016), and here we generalize this result to other measures. In percentage terms, Arctic sea-ice coverage and thickness have trended downward at a similar rate, with SIE, SIA, and SIT falling by about 50% over the sample. When combined, these declines account for the 75% drop in SIV over the sample. In the shaded regions of Figure 2, the carbon trends are extrapolated beyond the historical sample until the arrival of an ice-free Arctic. The unconstrained univariate column reveals that projected carbon at first IFA is higher for SIE than for SIA, slightly lower for SIT , and notably lower for SIV . Comparing the univariate and constrained bivariate columns of Figure 2, one sees that imposition of the simultaneous first-IFA constraint changes the slopes of the fitted and extrapolated carbon trends. This is consistent with the fact that for each pair of indicators, the common constrained first-IFA carbon levels lie between the two unconstrained first-IFA carbon levels. The qualitative insights from Figure 2 are detailed in Table 1. The estimated b's are all negative and highly statistically significant, and all R 2 values are above .80. The univariate projected first-IFA concentrations for SIA, SIE, SIT , and SIV are 488 ppm, 514 ppm, 484 ppm, and 437 ppm, respectively. When SIA is modeled jointly with SIE, SIT , and SIV , the bivariate constrained projected first-IFA concentrations are 509 ppm, 488 ppm, and 445 ppm. That is, the bivariate first-IFA carbon levels are higher or lower than the univariate estimate depending on which additional indicator is blended with SIA. Note that the bivariate constrained blending of SIA with the other indicators reduces the range of the projected first-IFA concentrations; the univariate range is 514 ppm -437 ppm = 77 ppm, whereas the constrained bivariate range is 509 ppm -445 ppm = 64 ppm. While our bivariate common first-IFA constraint imposes the joint occurrence of a zero-ice event, the literature has often focused on the occurrence of a nearly ice-free Arctic, or NIFA, defined as an SIA or SIE level of 1 million km 2 or less. This definition reflects a view that certain Northern coastal regions, notably, the Canadian Arctic Archipelago, may retain small amounts of landfast sea ice even after the open Arctic Sea is ice free (Wang and Overland, 2009). However, the resistance of such landfast ice to melting remains an open research issue, and there is much uncertainty about how resilient such a circumscribed "Last Ice" refuge will be to further warming (Cooley et al., 2020;Mudryk et al., 2021;Schweiger et al., 2021). Reflecting this uncertainty, while we also consider the usual definition of a nearly ice-free Arctic (10 6 km 2 ) for SIA and SIE, we do not account for any possible structural shift in near-zero sea-ice dynamics. Accordingly, Table 1 reports the projected first-NIFA concentrations for SIA. For SIT and SIV , there are no similar definitions of nearly ice-free concentrations in the literature, so we concentrate on the occurrence of a first NIFA only in SIA. Of course, any first-NIFA concentration is lower than the associated first-IFA level, but the precise difference depends on the indicator(s) examined and model used (univariate vs. bivariate). The univariate projected SIA first-NIFA concentration is 464 ppm. It almost matches the constrained bivariate SIA+SIT SIA first-NIFA projection at 460 ppm, but is still much higher than that for SIA+SIV , and lower than that for SIA+SIE. The range of the first-NIFA concentration projections is 480 ppm -430 ppm = 50 ppm. The "carbon budget" is a popular concept used to quantify the remaining amount of carbon that can be emitted until some climate event is reached. Analogously, we consider a carbon budget as the difference between the current atmospheric CO 2 concentration and the concentration associated with the first occurrence of IFA or NIFA. Given a 2021 CO2C value of 417 ppm and an estimated concentration of 464 ppm at first NIFA (measured under the SIA univariate model), our estimate of the Arctic carbon budget available until the likely occurrence of a nearly ice-free Arctic is 47 ppm. This concentration carbon budget is slightly higher or lower for the other model pairs. Finally, although we have thus far discussed only point forecasts, we also provide simulationbased 90% interval forecasts in Figure 2, shown as shaded regions around the point fore-casts. 13 The intervals are generally quite narrow, despite their appropriate widening as the projection horizon lengthens. Note also that our models are constrained in expectation only; that is, their conditional mean functions are constrained but their shock variances are not, so that some intervals are always wider than others (compare, e.g., the SIA and SIV intervals).
There are of course many under-appreciated or unmodeled aspects of Arctic sea-ice diminution that may affect the accuracy of both our point and interval forecasts. Some of these may speed or slow sea-ice loss. They include exceptional internal variation, which may result in a warming hiatus (Miller and Nam, 2020) that slows loss, and tipping points associated with methane release from permafrost melting (Chylek et al., 2022), which could promote an even steeper decline in Arctic sea ice. The monotonically increasing CO 2 concentrations in the historical data and SSP projected scenarios of Figure 1 allow us to translate carbon-trend representations like (3) (in ice-carbon space) into time-trend representations (in ice-time space). In particular, the one-to-one SSP scenario mappings allow us to translate the first IFA/NIFA carbon level projections into projected first IFA/NIFA years, conditional on the SSP scenario. Moreover, the curvature of the SSP scenarios produces a curvature in the time trends implied by the linear carbon trends. Because SIA is linearly related to CO2C and the SSP scenarios used to convert carbon trends into time trends are nonlinear, the implied SIA time trend will be nonlinear. For example, with the CO 2 concentration in the SSP 7.0 scenario increasing over time at an increasing rate, the implied sea-ice time trend will be decreasing at an increasing rate. Of course, translating the carbon-trend forecasts into ice-time space requires a second layer of conditioning relative to the carbon-trend forecasts in ice-carbon space insofar as inverse SSP scenarios are used to infer ice-time profiles, including NIFA and IFA times. That is, the timing of the first NIFA predictions from a given estimated carbon-trend model is conditional on the assumption that future emissions follow a certain SSP path. The implied fitted and forecasted nonlinear time trends are shown in Figure 3, the format of which parallels that of Figure 2. The implied time trends are obtained by running the fitted carbon trends of Figure 2 through the (inverse) concentration schedule SSP3 7.0 of Figure 1. Accordingly, each concentration in ice-carbon space yields a corresponding time in ice-time space. For example, as discussed earlier, our univariate estimate of concentration at first NIFA is 464 ppm. As seen by following the upper dashed horizontal line in Figure 1 (and as recorded in Table 1), a concentration of 464 ppm translates under SSP3 7.0 into a first NIFA year of 2033. Table 1 contains implied projected first-IFA and first-NIFA years. The univariate projected SIA first-IFA year is 2039 and the corresponding bivariate constrained projected first-IFA years are 2044, 2039, and 2028, depending on whether SIA is modeled jointly with SIE, SIT , or SIV , respectively. The range of univariate projected first-IFA years is 2045 -2026 = 19 years, and the constrained bivariate range is 2044 -2028 = 16 years. Table 1 also reports projected SIA first-NIFA years. The univariate projected SIA first-NIFA year is 2033. This date is close to the constrained bivariate projected SIA+SIT first-NIFA year and ranges between those for SIA+SIV and SIA+SIE, respectively. The range of the projected first-NIFA years is 2037 -2024 = 13 years. Across all univariate and bivariate specifications, the average SIA first-NIFA data falls between 2031 (average bivariate) and 2033 (univariate), so the occurrence of a nearly ice-free Arctic is clearly projected to occur fairly soon. Finally, as we will describe in detail below, our first-IFA and first-NIFA dates are robust to use of low, medium, or high concentration growth scenarios. We have emphasized the medium (SSP3 7.0) scenario as a plausible baseline, but low (SSP2 4.5) and high (SSP5 8.5) scenarios provide similar dates, because first NIFA occurs relatively rapidly before the carbon levels in the various scenarios can diverge. This can be seen in Figure 1, where the horizontal CO 2 concentration line at 464 ppm cuts the medium and high SSP schedules at nearly-identical times, and the low SSP schedule just a few years later. That is, our results indicate that global adoption of a somewhat lower emissions path would likely delay the arrival of a seasonally ice-free Arctic by only a few years. 4 Quadratic Time-Trend Models, Fits, and Forecasts

In the previous section, we focused on the connection between higher atmospheric CO 2 concentration and melting Arctic ice, implicitly recognizing that higher greenhouse gas concentrations raise global surface temperatures and lead to melting ice. Remarkably, a simple linear regression appeared adequate, based on the past forty years, for forecasting CO2C at first IFA and NIFA. We then used a leading SSP concentration scenario to convert the forecasted sea-ice paths from ice-carbon space to ice-time space. That is, we ultimately forecasted the path and pattern of SIA diminution over time, and in particular, the years of first IFA and first NIFA. The nonlinear SSP path converted the linear carbon trend into a nonlinear implied time trend, decreasing at an increasing rate. As a first and major robustness check, we now compare the nonlinear time trends obtained from the ice-carbon relationship to deterministic polynomial trends fit directly to SIA, as in Diebold and Rudebusch (2022). In particular, we allow for the relevant ice-time nonlinearity by allowing for quadratic time effects. In the univariate approach, we use

where SIA t is sea-ice area, T IM E is a time trend dummy, and ϵ t represents deviations from trend. Similarly, in the bivariate approach we use

where x t is SIA t and y t is one of the other sea-ice indicators (SIE t , SIT t , SIV t ). Linearity of course emerges as a special case. In our flagship constrained bivariate approach, we move to the equivalent root form of the quadratic system,

where the roots of the x trend are 1/γ x and α x /β x , and the roots of the y trend are 1/γ y and α y /β y .
Imposing the common first-IFA constraint (i.e., both trends reach zero simultaneously) amounts to constraining the two quadratics to have a common root, γ x = γ y = γ, which can immediately be done in a joint estimation of the two quadratic trends. 14The fitted direct quadratic time trends are shown in Figure 4. The point forecast paths are generally similar to those of the carbon-implied time trends in Figure 3. The interval forecast paths are wider, however, than for carbon-implied time trends. This is to be expected because the carbon-implied time trends of Figure 3 assume the amount of curvature via an assumed SSP path on which they condition, whereas the direct time trends of Figure 4 estimate the curvature, which is challenging. The extra parameter estimation uncertainty then translates into wider projection intervals. Additionally, note that only part of a time trend is constrained by imposing identical IFA times (one root is left free), and, as a result, different indicators approach IFA at different speeds. The different curvatures produce different asymmetries in the IFA intervals, with flatter trends (like SIA) promoting asymmetry and steeper trends (like SIV ) promoting symmetry. Put differently, the carbon-implied trends provide forecasts conditional on a linear carbontrend model and a particular assumed SSP scenario, whereas the directly-fitted quadratic time trends do not, instead (implicitly) conditioning simply on continued "business as usual". Neither approach is necessarily "better" than the other, and indeed they both address the same question ("What is the likely future path of Arctic sea ice? "), but the perspective and methods used to answer the question differ across the approaches. Given the common question but different approaches to answering it, it is certainly of interest to compare the answers, and reassuring to see that they cohere closely under a leading carbon scenario. Parameter estimates for both the unconstrained univariate and constrained bivariate direct quadratic trends are highly statistically significant, and all R 2 's are above 0.80. The key bivariate projected first-NIFA years are 2043, 2037, and 2036, for SIA+SIE, SIA+SIT , and SIA+SIV , respectively. The projected first-NIFA year range is 2043 -2036 = 7 years. The mean projected first-NIFA year is 2039, about 8 years after NIFA projected by the carbon-trends estimation. All told, the point forecasts of first-NIFA years from the carbon-trend models and the direct time-trend models agree quite closely, in terms of NIFA arriving very soon, but the carbon-trend models all feature somewhat earlier NIFA. We now proceed to provide a more complete characterization of the uncertainty surrounding the first-NIFA point forecasts, via a full accounting of first-NIFA probability distributions. Here we maintain our focus on first-NIFA years, but we build up approximations to their full probability distributions. Unlike our earlier first-NIFA point projections, which are based simply on extrapolated deterministic trends, here we account for random variation, which lets us learn not only about central tendencies of the first-NIFA distributions (mean, median, mode), but also other moments and related statistics (standard deviation, skewness, kurtosis, left-and right-tail percentiles, etc.) Results appear in  1. For example, for the leading baseline case, carbon trend SIA+SIT using the mid-range concentration scenario SSP3 7.0, the distribution median is 2031, a bit earlier than the 2033 reported in Table 1. The slightly earlier central tendencies of NIFA arrival in Table 3 occur because, when simulating sample pathways, the addition of stochastic shocks to a declining deterministic trend, particularly one concave to the origin (such as ours), might shorten the time necessary to achieve a given threshold. The random variation raises the possibility of a transitory fall below the threshold even before the deterministic trend achieves it on a sustained basis. 15 According to the simulation, stochastic shocks are more likely to result in early NIFA years (79.85%) than later NIFA years (9.51%). The discrepancy between the timing of the deterministic first NIFA year and the mean (or median or mode) of the stochastic first NIFA year depends on the curvature of the trend and the variance of the shocks. A steeper slope, a more concave trend curve, or a smaller error variance will reduce the timing difference. The last two columns of Table 2 characterize the right tails of the NIFA arrival distributions. Consider again the baseline case: carbon trend SIA+SIT with medium concentration scenario SSP3 7.0. The "80%" column of Table 2 reports an eighty percent chance of first NIFA by 2033. Even the most extreme late arrival (moving up and down the 80% column) is fairly soon (2057). Returning to the baseline case but increasing the confidence level from 80% all the way to 95% increases the NIFA arrival year by only one year -from 2033 to 2034. Alternatively, maintaining the 80% confidence level but moving to the optimistic lower concentration scenario SSP2 4.5 again postpones NIFA arrival by (only) five years. In Section 4, we checked the robustness of our carbon-trend results by comparing them to directly fitted time trends. In this section, we implement a second major robustness check by providing and comparing carbon-trend results when carbon is measured as cumulative CO 2 emissions since 1850 (CO2E ) rather than atmospheric CO 2 concentration (CO2C ). 16The emissions-based carbon-trend results are quite similar to the earlier results based on concentration. A linear carbon trend for emissions appears to fit well for each measure of Arctic sea ice. The projected cumulative emissions at first SIA IFA is generally on the order of 2842 to 4179 gt. The corresponding CO2E for first NIFA (univariate SIA) is 3275 gt, which when compared to the current cumulative emissions -projected for 2021 under SSP3 7.0 -of about 2499 gt puts the carbon budget at around 776 gt for reaching first NIFA. The implied projected first-IFA and first-NIFA years from the emissions-based carbontrend regressions using the SSP3 7.0 scenario are also quite similar to the earlier results based on concentration.
The univariate projected SIA first-IFA year is 2043, and the corresponding bivariate constrained projected first-IFA years are 2049, 2043, and 2031, depending on whether SIA is modeled jointly with SIE, SIT , or SIV , respectively. These are several years later than the corresponding first-IFA years obtained using atmospheric concentration.
The implied emissions-based first-NIFA probability distributions appear in Table 3. They are nearly identical in shape to the earlier-reported concentration-based first-NIFA distributions, with one difference: the emissions-based distribution is shifted slightly rightward relative to the concentration-based distribution (i.e., shifted toward later first-NIFA years). Using SSP3 7.0 and SIA+SIT , for example, produces an emissions-based median first NIFA of 2033, for example, in contrast to our earlier concentration-based 2031. The right-shifting of the emissions-based distribution makes its left tail thinner, so its 80th and 95th percentiles increase by even more relative to the concentration-based distribution. Moving from concentration to emissions increases the 80th percentile from 2033 to 2036 (again using SSP3 7.0 and SIA+SIT ) and increases the 95th percentile from 2034 to 2038. We have constructed projections of September Arctic sea ice using a variety of specifications. We prefer bivariate rather than univariate projections, because they enlarge the underlying information set in two ways. First, the bivariate projections blend the information in the history of SIA with that in the histories of other sea-ice indicators, such as SIE, SIT , and SIV , which may improve projection accuracy. Second, they provide regularization by enabling us to constrain first IFA to coincide across sea-ice indicators, thereby imposing appropriate geophysical constraints on otherwise flexible statistical models. Importantly, our analysis produces full probability distribution projections of Arctic sea ice and, in particular, distributions of first-NIFA years based on constrained bivariate models. Our empirical results reconcile linear carbon trends with quadratic time trends and provide constrained multivariate support for an early probabilistic September disappearance of sea ice. Notz and Stroeve (2016) and others have also considered carbon regressions of the kind we use, but again, we generalize them to a multivariate setting, impose equal-IFA constraints, and construct density forecasts. Using atmospheric CO 2 concentration data, the median of our preferred first-NIFA distribution is 2031, with 80% probability by 2033 and 95% probability by 2034. Using cumulative CO 2 emissions data, the median is 2033, with 80% probability by 2036 and 95% probability by 2038. Therefore, it is very likely that the September Arctic will be nearly ice-free at some time before the middle of the next decade. Our results are largely robust to the modeling strategy (univariate or multivariate, carbon trend or direct time trend), the sea-ice indicator variable blended with SIA (SIE, SIT , SIV ), concentration path assumptions (high SSP5 8.5, medium SSP3 7.0, low SSP2 4.5), CO 2 measures (atmospheric concentration, cumulative emissions), and statistical confidence levels (80%, 95%). Note in particular that although both this paper and Diebold and Rudebusch (2022) agree that Arctic sea ice will disappear on a seasonal basis within a couple of decades, the two papers use very different information sets and methods, and the fact that they agree is itself an important result. Hence, unless we somehow manage to follow something like the extremely low SSP1 2.6 concentration path, there is no escaping the sharp bottom line: The Arctic will become seasonally ice free very soon -most likely by the midto late 2030's. In closing, we note that although this paper is about refining statistical sea-ice projections, not about comparing statistical and climate-model projections (unlike our earlier paper, Diebold and Rudebusch (2022)), it nevertheless provides indirect evidence supporting the speculations of, for example, Stroeve et al. (2007) and Stroeve et al. (2012) that the slow decline in Arctic sea ice projected by climate models reflects an underestimation of the effects of GHGs. Our forecasts are indeed generally quite different from those of typical climate models; we predict an earlier first NIFA (in the 2030's), whereas the latest CMIP6 coupled-climate models tend to obtain first NIFA around mid-century. Hence, as regards future research, it will be of interest to assess whether the linear relationship between Arctic sea ice and carbon dioxide, which we have documented in the observational record and used to make probabilistic assessments of NIFA arrival, is also present in simulated paths from large-scale dynamical climate models. If so, it will be of great interest to assess whether the "b" parameters (as defined by our equation ( 3)) embedded in various climate models are of sufficient magnitude to achieve consistency with the rapid observed historical Arctic sea-ice decline, and our projected continued rapid decline. Diebold and Rudebusch (2023) take initial steps in that direction. For the sample from 1979 to 2021, we consider four measures of September Arctic sea ice: area, extent, thickness, and volume. Area data are from the National Snow and Ice Data Center (Sea Ice Index monthly dataset, Version 3, Dataset ID G02135, https://nsidc.org/data/G02135/versions/3). Area is from 30.98N, measured in 10 6 km 2 . The satellites miss the "pole hole", and the published SIA data exclude it, implicitly assuming that the pole hole has zero ice (c=0). A better approximation is to assume that the pole hole has full ice (c=1). Hence we first fill the pole hole by adding 1.19 × 10 6 km 2 to SIA from sample start through July 1987, 0.31 × 10 6 km 2 from September 1987 through December 2007, and 0.029 × 10 6 km 2 from January 2008 to present.17

Extent data are from the National Snow and Ice Data Center (Sea Ice Index monthly dataset, Version 3, Dataset ID G02135, https://nsidc.org/data/G02135/versions/3). Extent is from 30.98N, measured in 10 6 km 2 . Thickness and volume date are from PIOMAS at the Polar Science Center, http://psc. apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly/data/. Thickness is from 49N, measured in m, where greater than 0.15 m. Volume is from 49N, measured in 10 3 km 3 . Volume data are published monthly. Thickness data are published on daily (file PIOMAS.thick.daily.1979.2020.Current.v2.1-2.dat), and we transform the daily data into monthly averages. The CO 2 measures are annual September atmospheric concentration and cumulative emissions. The historical

Cumulative CO 2 emissions data are taken from Rogelj et al. (2021) and the corresponding IPCC report (Allan et al., 2022). The actual dataset (Rogelj et al., 2021) is at https: //data.ceda.ac.uk/badc/ar6_wg1/data/spm/spm_10/v20210809.
From the latter, we use the following files:

Our bivariate bootstrap projection interval simulation procedure, of which our univariate procedure is a special case, proceeds as follows at bootstrap replication i. First we estimate the model and collect residuals. Second, we draw pairs of those residuals, to preserve their cross-correlation, and we add them back to the in-sample fitted values. Third, we re-estimate the model on this synthetic data and use the re-estimated model to project out-of-sample. Finally, we draw pairs of residuals and add them to the forecasted paths to get the simulated path realizations. We repeat this for i = 1, 2, ..., 1000. We then sort the simulated path values at each date and obtain the 90% intervals shown in Figures 3 and4, with left and right endpoints given by the fifth and ninety-fifth percentiles, respectively, of the simulated path values. Our bivariate bootstrap first NIFA simulation proceeds identically, with one extra step. At the end, once the simulated path realizations are in hand, we calculate the corresponding simulated first NIFA realizations, NIFA i , i = 1, 2, ..., 1000. We then obtain the first NIFA densities shown in Figures 3 and4 by applying a nonparametric kernel density estimator to NIFA i , i = 1, 2, ..., 1000. One might naturally wonder whether our bootstrap should account for dependence, which it does not. Fortunately, however, the residuals in our annual models of September Arctic sea ice display negligible dependence. Indeed ρ(1) = 0.06 for both linear carbon-trend residuals and quadratic time-trend residuals. (The key is to recognize that although for obvious reasons September shocks are likely to be highly correlated with those of earlier or later months like August or October, they are much less likely to be highly correlated with those of earlier or later Septembers.) Notes

We follow the usual convention and define NIFA as sea-ice coverage of less than 10 6 km

as described below. See Appendix A for data sources and details. See https://www.ncdc.noaa.gov/monitoring-references/dyk/arctic-ice for details. See, for example,Bunzel et al. (2018),Chevallier et al. (2017),Shalina et al. (2020b),Notz andSIMIP (2020), andZygmuntowska et al. (2014). Although other greenhouse gases, such as methane, are also responsible for global warming, CO 2 is by far the most important. Moreover, exploration of a broad aggregation of greenhouse gases that takes into account the varying global warming impacts and duration of different gases -so-called "CO 2 equivalent"produced results similar to those reported below. See Appendix A for details. See, for example,Arias (2022). Some researchers measure CO 2 as atmospheric concentration and others as cumulative emissions. As discussed above, we employ the former; however, as shown in section 6, our results are robust to use of cumulative emissions. As regards possible endogeneity of CO2C, note that we are not attempting to estimate a structural coefficient. Rather, we are simply estimating a reduced-form best linear predictor under quadratic loss, for which least squares in levels is always consistent, regardless of the endogeneity/exogeneity status of the regressors, and regardless of whether the series are stationary, deterministically trending, stochastically trending (integrated), or cointegrated(Sims et al., 1990). Of course, the bivariate first-IFA-constrained modeling approach can be generalized to more than two indicators. Indeed, we have applied it to 3-and 4-variable sets of our sea-ice indicators and obtained similar empirical results to those reported below. Subtleties also arise with 3-and 4-variable systems, such as additional potentially relevant restrictions among indicators, as when, for example, SIV is at least the approximate product of SIE and SIT . The basic theoretical idea underlying linearity is the insight that climate system response fluxes may be approximated linearly, while still reproducing many aspects of the historical climate record -including but not at all limited to carbon-trend linearity, as in the classic work ofRaupach (2013). The interval forecasts are obtained from 1000 bootstrap simulations, and they account for both parameter estimation uncertainty and disturbance uncertainty. For details, see Appendix B. That is, in the joint regression, a i = α i , b i = -β i -α i γ, and c i = β i γ, for i = x, y. Consider, for example, a gentle deterministic trend that slowly approaches zero from above. This series will continue to run just above zero for a while before hitting zero. Adding random shocks to that trend, on the other hand, will almost certainly push the series below zero considerably sooner. We call this the "principle of stochastic precedence". The emissions data are based onRogelj et al. (2021) and served as input to the IPCC AR6 report(Allan et al., 2022). The emissions-based regressions use the historical CO 2 data from 1979 to 2019 and two additional years of data from the SSP3 7.0 scenario to complete the sample. See Appendix A for details. See the NSIDC Sea Ice Index Version 3 User Guide, https://nsidc.org/data/G02135/versions/3, p. 17. See: https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt. The data file was created on 11/05/2021 at 10:28:56. The data were retrieved on 11/23/2021. Our calculations are based on the de-seasonalized data column.