[simone-silvestri]

This PR adds a penetrative radiation correction to the convecting path of CATKE. In practice the convective velocity scale $w_\bigstar$ should depend on the external bouyancy forcing which, with penetrative radiation, is not constrained only on the surface but penetrates inside the mixed layer. Therefore, we should account for the integral of this forcing within the mixed layer. We also need to compute the convective length scale, which in CATKE is a measure of the mixed layer depth, and in the case of penetrative radiation requires an inversion to compute. Here it is implemented with 3 Newton iterations which (claude) verified is enough for any practical case.

To be more detailed, in a convection event, CATKE need a measure of $w_\bigstar$ and a measure of the convective mixing length. The latter is just an estimate of the mixed layer so I will call it $h$.
As per Wagner et al 2025, the assumed relation is

$$w_\bigstar^3 \sim h J^b$$

which comes from the assumption that the turbulent flux $w^\prime b^\prime$ is linear inside the mixed layer ($b$ is uniform) and $J^b$ is the (positive) buoyancy flux at the surface.

Now, if we introduce radiation, the balance above changes Let's start with the buoyancy budget

$$\partial_t b = -\partial_z F = -\partial_z (\overline{w'b'} + F^r )$$

where

$$\overline{w'b'}(0) = J^{bc}, \qquad F(0) = J^{bc} + J^r = J^b .$$

Note that $J^r$ is always negative as the solar radiation is the only penetrating radiation and $|J^{bc}| > |J^r|$ with $J^{bc}$ positive since we are in a convecting regime (solar radiation at the surface reduces convection).

Again, we assume that (1) the layer is well mixed ($b$ is uniform in $-h < z < 0$) then $F$ should be linear in $z$, and (2) there is no entrainment at the bottom of the mixed layer (same assumption of base CATKE) which means that $\overline{w'b'}(-h) = 0$ and the flux at the bottom of the mixed layer is just the penetrating radiation

$$F(-h) = J^r \mathcal{T}(- h)$$

where $\mathcal{T}$ is the trasmitted fraction of the (surface) penetrating radiation $J^r$.
Then a linear $F$ is easily computed as

$$F(z) = J^b + (J^b - J^r \mathcal{T}(- h)) \frac{z}{h}$$

and the turbulent flux is therefore

$$\overline{w'b'}(z) = F(z) - F^r(z) = J^b + (J^b - J^r \mathcal{T}(-h)) \frac{z}{h} - J^r \mathcal{T}(- z)$$

We then integrate in the layer the turbulent flux to get

$$\int_{-h}^0 \overline{w'b'} dz = w_\bigstar = \frac{h}{2}(J^b + J^r \mathcal{T}(-h)) - J^r \int_{-h}^0 \mathcal{T} dz$$

we normalize by 2 to keep the same scaling of the original CATKE implementation with no radiation.

In the limit of radiation absorbed at the surface $\mathcal{T}(-h) = 0$ and $\int_{-h}^0\mathcal{T}dz = 0$ returns back the previous formulation and in the limit of radiation passing through ($\mathcal{T}(h) = 1$ and $\int_{-h}^0 \mathcal{T} dz = h$) we recover $w_\bigstar^3 = h(J^b - J^r) = hJ^{bc}$ meaning that radiation is not accounted for in the mixed layer buoyancy budget at all.

Now, this formulation is a bit more complicated and requires a newton solve to find $h$ which is out convective mixing length. cc @rafferrari

@glwagner I would need your expertise here to understand if this implementation is sensible or not.

[glwagner]

How does this interact with single column mode / ensemble calibration? (thinking about this a bit, perhaps it would be a single specification for all columns?)

Perhaps the more important question is how this would work if the absorption depends on a prognostic variable like chlorophyll. It looks to me like the type of penetrative_radiation somehow has to be able to compute things like transmitted_thickness?

A separate but related comment is whether this better belongs all the way up in the model struct itself. For example in Breeze:

https://github.com/NumericalEarth/Breeze.jl/blob/f33703ad9424c73287f810f35a35a52643f581fc/src/AtmosphereModels/atmosphere_model.jl#L73

I think the main question is whether a bespoke interface for radiative fluxes would be used by any other closure except CATKE.

Another option is to generalize (further....) this CATKE implementation to incorporate all temperature (thus energy) forcings. In that case one could take a more minimal approach with the closure option itself, something like radiation_adjustment=true or something...

[simone-silvestri]

CATKE implementation to incorporate all temperature (thus energy) forcings.

I am testing this generalization, but it requires a bit of a complication. If we do not have a closed-form vertical integral apparently we need to perform a complicated bisection on the column to find h which requires an extra memory field allocation and an additional kernel.

[simone-silvestri]

I think the main question is whether a bespoke interface for radiative fluxes would be used by any other closure except CATKE.

Maybe we can generalize the radiation interface all closures? This could probably be done easily. How the closure ingests this interface is then closure-dependent.

[rafferrari]

The last equation should have w*^3 rather than w^*. Other than that, we should generalize the expression to allow for entrainment. As for whether the formulation should be extended to apply to atmospheres, I am not sure it would be very useful, because we implicitly assume the flux is destabilizing throughout the whole water column. This would not apply in an atmosphere where the radiative flux can be non-monotonic.

[simone-silvestri]

I have moved the interface to the parent module TurbulenceClosures so possibly other closures can define it and use it.

[glwagner]

I have moved the interface to the parent module TurbulenceClosures so possibly other closures can define it and use it.

I believe the interface belongs in Forcings. Then if you add a "radiative" forcing, it can be automatically picked up by a closure

[codecov]

Codecov Report

❌ Patch coverage is 88.77551% with 11 lines in your changes missing coverage. Please review.
✅ Project coverage is 73.60%. Comparing base (f3654c9) to head (188562c).
⚠️ Report is 2 commits behind head on main.

Files with missing lines	Patch %	Lines
...erticalDiffusivities/catke_vertical_diffusivity.jl	81.25%	6 Missing ⚠️
...EBasedVerticalDiffusivities/catke_mixing_length.jl	86.36%	3 Missing ⚠️
...asedVerticalDiffusivities/penetrative_radiation.jl	95.83%	1 Missing ⚠️
...erticalDiffusivities/tke_top_boundary_condition.jl	75.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #6099      +/-   ##
==========================================
- Coverage   74.72%   73.60%   -1.12%     
==========================================
  Files         443      444       +1     
  Lines       28008    28055      +47     
==========================================
- Hits        20928    20651     -277     
- Misses       7080     7404     +324     

Flag	Coverage Δ
buildkite	70.74% <100.00%> (+0.06%)	⬆️
distributed/output_combining	19.89% <5.15%> (-0.03%)	⬇️
distributed/tripolar	19.89% <5.15%> (-0.03%)	⬇️
julia	70.74% <100.00%> (+0.06%)	⬆️
makie	?
metal	16.30% <47.19%> (+0.04%)	⬆️
reactant	?
sharding	?

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[navidcy]

The error in the CI will heal after #6412 is merged.