simone-silvestri · Nov 6, 2026 · Sep 28, 2026 · Sep 28, 2026 · Sep 28, 2026 · Sep 28, 2026
diff --git a/docs/src/models/turbulence_closures.md b/docs/src/models/turbulence_closures.md
@@ -219,6 +219,7 @@ CATKEVerticalDiffusivity{VerticallyImplicitTimeDiscretization}
 ├── negative_tke_time_scale: 60.0
 ├── minimum_convective_buoyancy_flux: 1.0e-11
 ├── tke_time_step: nothing
+├── penetrative_radiation: nothing
 ├── mixing_length: TKEBasedVerticalDiffusivities.CATKEMixingLength
 │   ├── Cˢ:   1.131
 │   ├── Cᵇ:   0.967

diff --git a/src/Forcings/Forcings.jl b/src/Forcings/Forcings.jl
@@ -6,7 +6,11 @@ using Adapt: Adapt, adapt
 using DocStringExtensions: TYPEDSIGNATURES
 
 using Oceananigans: Oceananigans
+using Oceananigans.BuoyancyFormulations: BuoyancyForce, BuoyancyTracerFormulation, SeawaterBuoyancy,
+                                         get_temperature_and_salinity, thermal_expansionᶜᶜᶠ
 using Oceananigans.Fields: field, location
+using Oceananigans.Grids: Center, Face, inactive_cell, znode
+using Oceananigans.Operators: ∂zᶜᶜᶜ
 using Oceananigans.OutputReaders: FlavorOfFTS
 using Oceananigans.Units: Time
 using Oceananigans.Utils: named_tuple
@@ -21,5 +25,6 @@ include("forcing.jl")
 include("model_forcing.jl")
 include("compute_forcing.jl")
 include("tidal_forcing.jl")
+include("penetrative_radiation.jl")
 
 end # module
diff --git a/src/Forcings/penetrative_radiation.jl b/src/Forcings/penetrative_radiation.jl
@@ -0,0 +1,80 @@
+#####
+##### Interface for radiation absorbed below the surface
+#####
+
+"""
+    AbstractRadiativeForcing
+
+Forcing by a radiative flux that enters through the surface and is absorbed below it, heating the temperature, or the
+buoyancy of a `BuoyancyTracer`, with the convergence of the flux. Subtypes extend `surface_radiative_flux`,
+`transmitted_fraction`, `transmitted_fraction_derivative` and `transmitted_thickness`.
+"""
+abstract type AbstractRadiativeForcing end
+
+"""
+    penetrative_radiation(forcing)
+
+Return the radiation absorbed below the surface that `forcing` applies, or `nothing`.
+"""
+penetrative_radiation(forcing) = nothing
+penetrative_radiation(radiation::AbstractRadiativeForcing) = radiation
+
+function penetrative_radiation(forcings::Union{Tuple, NamedTuple})
+    for forcing in forcings
+        radiation = penetrative_radiation(forcing)
+        isnothing(radiation) || return radiation
+    end
+    return nothing
+end
+
+"""
+    surface_radiative_flux(i, j, grid, radiation)
+
+Return the flux of the forced tracer that `radiation` carries through the surface of column `(i, j)`, positive upward.
+"""
+function surface_radiative_flux end
+
+"""
+    transmitted_fraction(i, j, grid, radiation, d)
+
+Return the fraction of the surface radiative flux of `radiation` that reaches the depth `d` in column `(i, j)`.
+"""
+function transmitted_fraction end
+
+"""
+    transmitted_fraction_derivative(i, j, grid, radiation, d)
+
+Return the derivative of `transmitted_fraction` with respect to the depth `d` in column `(i, j)`.
+"""
+function transmitted_fraction_derivative end
+
+"""
+    transmitted_thickness(i, j, grid, radiation, h)
+
+Return the integral of `transmitted_fraction` from the surface to the depth `h` in column `(i, j)`.
+"""
+function transmitted_thickness end
+
+# Flux through the bottom face of cell k, positive upward; the bottom cell absorbs what reaches it
+@inline function radiative_flux(i, j, k, grid, radiation)
+    d = znode(i, j, size(grid, 3) + 1, grid, Center(), Center(), Face()) - znode(i, j, k, grid, Center(), Center(), Face())
+    J = surface_radiative_flux(i, j, grid, radiation) * transmitted_fraction(i, j, grid, radiation, d)
+    return ifelse(inactive_cell(i, j, k - 1, grid), zero(J), J)
+end
+
+@inline (radiation::AbstractRadiativeForcing)(i, j, k, grid, clock, fields) = - ∂zᶜᶜᶜ(i, j, k, grid, radiative_flux, radiation)
+
+"""
+    surface_radiative_buoyancy_flux(i, j, grid, radiation, buoyancy, fields)
+
+Return the buoyancy flux of `radiation` through the surface of column `(i, j)`, positive upward.
+"""
+@inline surface_radiative_buoyancy_flux(i, j, grid, ::Nothing, buoyancy, fields) = zero(grid)
+@inline surface_radiative_buoyancy_flux(i, j, grid, radiation::AbstractRadiativeForcing, ::BuoyancyTracerFormulation, fields) = surface_radiative_flux(i, j, grid, radiation)
+
+@inline function surface_radiative_buoyancy_flux(i, j, grid, radiation::AbstractRadiativeForcing, buoyancy::BuoyancyForce{<:SeawaterBuoyancy}, fields)
+    b = buoyancy.formulation
+    T, S = get_temperature_and_salinity(b, fields)
+    α = thermal_expansionᶜᶜᶠ(i, j, size(grid, 3) + 1, grid, b.equation_of_state, T, S)
+    return b.gravitational_acceleration * α * surface_radiative_flux(i, j, grid, radiation)
+end
diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl
@@ -8,7 +8,7 @@ using Oceananigans.BoundaryConditions: FieldBoundaryConditions, needs_implicit_s
 using Oceananigans.BuoyancyFormulations: materialize_buoyancy, validate_buoyancy
 using Oceananigans.DistributedComputations: Distributed
 using Oceananigans.Fields: CenterField, Field, TracerFields, ZeroField, tracernames
-using Oceananigans.Forcings: model_forcing
+using Oceananigans.Forcings: model_forcing, penetrative_radiation
 using Oceananigans.Grids: AbstractHorizontallyCurvilinearGrid, Center, Face, MutableVerticalDiscretization,
                           architecture, halo_size
 using Oceananigans.ImmersedBoundaries: GridFittedBoundary, ImmersedBoundaryGrid
@@ -18,7 +18,8 @@ using Oceananigans.TimeSteppers: AbstractLagrangianParticles, Clock, TimeStepper
                                  time_discretization, update_state!
 using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization, add_closure_specific_boundary_conditions,
                                        build_closure_fields, closure_required_tracers, implicit_diffusion_solver,
-                                       initialize_closure_fields!, validate_closure, with_tracers
+                                       initialize_closure_fields!, validate_closure, with_penetrative_radiation,
+                                       with_tracers
 using Oceananigans.Utils: tupleit
 
 PressureField(grid) = (; pHY′ = CenterField(grid))
@@ -338,6 +339,7 @@ function build_hydrostatic_free_surface_model(grid, ::Val{tracer_names}, timeste
     # constructed inside tendency kernels, ensuring field dependency indices
     # are consistent between materialization and tendency computation.
     model_fields = merge(hydrostatic_fields(velocities, free_surface, tracers), auxiliary_fields)
+    closure = with_penetrative_radiation(closure, penetrative_radiation(forcing))
     forcing = model_forcing(forcing, model_fields, prognostic_fields)
     transport_velocities = transport_velocity_fields(velocities)
 

diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl
@@ -7,7 +7,7 @@ using Oceananigans.BoundaryConditions: MixedBoundaryCondition, needs_implicit_so
 using Oceananigans.BuoyancyFormulations: materialize_buoyancy, validate_buoyancy
 using Oceananigans.DistributedComputations: Distributed
 using Oceananigans.Fields: CenterField, Field, VelocityFields, ZFaceField, tracernames
-using Oceananigans.Forcings: model_forcing
+using Oceananigans.Forcings: model_forcing, penetrative_radiation
 using Oceananigans.Grids: architecture, topology
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid
 using Oceananigans.Models: AbstractModel, extract_boundary_conditions, materialize_free_surface, materialize_tracers,
@@ -17,7 +17,7 @@ using Oceananigans.TimeSteppers: AbstractLagrangianParticles, Clock, TimeStepper
                                  time_discretization, update_state!
 using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization, build_closure_fields,
                                        implicit_diffusion_solver, initialize_closure_fields!, validate_closure,
-                                       with_tracers
+                                       with_penetrative_radiation, with_tracers
 using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: FlavorOfCATKE
 using Oceananigans.Utils: tupleit
 
@@ -349,6 +349,7 @@ function build_nonhydrostatic_model(grid, ::Val{tracer_names}, timestepper,
     timestepper = TimeStepper(timestepper, grid, prognostic_fields; implicit_solver)
 
     # Materialize forcing for model tracer and velocity fields.
+    closure = with_penetrative_radiation(closure, penetrative_radiation(forcing))
     forcing = model_forcing(forcing, model_fields, prognostic_fields)
 
     # Initialize boundary transport container

diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl
@@ -283,8 +283,8 @@ include("Simulations/Simulations.jl")
 include("Coriolis/Coriolis.jl")
 include("BuoyancyFormulations/BuoyancyFormulations.jl")
 include("StokesDrifts.jl")
-include("TurbulenceClosures/TurbulenceClosures.jl")
 include("Forcings/Forcings.jl")
+include("TurbulenceClosures/TurbulenceClosures.jl")
 include("Biogeochemistry.jl")
 
 # TODO: move above

diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl
@@ -148,6 +148,18 @@ function hydrostatic_turbulent_kinetic_energy_tendency end
 function buoyancy_force end # to be defined by models that want to use TurbulenceClosures
 function buoyancy_tracers end # to be defined by models that want to use TurbulenceClosures
 
+#####
+##### Interface for radiation absorbed below the surface
+#####
+
+"""
+    with_penetrative_radiation(closure, radiation)
+
+Return `closure` with `radiation` handed to every closure in it that adjusts for radiation absorbed below the surface.
+"""
+with_penetrative_radiation(closure, radiation) = closure
+with_penetrative_radiation(closures::Tuple, radiation) = map(closure -> with_penetrative_radiation(closure, radiation), closures)
+
 #####
 ##### Fallback: flux = 0
 #####

diff --git a/...ce_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl b/...ce_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl
@@ -15,13 +15,15 @@ using Oceananigans.BoundaryConditions: BoundaryConditions, DefaultBoundaryCondit
 using Oceananigans.BuoyancyFormulations: BuoyancyForce, BuoyancyTracer, SalinitySeawaterBuoyancy, SeawaterBuoyancy,
                                          TemperatureSeawaterBuoyancy, top_buoyancy_flux, ∂z_b
 using Oceananigans.Fields: CenterField, XFaceField, YFaceField, ZFaceField, ZeroField
+using Oceananigans.Forcings: surface_radiative_buoyancy_flux, transmitted_fraction, transmitted_fraction_derivative,
+                             transmitted_thickness
 using Oceananigans.Grids: Center, Face, inactive_cell, inactive_node, peripheral_node, static_column_depthᶜᶜᵃ
 using Oceananigans.Operators: Δzᶜᶜᶜ, Δzᶜᶠᶠ, Δzᶠᶜᶠ, Δz⁻¹ᶜᶜᶜ, Δz⁻¹ᶜᶠᶜ, Δz⁻¹ᶠᶜᶜ, ℑxᶜᵃᵃ, ℑxᶠᵃᵃ, ℑyᵃᶜᵃ, ℑyᵃᶠᵃ, ℑzᵃᵃᶜ, ℑzᵃᵃᶠ,
                               ∂zᶜᶠᶠ, ∂zᶠᶜᶠ
 using Oceananigans.TimeSteppers: TimeSteppers, time_discretization
 using Oceananigans.TurbulenceClosures: TurbulenceClosures, AbstractScalarDiffusivity, VerticalFormulation,
                                        VerticallyImplicitTimeDiscretization, buoyancy_force, buoyancy_tracers,
-                                       dissipation, getclosure, shear_production
+                                       dissipation, getclosure, shear_production, with_penetrative_radiation
 using Oceananigans.Utils: Utils, f32_safe_cbrt, launch!, named_tuple, prettysummary
 
 const c = Center()
@@ -145,6 +147,7 @@ get_top_tracer_bcs(::SeawaterBuoyancy, tracers) = (T = tracers.T.boundary_condit
 get_top_tracer_bcs(::TemperatureSeawaterBuoyancy, tracers) = (; T = tracers.T.boundary_conditions.top)
 get_top_tracer_bcs(::SalinitySeawaterBuoyancy, tracers)    = (; S = tracers.S.boundary_conditions.top)
 
+include("penetrative_radiation.jl")
 include("tke_top_boundary_condition.jl")
 include("catke_vertical_diffusivity.jl")
 include("catke_mixing_length.jl")

diff --git a/...osures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_equation.jl b/...osures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_equation.jl
@@ -22,14 +22,13 @@ end
 @inline TurbulenceClosures.dissipation(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, args...) = zero(grid)
 
 @inline function dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure::FlavorOfCATKE, velocities, tracers,
-                                             buoyancy, surface_buoyancy_flux)
+                                             buoyancy, closure_fields)
 
     # Convective dissipation length
     Cᶜ = closure.turbulent_kinetic_energy_equation.CᶜD
     Cᵉ = closure.turbulent_kinetic_energy_equation.CᵉD
     Cˢᵖ = closure.mixing_length.Cˢᵖ
-    Jᵇ = surface_buoyancy_flux
-    ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, Jᵇ)
+    ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, closure_fields)
 
     # "Stable" dissipation length
     Cˡᵒ = closure.turbulent_kinetic_energy_equation.CˡᵒD
@@ -51,7 +50,7 @@ end
 @inline function dissipation_rate(i, j, k, grid, closure::FlavorOfCATKE,
                                   velocities, tracers, buoyancy, closure_fields)
 
-    ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, closure_fields.Jᵇ)
+    ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, closure_fields)
     e = tracers.e
     FT = eltype(grid)
     eᵢ = @inbounds e[i, j, k]
@@ -96,7 +95,7 @@ end
                                tke::CATKEEquation, closure::CATKEVD,
                                buoyancy, top_tracer_bcs, top_velocity_bcs)
 
-    wΔ³ = top_convective_turbulent_velocity_cubed(i, j, grid, clock, fields, buoyancy, top_tracer_bcs)
+    wΔ³ = top_convective_turbulent_velocity_cubed(i, j, grid, clock, fields, buoyancy, top_tracer_bcs, closure.penetrative_radiation)
     u★ = friction_velocity(i, j, grid, clock, fields, top_velocity_bcs)
 
     Cᵂu★ = tke.Cᵂu★