Unphysical subsurface temperatures simulated in glacierized grid cells using Noah-MP

[emilycollier]

Model Version: 4.1.2 and 4.2

Type: Bug

Source file: phys/module_sf_noahmp_glacier.F

Description: Subsurface temperatures exceeding the melting point are simulated by Noah-MP (opt_gla = 1, opt_stc = 1, opt_alb =2) in glacierized grid cells for a simulation of the Alps during summer.

The figure below shows output at an example glacierized grid cell with zero vegetation for a test case in August 2018 with a clean build of WRF 4.2. There is no snow in this grid cell after the initial 0.1 m melts away at timestep ~4,400 (dt = 9s). The problem first appears after ~27,700 timesteps in the uppermost subsurface level when: (i) all subsurface levels have warmed to the melting point from the initial ice-column temperature of 263.15 K (panel a) (ii) all ice has melted in the uppermost level (panel b)

At this point, TSLB(z=1) is not reset to the melting point in the PHASECHANGEGLACIER subroutine at L1790 to L1805, because MICE(1) == 0 and therefore IMELT(1) == 0 (panel d): ! Calculate the energy surplus and loss for melting and freezing DO J = 1,NSOIL IF (IMELT(J) > 0) THEN HM(J) = (STC(J)-TFRZ)/FACT(J) STC(J) = TFRZ_

However, excess temperature is also not used to melt ice in lower layers at every timestep further in the same subroutine at L1914 to L1941, because the amount of potentially melted ice XM(1) is initially small (e.g., 2.5764965E-03 for TSLB(1)=273.1704) and has to grow over 20+ timesteps to reach the specified threshold of 0.1 kg/m2 (panel e): ! NOW REMOVE EXCESS HEAT BY MELTING ICE IF (ANY(STC(1:4) > TFRZ) .AND. ANY(MICE(1:4) > 0.)) THEN DO J = 1,NSOIL
IF ( STC(J) > TFRZ ) THEN
HEATR(J) = (STC(J)-TFRZ)/FACT(J) XM(J) = HEATR(J)DT/HFUS DO K = 1,NSOIL IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.1) THEN IF (MICE(K) > XM(J)) THEN ! LAYER ABSORBS ALL MICE(K) = MICE(K) - XM(J) XMF = XMF + HFUS XM(J)/DT STC(K) = TFRZ

Once these conditions are met, TSLB(z=1) in glacierized grid cells repeatedly increases up to ~273.6 K before being reset. The atmosphere doesn’t “see” the unphysical values because the surface temperature in these cells is constrained at the melting point elsewhere, e.g., for opt_stc = 1, in GLACIER_FLUX if there is still any soil ice or overlying snow.

However, the apparent stability of the TSLB fluctuations depends on an unrealistically high ice albedo of 0.8. If the ice albedo is decreased to a more reasonable value (e.g., 0.3), which provides more accurate feedbacks to the atmosphere, the entire ice column eventually melts. At this point, TSLB begins to increase rapidly, exceeding ~800 K in a couple of days. Here is an example from another grid point (I=60, J=1; write freq = 30 mins):

For this gridpoint at least, TSLB (z=1) appears to grow due to a positive ground heat flux, because the negative flux that is computed with respect to TSLB(z=1) in the loop in GLACIER_FLUX is replaced by a positive one when the fluxes are re-evaluated after TGB is reset to the melting point at L1085-L1100.

In addition, precision errors in MICE at z=2 and z=3 might contribute, since instead of all ice melting completely in these levels, it is reduced to a constant minimum value of ~O(10e-5), at least when checked at L1772-L1786. This may mean that any checks on soil ice being greater than zero are satisfied, for example for the TGB reset.

Code changes tested: Setting the threshold for melting or refreezing to zero in PHASECHANGE_GLACIER removes the unphysical values in TSLB(z=1) at all glacierized grid points while total soil ice is non-zero (XM_THRESH in Figure 3, which shows output at I=60, J=1). However, I’m not sure if those thresholds are there to prevent another instability.

L1922 < IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.1) THEN > IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.) THEN L1951 < IF (J .NE. K .AND. MLIQ(K) > 0. .AND. XM(J) < -0.1) THEN > IF (J .NE. K .AND. MLIQ(K) > 0. .AND. XM(J) < 0.) THEN

There is also still the issue of how best to handle the case where all glacier ice will melt... enforce a minimum ice fraction per level and include drainage? Remove the surface temperature constraint once the top level has completely melted? These kind of changes are more involved and I guess would need developer input.

Best, Emily

[davegill]

@emilycollier @weiwangncar @dudhia Emily, You win Best PR Ever Award!

[dudhia]

Agree with Dave. Let's add @barlage in case he is not following this.

[emilycollier]

Probably TMI, but thanks ;)

Apologies for my late reply. Before making a pull request -- if that would be welcome -- I'd like to have your input on the suggested fix. I checked setting the XM threshold in phys/module_sf_noahmp_glacier.F to zero only in the block for removing excess heat (SOL1), which is actually the only change that is required to remove the unphysical subsurface glacier temperatures, and in both this block and the one for removing excess cold (SOL2): SOL1: 1922c1922 < IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.) THEN > IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.1) THEN SOL2: 1922c1922 < IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.) THEN > IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.1) THEN 1951c1951 < IF (J .NE. K .AND. MLIQ(K) > 0. .AND. XM(J) < 0.) THEN > IF (J .NE. K .AND. MLIQ(K) > 0. .AND. XM(J) < -0.1) THEN

There are some pretty big impacts from changing the excess cold threshold. For the first 34 history writes at two-hourly frequency, only TSLB below z=2 in glacier grid cells differs between the two simulations (warmer in SOL2 at z=2 and cooler in z = 3 & 4). At write 35, TSLB(z=1) differs at a single grid point, and therefore so do TGB, GRDFLX and other surface energy vars.

By the next write two hours later, TGB amongst many vars is impacted nearly domain wide, with anomalies of 1-2K after a few more hours.

So the code changes have a large impact and it would be great to have your input on which change, if any, would be best to contribute.

[dudhia]

We should include @barlage for his thoughts on this as the NoahMP code owner. Jimy

On Mon, May 25, 2020 at 6:44 AM Emily Collier notifications@github.com wrote:

Probably TMI, but thanks ;)

Apologies for my late reply. Before making a pull request -- if that would be welcome -- I'd like to have your input on the suggested fix. I checked setting the XM threshold in phys/module_sf_noahmp_glacier.F to zero only in the block for removing excess heat (SOL1), which is actually the only change that is required to remove the unphysical subsurface glacier temperatures, and in both this block and the one for removing excess cold (SOL2): SOL1: 1922c1922 < IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.) THEN

IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.1) THEN SOL2: 1922c1922 < IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.) THEN IF (J .NE. K .AND. MICE(K) > 0. .AND. XM(J) > 0.1) THEN 1951c1951 < IF (J .NE. K .AND. MLIQ(K) > 0. .AND. XM(J) < 0.) THEN IF (J .NE. K .AND. MLIQ(K) > 0. .AND. XM(J) < -0.1) THEN

There are some pretty big impacts from changing the excess cold threshold. For the first 34 history writes at two-hourly frequency, only TSLB below z=2 in glacier grid cells differs between the two simulations (warmer in SOL2 at z=2 and cooler in z = 3 & 4). At write 35, TSLB(z=1) differs at a single grid point, and therefore so do TGB, GRDFLX and other surface energy vars. [image: Screenshot 2020-05-25 at 14 25 51] https://user-images.githubusercontent.com/51397205/82812769-01c11180-9e94-11ea-8ddb-e09f49d58216.jpg

By the next write two hours later, TGB amongst many vars is impacted nearly domain wide, with anomalies of 1-2K after a few more hours. [image: Screenshot 2020-05-25 at 14 26 08] https://user-images.githubusercontent.com/51397205/82812912-4c428e00-9e94-11ea-9c20-34b483613e5c.jpg [image: Screenshot 2020-05-25 at 14 26 26] https://user-images.githubusercontent.com/51397205/82813192-de4a9680-9e94-11ea-9cfd-855dfa153ac0.jpg

So the code changes have a large impact and it would be great to have your input on which change, if any, would be best to contribute.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/wrf-model/WRF/issues/1185#issuecomment-633554215, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIZ77EGV5YUNNEUCCVDW43RTJR2HANCNFSM4MVSLZNA .

[MadleneP]

Hello!

I found a similar issue in my simulations over the Alps, although I did not change the albedo. I am using WRF v4.3 to downscale re-analysis data (20CRv3) to a 2 km resolution in a nested domain over an area over the Alps in hourly resolution. As land surface model, I am using NoahMP with opt_gla=1. My simulation starts on 15-08-1849 and it is crashing on 26-09-1849 with the following error:

ERRENG = 1.16767883E-02 108 62 59.4008 157.194911049.5908-8391.0439-2756.3530 -------------- FATAL CALLED --------------- FATAL CALLED FROM FILE: LINE: 3004 Energy budget problem in NOAHMP GLACIER

The grid cell in which the model crashes has a land category of snow and ice and an altitude of 2362 m. After a closer look, and after outputing in 1 minute resolution, I found highly unrealistic soil temperatures as well, of up to +1209°C in the upper level. The other soil levels are affected as well: 2nd layer has temperatures of up to 512°C, 3rd layer up to 56°C, and 4th layer up to 3.3°C. This deviation from the normal range can be seen at three different points in time. However, the model seems to recover after each of these jumps up until roughly two model-days (on 23-09-1849 at 22:00) before the crash on 26-09-1849 at 05:05.

The high soil temperatures coincide with highly unrealistic negative sensible and latent heat fluxes of close to -1500 W/m² and positive ground heat fluxes of up to +2000 W/m².

However, starting 5 minutes before the crash, the sensible heat flux jumps from –990 W/m² to +46155 W/m² and the ground heat flux to -7444 W/m², at which point the model crashes with the energy budget problem in the glacier module of the NoahMP.

There are also jumps in the skin temperature from 0.2° to 75°C from one minute to the next, as well as in the 2 m temperature from 0.9° to 37°C. A sudden drop in the 2 m relative humidity from 99.7 to 8.9%. A jump in the surface upward longwave flux from 315 to 817 W/m² 1 minute after the jump seen in other variables. The skin temperature looks strange already 2 model-days before the crash, when it doesn't show any day-night cycle anymore. Other variables, such as 10 m wind speed and surface pressure are affected as well. So the model seems to behave at specific grid cells in an erratic way.

I tested also gla_opt=2. The model did not crash anymore at 26-09-1849. However, the soil temperatures in the upper level do not look as I expect them to look. The diurnal cycle looks too strong for a glacierized grid cell, I assume, and there are times when temperatures are completely flat. The grid cell where the model used to crash when using opt_gla=1 is represented by the red line.

I would be very grateful if you could give some suggestions on how to overcome this issue. Is there a way to switch off the glacier module, assuming that the source of the error is indeed the glacier module? Are there any chances to fix this problem in the next WRF version?

Many thanks in advance, Madlene