The global mean sea level rise predicted by its causative budget components during 2018 – 2050

ABSTRACT This study establishes a predictive empirical model, the first of its kind, which is innately a cause-and-effect representation of the observed global mean sea level over time. The model uses the trends of its observed budget components, including the temporal variations in mass of glaciers, Greenland, Antarctic ice sheets, terrestrial water storage, steric effect, and astronomical forcings of luni-solar origin, as independent variables in representing global mean sea level anomalies, namely, the dependent variable. The model parameters are estimated using monthly globally averaged satellite altimetry measurements and the yearly rates of the budget components during 1993–2018 as a priori information. Prospective monthly global mean sea level anomalies are then quantified and tabulated together with their root mean square error of prediction at 5% significance level for the period 2018–2050.


Introduction
Ocean thermal expansion and halosteric contraction, ice sheet, mountain glacier and ice cap mass balance, changes in terrestrial surface and groundwater storage, anthropogenic impoundment of waters in dams, and to a much lesser importance, permafrost degradation, snow and atmospheric water vapour change, are major geophysical contributors to the present-day Global Mean Sea Level, (GMSL), variations (WCRP Global Sea Level Budget Group, 2018). They are the components of what is called the GMSL budget that enable quantification of the budget closure/misclosure under the conservation of water mass on Earth.
A recent report by the WCRP Global Sea Level Budget Group (2018) elaborates exhaustively the GMSL budget for the period 1993 to 2018 and its components using recent literature on sea level science including the US Fourth Climate Assessment Report (USGCRP, 2018), and the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (IPCC, 2014).
The practice of evaluating GMSL budget is to quantify its geophysical contributors (budget components) from various sources and to calculate the misclosure of the budget reflecting deviation from zero during a period for which all time series data are available. WCRP Global Sea Level Budget Group (2018) calculated this misclosure using the linear rates of changes (trends) of the budget components. A recent study by İz et al. (2019) adjusted the misclosure of the reported linear trends of the budget components and their uncertainties during the period of 2005-2015 under the GMSL budget closure constraint. An alternative approach by H.B and Shum (2020a) achieved the budget closure by adjusting the linear trends (velocities) of the GMSL budget components simultaneously during the same period. Another study by H.B and Shum (2020a) adjusted GMSL budget during the period 2005-2015 with emphasis on closing the budget on a year-by-year basis as opposed to using linear trends of the GMSL budget components. The adjustment also accounts for the effect of snow, water vapour, and permafrost mass components as a lump sum. A more recent study by Iz (2022) detected periodicities of luni-solar forcings in the yearly misclosures of the global mean sea level budget during 1900-2018 as another GMSL budget component of astronomical origin.
The synergy among the budget components suggests their use to establish a cause-and-effect GMSL model. The observed variations in GMSL rise as the dependent variable can be modelled as a function of the remaining observed budget components as independent variables. The model can be used to assess present-day GMSL rise as well as for its prediction, which will be demonstrated in this study In the following sections, this cause-and-effect representation of the GMSL is presented. The subsequent sections start with the description of the observed time series of the GMSL budget components during 1993-2018, followed by a preliminary assessment of the GMSL budget closure. The observed GMSL anomalies are then modelled as a function of the budget components' velocities (trends) and evaluated using globally averaged monthly Satellite Altimetry (SA) data and the velocities of the remaining budget components. Prospective monthly global mean sea level anomalies are then quantified and tabulated together with their root mean square error of prediction at five percent significance level for the period 2018-2050.

Data on yearly time series for the global mean sea level budget components
The GMSL budget data used in this study was compiled by Frederikse et al. (2020), which relies on various earlier works cited in their manuscript and will not be repeated here. Figure 2 through Figure 6 display the yearly variations of the GMSL budget components for the period 1990-2018 for a total of 26 yearly averaged values for each one of the GMSL budget components together with their 1σ confidence intervals (CI). The yearly averaged time series are expressed in mm of equivalent sea level during this period.

The GMSL sea level budget
Following the narrative in WCRP Global Sea Level Budget Group (2018), İz et al. (2019), the GMSL budget closure is formulated as follows, GMSLðtÞ À GMSLðtÞ STERIC À MðtÞ GLACIERS À À MðtÞ GREENLAND À MðtÞ ANTARCTICA À MðtÞ TWS ¼ 0 where the components of the above expression represent time-dependent mass contributors (barystatic) to global mean sea level GMSL t ð Þ by the mass of glaciers,MðtÞ GLACIERS , Greenland, MðtÞ GREENLAND , Antarctic ice sheets, MðtÞ ANTARCTICA , and terrestrial water storage (TWS), MðtÞ TWS , respectively. The steric component, GMSLðtÞ STERIC , refers to the contributions of ocean thermal expansion and salinity to sea-level change, respectively (ibid.).
The GMSL budget misclosure, W t , can be quantified year by year using annually and globally averaged sea level height anomalies of each component, h t , i.e.
where superscript t is the epoch of the yearly averaged anomaly. The yearly misclosures are the lump sum effect of the potential biases and noise in the yearly averaged time series of the budget components as well as any other unknown contributors to the GMSL budget ( Figure 7). The same figure also displays the standard errors (SE) of the misclosures on a yearly basis calculated from the uncertainties of the budget components using variance propagation through Eqn.
(2). Figure 8 reveals that the variability of the misclosures is large with Root Sum of Squares (RSS) of misclosures RSS(W t ) ¼ 209mm: Their yearly fluctuations over time overlap with the contributions of either Antarctic, Glaciers, or Greenland time series or their sums. In the following section, the yearly fluctuations of the budget components are modelled and estimated together with their uncertainties using Generalised Least Squares (GLS).

Yearly velocities of GMSL budget components
The time progressions of yearly averaged observed anomalies for each of the budget components exhibit preponderantly a linear trend (Figure 1-6). They can be modelled as follows, The above observation equations include two unknown parameters for each component; their datums, h t0 COMPONENT that are defined in the middle of the series, t 0 with Δt Δt À t 0 , and their linear velocities, v COMPONENT . The random variations that are not explained by the linear trend parameters are represented by ε t COMPONENT for which their uncertainties are available as shown in Figure  The parameters and their standard errors for each of the budget components are estimated using GLS solutions together with their Adjusted R 2 (Table 1). Except TWS, the trend models explain over 90 percent of the variations in the yearly anomalies for the budget components. The estimated velocities enable the quantification of the budget velocity misclosure through the time derivative of Eqn. (2). The misclosure is not statistically significant at 2σ level. Nonetheless, its estimated magnitude is large, which can be attributed to the limited power of each model using only 26 yearly averaged values during 1993-2018 with 24 degrees of freedom. A recent study by this investigator (H.B, 2022) using yearly data W t reported a misclosure rate much smaller, 0.08 ± 0.00 mm/yr. using longer time series . Another contributor to the large misclosure is the unmodeled episodic variations in TWS, which are more pronounced during the satellite altimetry era. In addition, the same earlier study by Iz (2022) has also detected signatures of luni-solar forcings with periods 18.6 and 11.1 yr. in GMSL budget misclosures, which are present in GMSL observations. The presence of such variations of astronomical origin was also shown at globally distributed tide gauge (TG) measurements (Iz 2014). The omission of their effects biases the GMSL velocity estimates resulting in a larger misclosure rate.
The velocity estimates and their uncertainties for the budget components are the prerequisites for the      causal model for the GMSL predictions, which is discussed in the following section.

Modelling GMSL anomalies as a function of its causal budget components
The water conservation in the climate system enabling the time-dependent sea level budget equation given by Eqn. (1) can be expressed for a GMSL model as a function of the anomalies of its budget components as, The above expression also includes periodic variations due to the lunisolar forcing, 18.6 and 11.1 yr. periods, respectively, detected in TG and GMSL anomalies (Iz, 2014(Iz, , 2022. Note that the presence of a uniform acceleration is ambiguous during the SA era (1993-2022) as reported by İz and Shum (2020b).
In this model, the velocities of the budget components cannot be estimated because they are co-linear, sharing the same Δt. Nonetheless, the linear rate of change of the GMSL together with their standard errors can still be estimated introducing the trends of its budget components shown in Table 1 as a priori information as follows, Introducing the above stochastic constraints on the budget components' velocities as a priori information enables a proper GLS solution to the rank deficient observation equations based on Eqn.(4). This is the topic of the following sections.

Statistical model and the GMLS during SA era
The observation equations for Eqn. (4) and the a priori estimated linear trends of budget component Eqn. (5) can be cast together in the following matrix/vector format, In this expression, y GMSL is the n GMSL � 1 vector of GMSL observations. The n υ � 1 vector of observations y υ refers to the previously estimated velocities.
The A GMSL is the n GMSL design matrix for the SA component. The n υ matrix A υ , belongs to the velocities of the remaining budget components and x v is the m � 1 vector of unknown parameters that appear in Eqn. (4).
In the above representation, there are two sets of disturbances. The first set are the disturbances for the component velocities, ε v , excluding GMSL series, that are estimated in the previous sections GLS solutions. The estimated velocities are assumed to be independent of each other. 1 Their uncertainties σ v will be represented by a 5 × 5 diagonal matrix Σ υ .
The other set of disturbances, ε GMSL , are related to the SA altimetry data. They are represented by an n GMSL � 1 vector, ε GMSL . Although this investigator did not have access to high-resolution time series for all the budget components, globally averaged monthly satellite altimetry measurements shown in Figure 9 were made accessible online by GSFC (2017). The availability of 309 monthly averages, compared to their 26 yearly averages, provides more statistical power and precision in the GLS solution. Meanwhile, the lack of monthly measurements for the remaining budget components is not an impediment for the formulation, which is designed to use the a priori velocity estimates instead of observations. The velocity estimates were robust despite the limited number of yearly values thanks to their nearly linear progressions in time.  The marginal effect of increased statistical power, however, is costly. Monthly GMSL measurements exhibit high first order autocorrelations, AR (1), as reported by İz and Shum (2020b). The autoregressive property of the disturbances, ε GMSL t ; may be attributed to the averaging/smoothing of the globally averaged SA data. An AR(1) process may also be caused by unmodeled systematic variations of physical origin. Therefore, modelling AR(1) effect is an integral part of the model. High autocorrelation of disturbances causes random excursions in the residual series that may be triggered by what is called in econometrics as random innovations/ shocks and has the propensity to confound the interpretation of the results as shown in İz and Shum (2020b).
The disturbances of monthly averaged GMSL measurements can be represented as, In this expression, À 1 � ρ � 1 is the unknown autocorrelation coefficient of the AR(1) process. The stochastic processes for the random noise u GMSL t , have the following assumed distributional properties, They are identically and independently distributed, i.i. d. This expression together with Eqn. (7) gives, The corresponding V/C matrix denoted by Σ GMSL is given by, For the monthly averaged globally averaged SA measurements, the estimated AR(1) correlation coefficient is ρ ffi 0:94 (İz & Shum, 2020b) with an assumed a priori variance of unit weight σ 2 ¼ 1.
The following expressions are obtained by manipulating the Normal Equations of Eqn.(6) of the model parameters b x GMSL and their m V/C matrix, Σx, The parameters estimated using the above relationships are listed in Table 2. The adjusted velocity    estimates for the budget components are not significantly different than those obtained from separate a priori solutions. The amplitudes of the periodic components are also large in magnitude. Nonetheless, the nodal amplitude is a borderline value at α ¼ 0:05% significance level. The first order autocorrelated residuals, AR(1), and the i.i.d distributed random noise are shown in Figure  10 and Figure 11 respectively. The model explains 98.3% of the variation in monthly GMSL anomalies, which is an impressively good fit thanks to the inclusion of the periodic GMSL variations of Luni-solar origin into the model. Therefore, the model is well suited for GMSL predictions, which is the topic of the following section.

GMSL prediction
This section starts with a summary presentation of the statistical method to predict monthly GMSL anomalies known as Best Linear Unbiased Prediction (BLUP). Detailed information about the method can be found in Toutenburg (1982). An application of the method to sea level studies appears in . Subsequently, GMSL rise will be predicted for the period 2021-2050 together with their prediction intervals.
Denoting the subscript ' * ' to reference predicted values, n � is the number of observations to be predicted. If these observations are generated by the same linear model given by Eq. (7), then they can be expressed as follows, In this expression, y � is the n � � 1 vector of observations to be interpolated, extrapolated, predicted, and A � is the n � design matrix. The statistical properties of the disturbances are carried out over the predictions with ε � , which is the n � � 1 vector of disturbances. Furthermore, it is assumed that the past and future disturbances are uncorrelated, i.e. E εε T The predictor of y � is given by its GLS estimator (Toutenburg, 1982, p. 135 In this expression, ŷ � is the n � � 1 vector of predicted observations. Further manipulation of the above expression with the above statistical properties of the model gives, where b x is an unbiased estimate of x.
In practice however, one needs information about how much predicted values would differ from their actual values, i.e. b y � À y � . To this end, the corresponding statistic is given by its mean square error of prediction, MSEP, which is defined as follows, Prediction Intervals, PI, are different than the Confidence Intervals, CI, where the latter involves the expected value of the predicted value, i.e. its theoretical mean as if there are there are multitude of realisations of anomalies at a given epoch, which is of course not relevant in practice compared to predicting the actual realisation of the GMSL anomalies.
If the disturbances follow a normal distribution, i.e. εÑ 0; Σ y À � , then the corresponding PI, for each predicted value are computed from In these expressions, z 1À α;DF is the predefined zscore, DF is the degrees of freedom, and α is the significance level. The error of the predicted observation is denoted by its root mean square error of prediction, rmsep. The GMSL rise for the period 2018-2050 was predicted using the above BLUP. Figure 12 displays their progression with noticeable contributions of the conflated luni-solar effects with natural variations at global scale. The predicted anomalies are also tabulated in Table 3. The predictions are precise despite the use of rmsep instead of CI that accounts the variability in the observed GMSL rather than its expected value.

Conclusion
The model presented in this study is a cause-and-effect representation of the GMSL and effectively explains 98% of the monthly variations in globally observed SA measurements with an optimal number of parameters. Hence, it lends itself naturally predicting GMSL variations (İz & Shum, 2020b).
In this study, no comparison was made with any projections because the underlying assumptions of the predictions of this study are binary and should not be interpreted as possibilities contrary to projections. The tabulated prospective GMSL anomalies together with their rmsep are to serve as two mutually exclusive scenarios. If the smooth and steady contributions of the budget components continue in the future, which is likely to be the case based on their history , then the predictions will be accurate, precise, and useful for risk assessments at global scale during 2018-2050. If not, then the predicted GMSL anomalies will provide quantitative baseline values (observedpredicted), to assess the impact of the interventions put in place to mitigate climate change at global scale during this period. The resolution of the predictions can be further improved by introducing the impact of episodic large-scale events such as ENSO (El Niño and Southern Oscillation) on sea-level rise (Cazenave et al., 2012) into the model.

Notes
1. This assumption is correct in the sense that the budget components are not observationally correlated. Nonetheless, there may still be a cross-correlation because one component may be influencing the other with a lag, such as steric changes impacting melting of Greenland. 2. Alternative derivations are given in Toutenburg (1982) for the case of correlated disturbances if this relationship is known.