error in calculating water balance using hperwat function

[AtrCheema]

I am using hperwat function to simulate processes in a simple PERLND and calculate water balance using following lines of code

import pandas as pd
import numpy as np

np.random.seed(seed=77)

from perwat import pwater

tindex = pd.date_range('20110101', '20111210', freq='H')

general = {
    'sim_len': len(tindex),
    'tindex': tindex,
    'sim_delt':60  # number of minutes in time step
}

ui = {
    'CSNOFG': 0,
    'RTOPFG': 1,  #
    'UZFG': 0,   # algorithm used with `0` causes a lot of surface runoff (thus higher pero) at the expense of uzet
    'VCSFG': 0,
    'VLEFG': 0,
    'IFFCFG': 0,
    'IFRDFG': 0,

    'FOREST': 0.0,
    'LZSN': 0.25,
    'INFILT': 1.0,  #  index to the infiltration capacity of the soil.
    'LSUR': 1.0,
    'SLSUR': 1.5,
    'KVARY': 1.0, # affects the behavior of groundwater recession flow, enabling it to be non-exponential in its decay with time.
    'AGWRC': 0.4, # basic groundwater recession rate if KVARY is zero and there is no inflow to groundwater; AGWRC is defined as the rate of flow today divided by the rate of flow yesterday.

    'PETMAX': 4.0,
    'PETMIN': 1.7,
    'INFEXP': 2.0,
    'INFILD': 2.0,
    'DEEPFR': 0.1,  # fraction of groundwater inflow which will enter deep (inactive) groundwater, and, thus, be lost from the system as it is defined in HSPF
    'BASETP': 0.7,  # fraction of remaining potential E-T which can be satisfied from baseflow (groundwater outflow), if enough is available.
    'AGWETP': 0.5,

    'FZG': 0.0394,
    'FZGL': 0.1,
    'LTINFw': 1.0,
    'LTINFb': 0.0,

    'CEPS': 0.0,
    'CEPSC': 0.0,  # interception storage capacity for this canopy layer. This parameter is only used if VCSFG = 0
    'SURS': 0.0,  #  initial interception storage for this canopy layer.
    'UZS': 0.0,
    'UZSN': 0.025,
    'IFWS': 0.0,  # initial interflow storage.
    'LZS': 0.025,  # initial lower zone storage
    'AGWS': 0.0,
    'GWVS': 0.0,

    'NSUR': 0.1,
    'INTFW': 0.0,
    'IRC': 0.25,
    'LZETP': 0.0,

    'VIFWFG': 0,
    'VIRCFG': 0.1,
    'VNNFG': 0,
    'VUZFG': 0,

    'HWTFG': False

}

prec = np.random.random(len(tindex))
ts = {
    'PREC': prec,
    'PETINP': np.add(prec, 0.5)   # Input potential E-T
}

errorV, errorM = pwater(general, ui, ts)

ifwi = np.sum(ts['IFWI'])
deepfr = np.sum(ts['DEEPFR'])
pet = np.sum(ts['PET'])

# Input
_in = np.sum(ts['SUPY'])

# evapotranspiration
cepe = np.sum(ts['CEPE'])
uzet = np.sum(ts['UZET'])
lzet = np.sum(ts['LZET'])
agwet = np.sum(ts['AGWET'])
baset = np.sum(ts['BASET'])
_taet = cepe +  uzet + lzet + agwet + baset
taet = np.sum(ts['TAET'])
d_et = taet-_taet
if -1e-4 > d_et>1e-4:
    print('Problem in ET balance')
print('{:<10} {:<10} {:<10} {:<10} {:<10} {:<10}'.format('cepe', 'uzet', 'lzet', 'agwet', 'baset', 'taet'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f}'.format(cepe, uzet, lzet, agwet, baset, taet))

# Outflow
suro = np.sum(ts['SURO'])
ifwo = np.sum(ts['IFWO'])
agwo = np.sum(ts['AGWO'])
igwi = np.sum(ts['IGWI'])
pero = np.sum(ts['PERO'])
_pero = ifwo + agwo + igwi + suro
d_pero = pero-_pero
print('')
if -1e-4 > d_pero or d_pero>1e-4:
    print('Problem in balance of outflow from PERLND')
print('{:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6}'.format('ifwo', 'agwo', 'igwi', 'suro', 'pero'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f}'.format(ifwo, agwo, igwi, suro, pero))

# Total storage
ceps = np.sum(ts['CEPS'])   # Interception storage (for each
surs = np.sum(ts['SURS'])   # Surface (overland flow) storage
uzs = np.sum(ts['UZS'])
ifws = np.sum(ts['IFWS'])
lzs = np.sum(ts['LZS'])
agws = np.sum(ts['AGWS'])   # Active groundwater storage
tgws = np.sum(ts['TGWS'])   # Total groundwater storage
pers = np.sum(ts['PERS'])   # Total water stored in the PLS
_pers = ceps + surs + uzs + ifws + lzs + agws + tgws
d_pers = pers-_pers
print('')
if -5 > d_pers > 5:
    print('Probelem in Storage balance')
print('{:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6}'
      .format('ceps', 'surs', 'uzs', 'ifws', 'lzs', 'agws', 'tgws', 'pers'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} '
      .format(ceps, surs, uzs, ifws, lzs, agws, tgws, pers))

print('\nTOTAL WATER BALANCE')
print('{:<18} {:<18} {:<18}'.format('Preciptation', 'Evapotranspiration', 'Total Outflow'))
print('{:<18.3f} {:<18.3f} {:<18.3f} '.format(_in, taet, pero))

I am getting following outputs from this code

cepe       uzet       lzet       agwet      baset      taet      
0.000      465.502    0.000      1018.299   0.000      1483.801  
Problem in balance of outflow from PERLND
ifwo       agwo       igwi       suro       pero      
0.000      0.000      53.595     3073.925   3073.925  
ceps       surs       uzs        ifws       lzs        agws       tgws       pers      
0.000      0.000      0.044      0.000      7885.728   0.000      0.000      7885.772   

TOTAL WATER BALANCE
Preciptation       Evapotranspiration Total Outflow     
4076.324           1483.801           3073.925  

which is not correct. Can you tell me where I am making mistake or is this the error in code?

[PaulDudaRESPEC]

I’ve been considering your code above. I think HSPF considers PERO to be the sum of AGWO, IFWO and SURO only – but you also add IGWI, the deep groundwater losses. So I don't think this is an error in the HSPSquared hperwat; it's just an issue with how you are building your water balance summary.

[AtrCheema]

@PaulDudaRESPEC Thank you very much for response. 2 questions. Ignoring IGWI, sum of PERO and TAET is still more than precipitation. This must not come from storage as I am starting all storages with 0.0 or almost 0.0.

If deep percolation is ignored from water balance, how does model maintains water balance, I mean from where this deep percolation is considered to be coming from. I believe deep percolation comes from groundwater storage.

I will really appreciate your response.

[PaulDudaRESPEC]

I suggest you try running with DEEPFR = 0.0, and then comparing the results with these with DEEPFR = 0.1. I think that will help clarify what is happening.

[AtrCheema]

@PaulDudaRESPEC Keeping DEEPFR=0.0 makes IGWI=0.0 while making it 0.1 increases IGWI which make quite sense but total water balance still remains same. Following are the results using same code as above except changing the value of DEEPFR

DEEPFR=0.0

cepe       uzet       lzet       agwet      baset      taet      
0.000      465.502    0.000      1071.894   0.000      1537.396  
ifwo       agwo       igwi       suro       pero      
0.000      0.000      0.000      3073.925   3073.925  
ceps       surs       uzs        ifws       lzs        agws       tgws       pers      
0.000      0.000      0.044      0.000      7885.728   0.000      0.000      7885.772   
TOTAL WATER BALANCE
Preciptation       Evapotranspiration Total Outflow     
4076.324           1537.396           3073.925 

DEEPFR =0.1

cepe       uzet       lzet       agwet      baset      taet      
0.000      465.502    0.000      1018.299   0.000      1483.801  
Problem in balance of outflow from PERLND
ifwo       agwo       igwi       suro       pero      
0.000      0.000      53.595     3073.925   3073.925  
ceps       surs       uzs        ifws       lzs        agws       tgws       pers      
0.000      0.000      0.044      0.000      7885.728   0.000      0.000      7885.772   
TOTAL WATER BALANCE
Preciptation       Evapotranspiration Total Outflow     
4076.324           1483.801           3073.925  

[PaulDudaRESPEC]

Since the storage terms are state variables, I suggest writing out the values for surs, uzs, lzs, ... pers after the last time step of the simulation (instead of summing them). Then comparing to the initial values you can see how much each storage term has changed over the course of the simulation.

[AtrCheema]

@PaulDudaRESPEC I really appreciate your time. I calculated storages and changes in storages but since I am using zero initial storages, there is increase in storages after simulation which means some part of SUPY must also account for this increase in storages. The code to calculate is following

import pandas as pd
import numpy as np

np.random.seed(seed=77)

from perwat import pwater  # import pwater function here

tindex = pd.date_range('20110101', '20111210', freq='H')

general = {
    'sim_len': len(tindex),
    'tindex': tindex,
    'sim_delt':60  # number of minutes in time step
}

ui = {
    'CSNOFG': 0,
    'RTOPFG': 1,  #
    'UZFG': 0,   # algorithm used with `0` causes a lot of surface runoff (thus higher pero) at the expense of uzet
    'VCSFG': 0,
    'VLEFG': 0,
    'IFFCFG': 0,
    'IFRDFG': 0,

    'FOREST': 0.0,
    'LZSN': 0.25,
    'INFILT': 1.0,  #  index to the infiltration capacity of the soil.
    'LSUR': 1.0,
    'SLSUR': 1.5,
    'KVARY': 1.0, # affects the behavior of groundwater recession flow, enabling it to be non-exponential in its decay with time.
    'AGWRC': 0.4, # basic groundwater recession rate if KVARY is zero and there is no inflow to groundwater; AGWRC is defined as the rate of flow today divided by the rate of flow yesterday.

    'PETMAX': 4.0,
    'PETMIN': 1.7,
    'INFEXP': 2.0,
    'INFILD': 2.0,
    'DEEPFR': 0.0,  # fraction of groundwater inflow which will enter deep (inactive) groundwater, and, thus, be lost from the system as it is defined in HSPF
    'BASETP': 0.7,  # fraction of remaining potential E-T which can be satisfied from baseflow (groundwater outflow), if enough is available.
    'AGWETP': 0.5,

    'FZG': 0.0394,
    'FZGL': 0.1,
    'LTINFw': 1.0,
    'LTINFb': 0.0,

    'CEPS': 0.0,
    'CEPSC': 0.0,  # interception storage capacity for this canopy layer. This parameter is only used if VCSFG = 0
    'SURS': 0.0,  #  initial interception storage for this canopy layer.
    'UZS': 0.0,
    'UZSN': 0.025,
    'IFWS': 0.0,  # initial interflow storage.
    'LZS': 0.025,  # initial lower zone storage
    'AGWS': 0.0,
    'GWVS': 0.0,

    'NSUR': 0.1,
    'INTFW': 0.0,
    'IRC': 0.25,
    'LZETP': 0.0,

    'VIFWFG': 0,
    'VIRCFG': 0.1,
    'VNNFG': 0,
    'VUZFG': 0,

    'HWTFG': False

}

prec = np.random.random(len(tindex))
ts = {
    'PREC': prec,
    'PETINP': np.add(prec, 0.5)   # Input potential E-T
}

errorV, errorM = pwater(general, ui, ts)

ifwi = np.sum(ts['IFWI'])
deepfr = np.sum(ts['DEEPFR'])
pet = np.sum(ts['PET'])
infil = np.sum(ts['INFIL'])   # mm/ivld  infilteration to the soil
# cepo= np.sum(ts['CEPO'])

# Input
_in = np.sum(ts['SUPY'])

# evapotranspiration
cepe = np.sum(ts['CEPE'])
uzet = np.sum(ts['UZET'])
lzet = np.sum(ts['LZET'])
agwet = np.sum(ts['AGWET'])
baset = np.sum(ts['BASET'])
_taet = cepe +  uzet + lzet + agwet + baset
taet = np.sum(ts['TAET'])
d_et = taet-_taet
print('\nEvapotranspiration')
if -1e-4 > d_et>1e-4:
    print('Problem in ET balance')
print('{:<10} {:<10} {:<10} {:<10} {:<10} {:<10}'.format('cepe', 'uzet', 'lzet', 'agwet', 'baset', 'taet'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f}'.format(cepe, uzet, lzet, agwet, baset, taet))

# Outflow
suro = np.sum(ts['SURO'])
ifwo = np.sum(ts['IFWO'])
agwo = np.sum(ts['AGWO'])
igwi = np.sum(ts['IGWI'])
pero = np.sum(ts['PERO'])
_pero = ifwo + agwo + suro
d_pero = pero-_pero
print('\n Outflows')
if -1e-4 > d_pero or d_pero>1e-4:
    print('Problem in balance of outflow from PERLND')
print('{:<10.6} {:<10.6} {:<10.6} {:<10.6} '.format('ifwo', 'agwo',  'suro', 'pero'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f}'.format(ifwo, agwo, suro, pero))

# Total storage
ceps = np.sum(ts['CEPS'])   # Interception storage (for each
surs = np.sum(ts['SURS'])   # Surface (overland flow) storage
uzs = np.sum(ts['UZS'])
ifws = np.sum(ts['IFWS'])
lzs = np.sum(ts['LZS'])
agws = np.sum(ts['AGWS'])   # Active groundwater storage
tgws = np.sum(ts['TGWS'])   # Total groundwater storage
pers = np.sum(ts['PERS'])   # Total water stored in the PLS
_pers = ceps + surs + uzs + ifws + lzs + agws + tgws
d_pers = pers-_pers
print('\nSTORAGES')
if -5 > d_pers > 5:
    print('Probelem in Storage balance')
print('{:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6}'
      .format('ceps', 'surs', 'uzs', 'ifws', 'lzs', 'agws', 'tgws', 'pers'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} '
      .format(ceps, surs, uzs, ifws, lzs, agws, tgws, pers))

# change in storage
d_ceps = ceps - ui['CEPS']
d_surs = surs - ui['SURS']
d_uzs = uzs - ui['UZS']
d_ifws = ifws - ui['IFWS']
d_lzs = lzs - ui['LZS']
d_agws = agws - ui['AGWS']
print('\nChanges in Storages \n{:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6} {:<10.6}'
      .format('d_ceps', 'd_surs', 'd_uzs', 'd_ifws', 'd_lzs', 'd_agws'))
print('{:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f} {:<10.3f}'
      .format(d_ceps, d_surs, d_uzs, d_ifws, d_lzs, d_agws))

print('\nTOTAL WATER BALANCE')
print('{:<18} {:<18} {:<18}'.format('Preciptation', 'Evapotranspiration', 'Total Outflow'))
print('{:<18.3f} {:<18.3f} {:<18.3f} '.format(_in, taet, pero))

The output from this code is

Evapotranspiration

cepe       uzet       lzet       agwet      baset      taet      
0.000      465.502    0.000      1071.894   0.000      1537.396  

Outflows

ifwo       agwo       suro       pero       
0.000      0.000      3073.925   3073.925  

STORAGES

ceps       surs       uzs        ifws       lzs        agws       tgws       pers      
0.000      0.000      0.044      0.000      7885.728   0.000      0.000      7885.772   

Changes in Storages

d_ceps     d_surs     d_uzs      d_ifws     d_lzs      d_agws    
0.000      0.000      0.044      0.000      7885.703   0.000     

TOTAL WATER BALANCE

Preciptation       Evapotranspiration Total Outflow     
4076.324           1537.396           3073.925 

[PaulDudaRESPEC]

Are you running the official pwater code, or have you modified that code somehow? I'm suspicious that you are getting more agwet than I'd expect.

[AtrCheema]

@PaulDudaRESPEC Yes that was the case. I 'somehow' made some change in code which led to oeverestimation of agwet.

Thank you very much for your time.