-
Notifications
You must be signed in to change notification settings - Fork 680
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.
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
Unphysical subsurface temperatures simulated in glacierized grid cells using Noah-MP #1185
Comments
@emilycollier @weiwangncar @dudhia |
Agree with Dave. Let's add @barlage in case he is not following this. |
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 ***@***.***> 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
<#1185 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AEIZ77EGV5YUNNEUCCVDW43RTJR2HANCNFSM4MVSLZNA>
.
|
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
|
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 PHASECHANGE_GLACIER 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
The text was updated successfully, but these errors were encountered: