Another simple radiation scheme

[milankl]

Andrew Williams @AndrewILWilliams suggested

Another simple radiation scheme which might be worth trying, and which more faithfully mimics realistic radiative cooling profiles is that of Jeevanjee and Romps (2018), and recently used in Sec 2.3 of Jake Seeley's paper (https://www.jacobtseeley.com/files/seeley2023.pdf) and Nadir's 2022 paper (https://nadirjeevanjee.com/papers/19iris.pdf -- Eq. 2).

The basic insight is that, if F is the net upwards radiative flux, the flux divergence with respect to temperature doesn't vary that much. That is, \partial{T} F is approximately a fixed function of temperature (see Eq. 1 of the Seeley paper). This means one can get radiative heating rates by calculating "\partial{T} F" (which is just a fixed function of your local gridbox temperature) and then transforming to pressure coordinates by multiplying by dT/dP (or whatever vertical coordinate SpeedyWeather uses). Then you can get your K/s heating rates.

[milankl]

@AndrewILWilliams maybe you can help me understand that, it says in Seeley, 2023

With $T_t = 200K$ and $\alpha = 0.025 W/m^2/K$, I'd discretize this with $F_k$ being the upward flux between layer $k$ and $k+1$ so that we'd have

$$Fk = F{k+1} + (Tk - T{k+1}) \alpha (Tt - T{k+1})$$

For $N=3$ vertical levels with $T_1 = 280 K, T_2 = 290 K, T_3 = 300 K$ and sea/land surface temperature $T_4 = 300K$ we'd have

$$ F_N = 0$$

because there's no temperature difference between the surface and the lowermost layer but

$$F_2 = F_N + 25 W/m^2K = 25 W/m^2K$$

so the lowermost layer doesn't receive any net radiative flux, but loses it to the layer above. Then

$$F_1 = F_2 + 22.5 W/m^2K = 47.5 W/m^2 K$$

so that layer 1 receives 25 from below but loses 47.5 to above, net losing 22.5K, so a bit less than the layer below (makes sense to me) in order to close this you'd either need to define an outer space temperature $T_0$ or, and I think that's what they mean in the paper to not define and longwave outgoing radiation but a simple relaxation term within that layer. Did I get this right? Because what you said sounds like one should use a prescribed dT/dp profile, which I'd rather not use because it's too much freedom and it doesn't seem to be used in their paper either.

[AndrewILWilliams]

Aaah, I think I misunderstood how Jake implemented it. I thought that you could just use the local temperature to get $\partial{T} F$, then calculate local $dT/dp$ on the model grid and use that to get $(g/c{p}) * \partial_{p} F$, which would be the radiative heating rate to impose (in a model with a pressure-based vertical coordinate, but maybe different for speedy). I haven't thought this through / tried to implement it myself, so apologies if I'm missing something obvious!

It seems like Jake is doing something more similar to what you've suggested, and the 'closure' is through relaxing the stratosphere to a prescribed $T_{strat}$.

[milankl]

in a model with a pressure-based vertical coordinate, but maybe different for speedy

We use sigma coordinates so $\sigma = \frac{p}{p_s}$, fraction of the surface pressure, that's alwasy easily to translate. But I thought you'd just use the vertical change in temperature, i.e. from one layer to the next directly - as they state that \" $\partial / \partial T$ is a vertical derivative". You can probably multiply both sides with $dT/dp$ to have $\partial F / \partial p$ on the left-hand side, but you'd need to do that with $\alpha(T - T_t)$ too and that looks like a detour to me?

[milankl]

Is Jack Seeley on github and we could mention him here? Otherwise I could just implement that, it looks easy enough and see what happens 😆

[jtseeley]

Hi folks! Cool idea to implement this scheme. I am not sure what makes the most sense for your model, but what Andrew suggested (multiplying by $dT/dp$ to convert from $\partial_T F$ to $\partial_p F$) makes sense to me. The cloud-resolving model I used (DAM) actually needed flux divergences in W/m $^3$, so I did the following:

1. compute the domain-mean temperature profile (not necessary for a single column)
2. use that temperature profile to compute the profile of $F\mathrm{net}$ on the model vertical grid. I did this by analytically integrating the flux divergence in temperature coordinates (equation 1 from my 2023 paper), which yields $F\mathrm{net} = -\frac{\alpha}{2}(T-T_t)^2$ (up to a constant of integration, which can be set to 0)
3. compute the flux divergence in W/m $^3$ using a finite difference of the $F_\mathrm{net}$ computed above and the model's $dz$ grid
4. heating rates at altitudes above where $T=T_t$ are calculated using the Newtonian relaxation to $T_t$ on a prescribed timescale

[jtseeley]

Happy to consult further if what I said is not clear!

[AndrewILWilliams]

I didn't know you were on Github, @jtseeley ! This is great though, I wasn't entirely sure how you'd implemented it.

[milankl]

Thanks @jtseeley for the input. I'm still confused why I should go into pressure coordinates if the equation doesn't depend on it. I've interpreted the equation maybe more like $\partial F = -\partial T \alpha (T - T_t)$ meaning that there's a change of the flux $\partial F$ between two vertical layers (but their distance doesn't matter) if there's a change in temperature $\partial T$ between them. This makes intuitively sense to me as there's no net radiative flux between two layers of the same temperature. Furthermore because of the minus sign (for $T - T_t > 0$ in the troposphere) it means that if the temperature increases from one layer to the next the radiative flux goes in the opposite direction as the warmer layer radiates more towards the colder layer than the colder layer radiates to the warmer layer. Also makes sense to me. Particularly that the distance between the two layers doesn't actually matter because you assume them to be infinitely large horizontal discs, so that there's no radial decay in intensity. That's why I find it weird to introduce a dependence on the vertical coordinate $dp$ (or its change). Following Andrew's suggestion would work, I don't disagree but then you'd need to compute a $dT/dp$ profile for the right-hand side which you could do for the local column but then the $dp$'s would drop out anyway or you could use some reference profile (e.g. a global or domain-wide average), which, however seems to defeat the point that this vertical radiative flux should only depend on the vertical profiles of that column.

[milankl]

up to a constant of integration, which can be set to 0

I believe one would want to include some land or sea surface temperature into that equation as outlined above. Otherwise if you set that constant of integration to 0 this would be equivalent of having no radiative flux into/out of the column from the surface as well as at the top. This doesn't make sense to me as the lowermost model layer should receive a radiative flux from the surface too (but whether that flux cools/warms the land or sea is then more a question of one's land and ocean model). This would then be equivalent of prescribing the boundary condition flux as that same flux but using land/sea surface temperature for $\partial T$ and $T$. For the top of the atmosphere such a boundary condition would then need to make an assumption about the temperature of "space" or I guess we can talk about the stratosphere here assuming it's not explicitly represented. I see that you don't necessarily want to use 0K here because you want to keep the uppermost level close to the $T_t$ temperature hence you introduced that tropopause relaxation, which makes sense to me. But that sounds to me as if you're essentially setting $F_0 = 0$ (the flux between uppermost layer and stratosphere/space (equivalent to a constant temperature above your uppermost layer). And then you prevent any drift in your tropopause by that relaxation term.

Please tell me if I'm misunderstanding the meaning of that equation!

[AndrewILWilliams]

Following Andrew's suggestion would work, I don't disagree but then you'd need to compute a dT/dp profile for the right-hand side which you could do for the local column but then the dp's would drop out anyway or you could use some reference profile (e.g. a global or domain-wide average), which, however seems to defeat the point that this vertical radiative flux should only depend on the vertical profiles of that column.

Hi Milan. I'm a bit confused, why would the dps drop out? You have an equation for dF/dT, but I assume the model needs K/s heating rates, in which case you need info on the atmospheric mass (which can be obtained by multiplying by the local dT/dp). I don't think pressure cancels out anywhere? I may have misunderstood though.

[jtseeley]

I would only advise using this parameterization of the net flux to get local (layerwise) radiative heating rates in the column. It is not meant as a full radiative transfer emulator — i.e., no need to worry about the numerical value of net flux into/out of the column, or to tie that boundary condition to the underlying land/ocean model. All that matters is the vertical derivative, so you can offset the whole net flux profile by any integration constant you want and you will obtain the same heating rates! I hope that's clear.

[jtseeley]

And I agree with Andrew regarding the conversion to K/s heating rates.

[milankl]

Great, just implemented that with #502!

I'm a bit confused, why would the dps drop out? You have an equation for dF/dT, but I assume the model needs K/s heating rates, in which case you need info on the atmospheric mass (which can be obtained by multiplying by the local dT/dp). I don't think pressure cancels out anywhere? I may have misunderstood though.

@AndrewILWilliams In SpeedyWeather you can either set a tendency per layer $dT/dt$ or an upward/downward flux $F$. Because the original formulation here is as flux I'm just using that. You are right that there's a $\frac{1}{dp} \approx \frac{1}{\Delta p}$ multiplied to the absorbed flux in a given layer of width $\Delta p$ (also to get indeed heating rates in $K/s$) but I don't need to compute/assume any temperature profile $dT/dp$ to get the fluxes $F$ -- the temperature differences $dT$ between two layers is enough (which eventually becomes a local, instantaenous $dT/dp$ when the flux divergence is calculated)

It is not meant as a full radiative transfer emulator — i.e., no need to worry about the numerical value of net flux into/out of the column, or to tie that boundary condition to the underlying land/ocean model. All that matters is the vertical derivative, so you can offset the whole net flux profile by any integration constant you want and you will obtain the same heating rates!

@jtseeley Yep, exactly I'm setting the flux from the uppermost layer into outer space to 0 as well as the flux into the lowermost layer from the surface. Consequently we have a divergence of a flux that's zero integrated over the column which only redistributes temperature from warmer layers to colder layers. I really only have to do

https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/027b3979665d8376021280a3519dda3c5a990518/src/physics/longwave_radiation.jl#L101-L111

and can luckily skip most of your points 1-3. However, I had to adjust the 6h tropopause relaxation time to 24h as 6h causes numerical instabilities at lower resolutions (T31, 400km global). But I guess that time scale is really only to prevent the temperature flux piling up in the uppermost layer and a 24h time scale seems to do that equally well.

Also, in the version of your 2023 paper here $\alpha$ in Table 1 has the wrong units, should be $W/m^2/K^2$ (as you also write in the text) not $W/m^2/K$ 😉

[jtseeley]

Whoops, good catch on the typo in the table! Thanks, Milan.

On Sun, Mar 24, 2024 at 3:02 AM Milan Klöwer @.***> wrote:

Great, just implemented that with #502 https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/502!

I'm a bit confused, why would the dps drop out? You have an equation for dF/dT, but I assume the model needs K/s heating rates, in which case you need info on the atmospheric mass (which can be obtained by multiplying by the local dT/dp). I don't think pressure cancels out anywhere? I may have misunderstood though.

@AndrewILWilliams https://github.com/AndrewILWilliams In SpeedyWeather you can either set a tendency per layer $dT/dt$ or an upward/downward flux $F$. Because the original formulation here is as flux I'm just using that. You are right that there's a $\frac{1}{dp} \approx \frac{1}{\Delta p}$ multiplied to the absorbed flux in a given layer of width $\Delta p$ (also to get indeed heating rates in $K/s$) but I don't need to compute/assume any temperature profile $dT/dp$ to get the fluxes $F$ -- the temperature differences $dT$ between two layers is enough.

It is not meant as a full radiative transfer emulator — i.e., no need to worry about the numerical value of net flux into/out of the column, or to tie that boundary condition to the underlying land/ocean model. All that matters is the vertical derivative, so you can offset the whole net flux profile by any integration constant you want and you will obtain the same heating rates!

@jtseeley https://github.com/jtseeley Yep, exactly I'm setting the flux from the uppermost layer into outer space to 0 as well as the flux into the lowermost layer from the surface. Consequently we have a divergence of a flux that's zero integrated over the column which only redistributes temperature from warmer layers to colder layers. I really only have to do

https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/027b3979665d8376021280a3519dda3c5a990518/src/physics/longwave_radiation.jl#L101-L111

and can luckily skip most of your points 1-3. However, I had to adjust the 6h tropopause relaxation time to 24h as 6h causes numerical instabilities at lower resolutions (T31, 400km global). But I guess that time scale is really only to prevent the temperature flux piling up in the uppermost layer and a 24h time scale seems to do that equally well.

Also, in the version of your 2023 here https://www.jacobtseeley.com/files/seeley2023.pdf $\alpha$ in Table 1 has the wrong units, should be $W/m^2/K^2$ (as you also write in the text) not $W/m^2/K$ 😉

— Reply to this email directly, view it on GitHub https://github.com/SpeedyWeather/SpeedyWeather.jl/issues/456#issuecomment-2016715913, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABHONRNOZ3G36C5BLECUQUDYZZ263AVCNFSM6AAAAABDW3O3X2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAMJWG4YTKOJRGM . You are receiving this because you were mentioned.Message ID: @.***>