3.1 Temporal model for catchment-averaged annual SMB

We fit a generative statistical model for catchment-averaged SMB using an approach adapted from the work of Hu and Castruccio (2021) on other climate fields. We define the n-dimensional vector M(t) to be the catchment-averaged SMB in each of n catchments at time t, and we assume that it can be described by an additive model with a temporal variability vector μ(t) and a noise term vector ϵ(t) of the form

In the example case we present here, the vectors M,μ, and ε have one entry for each of the n=260 catchments at time t, and we evaluate each at a total of m=30 time steps.

The temporal trend μ(t) includes historical mean SMB for each catchment β₀ and the forcing variable f(t) with a linear coefficient β₁. The forcing variable, f(t), can be an external process which causes slow changes in SMB, such as atmospheric temperature (f(t)=T_A(t)) or simply a prescribed dependence on time (e.g., f(t)=t). Finally, Eq. (1b) includes autoregressive terms up to order p contained in the diagonal matrices ${\mathbf{\Phi }}_i,i=\mathrm{1},\mathrm{\dots },p$. The temporal trend as written would thus approximate an autoregressive process of order p, AR(p). Section 3.2 discusses how we identified AR(p) as the best type of temporal model for this application. At this stage, fitting temporal models to annually aggregated time series, we exclude seasonal terms from the temporal trend μ(t); seasonality is incorporated deterministically during the downscaling process described in Sect. 3.6. All stochasticity in this generation technique enters through interannual variability.

The noise term ϵ(t) is assumed to be independent, identically distributed in time, and from an n-dimensional normal distribution with mean of zero and covariance matrix Σ. As we describe in Sect. 3.3, spatial correlations between catchments are captured in ϵ(t).

3.2 Selecting candidate model type and order

We tested several model types in search of the most appropriate way to represent interannual SMB variability in Eq. (1b). Three criteria inform our selection of candidate stochastic model types for temporal SMB variability. First, we would like our candidate temporal models to capture the timescales of variability apparent in the data based on standard statistical methods that are likely to be familiar to glaciologists. Second, we would like our methods to build on existing open-source software such that other researchers can test and apply our work. For that reason, we prioritize models with existing fitting routines in Python or R. Finally, we would like to be able to compare our findings to those of King and Watson (2020) for the Antarctic Ice Sheet, so we prioritize temporal model families that those authors also tested.

These criteria guide our investigation of three common types of temporal models. All temporal models we test belong to the autoregressive fractionally integrated moving-average (ARFIMA) family of models. The first type, order-p autoregressive AR(p) models, is the simplest of the ARFIMA family. They assume that the value of SMB at time t depends linearly on values of SMB at times $(t-\mathrm{1},t-\mathrm{2},\mathrm{\dots },t-p)$. An AR(0) model is equivalent to a white noise model scaled to the data. The second type, ARIMA models of order $(p,d,q)$, includes order-p autoregressive terms applied to a series that has been differenced d times to reach stationarity, as well as dependence on a weighted moving average of the past q residual noise terms. Finally, general ARFIMA models are similar to ARIMA models but allow non-integer values for d, accounting for “long memory” in the time series. To avoid confusion, we henceforth use ARFIMA to refer only to ARFIMA models that do include non-integer differencing d, and we refer to the special cases ARIMA and AR(p) by their own names. King and Watson (2020) tested AR(p) and ARIMA models; they also tested generalized Gauss–Markov models, for which we were unable to find an open-source fitting routine, but which are very similar to ARFIMA models of order $(p,d,\mathrm{0})$.

For each catchment, we estimate β₀ as the 1980–2012 mean and remove it from the series. We then use conditional maximum likelihood (ordinary least squares) to optimize values of β₁ and the remaining parameters of Eq. (1b) associated with each candidate model type (AR, ARIMA, ARFIMA) over a range of orders $(p,d,q)$. We perform the model fitting with built-in functions from the Python package statsmodels v0.12.2 (Seabold and Perktold, 2010): $\mathtt{statsmodels}.\mathtt{tsa}.\mathtt{ar}\mathtt{\_}\mathtt{model}.\mathtt{AutoReg}$ and $\mathtt{statsmodels}.\mathtt{tsa}.\mathtt{arima}.\mathtt{model}.\mathtt{ARIMA}$. In each case, we assume a linear dependence on time, β₁f(t)=β₁t in Eq. (1b). Statsmodels does not include a built-in function to fit ARFIMA models, so we apply fractional differencing following Kuttruf (2019) and subsequently test ARIMA$(p,\mathrm{0},q)$ with the built-in function.

We analyze the Bayesian information criterion (BIC) as returned by the statsmodels built-in function for the temporal models fit to each catchment series. The BIC is given by

where ℓ is the log-likelihood function of the given temporal model on the data, T is the number of observations, and df is the number of degrees of freedom in the generator. Minimizing the BIC balances a maximization of log-likelihood ℓ – the probability that a stochastic generator of this type could have produced the data series from the process model being fit – with a penalty for excess parameters (overfitting). We select the temporal model with the lowest BIC for each catchment for each SMB process model. We analyze the preferred temporal model types across all catchment–model pairs to identify the most suitable class of temporal models. We chose to select the minimum BIC to encourage computationally cheap models with fewer parameters (as in King and Watson, 2020); we note that statsmodels also returns other common metrics of model fit such as the Akaike information criterion, which could be selected by users with other priorities.

Figure 2Autocorrelation (a) and partial autocorrelation (b) functions for the annual SMB computed by seven GrSMBMIP process models for Kangerlussuaq Glacier for 1980–2012. Lollipop markers indicate ACF or PACF values; the blue shaded region indicates the 95 % confidence interval around zero. Values significantly different from zero therefore appear as blue points outside the shaded region.

Download

To decide the range of orders $(p,d,q)$ to test in our model fitting, we use the autocorrelation and partial autocorrelation functions (ACF and PACF, respectively) to target the relevant timescales of variability. In a purely autoregressive process, the number of values significantly different from 0 before the first nonsignificant value in the PACF would indicate the AR order p. In a purely moving-average process, the number of values significantly different from 0 before the first nonsignificant value in the ACF would indicate the MA order q. These metrics cannot be used to determine the order $(p,d,q)$ of a more general ARFIMA process, but we use them as qualitative indicators of an appropriate range for testing. The ACF and PACF values differ per process model and per catchment; an example for the Kangerlussuaq Glacier annual SMB is shown in Fig. 2. In that example, significant autocorrelation is apparent at a lag time of 4 years for several process models, though with several previous values not significantly different from zero; the partial autocorrelation is not significant for any lag shown. Ice-core-derived ACF and PACF show significant values at timescales of up to 5 years, tapering to values not significantly different from 0 at longer timescales (Fig. B1). The combination of evidence from ice cores and from process model ACF and PACF in multiple catchments suggests several lags ≤5 years with significant ACF or PACF. We therefore choose to test values of p and q from 0 to 5. We determine the order of differencing required to reach stationarity, d, using augmented Dickey–Fuller and KPSS tests of stationarity on each catchment time series. Both tests agreed that the de-meaned catchment average time series were stationary, so d=0 should be appropriate, but for completeness we also tested d=0.5 and d=1. Among the range of values $(p,d,q)$ tested, we select the best-fit model as the one with the lowest BIC. We note that comparing the BIC of model fits among temporal models of different orders requires a consistent base dataset and fitting method (for example, the same software package and optimization scheme for all models). We computed the BIC using statsmodels built-in functions, setting the optional argument 𝚑𝚘𝚕𝚍_𝚋𝚊𝚌𝚔=max($p,d,q$) to ensure that lower-order models were fit to training data series of the same length as higher-order models.

Figure 3b shows example best fits for four model types and their BIC (see legend). The best fit to the NHM-SMAP example data in Fig. 3 is white noise with a trend (AR(0)). Both AR(p) models shown have lower BIC than the more complicated ARIMA and ARFIMA models. In this example, the best-fit AR(0) and ARFIMA series capture a trend with little other temporal information, while the AR(5) and ARIMA(1,0,1) series capture larger temporal variability with the expense of added parameters. We note that the models capturing only a trend in Fig. 3b will still generate stochastic series with temporal variability; the distinction is that almost all of the temporal variability in the final generator will come from the spatial noise generation process described in the next section.

In every catchment and SMB model tested (1820 catchment–model pair time series tested), AR models were the most suitable. There were no basin–model pairs where ARIMA or ARFIMA fits were preferred to AR(p) fits. Further, white noise with a trend (AR(0)) was preferred to any higher-order statistical fit for catchment-aggregated SMB in most basins. Each process model had some basins where higher-order AR(p) models were preferred (Fig. A1).

The example we present below allows the order of the autoregressive model fit to vary by basin. Users may decide to keep that flexibility, which adds some complication in storing the model parameters, or they may opt for the simplest AR(0) fit for every basin and allow residual variability to be captured in the spatial noise generation (Sect. 3.3) and downscaling (Sect. 3.6).

Figure 3Catchment mean SMB from seven Greenland-wide models. (a) Time series at monthly scale, as originally presented in the model output data. (b) Time series summed to annual scale, with series from example best-fit stochastic generators overlaid. The Bayesian information criterion for each model's fit to an example process model (NHM-SMAP) is shown in the figure legend. Lower BIC values indicate more preferred models.

Download

3.3 Estimating SMB covariance between catchments

Thus far, we have described a method for fitting and generating time-varying SMB for individual catchments with no correlation beyond the catchment scale. However, SMB over the Greenland Ice Sheet may vary coherently at spatial scales beyond those of single outlet glacier catchments due to large-scale processes in atmospheric circulation (Lenaerts et al., 2019). Motivated by this physical intuition, we introduce spatially informed noise generation.

Following Hu and Castruccio (2021), we construct a matrix of variance Σ for catchment-level noise terms (Eq. 1c above) as

where D is the diagonal matrix of per-catchment standard deviations and C is the spatial correlation matrix among all catchments. Note that this formulation assumes a catchment-specific variance D, so the SMB is assumed to vary differently within each catchment, implicitly accounting for different catchment sizes. The spatial model C is defined on the inter-catchment correlation, which we assume not to depend on the catchment size.

The spatial correlation pattern may differ for different SMB process models, so we construct the matrix of variance for each SMB process model separately. We calculate the empirical correlation matrix $\widehat{\mathbf{C}}$, which is an approximation of C, from the residuals of per-catchment best-fit temporal models described in Sect. 3.2. We save each residual (length n=28, with 5 years held back from the 33-year training set to accommodate consistent fitting of AR orders up to p=5) as a row in a 260×28 matrix R, with one row for each catchment. The empirical correlation matrix $\widehat{\mathbf{C}}$ is then the 260×260 matrix of correlation coefficients of the residuals, which we compute using 𝚗𝚞𝚖𝚙𝚢.𝚌𝚘𝚛𝚛𝚌𝚘𝚎𝚏(𝚁). The empirical correlation matrix computed from ANICE-ITM output is shown in Fig. 4a.

Because the number of catchments we seek to simulate (m=260 for Greenland) is considerably larger than the number of data points used to train individual statistical models (33 years of catchment-aggregated SMB for each catchment), $\widehat{\mathbf{C}}$ is singular. Therefore, we must enforce a sparsity condition to reduce the influence of spurious information. We estimate a sparse correlation matrix Γ using the graphical lasso algorithm described in Friedman et al. (2007).

We apply the 𝙶𝚛𝚊𝚙𝚑𝚒𝚌𝚊𝚕𝙻𝚊𝚜𝚜𝚘𝙲𝚅 function from the Python package scikit-learn v0.24.2 (Pedregosa et al., 2011), which estimates a sparse correlation matrix Γ with the following formulation:

where K is the inverse correlation matrix and α is a positive regularization parameter. Higher values of α lead to sparser resulting matrices Γ. In our implementation, we allow 𝙶𝚛𝚊𝚙𝚑𝚒𝚌𝚊𝚕𝙻𝚊𝚜𝚜𝚘𝙲𝚅 to select the best value of α through cross-validation. Figure 4b shows the sparse correlation matrix resulting from applying this method to ANICE-ITM output.

Figure 4Illustration of inter-catchment covariance in ANICE-ITM output data. (a) The empirical correlation matrix $\widehat{C}$, computed as described in Sect. 3.3. (b) The sparse covariance matrix that results from applying 𝙶𝚛𝚊𝚙𝚑𝚒𝚌𝚊𝚕𝙻𝚊𝚜𝚜𝚘𝙲𝚅 to $\widehat{C}$. (c) The first row of the sparse covariance matrix (line 0 in panel b) translated to map view. Catchment 0 is Umiammakku Isbræ, indicated on the map with UMI and an arrow to its terminus. All panels share the color bar shown below panel (c).

Each row in the sparse correlation matrix Γ represents the correlation of a given catchment with each other catchment. Figure 4c translates the information in the first row of Γ to a map of Greenland. The first row represents catchment 0 in the Mouginot and Rignot (2019) dataset: Umiammakku Isbræ. Umiammakku has the strongest correlation with itself (dark red shading), moderate positive correlation (lighter red shading) with surrounding catchments and a few more distant catchments, and zero or slight negative correlation (light blue shading) with other catchments in Greenland. We note that these correlations are inferred from the process model data – ANICE-ITM output, in Fig. 4 – rather than imposed by physical intuition. As such, the precise structure of the spatial correlation matrix will depend on how the data are aggregated. We would expect slightly different spatial correlations if they were computed with monthly data or using different catchment outlines. Users must also remember that the spatial correlations shown in Fig. 4 are computed on the residuals of temporal model fits, not on the SMB series themselves.

3.4 Forward modeling

Finally, we generate a set of realizations of the forward stochastic generator. Each realization is the sum of an autoregressive component and a draw ϵ(t) from the normal distribution with spatial covariance, as described in Eq. (1a) and Sect. 3.3. We find the Cholesky decomposition Γ=LL^T of the sparse correlation matrix and use the lower triangular component to generate spatially informed noise. The draw ϵ_k(t) for the kth catchment is found by matrix multiplication:

where N_j is a random normal matrix of shape (m,Y) for m the number of catchments, Y is the number of years in the desired time series, and $\widehat{k}$ selects the kth row of the matrix.

We generate realizations of catchment mean SMB for an example catchment: Kangerlussuaq Glacier. Each realization is a single time series of catchment mean SMB with variability described by the stochastic generator. Figure 5 shows 10 realizations of Kangerlussuaq SMB from 1980–2050, with process model training data overlaid in black for 1980–2012. By inspection, the stochastic realizations (blue lines on Fig. 5) have variability of similar amplitude and timescale to the process model series. The 10 realizations, generated in a few seconds on a laptop, fill the expected range of uncertainty in annual SMB. We interpret these stochastic realizations to be an efficiently generated forcing for ensemble simulations of ice sheet change given SMB subject to internal variability.

Figure 5Forward simulation of Kangerlussuaq catchment SMB to 2050, with 1980–2010 mean removed, generated using an AR(4) model with spatially informed noise. The black line shows the results of the process model ANICE-ITM during the period simulated for GrSMBMIP. Blue lines are single realizations of the stochastic generator.

Download

3.5 Nonstationarity

The GrSMBMIP process model historical output we use as our example application did not exhibit nonstationarity, according to the KPSS and augmented Dickey–Fuller tests applied to each output series (Sect. 3.2). We therefore fit stochastic generators that were stationary by construction, assuming the underlying distribution of the data did not change over the period of simulation. We generated stochastic forward simulations as shown in Fig. 5 to illustrate the possibility of generating time series with consistent variability outside the training period. Those simulations fit a linear trend to the training data and assumed that the trend and amplitude of variability remained constant into the future. For scientific applications that study periods of varying climate – for example, glacial–interglacial periods or century-scale climate projections with anthropogenic forcing – it is expected that SMB time series would not be well fit by stationary models (Weirauch et al., 2008; Bintanja et al., 2020).

To fit a stochastic generator to time series with statistics (mean, trend, variance) varying over time, the user could subdivide the training data series into periods with stationary trends and variance. Piecewise linear trends could be computed on the sub-series and each series normalized by its variance to create a “z-score” time series. The stochastic temporal model could then be fit to the z-score time series as described above. The output of the stochastic generator would then be re-scaled by the variance in each period to produce ensembles of nonstationary SMB series. The best choice of break points to subset the training data will depend on the user's priorities and the time period being investigated; we do not pursue z-score re-scaling any further in this example. For more formal discussion of bias correction in the case of climate data whose distribution changes in time, we refer the interested reader to applied statistics literature, e.g., Zhang et al. (2021) and Poppick et al. (2016).

A user generating realizations of SMB at a particular location, or aggregated over some area, could use the method described up to this point. For example, this method could generate realizations of aggregated SMB to support detection of departures from background variability, as in Wouters et al. (2013). The next section describes how to downscale SMB from the catchment annual average to spatially extensive SMB fields at sub-annual timescales.

3.6 Elevation downscaling

To force an ice sheet model, we require a two-dimensional SMB field on the mesh of the model, rather than catchment-aggregated time series. We now apply a spatial and temporal downscaling approach to produce gridded SMB from the stochastically generated series at sub-annual time steps. The downscaling assumes that within each glacier catchment and for a given time of the seasonal cycle, the SMB variation within a catchment can be described by a piecewise linear function with respect to elevation. This downscaling recognizes that, particularly in our Greenland example, there is a strong seasonal cycle in SMB and that the spatial variations of SMB within a glacier catchment are mostly a function of elevation. As shown in Fig. 6, these assumptions are generally quite good for Greenland SMB, and they are reflected in other statistical downscaling approaches that have been previously applied in deterministic frameworks (Hanna et al., 2011; Wilton et al., 2017; Sellevold et al., 2019; Goelzer et al., 2020a). Further, the method generates fields with realistic spatiotemporal variability and elevation dependence, which can be embedded within an ice sheet model (e.g., Verjans et al., 2022) to capture the known feedback between ice sheet surface elevation change and SMB change (Edwards et al., 2014; Lenaerts et al., 2019).

For each point p in a given catchment, we need the surface elevation z(p) used to force the physical SMB model underlying our stochastic generator and the local SMB spatial anomaly,

where A(p,t) is the process model SMB at point p and time t and $\stackrel{\mathrm{‾}}{A}\left(t\right)$ is the catchment mean SMB computed from the same process model at time t. We group all local elevation–anomaly pairs by month – for example, all January values together and all June values together – and fit a piecewise linear mass balance gradient for each month τ:

where c₀ is the minimum SMB and the segment slopes $(c_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}})$ and break points (z₁,z₂) are free parameters optimized by BIC and AIC. In each catchment we thus have 12 functions Λ_τ, one for each month. The monthly mass balance gradients Λ_τ reintroduce seasonal variation. When taken as a function of ice sheet surface elevation z^* that could be evolving in time, they also allow feedback between surface mass balance and dynamically evolving ice sheet geometry.

Example fits for the Kangerlussuaq Glacier catchment, computed from ANICE-ITM output covering 1980–1985, are shown in Fig. 6 (and the same example is shown computed from RACMO data in Fig. A3). The left panels show the spatial anomaly in map view, with the terminus of the glacier to the southeast; the right panels show local SMB departure from the catchment mean as a function of ice surface elevation. The spatial pattern in the example data shows strong departures from the catchment mean throughout the lowest portion of the glacier. January SMB in the lower reaches tends to exceed the catchment mean (blue shading); July SMB in the same area tends to be much below the mean (dark red points). Higher elevations show less pronounced departures from the catchment mean (lighter shading).

Finally, we produce time series of monthly local mass balance a for each grid point p of the kth catchment:

where t is the time in months since the start of the series, M(t) is the annual catchment mean SMB generated by the stochastic generator, Λ_τ is the local SMB spatial anomaly for month τ as defined above, and $z^{*}(\mathit{p},t)$ is the local surface elevation at time t. The same principle could be adapted for training data provided at even finer temporal resolution, though a large training dataset may be needed to capture the relevant variability in sub-monthly SMB.

To illustrate the method, we applied the elevation-based downscaling to estimate local SMB series at two different points, distributed across elevation, in the Kangerlussuaq Glacier catchment. Figure 7b shows those time series. Blue lines are the stochastically generated SMB, downscaled to a single point in space; black lines are the process model output at the grid cell nearest to the selected point. The point represented in the bottom panel is near the terminus and shows large-amplitude seasonal and interannual variations in both the process model and the downscaled stochastic realizations. The stochastic realizations closely track the process model series, while also including interannual variability in winter and summer SMB that differs between realizations. The point in the top panel is in the accumulation area. For that point, the range among the stochastic realizations is wider than the apparent variability in the process model series. The seasonal cycle has approximately correct amplitude. We interpret the variability in the catchment-averaged SMB to be dominated by large-amplitude variation near the terminus (Figs. 6 and A3), which is then reflected in the stochastic generator fit to the process model series. We further discuss this overestimate of accumulation zone interannual variability in the next section.

Figure 6SMB downscaling in time and space, shown for two example months (rows): January and July. Left panels: SMB spatial anomaly (difference from catchment mean) for each point within the Kangerlussuaq Glacier catchment, based on the ANICE-ITM model contribution to GrSMBMIP. Right panels: SMB lapse rate with elevation for each month, deduced from the anomaly fields shown. Colored data points represent individual pixel values for 1980–1985. The color map is consistent across panels.

Download

Figure 7(a) Location of example points (blue stars and numbers) in the Kangerlussuaq Glacier catchment. Circle markers show SMBMIP grid points, colored by the seasonal range in mass balance at that location, computed as local mass balance in December minus local mass balance in July. (b) SMB time series scaled from catchment mean down to local (single grid point) values. The series in the upper panel are scaled to point 1, in the accumulation area; the series in the lower panel are scaled to point 2, near the terminus. As in previous plots, the black lines in each series are process model output (ANICE-ITM for the example case) and the blue lines are stochastic realizations. Series share x and y axes.

Download