2.1 Bayesian inference

Data assimilation, loosely the fusion of data and models, can be formalized as the application of Bayesian inference (Wikle and Berliner, 2007). As such, we begin the section by briefly reviewing the key ideas behind Bayesian inference that is at the core of most modern DA schemes, including those explored herein. We refer the reader to the comprehensive texts of Gelman et al. (2013), Särkkä and Svensson (2023), Murphy (2023) and Evensen et al. (2022) for a thorough Bayesian treatment from the perspectives of statistics, applied mathematics, machine learning and geophysical DA, respectively. A more thorough analysis of the connection between formal Bayesian inference and practical cryospheric DA is provided in Alonso-González et al. (2022).

In essence, Bayesian inference can be viewed as updating beliefs about some uncertain (also known as random) variables θ given some data $\mathit{y}\in {\mathbb{R}}^{N_{\mathrm{o}}}$. The uncertain variables in the vector θ can generally be of any form. Here, without loss of generality, we will restrict our attention to continuous model parameters $\mathit{\theta }\in {\mathbb{R}}^{N_{\mathrm{p}}}$. Belief updating can then be achieved by combining the basic sum and product rules of probability (Jaynes, 2003; MacKay, 2003) to obtain Bayes' rule

which states that the posterior belief after conditioning on the data, p(θ∣y), is proportional to the product of the likelihood, p(y∣θ), which is loosely speaking what the data tell us, and the prior, p(θ), which encodes what we believed about θ before considering the data. For the model evidence term p(y) in the denominator of Eq. (1) we again combine the product and sum rules to see that

where here and throughout this study the limits of integration are implicitly over the entire support of the integrand. From Eq. (2) the evidence is independent of the parameters θ and simply the integral of the numerator of Eq. (1) which ensures that the posterior integrates to one over its support. Thus, the evidence can be seen as just a normalizing constant, although it plays an important role as a marginal likelihood for higher levels of inference, since it is implicitly conditioned on the model (MacKay, 2003; Murphy, 2023). In this Bayesian framework, the probability densities p(⋅) encode beliefs (i.e. epistemic uncertainties) about their arguments. If the subjective nature of beliefs is unpalatable, it may help to imagine that these are the beliefs held by an abstract numerical agent with a probabilistic model of the world (Hennig et al., 2022). This probabilistic numerics perspective is instructive for the adaptive methods that we will present here.

In theory, by inspection of Eq. (1), performing Bayesian inference is simply (up to a normalizing constant) multiplying the prior and the likelihood. Naively then, we could just enumerate the posterior through an exhaustive grid approximation provided that we select a sufficiently fine discretization of the variables θ (MacKay, 2003). However, as soon as we have more than a couple of uncertain variables in θ and/or very informative data in y, this grid approach becomes intractable due to the so-called curse of dimensionality (Murphy, 2023). Recalling and extending the earlier analogy, the task of inferring a possibly multi-modal posterior distribution p(θ∣y) is akin to looking for an unknown number of needles in a multi-dimensional haystack: the larger the support (e.g. dimensionality) of the prior, the larger the haystack, and the more informative (numerous and/or accurate) the observations, the smaller the needles. It is this computational challenge that has spurred the development of sophisticated Bayesian inference algorithms, including ensemble-based (also known as Monte Carlo) algorithms that emerged from statistical physics (MacKay, 2003; Robert and Casella, 2004) and have, along with variational methods, enjoyed widespread adoption in geophysical DA (Reich and Cotter, 2015; Evensen et al., 2022). In this study, we restrict our attention to ensemble-based methods since these gradient-free methods are the most widely used in cryospheric DA due to their relative ease of use and comparatively robust uncertainty quantification (see Alonso-González et al., 2022, and references therein). At the same time, more research is warranted to explore the ever-growing plethora of inference algorithms (Evensen et al., 2022; Hennig et al., 2022; Murphy, 2023; Särkkä and Svensson, 2023), many of which remain relatively untested in cryospheric science and geoscience more generally.

2.2 Cryospheric inverse problems

To help concretize the above formalization of cryospheric DA, it is instructive to consider the perspective of inverse modeling (Sanz-Alonso et al., 2023). Let 𝒢(⋅) denote our forward (or data generating) model that maps from the uncertain snow model parameters to predicted (modeled) snow observations $\widehat{\mathit{y}}\in {\mathbb{R}}^{N_{\mathrm{o}}}$, i.e., $\widehat{\mathit{y}}=\mathcal{G}\left(\mathit{\theta }\right)$. In practice, the forward model combines an observation model for $\widehat{\mathit{y}}=\mathcal{H}\left(\mathit{x}\right)$ with a dynamical model for the full model state trajectory x=ℳ(ϕ) and a (bounding) parameter transformation step ϕ=𝒯(θ) (Gelman et al., 2013; Alonso-González et al., 2022; Pirk et al., 2022). Next, consider the typical case where we are given a set of noisy snow observations y which we assume to be related to some true snow observable y^⋆ through $\mathit{y}={\mathit{y}}^{\star }+\mathit{ϵ}$, where ϵ is the observation error. By making the usual strong constraint (or perfect model) assumption (Evensen et al., 2022) that the forward model maps perfectly (without error) from the parameters to the observed variables, we can define some true parameter set θ^⋆ to exist such that ${\mathit{y}}^{\star }=\mathcal{G}\left({\mathit{\theta }}^{\star }\right)$. In general, these parameters can include internal model parameters, initial conditions, and/or forcing terms. This corresponds to a widely used strong constraint version of the forcing formulation of cryospheric DA, wherein model states are completely determined by this parameter vector and the forward model (Alonso-González et al., 2022). We adopt this approach throughout without loss of generality, as the forcing formulation can accommodate a weak constraint with model error (Evensen et al., 2022).

Using the definition of the observation error, we can establish the following forward relationship between the noisy observations we are given and some presumed true parameters of interest

At least conceptually, the task in cryospheric data assimilation can now be cast as somehow inverting 𝒢(⋅) in Eq. (3) to solve for θ^⋆. Unfortunately, even in an idealized linear and noise-free (ϵ=0) case, this is typically an ill-posed problem in that an exact solution θ^⋆ may not exist or be unique (Sanz-Alonso et al., 2023). In the more challenging noisy and possibly non-linear practical settings the problem is always ill-posed since we invariably assimilate noisy observations where the observation error is uncertain.

Due to ill-posedness, the quest for a universally optimal (let alone exact) solution θ^⋆ is nonsensical since infinitely many solutions θ are typically admissible. A way forwards is to instead seek a probabilistic (i.e. Bayesian) solution to the inverse problem in Eq. (3) where we construct an observation error model by treating ϵ as an uncertain variable. For convenience, as is common practice in DA (Carrassi et al., 2018), we assume independent additive zero-mean Gaussian observation errors $\mathit{ϵ}\sim N(\mathbf{0},\mathbf{R})$ where R is an N_o×N_o diagonal observation error covariance matrix. This assumption helps formulate the likelihood p(y∣θ), the probability of obtaining the fixed observations y given that a parameter set θ is true, since if $\mathit{\theta }={\mathit{\theta }}^{\star }$ (due to the conditional ∣θ) then from Eq. (3) $\mathit{ϵ}=\mathit{y}-\mathcal{G}\left(\mathit{\theta }\right)$ is the residual so we apply the assumed observation error model to obtain a Gaussian likelihood (Evensen et al., 2022)

where $c_y=\det (\mathrm{2}\mathit{\pi }\mathbf{R}{)}^{-\mathrm{1}/\mathrm{2}}$ and $\widehat{\mathit{y}}=\mathcal{G}\left(\mathit{\theta }\right)$ denotes the predicted (i.e. modeled) observations given a particular parameter set θ whereby diagnosing the likelihood in Eq. (4) requires point-wise evaluations of the forward model to evaluate the residual. In accordance with the likelihood principle, p(y∣θ) should be viewed as a function of the uncertain parameters θ rather than the fixed (albeit noisy) observations y that we are assimilating (MacKay, 2003). If we now combine this likelihood with a regularizing prior p(θ) that encodes initial beliefs concerning the parameters, then in principle the full probabilistic solution to the inverse problem Eq. (3) is obtained by inferring the posterior through Bayes' rule in Eq. (1). As noted by Evensen et al. (2022) this process of Bayesian inference is just point-wise multiplication that does not in itself involve any explicit (matrix or function) inversion, but it nonetheless offers a probabilistic framework for solving general geophysical inverse problems (Sanz-Alonso et al., 2023). Although distinctions are sometimes made between DA and inverse modeling, the two fields are highly complementary and unified under the umbrella of Bayesian inference (Reich and Cotter, 2015; Evensen et al., 2022; Sanz-Alonso et al., 2023). Thereby, the process of cryospheric DA can be viewed as the solution to dynamical (time-varying) cryospheric inverse problems. The Bayesian dynamics then determines if one is solving a filtering or more general smoothing problem (Alonso-González et al., 2022; Särkkä and Svensson, 2023), either way the solutions are typically obtained via numerical approximation of Bayesian inference which is our focus.

2.3 Markov chain Monte Carlo

Markov Chain Monte Carlo (MCMC) methods, which are asymptotically exact, are widely considered the reference computational tool for approximating Bayesian posteriors in practice (Gelman et al., 2013; Murphy, 2023). Here we provide a brief overview of MCMC since we will use it as a gold-standard computational benchmark against which to gauge the performance of other cryospheric DA methods (Law and Stuart, 2012). As expounded in Robert and Casella (2004), MCMC originated with the seminal physics paper of Metropolis et al. (1953), whose results were later generalized to statistics by Hastings (1970), introducing the Random Walk Metropolis (RWM) algorithm that is the ancestor of modern MCMC methods, among which gradient-based methods such as Hamiltonian Monte Carlo (MacKay, 2003; Neal, 2011) are arguably the state-of-the-art. All MCMC schemes construct Markov chains so as to asymptotically sample from a target distribution of interest, where the natural choice for Bayesian DA is the posterior p(θ∣y). More specifically, MCMC algorithms take sequential Markovian (memoryless) steps through the parameter space and probabilistically either accept or reject the newly proposed step by comparing its posterior density to that of the current step, eventually obtaining samples from the posterior distribution in Eq. (1).

We base our implementation here on the aforementioned RWM algorithm which remains arguably the most archetypal MCMC method. In this approach, a new step ${\mathit{\theta }}_{i+\mathrm{1}}^{*}$ in the chain is proposed by randomly drawing from a proposal distribution $q({\mathit{\theta }}_{i+\mathrm{1}}^{*}\mid {\mathit{\theta }}_i)$ that is conditioned (usually by centering) on the current step θ_i. This proposed step is probabilistically accepted based on the value of the acceptance ratio

on the condition that the acceptance rule $u_{i+\mathrm{1}}\le \min ({\mathit{\alpha }}_{i+\mathrm{1}},\mathrm{1})$ holds, where $u_{i+\mathrm{1}}\sim U(\mathrm{0},\mathrm{1})$ is a realization of a random variable that is uniformly distributed between 0 and 1, otherwise it is rejected. The evidence p(y), which appears as a constant in the posterior density Eq. (1) for both the current and the proposed step, cancels out in the acceptance ratio, so we do not need to estimate this intractable quantity. This, together with asymptotic guarantees, helps explain the historical popularity of MCMC sampling. The more general form of Eq. (5) introduced by Hastings (1970) also involves proposal densities, but these cancel out in RWM and were hence omitted here. Upon acceptance, the chain moves to the proposed point in parameter space such that ${\mathit{\theta }}_{i+\mathrm{1}}={\mathit{\theta }}_{i+\mathrm{1}}^{*}$. Upon rejection, the chain stays at the current point ${\mathit{\theta }}_{i+\mathrm{1}}={\mathit{\theta }}_i$. Note that a proposed point will always be accepted if it has a higher posterior density since then ${\mathit{\alpha }}_{i+\mathrm{1}}>\mathrm{1}$ so the acceptance rule always holds. Moreover, the proposed point ${\mathit{\theta }}_{i+\mathrm{1}}^{*}$ will always have a non-zero probability of being accepted. This probabilistic acceptance rule helps ensure that the chain will asymptotically (i→∞) sample from the posterior. To sample from the posterior with MCMC, one needs to run the Markov chain for enough iterations to ensure that it mixes properly and converges to sampling from the posterior. In practice, although there are some diagnostics that can be used as a guide (Gelman et al., 2013), what constitutes enough iterations is uncertain a priori and can remain challenging to verify post hoc. As such, MCMC is typically run for a very large number (i.e., tens of thousands) of iterations to ensure convergence (Cleary et al., 2021). An initial part of the chain is discarded as a burn-in period to avoid the biasing effects of the chain initialization. Due to the Markov property, where the next step only depends on the current step, it is also clear that the subsequent iterations in the chain will be auto-correlated rather than independent. As such, the remaining part of the chain can be subsampled uniformly at random to obtain more independent samples from the posterior distribution.

The most standard RWM approach employs a multivariate normal (Gaussian) proposal of the form $q({\mathit{\varphi }}_{j+\mathrm{1}}^{*}\mid {\mathit{\varphi }}_j)=N({\mathit{\varphi }}_{j+\mathrm{1}}^{*}\mid {\mathit{\varphi }}_j,{\mathbf{\Sigma }}_q)$ where the mean is the current point ϕ_j and Σ_q is the proposal covariance matrix. To reduce the number of tuning parameters in the proposal, the latter can be made isotropic by using a scalar (constant diagonal) matrix of the form ${\mathbf{\Sigma }}_q={\mathit{\sigma }}_q^{\mathrm{2}}\mathbf{I}$ where ${\mathit{\sigma }}_q^{\mathrm{2}}$ is the proposal variance and I is an identity matrix. In practice, good mixing of the RWM method requires judicious hand-tuning of σ_q, and even then the use of an isotropic covariance matrix remains wasteful to efficiently explore higher dimensional parameter spaces given that the posterior is often anisotropic (MacKay, 2003). A relatively simple workaround is to employ an adaptive RWM method that automatically modifies a generally anisotropic proposal covariance on the fly using the history of the evolving Markov chain to obtain better mixing properties. Among the existing adaptive MCMC algorithms, here we chose to employ the robust adaptive Metropolis (RAM) method proposed by Vihola (2012) using the hyperparameters suggested therein with $N_{\mathrm{s}}=\mathrm{20}\times {\mathrm{10}}^{\mathrm{3}}$ steps and discarding the first 10  % as a burn-in phase (Murphy, 2023). The RAM method was chosen in-lieu of even more sophisticated MCMC methods such as Hamiltonian Monte Carlo (Neal, 2011) that may have mixed even faster since the RAM method is considerably easier to implement. Moreover, by running RAM for tens of thousands of iterations we can be fairly confident that it converges and represents a strong benchmark. A crucial point here is that it is likely that most if not all MCMC methods are too computationally expensive to deploy at scale (i.e., for large areas) in practical cryospheric DA, although further research into applying more sophisticated MCMC schemes (Murphy, 2023) is needed. Herein, MCMC in the form om RAM is used as a reference standard against which we gauge the performance of the more tractable ensemble-based cryospheric DA schemes (Law and Stuart, 2012), especially the adaptive particle smoother that is the focus of this study.

2.4 Ensemble Kalman methods

Ensemble Kalman methods, introduced by Evensen (1994) and described in detail in Evensen et al. (2022), are ensemble-based extensions of the classic Kalman methods (Jazwinski, 1970) which provide exact Bayesian inference methods for linear Gaussian models (Särkkä and Svensson, 2023). In such models the mapping from hidden states and/or parameters to observations is linear while both the prior and the likelihood are Gaussian. Ensemble Kalman methods help to relax these assumptions in the sense that approximate yet efficient inference is still possible when they are violated, making them applicable also to non-linear geophysical problems (Evensen et al., 2022). Although other non-linear variations on classic Kalman methods also exist (Särkkä and Svensson, 2023), the ease of implementation and the robust performance of ensemble Kalman methods have made them highly applicable for geophysical DA in general (Carrassi et al., 2018) as well as cryospheric DA in particular (see Alonso-González et al., 2022, and references therein). Crucially, the more recent development of iterative ensemble Kalman methods (Emerick and Reynolds, 2013; Garbuno-Inigo et al., 2020) have helped to enhance the capability of this class of DA schemes for highly non-linear and/or complex problems (e.g. Pirk et al., 2022, 2024; Keetz et al., 2025).

In practical cryospheric DA, several numerical experiments have previously shown that iterative ensemble Kalman methods can outperform basic particle methods while maintaining more robust posterior uncertainty quantification (Aalstad et al., 2018; Alonso-González et al., 2022). It is thus instructive to also include such experiments here to provide an additional benchmark for the performance of the proposed adaptive particle method. Rather than providing a high-fidelity benchmark like the MCMC experiments, the more computationally tractable iterative ensemble Kalman experiments should be seen as a more practical operational benchmark for cryospheric DA. For the sake of completeness, the rest of this section provides a brief overview on the implementation of iterative ensemble Kalman methods. We refer to Evensen et al. (2022) and Alonso-González et al. (2022) for more details on the theory of ensemble Kalman methods and their implementation for cryospheric DA, respectively.

Here we adopt a specific iterative ensemble Kalman method known as the ensemble smoother with multiple data assimilation (ES-MDA; Emerick and Reynolds, 2013) due to its relative ease of implementation and robust performance (Aalstad et al., 2018; Alonso-González et al., 2022, 2023). We note in passing that other promising variations on the iterative ensemble Kalman method exist (Garbuno-Inigo et al., 2020; Evensen et al., 2022) and are worthy of further investigation in cryospheric DA, but we do not expect their performance to differ markedly from the ES-MDA. The ES-MDA scheme is initialized by sampling an initial ensemble of $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{e}}$ parameter vectors from the prior ${\mathit{\theta }}_i^{\left(\mathrm{0}\right)}\sim p\left(\mathit{\theta }\right)$ and it then proceeds by cycling between a prediction and update step for $\mathrm{\ell }=\mathrm{0},\mathrm{\dots },N_{\mathrm{a}}$ iterations:

1. Run the forward model to obtain the predicted observations ${\widehat{\mathit{y}}}_i^{\left(\mathrm{\ell }\right)}=\mathcal{G}\left({\mathit{\theta }}_i^{\left(\mathrm{\ell }\right)}\right)$ from the state x (see Sect. 2.2).

2. If ℓ<N_a, perform a tempered ensemble Kalman update step

   $$\begin{array}{}\text{(6)}& {\mathit{\theta }}_i^{(\mathrm{\ell }+\mathrm{1})}={\mathit{\theta }}_i^{\left(\mathrm{\ell }\right)}+{\mathbf{K}}^{\left(\mathrm{\ell }\right)}\left(\mathit{y}-{\widehat{\mathit{y}}}_i^{\left(\mathrm{\ell }\right)}\right),\end{array}$$

   where K^(ℓ) is the tempered ensemble Kalman gain for iteration ℓ obtained from (ensemble) covariance matrices (Evensen et al., 2022).

The steps above are implicitly carried out for all $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{e}}$ ensemble members, with N_a+1 iterations of the prediction step and N_a iterations of the update step. The update step itself is essentially free, so the total cost of the ES-MDA is $(N_{\mathrm{a}}+\mathrm{1})\times N_{\mathrm{e}}$ forward model simulations where it is possible to parallelize across the ensemble dimension in each iteration ℓ. Recall from Sect. 2.2 that several steps (observation, dynamics, transformation) are baked into the forward model 𝒢(⋅). Therein, it is the cost of running the dynamical model ℳ(⋅) to obtain the ensemble of hidden model states x of interest that completely dominates the computational burden of 𝒢(⋅) in step 1. So, while (single chain) MCMC costs tens of thousands of strictly sequential forward model runs, with a typical setting of N_a=4 (Aalstad et al., 2018; Pirk et al., 2022; Alonso-González et al., 2022) ES-MDA only incurs a computational cost of $N_{\mathrm{a}}+\mathrm{1}=\mathrm{5}$ iterations of an ensemble of N_e=100 parallelizable forward model runs. With N_a=1 the ES-MDA reverts to the original (non-iterative) ensemble smoother (ES) scheme proposed by van Leeuwen and Evensen (1996) which we also include here in the benchmarking of the new AdaPBS method. We emphasize that even this non-iterative ES with N_a=1 has a cost of $(N_{\mathrm{a}}+\mathrm{1})\times N_{\mathrm{e}}=\mathrm{200}$ model runs as it requires rerunning an ensemble of model simulations with the updated parameters to obtain posterior state predictions.

2.5 Particle methods

Particle methods (van Leeuwen, 2009; van Leeuwen et al., 2019), also known as sequential Monte Carlo (SMC) (Chopin and Papaspiliopoulos, 2020), rose to prominence with the work of Gordon et al. (1993) and Kitagawa (1996) around the same time as ensemble Kalman methods (Evensen, 1994). Moreover, similar inference methods have arguably independently been (re)discovered in geoscience (Beven and Binley, 1992; van Leeuwen and Evensen, 1996; Margulis et al., 2015). These particle methods also have roots back to the dawn of Monte Carlo methods in physics (Hammersley and Morton, 1954; Robert and Casella, 2004). The key inference mechanism that powers particle methods is importance sampling (MacKay, 2003), which weighs the importance of an ensemble of particles (ensemble members) θ_i according to their posterior probability density and the density of the proposal. By performing this sequentially in time and resampling the particles based on their weights (i.e., fitness) after each importance sampling step, we recover the algorithm known as sequential importance resampling (SIR; Smith and Gelfand, 1992; Doucet et al., 2000) at the core of particle methods used for Bayesian filtering and smoothing (Kitagawa, 1996) and static inference problems (Chopin, 2002). The cycling of prediction (mutation) and updating (selection) has clear connections with metaheuristic evolutionary methods (Campelo and Aranha, 2023) such as genetic algorithms (Holland, 1992) that can be mathematically formalized as particle methods (Del Moral, 2004). This perspective also provides a connection to other evolutionary-inspired Bayesian such as Shuffled Complex Evolution Metropolis (SCEM-UA; Vrugt et al., 2003) and Differential Evolution Adaptive Metropolis (DREAM; Vrugt et al., 2008), that remain state of the art optimizers and samplers, respectively, in hydrology. However, following Sörensen (2015) we emphasize that the adaptive particle method proposed herein and our playful title use evolutionary principles as a conceptual model for inspiration and not as a justification. Justification can be found in the literature on the convergence of general SMC methods (Chopin and Papaspiliopoulos, 2020) and more specifically adaptive multiple importance sampling (Marin et al., 2019). Having introduced the notion of particle methods, the rest of this section will outline the principle of importance sampling and how this is incorporated into SIR. Together, this basic SIR theory suffices to grasp vanilla particle methods related to the seminal “bootstrap” particle filter of Gordon et al. (1993). These basic particle methods have become popular approaches to DA in snow science (e.g. Leisenring and Moradkhani, 2011; Margulis et al., 2015; Charrois et al., 2016; Magnusson et al., 2017; Piazzi et al., 2018; Fiddes et al., 2019; Smyth et al., 2019; Alonso-González et al., 2021; Cluzet et al., 2021; Oberrauch et al., 2024; Girotto et al., 2024; Sun et al., 2025) and glaciology (Navari et al., 2016; Landmann et al., 2021; Cao et al., 2025).

2.6 Adaptive particle methods

The typical particle methods used in cryospheric DA (see Alonso-González et al., 2022, and references therein) can be considered to be basic particle methods in the spirit of the seminal bootstrap particle filter (Gordon et al., 1993) in that they do not leverage recent developments (van Leeuwen et al., 2019; Chopin and Papaspiliopoulos, 2020). They are basic in the sense that they use the simplest and most convenient choice of proposal distribution (van Leeuwen, 2009), namely using the prior as the proposal. In this case, i.e. inserting q(θ_i)=p(θ_i) in Eq. (12), the unnormalized weights simply correspond to likelihood evaluations ${\stackrel{\mathrm{̃}}{w}}_i=p(\mathit{y}\mid {\mathit{\theta }}_i)$. Thereby, the normalized weights in the weighted mean approximation of the posterior expectation in Eq. (14) just become normalized likelihoods of the form $w_i=p(\mathit{y}\mid {\mathit{\theta }}_i)/{\sum }_{k=\mathrm{1}}^{N_{\mathrm{e}}}p(\mathit{y}\mid {\mathit{\theta }}_k)$. This choice of proposal also makes it easy to implement particle filtering via SIR since one does not need direct access to the dynamic prior which is generally no longer available in closed form. Nonetheless, in terms of the importance sampling step itself the use of the prior as the proposal is generally highly suboptimal. This has motivated more sophisticated particle methods such as those involving tempering (Chopin and Papaspiliopoulos, 2020; Murphy, 2023), hybrid methods (Pirk et al., 2022; Särkkä and Svensson, 2023), and adaptive importance sampling methods (Bugallo et al., 2017). Our focus will be on the latter, although we note in passing that it is entirely possible to combine these approaches (e.g. Koblents and Míguez, 2015).

The idea behind adaptive importance sampling is fairly simple: apply importance sampling iteratively and use the results of the previous iteration to adapt the proposal for the next. The reason that this approach tends to work is also clear: after each iteration the weighted particle ensemble tends to be a better approximation of the posterior than the proposal and so we can use these particles to allow the proposals to evolve by adapting toward the posterior. As such, adaptive importance sampling leverages the power of iterations that iterative ensemble Kalman methods also exploit. At the same time, adaptive importance sampling adapts the proposal to the complexity of the target posterior at hand so that the cost of the algorithm is not set in stone. Thereby, these adaptive algorithms can converge rapidly via early stopping and only in the worst case will they run for a predefined maximum number of iterations. A plethora of adaptive importance sampling methods exist (Bugallo et al., 2017) and many of these are worthy of further investigation in cryospheric DA. The AMIS method that we chose to adopt was based on preliminary non-exhaustive testing of a subset of these methods and came out on top in terms of efficiency and performance.

2.6.2 Adaptive particle batch smoother

Recall that we seek to enhance a popular cryospheric DA algorithm known as the PBS that was introduced by Margulis et al. (2015) and has since been widely adopted for cryospheric reanalysis (e.g. Margulis et al., 2016; Navari et al., 2016; Cortés and Margulis, 2017; Fiddes et al., 2019; Alonso-González et al., 2021; Liu et al., 2021; Girotto et al., 2024; Cao et al., 2025; Sun et al., 2025) by embedding it within the powerful and more general adaptive framework of the AMIS algorithm (Cornuet et al., 2012). We will refer to this proposed new DA method as the adaptive PBS, or AdaPBS for short. In the current AdaPBS implementation, we immediately make one simplification compared to the general AMIS algorithm laid out above which is to fix the number of particles $N_{\mathrm{e}}^{\left(\mathrm{\ell }\right)}$ sampled in each iteration to a constant ensemble size, i.e. $N_{\mathrm{e}}^{\left(\mathrm{\ell }\right)}=N_{\mathrm{e}}$ for all iterations ℓ. The parsimonious choice of a constant ensemble size for all iterations reduces the number of hyperparameters that the user has to specify and facilitates comparisons with the conventional PBS with a given ensemble size. We still suspect that varying the size of the ensemble during the iterations may lead to performance improvements, which is a topic worthy of future research. In our case of a fixed ensemble size at each iteration, the ratio in the DM simplifies to $N_{\mathrm{e}}^{\left(j\right)}/N_{\mathrm{h}}^{\left(\mathrm{\ell }\right)}=\mathrm{1}/\mathrm{\ell }$ since $N_{\mathrm{h}}^{\left(\mathrm{\ell }\right)}=\mathrm{\ell }{N}_{\mathrm{e}}$ so each component of the DM proposal is equally weighted. As such, the typically suboptimal choice of the prior as the initial proposal, which is implicit in the PBS, with a mixture weight of $\mathrm{1}/\mathrm{\ell }$ becomes less influential as ℓ increases and the latter adaptive iterations together occupy a continuously increasing total mixture weight of $(\mathrm{\ell }-\mathrm{1})/\mathrm{\ell }\to \mathrm{1}$ as ℓ→∞.

For simplicity, but without loss of generality, in the current version of AdaPBS we employ multivariate normal (Gaussian) proposal distributions $q^{\left(\mathrm{\ell }\right)}\left(\mathit{\theta }\right)=N(\mathit{\theta }\mid {\mathit{\mu }}^{\left(\mathrm{\ell }\right)},{\mathbf{C}}^{\left(\mathrm{\ell }\right)})$ that are parametrized by hyperparameters in the form of a mean vector μ^(ℓ) and covariance matrix C^(ℓ). In principle, we could set these proposals to anything with at least the same support as the target posterior so we have great freedom in picking our proposal distribution. In practice, however, importance sampling works better if the proposal is similar to the target posterior distribution. As such, it is advantageous to update the hyperparameters of the successive proposal distributions using the current particle approximation of the posterior in the form of the ensemble ${\mathit{\theta }}_r^{\left(\mathrm{\ell }\right)}$ obtained in step 6 of AMIS. Recalling that this is an equally weighted ensemble, an obvious choice is to estimate these hyperparameters using ensemble statistics from the successive posterior approximations. Thus, in step 7 of AMIS the mean vector and covariance matrix for the Gaussian proposal in the next iteration ℓ+1 are simply set to the ensemble mean vector

$$\begin{array}{}\text{(24)}& {\mathit{\mu }}^{(\mathrm{\ell }+\mathrm{1})}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}}}{\mathbf{\Theta }}^{\left(\mathrm{\ell }\right)}\mathbf{1},\end{array}$$

and ensemble covariance matrix

$$\begin{array}{}\text{(25)}& {\mathbf{C}}^{(\mathrm{\ell }+\mathrm{1})}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}}}\left({\mathbf{\Theta }}^{\left(\mathrm{\ell }\right)}-{\mathit{\mu }}^{(\mathrm{\ell }+\mathrm{1})}{\mathbf{1}}^{\mathrm{T}}\right){\left({\mathbf{\Theta }}^{\left(\mathrm{\ell }\right)}-{\mathit{\mu }}^{(\mathrm{\ell }+\mathrm{1})}{\mathbf{1}}^{\mathrm{T}}\right)}^{\mathrm{T}},\end{array}$$

where the particle ensemble ${\mathit{\theta }}_r^{\left(\mathrm{\ell }\right)}$ is stored as column vectors in the N_p×N_e matrix Θ^(ℓ) and 1 is a N_e×1 column vector of ones.

If, using suitable transformations, we further require a prior $p\left(\mathit{\theta }\right)=N(\mathit{\theta }\mid \mathit{\mu },\mathbf{C})$ that is Gaussian and employing the usual Gaussian likelihood Eq. (4) of the PBS $p(\mathit{y}\mid \mathit{\theta })=N(\mathit{y}\mid \widehat{\mathit{y}},\mathbf{R})$ where the mean $\widehat{\mathit{y}}=\widehat{\mathit{y}}\left(\mathit{\theta }\right)$ is the predicted observation vector and R is the observation error covariance then the (unnormalized) target posterior is

$$\begin{array}{}\text{(26)}& p(\mathit{y}\mid \mathit{\theta })p\left(\mathit{\theta }\right)=A\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\varphi }\left(\mathit{\theta }\right)\right),\end{array}$$

where $A=\det \left(\mathrm{2}\mathit{\pi }\mathbf{R}{)}^{-\mathrm{1}/\mathrm{2}}\det \right(\mathrm{2}\mathit{\pi }\mathbf{C}{)}^{-\mathrm{1}/\mathrm{2}}$ is a constant that does not depend on θ and we have introduced the unnormalized and scaled negative log posterior ϕ (not to be confused with the uncertain bounded parameters φ=𝒯(θ)) defined as

$$\begin{array}{}\text{(27)}& \mathit{\varphi }\left(\mathit{\theta }\right)={\left[\mathit{y}-\widehat{\mathit{y}}\right]}^{\mathrm{T}}{\mathbf{R}}^{-\mathrm{1}}\left[\mathit{y}-\widehat{\mathit{y}}\right]+{\left[\mathit{\theta }-\mathit{\mu }\right]}^{\mathrm{T}}{\mathbf{C}}^{-\mathrm{1}}\left[\mathit{\theta }-\mathit{\mu }\right],\end{array}$$

for economy. We may construct a similar expression for the Gaussian proposal distributions $q^{\left(\mathrm{\ell }\right)}\left(\mathit{\theta }\right)=N(\mathit{\theta }\mid {\mathit{\mu }}^{\left(\mathrm{\ell }\right)},{\mathbf{C}}^{\left(\mathrm{\ell }\right)})$ of the form

$$\begin{array}{}\text{(28)}& q^{\left(\mathrm{\ell }\right)}\left(\mathit{\theta }\right)=c^{\left(\mathrm{\ell }\right)}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\psi }}^{\left(\mathrm{\ell }\right)}\left(\mathit{\theta }\right)\right),\end{array}$$

where we have defined the shorthand function

$$\begin{array}{}\text{(29)}& {\mathit{\psi }}^{\left(\mathrm{\ell }\right)}\left(\mathit{\theta }\right)={\left[\mathit{\theta }-{\mathit{\mu }}^{\left(\mathrm{\ell }\right)}\right]}^{\mathrm{T}}{{\mathbf{C}}^{\left(\mathrm{\ell }\right)}}^{-\mathrm{1}}\left[\mathit{\theta }-{\mathit{\mu }}^{\left(\mathrm{\ell }\right)}\right],\end{array}$$

where $c^{\left(\mathrm{\ell }\right)}=\det (\mathrm{2}\mathit{\pi }{\mathbf{C}}^{\left(\mathrm{\ell }\right)}{)}^{-\mathrm{1}/\mathrm{2}}$ is a normalizing constant that only depends on the sampling proposal covariance C^(ℓ). Inserting Eq. (29) in Eq. (17) with fixed N_e the Gaussian DM proposal at iteration ℓ has the following simple analytical form

$$\begin{array}{}\text{(30)}& {\mathit{\upsilon }}^{\left(\mathrm{\ell }\right)}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\ell }}}\sum _{j=\mathrm{1}}^{\mathrm{\ell }}{c}^{\left(j\right)}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\psi }}^{\left(j\right)}\left(\mathit{\theta }\right)\right),\end{array}$$

whereby it is readily verified that a numerically stable estimate of the logarithm of the unnormalized weights in Eq. (18) for the particle history ($k=\mathrm{1},\mathrm{\dots },\mathrm{\ell }$ and $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{e}}$) at iteration ℓ is given by

$$\begin{array}{}\text{(31)}& \log \left({\stackrel{\mathrm{̃}}{w}}_i^{\left(k\right)}\right)=\log \left(A\right)-\mathrm{0.5}\mathit{\varphi }\left({\mathit{\theta }}_i^{\left(k\right)}\right)-{\mathrm{LSE}}_j\left(a^{\left(j\right)}\left({\mathit{\theta }}_i^{\left(k\right)}\right)\right),\end{array}$$

where with some abuse of notation the log-sum-exp term LSE_j(⋅) related to the DM is

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& {\mathrm{LSE}}_j\left(a^{\left(j\right)}\right({\mathit{\theta }}_i^{\left(k\right)}\left)\right)=a^{(\star )}\left({\mathit{\theta }}_i^{\left(k\right)}\right)\\ & +\log \left(\sum _{j=\mathrm{1}}^{\mathrm{\ell }}\exp \left(a^{\left(j\right)}\left({\mathit{\theta }}_i^{\left(k\right)}\right)-a^{(\star )}\left({\mathit{\theta }}_i^{\left(k\right)}\right)\right)\right),\end{array}\end{array}$$

which is implicitly a function with $j=\mathrm{1},\mathrm{\dots },\mathrm{\ell }$ arguments

$$\begin{array}{}\text{(33)}& a^{\left(j\right)}\left({\mathit{\theta }}_i^{\left(k\right)}\right)=\log \left(c^{\left(j\right)}\right)-\log \left(\mathrm{\ell }\right)-\mathrm{0.5}{\mathit{\psi }}^{\left(j\right)}\left({\mathit{\theta }}_i^{\left(k\right)}\right),\end{array}$$

whose maximum across ℓ iterations is $a^{(\star )}\left({\mathit{\theta }}_i^{\left(k\right)}\right)={max}_j\left(a^{\left(j\right)}\left({\mathit{\theta }}_i^{\left(k\right)}\right)\right)$. We emphasize that to obtain the logarithm of the unnormalized weights Eq. (31) for the particle history at each iteration ℓ both the negative log posterior in Eq. (27) and the LSE form of the DM in Eq. (32) have to be evaluated for each of the ℓN_e particles. To obtain the logarithm of the self-normalized weights Eq. (19) we also need to use the logarithm of all the ℓN_e historical unnormalized weights from Eq. (31) to evaluate the logarithm of the evidence approximation Eq. (20) which is given by

$$\begin{array}{}\text{(34)}& \log \left(Z_{\mathit{\upsilon }}^{\left(\mathrm{\ell }\right)}\right)=-\log \left(\mathrm{\ell }{N}_{\mathrm{e}}\right)+{\mathrm{LSE}}_{k,i}\left(\log \left({\stackrel{\mathrm{̃}}{w}}_i^{\left(k\right)}\right)\right),\end{array}$$

where the new log-sum-exp term ${\mathrm{LSE}}_{k,i}(\cdot )$ is implicitly a function with ℓN_e arguments

$$\begin{array}{}\text{(35)}& \begin{array}{rl}& {\mathrm{LSE}}_{k,i}\left(\log \left({\stackrel{\mathrm{̃}}{w}}_i^{\left(k\right)}\right)\right)=\mathrm{\Omega }\\ & +\log \left(\sum _{k=\mathrm{1}}^{\mathrm{\ell }}\sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}\exp \left(\left(\log \right({\stackrel{\mathrm{̃}}{w}}_i^{\left(k\right)})-\mathrm{\Omega }\right)\right),\end{array}\end{array}$$

where $\mathrm{\Omega }={max}_{k,i}\left(\log \left({\stackrel{\mathrm{̃}}{w}}_i^{\left(k\right)}\right)\right)$. Note that Eq. (34) only needs to be evaluated once per iteration ℓ. Finally, combining Eqs. (31) and (34), the stable logarithm of the self-normalized weights Eq. (19) that we seek can now be computed as

$$\begin{array}{}\text{(36)}& \log \left(w_i^{\left(k\right)}\right)=\log \left({\stackrel{\mathrm{̃}}{w}}_i^{\left(k\right)}\right)-\log \left(\mathrm{\ell }{N}_{\mathrm{e}}\right)-\log \left(Z_{\mathit{\upsilon }}^{\left(\mathrm{\ell }\right)}\right),\end{array}$$

which we have expressed in terms of the log evidence approximation Eq. (34) since this is a useful additional output from AdaPBS for hierarchical Bayesian inference (Robert, 2007; Murphy, 2023). The self-normalized weights for AdaPBS are now obtained trivially by taking the exponential of Eq. (36). These weights can then be resampled N_e times to obtain an equally weighted ensemble of particles as described in step 6 of the AMIS algorihtm. Next, as described in step 7, upon convergence these N_e resampled particles are output as a particle approximation of the posterior. Otherwise we use these particles to design the Gaussian sampling proposal through Eqs. (24) and (25) for the next iteration ℓ+1.

The overall workflow of the AdaPBS algorithm is visualized in Fig. 1. As with the standard PBS, AdaPBS is initialized by sampling a particle ensemble from the prior, running this ensemble through the forward model, and computing self-normalized importance weights via Eq. (36). Ensemble diversity is then diagnosed using the effective sample size N_eff in Eq. (22). If N_eff exceeds a predefined threshold, the algorithm stops successfully, outputting this diverse ensemble as the posterior approximation. Otherwise, it adapts a new sampling proposal using clipping Eq. (23) and ensemble statistics via Eqs. (24) and (25). A new particle ensemble is drawn from this adapted sampling proposal, run through the forward model, and weights are obtained using Eq. (36), now incorporating the entire particle history via a mixture proposal Eq. (17). Adaptation continues until the ESS-based diversity threshold is met or a predefined maximum number of iterations is reached. Workflow of the adaptive particle batch smoother (AdaPBS), illustrating its extension (green) from the non-iterative PBS method (yellow). The algorithm iterates by adapting the proposal until the effective sample size (N_eff) reaches a predefined fraction of the ensemble size (N_e), determined by the diversity threshold (τ), or until a maximum number of iterations (ℓ_max) is reached. AdaPBS is mathematically equivalent to PBS if the diversity threshold is met in the first iteration (ℓ=1) or if the algorithm is constrained to a single iteration (ℓ_max=1).

Figure 1Flowchart showing the workflow in the adaptive particle batch smoother (AdaPBS) as an extension (green) of the non-iterative PBS method (yellow). The AdaPBS algorithm continues iterating by adapting the proposal until convergence, achieved when the effective sample size (N_eff) reaches at least a fraction of the ensemble size (N_e) determined by the desired diversity threshold (τ), or until a pre-defined maximum number of iterations (ℓ_max) is reached. Both in the special case of convergence in the first iteration (ℓ=1) and in the non-iterative case (ℓ_max=1), AdaPBS reduces to PBS.

Download

2.7 Probabilistic evaluation diagnostics

To evaluate the performance of the respective data assimilation schemes we employ two probabilistic evaluation diagnostics, namely the Continuous Ranked Probability Score (CRPS; Gneiting et al., 2005) and the Kullback–Leibler Divergence (KLD; Murphy, 2023). These probabilistic verification measures allow us to evaluate the performance of the entire approximate posterior distribution as represented by an ensemble of particles rather than just point estimates such as the posterior mean.

The CRPS is a strictly proper scoring rule defined by (Gneiting and Raftery, 2007)

where $P\left(x\right)={\int }_{\mathrm{\infty }}^{x}p\left(\mathit{\chi }\right)\mathrm{d}\mathit{\chi }$ is the cumulative (posterior or prior) probability distribution of the uncertain variable or parameter x of interest with density p(x) that we are inferring, x^⋆ is the reference truth value, and H(⋅) is the Heaviside function which is 1 for positive arguments (here $x>x^{\star }$) and 0 otherwise. The CRPS is non-negative and inherits the same units as the variable x of interest. In the special case of a deterministic forecast, where P(x) is also a Heaviside function, the CRPS reduces to the absolute error. As such, the CRPS is negatively oriented with the best possible value being 0 indicating perfect agreement with the reference. In the general non-degenerate probabilistic setting the CRPS measures both the accuracy (goodness of fit) and the precision (sharpness) of the distribution p(x) relative to the reference x^⋆, making it an apt choice for evaluating ensemble-based data assimilation experiments in cryospheric science (Piazzi et al., 2018; Cao et al., 2025). In practice, in line with previous studies (Alonso-González et al., 2023; Mazzolini et al., 2025; Alonso-González et al., 2026), we estimate the CRPS for probabilistic predictions obtained via ensemble-based DA schemes in MuSA by assuming that the marginal predictive prior and posterior snow depth distributions follow a normal distribution. Thereby, we only need to store the mean and standard deviation of the ensemble predictions of each scheme which can be plugged into the simple analytical expression for the Gaussian CRPS given by Gneiting et al. (2005).

The CRPS is used to evaluate the performance of the ensemble of snowpack states obtained in MuSA by comparing them to assimilated or independent observations. In a similar vein, we employ the KLD as a probabilistic evaluation diagnostic that quantifies how close the approximate posterior distributions over parameters are to the reference posterior obtained through MCMC simulation using the RAM method. In particular, we use the so-called reverse KLD defined by (Murphy, 2023)

where q(θ|y) is an approximating distribution of the target posterior distribution p(θ|y). The KLD, also known as relative entropy (MacKay, 2003), is a dimensionless measure of the distance between the two distributions that it takes as arguments. As with the CRPS, it is negatively oriented and non-negative with the best possible value being 0 in the case that the input distributions used as arguments are equal. It happens to be asymmetric in these arguments, and this so-called reverse form of the KLD is widely used as an objective function in variational Bayesian inference which recasts inference as an optimization problem over tractable variational distributions q(θ|y) used to approximate the intractable posterior p(θ|y) (Murphy, 2023). While this use of reverse KLD as an objective function motivated our use of this metric, we use it in the slightly different context of evaluation rather than optimization. In particular, we use the reference posterior samples simulated via MCMC using RAM to define the target posterior distribution p(θ|y) and we use the reverse KLD to measure how close the more tractable approximations q(θ|y) from the ensemble-based DA schemes, including the AdaPBS, are to this target. To simplify the calculation of the KLD, especially given that we only have samples from the distributions p and q, we focus on the reverse KLD of the marginal distributions of the perturbation parameters θ in the vector θ in transformed space as introduced in Sect. 3.2. Note that the KLD is invariant to transformations (Murphy, 2023), and by sticking to the transformed unbounded space we can better approximate these marginals as Gaussian. Assuming Gaussianity, the marginal reverse KLD has a simple analytical solution of the form (see supplementary material 5.1.2 of Murphy, 2023)

where μ_q and σ_q are the mean and standard deviation, respectively, of the ensemble-based approximation ${\mathit{\theta }}_i\sim q\left(\mathit{\theta }\right|\mathit{y})$. In addition to computing reverse KLD of the various ensemble-based DA schemes in MuSA we also compute the reverse KLD of the marginal prior, i.e. setting $q\left(\mathit{\theta }\right|\mathit{y})=p(\mathit{\theta })$, as an additional benchmark to help contextualize this distance measure. For all these reverse KLD calculations, we use the MCMC samples obtained via RAM (Sect. 2.3) to define the reference marginal posterior distribution ${\mathit{\theta }}_j\sim p\left(\mathit{\theta }\right|\mathit{y})$ with reference sample mean μ_p and standard deviation σ_p.

3.1 Data assimilation framework and experimental design

All experiments presented in this study were developed using the Multiple Snow Data Assimilation System (MuSA) (Alonso-González et al., 2022). MuSA is a versatile, open-source data assimilation tool designed to integrate diverse snowpack observations with an ensemble of numerical snowpack simulations. Its modular architecture allows users to implement different algorithms and snowpack models. In addition to the set of algorithms already available in MuSA, we have implemented three new algorithms to develop the experiments proposed in this work. This includes the aforementioned AdaPBS algorithm and two MCMC algorithms, namely RWM and RAM, where RAM is used as a sta. The updated MuSA code is released as a new version of the tool (Alonso-González and Aalstad, 2025). In order to provide an intuitive understanding of the behavior of the AdaPBS, we have performed several experiments that compare its performance to that of previously proposed ensemble-based cryospheric DA algorithms.

3.2 Experiment 1: Assimilating drone-based snow depth in a temperature index model

First, we assimilated snow depth observations in an ensemble of simulations generated by a simple temperature index model (Hock, 2003) implemented in MuSA. We emphasize that the primary objective of these experiments is to assess the performance of the AdaPBS algorithm by exploiting a simpler yet computationally efficient model and not to achieve the most physically complex simulations possible. Moreover, in part due to their simplicity in terms of input data and few calibration parameters, temperature index models are being used effectively in operational settings (Lussana et al., 2018) as well as for hemispheric-scale snow reanalysis (Elias Chereque et al., 2024) and global glacier modeling (Rounce et al., 2023). Our temperature index forward modeling implementation relies on the strong assumption of constant density, fixed at a climatological value of 300 kg m⁻³, along with a seasonally constant hourly temperature index factor a set to 0.1375 mm h⁻¹ K⁻¹. Typically a would also be treated as an uncertain parameter (e.g. Rounce et al., 2023), but to facilitate the comparison of schemes and visualizing the results we chose to fix it here. Thus, the near surface air temperature T_n [K] dependent snowmelt rate M_n [mm h⁻¹] at each hourly timestep n is computed as follows:

where the melting temperature T₀=273.15 [K] is set to the freezing point of water. Often this is treated as a calibration parameter to account for errors in the air temperature forcing, but we do this more explicitly by perturbing the air temperature field itself with an additive bias parameter b that we seek to infer. By combining the snowmelt rate with the snowfall rate S_n [mm h⁻¹], the hourly SWE D_n [mm] is updated via

where c is the precipitation perturbation parameter that we also seek to infer. Note that the air temperature bias b also affects the snowfall rate S_n directly because the precipitation phase is diagnosed using the perturbed air temperature T_n+b.

The temperature index model DA experiments were conducted in the Izas experimental catchment, in the Central Spanish Pyrenees (Revuelto et al., 2017). The snow depth data were acquired through drone surveys utilizing structure-from-motion techniques (Revuelto et al., 2021). Meteorological forcing data in the form of air temperature and total precipitation (P), were generated using the Micromet downscaling tool (Liston and Elder, 2006b), driven by the ERA5 reanalysis (Hersbach et al., 2020). These experiments were conducted for a single grid cell from which we retrieve both the forcing and observations using a grid spacing of 5 m spatial resolution. The cell was located in a concavity within the basin and therefore exhibited significant snow accumulations resulting from snow redistribution that makes it particularly challenging to simulate for the temperature index model with standard forcing. These drone-acquired data and forcing datasets have previously been used in data assimilation studies (Alonso-González et al., 2022, 2023) Following these studies and the findings of Revuelto et al. (2021), we assume an observation error standard deviation of σ_y=0.2 m for the drone-based snow depth retrievals. The corresponding observation error variance ${\mathit{\sigma }}_y^{\mathrm{2}}$ is used to construct the diagonal observation error covariance matrix R in the likelihood Eq. (4). The ensemble of simulations was generated by perturbing the precipitation and temperature timeseries with the pseudo-randomly sampled constant in time prior perturbation parameters b and c covering the 2018/2019 snow season. The data assimilation window covered the whole snow season and, as such, our experiments were developed using batch smoothers, but we would still expect similar relative performance (i.e. ranking) of these cryospheric DA schemes when applied as filters.

For simplicity, and without loss of generality for the DA, we assume a priori parameter independence, such that the joint prior factorizes into the product of its marginal prior distributions. This assumption can be relaxed using additional background knowledge on parameter dependence (e.g., Pirk et al., 2022; Alonso-González et al., 2023). Yet, even with prior independence, correlation between parameters can be inferred in the posterior via the likelihood (Cao et al., 2025). The prior probability distribution for the additive air temperature bias is a normal distribution with mean μ_b=0 and standard deviation σ_b=1. For the prior on the multiplicative precipitation scaling, we use a strictly positive lognormal distribution with an associated normal mean μ_c=0.1 (so the actual median perturbation in model space is exp (μ_c)=1.11) and standard deviation σ_c=0.5. Note that we carry out inference directly on the unbounded space of the log-transformed precipitation perturbation parameter and then apply an exponential transform to map this back to the model parameter c following common practice (Gelman et al., 2013; Alonso-González et al., 2022). The location parameters of these prior distributions were chosen in a conservative manner, i.e., near unbiased, with a mean temperature bias of 0 and a precipitation scaling centered close to 1. This corresponds to quite a general case where, without additional background knowledge or data, the orientation of the forcing bias is not known a priori. The relatively large prior spread parameters (σ_c and σ_b) encode considerable initial epistemic uncertainty in the two forcing perturbation parameters.

Here we have compared the results of the approximate Bayesian inference using different DA algorithms, to evaluate and benchmark the performance of the novel AdaPBS algorithm. First, we have solved the data assimilation problem using adaptive MCMC using RAM (Vihola, 2012), which we consider a gold-standard albeit computationally costly reference. Taking inspiration from Emerick and Reynolds (2011), we initialized the Markov chain used in RAM at the approximate posterior mean obtained from ES-MDA with N_a=4 iterations and N_e=100 ensemble members. The idea is to accelerate MCMC by starting the chain within or close to the typical set of the target posterior distribution. The length of the Markov chain in the RAM algorithm was set to 20 000 and the initial 10  % of the samples were discarded as burn in. Discarding the burn in period avoids initialization artifacts whereas the relatively long chain gives the Markov Chain time to reach its target posterior distribution.

We then performed the same exercise using a Particle Batch Smoother (PBS), an ES, an ES-MDA and finally we ran the AdaPBS. The hyperparameters of the analysis (i.e. the N_e=100 number of ensemble members and the prior distributions) were the same in all cases. In the case of the AdaPBS, the N_eff limit was set to 30  %, which means that for our N_e=100 member ensemble at least 30 particles should show non-negligible weights before stopping the iterations. The results of each algorithm were compared with the samples from MCMC-RAM which we used as a reference posterior. The agreement between the approximate posteriors obtained via ensemble-based DA and the MCMC reference was quantified using the KLD as a probabilistic evaluation diagnostic as described in Sect. 2.7.

3.3 Experiment 2: Assimilating hourly ESM-SnowMIP snow depth data in FSM2

In the second set of experiments, we conducted several site-level simulations using a physics-based snow model of intermediate complexity, the Flexible Snow Model (FSM2; Essery et al., 2025). In these experiments, we assimilated hourly observations of snow depth at six different geographical locations using both AdaPBS and ES-MDA for comparison. For simplicity, we assume the same observation error standard deviation of σ_y=0.2 m as in Experiment 1. While this is likely to be an overestimate of individual measurement errors at these sites, inflating independent observation errors in our model helps implicitly account for the fact that the actual errors in these data are correlated in time, which reduces their information content (Evensen et al., 2022). Assimilating snow depth observations at different locations over different periods of time allowed us to obtain a general overview of the behavior of the algorithms in different climates and weather regimes. Observations were obtained from a database generated under the framework of the Earth System Model-Snow Model Intercomparison Project (ESM-SnowMIP; Ménard et al., 2019). Specifically, we conducted simulations at Reynolds Mountain East, Idaho, USA (RME); Sapporo, Japan (SAP); Senator Beck, Colorado, USA (SNB); Sodankylä, Finland (SOD); Swamp Angel, Colorado, USA (SWE) and Weissfluhjoch, Switzerland (WFJ). The meteorological forcing was obtained directly from the nearest grid cell of the ERA5 reanalysis (Hersbach et al., 2020) without applying topographic corrections. While site-level meteorological forcing was available from Ménard et al. (2019), we used the more uncertain, unadjusted global reanalysis data to provide a more challenging experimental setup that is generally applicable to any site with snow depth data. Specifically, each water year, we assimilated a large batch of $\mathrm{24}\times \mathrm{365}=\mathrm{8760}$ hourly snow depth observations.

To explore the performance of the algorithms in the case of a higher dimensional parameter space, we have generated 7 perturbation parameters to correct the forcing (air temperature, precipitation, relative humidity, surface pressure, wind speed, and incoming long- and short-wave radiations). On top of this, we have also included the 12 main tunable internal model parameters related to snow (not vegetation) in FSM2 (Essery et al., 2025) as listed in Table 1. The prior distributions for the temperature and precipitation parameters are the same as those for the previous experiment described in Sect. 3.2. For the new model parameters, we used a multiplicative logit-normal distribution bounded in the physical space between 0.8 and 1.2, defined by the moments of its underlying normal distribution μ=0 and σ=1. As such, we perturb these new parameters in the ±20  % range from their reference values in the physical space, ensuring that the combinations parameters remain in their physical bounds, maintaining the numerical stability of FSM2. We performed a dependent validation to compare the performance of AdaPBS with ES-MDA and PBS, by comparing the posterior snow simulations with the assimilated observations, using bias and root mean squared error (RMSE) as well as an uncertainty-aware metric in the form of the continuous ranked probability score (CRPS) described in Sect. 2.7. We removed the timesteps where both the simulation and observations had snow depths of zero in the computation of the error and CRPS to avoid artificially deflating the validation metrics by not considering the seasonality of the snow depth values. That is to say, we are primarily concerned with model performance in the periods in which the model or observations indicate the presence of a seasonal snowpack. It is only in these periods that the model is susceptible to snow commission or omission errors rather than the mostly trivially correct no-snow prediction outside the snow season.

Unlike Experiment 1 with the temperature index model, we did not run MCMC in this more complex experiment due to the high computational cost and marginal additional benefits to the analysis. As such, the KLD metric could not be used without an MCMC-based reference distribution. Instead, this experiment focused on comparing the predictive performance in state space of AdaPBS to that of ES-MDA, an advanced and competitive baseline for DA with more complex snow models (Alonso-González et al., 2022, 2026). In this higher-dimensional experiment, the entire 19-dimensional parameter space is too large for comprehensive and instructive visualization. An alternative could be to investigate the importance of each parameter by estimating the information gain (van Hove et al., 2026), but such parameter analysis strays beyond the focus of this experiment. Except for a handful of forcing parameters partly analyzed in Experiment 1, the majority of the 19 parameters are so-called nuisance parameters (Jaynes, 2003). They are not of primary interest beyond their effect on the state and their role in making the DA problem more challenging. We thus focus on the performance of AdaPBS and ES-MDA for estimating the snow depth state observed at the selected ESM-SnowMIP sites. The non-iterative PBS is also tested in this experiment to directly assess the impact of the adaptations within the iterative AdaPBS.

Table 1Tabular overview of the internal FSM2 parameters that we perturbed at ESM-SnowMIP sites. Columns from left to right: the parameter name in the FSM2 code, its default value, physical units (– is dimensionless), and a description adapted from Essery (2015).

Download Print Version | Download XLSX

4.1 Inter-comparison of algorithms using a temperature index model

Despite relying on identical priors, observations, observation error model, and forward model in the form of a temperature index snow model, the performance of the DA algorithms differed substantially. This conclusion aligns with the few previous publications that developed inter-comparison experiments on ensemble-based cryospheric data assimilation algorithms (Leisenring and Moradkhani, 2011; Margulis et al., 2015; Aalstad et al., 2018; Alonso-González et al., 2022; Cao et al., 2025) and highlights the importance of selecting the algorithm with careful consideration of the problem at hand by balancing accuracy and computational cost.

Given the highly informative nature of the drone-based snow depth observations, the posterior approximation obtained through MCMC samples via RAM showed a challenging, narrow, banana-shaped distribution (Fig. 2). This demonstrates that even in this relatively simple case – where in practice we are only updating 2 parameters (b and c) by assimilating just a handful (5) of snow depth observations – efficiently sampling from the posterior is nonetheless a challenging task. This is probably exacerbated in this case by the fact that we are using batch smoothers. On the one hand, such smoothers can provide more information than filters by assimilating the complete trajectory of observations in the water year all at once and propagating information backwards in time. On the other hand, the larger immediate information gain from batch smoothers can make posterior sampling more challenging than with sequential filters.

The original PBS algorithm of Margulis et al. (2015) has several benefits. It is relatively easy to implement, relies on few assumptions, and its computational cost is lower than that of other schemes, even when considering the same number of members in the ensemble since (unlike the ES) it does not require any reruns. These benefits do come at a cost, as PBS is prone to collapse in more challenging settings. This means that sometimes, especially with very informative (i.e. numerous and/or accurate) observations, the ensemble collapses and degeneracy ensues (Snyder et al., 2008; Morzfeld et al., 2017). Using a familiar analogy, this is like looking for smaller needles in a haystack. A similar result arises if the priors are relatively far (as measured, e.g., by a modified variant of the KLD in Sect. 2.7) from the posterior in the parameter space, such as in the case of overconfident and/or strongly biased priors. More generally, “far” is not just the distance between the prior and posterior mean, but whenever the overlap between the prior and the posterior is small, including when a broad and/or high-dimensional prior encompasses a narrow posterior. Following a similar analogy, this is akin to having a larger haystack. These needle and haystack effects can yield sub-optimal solutions with high Monte Carlo variance (Robert and Casella, 2004) that make it difficult, if not impossible, to accurately estimate uncertainty. So even if SNIS, which powers the PBS, is asymptotically unbiased (Rainforth et al., 2020), the PBS in particular and basic importance sampling in general may exhibit a large variance since the prior is implicitly used as the proposal, which is generally far from the posterior (MacKay, 2003). In the worst case, this will result in complete particle degeneracy with N_eff≃1. This was clearly the case in our particular experiment, as shown in Fig. 2, where the majority of the weight collapsed onto a single particle despite the fact that the range of the prior parameter distribution clearly encompassed the reference posterior MCMC samples obtained via the RAM algorithm. In the model state space, the fact that a single particle carries all the probability leads to a degenerate posterior snow depth ensemble that is not close to the assimilated observations and completely lacks well-calibrated uncertainty quantification.

In this case, the skewed prior on the precipitation perturbation combined with the challenging banana shape of the posterior are key challenges leading to degeneracy since very few particles were sampled in the range of the typical set of the posterior. The situation would improve if a more informative prior were used with a higher median precipitation perturbation parameter and even a positive correlation between the perturbations to better sample the banana shape. The problem, of course, is that these prior hyperparameter settings are rarely, if ever, known a priori, even from expert knowledge, and setting them directly based on the data leads to incoherent “double-dipping” via so-called empirical Bayes (Robert, 2007; Gelman et al., 2013; Murphy, 2023). Adaptive particle methods present a solution that, instead of incoherently changing the prior to avoid degeneracy, simply adapts the proposal towards the target posterior so as to perform more sample efficient approximate Bayesian inference that remains coherent.

It may seem surprising that importance sampling via PBS struggles to accurately capture the reference posterior in this apparently simple experiment with a two-dimensional parameter space. An alternative would be gradient-based optimization techniques, such as gradient descent or Newton's method, used in both machine learning (Murphy, 2023) and variational DA (Bannister, 2017), to identify the mode with relatively few iterations and find a Laplace approximation of the posterior (MacKay, 2003). However, as is often the case in practical DA in the absence of automatic differentiation (Murphy, 2023), we do not have access to gradients and are focused on the harder problem of “black-box” (i.e., gradient-free) inference that is naturally handled by ensemble-based methods (Sanz-Alonso et al., 2023). While similar low-dimensional calibration problems are routinely encountered in hydrology and can be tackled via methods such as generalized likelihood uncertainty estimation (GLUE; Beven and Binley, 1992), SCEM-UA optimization (Vrugt et al., 2003), or DREAM MCMC sampling (Vrugt et al., 2008), Nonetheless, these methods usually involve orders of magnitude more forward model evaluations than the N_e=100 particles used by PBS in this experiment. For example, Clark et al. (2006) demonstrated the need for a computationally intensive SODA approach, using SCEM-UA in an outer loop to optimize parameters, together with an EnKF inner loop for state estimation, to rigorously calibrate a simple two-parameter snow model. Furthermore, the widely used hydrological calibration technique GLUE (Beven and Binley, 1992), tailored to identify equifinality (Beven, 2006), is equivalent to PBS when using a Gaussian likelihood. As such, GLUE could only outperform PBS by increasing the ensemble size or modifying the likelihood, which are equally applicable to PBS. Both black-box optimization and inference are challenging problems where particle methods have proven to be highly competitive (Duan et al., 2023; Chopin and Papaspiliopoulos, 2020). In fact, these particle methods can be seen as formalizing highly successful meta-heuristic evolutionary algorithms that are widely used for optimization (Del Moral, 2004). The takeaway here is that particle methods can be improved either by increasing the ensemble size, at considerable computational expense given typically poor scaling (Snyder et al., 2008), or by adapting them for example using various iterative methods (Bugallo et al., 2017; Chopin and Papaspiliopoulos, 2020).

Figure 2Results from the particle batch smoother (PBS) assimilating drone-based snow depth observations in a temperature index model for a single cell at Izas in water year 2019: The model parameter space visualized as the temperature bias parameter b on the y-axis and the precipitation correction c on the x-axis with the prior ensemble of particles shown with red dots while 2 green stars with size proportional to weight indicate particles with non-negligible PBS weights (N_eff=1.26) and the blue banana-shaped distribution shows the reference posterior MCMC samples obtained via RAM (left panel); The model state space (right panel) showing the trajectory of snow depth with mean (solid line) ±1 standard deviation (shading) from the prior (orange), PBS posterior (green), and reference MCMC posterior (blue), along with the assimilated observations (yellow dots).

Download

The ES (van Leeuwen and Evensen, 1996), i.e. the batch smoother version of the EnKF (Evensen, 1994), has also been widely used in cryospheric DA particularly in reanalysis settings (Durand et al., 2008; Girotto et al., 2014). Although the EnKF has mainly been used to directly update snow model states (e.g. De Lannoy et al., 2012), the ES is better suited to update the forcing and internal parameters in a forcing formulation of the data assimilation problem (Evensen et al., 2022). Since the update moves the parameters in parameter space through the ensemble Kalman analysis step, with the ES it is necessary to generate an ensemble of snow model simulations twice, once with the prior parameters and once with the posterior parameters. This has the obvious advantage of maintaining model stability and physical consistency as long as the parameters remain within their physical bounds (Alonso-González et al., 2022). The latter boundedness can be easily controlled using transformation techniques such as Gaussian anamorphosis that also help to accommodate the Gaussian assumption (Bertino et al., 2003; Aalstad et al., 2018). However, with the ES this forcing formulation approach entails a non-parallelizable doubling of the computational cost (2N_e simulations) compared to solving the same problem with the same ensemble size using the PBS (N_e simulations). In practice, the fact that the ES is less prone to ensemble collapse allows one to reduce the ensemble size compared to PBS such that in some problems the ES may actually end up being less costly. This gives a lot of flexibility when configuring the assimilation system according to the available computational resources especially in spatio-temporal settings where localization methods are more easily implemented with ensemble Kalman schemes (Evensen et al., 2022). Despite these advantages, the linear assumption in combination with models used in cryospheric sciences, which are typically non-linear, can lead to highly suboptimal posterior approximations with poorly calibrated and even underconfident uncertainty estimates. In fact, this is what we see in Fig. 3 with slight and underconfident updates both in parameter and state space, in line with previous studies (Aalstad et al., 2018; Alonso-González et al., 2022).

Figure 3Analogous to Fig. 2 but for the ensemble smoother (ES) with the prior (top left) and posterior (top right) in the model parameter space and the trajectory of snow depth (bottom) in the model state space.

Download

The typical underlying linear assumption of ES, which is the main cause of suboptimal results, can be strongly violated in practical cryospheric DA problems. The iterative ES-MDA (Emerick and Reynolds, 2013) has proven to be a powerful cryospheric DA algorithm in the few previous inter-comparisons that have been carried out to date (Aalstad et al., 2018; Alonso-González et al., 2022). Iterations allow for a progressive movement of parameters towards the posterior via likelihood tempering (Chopin and Papaspiliopoulos, 2020; Murphy, 2023), which in practice limits the effects of model non-linearity as shown in Fig. 4. One limitation is that the number of iterations must be pre-set in the ES-MDA of Emerick and Reynolds (2013), making it an important hyperparameter that must be carefully adjusted balancing between computational cost and solution robustness. When implemented in a distributed cell-by-cell manner, this is problematic, as it is difficult to adjust the number of iterations cell by cell, opting in practice for a fixed number of iterations for the whole domain. Such a global “one size fits all” approach is very likely to result in a waste of computational resources through an unnecessary number of model realizations compared to a locally adaptive approach. In addition, the computational cost of the linear algebra operations needed for ES-MDA may become non-negligible in the case of assimilating a large number of observations.

Figure 4Analogous to Fig. 3 but for an iterative ensemble smoother (IES) in the form of the ensemble smoother with multiple data assimilation (ES-MDA). The top panels show the evolution from the prior to the posterior ensemble members (red) across the MDA iterations along with the reference banana-shaped MCMC posterior samples (blue). As before, the bottom panel shows the predicted snow depth in model state space, but with a better calibrated IES posterior (green). The zero-based iteration counter from 0 to 4 corresponds to ℓ in the algorithm description in Sect. 2.4 with N_a=4 iterations of the update step and $N_{\mathrm{a}}+\mathrm{1}=\mathrm{5}$ iterations of the prediction step.

Download

Here we introduce a new algorithm, namely the AdaPBS, as an effort to combine the advantages of the different algorithms presented above into a single tool. The AdaPBS has the potential to be a powerful iterative DA method, aspiring to sample from the posterior with performance similar to that of ES-MDA. Unlike the ES-MDA, AdaPBS does not rely on assumptions of linearity or Gaussianity, which facilitates its implementation, especially for users with limited experience in data assimilation. This also opens up for the possibility of using tailored likelihood functions such as those involving zero-inflation (Smith et al., 2010; Tang et al., 2023) that may be better suited to double-bounded cryospheric variables.

Figure 5Analogous to Fig. 4 but for the adaptive particle batch smoother (AdaPBS) showing the particle and effective sample size (ESS, i.e., N_eff) evolution across iterations. Note that the results after the initial iteration (top left) correspond exactly (ignoring Monte Carlo variance) to those that would be obtained from the PBS in Fig. 2. The zero-based iteration counter from 0 to 4 corresponds to ℓ−1 in the algorithm description in Sect. 2.6, so here the AdaPBS converged (here N_eff>30) in ℓ=5 iterations with exactly the same computational cost as the ES-MDA with $N_{\mathrm{a}}+\mathrm{1}=\mathrm{5}$ iterations of the prediction step.

Download

A clear computational advantage of the AdaPBS over ES-MDA is that it does not require pre-selecting a fixed number of iterations, enabling the use of early stopping strategies throughout the simulation domain. This allows the number of iterations to be automatically adapted for each cell in a spatio-temporally distributed model domain, depending on the difficulty of the local (or localized) inference problem being solved. We have used an early stopping criterion that stops the iterations when a threshold on the classical estimate of the ESS (N_eff in Eq. 22), i.e., a minimum number of ensemble members with considerable weight is reached, making it more robust to collapse than PBS. In this case, we used a relatively high diversity threshold of τ=0.3 for $N_{\mathrm{eff}}/N_{\mathrm{e}}$, which in practice resulted in the same number of iterations as ES-MDA. However, it should be noted that in this example, a respectable N_eff=22 was already achieved by iteration 3 as shown in Fig. 5, which may be sufficient for many applications.

Thus, a key hyperparameter to set in AdaPBS is the diversity threshold τ, which controls the target minimum effective sample size N_eff. Since τ is directly related to the desired performance of the algorithm, it is more interpretable for users than the fixed number of iterations N_a required in ES-MDA (Emerick and Reynolds, 2013). This makes AdaPBS an easy-to-implement algorithm, resistant to collapse, and potentially able to significantly reduce computational cost compared to ES-MDA when applied to large domains with many cells. Furthermore, even for very challenging cases with AdaPBS struggling to reach the threshold diversity τ, the algorithm will still terminate after a predefined maximum number of iterations ℓ_max which can be allocated depending on the computational budget at hand. Even in these challenging edge cases we have found that although the particle diversity and resulting uncertainty quantification may be worse than desired, the AdaPBS method will by construction still typically provide more accurate and reliable inference than the PBS.

In Table 2 we compare the KLD distance measures for the marginal posterior distributions obtained for temperature and precipitation parameters using the respective algorithms. This shows that although AdaPBS is not capable of reaching precisely the same low KLD values as the ES-MDA, it shows quite a comparable performance to this state-of-the art iterative ensemble Kalman method. In fact, the quantitative differences appear minimal when the distributions are compared graphically in Fig. 6. More generally, it is clear that the iterative methods, namely the AdaPBS and ES-MDA, are much closer to being able to match the performance of MCMC than their non-iterative counterparts, namely the PBS and ES. Crucially, the AdaPBS can adapt to the complexity of the problem at hand and can thus be guaranteed to incur a lower computational or in the worst case equal cost to the ES-MDA by setting ${\mathrm{\ell }}_{max}=N_{\mathrm{a}}+\mathrm{1}$. This cost difference is highlighted in the cost column of Table 2, where in the best case with a less complex target posterior the AdaPBS has the chance to trigger early stopping already after a single iteration (L=1) and computationally affordable as the PBS on which it is based. Moreover, like the PBS, the AdaPBS can handle more general probabilistic models than the Gaussian models assumed by the ES-MDA and other ensemble Kalman methods.

Table 2Performance of cryospheric DA algorithms in terms of inferential accuracy and cost. Accuracy is gauged in terms of how well the posterior approximations q from the algorithms match the reference posterior estimate p from RAM as measured by the Gaussian approximation of the marginal reverse KLD $D_{\mathrm{KL}}\left(q\right|\left|p\right)$ in Eq. (39) for air temperature bias b and the precipitation correction c in the temperature index model experiments. Computational cost is measured in terms of the number of forward model runs as multiples of the ensemble size N_e=100, ES-MDA iterations N_a=4, AdaPBS iterations L≤ℓ_max where typically ${\mathrm{\ell }}_{\max }=N_{\mathrm{a}}+\mathrm{1}$, and the number of Markov chain steps in RAM $N_{\mathrm{s}}=\mathrm{20}\times {\mathrm{10}}^{\mathrm{3}}$.

Download Print Version | Download XLSX

Figure 6Comparison of Gaussian kernel density estimates (KDEs) obtained via scipy.stats.gaussian_kde (Scott's rule bandwidth; Virtanen et al., 2020) of the marginal prior and marginal approximate posterior distributions of the parameters from the respective algorithms in the transformed (unbounded) space. The MCMC posterior (orange) from RAM is considered the gold-standard.

Download

It should also be noted that the estimated posterior uncertainty differs between ES-MDA and AdaPBS. The posterior standard deviation in ES-MDA is slightly higher compared to that of AdaPBS. At least in the present example, these differences are nonetheless relatively minor, especially in the model state space visualized here in terms of snow depth. While quantitative probabilistic evaluation using diagnostics such as KLD is helpful, qualitative visual comparisons to an MCMC reference posterior ensemble trajectory in state space can be equally instructive. Inspecting the snow depth trajectories in Figs. 2–5 shows that the patterns identified in parameter space carry over to snow depth in state space. Using the RAM-based MCMC posterior trajectory as a reference: PBS is biased and highly degenerate, ES is biased and highly under-confident, ES-MDA is unbiased yet slightly under-confident, whereas AdaPBS overlaps almost perfectly with the reference ensemble. To date, most snow DA studies focus primarily on the predictive performance of point estimates such as the posterior mean, median, or mode (e.g. Margulis et al., 2016; Aalstad et al., 2018; Alonso-González et al., 2021; Oberrauch et al., 2024). Nonetheless, as snow DA continues to transition towards fully Bayesian methods, an emphasis on rigorous uncertainty quantification will become increasingly important for probabilistic prediction (Reich and Cotter, 2015) and decision-making (Robert, 2007). As we have shown, benchmarking via MCMC reference posteriors (Law and Stuart, 2012) allows for direct qualitative and quantitative evaluation of more approximate posterior estimates obtained from tractable DA algorithms such as AdaPBS.

Despite the many advantages of the AdaPBS scheme that we have highlighted, the ES-MDA still retains a key advantage: There are more localization methods developed for ensemble Kalman-based algorithms (Evensen et al., 2022) than for those based on particle methods (Farchi and Bocquet, 2018). This facilitates the development of spatiotemporal assimilation initiatives (Alonso-González et al., 2023; Mazzolini et al., 2025; Alonso-González et al., 2026), as opposed to the purely temporal example presented here, which allows non-local observations in remote cells of a distributed model domain to be considered to update the local cell in question. In this way, information can be propagated in space, correcting areas of the domain even when they have no local observations, or integrating point-scale observations into distributed simulations. A future line of research with great potential will be the development of these spatio-temporal methods for AdaPBS and other particle-based methods, which remains at the frontier of current DA research. Moreover, in conjunction with the localization problem, it remains to be seen if adaptive particle methods can be scaled up to as high-dimensional problems as ensemble Kalman methods (Carrassi et al., 2018; Pirk et al., 2024). Nonetheless, the results from our inter-comparison of cryospheric DA schemes with a temperature index model show that these adaptive particle methods are likely to perform favorably in local cryospheric DA problems with low dimensional parameter spaces. In the ensuing experiments with FSM2, we push the model complexity, number of observations, and the dimensionality of the parameter space considerably to further test the AdaPBS.

4.2 Applying AdaPBS in FSM2 at SnowMIP sites

The large number of assimilated observations in these experiments with thousands of observations per annual DA window represents a challenge for any batch smoothing assimilation algorithm. On average, AdaPBS required 8 iterations per season to address this problem. In principle, such a high number of iterations in AdaPBS would be expected to result in a considerable increase in computational cost compared to the typical number of iterations used by ES-MDA in most cryospheric applications (which in our experience is around four) since a higher number of iterations involve a large number of sequential model realizations. However in this case, the computational cost of the linear algebra in ES-MDA scaled with the number of observations, resulting in a much higher execution time compared with AdaPBS. The site with the largest run times due to the longest simulation, namely WSJ with 6 water years, took 16 min on a single core with AdaPBS compared with 5 h using two cores for ES-MDA. Nonetheless, this conclusion should be interpreted with some caution, as both linear algebra operations and ensemble generation can be parallelized. There may be solutions to circumvent these issues depending on the problem at a hand related to the the parallelization scheme of the computing infrastructure, the hyperparamters (i.e. number of iterations for ES-MDA and diversity threshold for AdaPBS), and further optimizing linear algebra operations (Hennig et al., 2022). Moreover, previous experience suggests that a larger number of iterations of ES-MDA, potentially exceeding N_a=10, may also be necessary to ensure convergence based on earlier work applying the ES-MDA to higher-dimensional problems (Pirk et al., 2024). In a similar vein, the validation metrics of AdaPBS may improve if a larger N_eff threshold is selected, at the cost of increasing the computational expense. More generally, the large number of observations can be assimilated in a more tractable manner via Bayesian filtering (Särkkä and Svensson, 2023) and tempering (Chopin and Papaspiliopoulos, 2020) techniques. Combining adaptive particle methods with these techniques is likely to be a promising future research direction in cryospheric DA.

Figure 7Comparison of the performance of the ESA-MDA (top) and the AdaPBS (bottom) at the Weissfluhjoch (WFJ) ESM-SnowMIP site in the eastern Swiss Alps. The red and green lines with corresponding shading show the prior (“open-loop” in red) and posterior (“updated” in green) mean ±1 standard deviation along with the hourly assimilated snow depth observations (purple dots).

Download

The results of AdaPBS in this more challenging setting with higher-dimensional parameter and observation spaces are similar to those obtained by ES-MDA, as shown in Table 3. It is difficult to conclude confidently which of the two schemes works best, so we conclude tentatively that their performance is similar even in scenarios as complex as the one proposed. The biggest differences are found at SNB, both in terms of RMSE (ES-MDA: 0.21, AdaPBS: 0.44) and CRPS (ES-MDA: 0.15, AdaPBS: 0.29), despite both methods presenting a very similar bias for this site. While the AdaPBS has a tendency to either match or perform slightly worse than ES-MDA in terms of RMSE and CRPS, it generally exhibits a lower bias than ES-MDA particularly for SWA (ES-MDA: −0.16, AdaPBS: −0.06). However, the fact that there are notable improvements with AdaPBS compared to ES-MDA in only one of the sites (SWA), worse performance at one site (SNB), and nearly identical performance at the remaining sites, is clearly not sufficient grounds to claim that there is considerable differences in the performance of these schemes. The results in Table 3 also show that AdaPBS consistently outperforms PBS in this experiment for all metrics across all sites, with the exception of SOD where the low bias differs by only 0.01. Overall, when considering metrics averaged over all sites, use of AdaPBS leads to a 25 %, 41 %, and 60 % reduction in RMSE, CRPS, and bias, respectively, compared to PBS. A likely reason for the large gains from adaptation is that the non-iterative PBS struggles to capture the dense observations, as demonstrated by the average of eight iterations needed in AdaPBS. As a result, PBS often degenerates completely to a suboptimal solution. Although not explicitly evaluated here, a single ES update is expected to also lead to a suboptimal solution due to the combined effects of nonlinearity and the highly informative nature of the observations. This expectation is supported by Experiment 1 as well as previous experiments comparing the ES-MDA to the ES in snow data assimilation (Aalstad et al., 2018; Alonso-González et al., 2022) and other applications (Pirk et al., 2022, 2024; Keetz et al., 2025). It is thus encouraging to see the AdaPBS, a purely particle-based iterative batch smoother, be able to match the performance of the ES-MDA given that AdaPBS provides added flexibility in terms of adapting to the problem at hand through early stopping and is readily extended to non-Gaussian likelihoods (e.g. Smith et al., 2010). Figure 7 demonstrates the generally very good performance of both the ES-MDA and AdaPBS algorithms at WFJ in terms of matching the temporally dense snow depth observations. The same holds for most of the other ESM-SnowMIP sites visualized in the Supplement (see Fig. S1). However, suboptimal posterior predictions are obtained in several water years, especially with AdaPBS but also ES-MDA, at the SNB which turns out to be the most challenging among all the sites (see Fig. S3). Nonetheless, it should be noted that no explicit ERA5 forcing downscaling has been performed other than implicitly by the assimilation itself by inferring the forcing perturbation parameters. The resolution of ERA5 may not be sufficient to capture all the particularities of local observations even after the implicit downscaling at least with the current probabilistic model configuration.

Table 3Snow depth evaluation metrics in the form of the root mean square error (RMSE), mean continuous ranked probability score (CRPS), and mean error (bias), for the Prior, PBS, IES (ES-MDA), and AdaPBS across 6 ESM-SnowMIP sites.

Download Print Version | Download XLSX

4.3 Recommendations for applying and extending the AdaPBS

Broadening the discussion, we provide recommendations informed by the above snow DA experiments and recent applications of AdaPBS to permafrost (Willmes et al., 2025) and glaciers (Yang et al., 2026). In addition to ensemble size N_e, which is foundational for all ensemble-based DA, two key algorithm-specific hyperparameters control AdaPBS: the diversity threshold $\mathrm{0}<\mathit{\tau }<\mathrm{1}$ and the maximum number of iterations ℓ_max≥1 (Sect. 2.6). On the one hand, τ serves as a lower bound for the desired quality of the posterior approximation: an ensemble run that converges successfully will achieve an N_eff of at least τN_e. On the other hand, ℓ_max serves as an upper bound on the computational cost: AdaPBS will run for at most ℓ_max iterations with a total of ℓ_maxN_e model runs, but it may converge and terminate much earlier. This upper bound is only reached if AdaPBS converges in the final possible iteration or, in the worst case, if it fails to converge within the prescribed ℓ_max iterations. There is thus an inherent tradeoff in these hyperparameters: a higher threshold τ≃1 helps improve the particle approximation of the posterior while often requiring more iterations for convergence, demanding a larger ℓ_max with increased computational cost. At the opposite extreme, one could simply set ℓ_max=1 to minimize computational cost, but then AdaPBS is identical to the PBS and has no way to adapt and evolve beyond degeneracy when it occurs. A useful default is ℓ_max=5, mirroring the typical cost of the ES-MDA that uses N_a+1 ensemble iterations with N_a=4 assimilation cycles (Alonso-González et al., 2026), but where AdaPBS imposes fewer assumptions and can stop early upon convergence.

The optimal balance between these two hyperparameters depends on the problem at hand and the available computational budget. Setting the threshold in the range $\mathrm{0.1}<\mathit{\tau }<\mathrm{0.5}$ is generally enough to obtain a good posterior approximation, while the number of iterations depends on the problem. This number may vary considerably with the complexity of the problem and even within a given distributed multi-year simulation. Ensuring convergence across the board may require setting ℓ_max to 10 or more. In practice, it is often entirely acceptable to end up with a low percentage (say ≤10  %) of cases in which the algorithm failed to converge. We emphasize that failure to converge does not imply a complete failure of the algorithm, but rather that the particle approximation is not as diverse as desired for uncertainty quantification via (e.g.) the posterior variance. The ensemble may still provide a reasonable point estimate for the posterior mean and avoid complete degeneracy (N_eff=1). Since the AdaPBS is designed as a direct extension of the PBS, instances of failed convergence will correspond to cases when the PBS would be at least as degenerate.

Experiments 1 and 2 demonstrate the role of problem complexity. For the simpler temperature index model with two parameters, AdaPBS exceeded the desired threshold of τ=0.3 (i.e., N_eff=30 with N_e=100) with N_eff=84 in 5 iterations, providing a strong approximation that nearly converged already at 4 iterations (Fig. 5). In Experiment 2, the intermediate complexity snow model FSM2, a higher dimensional parameter space, and numerous observations conspired to make inference challenging (Snyder et al., 2008; Morzfeld et al., 2017). Following the “needle in a haystack” analogy for inferring the posterior: increased model complexity may result in an oddly shaped needle, more informative data makes the needle smaller, and higher dimensionality makes the haystack bigger. Thus, AdaPBS required, on average, nearly twice as many iterations (8) to converge with the same τ=0.3 as Experiment 1. For both experiments, the number of iterations could have been reduced by lowering the bar on diversity to say τ=0.1, which can suffice in operational settings focused on prediction using point estimates rather than uncertainty quantification (e.g. Oberrauch et al., 2024).

Valuable lessons can also be learned from the two other contemporaneous studies using AdaPBS for cryospheric DA. Willmes et al. (2025) used AdaPBS to assimilate Sentinel-2 and in-situ FSCA data into the detailed multiphysics cryospheric model CryoGrid (Westermann et al., 2023), inferring two parameters to constrain ground surface temperatures around Bayelva near Ny-Ålesund, Svalbard. To configure AdaPBS, the maximum number of iterations was set to ℓ_max=8, and the diversity threshold to τ=0.3. In this two-dimensional parameter space, across all grid cells, convergence was achieved within an average of ℓ=3 iterations, and only rarely did AdaPBS fail to converge. To contend with the relatively high computational cost of this CryoGrid configuration, parallel simulations were run with the ensemble size N_e=20 equal to the number of locally available cores. As such, AdaPBS could optimally leverage parallelization to avoid waiting for new particles to complete while benefiting from adaptation across iterations to avoid degeneracy. For example, the total computational cost of 5 iterations with this configuration is identical to running a non-adaptive PBS with a larger ensemble of N_e=100, as both would require 5 sequential batches of 20 parallel simulations to finish. This approach can be extended to spatially distributed simulations across large domains, where AdaPBS is, like all ensemble-based methods, highly parallelizable in the ensemble dimension. Furthermore, later adaptive iterations will be cheaper to run as some grid cells will have converged and can be skipped. It may also be possible to fold in new cells to replace converged cells, helping to maintain an even computational load. For very large scale applications, it is advantageous to jointly chunk simulations in both the ensemble and spatial dimensions, vectorizing within chunks and parallelizing across. Generally, multiplying the ℓ_max hyperparameter with typical model run times provides a useful worst-case estimate for wall-clock time for allocating computing resources.

Further lessons can be learned from Yang et al. (2026), who used AdaPBS to jointly calibrate 4 parameters controlling surface mass balance and frontal ablation in the coupled PyGEM-OGGM glacier models (Rounce et al., 2023; Maussion et al., 2019). Therein, the evolution of 71 marine-terminating glaciers on Svalbard was inferred for the period 2000–2019 by jointly assimilating glacier-wide mass balance and terminus-positions. Using a configuration with N_e=200, τ=0.1, and ℓ_max=10, AdaPBS converged successfully within 10 iterations for 87 % of the glaciers considered. In terms of performance, the posterior mean from AdaPBS considerably reduced root mean squared errors compared to the prior mean for both the converged and non-converged glaciers. By extending the maximum number of iterations to ℓ_max=100, AdaPBS also managed to converge for most of the remaining glaciers. This showcases both the need for selecting an appropriate AdaPBS configuration and that infrequent non-convergence need not adversely affect overall performance. Moreover, as suggested in Yang et al. (2026), methods such as Bayesian optimization (Hennig et al., 2022) can be used to automatically tune these hyperparameters.

The AdaPBS algorithm that we have presented relies on adaptive importance sampling (Bugallo et al., 2017), specifically AMIS (Cornuet et al., 2012), to extend the PBS (Margulis et al., 2015) by allowing it to evolve beyond collapse when degeneracy strikes. At the same time, adaptation is but one direction for improvement to make particle methods less prone to degeneracy. Likelihood tempering (Neal, 2001; Murphy, 2023), which also powers the ES-MDA (Emerick and Reynolds, 2013; Stordal and Elsheikh, 2015), inflates the observation error covariance and iteratively ensures a smoother transition from prior to posterior. Data tempering (Chopin, 2002; Murphy, 2023) randomly divides the observation vector into mini-batches and assimilates these sequentially, allowing particle methods to assimilate more data without degeneracy (e.g. van Hove et al., 2025, 2026). Particle rejuvenation methods such as jitter (Farchi and Bocquet, 2018) or the resample-move approach using MCMC (Gilks and Berzuini, 2001) also help avoid degeneracy and are frequently used together with tempering for advanced particle methods (Chopin and Papaspiliopoulos, 2020). Hybrid methods are also widely used (van Leeuwen et al., 2019), for example, employing the ES-MDA to design the proposal distribution (Pirk et al., 2022). Crucially, tempering, rejuvenation, and hybridization can all be used together with adaptation to design more robust particle methods for cryospheric DA. Furthermore, as shown theoretically in Elvira et al. (2019), particle methods relying on AMIS, such as AdaPBS, can be readily applied to the auxiliary particle filter, which may further improve operational snow hydrological forecasting systems (Oberrauch et al., 2024). Using adaptive particle methods for spatio-temporal DA (e.g. Mazzolini et al., 2025; Alonso-González et al., 2026), will require further progress on particle localization techniques (Farchi and Bocquet, 2018) for cryospheric DA by building on the efforts of Cluzet et al. (2021). As an extension of the PBS, AdaPBS has the potential to be directly applied to large-scale snow reanalysis efforts (Margulis et al., 2016; Alonso-González et al., 2021; Liu et al., 2021), where ℓ_max provides an upper bound on computational cost.