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