Hourly temperature calcs from daily max/min data

[MitchHeineman]

I noticed large temperature jumps in reviewing model performance during winter conditions with snow processes active. For models using daily max/min temperatures, I found simulated hourly temperatures changing by as much as 20°F at midnight. I compared the algorithm in settemp in climate.c with the corresponding SWMM4h code in wshed.for and found that they match. However, I also observed that the SWMM4h code was changed in 2005, when SWMM5 had already been released, and that the SWMM4h code does not conform with the algorithm's description in Appendix II of the SWMM4 manual. Both the SWMM3 and SWMM4 manuals state that “The interpolation is performed, using three different periods: 1) between the maximum of the previous day and the minimum of the present, 2) between the minimum and maximum of the present, and 3) between the maximum of the present and minimum of the following.” Temperature calculations for a given day should thus depend on the prior day’s maximum temperature, the current day’s max and min, and the following day’s minimum. However, the SWMM4h/SWMM5 algorithm depends on the prior day’s temperature range (not its maximum) and the current day’s max and min. It does not consider the following day at all. Ironically, the changes were meant to remediate the issue I’ve observed, as item 131 in 44guchng.txt says “Fix Runoff air temperature interpolation scheme in Sub. WSHED for snowmelt. Old scheme caused discontinuity at midnight. RED, 7/1/05.”

While hourly temperature datasets are more widely available today than in the past, it remains desirable for SWMM to include a temperature algorithm based on daily max/min temperatures. SWMM also optionally uses daily temperature as input for evaporation calcs but this issue is largely irrelevant there, as evaporation processes function on longer timescales than snow. I propose to add an option for calculating temperature based on the algorithm described in the SWMM3/4 manuals. I can draft code if that would help move the process along. Otherwise, I present here VBA adapted from code I wrote years ago that is a slight simplification of the older algorithm – whereas the SWMM methods place minimum temperature at sunrise and maximum three hours before sunset, this code places the minimum at 6 AM and the maximum at 3 PM.

The figure presented here compares temperatures computed in SWMM5 with those for the VBA algorithm using max/min data from Lawrence, Massachusetts (USW00094723). While “my” algorithm yields temperatures outside the range of the reported max/min late on Feb. 2, its changes are much smoother and overall average temperatures conform just as well with the average of the max/min dataset. The SWMM5 temps drop from 30°F to 12°F at midnight on Feb. 3, jump from 2°F to 16°F at midnight on Feb. 5, and generally don't look pretty. The SWMM4 method remains imperfect but yields better overall results.

    1 2023    2    1   30   16  
    1 2023    2    2   38   22  
    1 2023    2    3   33  -10  
    1 2023    2    4   16  -12  
    1 2023    2    5   51   16  
    1 2023    2    6   49   30  

Function CurrTemp(hiermax As Integer, todaymax As Integer, todaymin As Integer, morgenmin As Integer, hour As Integer) As Double
' Method in SWMM4 manual, Appendix II
' M. Heineman 1996 in xBASE; ported to Fortran 2004, VBA 2024
' Assume Min temp at 6 AM, Max temp at 3 PM
' Hour ranges from 1 to 24
Dim pi As Double
pi = 2# * WorksheetFunction.Acos(0#)
Select Case (hour)
    Case 1 To 5
        CurrTemp = hiermax + (todaymin - hiermax) * (1 - Cos(pi * (hour + 9) / 15)) * 0.5
    Case 6 To 15
        CurrTemp = todaymin + (todaymax - todaymin) * (1 - Cos(pi * (hour - 6) / 9)) * 0.5
    Case Else
        CurrTemp = todaymax + (morgenmin - todaymax) * (1 - Cos(pi * (hour - 15) / 15)) * 0.5
End Select
End Function

[cbuahin]

@MitchHeineman, this fix is simple enough that we may be able to add it to the v5.3.0 release that we are working on. Let me look into it further. Thanks for raising the issue and proposing the solution.

[cbuahin]

I should add that it will be good to understand how large the bound exceedances can get and their implications for the snowmelt processes. I don't think they will be large, but I plan to check to see.

[MitchHeineman]

To make sure I'm being clear - I'm not disputing the methodology whereby min temp is at sunrise and max 3 h before sunset. My concern is how temps are calculated from afternoon maximum until the following sunrise.

[MitchHeineman]

Further checking shows that I was wrong to state that SWMM44h and SWMM5 use the previous day's range. They instead use the prior day's max temp, conforming with the description in the SWMM 3 through 5 manuals. However, the issue that the subsequent day's minimum temperature is ignored has been present since at least the 1990s, based on this 1994 SWMM4 code:

C******** COMPUTE HOURLY TEMPERATURE USING SINSUOIDAL INTERPOLATION.
C=======================================================================
                     DECL  = 0.40928*COS(0.017202*(172.0-DAYNO))
                     HRANG = 3.8197*ACOS(-TAN(DECL)*ANGLAT)
                     HRSR  = 12.0 - HRANG + DTLONG + 0.5
                     HRSS  = 12.0 + HRANG + DTLONG - 2.5
                     DHRDY =  HRSR - HRSS
                     DYDIF =  24.0 + HRSR - HRSS
                     HRDAY = (HRSR + HRSS)/2.0
                     TAVE  = (TMAX + TMIN) / 2.0
                     TRNG  = (TMAX - TMIN) / 2.0
                     TRNG1 = (TMAX1 - TMIN)
                     ENDIF
C#### RED (WCH), 1/31/94.  SHOULD HAVE ENDIF HERE, OTHERWISE
C####   ONLY ONE AIR TEMPERATURE PER DAY.  BUT ALSO NEED ADDITIONAL
C####   IF-STMT HERE TO PERFORM FOLLOWING ONLY FOR ISNOW > 1.
      IF(ISNOW.GT.1) THEN
                     HR    = TIMDAY/3600.0
                     IF(HR.LT.HRSR) TA = TMIN +
     +                  TRNG1 * SIN(3.14159 / DYDIF * (HRSR - HR))
                     IF(HR.GE.HRSR.AND.HR.LE.HRSS) TA = TAVE
     +                 + TRNG * SIN(3.14159 / DHRDY * (HRDAY - HR))
                     IF(HR.GT.HRSS) TA  = TMAX -
     +                  TRNG * SIN(3.14159 / DYDIF *(HR - HRSS))
                     ENDIF