How to properly diagnose T2m and Q2m with NoahMP?

[HelinWei-NOAA]

Unlike Noah LSM, NoahMP uses two-tile (vegetation and bare soil) approach and makes T2m and Q2m diagnosis more complicated. There are several ways to diagnose T2m and Q2m.

(1) NoahMP produces T2m and Q2m separately over vegetation and bare soil, then area averaged based on the fveg (maximum GVF) is used to get the corresponding values over the model grid.
(2) NoahMP produces aggregated Tsurf and Qsurf, and sfc_diag.f will compute T2m and Q2m based on stability outputs from third sfc_diff call. However, there are different approaches to get aggregated Tsurf and Qsurf:
(a) area averaged based on fveg,
(b) use aggregated SH and LH with the surface exchange coefficient from third stability call

Which one should be used should depend on the consistency between the fluxes and the T/Q at the first model level and surface.

From our tests, the differences can be significant with the different approaches. The current approach is (2)(a) (Tsurf/Qsurf is from area averaged). It is found that Td2m based on Qsurf is way too high over the bare soil.
Td2m over bare soil
P8ALSMZ Td2m is too high. For the other two cases (P7C and P8A) using NoahMP, Td2m seems okay. In both cases, (2)(b) is used for Qsurf. But for other regions, we found Td2m is too low and not consistent with the surface cape pattern.

Therefore the issue here is we need to determine which approach we should use and why Qsurf is so high over bare soil with method (2)(a).

Also, for (1), is the internal NoahMP method for determining T2m and Q2m for sfcdif options 3 (global) and 4 (regional) correct?

[HelinWei-NOAA]

@barlage Can you add Weizhong?

[HelinWei-NOAA]

@wzzheng90 Can you see this thread?

[barlage]

@HelinWei-NOAA I can't. Is this because Weizhong does not have a fork of the repository or is not a collaborator on your fork?

[HelinWei-NOAA]

Using Mike's UFS land driver to test Qsurf difference with 2b and 2a over SCAN sites

[wzzheng90]

Now I can see it.

[HelinWei-NOAA]

over vegetation:qsfc = (0.622eah)/(sfcprs-0.378eah)
over bared soil:qsfc = 0.622*(estg * rhsur)/(psfc-0.378*(estg*rhsur))

Why use psfc (pressure at lowest model layer) not use surface pressure over bared soil?

[HelinWei-NOAA]

It looks like sfcprs=psfc, both surface pressure

over vegetation:qsfc = (0.622_eah)/(sfcprs-0.378_eah) over bared soil:qsfc = 0.622*(estg * rhsur)/(psfc-0.378*(estg*rhsur))

Why use psfc (pressure at lowest model layer) not use surface pressure over bared soil?

[HelinWei-NOAA]

I add some calculations after the third stability call in sfc_noahmp_drv.F90 to calculate 2m T/Q

  ch2       (i)   = 0.4*0.4/(fm1(i)*fh21(i))
  ch2       (i)   = max(ch2(i), 0.00001/zf(i))
  snwdph    (i)   = snow_depth * 1000.0       ! convert from m to mm; wait after the stability call

! qsurf (i) = q1(i) + evap(i)/(con_hvapdensitych(i)wind(i))
t2mmp(i) = tsurf(i) + hflx(i)/(con_hvapdensitych2(i)wind(i))
q2mp(i) = qsurf(i) + evap(i)/(con_hvapdensitych2(i)*wind(i))

Now I can compare 2m T/Q between method 1 and 2(a). TEST1 is method 1 and TEST4 is 2(a)
You can see the difference of T2m is small but huge for Q2m

@barlage Do you have the observation to compare?

[wzzheng90]

T2m (or q2m) is usually interpolated from Tskin and Tlev1 (or qs and qlev1) with the stability. Is there any specific reason to use the surface fluxes to derive them?

[HelinWei-NOAA]

Add test2 for method 2(b)

  tsurf     (i)   = t1(i) + hflx(i)/(con_hvap*density*ch(i)*wind(i))
  qsurf     (i)   = q1(i) + evap(i)/(con_hvap*density*ch(i)*wind(i))
  t2mmp(i)        = tsurf(i) - hflx(i)/(con_hvap*density*ch2(i)*wind(i))
  q2mp(i)        = qsurf(i) - evap(i)/(con_hvap*density*ch2(i)*wind(i))

This time even T2m has some large difference, something is weird

[HelinWei-NOAA]

In theory they should be very close. If we can add sfc_diag.f to ufs land driver, then we can get the exact same 2m T/Q as UFS.

T2m (or q2m) is usually interpolated from Tskin and Tlev1 (or qs and qlev1) with the stability. Is there any specific reason to use the surface fluxes to derive them?

[wzzheng90]

The surface fluxes come from the tiles and ch(2) and some others come from the 3rd call. There would be quite large inconsistent.

[HelinWei-NOAA]

That is what the current model doing and we just make sure the flux is consistent before and after the third call.

The surface fluxes come from the tiles and ch(2) and some others come from the 3rd call. There would be quite large inconsistent.

[RongqianYang-NOAA]

'fveg' does not change with the current option (max). Unlike in Noah, the current option for 'opt_trs' (z0m/z0h) , where 'fveg' is involved, does not make too much difference. That will cause many problems over sparse canopy, especially when/if 'fveg' is used in the verification.

[HelinWei-NOAA]

To rule out the cause by fveg=shdmax, I have repeated test1, test2, test4 but dveg=1 which uses shdfac (gvf). The scenario is very similar. So the most likely reason is LH over bared soil is too high.

[HelinWei-NOAA]

Notice that the blue lines (test2 and test21) for both dveg=1 and dveg=4 are similar. For the other two tests, dveg=1 has higher values.

[HelinWei-NOAA]

Two critical inputs for the third stability call are the composited z0h and z0. How to composite them can make quite difference. opt_trs=4 is used for all of my previous tests. I ran 3 more tests (test3, test 5, test6) with 2(a) but using option 1, 2,3 for thermal roughness. Here you can see the difference. The reason we have larger difference on Q is because of the large difference on z0h.

[RongqianYang-NOAA]

First, the previous t2mmp and q2mmp calculations should be "tsurf/qsurf" - (minus) ... not + (plus)?.
Second, replace the z0m/z0h in the driver call (third) with the original z0m/z0h (used for Noah), use only Tskin from Noah MP, use the similar method for z0m/z0h for the tiles.

[RongqianYang-NOAA]

Which 'fveg' option is used for tests 4, 5, and 6? Still 'shdmax'? The differences in thermal roughness length among the options are very large indeed. All the options expect a variable 'fveg' - not max GVF - otherwise, there will be a problem. A invariable 'fveg' looks similar to the constant LAI(=3) in Noah. The variable table LAI in Noah MP is not used in the formulations.

[HelinWei-NOAA]

You are right about the sign. I have updated some plots with this correction. Now method 2(b) is more like an outlier and the difference because of zoh becomes smaller.

First, the previous t2mmp and q2mmp calculations should be "tsurf/qsurf" - (minus) ... not + (plus)?. Second, replace the z0m/z0h in the driver call (third) with the original z0m/z0h (used for Noah), use only Tskin from Noah MP, use the similar method for z0m/z0h for the tiles.

[HelinWei-NOAA]

shdmax is used for tests 4, 5 ,6

Which 'fveg' option is used for tests 4, 5, and 6? Still 'shdmax'? The differences in thermal roughness length among the options are very large indeed. All the options expect a variable 'fveg' - not max GVF - otherwise, there will be a problem. A invariable 'fveg' looks similar to the constant LAI(=3) in Noah. The variable table LAI in Noah MP is not used in the formulations.

[HelinWei-NOAA]

It seems to me that using t1/q1 to derive tsurf/qusrf is not a good way in method 2(b) because in many situations the first model level is not inside the surface layer. The similarity theory (constant flux) can't be applied between the surface and the first model level.

[RongqianYang-NOAA]

It is believed that the tiled approach can generate better fluxes. Since we can't average stability (even the cm and ch are available for both tiles in the code - after vege_flux and bare_flux calls - they are not used), it would be a good idea if we have a reliable cm/ch for the grid (from the 3rd call), it depends also on the Tskin formulation (in addition to the two roughness lengths). 2(b) relies on the accuracy of tskin, if the z0m/z0h performs ok with Noah, i.e., 2(b) puts more trust on the fluxes and the sfc_diff.f call from the atmos. model. The sfc_diff.f needs a variable 'fveg'.

[HelinWei-NOAA]

opt_trs=2, test8 zohtotal weighed by shdfac, test9 flux-based weighted
lweight=(abs(evc)+abs(tr)+abs(evg))/(abs(evc)+abs(tr)+abs(evg)+abs(evb))

result

[barlage]

some heat flux consistency checks

https://docs.google.com/presentation/d/1ryasx2V2EIjZG1YTU2CBOgLwFeL4y2VOtN-HraGok-M