A new module calculateWaterBalance is proposed in this branch. This module has two tasks:

• Task 1: Calculate water balance and check closure
• Task 2 (optional): Correct water balance error

Task 1

Water balance concept:
total inflow = total outflow + residual
total inflow = precipitation
total outflow = evapotranspiration + runoff + bottom boundary flux + delta storage
bottom boundary flux = flux at bottom boundary (in case of no groundwater coupling)
bottom boundary flux = recharge (in case of groundwater coupling)
delta storage = sum[ (theta(t+1) - theta(t)) / timestep * soil thickness]
residual = total inflow - total outflow
error = residual / mean([total inflow + total outflow]) * 100%

Task 1 is done in lines 54-90

Task 2 (optional)

The user can decide to activate this option by adding one line in the config file closeWaterBalance=1. This line is read by the function read_config. If that line is not there, Task 2 is not activated. If the water balance error is large, a correction needs to be made.

Correction concept
To close the water balance, one of the water fluxes (ET, runoff, bottom boundary flux, or delta storage) must be adjusted. Since all fluxes depend on soil moisture, the delta storage is chosen for correction. Why? Because 1) the adjustment is distributed across the entire soil moisture profile, so each soil layer undergoes only a minor modification, minimizing potential disturbances, and 2) the matric potential is a state variable which is solved in the Richards equation adopted in STEMMUS, then this calculated matric potential is converted to soil moisture using the Van Genuchten equation, so the calcualted soil moisture has some uncertainty due to the parameters of Van Genuchten. This approach prevents large changes in other fluxes that could occur if they were corrected directly. The correction in the soil moisture profile is done through a loop of 30 iterations, which does the following:

• Calculate the correction in the soil moisture profile from the residual (line
• Apply a constraint to that correction: max change in soil moisture (to avoid dramatic changes) (line 122)
• Add the correction value to the soil moisture profile (lines 129-132)
• Ensure corrected soil moisture values remain within residual and saturated water content limits (line 135)
• Recalculate the correction value of soil moisture after checking residual and saturation bounds (line 138)
• Recalculate corrected delta storage (based on the corrected soil moisture profile) (lines 141-143)
• Calculate corrected total inflow, total outflow, residual and error (lines 145-151)

The water balance module is placed in the main function within the loop of STEMMUS (here) and also runs after the last update of soil moisture here

Outputs of all fluxes (lines 159-175) are written to a new CSV file named waterBalance (whether the correction is switched on or off)

Additional related changes

• Calculation of Dunnian runoff is moved from loadForcingData to calculateBoundaryConditions to account properly for the snowmelt as part of total precipitation (which was discussed in a former issue and added to the code -> here). Now the two runoff components (Hortoian and Dunnian) are in one function (calculateBoundaryConditions), which i think is better