Estimating water balance from .SAD file

[Carlos Alberto Arnillas]

Hello

I am trying to estimate the water balance (as printed at the end of the site***.OUT file) using the .SAD files. My codes (using R) are below. As can be noticed, no value matches with the ones reported in the .OUT file, not even the before/after values. What am I doing wrong? (In yellow the ones that are a close match).

Thanks in advance

Carlos Alberto

sa <- sao[,by=SAID,.(WSA=unique(WSAha))] # Table with the area per subarea

sad[,date := as.IDate(sprintf("%4d-%2d-%2d", Y, M, D))] # .SAD table

# focusing on the water components

w_sa <- unique(sad[,.(SAID, date, Y, M, D, PRCP, ET, Q, QIN, PRK,

                   DPRK, WYLD, SSF, RSSF, SSFI, QRF, SW, RZSW, GWST)])

# Estimating the total water flows

w_sa[sa, on="SAID", .(PRCP=sum(PRCP*WSA*10),  # Precipitation

                      ET=sum(ET*WSA*10),      # Evapotranspiration

                      Q=sum(Q*WSA*10),        # Flow

                      Q13=sum(Q*WSA*(SAID==13)*10), # Flow from the outflow only

                      PRK=sum(PRK*WSA*10),    # Percolation

                      QRF=sum(QRF*WSA*10),    # Quick return flow

                      WYLD=sum(WYLD*WSA*10),  # Water yield

                      DPRK=sum(DPRK*WSA*10))] # Deep percolation

# Results:

#        PRCP       ET       Q     Q13     PRK     QRF    WYLD    DPRK

# 1: 35377089 18319300 8526136 55004.1 8363367 609.715 8739426 8234800

# Values before and after

w_sa[date %in% range(w_sa$date)][sa, on="SAID"][, by=date, 

                                                .(n=.N, SW=sum(SW*WSA*10), RZSW=sum(RZSW*WSA*10), GWST=sum(GWST*WSA*10))]

         date  n         SW       RZSW     GWST

1: 2002-01-01 28   1135.431   1135.431 28791490

2: 2015-12-31 28 339030.372 339030.372 10757589

# From SITE13.OUT

# TOTAL WATER BALANCE (m3)

# 

# PER =-0.740862D+01  DF  =-0.653582D+05  BSW = 0.642210D+06  

# BGWS= 0.125815D+06  BRSV= 0.000000D+00  BSNO= 0.000000D+00

# BSWL= 0.000000D+00  PCP = 0.353898D+08  WYLD= 0.867023D+07 

# DPRK= 0.823769D+07  ET  = 0.183257D+08  IRG = 0.000000D+00

# QIN = 0.000000D+00  PSOQ= 0.000000D+00  EVRT= 0.860678D+04 

# RSIR= 0.000000E+00  WLIR= 0.000000E+00  IRDL= 0.000000D+00

# FSW = 0.882192D+06  FGWS= 0.427163D+05  FRSV= 0.000000D+00 

# FSNO= 0.560344D+05  FSWL= 0.000000D+00

[Luca]

Carlos,

Which version of APEX are you using and how many subareas are included in the watershed you are simulating?

I did a quick test simulating a single subarea and using Excel to aggregate daily values of precipitations and evapotranspiration and I got a perfect match for the precipitation while the evapotranspiration was off by 4.4 m3.

Best,

Luca

[Carlos Alberto Arnillas]

Hello Luca

I am using the version that comes with ArcAPEX, and working with 28 subareas.

What do you think causes the mismatch in evapotranspiration? I know it is not large, but wondering if the same problem could became larger when weights have to be included.

(please, let me know if you want me to post my input/output text files)

Carlos Alberto

[Luca]

Carlos,

I spent some time with your files and did some tests. I calculated PRCP and ET with Excel and then with R. I got the same results with both methods and they were different from the values you have calculated.

I have calculated a PRCP of  35,389,746 m3, and an ET of 18,325,660 m3wich are quite close to the values reported in the OUT file. I think the difference is due to rounding when the daily values are printed in the SAD file.

Here is the code I used in R

library(readr)
library(dplyr)

# Reading a csv file with subarea ID and size

wsa <- read.csv(file = "D:/Debug run APEX/APEX forum Carlos Alberto water balance TxtInOut/SA_WSA.csv",
                header = T, sep = ",", dec = ".")

# Reading the SAD file
sadFile <- read_delim(file = "SITE13.SAD", delim = " ", skip = 9, col_names = T, trim_ws = T)

# Working with the SAD file to get the cumulative values for precipitation and evapotranspiration
wkSAD <- sadFile %>%
    select(., SA, ID, Y, M, D, PRCP, ET) %>%
    group_by(., SA) %>%
    summarise(., PRCP = sum(PRCP), ET = sum(ET)) %>%
    left_join(., wsa, by = 'SA') %>%
    mutate(., PRCPm3 = PRCP * WSA * 10,
           ETm3 = ET * WSA * 10)

Try to revise your code to see if you can get better results and let me know if you want me to check any other variables.

Luca

[Carlos Alberto Arnillas]

Hello Luca

Thanks for your help. I reviewed the code, and found that the mismatch was caused because I was using the area reported in the SAO file, which is a rounded version of the area reported in the original SUB file.

Now I am trying to follow Rachel Mason et al.'s paper to check the overall water balance. Is it possible to estimate all the watershed summary stats from the SAD file or other daily files? If so, what would be the formulas to do it? (or where can I find them). So far, I only managed to find how to estimate PCP, ET and DPRK. 
Thanks again!

Carlos Alberto

[Luca]

Carlos,

It should be possible to estimate the balance using the daily data provided in the SAD file. I have never done it but, all the aggregated data (monthly, annual, and total) values are calculated starting from the daily values used during the simulation. As I said, I have never attempted it before but every calculation should be straightforward.

I don't have any documentation on how to aggregate the daily value of each variable but, try a simple aggregation first and let me know if you see any unusual results and I will be happy to check how that aggregated value is calculated in the model.

Let me know if you need more information.

Luca

[Carlos Alberto Arnillas]

Thanks Luca

Here I'm trying to link Rachel et al.'s definitions of the water balance with the equations that I am using, the results I got, and the abbreviations that I think correspond to the .OUT file. I tried to add some colours for the values that I was able to match.

These three ones are straightforward because in all them the balance does not have to account for water movement between subareas, and match the values reported.

PCP:    Precipitation            sum(PRCP*abs(WSA)*10)    35389746   PCP

ET:     Evapotranspiration       sum(ET*abs(WSA)*10)      18325660   ET
DPRK:   Deep percolation         sum(DPRK*abs(WSA)*10)     8237695   DPRK

The before/after numbers should be straightforward,  but I can't match my results with any of the ones provided by APEX (see below), except for snow

Balance (before/after)

BSW:    Shallow water (before)   sum(SW*abs(WSA)*10) on the first date      1135.737    BSW, BGWS?

FSW:    Shallow water (final)    sum(SW*abs(WSA)*10) on the last date     339153.159    FSW, FGWS?

BSNO:   Snow (final)             sum(SNO*abs(WSA)*10) on the first date        0

FSNO:   Snow (final)             sum(SNO*abs(WSA)*10) on the last date     56034.35

QIN = IRG = 0 is fine, as there is no input of water into the watershed besides precipitation.  And with BRSV = FRSV = 0 because there are no reservoirs in the area.

No clue how to estimate the other parameters, like surface runoff, (SSF) percolation (PRK), subsurface outflow (SSO), or quick return flow (QRF). All exist —except SSO— in the .SAD file, but as they represent flow between subareas, they are probably meaningless. 

When I focus on the water yield (WYLD) of the subarea in the outlet, I got 

WYLD13:  sum(WYLD*abs(WSA)*(SA==13)*10)     33421.93

And adding the water yield from each subarea I got:

WYLD13:  sum(WYLD*abs(WSA)*(SA==13)*10)   8742525 

The former is very far from the PRCP - ET - DPRK = 8826391. The second suppose that the water is flowing from a field to a channel and that the channel water is not considered as an input for a subarea (only to its channel).

In bold, the variables that I think I figure out. The others, are unclear.

# From SITE13.OUT
# TOTAL WATER BALANCE (m3)
#
# PER =-0.740862D+01 DF  =-0.653582D+05 BSW = 0.642210D+06
# BGWS= 0.125815D+06 BRSV= 0.000000D+00 BSNO= 0.000000D+00
# BSWL= 0.000000D+00 PCP = 0.353898D+08 WYLD= 0.867023D+07
# DPRK= 0.823769D+07 ET  = 0.183257D+08 IRG = 0.000000D+00
# QIN = 0.000000D+00 PSOQ= 0.000000D+00 EVRT= 0.860678D+04
# RSIR= 0.000000E+00 WLIR= 0.000000E+00 IRDL= 0.000000D+00
# FSW = 0.882192D+06 FGWS= 0.427163D+05 FRSV= 0.000000D+00
# FSNO= 0.560344D+05 FSWL= 0.000000D+00

Thanks again for all your help.

Carlos Alberto

[Luca]

Carlos,

This time I didn't have time to try to calculate the values from the daily data like I did before but, here is how the model calculates some of the variables included in the final water balance.

BSW is the soil water content at the beginning of the simulation reported in m3. It is calculated as the sum of the soil water content of each subarea at the beginning of the simulation (not even at the end of day one). I think you can get this information from the OUT file (at the very beginning of it) and from the DCN file where the first set of data are referred to the initial conditions.

FSW is the soil water content at the end of the simulation reported in m3. It is calculated as the sum of the soil water content of each subarea on the last day of the simulation. You will find the daily data on the DCN fil.

BGWS is the groundwater storage at the beginning of the simulation in m3. It is calculated as the groundwater stored in each subarea at the beginning of the simulation. Probably, the value on the first day of the SAD file will give you a good approximation of the initial value.

FGWS is the groundwater storage at the end of the simulation in m3. It is calculated as the groundwater stored in each subarea on the last day of the simulation. You can use the last value reported in the SAD file.

EVRT is the evaporation from the flow from all the subareas included in the watershed. You should be able to print it in the SAD file adding the variable 147 in the PRINT file (let me know if you need more information on this).

WYLD is the sum of the water yield from the subarea where the outlet of the watershed is located (or this is what I can understand looking at the code of the model). You should be able to get it with the daily values. The model calculates it accumulating the daily values in several steps so, I am wondering if the rather large difference can be caused by rounding...

Let me know what you find and if you need more information.

Luca