Water Balance Terms in 1D Noah-MP

[mcthreems]

Hello, I am trying to understand how the water balance terms work in the 1D (single point) Noah-MP model. I have identified 11 output variables that seem like they should represent the full water budget of the column for a run using the default model settings, these being: 'RAINRATE', 'ECAN', 'ETRAN', 'EDIR', 'UGDRNOFF', 'SFCRNOFF', 'CANICE', 'CANLIQ', 'SNEQV', 'SOIL_M'. I converted each output into mm/day (including converting to a flux from an accumulation or storage term for the latter 6), but there appears to be a missing term(s) as the budget from these outputs is not 0. Are there additional outputs I need to include in this calculation?

I also want to ask whether there actually is a true water balance in this model that I can calculate using the output variables?

[cenlinhe]

Here is the source code where the model checks water balance/residual every model timestep (if the residual is large, the model will crash): https://github.com/NCAR/noahmp/blob/master/src/BalanceErrorCheckMod.F90#L58
The 11 output variables you identified are the major water balance terms used if you do not activate irrigation, tile drainage etc. (special processes). Please also double check the unit of each term (particularly the runoff). For example, 'UGDRNOFF', 'SFCRNOFF' are accumulated quantities, whereas 'RAINRATE', 'ECAN', 'ETRAN', 'EDIR', are flux terms. 'CANICE', 'CANLIQ', 'SNEQV', 'SOIL_M' are state quantities.

[mcthreems]

Thanks for directing me to this code section. I reproduced the method from this code to calculate the water balance but my results were the same. I did make sure to be careful of the different output units in the calculation. The imbalance I am seeing is very close in magnitude to the precipitation which is strange to me even just on its own. My outputs are hourly while the Noah timestep is 30 minutes, would this create an issue when trying to calculate the water balance?

[cenlinhe]

This is indeed weird. hourly output should be fine. Could you please show your accumulated timeseries of each budget term and residual?

[mcthreems]

Here are the accumulations over each day for all components. I calculated the residual for this figure using the method from the code section you sent. The simulation itself uses the default single-point Noah-MP setup with NLDAS2 forcing.

[cenlinhe]

Your residual curve indeed looks the same as precipitation curve. Right now, I am not quite sure why. I will try to discuss with other colleagues.
@tslin2 @RAbolafiaRosenzweig @CharlesZheZhang When you have time, could you please take a look at this issue?

[RAbolafiaRosenzweig]

Can you please share the script you use to calculate the residual error? From a single pixel can you please also share a small sample (e.g., 24-50) time steps for a table showing the hourly values of each WB term and the calculated residual error?

[mcthreems]

Here is the table of sample timestep data plus the code used to calculate the values within:

rawdat = '/zdata/data/matt_data/hrldas/hrldas/run/inputs/exp_setups/NLDAS2_38.0_-100.0_dx0.5_dy0.5/out_NLDAS2_gw3_conus_def_nrdef/2000010101.LDASOUT_DOMAIN1.loop0001'
varnames = ['SOIL_M','RAINRATE','ECAN','ETRAN','EDIR','UGDRNOFF','SFCRNOFF','CANICE','CANLIQ','SNEQV']

tabarr = []
with netcdf(rawdat) as ncf:
for varname in varnames:
idat = np.squeeze(np.asarray(ncf.variables[varname]))
if varname not in ['UGDRNOFF','SFCRNOFF']:
idat_sub = np.copy(idat[1000:1050])
else:
idat_sub = np.copy(idat[999:1050])
accvec = np.zeros((50,))
if varname == 'SOIL_M':
idat_sub[:,0] = idat_sub[:,0]0.11000. #convert volumetric moisture to mm equivalent using layer thicknesses
idat_sub[:,1] = idat_sub[:,1]0.31000.
idat_sub[:,2] = idat_sub[:,2]0.61000.
idat_sub[:,3] = idat_sub[:,3]*1.*1000.
accvec[1:] = np.cumsum(np.diff(np.sum(idat_sub,axis=-1))) #sum soil layer values, then take difference before doing cumsum since it's storage
elif varname in ['CANICE','CANLIQ','SNEQV']:
accvec[1:] = np.cumsum(np.diff(idat_sub)) #take difference before doing cumsum since it's storage
elif varname == 'RAINRATE':
accvec[:] = np.cumsum(idat_sub) #precip is mm/timestep, can do cumsum directly
elif varname in ['ECAN','ETRAN','EDIR']:
accvec[:] = np.cumsum(60.*idat_sub) #evap is mm/s, convert to mm/hr to match timestep before doing cumsum
else:
accvec[:] = np.copy(idat_sub[1:] - idat_sub[0]) #runoff is already an accumulation, so just remove accumulation from before subset

    tabarr.append(np.copy(accvec))

resvec = tabarr[0] + tabarr[7] + tabarr[8] + tabarr[9] - (tabarr[1] - tabarr[2] - tabarr[3] - tabarr[4] - tabarr[5] - tabarr[6]) #calculate residual
tabarr.append(np.copy(resvec))
tabarr = np.asarray(tabarr)

[RAbolafiaRosenzweig]

Thanks for sharing here. It is hard to follow on the GitHub interface. Can you please attach your script file, netcdf data file, and summary table (e.g., as an excel file)?

[mcthreems]

I've attached the script file and csv output file, but the nc data file was too large. I have made that file available at this link, let me know if it doesn't work: https://www.icloud.com/iclouddrive/0e4sCD-iX1F44LvhZZSJ45zRg#2000010101.LDASOUT%5FDOMAIN1.loop0001
water_bal_test_csv.csv
water_bal_test_table.py

[RAbolafiaRosenzweig]

It looks like the residual error is primarily from inconsistent units in the precipitation. Namely, the output file you provided notes rainrate in units of mm/timestep. Because you are running the model at a 30 minute time step, and the other fluxes for your WB error calculation are hourly, you will need to multiply PRCP by 2.

[mcthreems]

Yes, this seems to be it. After applying this correction I do still get some accumulating residual over time, but much less (now only about 40mm total accumulated residual over the 25-year run). I'll mark this as the answer for my main question, but I am also curious about whether the balance terms change when using the TOPMODEL with aquifer option (runoff option #1 in the namelist). My understanding for this scheme is the subsurface runoff term is taken out of the aquifer, and water first drains from the lowest soil layer into the aquifer. Is this understanding correct or is it missing something?

[cenlinhe]

The aquifer water is accounted for in TOPMODEL case. The remaining residual may be due to the soil moisture solver inaccuracy and/or too coarse vertical soil layer configuration. We will discuss with key model developers to see how we could improve the solver to further reduce residual.

[tslin2]

you may try to use pFUnit to discover the water imbalance
https://github.com/Goddard-Fortran-Ecosystem/pFUnit

[RAbolafiaRosenzweig]

The remaining residual water balance error appears to stem from a mismatch between the model time step (likely 30 minutes) and the output time step (hourly). The output variables used for the water balance calculation represent instantaneous values rather than time-accumulated or averaged quantities. As a result, at some output times the instantaneous fluxes at the end of the output interval do not fully correspond to the change in storage over that interval. Future model updates will aim address this by providing accumulated flux outputs. In the interim, achieving precise water balance closure from outputs at all times requires that the output time step match the native model time step.