Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

additional factor when converting surface heat fluxes from dry to moist theta #1259

Merged
merged 4 commits into from
Dec 29, 2020
Merged

additional factor when converting surface heat fluxes from dry to moist theta #1259

merged 4 commits into from
Dec 29, 2020

Conversation

matzegoebel
Copy link
Contributor

@matzegoebel matzegoebel commented Jul 22, 2020 •
edited by dudhia
Loading

To close the dry theta budget at the lowest grid point when use_theta_m is enabled, the surface heat fluxes, and in the case of km_opt=5 also the non-local heat flux, are multiplied with 1+Rv/Rd*qv.

TYPE: bug fix

KEYWORDS: surface heat flux, use_theta_m, moist theta, non-local flux, 3D TKE

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES: The conversion of surface heat fluxes from dry to moist theta in WRF is inconsistent. This leads to an unclosed dry theta budget at the lowest grid point. The surface sensible heat flux in the routines vertical_diffusion_2 and vertical_diffusion_implicit should be multiplied with as is derived below.
Using Reynold's decomposition and averaging we can show for the perturbation moist theta:

with water vapor mixing ratio q and

This yields for the converted surface heat flux:

The triple correlation in the equation above cannot be easily estimated in WRF, but judging from my tests it probably has a negligible impact.
Instead of the given equation, WRF uses:

We therefore multiplied the surface heat flux with , when use_theta_m=1.
For km_opt=5, also the contribution from the non-local heat flux is multiplied with a correction factor.

LIST OF MODIFIED FILES: dyn_em/module_diffusion_em.F

TESTS CONDUCTED: checked the dry theta budget of a simulation with use_theta_m=1: simple idealized LES once with model-computed surface heat flux and once with prescribed heat fluxes, RRTMG radiation, no microphysics, zero eddy diffusivities. Changed the code to output advection, SGS diffusion and radiation tendencies for dry theta (regardless of use_theta_m and without changing the tendencies used by the model). These forcing tendencies are averaged online over 30 minutes. The sum of the forcing terms is plotted against the actual tendency of dry theta (change of instantaneous dry theta between the half-hourly time stamps divided by 30 minutes) for all gridpoints and time stamps. All tendency terms are divided by the dry air mass to obtain more familiar units. The modified code fixes the issue of an unclosed dry theta budget at the lowest gridpoint (see attached figures). The simulations were repeated for km_opt=5.

RELEASE NOTE: To close the dry theta budget at the lowest grid point when use_theta_m is enabled, the surface heat fluxes, and in the case of km_opt=5 also the non-local heat flux, are multiplied with 1+Rv/Rd*qv.

test_fluxes_flat_thd_3rd+5th
test_fluxes_flat_thdm_3rd+5th
test_fluxes_flat_fixed-sfc-fluxes_thdm_3rd+5th

@davegill
Copy link
Contributor

@skamaroc @dudhia @weiwangncar
Folks,
What is the best way for us to include Joe in this conversation? Does he do github?

@weiwangncar
Copy link
Collaborator

@davegill Send him an email, and he can view the PR.

@dudhia
Copy link
Collaborator

dudhia commented Jul 22, 2020

I checked and it is the same as Bill suggested.
This is a correct equation fix in my view. Will check code fix now.

@davegill
Copy link
Contributor

@MATZE77 @dudhia @weiwangncar
Matthias,
You have added two IF tests inside of DO loops. Could you re-write the code to put the IF tests outside of the DO loops?

@matzegoebel
Copy link
Contributor Author

@davegill
there is no additional if test, just two additional ELSE clauses compared with the original code, so I don't quite get what you mean. Can you elaborate?

@dudhia
Copy link
Collaborator

dudhia commented Jul 23, 2020

This looks good to me. Especially appreciate hearing about the conservation aspect of this fix.
Someone should run isfflx=0 with a fixed heat flux and moisture case to fully check the change.
It also applies to the km_opt=5 3d tke scheme.

@davegill davegill changed the base branch from release-v4.2.1 to release-v4.2.2 July 23, 2020 23:03
@davegill
Copy link
Contributor

davegill commented Jul 31, 2020 •
edited
Loading

@dudhia @weiwangncar @mgduda
Folks,
Take a look at the text of this PR commit message. I did not know that you can insert this text to manufacture equations on the fly?
Example 1: https://latex.codecogs.com/gif.latex?1-\frac{1}{2}+\frac{1}{3}...
becomes

Example 2: https://latex.codecogs.com/svg.latex?1-\frac{1}{2}+\frac{1}{3}...
becomes

Example 3: https://latex.codecogs.com/png.latex?1-\frac{1}{2}+\frac{1}{3}...
becomes

@weiwangncar
Copy link
Collaborator

Pretty cool!

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Aug 10, 2020 •
edited
Loading

@dudhia
What do you mean with running a test case to fully check the change? What exactly do you want to have tested/checked?
I could run the same test case as described in the initial PR message but for all combinations of isfflx in {0,1} and km_opt in {2,5}. Would that suffice?

@dudhia
Copy link
Collaborator

dudhia commented Aug 10, 2020 via email

@matzegoebel
Copy link
Contributor Author

So now I repeated the runs described in the initial PR message with isfflx in {1,2} (radiation scheme vs. specified heat flux} and km_opt in {2,5} once with the original code and once with the modified code.
org_hfx_2_None_3rd+5th
org_hfx_2_0 15_3rd+5th
org_hfx_5_None_3rd+5th
org_hfx_5_0 15_3rd+5th
fixed_hfx_2_None_3rd+5th
fixed_hfx_2_0 15_3rd+5th
fixed_hfx_5_None_3rd+5th
fixed_hfx_5_0 15_3rd+5th

@matzegoebel
Copy link
Contributor Author

For km_opt=5 and isfflx=2, there is still some inconsistency visible close to the surface. In my opinion, the non-local flux also needs to be converted to theta_m. I added the conversion factor in the latest commit. That fixes the issue as can be seen here:
fixed_hfx_nlflux_5_None_3rd+5th
fixed_hfx_nlflux_5_0 15_3rd+5th

1 similar comment
@matzegoebel
Copy link
Contributor Author

For km_opt=5 and isfflx=2, there is still some inconsistency visible close to the surface. In my opinion, the non-local flux also needs to be converted to theta_m. I added the conversion factor in the latest commit. That fixes the issue as can be seen here:
fixed_hfx_nlflux_5_None_3rd+5th
fixed_hfx_nlflux_5_0 15_3rd+5th

dudhia
dudhia reviewed Aug 11, 2020
View reviewed changes
dyn_em/module_diffusion_em.F Outdated
@@ -8044,6 +8050,9 @@ SUBROUTINE vertical_diffusion_implicit(ru_tendf, rv_tendf, rw_tendf, rt_tendf,&
DO i = i_start, i_end
nlflux_rho(i,k,j) = nlflux(i,k,j) * &
(fnm(k)*rho (i,k,j)+fnp(k)*rho (i,k-1,j))
IF(config_flags%use_theta_m == 1)THEN
nlflux_rho(i,k,j) = nlflux_rho(i,k,j)*(1.+rvovrd*moist(i,k,j,P_QV))
ENDIF
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure about this. Check if this can be done instead with the line that sets sflux(i,j) based on hfx.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be equivalent. When looking at how the non-local flux is computed, you see that all terms that contribute are proportional to sflux (via the variables amf1, amf2, bmf2, hfxpbl):

WRF/dyn_em/module_diffusion_em.F

Lines 5206 to 5226 in 6d0090f

sflux0 = (a11+a12*sfcfracn)*sflux(i,j)
snlflux0 = nlfrac*sflux0
amf1 = snlflux0/sfcfracn
amf2 = -snlflux0/(mlfrac-sfcfracn)
bmf2 = -mlfrac*amf2
amf3 = snlflux0*enlfrac2(i,j)/deltaoh(i,j)
bmf3 = -amf3*mlfrac
hfxpbl(i,j) = amf3+bmf3
pth1 = pthnl(delxy,hpbl(i,j))
hfxpbl(i,j) = hfxpbl(i,j)*pth1
DO k = kts, ktf
zfacmf(i,k,j) = max((zq(i,k+1,j)/hpbl(i,j)),zfmin)
IF(pblflg(i,j).and.k.LT.kpbl(i,j)) THEN
IF(zfacmf(i,k,j).LE.sfcfracn) THEN
nlflux(i,k,j) = amf1*zfacmf(i,k,j)
ELSE IF (zfacmf(i,k,j).LE.mlfrac) THEN
nlflux(i,k,j) = amf2*zfacmf(i,k,j)+bmf2
ENDIF
nlflux(i,k,j) = nlflux(i,k,j) + hfxpbl(i,j)*exp(-entfacmf(i,k,j))
nlflux(i,k,j) = nlflux(i,k,j)*pth1

I thought it would be more intuitive to view the non-local flux as being computed for dry theta and then convert it via the same equation derived in the initial PR message:

The water vapor flux term is dropped because there is no non-local water vapor flux in this scheme.

But if you still think it would be better/more readable to convert sflux instead, I can also do that.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a general comment, but I will put one here too.
The nlflux terms are proportional to the surface heat flux, so it should be possible to just modify the surface heat flux the same way you did for other options. The difference from this code is that the correcting q term should be the surface one.

@dudhia
Copy link
Collaborator

dudhia commented Aug 11, 2020 via email

@dudhia
Copy link
Collaborator

dudhia commented Aug 12, 2020 via email

@matzegoebel
Copy link
Contributor Author

ok, I get your point that the change would then be more similar to my other changes. That would however change the grid variable nlflux (in contrast to my change, which is only local). I think it might be better if that potential output variable would be the same regardless of whether theta_m is used or not. In addition, grid%nlflux is also used in the subroutine tke_bouyancy. I think that in the bouyancy tke production term, the non-local heat flux should not include the moisture factor. If we already multiplied the term before, we would need to divide it here again, which is somehow inelegant. What do you think?

@dudhia
Copy link
Collaborator

dudhia commented Aug 17, 2020 via email

@matzegoebel
Copy link
Contributor Author

ah, now I understand that there is a difference between the two approaches. You suggest using the surface q for the correction factor at all levels, while I implemented the k-level q.
In my opinion, taking the k-level q is correct. This modification is different from my other modifications because it applies to the non-local flux at all levels not only at the surface. So we cannot use an analogy argument.
I took again a look at the Zhang 2018 paper. They write the equations for the non-local flux of dry potential temperature. To convert this flux to moist potential temperature we need to multiply by a correction factor using the k-level q as in the equation that I wrote above:

How the dry non-local flux is calculated, i.e. that is proportional to the surface heat flux, should be irrelevant for the conversion.
We could also involve the developers of this package to come to a decision.

@dudhia
Copy link
Collaborator

dudhia commented Aug 25, 2020 via email

@matzegoebel
Copy link
Contributor Author

Why do you say, that internally at other levels theta_m is the variable already? In the non-local flux subroutine, I can't see any connection to theta_m, neither at the surface nor higher up. The routine uses the input variable th_phy, which is dry theta.
The method uses a lot of constant parameters that are derived from a coarse-grained LES. These parameters are tied to the dry theta flux. If we just convert the surface flux, these parameters would not fit anymore. Instead, I suggest converting the whole non-local flux after computation, which requires the k-level q. But maybe we need a third opinion on this?

@dudhia
Copy link
Collaborator

dudhia commented Sep 14, 2020 via email

@matzegoebel
Copy link
Contributor Author

yes, var_mix corresponds to theta_m. But that variable is only used in the local part (Deardorff type) of the scheme and has nothing to do with the non-local part. The non-local flux is computed separately and its divergence is added as a source term to the implicit solution in line 8107:

WRF/dyn_em/module_diffusion_em.F

Lines 8102 to 8108 in f311cd5

DO k = kts+1, ktf-1
rdz_w = -g/dnw(k)/(c1h(k)*mu(i,j) + c2h(k))
a(k) = -rdz_w*rdz(i,k,j)*xkxavg(i,k,j)*dt
b(k) = 1.+rdz_w*(rdz(i,k+1,j)*xkxavg(i,k+1,j)+rdz(i,k,j)*xkxavg(i,k,j))*dt
c(k) = -rdz_w*rdz(i,k+1,j)*xkxavg(i,k+1,j)*dt
d(k) = -rdz_w*(nlflux_rho(i,k+1,j)-nlflux_rho(i,k,j))*dt + var_mix(i,k,j)
ENDDO

In my opinion, the var_mix in line 8107 belongs to the local part of the scheme. So maybe that's confusing, because it stands in the same line as the non-local source term. The non-local flux in line 8107 nevertheless corresponds to a dry theta flux at all vertical levels.

@dudhia
Copy link
Collaborator

dudhia commented Sep 14, 2020 via email

@matzegoebel
Copy link
Contributor Author

I don't think so, because the change in my line 8054 only affects the non-local flux which is only taken into account for k > kts, while the change in line 8107 converts the surface heat flux, so for k=kts. These are two separate issues.

@dudhia
Copy link
Collaborator

dudhia commented Sep 14, 2020 via email

@serafis3
Copy link

Dear all, I'm the PI of @matzegoebel's project in Innsbruck and I'm jumping in because we had a few discussions about this issue recently. I don't have the time to dig deep into the WRF code, but I'm hoping to help with some considerations on general principles. So...

(1) At any vertical level k, the following flux decomposition holds true (neglecting third-order moments):

(2) If now we want the turbulent flux to be parameterized with a local plus a non-local component, we simply add the non-local fluxes where appropriate on the RHS:

(3) The above can be restated as:

I guess this is what @matzegoebel suggested back on 12 August. I interpret the equation as follows: Regardless of how the non-local contributions to turbulent flux are modelled, consistency with the local tendency of theta_m requires them to be weighted with the local values of theta and q. Does this sound convincing?

I do see the uncomfortable aspect of the thing: Even if the non-local term remains proportional to the surface fluxes, now its impact is not constant throughout the BL. This is a side-effect of reconstructing the fluxes of theta_m from those of theta and q.

@dudhia
Copy link
Collaborator

dudhia commented Sep 21, 2020 via email

@serafis3
Copy link

serafis3 commented Sep 22, 2020 •
edited
Loading

Thanks for getting back! Sure, there is no doubt about the purpose and origin of the nonlocal terms. To clarify my thought: In Eq (3) above, the nonlocal fluxes of (dry) theta and q (if existing) remain just as they are in the Shanghai Met Service scheme (kmopt=5). But there is an open issue with theta_m conservation. A closure such as:

is fine with theta in a dry atmosphere. But if one wishes to conserve theta_m in a moist atmosphere, it seems the closure should become:

which expands to:

This doesn't change the parameterization of the nonlocal fluxes at all, it only adds two terms to account for the influence of moisture on theta_m. If there is no need to conserve theta_m, those two terms can be left out... it can be a legitimate modelling choice, if the consequences are clear.

I'm not sure if the authors of the SMS scheme considered the case of moist air. It doesn't seem so. Yes, getting their feedback would be precious!

Stefano

@matzegoebel
Copy link
Contributor Author

@serafis3, thanks for your support. I totally agree.
In order to make equations be displayed, the URL has to be embedded in an image command like this:
<img src="https://latex.codecogs.com/gif.latex?\overline{w'\theta_{m}'}" />

@dudhia
Copy link
Collaborator

dudhia commented Sep 22, 2020 via email

@matzegoebel
Copy link
Contributor Author

I just pushed commit 654069f which leads to an even better conservation of dry theta. Instead of multiplying the non-local flux in km_opt=5 directly with q(k), the vertical derivative of nlflux is first calculated and then multiplied with q(k). This is more consistent, because nlflux is on w-levels, while q is on mass levels, as is the vertical derivative of nlflux. I did a few simple CBL simulations to show the difference between not modifying the non-local flux at all, using q(k) in the conversion factor, and using q(kts) in the conversion factor, as Jimy suggests. The following figure shows the RMSE of comparing the sum of all forcings with the actual dry theta tendency. Each lines represents an average over 10 minutes, x, y and over 4 ensemble simulations with perturbed initial conditions.

rmse_nlflux_conversion_tmean

Using q(k) in the conversion, seems to lead to a slightly lower error at most levels, especially at the beginning of the simulation,
The results at later times need be taken with care, since the different simulations are not exactly equivalent due to the different conversion factors. Therefore, I argue that the q(k) version should be the correct one (in addition to the theoretical arguments that we presented earlier).
In any way, both versions are much better than not modifying the non-local flux (the green line in the left plot actually extends to 5e-5, thus an order of magnitude worse).

@dudhia
Copy link
Collaborator

dudhia commented Oct 12, 2020 via email

@matzegoebel
Copy link
Contributor Author

yes that's right, in line 8110

@davegill
Copy link
Contributor

@dudhia @weiwangncar
What is the status of this PR? Are we waiting for external comments from anyone?

@dudhia
Copy link
Collaborator

dudhia commented Oct 23, 2020 via email

@dudhia
Copy link
Collaborator

dudhia commented Dec 16, 2020

@matzegoebel I want to finalize for 4.2.2 bug-fix which is about ready.
Can we do it my way for km_opt=5? Which is to use kts for qv. Thanks.

@dudhia
Copy link
Collaborator

dudhia commented Dec 16, 2020

In fact, it seems like it is already this way, so I will go ahead and approve soon after you confirm it is ready now.

@matzegoebel
Copy link
Contributor Author

ok

@dudhia
Copy link
Collaborator

dudhia commented Dec 16, 2020 via email

@davegill
Copy link
Contributor

@matzegoebel @dudhia
Jimy and Matthias,
We will do the merge.

@matzegoebel
Copy link
Contributor Author

I just see that I have a bunch of equations and plots in the commit message. Should remove those?

@davegill
Copy link
Contributor

@matzegoebel
If those equations and figures are still valid, they can certainly remain.

@weiwangncar weiwangncar merged commit 61029b1 into wrf-model:release-v4.2.2 Dec 29, 2020
@dudhia
Copy link
Collaborator

dudhia commented Dec 29, 2020 via email

@matzegoebel
Copy link
Contributor Author

matzegoebel commented Dec 30, 2020 •
edited
Loading

@dudhia

It worked, but now I realize the change I wanted, k-> kts, is not in there yet.

Oh yes you're right. Sorry, I overlooked this.

We also had Xu Zhang suggest he preferred it to be in.

He changed his mind several times during the email discussion, but in his latest mail on October 26th he wrote "Hi Matthias, I totally agree with your opinion in the view of theory." This refers to the discussion we had about whether to use q(k) or q(kts). So I guess he actually supports my view. But we can also ask him again for a final statement, if you're not convinced.

@dudhia
Copy link
Collaborator

dudhia commented Dec 30, 2020 via email

vlakshmanan-scala pushed a commit to scala-computing/WRF that referenced this pull request Apr 4, 2024
@matzegoebel @vlakshmanan-scala
Add factor when converting surface heat fluxes from dry to moist theta (
fb22b5a 
wrf-model#1259)

TYPE: bug fix

KEYWORDS: surface heat flux, use_theta_m, moist theta, non-local flux, 3D TKE

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES: The conversion of surface heat fluxes from dry to moist theta in WRF is inconsistent. This leads to an unclosed dry theta budget at the lowest grid point. The surface sensible heat flux in the routines `vertical_diffusion_2` and `vertical_diffusion_implicit` should be multiplied with <img src="https://latex.codecogs.com/svg.latex?1+\frac{R_{v}}{R_{d}}q"  /> as is derived below.
Using Reynold's decomposition and averaging we can show for the perturbation moist theta:
<img src="https://latex.codecogs.com/svg.latex?\theta_{m}'=(\theta(1+q\alpha))'=\theta(1+q\alpha)-\overline{\theta(1+q\alpha)}=\theta-\overline{\theta}+\alpha(q\theta-\overline{q\theta})=\\=\theta'+\alpha\left[(\overline{\theta}+\theta')(\overline{q}+q')-\overline{\theta}\overline{q}-\overline{\theta'q'}\right]=\theta'(1+\alpha\overline{q})+\alpha\left[\overline{\theta}q'+\theta'q'-\overline{\theta'q'}\right]"  /> 
with water vapor mixing ratio q and
 <img src="https://latex.codecogs.com/svg.latex?\alpha=\frac{R_{v}}{R_{d}}\approx1.61"  />

This yields for the converted surface heat flux:
<img src="https://latex.codecogs.com/svg.latex?\overline{w'\theta_{m}'}_0=(1+\alpha\overline{q}_0)\overline{w'\theta'}_0+\alpha\overline{\theta}_0\,\overline{w'q'}_0+\alpha\overline{w'\theta'q'}_0\approx(1+\alpha\overline{q}_0)\overline{w'\theta'}_0+\alpha\overline{\theta}_0\,\overline{w'q'}_0"  />
The triple correlation in the equation above cannot be easily estimated in WRF, but judging from my tests it probably has a negligible impact.
Instead of the given equation, WRF uses:
<img src="https://latex.codecogs.com/svg.latex?\overline{w'\theta_{m}'}_0=\overline{w'\theta'}_0+\alpha\overline{\theta}_0\overline{w'q'}_0 "  />
We therefore multiplied the surface heat flux with <img src="https://latex.codecogs.com/svg.latex?1+\alpha\overline{q}_0" />, when `use_theta_m=1`.
For `km_opt=5`, also the contribution from the non-local heat flux is multiplied with a correction factor.

LIST OF MODIFIED FILES: dyn_em/module_diffusion_em.F 

TESTS CONDUCTED: checked the dry theta budget of a simulation with use_theta_m=1: simple idealized LES once with model-computed surface heat flux and once with prescribed heat fluxes, RRTMG radiation, no microphysics, zero eddy diffusivities. Changed the code to output advection, SGS diffusion and radiation tendencies for dry theta (regardless of use_theta_m and without changing the tendencies used by the model). These forcing tendencies are averaged online over 30 minutes. The sum of the forcing terms is plotted against the actual tendency of dry theta  (change of instantaneous dry theta between the half-hourly time stamps divided by 30 minutes) for all gridpoints and time stamps. All tendency terms are divided by the dry air mass to obtain more familiar units. The modified code fixes the issue of an unclosed dry theta budget at the lowest gridpoint (see attached figures). The simulations were repeated for km_opt=5.

RELEASE NOTE: To close the dry theta budget at the lowest grid point when use_theta_m is enabled, the surface heat fluxes, and in the case of km_opt=5 also the non-local heat flux, are multiplied with 1+Rv/Rd*qv.

![test_fluxes_flat_thd_3rd+5th](https://user-images.githubusercontent.com/17001470/88155314-389f7380-cc08-11ea-9dac-f7ec6cceeaed.png)
![test_fluxes_flat_thdm_3rd+5th](https://user-images.githubusercontent.com/17001470/88155312-3806dd00-cc08-11ea-8226-4ed46c5c556e.png)
![test_fluxes_flat_fixed-sfc-fluxes_thdm_3rd+5th](https://user-images.githubusercontent.com/17001470/88155315-39380a00-cc08-11ea-8cbd-6b9fd6217c1c.png)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@davegill davegill davegill left review comments

@dudhia dudhia dudhia approved these changes

Assignees
No one assigned
Labels
bug Dynamics ARW
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
5 participants
@matzegoebel @davegill @weiwangncar @dudhia @serafis3