diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl old mode 100755 new mode 100644 index 58b3e7849..6570a7652 --- a/examples/meridional_heat_transport_ecco.jl +++ b/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) diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl old mode 100644 new mode 100755 diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 3c7783ae0..8db2d120b 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,7 +1,7 @@ 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, @@ -9,10 +9,10 @@ export frazil_temperature_flux, net_ocean_temperature_flux, sea_ice_ocean_temper 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 diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 6bc6338db..a4c583fd9 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/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 diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 1421eadc5..b7e9886cb 100644 --- a/src/NumericalEarth.jl +++ b/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ᶜᶠᵃ