ekalmX — Covariance Localization & Ensemble Inflation¶

Subject: Gap analysis of covariance localization methods and ensemble inflation strategies for the ensemble Kalman filter stack.

Date: 2026-04-03

Part A — Covariance Localization¶

A.1 Why Localize¶

With a finite ensemble of size \(N_e\), the sample covariance \(P = \frac{1}{N_e - 1} X'^T X'\) has rank at most \(N_e - 1\), while the true covariance may be full rank. The deficit manifests as spurious long-range correlations — the sample covariance assigns non-zero correlation between physically distant variables purely due to sampling noise. These spurious correlations cause the Kalman gain to draw information from distant, uninformative observations, leading to filter instability and eventual divergence.

Localization suppresses these artifacts by tapering covariance entries as a function of physical distance, zeroing out correlations beyond a characteristic radius.

A.2 Localization Framework¶

Two complementary approaches exist:

B-localization (covariance tapering). The localized covariance is formed via the Schur (Hadamard / element-wise) product:

\[P_{\text{loc}} = \rho \circ P\]

where \(\rho_{ij} = \rho(d_{ij} / r)\) is a positive-definite taper function of the distance \(d_{ij}\) between grid points \(i\) and \(j\), with localization radius \(r\). The Schur product of two PSD matrices is PSD (Schur product theorem), so \(P_{\text{loc}}\) remains a valid covariance. Applied to the gain:

\[K_{\text{loc}} = (\rho \circ C^{xH})\,(\rho \circ C^{HH} + R)^{-1}\]

R-localization (observation-space). Instead of tapering the covariance, inflate the observation noise locally:

\[R_{\text{loc}} = \text{diag}(1 / \rho_i) \circ R\]

where \(\rho_i\) is the taper weight for observation \(i\) relative to the analysis grid point. Observations far from the analysis point receive inflated noise and contribute less to the update. This is the approach used by the LETKF — each grid point sees only local observations with distance-weighted noise inflation.

A.2.1 Geographic localization¶

For real-world atmospheric and oceanic applications the distance metric is geographic distance (great-circle or projected), not grid index. filterax ships a GeoLocalizer for this case (Decision D14):

class GeoLocalizer(AbstractLocalizer):
    """Gaspari-Cohn (or other taper) over geographic distance."""

    coords: Float[Array, "N_x 2"]            # lon/lat of state grid points (static)
    frame: "LocalFrame" = eqx.field(static=True)  # CRS + projection origin, built with pyproj
    radius_km: float = eqx.field(static=True)
    taper: str = eqx.field(static=True, default="gaspari_cohn")

LocalFrame is built once with pyproj (outside jax.jit) and held as static metadata. Pairwise distances are precomputed at construction. The JIT path is just taper(distances / radius_km). pyproj is therefore a soft optional dependency — the localizer itself remains JAX-pure.

GeoLocalizer pairs with patcher-based LETKF (see A.2.2) for large-domain assimilation: the patcher slices coords per patch; GeoLocalizer applies the taper inside each patch.

A.2.2 Patcher-localized LETKF¶

Continental- or basin-scale domains do not fit a full ensemble × full state in device memory. filterax addresses this with patch decomposition (Decision D16):

Wave 2 L0 primitives: create_patches, assign_obs_to_patches, blend_patches.
Wave 4 L2 model: LocalEnKF / patcher-LETKF — runs an independent LETKF analysis per patch, blends overlapping regions.
Wave 4 extension point: AbstractPatcher protocol so future geostack patchers (geotoolz.patch) can be plugged in without changing filterax-side code.

class AbstractPatcher(eqx.Module):
    """Decompose a spatial domain into overlapping patches and stitch back."""

    @abc.abstractmethod
    def patches(self, state) -> Iterable["Patch"]: ...

    @abc.abstractmethod
    def stitch(self, patches: Iterable["Patch"], shape) -> Array: ...

filterax ships an in-house implementation in Wave 2; it migrates to consuming geotoolz.patch (still incubating in geostack) once that surface stabilizes.

A.3 Gap Catalog¶

Gap 1: Gaspari-Cohn Taper¶

Domain: Any gridded or unstructured spatial domain with a distance metric.

Math (Gaspari & Cohn 1999): A 5th-order piecewise polynomial on \(z = d/r \in [0, 2]\):

\[\rho(z) = \begin{cases} -\frac{1}{4}z^5 + \frac{1}{2}z^4 + \frac{5}{8}z^3 - \frac{5}{3}z^2 + 1, & 0 \leq z \leq 1 \\ \frac{1}{12}z^5 - \frac{1}{2}z^4 + \frac{5}{8}z^3 + \frac{5}{3}z^2 - 5z + 4 - \frac{2}{3z}, & 1 < z \leq 2 \\ 0, & z > 2 \end{cases}\]

Properties: Compactly supported (zero beyond \(2r\)). \(C^2\) smooth at all transition points (\(z = 0, 1, 2\)). Positive definite — the localized covariance is guaranteed PSD. The gold standard in operational NWP and ocean DA.

def gaspari_cohn(
    distances: Float[Array, "..."],
    radius: float,
) -> Float[Array, "..."]:
    """Gaspari-Cohn 5th-order piecewise polynomial taper.

    Compactly supported: zero beyond 2 x radius.
    Smooth (C2) at all transition points.
    """

Ref: Gaspari, G. & Cohn, S. E. (1999). Construction of correlation functions in two and three dimensions. Quart. J. Roy. Meteor. Soc., 125, 723-757.

Impl ref: filterX/api/primitives.md §localization

Gap 2: Gaussian Taper¶

Domain: General. Useful when compact support is not required.

Math:

\[\rho(d) = \exp\!\left(-\frac{d^2}{2r^2}\right)\]

Properties: Infinitely differentiable (\(C^\infty\)). Positive definite. Not compactly supported — decays exponentially but never reaches zero. Simpler than Gaspari-Cohn but retains weak long-range correlations.

def gaussian_taper(
    distances: Float[Array, "..."],
    radius: float,
) -> Float[Array, "..."]:
    """Gaussian taper: exp(-d^2 / (2 radius^2))

    Not compactly supported — decays but never reaches zero.
    """

Ref: Standard. See also Gaspari & Cohn (1999) for comparison.

Impl ref: filterX/api/primitives.md §localization

Gap 3: Hard Cutoff¶

Domain: Testing and prototyping only.

Math:

\[\rho(d) = \mathbf{1}\{d \leq r\}\]

Properties: Discontinuous at \(d = r\). Not positive definite in general (the resulting matrix may lose PSD-ness). Can cause filter instability due to the sharp boundary. Useful only as a baseline or for debugging.

def hard_cutoff(
    distances: Float[Array, "..."],
    radius: float,
) -> Float[Array, "..."]:
    """Binary cutoff: 1 if d <= radius, 0 otherwise.

    Discontinuous — can cause filter instability. Testing only.
    """

Impl ref: filterX/api/primitives.md §localization

Gap 4: SOAR — Second-Order Auto-Regressive¶

Domain: General. Common in meteorological background error modelling.

Math:

\[\rho(d) = \left(1 + \frac{d}{r}\right) \exp\!\left(-\frac{d}{r}\right)\]

Properties: \(C^1\) smooth (first derivative continuous, second is not). Positive definite. Approximates compact support — decays to \(< 0.01\) by \(d \approx 5r\). Often used as a correlation model in background error covariance (e.g., Met Office VAR). Simpler than Gaspari-Cohn but less smooth.

def soar_taper(
    distances: Float[Array, "..."],
    radius: float,
) -> Float[Array, "..."]:
    """Second-Order Auto-Regressive taper.

    rho = (1 + d/r) exp(-d/r)

    C1 smooth, positive definite, approximate compact support.
    """

Ref: Thiebaux, H. J. & Pedder, M. A. (1987). Spatial Objective Analysis. Academic Press.

Gap 5: Adaptive Localization¶

Domain: General. Eliminates manual tuning of the localization radius.

Math (Anderson 2007, 2012): The optimal localization radius for each state-observation pair is estimated from innovation statistics. For each variable \(x_i\) and observation \(y_j\), the sample correlation \(r_{ij} = C^{xH}_{ij} / (\sigma_{x_i} \sigma_{y_j})\) is compared to its sampling uncertainty:

\[\text{se}(r_{ij}) \approx \frac{1 - r_{ij}^2}{\sqrt{N_e - 2}}\]

If \(|r_{ij}| < \text{se}(r_{ij})\), the correlation is not statistically significant and is set to zero. The effective localization radius adapts per-variable — strongly correlated pairs retain large radii, weakly correlated pairs are suppressed aggressively.

Properties: Data-driven, no manual radius selection. Requires multiple cycles to estimate correlations reliably. More complex to implement (per-pair significance testing). Can be combined with a Gaspari-Cohn base taper.

def adaptive_localization(
    particles: Float[Array, "N_e N_x"],
    obs_particles: Float[Array, "N_e N_y"],
    significance_level: float = 0.05,
) -> Float[Array, "N_x N_y"]:
    """Adaptive localization (Anderson 2007, 2012).

    Estimates per-variable localization weights from the statistical
    significance of ensemble correlations. Correlations below the
    sampling noise floor are zeroed out.
    """

Ref: Anderson, J. L. (2007). Exploring the need for localization in ensemble data assimilation using a hierarchical ensemble filter. Physica D, 230, 99-111. Anderson, J. L. (2012). Localization and sampling error correction in ensemble Kalman filter data assimilation. Mon. Wea. Rev., 140, 2359-2371.

Gap 6: Domain Localization (for LETKF)¶

Domain: Gridded domains with the Local Ensemble Transform Kalman Filter.

Math: Domain localization is not a taper function applied to the covariance. Instead, for each analysis grid point \(i\):

Select all observations within radius \(r\) of grid point \(i\): \(\mathcal{O}_i = \{j : d(i, j) \leq r\}\).
Apply distance-dependent observation-noise inflation to the selected observations: \(R_{\text{loc},jj} = R_{jj} / \rho(d_{ij})\).
Run a local ETKF analysis using only \(\mathcal{O}_i\) — a small \(N_e \times |\mathcal{O}_i|\) problem.

The localization is implicit in the observation selection and noise weighting. No global taper matrix is ever formed, making this approach memory-efficient for large state spaces.

Properties: Embarrassingly parallel over grid points. Natural fit for LETKF. The local analyses are independent and can run on separate cores/devices. Observation selection uses a spatial index (KD-tree, ball tree) for efficiency.

def select_local_obs(
    grid_point: Float[Array, "D"],
    obs_coords: Float[Array, "N_y D"],
    radius: float,
) -> Int[Array, "N_local"]:
    """Select observation indices within localization radius."""

def inflate_local_obs_noise(
    obs_noise_local: Float[Array, "N_local N_local"],
    distances: Float[Array, "N_local"],
    taper_fn: Callable = gaspari_cohn,
    radius: float = 1.0,
) -> Float[Array, "N_local N_local"]:
    """R-localization: inflate observation noise by 1/rho(d).

    R_loc = diag(1 / rho(d)) @ R_local

    Observations near the grid point keep their original noise;
    distant observations receive inflated noise and contribute less.
    """

Ref: Hunt, B. R., Kostelich, E. J., & Szunyogh, I. (2007). Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter. Physica D, 230, 112-126.

Part B — Ensemble Inflation¶

B.1 Why Inflate¶

Ensemble Kalman filters systematically underestimate posterior uncertainty due to three compounding effects:

Finite ensemble size. With \(N_e\) members, the sample covariance underestimates the true variance — the ensemble cannot represent all directions of uncertainty.
Imperfect model. Model error introduces unrepresented uncertainty that the ensemble does not account for.
Localization side effects. Covariance tapering reduces the effective rank of the covariance further, compounding the underestimation.

Without inflation, the ensemble spread collapses over repeated assimilation cycles. The filter becomes overconfident, assigns excessive weight to the forecast, and rejects observations — a failure mode called filter divergence.

B.2 Gap Catalog¶

Gap 7: Multiplicative Inflation¶

Domain: General. The simplest and most common inflation method.

Math:

\[x^{(j)}_{\text{inflated}} = \bar{x} + \lambda (x^{(j)} - \bar{x}), \qquad \lambda > 1\]

Equivalently, \(X'_{\text{inflated}} = \lambda X'\). The inflated covariance is \(P_{\text{inflated}} = \lambda^2 P\).

Properties: One scalar parameter. Applied before the analysis step (prior inflation) or after (posterior inflation). Typical values: \(\lambda \in [1.01, 1.10]\). Does not account for spatially varying model error.

def inflate_multiplicative(
    particles: Float[Array, "N_e N_x"],
    factor: float,
) -> Float[Array, "N_e N_x"]:
    """Multiplicative inflation: X'_inflated = factor * X'

    Scales anomalies around the ensemble mean.
    factor > 1 increases spread, factor = 1 is identity.
    """

Ref: Anderson, J. L. & Anderson, S. L. (1999). A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. Mon. Wea. Rev., 127, 2741-2758.

Impl ref: filterX/api/primitives.md §inflation

Gap 8: RTPS — Relaxation to Prior Spread¶

Domain: General. Widely used in operational NWP (e.g., NCEP GFS).

Math (Whitaker & Hamill 2012): Relaxation to Prior Spread rescales the analysis anomalies so that the per-variable standard deviation interpolates between the analysis and forecast values:

\[\sigma_{\text{relaxed},i} = (1 - \alpha)\,\sigma^a_i + \alpha\,\sigma^f_i\]

\[x^{(j)}_{\text{relaxed}} = \bar{x}^a + \frac{\sigma_{\text{relaxed},i}}{\sigma^a_i} \cdot (x^{(j)}_a - \bar{x}^a)\]

where \(\sigma^a_i\) and \(\sigma^f_i\) are the per-variable analysis and forecast ensemble standard deviations. With \(\alpha = 0\) the analysis is unchanged; with \(\alpha = 1\) the forecast spread is fully restored.

Properties: Per-variable inflation — spatially adaptive. Applied after the analysis step. Preserves the analysis mean. One scalar parameter \(\alpha\), but the effective inflation factor varies across variables.

def inflate_rtps(
    analysis_particles: Float[Array, "N_e N_x"],
    forecast_particles: Float[Array, "N_e N_x"],
    alpha: float,
) -> Float[Array, "N_e N_x"]:
    """Relaxation to Prior Spread (Whitaker & Hamill 2012).

    sigma_relaxed = (1 - alpha) sigma_analysis + alpha sigma_forecast

    Parameters
    ----------
    alpha : relaxation coefficient in [0, 1].
        0 = pure analysis, 1 = restore full forecast spread.
    """

Ref: Whitaker, J. S. & Hamill, T. M. (2012). Evaluating methods to account for system errors in ensemble data assimilation. Mon. Wea. Rev., 140, 3078-3089.

Impl ref: filterX/api/primitives.md §inflation

Gap 9: RTPP — Relaxation to Prior Perturbations¶

Domain: General. Commonly paired with ETKF/LETKF.

Math (Zhang et al. 2004): Directly interpolates the analysis perturbations toward the forecast perturbations:

\[X'_{\text{relaxed}} = (1 - \alpha)\,X'_a + \alpha\,X'_f\]

Unlike RTPS which operates on per-variable spread (a scalar per variable), RTPP operates on the full perturbation matrix and thus preserves inter-variable correlation structure from the forecast.

Properties: One scalar parameter. Simpler than RTPS — no per-variable standard deviation computation. Modifies correlation structure (RTPS does not). Applied after the analysis step.

def inflate_rtpp(
    analysis_particles: Float[Array, "N_e N_x"],
    forecast_particles: Float[Array, "N_e N_x"],
    alpha: float,
) -> Float[Array, "N_e N_x"]:
    """Relaxation to Prior Perturbations (Zhang et al. 2004).

    X'_relaxed = (1 - alpha) X'_analysis + alpha X'_forecast

    Directly interpolates perturbations (not spread).
    """

Ref: Zhang, F., Snyder, C., & Sun, J. (2004). Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter. Mon. Wea. Rev., 132, 1238-1253.

Impl ref: filterX/api/primitives.md §inflation

Gap 10: Additive Inflation¶

Domain: General. Explicitly models model error.

Math:

\[X'_{\text{inflated}} = X' + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, Q_{\text{add}})\]

Random noise drawn from a prescribed model-error covariance \(Q_{\text{add}}\) is added to each ensemble member's perturbations. This directly injects additional spread in directions specified by \(Q_{\text{add}}\).

Properties: Requires specifying \(Q_{\text{add}}\) — a non-trivial modelling choice. Can target specific variables or spatial scales. Common sources for \(Q_{\text{add}}\): climatological variability, lagged forecast differences, or stochastic physics tendencies.

def inflate_additive(
    key: PRNGKey,
    particles: Float[Array, "N_e N_x"],
    noise_cov: AbstractLinearOperator,
) -> Float[Array, "N_e N_x"]:
    """Additive inflation: X' <- X' + epsilon, epsilon ~ N(0, Q_add).

    Adds stochastic model-error perturbations to the ensemble.

    Parameters
    ----------
    noise_cov : model-error covariance Q_add.
        Sampling delegates to gaussx.
    """

Ref: Hamill, T. M. & Whitaker, J. S. (2005). Accounting for the error due to unresolved scales in ensemble data assimilation: A comparison of different approaches. Mon. Wea. Rev., 133, 3132-3147.

Gap 11: Adaptive Inflation¶

Domain: General. Eliminates manual tuning of the inflation factor.

Math (Anderson 2007, 2009): The optimal inflation factor is estimated from the consistency of innovations. The expected innovation variance is:

\[\text{var}(v) = C^{HH} + R\]

If the observed innovation variance exceeds this (innovations are too large), the ensemble is under-dispersive and inflation should increase. Anderson (2009) frames this as a Bayesian estimation problem: the inflation factor \(\lambda\) is treated as a random variable with a prior \(\mathcal{N}(\lambda_{\text{prior}}, \sigma^2_\lambda)\), updated each cycle via:

\[\lambda_{\text{post}} = \frac{\sigma^2_\lambda \cdot \hat{\lambda}_{\text{obs}} + \sigma^2_{\text{obs}} \cdot \lambda_{\text{prior}}}{\sigma^2_\lambda + \sigma^2_{\text{obs}}}\]

where \(\hat{\lambda}_{\text{obs}}\) is the inflation implied by the observed innovation variance and \(\sigma^2_{\text{obs}}\) is its uncertainty. The posterior becomes the prior for the next cycle.

Properties: Self-tuning — no manual \(\lambda\) selection. Can be applied per-variable (spatially varying inflation). Requires multiple cycles to converge. Needs careful initialization of the prior.

def inflate_adaptive(
    particles: Float[Array, "N_e N_x"],
    obs: Float[Array, "N_y"],
    obs_op: Callable,
    obs_noise: AbstractLinearOperator,
    prior_mean: float = 1.0,
    prior_var: float = 0.01,
    min_factor: float = 1.0,
    max_factor: float = 1.2,
) -> tuple[Float[Array, "N_e N_x"], float, float]:
    """Adaptive multiplicative inflation (Anderson 2007, 2009).

    Estimates optimal inflation from innovation consistency:
    var(innovation) should equal C^{HH} + R.

    Returns inflated particles and updated (mean, variance) of the
    inflation factor for use as prior in the next cycle.
    """

Ref: Anderson, J. L. (2007). An adaptive covariance inflation error correction algorithm for ensemble filters. Mon. Wea. Rev., 135, 1286-1299. Anderson, J. L. (2009). Spatially and temporally varying adaptive covariance inflation for ensemble filters. Tellus A, 61, 72-83.

Impl ref: filterX/api/primitives.md §inflation

Gap 12: Ledoit-Wolf Shrinkage¶

Domain: General. Optimal covariance regularization from random matrix theory.

Math (Ledoit & Wolf 2004): The shrinkage estimator replaces the sample covariance with a convex combination toward a structured target:

\[P_{\text{shrunk}} = (1 - \lambda^*)\,P_{\text{sample}} + \lambda^*\,\mu I\]

where \(\mu = \text{tr}(P_{\text{sample}}) / N_x\) (the average eigenvalue) and \(\lambda^*\) is the oracle-approximating shrinkage intensity:

\[\lambda^* = \frac{\sum_{i,j} \text{Var}(P_{ij})}{\sum_{i,j} (P_{ij} - \mu \delta_{ij})^2}\]

minimizing \(\mathbb{E}\|P_{\text{shrunk}} - P_{\text{true}}\|_F^2\).

Properties: Not inflation per se, but serves a similar purpose — regularizes the sample covariance by pulling eigenvalues toward their mean. Guaranteed PSD. Analytically optimal (no tuning parameter). Particularly effective when \(N_e \ll N_x\). Can be applied as a covariance correction step before or instead of multiplicative inflation.

def ledoit_wolf_shrinkage(
    particles: Float[Array, "N_e N_x"],
) -> tuple[Float[Array, "N_x N_x"], float]:
    """Ledoit-Wolf optimal shrinkage of the sample covariance.

    P_shrunk = (1 - lambda*) P_sample + lambda* mu I

    Returns the shrunk covariance and the optimal shrinkage intensity.
    """

Ref: Ledoit, O. & Wolf, M. (2004). A well-conditioned estimator for large-dimensional covariance matrices. J. Multivariate Anal., 88, 365-411.

Section C — Shared Infrastructure & References¶

C.1 Shared Infrastructure¶

All localization and inflation primitives share common infrastructure and compose with the ensemble Kalman filter pipeline:

Component	Source	Notes
Distance computation	`filterx.primitives.distances`	Pairwise distances for taper evaluation
Taper function dispatch	`filterx.primitives.localization`	`gaspari_cohn`, `gaussian_taper`, `hard_cutoff`, `soar_taper`
Hadamard product	`filterx.primitives.localize`	\(P_{\text{loc}} = \rho \circ P\) — works on dense and structured operators
Ensemble statistics	`filterx.primitives.ensemble_mean`, `ensemble_anomalies`	Shared by inflation methods
Low-rank covariance	`gaussx.operators.LowRankUpdate`	Ensemble covariance as rank-\((N_e-1)\) operator
Structured solve	`gaussx.primitives.solve`	Woodbury solve for localized/inflated gain
Observation noise	`gaussx.operators`	`DiagonalLinearOperator` for R, R-localization via `diag(1/rho)`
Spatial indexing	`filterx.primitives.select_local_obs`	KD-tree / ball-tree for domain localization
PRNG	`jax.random`	Deterministic sampling for additive inflation

C.2 Composition Patterns¶

Localization and inflation are typically composed in a specific order within the assimilation cycle:

Prior inflation (multiplicative or adaptive) — before the analysis step.
Localized gain computation (B-localization or domain localization).
Analysis update — apply the localized gain.
Posterior inflation (RTPS, RTPP, or additive) — after the analysis step.

C.3 References¶

Gaspari, G. & Cohn, S. E. (1999). Construction of correlation functions in two and three dimensions. Quart. J. Roy. Meteor. Soc., 125, 723-757.
Anderson, J. L. & Anderson, S. L. (1999). A Monte Carlo implementation of the nonlinear filtering problem. Mon. Wea. Rev., 127, 2741-2758.
Zhang, F., Snyder, C., & Sun, J. (2004). Impacts of initial estimate and observation availability on convective-scale data assimilation. Mon. Wea. Rev., 132, 1238-1253.
Hamill, T. M. & Whitaker, J. S. (2005). Accounting for the error due to unresolved scales in ensemble data assimilation. Mon. Wea. Rev., 133, 3132-3147.
Anderson, J. L. (2007). An adaptive covariance inflation error correction algorithm for ensemble filters. Mon. Wea. Rev., 135, 1286-1299.
Hunt, B. R., Kostelich, E. J., & Szunyogh, I. (2007). Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter. Physica D, 230, 112-126.
Anderson, J. L. (2009). Spatially and temporally varying adaptive covariance inflation for ensemble filters. Tellus A, 61, 72-83.
Whitaker, J. S. & Hamill, T. M. (2012). Evaluating methods to account for system errors in ensemble data assimilation. Mon. Wea. Rev., 140, 3078-3089.
Anderson, J. L. (2012). Localization and sampling error correction in ensemble Kalman filter data assimilation. Mon. Wea. Rev., 140, 2359-2371.
Ledoit, O. & Wolf, M. (2004). A well-conditioned estimator for large-dimensional covariance matrices. J. Multivariate Anal., 88, 365-411.
Thiebaux, H. J. & Pedder, M. A. (1987). Spatial Objective Analysis. Academic Press.
Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res., 99, 10143-10162.
Vetra-Carvalho, S., et al. (2018). State-of-the-art stochastic data assimilation methods for high-dimensional non-Gaussian problems. Tellus A, 70, 1445364.