Weak-Constrained 4DVar

Strong-constraint 4DVar assumes the forward model is exact. In reality models carry biases, omit physics, and accumulate discretisation error. Weak-constraint 4DVar addresses this by treating model error itself as a control variable Trémolet, 2006: the dynamics are allowed to be wrong by an amount $\boldsymbol{\eta}_t$ at each step, and the analysis estimates those errors alongside the initial state.

Notation follows the strong-constraint note: state $\boldsymbol{u}_t$, background $\boldsymbol{u}_b$, covariances $\mathbf{B}$ / $\mathbf{R}_t$, free-run dynamics $M_t^{\text{free}}$, observation operator $H_t$ — plus the model-error $\boldsymbol{\eta}_t$ with covariance $\mathbf{Q}_t$.

Cost Function¶

The state now evolves under noisy dynamics,

so the control vector expands from $\boldsymbol{u}_0$ to the augmented $(\boldsymbol{u}_0, \boldsymbol{\eta}_1, \ldots, \boldsymbol{\eta}_T)$. The cost gains a third, model-error term:

with the trajectory recovered by recursively stepping (1).

Why the Model-Error Term Matters¶

Letting $\boldsymbol{\eta}_t$ absorb per-step errors improves the fit at every time without distorting $\boldsymbol{u}_0$. It fixes two failure modes of the strong constraint:

Implementation¶

The augmented control $(\boldsymbol{u}_0, \boldsymbol{\eta})$ is a PyTree; the gradient over it comes from autodiff, and a general optimiser (optimistix) handles the minimisation.

import jax
import jax.numpy as jnp
import einx
from jaxtyping import Array, Float

def weak_4dvar_cost(
    u0:  Float[Array, "N"],       # initial state
    eta: Float[Array, "T N"],     # model-error trajectory η₁..η_T
    ub:  Float[Array, "N"],       # background
    ys:  Float[Array, "T M"],     # observations
    M,                            # free-run step:  u -> M^free(u)
    H,                            # observation operator:  u -> H(u)
    apply_B_inv,                  # v -> B⁻¹ v
    apply_R_inv,                  # v -> R⁻¹ v
    apply_Q_inv,                  # v -> Q⁻¹ v
) -> Float[Array, ""]:
    # roll the noisy trajectory:  u_t = M(u_{t-1}) + η_t
    def step(u_prev, eta_t):
        u_t = M(u_prev) + eta_t
        return u_t, u_t
    _, traj = jax.lax.scan(step, u0, eta)                       # (T, N)

    db  = u0 - ub
    Jb  = 0.5 * einx.dot("N, N ->", db, apply_B_inv(db))        # background
    res = ys - jax.vmap(H)(traj)                                # (T, M)
    Jo  = 0.5 * jnp.sum(jax.vmap(lambda r: einx.dot("M, M ->", r, apply_R_inv(r)))(res))
    Jq  = 0.5 * jnp.sum(jax.vmap(lambda e: einx.dot("N, N ->", e, apply_Q_inv(e)))(eta))
    return Jb + Jo + Jq

# gradient over the augmented control (u0, η)
grad_J = jax.grad(weak_4dvar_cost, argnums=(0, 1))

Choosing the model-error covariance $\mathbf{Q}$¶

• Diagonal — white Gaussian model noise, $\mathbf{Q}_t = \operatorname{diag}(\boldsymbol{\sigma}^2)$. Simple, rarely realistic.

• Matérn (spatial) — a structured covariance for spatially correlated model error (FFT-based matvecs).

• Markovian (time-correlated) — an Ornstein–Uhlenbeck $\mathbf{Q}_t$ with correlation timescale $\tau$; error decorrelates over $\sim\tau$ steps, with physical motivation.

Computational Cost¶

The control dimension grows with the window length $T$,

For $T = 100$ and a 2D field with $D_u = 10^4$, that is $\sim 10^6$ control variables instead of 10⁴ — a large jump in the minimisation cost. Two mitigations:

• Low-rank model error — restrict $\boldsymbol{\eta}_t$ to a low-dimensional subspace (POD modes, a learned latent space) via a custom reparameterisation.

• Shorter sub-windows + cycling — run weak-constraint 4DVar over short sub-windows and chain them with a smoother cycle (the operational ECMWF pattern).

Limiting Cases¶

The model-error covariance interpolates a whole spectrum of methods:

That strong-constraint is the tight-$\mathbf{Q}$ limit is the cleanest way to see why weak 4DVar can only help (at the cost of more control variables) when the model genuinely has error.

Posterior¶

The Laplace posterior is Gaussian over the augmented control $(\boldsymbol{u}_0, \boldsymbol{\eta})$, with a block covariance over that joint space. Marginalising to the state trajectory $\{\boldsymbol{u}_t\}_{t=0}^{T}$ means applying the linearised forward at each time and propagating the covariance — implementation-heavy and rarely needed. When full trajectory uncertainty matters, an ensemble smoother is usually preferable.

When Weak-Constraint 4DVar Is the Right Answer¶