[Paper Review] Split Newton in Two and Free-Surface Becomes Solvable — Casulli–Zanolli (2012)

$V(\eta) + T\eta = \mathbf{b}$. One line. The first time I saw it, I thought a single Newton step would close the problem. Then I hit the edge of the well — cells whose water depth crosses zero. The Jacobian goes silently to zero there. Same cell, wet at one timestep, dry at the next. That is where every Newton-based free-surface code limps. What Casulli and Zanolli did in 2012 was cut Newton itself in two. The outer iterate rises from below; the inner iterate falls from above; they squeeze the solution from both sides, and monotone convergence is guaranteed.

Drop T to 0.05 and push b past 1.5. The scalar nested Newton tightens from above and from below at the same time.

At T = 0.05 the matrix $T+P-Q$ is nearly singular, and the iteration still stops within five steps. Past b = 1.5 the root sits on the flat tail of $V$, where $\eta > 1$, and the inner loop halts in one step.

Paper info#

• Authors: Vincenzo Casulli, Paola Zanolli
• Journal: Journal of Computational and Applied Mathematics, vol. 236 (2012), pp. 3937–3947
• DOI: 10.1016/j.cam.2012.04.025
• Keywords: mildly nonlinear systems / wetting and drying / free-surface hydrodynamics / Richards equation / confined–unconfined aquifers

Q1. Why does Newton die so often on free surfaces#

The problem looks innocent.

Here $\eta \in \mathbb{R}^N$ is the piezometric head per cell, $T$ is a symmetric positive semi-definite sparse matrix (a discrete diffusion stencil), and $\mathbf{V}$ is a diagonal nonlinear function — the per-cell volume $V_i(\eta_i) = \int_{-\infty}^{\eta_i} a_i(z)\,dz$. The function $a_i$ only needs to have bounded variation (discontinuities allowed). Three canonical free-surface examples:

• wetting/drying: $V(\eta) = \max(0, \eta)$. Dry cells have $V'=0$.
• confined–unconfined aquifer: $V(\eta) = \max(0, \min(c-(-h), \eta-(-h)))$. Two kinks.
• Richards equation in mixed form: $V(\eta)$ depends on soil, smooth but flat over wide ranges.

All three share one trait: $V'(\eta)$ vanishes on whole intervals. Standard Newton uses $J = V'(\eta^{(k)}) + T$. When $V'=0$ and the corresponding row of $T$ couples weakly (or null-space of $T$ is touched), $J$ turns singular. Unless the initial guess is close enough to the truth, the iteration blows up — that is the 30-year-old plague of free-surface codes.

Q2. Slice $\mathbf{V}(\eta)$ in two — Jordan decomposition#

The first trick is splitting the function. Bounded variation of $a_i$ means it is the difference of two non-decreasing functions (Jordan decomposition).

Automatically $V_i$ splits too.

Both $V_1$ and $V_2$ are non-decreasing and bounded. For wetting/drying, $p(\eta)=1\,[\eta \ge 0]$, $q\equiv 0$. For confined–unconfined, $p(\eta)=1\,[\eta \ge -h]$, $q(\eta)=1\,[\eta > c]$. Once split, each piece is monotone — that one sentence is the seed of everything else.

Q3. Two linearizations, one inside, one outside#

We solve

The outer loop linearizes $V_2$ around $\eta^{n-1}$.

with $Q^{n-1} = \mathrm{diag}(q_i(\eta_i^{n-1}))$. Start at $\eta^0 \le \ell$ (the lower fence). This is still nonlinear, so the inner loop linearizes $V_1$:

starting at $\eta^{n,0} \ge u$ (the upper fence). The fences $\ell, u$ are chosen so that $T+P-Q$ is a Stieltjes matrix (symmetric M-matrix). The reward:

• Inner iterates are monotone decreasing from above ($\eta^{n,m+1} \le \eta^{n,m}$).
• Outer iterates are monotone increasing from below ($\eta^{n+1} \ge \eta^n$).
• They squeeze toward the unique solution.

Convergence is proven by inequalities, not floating-point arguments. The famous "Newton converges only near the truth" caveat disappears. The five steps and the single step on the slider above are those inequalities, counted. Once the inner matrix stays Stieltjes, the solve is the same SPD problem as a V-cycle that clears the long waves Jacobi leaves behind.

Q4. A pumped well in 1D — Python implementation#

Section 6 of the paper studies a confined–unconfined aquifer in 2D. I'll shrink it to 1D and code it from scratch. The toy problem is paraboloid bottom/ceiling plus a centered pump — over days, the phreatic region expands and then dry ground emerges. None of the standard CFD toy problems (Sod, Karman, blast) match, so I keep the paper's own setup.

import numpy as np
 
# Per-cell volume: V_i(eta) = a0 * clamp(eta, -h_i, c_i) (shifted by -h_i)
def cell_volume(eta, h, c, a0=0.3):
    return a0 * (np.clip(eta, -h, c) - (-h))
 
# Diagonal Jacobians from the Jordan split
def p_diag(eta, h, a0=0.3):
    # V1 = a0 * max(0, eta - (-h)) -> derivative
    return np.where(eta >= -h, a0, 0.0)
 
def q_diag(eta, c, a0=0.3):
    # V2 = a0 * max(0, eta - c) -> derivative
    return np.where(eta > c, a0, 0.0)
 
def nested_newton_step(eta_prev, h, c, T, dt, phi, tol=1e-6):
    """One time step of the monotone nested Newton solve.
 
    Solves V(eta_new) + T eta_new = V(eta_prev) + dt*phi,
    where T is the SPD matrix from -dt*Laplacian.
    """
    N = len(eta_prev)
    V_prev = cell_volume(eta_prev, h, c)
    rhs0 = V_prev + dt * phi
 
    # Outer start: well below the lower fence (every cell dry)
    eta = np.full(N, -h.max() - 1.0)
    for n_out in range(20):
        Q = np.diag(q_diag(eta, c))
        V2 = np.where(eta > c, 0.3 * (eta - c), 0.0)
        d = rhs0 + V2 - Q @ eta
 
        # Inner start: well above the upper fence (every cell pressurized)
        eta_in = np.full(N, c.max() + 1.0)
        for m_in in range(30):
            P = np.diag(p_diag(eta_in, h))
            V1 = np.where(eta_in > -h, 0.3 * (eta_in - (-h)), 0.0)
            f = P @ eta_in - V1 + d
            A = T + P - Q                     # Stieltjes (symmetric M-matrix)
            eta_in = np.linalg.solve(A, f)    # 1D: dense; 2D: PCG
            r_in = V1 + (T - Q) @ eta_in - d
            if np.max(np.abs(r_in)) < tol:
                break
        eta = eta_in
        # Outer residual
        V_cur = cell_volume(eta, h, c)
        r_out = V_cur + T @ eta - rhs0
        if np.max(np.abs(r_out)) < tol:
            return eta, n_out + 1, m_in + 1
    return eta, 20, m_in + 1
 
# 1D well (paraboloid cross-section)
N = 41
L = 1000.0
x = np.linspace(0, L, N)
h_arr = 10 * np.maximum(0, 1 - ((x - L/2) / (L/2))**2)
c_arr = h_arr.copy()
 
# Diffusion matrix T from dt*kappa/dx^2 * Laplacian (zero-flux BC)
dx = L / (N - 1)
dt = 86400.0     # 1 day
kappa = 1.0      # hydraulic conductivity
alpha = dt * kappa / dx**2
T = -alpha * np.eye(N, k=-1) + 2*alpha*np.eye(N) - alpha*np.eye(N, k=1)
T[0, 0] = alpha; T[0, 1] = -alpha
T[-1, -1] = alpha; T[-1, -2] = -alpha
 
# Initial state: pressurized at the ceiling everywhere
eta = c_arr.copy()
phi = np.zeros(N); phi[N//2] = -8e-4   # central pump (sink)
 
for day in range(10):
    eta, n_out, m_in = nested_newton_step(eta, h_arr, c_arr, T, dt, phi)
    pre = int(np.sum(eta >= c_arr - 1e-6))
    dry = int(np.sum(eta <= -h_arr + 1e-6))
    print(f"day {day+1}: pre={pre} dry={dry} outer={n_out} inner_last={m_in}")

A typical run prints:

day 1: pre=39 dry=0 outer=5 inner_last=5
day 2: pre=37 dry=0 outer=4 inner_last=4
...
day 9: pre=0 dry=8 outer=1 inner_last=5
day 10: pre=0 dry=12 outer=1 inner_last=5

Raise the pump and step +1 day. The pre and dry counts from the printout show up as a patch on the section.

Above 0.0025, ten days open a black dry pocket in the middle. The active-cell count changes every day. The nested loop does not rebuild the Jacobian for that.

Outer averages 3 iterations, inner about 5. That matches Table 1 of the paper qualitatively. The most striking part is the pressurized-to-phreatic transition: both loops frequently halt in a single iteration (Remark 6 in the paper: if $u \le \tilde\eta^{1,1} \le \ell$, you are already at the root).

Q5. Evolve a 2D cross-section interactively#

Read the ten-day log by cell counts only. Day 1: 39 pressurized cells out of 41, 0 dry, outer 5, last inner 5. Day 2: pressurized drops to 37 and outer to 4. Days the printout skips stay unread. Day 9: pressurized 0, dry 8, outer 1, inner 5. Day 10: dry cells rise to 12 and the outer count stays at 1. The pump is $\phi=-8\times 10^{-4}$ on the center cell, the step is 86400, and length 1000 is split into 41 cells. $a_0=0.3$ enters $P$ only on the wet interval.

While the head sits on the ceiling $c$, $q$ is on and $Q$ equals $a_0$. Below $c$, $Q=0$. Below the floor $-h$, $P$ is 0 as well. The outer count falling from 5 to 1 on day 9 is the signal that the run has moved toward Remark 6. The inner count is still 5. Dry cells add rows with $P-Q=0$, and the number of inner linear solves does not shrink. Three outer passes times about five inner passes is still 15 to 25 linear systems, each one Stieltjes.

What the algorithm cannot do — one critical paragraph#

It is not a miracle cure. A few caveats:

1. Cells with $P-Q = 0$. When a cell is dry and the corresponding row of $T$ vanishes, the matrix loses the T2 irreducibility hypothesis. The paper recommends excluding dry cells from the system, but the active-set bookkeeping is non-trivial in practice.
2. It is not fast. Per time step you do outer × inner = 15 to 25 linear solves on average. Plain implicit Euler with CG might be faster when it converges at all. The price you pay is zero divergence events.
3. No direct extension to compressible Euler. The technique is bound to mildly nonlinear systems with diagonal $V$ of bounded variation. Riemann-solver nonlinearities do not factor this way.
4. Inner solver still matters. Each inner system is SPD, but how you precondition it is a separate art. The paper recommends PCG; multigrid is usually a step up for large grids.

Reproducibility scorecard#

• Math (§2–§5): self-contained. With pen and paper you can rebuild every inequality. A.
• Algorithm boxes (1 and 2): pseudocode is complete. A.
• Numerical example (§6): grid spacing, time step, all coefficients spelled out. My 30-line Python reproduces the paper qualitatively. A.
• Practical landing: OpenFOAM lacks a first-class wetting/drying module, but meshFlow/overFlowSolver variants could host it. Ocean/estuary codes (SCHISM, TELEMAC) already use related ideas.

When to drop the dry cells#

Rows with $P-Q=0$ break the T2 hypothesis. The paper says to remove dry cells. Day 10 has 12 dry cells out of 41. Remove them and the active set changes every day. Leave them and zero rows remain. Nested Newton keeps the inequality convergence either way. What changes is the size of the linear system.

Finite termination in Brugnano and Casulli (2008/2009) is a piecewise-linear statement. Smooth $a(\eta)$, as in Richards, belongs with Casulli and Zanolli (2010). A free surface made by waves lives in the depth that opens a Kelvin wedge. Here the free surface is a dry cell.