Heat diffusion in the ocean

The ocean moves heat around, and the simplest way it does so is diffusion: heat flows from warm to cold, smoothing out differences. Drop a patch of warm water into a cold sea and it spreads, flattens, and fades. The rule governing it is the diffusion equation, one of the most important partial differential equations in all of science:

$$\frac{\partial T}{\partial t} = D\,\frac{\partial^2 T}{\partial x^2}.$$

We can't solve a PDE symbolically in general, so we discretize: chop space into a grid, approximate the curvature ∂²T/∂x² by differences between neighbours, and step time forward. Below, a warm block diffuses along a line — set the diffusivity and the elapsed time and watch it spread, computed live on a grid in WebAssembly.

begin
    using PlutoUI, WasmMakie
end

diffusivity = 120 ÷100

time steps elapsed = 400

let
    # finite-difference heat equation on a 1-D grid. The bonds (diffusivity, time) become
    # locals before the loops; the update uses only the local D, so the nested time/space
    # stepping stays cleanly wasm-compilable.
    nx = 61
    D = Float64(diffi) / 100.0
    dt = 0.2
    xs = Vector{Float64}(undef, nx)
    T0 = Vector{Float64}(undef, nx)
    T = Vector{Float64}(undef, nx)
    Tn = Vector{Float64}(undef, nx)
    Tmid = Vector{Float64}(undef, nx)
    for i in 1:nx
        xs[i] = Float64(i)
        v = 0.0
        if i > 20
            if i < 42
                v = 1.0          # a warm block in the middle third
            end
        end
        T0[i] = v
        T[i] = v
        Tmid[i] = v
    end
    half = div(nsteps, 2)
    for step in 1:nsteps
        for i in 2:(nx - 1)
            Tn[i] = T[i] + D * dt * (T[i + 1] - 2.0 * T[i] + T[i - 1])
        end
        Tn[1] = 0.0
        Tn[nx] = 0.0
        for i in 1:nx
            T[i] = Tn[i]
        end
        if step == half
            for i in 1:nx
                Tmid[i] = T[i]
            end
        end
    end

    fig = Figure(size = (620, 340))
    ax = Axis(fig[1, 1])
    lines!(ax, xs, T0)      # the initial warm block
    lines!(ax, xs, Tmid)    # halfway through
    lines!(ax, xs, T)       # fully spread out
    fig
end

0.551

diff_peak = let
    # the final center temperature in a SEPARATE bond-dependent cell so the markdown
    # below interpolates it live (values inside a markdown cell's own `let` bake to the
    # slider defaults).
    nx = 61
    D = Float64(diffi) / 100.0
    dt = 0.2
    T = Vector{Float64}(undef, nx)
    Tn = Vector{Float64}(undef, nx)
    for i in 1:nx
        v = 0.0
        if i > 20
            if i < 42
                v = 1.0
            end
        end
        T[i] = v
    end
    for step in 1:nsteps
        for i in 2:(nx - 1)
            Tn[i] = T[i] + D * dt * (T[i + 1] - 2.0 * T[i] + T[i - 1])
        end
        Tn[1] = 0.0
        Tn[nx] = 0.0
        for i in 1:nx
            T[i] = Tn[i]
        end
    end
    floor(T[31] * 1000.0) / 1000.0
end

The peak (center) temperature has fallen to 0.551 from its starting value of 1.0 as the heat spread outward. Raise the diffusivity or let more time pass and the block flattens faster – diffusion never piles heat up, it only ever smooths it out. Push diffusivity past about 2.5 and the simple scheme goes unstable (the numbers blow up) – a first taste of why choosing the time step matters in PDE solvers.

From a grid to a climate

This little solver is the seed of a real ocean model. Add a flow that carries heat along (advection) and you get the advection-diffusion equation that moves warmth from the tropics toward the poles. Do it on a 2-D or 3-D grid, couple it to the atmosphere, and you have the circulation models behind climate projections.

The technique you just used — replace derivatives with differences between grid points, then march forward in time — is how essentially every PDE in physics and engineering gets solved on a computer. And the stability catch you can trigger with the slider is a real and famous constraint (the CFL condition) that every such solver must respect.

Appendix

A hand-written finite-difference solver for the 1-D diffusion equation, marched forward on a 61-point grid and drawn with WasmMakie — entirely in-browser WebAssembly, no Plots.jl and no PDE library. The method and the stability limit are the standard ones from the MIT advection-diffusion lecture.