Add meridional heat transport diagnostic via ocean heat uptake method

examples/meridional_heat_transport_ecco.jl

@@ -15,7 +15,7 @@ Nz = 50
 depth = 5000meters
 z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4)
 
-underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z)
+# underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z)
 underlying_grid = LatitudeLongitudeGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z, longitude = (0, 360), latitude = (-80, 80))
 bottom_height = regrid_bathymetry(underlying_grid;
                                   minimum_depth = 10,
@@ -77,9 +77,10 @@ end
 # And add it as a callback to the simulation.
 add_callback!(simulation, progress, IterationInterval(200))
 
-mht = Field(meridional_heat_transport(esm))
+mht_vT = Field(meridional_heat_transport(esm, MeridionalFluxMethod()))
+mht_OHC = Field(meridional_heat_transport(esm, TendencyMethod()))
 
-ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht);
+ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_vT, mht_OHC);
                                         schedule = TimeInterval(3hours),
                                         filename = "ocean_one_degree_mht",
                                         overwrite_existing = true)
@@ -90,22 +91,26 @@ run!(simulation)
 
 using Oceananigans
 
-mht  = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht"; backend = OnDisk())
+mht_vT  = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_vT"; backend = OnDisk())
+mht_OHC = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht_OHC"; backend = OnDisk())
 
-times = mht.times
+times = mht_vT.times
 Nt = length(times)
 
-grid = mht.grid
-Ny = size(mht.grid, 2)
+grid = mht_vT.grid
+Ny = size(mht_vT.grid, 2)
 
-mht_mean  = deepcopy(mht[1][1, :, 1])
+mht_vT_mean  = deepcopy(mht_vT[1][1, :, 1])
+mht_OHC_mean = deepcopy(mht_OHC[1][1, :, 1])
 
 for iter in 1:Nt
     @info "iteration $iter out of $Nt"
-    mht_mean  +=  mht[iter][1, :, 1]
+    mht_vT_mean  +=  mht_vT[iter][1, :, 1]
+    mht_OHC_mean += mht_OHC[iter][1, :, 1]
 end
 
-@. mht_mean = mht_mean / Nt
+@. mht_OHC_mean = mht_OHC_mean / Nt
+@. mht_vT_mean = mht_vT_mean / Nt
 
 using CairoMakie
 
@@ -114,6 +119,9 @@ ax = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)")
 
 φ = φnodes(grid, Face())
 
-lines!(ax, φ, mht_mean[1:Ny+1]  / 1e15, linewidth=4)
+lines!(ax, φ, mht_vT_mean[1:Ny+1]  / 1e15, linewidth=4, label="via vT")
+lines!(ax, φ, mht_OHC_mean[1:Ny+1] / 1e15, linewidth=4, label="via OHC")
+Legend(fig[2, :], ax, orientation=:horizontal)
+Label(fig[0, :], "Meridional heat transport", fontsize=16, tellwidth=false)
 
 save("mht.png", fig)

src/Diagnostics/Diagnostics.jl

@@ -1,18 +1,18 @@
 module Diagnostics
 
 export MixedLayerDepthField, MixedLayerDepthOperand
-export meridional_heat_transport
+export meridional_heat_transport, TendencyMethod, MeridionalFluxMethod
 export frazil_temperature_flux, net_ocean_temperature_flux, sea_ice_ocean_temperature_flux, atmosphere_ocean_temperature_flux,
        frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
        net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux,
        net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux
 
 using Oceananigans
 using Oceananigans.Architectures: architecture
-using Oceananigans.Models: buoyancy_operation
-using Oceananigans.Grids: new_data, inactive_cell, znode
 using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_regions!
 using Oceananigans.Fields: FieldStatus
+using Oceananigans.Grids: new_data, inactive_cell, znode
+using Oceananigans.Models: buoyancy_operation
 using Oceananigans.Utils: launch!
 using KernelAbstractions: @index, @kernel
 using NumericalEarth.EarthSystemModels: EarthSystemModel

src/Diagnostics/meridional_heat_transport.jl

@@ -1,20 +1,16 @@
+using Oceananigans.TimeSteppers: SplitRungeKuttaTimeStepper, QuasiAdamsBashforth2TimeStepper
 using ..EarthSystemModels: EarthSystemModel, reference_density, heat_capacity
 
+struct MeridionalFluxMethod end
+struct TendencyMethod end
+
 """
-    meridional_heat_transport(esm::EarthSystemModel;
+    meridional_heat_transport(esm::EarthSystemModel, MeridionalFluxMethod();
                               reference_temperature = 0)
 
-Return the meridional heat transport for the coupled `esm::EarthSystemModel` by computing
-the meridional heat flux.
-
-The meridional heat transport is computed via:
-
-```math
-\\mathrm{MHT} ≡ ρᵒᶜ cᵒᶜ ∫  v (T - T_{\\rm ref}) \\, \\mathrm{d}x \\, \\mathrm{d}z
-```
-
-Above, ``T_{\\rm ref}`` is a reference temperature and ``ρᵒᶜ`` and ``cᵒᶜ`` are the
-ocean reference density and specific heat capacity respectively.
+Return the meridional heat transport for the coupled `esm::EarthSystemModel` using either
+two methods: either by directly computing the meridional heat flux or indirectly using
+the total ocean heat content and the ocean heat uptake.
 
 !!! warning "Only works on LatitudeLongitudeGrid"
 
@@ -26,11 +22,93 @@ Arguments
 
 * `esm`: An EarthSystemModel.
 
+* The method for the computation. Available options are: `MeridionalFluxMethod()` (default)
+  and `TendencyMethod()`.
+
+  `MeridionalFluxMethod()` computes the meridional heat transport directly by summing
+  the meridional heat flux; `TendencyMethod()` computes the meridional heat
+  transport indirectly using the ocean heat content.
+
+  1. For `MeridionalFluxMethod()`, the meridional heat transport is computed via:
+
+     ```math
+     \\mathrm{MHT} ≡ ρᵒᶜ cᵒᶜ ∫  v (T - T_{\\rm ref}) \\, \\mathrm{d}x \\, \\mathrm{d}z
+     ```
+
+     Above, ``T_{\\rm ref}`` is a reference temperature and ``ρᵒᶜ`` and ``cᵒᶜ`` are the
+     ocean reference density and specific heat capacity respectively.
+
+  2. For `TendencyMethod()` we have:
+
+     Let ``T`` be three-dimensional (potential) temperature, ``ρᵒᶜ`` the ocean reference
+     density, ``cᵒᶜ`` the specific heat capacity, ``H`` the resting depth, and ``η`` the
+     free-surface elevation.
+
+     The column heat content per unit horizontal area (units of J m⁻²) is:
+
+     ```math
+     ℋ = ρᵒᶜ cᵒᶜ ∫_{-H}^{η} T \\, \\mathrm{d}z
+     ```
+
+     and its time-derivative is a heat flux, ``𝒬 ≡ ∂ℋ/∂t``.
+
+     The vertically-integrated heat budget is
+
+     ```math
+     𝒬 = \\frac{∂ℋ}{∂t} = - \\boldsymbol{∇}_h \\cdot \\boldsymbol{F}_h - 𝒬ᵃᵒ_{\\rm net} + ℛ
+     ```
+
+     where
+
+     * ``\\boldsymbol{F}_h`` is the depth-integrated horizontal heat flux vector (units W m⁻¹),
+       that includes advection and any parameterized lateral/eddy fluxes,
+     * ``𝒬ᵃᵒ_{\\rm net}`` is the [net ocean surface heat flux](@ref NumericalEarth.Diagnostics.net_ocean_heat_flux)
+       (units W m⁻²), and
+     * ``ℛ`` is the residual sources/sinks and non-closed terms (e.g. numerics, unaccounted
+       physics, mass/volume effects).
+
+     The total ocean heat content (OHC) South of latitude ``φ`` is:
+
+     ```math
+     \\mathrm{OHC}_S(φ, t) ≡ ∫_{A(φ)} ℋ \\, \\mathrm{d}A
+     ```
+
+     where ``A(φ)`` is the ocean area South of latitude ``φ``.
+
+     Integrating the vertically-integrated heat budget over ``A(φ)`` and using the divergence
+     theorem we get
+
+     ```math
+     \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S(φ, t) =
+        - ∮_{∂A(φ)} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ
+        - ∫_{A(φ)} 𝒬ᵃᵒ_{\\rm net} \\, \\mathrm{d}A
+        + ∫_{A(φ)} ℛ \\, \\mathrm{d}A
+     ```
+
+     where ``∂A(φ)`` is the boundary of ``A(φ)``, i.e., the circle at latitude ``φ``.
+
+     The northward meridional heat transport across latitude ``φ`` is
+
+     ```math
+     \\mathrm{MHT}(φ, t) ≡ ∮_{\\mathrm{lat}=φ} \\boldsymbol{F}_h \\cdot \\hat{\\boldsymbol{n}} \\, \\mathrm{d}ℓ
+     ```
+
+     with the understanding that ``\\mathrm{MHT} > 0`` implies northward heat transport.
+
+     Ignoring the residual ``ℛ``, the OHC-based diagnostic relation is
+
+     ```math
+     \\begin{align*}
+         \\mathrm{MHT} & = - ∫_{A(φ)} 𝒬ᵃᵒ_{\\rm net} \\, \\mathrm{d}A
+                           - \\frac{\\mathrm{d}}{\\mathrm{d}t} \\, \\mathrm{OHC}_S \\\\
+                       & = - ∫_{A(φ)} \\left( 𝒬ᵃᵒ_{\\rm net} + 𝒬 \\right) \\, \\mathrm{d}A
+     \\end{align*}
+     ```
 
 Keyword Arguments
 =================
 
-* `reference_temperature`: The reference temperature (in ᵒC) used for the calculation; default: 0 ᵒC.
+* `reference_temperature`: The reference temperature (in ᵒC) used for `MeridionalFluxMethod()`; default: 0 ᵒC.
 
   !!! info "Reference temperature"
 
@@ -72,7 +150,8 @@ Integral of BinaryOperation at (Center, Face, Center) over dims (1, 3)
     └── grid: 4×5×2 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×2 halo
 ```
 """
-function meridional_heat_transport(esm::EarthSystemModel; reference_temperature=0)
+function meridional_heat_transport(esm::EarthSystemModel, method=MeridionalFluxMethod();
+                                   reference_temperature=0)
 
     grid = esm.ocean.model.grid
 
@@ -81,15 +160,48 @@ function meridional_heat_transport(esm::EarthSystemModel; reference_temperature=
     grid isa OrthogonalSphericalShellGrid &&
         throw(ArgumentError("meridional_heat_transport diagnostic does not work on OrthogonalSphericalShellGrid at the moment; use LatitudeLongitudeGrid."))
 
-    FT = eltype(esm)
-    reference_temperature = convert(FT, reference_temperature)
+    if method isa MeridionalFluxMethod
+        FT = eltype(esm)
+        reference_temperature = convert(FT, reference_temperature)
+        return meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature)
+    elseif method isa TendencyMethod
+        @warn """
+        Computing the meridional heat transport via the TendencyMethod
+                 does not ensure the heat budget is closed!
+                 If you require a closed heat budget, then use the MeridionalFluxMethod.
+        """
+        return meridional_heat_transport_via_ocean_heat_content(esm)
+    else
+        throw(ArgumentError("Unknown method $(method); choose either MeridionalFluxMethod() or TendencyMethod()."))
+    end
+end
 
+function meridional_heat_transport_via_meridional_heat_flux(esm; reference_temperature)
     ρᵒᶜ = reference_density(esm.ocean)
     cᵒᶜ = heat_capacity(esm.ocean)
-
     T = esm.ocean.model.tracers.T
     v = esm.ocean.model.velocities.v
-
     MHT = Integral(ρᵒᶜ * cᵒᶜ * v * (T - reference_temperature), dims=(1, 3))
     return MHT
 end
+
+function meridional_heat_transport_via_ocean_heat_content(esm)
+    ρᵒᶜ = reference_density(esm.ocean)
+    cᵒᶜ = heat_capacity(esm.ocean)
+    ∂t_T = temperature_tendency(esm.ocean.model.timestepper)
+    𝒬ᵃᵒₙₑₜ = net_ocean_heat_flux(esm) |> Field
+
+    𝒬 = Integral(ρᵒᶜ * cᵒᶜ * ∂t_T, dims=3) |> Field
+    ∫Σ𝒬_dx = Integral(𝒬ᵃᵒₙₑₜ + 𝒬, dims=1) |> Field
+    MHT = CumulativeIntegral(- ∫Σ𝒬_dx, dims=2)
+    return MHT
+end
+
+temperature_tendency(timestepper::SplitRungeKuttaTimeStepper) = timestepper.Gⁿ.T
+
+function temperature_tendency(timestepper::QuasiAdamsBashforth2TimeStepper)
+    Gᵀⁿ = timestepper.Gⁿ.T
+    Gᵀ⁻ = timestepper.G⁻.T
+    χ = timestepper.χ
+    return (3/2 + χ) * Gᵀⁿ - (1/2 + χ) * Gᵀ⁻
+end

src/NumericalEarth.jl

@@ -56,7 +56,7 @@ export
     frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
     net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux,
     net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux,
-    meridional_heat_transport
+    meridional_heat_transport, TendencyMethod, MeridionalFluxMethod
 
 using Oceananigans
 using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ