Probabilistic calibration of glacier projections using data assimilation

The 2010 Eyjafjallajökull volcano eruption had serious consequences to civil aviation. This has initiated a lot of research on volcanic ash forecasting in recent years. For forecasting the volcanic ash transport after eruption onset, a volcanic ash transport and diffusion model (VATDM) needs to be run with Eruption Source Parameters (ESPs) such as plume height and mass eruption rate as input, and with data assimilation techniques to continuously improve the forecast. Reliable and accurate ash measurements are crucial for providing successful ash clouds advices. In the first<br/>phase of this research work, simulated aircraft-based volcanic ash measurements, will be assimilated into a transport model to identify the potential benefit of this kind of observations in an assimilation system. The results show that assimilating aircraft-based measurements can improve the state of ash clouds, and can provide an improved forecast. We also show that for an advice on the aeroplane flying level, aircraft-based measurements should preferably be taken at this level. Furthermore<br/>it is shown that in order to make an acceptable advice for aviation decision makers, accurate knowledge about uncertainties of ESPs and measurements is of great importance.<br/><br/>The forecast accuracy of distal volcanic ash clouds is important for providing valid aviation advice during volcanic ash eruptions. However, because the distal part of a volcanic ash plume is far from the volcano, the influence of eruption information on this part becomes rather indirect and uncertain, resulting in inaccurate volcanic ash forecasts in these distal areas. In this thesis, we use real-life aircraft in situ observations, measured in the North-West part of Germany during the 2010 Eyjafjallajökull eruption, in an ensemble-based data assimilation system to investigate the potential improvement on the forecast accuracy with regard to the distal volcanic ash plume. We show that the error of the analyzed volcanic ash state can be significantly reduced by assimilating real-life in situ measurements. After assimilation, it is shown that the model-based aviation advice for Germany, the Netherlands and Luxembourg can be improved. We suggest that with suitable aircrafts measuring once per day across the distal volcanic ash plume, the description and prediction of volcanic ash clouds in these areas can be improved significantly.<br/><br/>Among the data assimilation approaches, the ensemble Kalman filter (EnKF) is a well-known and popular method. A proper covariance localization strategy in the analysis step of EnKF is essential for reducing spurious covariances caused by the finite ensemble size, as shown for this application for assimilation of aircraft in situ measurements. After analyzing the characteristics of the physical forecast error covariances, we present a two-way tracking approach to define the localization matrix<br/>for covariance localization. The result shows that the Two-way-tracking Localized EnKF (TL-EnKF) effectively maintains the correctly specified physical covariances and largely reduces the spurious ones. The computational cost of TL-EnKF is also evaluated and is shown to be advantageous for both serial and parallel implementations. Compared to the commonly used distance-based covariance localization, the two-way tracking approach is shown to be more suitable. In addition, the covariance inflation approach is verified as an additional improvement to TL-EnKF to achieve more accurate results.<br/><br/>A timely prediction requires that the computations of the data assimilation system can be performed quickly (at least than the Wall-clock). We therefore investigate strategies for accelerating the data assimilation algorithm. Based on evaluations of the computational time, the analysis step of the assimilation turns out to be the most expensive part. After a study on the characteristics of the ensemble ash state, we propose a mask-state algorithm which records the sparsity information of the full ensemble state matrix and transforms the full matrix into a relatively small one. This will reduce the computational cost in the analysis step. Experimental results show the mask-state algorithm significantly speeds up the analysis step. Subsequently, the total amount of computing time for volcanic ash data assimilation is reduced to an acceptable level. The mask-state algorithm is generic and thus can be embedded in any ensemble-based data assimilation framework. Moreover, ensemble-based data assimilation with the mask-state algorithm is promising and flexible, because it implements exactly the standard data assimilation without any approximation and it realizes the satisfying performance without any change of the full model.<br/><br/>Infrared satellite measurements of volcanic ash mass loadings are often used as<br/>input observations for the assimilation scheme. However, these satellite-retrieved<br/>data are often two-dimensional (2D), and cannot easily be combined with a three-dimensional (3D) volcanic ash model to improve the volcanic ash state. By integrating available data including ash mass loadings, cloud top heights and thickness information, we propose a satellite observational operator (SOO) that translates satellite-retrieved 2D volcanic ash mass loadings to 3D concentrations at the top layer of the ash cloud. Ensemble-based data assimilation is used to assimilate the extracted measurements of ash concentrations. The results show that satellite data assimilation can force the volcanic ash state to match the satellite observations, and that it improves the forecast of the ash state. Comparison with highly accurate aircraft in situ measurements shows that the effective duration of the improved volcanic ash forecasts is about a half day.

Bayesian state estimation of a dynamical system utilising a stream of noisy measurements is important in many geophysical and engineering applications. In these cases, nonlinearities, high (or infinite) dimensionality of the state space, and sparse observations pose key challenges for deriving efficient and accurate data assimilation (DA) techniques. A number of DA algorithms used commonly in practice, such as the Ensemble Kalman Filter (EnKF) or Ensemble Square Root Kalman Filter (EnSRKF), suffer from serious drawbacks such as catastrophic filter divergence and filter instability. The analysis of stability and accuracy of these DA schemes has thus far focused either on finite-dimensional dynamics, or on complete observations in case of infinite-dimensional dissipative systems. We develop a unified framework for the analysis of several well-known and empirically efficient DA techniques derived from various Gaussian approximations of the Bayesian filtering schemes for geophysical-type dissipative dynamics with quadratic nonlinearities. We establish rigorous results on (time-asymptotic) accuracy and stability of these algorithms with general covariance and observation operators. The accuracy and stability results for EnKF and EnSRKF for dissipative PDEs are, to the best of our knowledge, completely new in this general setting. It turns out that a hitherto unexploited cancellation property involving the ensemble covariance and observation operators and the concept of covariance localization in conjunction with covariance inflation play a pivotal role in the accuracy and stability for EnKF and EnSRKF. Our approach also elucidates the links, via determining functionals, between the approximate-Bayesian and control-theoretic approaches to DA. We consider the ‘model’ dynamics governed by the two-dimensional incompressible Navier–Stokes equations (NSEs) and observations given by noisy measurements of averaged volume elements or spectral/modal observations of the velocity field. In this setup, several continuous-time DA techniques, namely the so-called 3DVar, EnKF and EnSRKF reduce to a stochastically forced NSEs. For the first time, we derive conditions for accuracy and stability of EnKF and EnSRKF. The derived bounds are given for the limit supremum of the expected value of the L 2 norm and of the H 1 Sobolev norm of the difference between the approximating solution and the actual solution as the time tends to infinity. Moreover, our analysis reveals an interplay between the resolution of the observations (roughly, the ‘richness’ of the observation space) associated with the observation operator underlying the DA algorithms and covariance inflation and localization which are employed in practice for improved filter performance.

Abstract. The growing availability of glacier observations poses a challenge for models to integrate this heterogeneous information in a dynamically consistent way. At the same time, estimates of current glacier volume and area remain uncertain, as many global inventories and thickness datasets date back to the early 2000s. We present the Open Global Glacier Data Assimilation Framework, named AGILE, a time-dependent variational method inspired by 4D-Var data assimilation. AGILE is built on a reimplementation of the OGGM flowline glacier evolution model in PyTorch, enabling full differentiability through automatic differentiation (AD). We test AGILE v0.1 in a series of idealized experiments designed to reflect common initialization and calibration scenarios in global glacier modeling. The goal is to recover glacier bed topography and distributed ice volume in 2020 through transient calibration, based on dynamical simulations starting in 1980. In these experiments, we assume a perfectly known mass balance and fixed ice dynamics parameters. While this setup simplifies real-world complexity, it allows us to isolate and evaluate the core functionality of the approach. Our results show that AGILE efficiently optimizes multiple control variables by leveraging AD-derived gradients, requiring only a few iterations to substantially improve upon initial guesses. We also examine the potential to reconstruct earlier glacier states (e.g., in 1980) without direct observations and find that this is fundamentally limited because glacier dynamics are governed by a diffusion equation, which leads to a loss of information about past states over time, even in an idealized setting. Overall, our experiments demonstrate AGILE’s potential as a flexible and efficient data assimilation framework. Its ability to integrate diverse datasets in a dynamically consistent manner makes it a promising tool for future real-world glacier modeling applications.

Global hydrological models have become a valuable tool for a range of global impact studies related to water resources. However, glacier parameterization is often simplistic or non-existent in global hydrological models. By contrast, global glacier models do represent complex glacier dynamics and glacier evolution, and as such hold the promise of better resolving glacier runoff estimates. In this study, we test the hypothesis that coupling a global glacier model with a global hydrological model leads to a more realistic glacier representation and consequently an improved runoff prediction in the global hydrological model. To this end, the Global Glacier Evolution Model (GloGEM) is coupled with the global hydrological model PCR-GLOBWB 2 using the eWaterCycle platform. For the period 2001–2012, the coupled model is evaluated against the uncoupled PCR-GLOBWB 2 in 25 large-scale (>50.000 km2) glacierized basins. The coupled model produces higher runoff estimates across all basins and throughout the melt season. In summer, the runoff differences range from 0.07 % for weakly glacier-influenced basins to 252 % for strongly glacier-influenced basins. The difference can primarily be explained by PCR- GLOBWB 2 not accounting for glacier flow and glacier mass loss, thereby causing an underestimation of glacier runoff. The coupled model performs better in reproducing basin runoff observations mostly in strongly glacier-influenced basins, which is where the coupling has the most impact. This study underlines the importance of glacier representation in global hydrological models and demonstrates the potential of coupling a global hydrological model with a global glacier model for better glacier representation and runoff predictions in glacierized basins.

Model and data reduction for data assimilation: Particle filters employing projected forecasts and data with application to a shallow water model

Due to climate warming, the glaciers of the Tibetan Plateau have experienced rapid mass loss and the patterns of glacier changes have exhibited high spatiotemporal heterogeneity, especially in interior areas. As the largest ice field within the Tibetan Plateau, the Puruogangri Ice Field has attracted a lot of attention from the scientific community. However, relevant studies that are based on satellite data have mainly focused on a few periods from 2000–2016. Long-term and multiperiod observations remain to be conducted. To this end, we estimated the changes in the glacier area and mass balance of the Puruogangri Ice Field over five subperiods between 1975 and 2021, based on multisource remote sensing data. Specifically, we employed KH-9 and Landsat images to estimate the area change from 1975 to 2021 using the band ratio method. Subsequently, based on KH-9 DEM, SRTM DEM, Copernicus DEM and ZY-3 DEM data, we evaluated the glacier elevation changes and mass balance over five subperiods during 1975–2021. The results showed that the total glacier area decreased from 427.44 ± 12.43 km2 to 387.87 ± 11.02 km2 from 1975 to 2021, with a decrease rate of 0.86 km2 a−1. The rate of area change at a decade timescale was −0.74 km2 a−1 (2000–2012) and −1.00 km2 a−1 (2012–2021). Furthermore, the rates at a multiyear timescale were −1.23 km2 a−1, −1.83 km2 a−1 and −0.42 km2 a−1 for 2012–2015, 2015–2017 and 2017–2021, respectively. In terms of the glacier mass balance, the region-wide results at a two-decade timescale were −0.23 ± 0.02 m w.e. a−1 for 1975–2000 and −0.29 ± 0.02 m w.e. a−1 for 2000–2021, indicating a sustained and relatively stable mass loss over the past nearly five decades. After 2000, the loss rate at a decade timescale was −0.04 ± 0.01 m w.e. a−1 for 2000–2012 and −0.17 ± 0.01 m w.e. a−1 for 2012–2021, indicating an increasing loss rate over recent decades. It was further found that the mass loss rate was −0.12 ± 0.02 m w.e. a−1 for 2012–2015, −0.03 ± 0.01 m w.e. a−1 for 2015–2017 and −0.40 ± 0.03 m w.e. a−1 for 2017–2021. These results indicated that a significant portion of the glacier mass loss mainly occurred after 2017. According to our analysis of the meteorological measurements in nearby regions, the trends of precipitation and the average annual air temperature both increased. Combining these findings with the results of the glacier changes implied that the glacier changes seemed to be more sensitive to temperature increase in this region. Overall, our results improved our understanding of the status of glacier changes and their reaction to climate change in the central Tibetan Plateau.

Data assimilation techniques, developed in the past two decades mainly for weather prediction, produce better forecasts by taking advantage of both theoretical/numerical models and real-time observations. In this paper, we explore the possibility of applying the four-dimensional variational data assimilation (4D-VAR) technique to the prediction of solar flares. We do so in the context of a continuous version of the classical cellular-automaton-based self-organized critical avalanche models of solar flares introduced by Lu and Hamilton (Astrophys. J. 380, L89, 1991). Such models, although a priori far removed from the physics of magnetic reconnection and magnetohydrodynamical evolution of coronal structures, nonetheless reproduce quite well the observed statistical distribution of flare characteristics. We report here on a large set of data assimilation runs on synthetic energy release time series. Our results indicate that, despite the unpredictable (and unobservable) stochastic nature of the driving/triggering mechanism within the avalanche model, 4D-VAR succeeds in producing optimal initial conditions that reproduce adequately the time series of energy released by avalanches and flares. This is an essential first step toward forecasting real flares.

Studies of the global sea-level budget (SLB) and the global ocean-mass budget (OMB) are essential to assess the reliability of our knowledge of sea-level change and its contributions. Here we present datasets for times series of the SLB and OMB elements developed in the framework of ESA's Climate Change Initiative. We use these datasets to assess the SLB and the OMB simultaneously, utilising a consistent framework of uncertainty characterisation. The time series, given at monthly sampling, include global mean sea-level (GMSL) anomalies from satellite altimetry; the global mean steric component from Argo drifter data with incorporation of sea surface temperature data; the ocean mass component from Gravity Recovery and Climate Experiment (GRACE) satellite gravimetry; the contribution from global glacier mass changes assessed by a global glacier model; the contribution from Greenland Ice Sheet and Antarctic Ice Sheet mass changes, assessed from satellite radar altimetry and from GRACE; and the contribution from land water storage anomalies assessed by the WaterGAP global hydrological model. Over the period Jan 1993–Dec 2016 (P1, covered by the satellite altimetry records), the mean rate (linear trend) of GMSL is 3.05 ± 0.24 mm yr−1. The steric component is 1.15 ± 0.12 mm yr−1 (38 % of the GMSL trend) and the mass component is 1.75 ± 0.12 mm yr−1 (57 %). The mass component includes 0.64 ± 0.03 mm yr−1 (21 % of the GMSL trend) from glaciers outside Greenland and Antarctica, 0.60 ± 0.04 mm yr−1 (20 %) from Greenland, 0.19 ± 0.04 mm yr−1 (6 %) from Antarctica, and 0.32 ± 0.10 mm yr−1 (10 %) from changes of land water storage. In the period Jan 2003–Aug 2016 (P2, covered by GRACE and the Argo drifter system), GMSL rise is higher than in P1 at 3.64 ± 0.26 mm yr−1. This is due to an increase of the mass contributions (now about 2.22 ± 0.15 mm yr−1, 61 % of the GMSL trend), with the largest increase contributed from Greenland. The SLB of linear trends is closed for P1 and P2, that is, the GMSL trend agrees with the sum of the steric and mass components within their combined uncertainties. The OMB budget, which can be evaluated only for P2, is also closed, that is, the GRACE-based ocean-mass trend agrees with the sum of assessed mass contributions within uncertainties. Combined uncertainties (1-sigma) of the elements involved in the budgets are between 0.26 and 0.40 mm yr−1, about 10 % of GMSL rise. Interannual variations that overlie the long-term trends are coherently represented by the elements of the SLB and the OMB. Even at the level of monthly anomalies the budgets are closed within uncertainties, while also indicating possible origins of remaining misclosures.

<p>Studies of the global sea-level budget (SLB) and ocean-mass budget (OMB) are essential to assess the reliability of our knowledge of sea-level change and its contributors. The SLB is considered closed if the observed sea-level change agrees with the sum of independently assessed steric and mass contributions. The OMB is considered closed if the observed ocean-mass change is compatible with the sum of assessed mass contributions. </p><p>Here we present results from the Sea-Level Budget Closure (SLBC_cci) project conducted in the framework of ESA’s Climate Change Initiative (CCI). We used data products from CCI projects as well as newly-developed products based on CCI products and on additional data sources. Our focus on products developed in the same framework allowed us to exercise a consistent uncertainty characterisation and its propagation to the budget closure analyses, where the SLB and the OMB are assessed simultaneously. </p><p>We present time series of global mean sea-level changes from satellite altimetry; new time series of the global mean steric component generated from Argo drifter data with incorporation of sea surface temperature data; time series of ocean-mass change derived from GRACE satellite gravimetry; time series of global glacier mass change from a global glacier model; time series of mass changes of the Greenland Ice Sheet and the Antarctic Ice Sheet both from satellite radar altimetry and from GRACE; as well as time series of land water storage change from the WaterGAP global hydrological model. Our budget analyses address the periods 1993–2016 (covered by the satellite altimetry records) and 2003–2016 (covered by GRACE and the Argo drifter system). In terms of the mean rates of change (linear trends), the SLB is closed within uncertainties for both periods, and the OMB, assessable for 2003–2016 only, is also closed within uncertainties. Uncertainties (1-sigma) arising from the combined uncertainties of the elements of the different budgets considered are between 0.26 mm/yr and 0.40 mm/yr, that is, on the order of 10% of the magnitude of global mean sea-level rise, which is 3.05 ± 0.24 mm/yr and 3.65 ± 0.26 mm/yr for 1993-2016 and 2003-2016, respectively. We also assessed the budgets on a monthly time series basis. The statistics of monthly misclosure agrees with the combined uncertainties of the budget elements, which amount to typically 2-3 mm for the 2003–2016 period. We discuss possible origins of the residual misclosure.</p>

Ensemble simulations have been performed with a climate model constrained to follow temperature histories obtained from a recent compilation of 56 well‐calibrated surface temperature proxy records, using a new data assimilation technique. First, we demonstrate that the data assimilation technique provides a faithful representation in the Northern Hemisphere of the signal recorded by the climate proxies at both the regional and gridbox scales. Second, by varying the external forcing, the parameters of the data assimilation method, and the parameters controlling the equilibrium climate sensitivity of the climate model, we demonstrate that the uncertainty in model results is much lower in simulations using data assimilation than in those without it. This observation implies that the data assimilation, using a set of 56 proxies, is providing an efficient and robust constraint on the simulated climate variability over the past centuries. At the hemispheric and continental scales, the model reconstructions using data assimilation are in good agreement with both the instrumental record of the past 150 years and reconstructions of climate in past centuries derived from the application of traditional statistical approaches to networks of proxy data. This increases the confidence in both the data assimilation and traditional statistical approaches. Our data assimilation method, however, is unable to provide a reliable reconstruction over the North Atlantic Ocean, which we attribute to the paucity of proxy data in that region.

Challenges for data assimilation – from convective‐scale to climate

<p>Ongoing global glacier retreat leads to sea-level rise and changes in regional freshwater availability. For an adequate adaptation to these changes, knowledge about the ice volume and the current dynamic state of glaciers is crucial. At regional to global scales, sparse observations made the dynamic state of glaciers very difficult to assess. Thanks to recent advances in global geodetic mass-balance and velocity assessments, new ways to initialize numerical models and ice thickness estimation emerge. In this contribution, we present the COst Minimization Bed INvErsion model (COMBINE), which aims to be a cheap, flexible global data assimilation and inversion method. COMBINE uses an existing numerical model of glacier evolution (the Open Global Glacier Model, OGGM) rewritten in the machine learning framework PyTorch. This makes the model fully differentiable and allows to iteratively minimize a cost function penalizing mismatch to observations. Thanks to the flexible nature of automatic differentiation, various observational sources distributed in time can be considered (e.g. surface elevation and area changes, ice velocities). No assumption about the dynamic glacier state is needed, releasing the equilibrium assumption often required for large scale ice volume computations. In this contribution, we will demonstrate the capabilities of COMBINE in several idealized and real-world applications, and discuss its added value and upcoming challenges for operational application.</p>

Slip at the ice-bed interface (basal motion) dominates the flow of many glaciers, and it is uncertain whether this velocity component will increase or slow in a warmer world. Past results from an idealized flowline glacier model show that declining basal motion induces a two-phase response that initially accelerates glacier retreat in a warming climate on a multidecadal timescale but lessens centennial-scale retreat and mass loss. In the present work, we utilize existing field-collected and remotely-sensed constraints on ice thickness, ice surface velocity, and the change in each of these terms to constrain the current rate of basal motion and its change over the past ~40 years. We focus on the ~1500 global glaciers with higher density of field-based ice thickness measurements in the GlaThiDa dataset (>18 measurements points per glacier). Utilizing these ice thickness and surface velocity constraints, we employ a flow model to estimate the rate of basal motion as the residual between observed surface velocity and modeled ice deformation. We first estimate the contribution of varying basal motion to observed changes in surface velocity across the study glaciers.  We then estimate these glaciers’ retreat and thinning responses to changing velocity and compare these with the magnitudes expected from atmospheric warming, constrained by published point measurements, mass balance models, and snowline observations. These results will constrain the extent to which evolving ice dynamics have amplified or mitigated the response of global glaciers to climate change over past decades. Further, this knowledge will provide insight into the potential importance of varying basal motion on projections of future glacier change, with implications for global sea level rise as well as local water resource and ecosystem management.

A better carbon-water flux simulation in multiple vegetation types by data assimilation

The water balance of high-alpine glacierized catchments is largely dominated by snow and ice processes. When modelling the hydrological response of such catchments, a reliable representation of snow/ice accumulation and melt should be ensured, especially when studying the effects of climate change. Even though numerous state-of-the-art hydrological models are able to adequately represent the contribution of snow melt into the total runoff with the use of complex approaches (e.g. energy balance models), glacier dynamics are still based on conceptual or empirical methods, which exhibit some limitations compared to more sophisticated models (e.g. explicit ice-flow dynamics).The Water Flow and Balance Simulation Model (WaSiM) is a process-based hydrological model that includes an empirical volume-area scaling approach for describing the glacier’s evolution. Although acceptable estimates can be obtained with this approach, an integration to a more complex glacier representation is still missing. For this reason, a coupling scheme between WaSiM and the Open Global Glacier Model (OGGM) is developed, hence accounting for explicit ice-flow dynamics.The workflow consists mainly on three steps: i) a first WaSiM run to obtain monthly values of temperature and precipitation that serve as input for the ii) second step, which is running OGGM. Finally, iii) a dynamic model run of WaSiM with the updated output from OGGM (annual glacier outlines and ice thickness) is performed. Within this last step, the glacier’s volume internally calculated by WaSiM (i.e. with the VA-scaling approach) is replaced by OGGM’s output, while performing a simultaneous multi-data set automatic calibration. In this calibration, only WaSiM parameters are adjusted and simulation results are compared against glacier mass balances (OGGM) and observed runoff. The performance of the calibration is then evaluated in terms of a weighted multi-objective function. Although the best fit between observed and simulated runoff is achieved when considering only runoff observations (single-data calibration), glacier components are better represented when calibrating the coupled model with the multi-data set (i.e. also including glacier mass balances). Therefore, a trade-off is made between general model performance and accurate runoff prediction. This coupling scheme is aimed for hydrological modellers with no additional expertise on glacier modelling, since OGGM is set up according to its default parameters. Finally, it could serve as a tool not only to predict the hydrological response of any glacierized catchment (even without any available glacier data), but also to make predictions under future climate projections with a more reliable representation of glaciers.