Critical regression in NOAH-MP snow physics between v4.2 and v4.7 - persistent cold pools

[meteoadriatic]

Problem Description: WRF v4.3+ with NOAH-MP LSM produces widespread unrealistic surface temperature cold pools (5-8K below surroundings) that do not appear in v4.2 using identical configuration and initial conditions (as already others found to be the case - see here). The cold pools correlate with areas of snow cover, including trace snow cover (<1mm snow water equivalent).

This is already reported several times, for example here, here, here and especially already mentioned thread. I made some research and want to report my findings, as I think this is important to resolve.

Key Findings:

1. Version comparison confirms regression:
   • v4.2.2: Trace snow melts by mid-morning; surface temperatures reach 282-285K at noon; spatially smooth temperature fields
   • v4.7: Trace snow persists through noon; surface temperatures stuck at 275-278K; distinct cold pool patterns remain
2. Bare soil test isolates vegetation coupling:
   • Setting VEGFRA=0 in wrfinput eliminates cold pools in v4.7
   • Problem seem to be specific to vegetated grid points
   • Suggests issue in vegetation-ground heat flux coupling with (trace?) snow

Diagnostic evidence at affected example gridpoint:

• SNEQV persists at 0.75mm through noon in v4.7 (vs 0.0mm in v4.2.2), in sunny autumn day
• Ground heat flux under vegetation (GHV) up to 100 W/m² and more into cold soil
• STC(1) stuck up to 5K or even more below surface temperature
• Extremely slow surface warming rate 0.12 K/hour (less than half the rate of bare soil)
  

Critical diagnostic finding: Thermodynamically impossible temperature profile: Diagnostics show STC(1) becomes colder than both the surface and the layer below it, one random example at affected gridpoint with trace snow amount:
TG (surface) = 280K
STC(1) (top soil) = 276K ← coldest layer - impossible
STC(2) (deeper) = 279K

Please investigate. It's not just "slower warming", compared to v4.2 (or NOAH!), it's fundamentally broken physics creating impossible temperature gradients.

I can provide assistance as much as possible, but alone was not able to find good solution.

Workaround only: reduce CBIOM in MPTABLE to 0.002 (seems to be correct value instead of 0.02), gets better, but does not resolve the issue completely.
Thank you,
Ivan

 

[meteoadriatic]

Few plots from 24h simulation, initialized at 2025-10-03 at 12 UTC using ECMWF input data.

2025-10-04 12 UTC (24h lead time). 14 hours local time, no cloud cover:

TG (range 275-305 K)


TSLB(0) (range the same as above) - note, the cold pool here is generally colder than TG and TSLB(1)


TSLB(1) (range the same as above) - this one is warmer than TSLB(0)


SNOWH (range: 0-0.18 m)


SWDOWN (range: 16-708 W/m2)



I can provide plots with v4.2 if needed, or whatever else you request.

 

[meteoadriatic]

Smoking gun (however, present in v4.2, too, where this issue is not apparent):

! adjusting snow surface temperature
     IF(OPT_STC == 2) THEN
      IF (SNOWH > 0.05 .AND. TG > TFRZ) THEN
        TGV = TFRZ
        TGB = TFRZ
          IF (VEG .AND. FVEG > 0) THEN
             TG    = FVEG * TGV       + (1.0 - FVEG) * TGB
             TS    = FVEG * TV        + (1.0 - FVEG) * TGB
          ELSE
             TG    = TGB
             TS    = TGB
          END IF
      END IF
     END IF

Similar for OPT_STC ==1 and OPT_STC == 3... why do we do this? OK This is valid only for SNOWH > 0.05 m, but I saw now change if I modify this threshold. But even more suspecting, in PHASECHANGE subroutine:

    DO J = ISNOW+1,NSOIL
         IF (MICE(J) > 0.0 .AND. STC(J) >= TFRZ) THEN  !melting
             IMELT(J) = 1
         ENDIF
         IF (MLIQ(J) > SUPERCOOL(J) .AND. STC(J) < TFRZ) THEN !refreezing
             IMELT(J) = 2
         ENDIF

         ! If snow exists, but its thickness is not enough to create a layer
         IF (ISNOW == 0 .AND. SNEQV > 0.0 .AND. J == 1) THEN
             IF (STC(J) >= TFRZ) THEN
                IMELT(J) = 1
             ENDIF
         ENDIF
    ENDDO

! Calculate the energy surplus and loss for melting and freezing

    DO J = ISNOW+1,NSOIL
         IF (IMELT(J) > 0) THEN
             HM(J) = (STC(J)-TFRZ)/FACT(J)
             STC(J) = TFRZ
         ENDIF

If I understand it well, this forces ground temperature to 0°C if snow is melting, for the traces of snow without layer (SNEQV > 0.0 .and. ISNOW == 0). Why we are doing that? I do understand the motivation from physical standpoint of melting ice, but I'm not sure that implemenation of this is correct as here we are resetting ground temperature, not snow temperature. I think this is a big problem, but I'm not convinced that the problem of cold pools is caused (only) by this.

There are more suspicios parts of the code, but I will stop here, waiting for someone's response.
Ivan

 

[Ming Chen]

Hi Ivan,

Thank you for the detailed description of the snow issue in NoahMP. In your test case, which data did you use to initialize WRF? We know that unrealistic snow depth in ERA5 may lead to severe cold bias later in the model simulation.

While the initial condition could be an issue, I agree with you that there might exist systemic bias in the land model. I will talk to our expert Cenlin and hopefully we can find the reason behind the cold bias. We will get back to you once we know for sure what is wrong.

 

[cenlinhe]

Smoking gun (however, present in v4.2, too, where this issue is not apparent):

! adjusting snow surface temperature
     IF(OPT_STC == 2) THEN
      IF (SNOWH > 0.05 .AND. TG > TFRZ) THEN
        TGV = TFRZ
        TGB = TFRZ
          IF (VEG .AND. FVEG > 0) THEN
             TG    = FVEG * TGV       + (1.0 - FVEG) * TGB
             TS    = FVEG * TV        + (1.0 - FVEG) * TGB
          ELSE
             TG    = TGB
             TS    = TGB
          END IF
      END IF
     END IF

Similar for OPT_STC ==1 and OPT_STC == 3... why do we do this? OK This is valid only for SNOWH > 0.05 m, but I saw now change if I modify this threshold. But even more suspecting, in PHASECHANGE subroutine:

    DO J = ISNOW+1,NSOIL
         IF (MICE(J) > 0.0 .AND. STC(J) >= TFRZ) THEN  !melting
             IMELT(J) = 1
         ENDIF
         IF (MLIQ(J) > SUPERCOOL(J) .AND. STC(J) < TFRZ) THEN !refreezing
             IMELT(J) = 2
         ENDIF

         ! If snow exists, but its thickness is not enough to create a layer
         IF (ISNOW == 0 .AND. SNEQV > 0.0 .AND. J == 1) THEN
             IF (STC(J) >= TFRZ) THEN
                IMELT(J) = 1
             ENDIF
         ENDIF
    ENDDO

! Calculate the energy surplus and loss for melting and freezing

    DO J = ISNOW+1,NSOIL
         IF (IMELT(J) > 0) THEN
             HM(J) = (STC(J)-TFRZ)/FACT(J)
             STC(J) = TFRZ
         ENDIF

If I understand it well, this forces ground temperature to 0°C if snow is melting, for the traces of snow without layer (SNEQV > 0.0 .and. ISNOW == 0). Why we are doing that? I do understand the motivation from physical standpoint of melting ice, but I'm not sure that implemenation of this is correct as here we are resetting ground temperature, not snow temperature. I think this is a big problem, but I'm not convinced that the problem of cold pools is caused (only) by this.

There are more suspicios parts of the code, but I will stop here, waiting for someone's response.
Ivan

This treatment is to ensure that the temperature is not above freezing point when there is still snow on the ground. This treatment is not the cause for the cold bias.

 

[cenlinhe]

The cold bias is a long-standing issue in WRF/Noah-MP. We have a paper in review that systematically investigated this issue and found that overestimated snow cover in the model is only contributing to 20% of the cold bias, while there must be some other dominant processes causing the cold bias. This is still an unresolved issue. The cause could potentially be land-atmosphere coupling (might be vegetation-related) or other uncertainties in atmospheric processes and/or snow and soil thermal conductivity. Here are some possible fixes related to land model to reduce the cold bias:
1. reduce CBIOM to a very small values (e.g., 10^-5);
2. increase snow thermal conductivity (TKSNO; you can use the constant formulation or apply a multiplier at this code line): noahmp/src/module_sf_noahmplsm.F at e5c0859874407859936739e8be8741f9aed369ee · NCAR/noahmp
3. increase SCFFAC and/or MFSNO parameter values in MPTABLE for specific vegetation types in your modeling domain.

 

[meteoadriatic]

Hi Ivan,

Thank you for the detailed description of the snow issue in NoahMP. In your test case, which data did you use to initialize WRF? We know that unrealistic snow depth in ERA5 may lead to severe cold bias later in the model simulation.

While the initial condition could be an issue, I agree with you that there might exist systemic bias in the land model. I will talk to our expert Cenlin and hopefully we can find the reason behind the cold bias. We will get back to you once we know for sure what is wrong.

Hello, thank you. In this particular case initial snow height = 0, and those 0-18 cm were fallen during the simulation period. The initial soil temperature that came from forcing dataset is correct/smooth, but during the simulation it went completely wrong for the areas where snow got accumulated. In other areas temperature profiles stay correct to the end of simulation.

The cold bias is a long-standing issue in WRF/Noah-MP. We have a paper in review that systematically investigated this issue and found that overestimated snow cover in the model is only contributing to 20% of the cold bias, while there must be some other dominant processes causing the cold bias. This is still an unresolved issue. The cause could potentially be land-atmosphere coupling (might be vegetation-related) or other uncertainties in atmospheric processes and/or snow and soil thermal conductivity. Here are some possible fixes related to land model to reduce the cold bias:
1. reduce CBIOM to a very small values (e.g., 10^-5);
2. increase snow thermal conductivity (TKSNO; you can use the constant formulation or apply a multiplier at this code line): noahmp/src/module_sf_noahmplsm.F at e5c0859874407859936739e8be8741f9aed369ee · NCAR/noahmp
3. increase SCFFAC and/or MFSNO parameter values in MPTABLE for specific vegetation types in your modeling domain.

Thank you, but the most interesting part is that WRF 4.2 does not have this problem. Something changed in 4.3 and I can't figure out what. I do understand that there are ways to force the model toward better solution, but those seem to be just a bandaid, not a fix. I really hope that there is a way to find what happened in v4.3 and to code a proper fix.

Ivan

 

[Ming Chen]

Hi Ivan,

I just double checked all the released version of WRF from V4.2 to V4.7. Except for a few bug fixes in Noah-MP for snow combination, vegetation fraction scaling, and urban ground heat flux sign, I couldn't find any significant changes.

Please confirm that the large cold bias occurs in WRFv4.3+, but results of WRFv4.2 is fine. Would you please run the same case using WRFV4.2.1 and WRFv4.2.2, and let me know whether these two versions also give you reasonable results? By these tests, I would like to know for sure that, starting from which version the large cold bias occurred. This information will be helpful for us to trace back and find possible bugs.

At present the cold bias issue you reported remains a mystery for us. I am actually not 1000% sure it is related to snow physics in NoahMP. There might have some other bugs in the model we are not aware of. We will look into this issue and get back to you if we figure out what is wrong.

At the same time, please keep us updated if you have any findings. Thanks in advance.

 

[meteoadriatic]

First of all, I need to apologize for the incorrect statement made above in the introductory post. I blindly trusted another forum member's post where he stated that the issue appeared with WRF 4.3. It is not the case. The issue appeared with WRF 4.4 and after, and v4.4 has seen a lot of NOAH-MP changes!

Here is an example of a result matrix of the same sinoptic situation as in my introductory posts - same date, same hour. The only difference is that this time the domain is shifted a bit toward the south-east, to center over the problem area. All of the runs in a matrix below, are done with the exactly the same namelists.

Interestingly, the issue seems to be somewhere at or above the surface (related to the vegatation cover?) as we see most unrealistic anomalies in TSK and T2 fields. Moving deeper into the soil [TG, TSLB(0), as well as even deeper down (not shown)], the cold pool gets less pronounced, relative to the unaffected versions (4.2.2, 4.3, and I believe, before that).

No wonder we were not able to find the cause of issue in v4.3. I will dig down into changes made in v4.4 later, but given the excessive amount of those in that version, I'm not confident that I will be able to detect the problem on my own. I will appreciate if someone with more NOAH-MP knowledge (@cenlinhe) can take a look into issue, with regard to these new important findings.

Also please let me know if you want any other tests to be ran here.

Thank you
Ivan

 

[haiqingsong]

First of all, I need to apologize for the incorrect statement made above in the introductory post. I blindly trusted another forum member's post where he stated that the issue appeared with WRF 4.3. It is not the case. The issue appeared with WRF 4.4 and after, and v4.4 has seen a lot of NOAH-MP changes!

Here is an example of a result matrix of the same sinoptic situation as in my introductory posts - same date, same hour. The only difference is that this time the domain is shifted a bit toward the south-east, to center over the problem area. All of the runs in a matrix below, are done with the exactly the same namelists.

Interestingly, the issue seems to be somewhere at or above the surface (related to the vegatation cover?) as we see most unrealistic anomalies in TSK and T2 fields. Moving deeper into the soil [TG, TSLB(0), as well as even deeper down (not shown)], the cold pool gets less pronounced, relative to the unaffected versions (4.2.2, 4.3, and I believe, before that).

No wonder we were not able to find the cause of issue in v4.3. I will dig down into changes made in v4.4 later, but given the excessive amount of those in that version, I'm not confident that I will be able to detect the problem on my own. I will appreciate if someone with more NOAH-MP knowledge (@cenlinhe) can take a look into issue, with regard to these new important findings.

Also please let me know if you want any other tests to be ran here.

Thank you
Ivan

View attachment 19333

Many changes in WRF4.4 about NOAHMP：Sync with NoahMP Github version with all NoahMP updates since v4.3 (#… · wrf-model/WRF@a82ce24
Canopy heat storage in canopy temperature calculation was added.

 

[cenlinhe]

Yes, we added a canopy heat storage treatment in WRF4.4, which may be the cause for the cold bias occuring since WRF4.4. Adjusting CBIOM parameter to a very small number is essentially to reduce or remove the canopy heat storage effect. This effect currently has a large uncertainty without direct observational constraints as described in other Q&A threads in WRF forum (mentioned by Ivan above). So I would suggest reducing the value of CBIOM to a very small value to remove this uncertainty canopy heat storage effect. In our newly offline Noah-MP version, we adjust the default CBIOM value down already based on many tests showing temperature bias with the original CBIOM value (0.02).

 

[vaakos]

Dear all,

I am writing here because this issue appears to be identical to that in this other forum post.

In that post, I wrote the following:

Dear all,

I'm writing a follow-up on this topic because the issue still persists in the latest model version, WRFV4.7.0.

The problem remains the formation of persistent, steady cold patches that appear near the surface in association with extensive cloud cover. These patches emerge after a few days of simulation.

I understand that WRF is primarily designed for short-term numerical weather prediction. However, this is not simply a case of the model developing biases over long-term integration. Rather, the phenomenon appears completely unphysical and may indicate a potential code error. Given the widespread use of WRF as a regional climate model, I believe this issue warrants thorough investigation.

To support this, I’ve uploaded output from the first few days of my simulation. In the hourly_output_d02_* files, the variables T2 and TSK clearly show unphysical values (e.g., over the western parts of the domain) as early as the second day of integration.

I’ve also included the namelist.input and met_em* files to facilitate reproducibility. The initial and boundary conditions are from ERA5.

Please follow this link to access the files.

I hope we can eventually find a solution to this.

Best regards,
Akos

I have re-uploaded the files referenced in that post to Box storage (filename: WRFV4.7.0_cold_patch_issue.tar.gz).

I sincerely hope this issue will be investigated and resolved soon.

Best regards,
Akos

 

[meteoadriatic]

OK thank you all, after detailed comparative testing it seems that introducing canopy heat storage effects in v4.4 did the regression in simulation accuracy.

Indeed, going very low with CBIOM parameter effectively reverses these changes. Going 10x lower as suggested in few places, that is from 0.02 to 0.002 is not enough to fully reverse the effect, but 0.001 almost is. I don't know what would the physically correct thing to put there be, but 0.001 seem to work reasonably.

Here is direct comparison, T2, range 273-293 K fixed on all plots.
Left: v4.3
Middle: v4.7.1 with CBIOM = 0.02 (currently default value)
Right: v4.7.1 with CBIOM = 0.001

Thank you all who participated. Should we create a pull request to set a better default value for next model update?

Ivan

 

[vaakos]

I do not currently have the data, but I remember seeing the cold patch issue in my test runs at convection-permitting resolutions, even when setting the CBIOM value as low as 0.0001, using specific physical parameterizations (e.g., microphysics and icloud). These cold patches were, however, very limited in extent and spatially isolated. Therefore, I’m not sure that reducing the CBIOM value is a generally applicable solution for all model configurations, and I think a thorough, physically consistent investigation should be carried out.

 

[cenlinhe]

This cold bias problem is not only a land model issue. It is also related to atmospheric parameterizations. It is still an unresolved issue and many groups are investigating this issue right now from different aspects. For now, we will submit a WRF PR to update the default CBIOM value as a temporary fix due to its large uncertainty. As more and more research is available related to this issue, we will add more fixes/solutions. Thank you all for looking into this issue.

 

[haiqingsong]

I do not currently have the data, but I remember seeing the cold patch issue in my test runs at convection-permitting resolutions, even when setting the CBIOM value as low as 0.0001, using specific physical parameterizations (e.g., microphysics and icloud). These cold patches were, however, very limited in extent and spatially isolated. Therefore, I’m not sure that reducing the CBIOM value is a generally applicable solution for all model configurations, and I think a thorough, physically consistent investigation should be carried out.

This is incredible.
Which microphysics scheme do you use?

 

[vaakos]

There were cases where switching from mp_physics=8 to mp_physics=10 brought back the cold patch problem (WRFV4.7.0).

 

[Ming Chen]

Dear All,
Thank you so much for the detailed information regarding the cold bias issue starting from WRFV4.4. We will look into this problem and keep you updated of any progress.