NWP | WRF-DA | Outer vs. Inner Loops | Incremental approach | Dual-resolution

WRF-DA

In WRF Data Assimilation (WRF-DA), particularly for 3DVar (Three-Dimensional Variational) and 4DVar, the use of "loops" is a strategy to handle the inherent non-linearity of the atmosphere while keeping the math computationally feasible.

Think of it as a "refinement" process: the outer loop handles the complex, messy physics, while the inner loop does the heavy lifting of finding the mathematical best fit.

1. The Outer Loop (Non-linear)

The Outer Loop is where the actual "real-world" physics happen. It manages the non-linear aspects of the model.

  • Function: It computes the Innovation Vector (the difference between the observations and the model's first guess).
  • Physics: It uses the full, non-linear version of the observation operators (functions that convert model variables like temperature/pressure into what a satellite or radar actually sees).
  • Updates: After the inner loop finishes its work, the outer loop takes that result, updates the "First Guess" (background state), and starts a new iteration if necessary.
  • Goal: To ensure the final analysis doesn't drift away from physical reality due to non-linear errors.

2. The Inner Loop (Linearized)

The Inner Loop is where the actual optimization (minimization) occurs. This is the most computationally expensive part, but it is simplified to be fast.

  • Function: It minimizes the Cost Function (\(J\)).
  • Simplification: To stay efficient, the inner loop uses Linearized versions of the model and observation operators (Tangent Linear and Adjoint models).
  • Iteration: It runs many times (often 50–100+ iterations) using an algorithm like the Conjugate Gradient method to find the "valley" or the minimum error point.
  • Resolution: Often, to save time, the inner loop is run at a lower resolution than the outer loop.

3. How They Work Together

The relationship is a nested cycle. You can visualize it as a person trying to find the lowest point in a rugged mountain range:

  1. Outer Loop: You stand at a point, look at your map (observations), and realize you are 500 meters off-track. You simplify the map of your immediate surroundings into a flat, tilted plane.
  2. Inner Loop: You run quickly down that "flat plane" (linearized path) until you reach what you think is the bottom.
  3. Outer Loop: You stop, look at the real, rugged map again, and realize that because the mountain is curvy, your "flat plane" math was slightly off. You update your position and start again.

Key Differences at a Glance

Feature Outer Loop Inner Loop
Model Type Full Non-linear Tangent Linear / Adjoint (Linearized)
Resolution High (Final Output Res) Often Lower (to save power)
Purpose Calculate Innovation (\(y - Hx\)) Minimize Cost Function (\(J\))
Iterations Few (usually 1–3) Many (dozens or hundreds)

Incremental formulation in detail mathematic equation form

To understand the outer and inner loops mathematically, we must look at the incremental formulation of the variational cost function. This approach is used because the original cost function is non-linear and difficult to minimize directly.

1. The Original (Non-linear) Cost Function

The goal of 3DVar/4DVar is to minimize the cost function \(J\) with respect to the state vector \(x\):

\[J(x) = \frac{1}{2}(x - x_b)^T B^{-1} (x - x_b) + \frac{1}{2}(y - H(x))^T R^{-1} (y - H(x))\]

Where:

  • \(x_b\): Background state (first guess).
  • \(B\): Background error covariance matrix.
  • \(y\): Observation vector.
  • \(H\): Non-linear observation operator.
  • \(R\): Observation error covariance matrix.

2. The Outer Loop: Linearization Point

The Outer Loop (indexed by \(n\)) handles the non-linearity by linearizing the problem around a "reference state" \(x^n\).

  1. Innovation Calculation: It computes the difference between actual observations and the model's non-linear projection: \[d^n = y - H(x^n)\]
  2. State Update: The state for the next outer loop is updated using the increment \(\delta x\) found in the inner loop: \[x^{n+1} = x^n + \delta x\]

At the very first outer loop (\(n=0\)), \(x^0 = x_b\).

3. The Inner Loop: Minimizing the Increment

To make the math solvable, we define the Analysis Increment \(\delta x = x - x^n\). We then perform a Taylor expansion on the non-linear operator \(H(x)\) around the reference state \(x^n\):

\[H(x^n + \delta x) \approx H(x^n) + \mathbf{H} \delta x\]

Where \(\mathbf{H}\) is the Tangent Linear (TL) version of the observation operator (the Jacobian matrix).

The Incremental Cost Function

The inner loop minimizes this simplified, quadratic cost function \(J(\delta x)\):

\[J(\delta x) = \frac{1}{2}(\delta x - [x_b - x^n])^T B^{-1} (\delta x - [x_b - x^n]) + \frac{1}{2}(\mathbf{H} \delta x - d^n)^T R^{-1} (\mathbf{H} \delta x - d^n)\]

The Gradient

The inner loop uses the Adjoint (\(\mathbf{H}^T\)) to calculate the gradient of the cost function, which tells the algorithm which way to "roll" toward the minimum:

\[\nabla J(\delta x) = B^{-1}(\delta x - [x_b - x^n]) + \mathbf{H}^T R^{-1} (\mathbf{H} \delta x - d^n)\]

The inner loop iterates (usually using a Conjugate Gradient or Lanczos method) until \(\nabla J(\delta x)\) is sufficiently small.

4. The Summary Cycle

The full mathematical workflow looks like this:

  1. Outer Loop (Start): Define \(x^n\). Calculate the non-linear innovation \(d^n = y - H(x^n)\).
  2. Inner Loop: * Start with \(\delta x = 0\).
    • Solve for \(\delta x\) that minimizes the quadratic \(J(\delta x)\) using the linearized \(\mathbf{H}\) and its adjoint \(\mathbf{H}^T\).
    • This is "inner" because it repeats many times for a single \(d^n\).
  3. Outer Loop (End): Update the state \(x^{n+1} = x^n + \delta x\).
  4. Repeat: If you have multiple outer loops, go back to Step 1 with the new \(x^{n+1}\) to re-linearize and get a more accurate \(H(x)\).

Papers

4D-Var for numerical weather prediction

Schematic representation of the way incremental 4D-Var solves the nonlinear minimisation problem. The full nonlinear cost function (thick blue line) is minimised through the minimisation of successive quadratic cost functions (thin black lines) obtained by linearisation of the cost function around successively more accurate guess fields (Plot courtesy of Sebastien Massart, European Centre for Medium-Range Weather Forecasts).

Schematic representation of the way incremental 4D-Var solves the nonlinear minimisation problem. The full nonlinear cost function (thick blue line) is minimised through the minimisation of successive quadratic cost functions (thin black lines) obtained by linearisation of the cost function around successively more accurate guess fields (Plot courtesy of Sebastien Massart, European Centre for Medium-Range Weather Forecasts).
Schematic representation of the way incremental 4D-Var solves the nonlinear minimisation problem. The full nonlinear cost function (thick blue line) is minimised through the minimisation of successive quadratic cost functions (thin black lines) obtained by linearisation of the cost function around successively more accurate guess fields (Plot courtesy of Sebastien Massart, European Centre for Medium-Range Weather Forecasts).

Dual-resolution

To configure WRF-DA for a multi-resolution or dual-resolution approach, you must leverage the incremental formulation capabilities in the namelist.davar. This allows the inner loop (minimization) to run on a coarser grid, which is then interpolated back to the high-resolution outer loop grid.

Here is the technical breakdown of the configuration.

Representation of the incremental 4D-Var (courtesy of Y. Tremolet).
Representation of the incremental 4D-Var (courtesy of Y. Tremolet).

1. Key Namelist Parameters

In your namelist.davar, the primary switches for resolution control are found in the &wrfvar1, &wrfvar7, and &wrfvar11 blocks.

Parameter Recommended Setting Description
multi_inc 1 or 2 Enables the multi-incremental approach.
low_res_increment .true. Explicitly tells the inner loop to use a reduced grid.
analysis_type "3DVAR" or "4DVAR" Works for both, but most effective for memory-heavy 4DVar.
ntmax 50, 50 Number of inner loop iterations for each outer loop.

2. Setting Up the Resolution Ratio

The WRF-DA system needs to know the relationship between your high-res grid (the outer loop) and your low-res grid (the inner loop). This is handled by defining the grid dimensions and the staggering ratio.

The &wrfvar11 Block

This is where the magic happens for resolution scaling:

1
2
3
4
5
6
7
8
9
10
11
12
&wrfvar11
  cv_options = 3,                      ! Using Control Variable option 3 (usually)
  check_max_iv = .true.,
  ! Resolution scaling:
  ! If outer loop is 3km and you want inner loop at 9km:
  ide_inner = 100,                     ! Number of grid points in X for inner loop
  jde_inner = 100,                     ! Number of grid points in Y for inner loop
  its_inner = 1,
  ite_inner = 100,
  jts_inner = 1,
  jte_inner = 100,
/

Note: The ide_inner and jde_inner should correspond to a grid that matches the domain extent but with a larger \(dx\) and \(dy\). WRF-DA will handle the horizontal interpolation of the increments (\(\delta x\)) back to the high-resolution grid using the internal dual-resolution operators.

3. The Multi-Outer Loop Configuration

To get the benefit of the "sharp and physically consistent" map you mentioned, you should run at least two outer loops. This allows the first inner loop to find a rough solution and the second outer loop to re-linearize the physics at high resolution.

1
2
3
4
&wrfvar7
  max_ext_its = 2,                     ! Number of Outer Loops
  ! Other minimization settings...
/

4. Operational Workflow Requirements

To successfully run in this mode, your input files must be prepared correctly:

  1. Background Error Covariance (\(B\)-matrix): The \(B\)-matrix (generated via gen_be) must be compatible with the inner loop resolution. If your inner loop is at 9km, you should provide a BE file statistics generated for a 9km grid.
  2. Background File (fg): This remains your high-resolution WRF input file (the first guess). The outer loop reads this to compute \(H(x)\).
  3. Observation Files (ob.ascii or bufr): These are independent of resolution, but ensure your thinning parameters in the namelist are appropriate for the inner loop resolution to avoid aliasing.

5. Performance Strategy

Running the inner loop at a lower resolution (e.g., a 3:1 ratio, like 3km outer/9km inner) typically results in:

  • Memory reduction: Approximately 90% reduction in the size of the Hessian matrix and background error calculations for the inner loop.
  • Speed: Minimization converges faster because the "state space" it is searching through is significantly smaller.

In detail with mathematic equation and a concrete example.

To understand the incremental formulation, we have to look at why the "standard" 3DVar equation is a nightmare for computers. The standard form requires calculating the gradient of the non-linear observation operator \(H(x)\) at every single step of the minimization.

The Incremental Formulation bypasses this by shifting the focus from the full state \(x\) to a small change, the increment \(\delta x\).

1. The Mathematical Derivation

The Total State

We define the state \(x\) as the sum of a reference state \(x^n\) (from the outer loop) and a small perturbation \(\delta x\): \[x = x^n + \delta x\]

Linearizing the Observation Operator

We use a Taylor expansion to linearize the non-linear operator \(H\) around the reference state \(x^n\): \[H(x^n + \delta x) \approx H(x^n) + \mathbf{H} \delta x\] Here, \(\mathbf{H}\) is the Tangent Linear (TL) model.

The Incremental Cost Function

Substituting these into the standard 3DVar equation, we get the version used in the Inner Loop:

\[J(\delta x) = \frac{1}{2}(\delta x - \Delta x_b^n)^T \mathbf{B}^{-1} (\delta x - \Delta x_b^n) + \frac{1}{2}(\mathbf{H} \delta x - d^n)^T \mathbf{R}^{-1} (\mathbf{H} \delta x - d^n)\]

Where:

  • \(\Delta x_b^n = x_b - x^n\): The distance between the background (first guess) and the current outer-loop reference state.
  • \(d^n = y - H(x^n)\): The Innovation Vector (the "mismatch" calculated in the outer loop).

2. A Concrete Example: Assimilating Temperature

Imagine you are trying to find the temperature (\(T\)) at a specific grid point.

Initial Conditions (Outer Loop 1)

  • Background (\(x_b\)): Your model says it is \(20^\circ\text{C}\).
  • Observation (\(y\)): A thermometer nearby measures \(23^\circ\text{C}\).
  • Current State (\(x^0\)): We start with the background, so \(x^0 = 20^\circ\text{C}\).
  • Innovation (\(d^0\)): \(y - H(x^0) = 23 - 20 = \mathbf{3^\circ\text{C}}\).

The Inner Loop (Optimization)

The inner loop doesn't know about \(20^\circ\text{C}\) or \(23^\circ\text{C}\). It only sees:

  • The "gap" it needs to close: \(3^\circ\text{C}\).
  • The error weights (\(B\) and \(R\)).

The inner loop iterates to find the best \(\delta x\). Let's say the math determines that since the thermometer is very reliable, we should move \(+2.5^\circ\text{C}\). * \(\delta x = +2.5^\circ\text{C}\)

Updating the State (Outer Loop 2)

Now we update our reference state: \[x^1 = x^0 + \delta x = 20 + 2.5 = \mathbf{22.5^\circ\text{C}}\]

If we run a second outer loop, the system re-calculates the physics (like air density or humidity impacts on temperature) at \(22.5^\circ\text{C}\). It realizes the innovation is now only \(0.5^\circ\text{C}\) and runs the inner loop again to refine that last half-degree.

3. Why this enables Multi-Resolution

The beauty of the equation \(J(\delta x)\) is that \(\delta x\) does not have to be the same size as \(x\).

  • Outer Loop: \(H(x^n)\) is calculated on your 3km grid. It produces a 3km map of innovations (\(d^n\)).
  • Inner Loop: We "downsample" that \(d^n\) to a 9km grid. We solve for \(\delta x\) on this 9km grid. Because a 9km grid has \(1/9\text{th}\) the points of a 3km grid, the matrices \(\mathbf{B}\) and \(\mathbf{H}\) shrink massively.
  • Final Step: Once the inner loop finds the 9km \(\delta x\), WRF-DA uses an interpolator to "up-sample" that \(+2.5^\circ\text{C}\) change back onto the 3km grid and adds it to the original \(x^n\).

4. Summary of the Flow

  1. Outer Loop (High-Res): Does the "Real Physics." It computes \(d^n = y - H(x^n)\).
  2. Inner Loop (Low-Res): Does the "Fast Math." It minimizes \(J(\delta x)\) to find the smoothest, most likely correction.
  3. Update (High-Res): Adds the correction back to the high-res model.

This is exactly how operational centers like ECMWF or NCEP handle global data—they couldn't possibly run the inner loop optimization at the full resolution of the forecast model without running out of memory.

Can we make the background covariance matrix by different resolution mesh and apply into any WRF?

The short answer is yes, but with a very important technical caveat: The resolution of the Background Error Covariance matrix (\(B\)) must match the resolution of the Inner Loop.

In WRF-DA, the \(B\) matrix is not just a pile of data; it defines the spatial correlation scales (how far an observation's influence spreads). If your \(B\) matrix is built for a 9km grid, the math expects a 9km grid.

1. The Resolution "Handshake"

In a multi-resolution setup, the resolutions must align as follows:

  1. Outer Loop: High-resolution (e.g., 3km). Uses fg (first guess) file at 3km.
  2. Inner Loop: Low-resolution (e.g., 9km).
  3. BE File (\(B\) Matrix): Must be 9km.

When you run the inner loop, WRF-DA performs the minimization in "Control Variable" space. If you provide a \(B\) matrix that was generated at 3km but tell the namelist to run the inner loop at 9km, the dimensions of the matrices won't match, and the program will crash with a "dimension mismatch" error during the initialization of the background error.

2. Can you apply a BE file to "any" WRF?

Technically, you can use a BE file from one domain on another, but it is not scientifically recommended unless the resolutions and physical characteristics match.

The "Same Resolution" Rule

If you have a BE file generated for a 10km grid over Europe, you can use it for a 10km grid over North America. The code will run because the grid spacing (\(dx, dy\)) and the number of vertical levels match.

The "Different Resolution" Problem

If you try to use a 30km BE file for a 3km WRF run:

  • The Math: It will fail because the covariance lengths (the distance in grid points) will be completely wrong.
  • The Physics: A 30km BE file assumes "synoptic-scale" errors. A 3km run has "convective-scale" errors. Using the wrong BE file will result in an analysis that is either too "smeared out" or physically inconsistent, leading to model instability (crashes) when the WRF forecast starts.

3. How to make a BE for a different resolution

If you want to use the multi-resolution approach (e.g., 3km Outer / 9km Inner), you must generate the BE file specifically for the Inner Loop resolution using the Gen_BE utility.

The Gen_BE Process:

  1. Run a "Training" Period: Run WRF at your Inner Loop resolution (e.g., 9km) for about 15–30 days, saving forecasts at 12h and 24h intervals (the "NMC Method").
  2. Extract Differences: Use gen_be to calculate the differences between these forecasts.
  3. Generate Matrix: The resulting be.dat will contain the statistics (variances and correlation lengths) specific to that 9km grid.

4. The Namelist Mapping

Once you have your BE file for the lower resolution, you map it in namelist.davar:

1
2
3
4
5
6
7
8
&wrfvar4
  be_filename = 'be.dat_9km',          ! The BE file matched to inner-loop res
/
&wrfvar11
  low_res_increment = .true.,
  ide_inner = 100,                     ! Dimensions of the 9km grid
  jde_inner = 100,
/

Summary Recommendation

To maintain the "sharpness" you mentioned earlier:

  • Generate your BE at the coarse resolution you intend to use for your inner loop.
  • Use at least 2 outer loops. This ensures the high-resolution physics in the outer loop can "correct" any spatial smoothing introduced by the lower-resolution \(B\) matrix in the inner loop.

Does outer loop's dimension match the dimension of B matrix? if not, is it valid?

The short answer is no, the dimensions do not have to match, and yes, it is perfectly valid—provided you are using the incremental formulation correctly.

In fact, the primary reason the incremental formulation exists is to allow the Outer Loop and the B-Matrix (via the Inner Loop) to exist at different resolutions.

1. The Dimensionality Breakdown

To understand why this is valid, look at the "spaces" where the math happens:

  • Analysis Space (Outer Loop): This is your high-resolution grid (e.g., \(1000 \times 1000\) points). The state vector \(x\) lives here.
  • Increment Space (Inner Loop): This is your lower-resolution grid (e.g., \(333 \times 333\) points). The increment \(\delta x\) lives here.
  • B-Matrix: This matrix represents the covariance of errors in the Increment Space. Therefore, the dimensions of \(B\) must match the dimensions of the Inner Loop, not the Outer Loop.

2. How the "Mismatch" is Validated Mathematically

If the dimensions don't match, how does the equation \(x^{n+1} = x^n + \delta x\) work? WRF-DA uses Projection Operators (Interpolation) to bridge the gap.

The math introduces two hidden operators, \(S\) (Slowing/Downsampling) and \(S^{-1}\) (Interpolating/Upsampling):

  1. Downsampling (Outer \(\to\) Inner): The innovation \(d^n = y - H(x^n)\) is calculated at high resolution. Before it enters the inner loop, it is smoothed or interpolated to the coarse grid: \[d_{low}^n = S(d^n)\]

  2. The Minimization: The inner loop solves for \(\delta x\) using the coarse \(B\) matrix. All matrices (\(B, H, R\)) now have matching "coarse" dimensions.

  3. Upsampling (Inner \(\to\) Outer): Once the optimal \(\delta x\) is found on the coarse grid, it is interpolated back to the high-resolution grid: \[x^{n+1} = x^n + S^{-1}(\delta x)\]

3. Why This Is Valid (and Often Better)

Using a \(B\) matrix with smaller dimensions than your outer loop is valid because of the "Information Density" principle:

  • The Physical Reality: Most background errors in meteorology are spatially correlated. A "bump" in pressure error at one point usually means a "bump" in the surrounding 20–50km.
  • The Efficiency: You don't need a \(3\text{km}\) resolution \(B\) matrix to represent a \(20\text{km}\) error correlation. By using a \(9\text{km}\) \(B\) matrix (Inner Loop), you are essentially filtering out "noise" and focusing on the scales where the data assimilation is most effective.

4. Technical Validity Requirements

For WRF-DA to accept this "mismatch" as valid, you must ensure:

  1. Coordinate Consistency: The coarse grid must cover the exact same geographic area as the fine grid.
  2. Namelist Mapping: You must explicitly define the inner loop dimensions in &wrfvar11 (e.g., ide_inner, jde_inner).
  3. The BE File: The be.dat file you provide must have been generated at that coarse resolution. If you provide a \(3\text{km}\) BE file but tell WRF-DA to run a \(9\text{km}\) inner loop, the code will stop because it cannot map a \(3\text{km}\) covariance structure onto a \(9\text{km}\) grid.

Summary Table

Component Dimension Purpose
Outer Loop \(N_{high}\) High-fidelity physics and final output.
Inner Loop \(N_{low}\) Mathematical optimization (minimization).
B-Matrix \(N_{low} \times N_{low}\) Matches the Inner Loop to define error scales.
Link Interpolation Moves data between \(N_{high}\) and \(N_{low}\).

The Verdict: It is not only valid; it is the industry standard for high-resolution forecasting. It prevents the \(B\) matrix from becoming so large that it crashes the computer's RAM.

NWP
#MPAS #WRF #NWP #GFS #FNL #ERA5 #DA #IFS #CMA
NWP | WRF-DA | Outer vs. Inner Loops | Incremental approach | Dual-resolution
https://waipangsze.github.io/2026/04/27/NWP-WRF-DA-Outer-vs-Inner-Loops/
Author
wpsze
Posted on
April 27, 2026
Updated on
April 27, 2026
Licensed under