Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

Negative Q2 values over very high elevations and snow


New member

I am running MPAS 7.3 on a 10 km mesh with NOAH (same version as in WRF) and I am facing negative Q2 values over very high elevations and snow covered areas in the southern Andes as described in the WRF forum in 2021. The ECMWF analysis used to initialize shows small values (0.1 g/kg) but no negative values of QVAPOR in the lowest model levels.

I made lots of tests and I figured out that when changing the sign of the "DEW" variable in module_sf_noahlsm.F, line 795 lots of negative Q2 values are gone.
For snow-free areas, e.g. in the NOPAC subroutine, DEW is set to -ETP1 while there is a comment "REINITIALIZE ETP1 TO ZERO" but I do not see this in the code.

I should add that at these grid cells the difference between TSK and the lowest model level temperature in the ECMWF analysis is in the range of 15 K (TSK << T3D(i,1,j))

I am absolutely not an expert in soil processes so I am curious to hear your opinion on that.
@ thomasuh

Thanks for the detailed description of this issue. I looked at the code module_sf_noahlsm.F, in which the value of upward ETA (latent heat flux at the surface) is negative. Positive value of DEW is equivalent to condensation (i.e., when evapotranspiration is negative).

Therefore, I believe the line 795:


is correct.

However, I am sorry that I don't have an answer why Q2 could be negative. I will talk to our expert and get back to you.
Last edited:
Below is the feedback from our expert:

"I took a look at this issue, but unfortunately I am not very familiar with the Noah process. My first impression is that it might be caused by the large near-surface temperature gradient but based on the Q2 calculation equation, there should not be negative values unless other variables in the Q2 equation have some unreasonable value. My suggestion would be checking the Q2 equation and related variables to see if they have reasonable values."

Personally I am really curious why Q2 could be negative and why surface could be so cold compared to the air above. I am suspicious that there might exist some issues in snow physics of Noah. Hopefully we can have time to debug .....

I will get back to you if I have any update.
@Ming Chen

Q2 is based on QSFC which is negative over these areas. At the same time, LH (and thus QFX) is also negative (about -10 W/m2). QSFC is calculated here (it is called Q1 in the NOAH code) from ETA_KINEMATIC. In case there is now snow (where also negative values of QSFC exist), ETA_KINEMATIC=ETA. However, the documentation at the beginning of the module says that ETA and ETA_KINEMATIC have different units.

The calculation of QSFC also involves "RCH" which is calculated here. As this is based on CH (calculated in the surface layer scheme), maybe there is an issue in the surface layer scheme in case a strong temperature gradient between TSK and T(i,1,j) exists. I tested both sf_monin_obukhov_rev and sf_monin_obukhov schemes, but with the same results.

What is also a bit unclear to me is this equation. I would assume that the first part is the Stefan–Boltzmann law, but then the constant "6.48E-08" seems to be not correct or this is a conversion factor?
Last edited:
@Ming Chen

Apparently the negative Q2 values are spin-up effects. For testing I just ran the simulation for a few minutes.
When extending the simulation duration to more than one hour, the negative Q2 values are getting less and less. I will try to run a simulation of about 4-6 hours and see what happens.
Apologize for the excitement that arose.
Many thanks for the update. We will take a look at the equation

RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0

and get back to you. Thanks again.
I have talked to our physics expert and Noah LSM expert. Unfortunately we haven't figured out yet the physical meaning of the constant " 6.48E-8" in the equation. I am sorry that I don't have an immediate answer to this question.
We will continue to pay attention to this issue and keep you updated if we found something.
Last edited:

For the equation to calculate RR, we tracked it down to an OSU 1d model by Ek and Mahrt (1991).

RGAS=287 and CP=1004 works out about right. Factor becomes 8/7.

Thereby the RR calculation is correct.
Hi Ming,

Thanks for the update.
After a couple of hours, grid cells with negative values of QSFC still remain over glacier grid cells at higher elevations. One thing I noticed in the glacier driver is the definition of ETA_KINEMATIC in the NOAH glacier subroutine where ETA_KINEMATIC is set to ETP while apparently it may have been ESNOW in earlier versions.
I've also seen that ETP3 is calculated here but is never used.
Apologize for being so persistent on this but I would like to understand what's going on or what's going wrong;)

I can upload a sample DIAG file to your Nextcloud so that you can have a look. However the file size is about 5.4 GB - is this file size still okay for upload?
Last edited:
Hi Ming,

I also stumbled over the calculation of "TH2" in NOAHDRV.
Lines 819+820 calculate the necessary factors for conversion to potential temperature based on PSFC (lowest model level pressure) and SFCPRS (middle level of the next layers. Line 821 calculate Theta on the lowest model level, however line 822 converts this back to "real" temperature "TH2".
In all the NOAH codes, "TH2" is defined as potential temperature at the lowest model level but instead it is an interpolated "real" temperature. This leads to a very small difference in the glacier code. This in turn could lead to negative values of "RAD" and thus too large negative values of ETP coming from the PENMAN subroutine.