CWatM Malawi Calibration summary

Malawi CWatM

Key considerations for Malawian context

Domestic water usage: Groundwater is the major source of domestic water in Malawi, around 75% of the population of Malawi are reliant on groundwater for their everyday needs. This is higher within rural areas and lower in urban areas, specific estimates vary but the BGS estimates ~82% of the population relying on groundwater for domestic purposes in rural areas and 20% in urban areas where gravity fed piped water supply systems are more common (https://earthwise.bgs.ac.uk/index.php/Hydrogeology_of_Malawi#:~:text=Malawi%20has%20abundant%20surface%20water,for%20irrigation%20by%20smallholder%20farmers.)

Groundwater: We have data of over 100,000 boreholes across Malawi which give an indication of where seasonal water shortages are experienced.

The main cities of Blantyre and Lilongwe have piped water supply systems.

Blantyre: The Blantrye Water Board extracts water from the Shire River around 40km from Blantyre (Walker’s Ferry), it is also 800m below Blantyre so requires pumping to the city. Blantyre Water Board currently processes 96 million litres of water/day. There is also 5 million litres conveyed to the city from Mudi Dam and an additional planned 20 million litres from the construction of a new source station at Likhubula in Mulanje. The water board currently supplies around 1.4 million people. (Source Blantrye Water Board).

Lilongwe: Water comes from the Lilongwe River, however, current output is insufficient for demand and there can be irregular water supply. There are two main treatment works supplying Lilongwe, able to supply 35,000 m3/day and 60,000 m3/day respectively.

Source:

https://www.water-alternatives.org/index.php/alldoc/articles/vol11/v11issue2/439-a11-2-6/file

https://www.diva-portal.org/smash/get/diva2:1264613/FULLTEXT01.pdf

Key considerations for CWatM: Groundwater still plays a significant role in domestic water usage within Malawi, however, piped surface water is a significant consideration within urban areas. Piped supply systems are struggling to keep up with demand.

Irrigation: Whilst a few irrigation schemes were developed in Malawi prior to 1960, these were limited. Between 1967 and 1982 16 irrigation schemes were built alongside the shores of Lake Malawi, the Lake Chilwa basin and the Lower Shire River Basin [@gwiyani-nkhoma2011].

The 2015 Government of Malawi Irrigation master plan stated that there is an estimated potential for 407,862 hectares of irrigated land in Malawi with only 104,634 having been developed having been developed by 2015. The area under irrigation is estimated to be growing since 2006 at approximately 6% per year (https://kulimamalawi.org/wp-content/uploads/2019/12/NIP_Final_Signed.pdf).

A 2019 report by the government of Malawi estimated that 116 249 hectares of land was under irrigation, with roughly half reported by commercial farmers and half by small-holder farmers (59, 655 hectares smallholder farmers and 56,594 hectares for the private/commercial farmers. The most common irrigation system for small holder farmers is gravity fed systems (47%), pumps (both treadle and motarized) are utilized by 43% of small holder farmers and watering cans are used by the remaining 10%. Source: https://agriculture.gov.mw/docs/Water%20Development.pdf

Around 80% of the population of Malawi are engaged in small-holder farming. It is quite common to use waste domestic water for small-scale watering/ informal irrigation. Capturing the level of small holder irrigation may be an important consideration.

Key takeaway for CWatM: Irrigation is increasing within Malawi and having a scenario of increased irrigation (to the maximum calculated capacity) would be helpful. About of half irrigation projects are small holder farmers and quite informal irrigation structures, this may mean that irrigation is currently underestimated in CWatM (if it accounts for formalized, commercial irrigation systems only).

Industrial water withdrawal: Industrial water usage currently represents a small fraction of total water withdrawal. Estimates are around 4% of total water consumption. This is largely groundwater with some estimates around 95% groundwater dependence for industry.

Sources:

https://knoema.com/atlas/Malawi/topics/Water/Water-Withdrawal/Industrial-water-withdrawal#:~:text=In%202020%2C%20industrial%20water%20withdrawal,average%20annual%20rate%20of%200.75%25.

https://centaur.reading.ac.uk/96540/17/Dossier_Country_Malawi.pdf

Lake Malawi considerations:

Main inputs: Rain, inflow from 4 main rivers

Outputs: Evaporation over lake, outflow into Shire River.

The outflow into the shire river is important for hydropower generation in Malawi (about 90% of Malawi’s hydropower generation is here).

https://link.springer.com/article/10.1007/s10661-020-08519-4

https://afjrd.org/index.php/jos/article/view/201

Calibration

See CWatM settings

Change path to location of ‘bigrun1’ file download

library(devtools)
#devtools::install_github("dof1985/cwatmRutils")
library(cwatmRutils)
library(raster)
library(ggplot2)
library(ncdf4)
#Define path root for run looking to analyse
pth <- "C:/Users/hinton/Documents/GitHub/MalawiLake/outputs/bigrun1"

Visualise discharge raster of CwatM output as a plot

fn <- "discharge_daily.nc"
monthly_discharge <- ncdf2raster(pth = sprintf("%s/%s", pth, fn),
                                 transpose = TRUE, temporal_fun = "mean", na.rm = TRUE)
plot(monthly_discharge)

Comparison to GRDC discharge data

#Load monthly discharge data for run from 1970 to 2010

fn <- "discharge_monthavg.nc"

monthly_discharge <- ncdf2raster(pth = sprintf("%s/%s", pth, fn),
                 transpose = TRUE, time = c(1,12*40), fun = mean, na.rm = TRUE)
monthly_discharge$date <- as.Date(monthly_discharge$time, origin = "1901-01-01")

# 
library(mondate)
start<- mondate("1971-01-01")
end<- mondate("2010-12-31")
as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))

  [1] "1971-01-01" "1971-02-01" "1971-03-01" "1971-04-01" "1971-05-01" "1971-06-01"
  [7] "1971-07-01" "1971-08-01" "1971-09-01" "1971-10-01" "1971-11-01" "1971-12-01"
 [13] "1972-01-01" "1972-02-01" "1972-03-01" "1972-04-01" "1972-05-01" "1972-06-01"
 [19] "1972-07-01" "1972-08-01" "1972-09-01" "1972-10-01" "1972-11-01" "1972-12-01"
 [25] "1973-01-01" "1973-02-01" "1973-03-01" "1973-04-01" "1973-05-01" "1973-06-01"
 [31] "1973-07-01" "1973-08-01" "1973-09-01" "1973-10-01" "1973-11-01" "1973-12-01"
 [37] "1974-01-01" "1974-02-01" "1974-03-01" "1974-04-01" "1974-05-01" "1974-06-01"
 [43] "1974-07-01" "1974-08-01" "1974-09-01" "1974-10-01" "1974-11-01" "1974-12-01"
 [49] "1975-01-01" "1975-02-01" "1975-03-01" "1975-04-01" "1975-05-01" "1975-06-01"
 [55] "1975-07-01" "1975-08-01" "1975-09-01" "1975-10-01" "1975-11-01" "1975-12-01"
 [61] "1976-01-01" "1976-02-01" "1976-03-01" "1976-04-01" "1976-05-01" "1976-06-01"
 [67] "1976-07-01" "1976-08-01" "1976-09-01" "1976-10-01" "1976-11-01" "1976-12-01"
 [73] "1977-01-01" "1977-02-01" "1977-03-01" "1977-04-01" "1977-05-01" "1977-06-01"
 [79] "1977-07-01" "1977-08-01" "1977-09-01" "1977-10-01" "1977-11-01" "1977-12-01"
 [85] "1978-01-01" "1978-02-01" "1978-03-01" "1978-04-01" "1978-05-01" "1978-06-01"
 [91] "1978-07-01" "1978-08-01" "1978-09-01" "1978-10-01" "1978-11-01" "1978-12-01"
 [97] "1979-01-01" "1979-02-01" "1979-03-01" "1979-04-01" "1979-05-01" "1979-06-01"
[103] "1979-07-01" "1979-08-01" "1979-09-01" "1979-10-01" "1979-11-01" "1979-12-01"
[109] "1980-01-01" "1980-02-01" "1980-03-01" "1980-04-01" "1980-05-01" "1980-06-01"
[115] "1980-07-01" "1980-08-01" "1980-09-01" "1980-10-01" "1980-11-01" "1980-12-01"
[121] "1981-01-01" "1981-02-01" "1981-03-01" "1981-04-01" "1981-05-01" "1981-06-01"
[127] "1981-07-01" "1981-08-01" "1981-09-01" "1981-10-01" "1981-11-01" "1981-12-01"
[133] "1982-01-01" "1982-02-01" "1982-03-01" "1982-04-01" "1982-05-01" "1982-06-01"
[139] "1982-07-01" "1982-08-01" "1982-09-01" "1982-10-01" "1982-11-01" "1982-12-01"
[145] "1983-01-01" "1983-02-01" "1983-03-01" "1983-04-01" "1983-05-01" "1983-06-01"
[151] "1983-07-01" "1983-08-01" "1983-09-01" "1983-10-01" "1983-11-01" "1983-12-01"
[157] "1984-01-01" "1984-02-01" "1984-03-01" "1984-04-01" "1984-05-01" "1984-06-01"
[163] "1984-07-01" "1984-08-01" "1984-09-01" "1984-10-01" "1984-11-01" "1984-12-01"
[169] "1985-01-01" "1985-02-01" "1985-03-01" "1985-04-01" "1985-05-01" "1985-06-01"
[175] "1985-07-01" "1985-08-01" "1985-09-01" "1985-10-01" "1985-11-01" "1985-12-01"
[181] "1986-01-01" "1986-02-01" "1986-03-01" "1986-04-01" "1986-05-01" "1986-06-01"
[187] "1986-07-01" "1986-08-01" "1986-09-01" "1986-10-01" "1986-11-01" "1986-12-01"
[193] "1987-01-01" "1987-02-01" "1987-03-01" "1987-04-01" "1987-05-01" "1987-06-01"
[199] "1987-07-01" "1987-08-01" "1987-09-01" "1987-10-01" "1987-11-01" "1987-12-01"
[205] "1988-01-01" "1988-02-01" "1988-03-01" "1988-04-01" "1988-05-01" "1988-06-01"
[211] "1988-07-01" "1988-08-01" "1988-09-01" "1988-10-01" "1988-11-01" "1988-12-01"
[217] "1989-01-01" "1989-02-01" "1989-03-01" "1989-04-01" "1989-05-01" "1989-06-01"
[223] "1989-07-01" "1989-08-01" "1989-09-01" "1989-10-01" "1989-11-01" "1989-12-01"
[229] "1990-01-01" "1990-02-01" "1990-03-01" "1990-04-01" "1990-05-01" "1990-06-01"
[235] "1990-07-01" "1990-08-01" "1990-09-01" "1990-10-01" "1990-11-01" "1990-12-01"
[241] "1991-01-01" "1991-02-01" "1991-03-01" "1991-04-01" "1991-05-01" "1991-06-01"
[247] "1991-07-01" "1991-08-01" "1991-09-01" "1991-10-01" "1991-11-01" "1991-12-01"
[253] "1992-01-01" "1992-02-01" "1992-03-01" "1992-04-01" "1992-05-01" "1992-06-01"
[259] "1992-07-01" "1992-08-01" "1992-09-01" "1992-10-01" "1992-11-01" "1992-12-01"
[265] "1993-01-01" "1993-02-01" "1993-03-01" "1993-04-01" "1993-05-01" "1993-06-01"
[271] "1993-07-01" "1993-08-01" "1993-09-01" "1993-10-01" "1993-11-01" "1993-12-01"
[277] "1994-01-01" "1994-02-01" "1994-03-01" "1994-04-01" "1994-05-01" "1994-06-01"
[283] "1994-07-01" "1994-08-01" "1994-09-01" "1994-10-01" "1994-11-01" "1994-12-01"
[289] "1995-01-01" "1995-02-01" "1995-03-01" "1995-04-01" "1995-05-01" "1995-06-01"
[295] "1995-07-01" "1995-08-01" "1995-09-01" "1995-10-01" "1995-11-01" "1995-12-01"
[301] "1996-01-01" "1996-02-01" "1996-03-01" "1996-04-01" "1996-05-01" "1996-06-01"
[307] "1996-07-01" "1996-08-01" "1996-09-01" "1996-10-01" "1996-11-01" "1996-12-01"
[313] "1997-01-01" "1997-02-01" "1997-03-01" "1997-04-01" "1997-05-01" "1997-06-01"
[319] "1997-07-01" "1997-08-01" "1997-09-01" "1997-10-01" "1997-11-01" "1997-12-01"
[325] "1998-01-01" "1998-02-01" "1998-03-01" "1998-04-01" "1998-05-01" "1998-06-01"
[331] "1998-07-01" "1998-08-01" "1998-09-01" "1998-10-01" "1998-11-01" "1998-12-01"
[337] "1999-01-01" "1999-02-01" "1999-03-01" "1999-04-01" "1999-05-01" "1999-06-01"
[343] "1999-07-01" "1999-08-01" "1999-09-01" "1999-10-01" "1999-11-01" "1999-12-01"
[349] "2000-01-01" "2000-02-01" "2000-03-01" "2000-04-01" "2000-05-01" "2000-06-01"
[355] "2000-07-01" "2000-08-01" "2000-09-01" "2000-10-01" "2000-11-01" "2000-12-01"
[361] "2001-01-01" "2001-02-01" "2001-03-01" "2001-04-01" "2001-05-01" "2001-06-01"
[367] "2001-07-01" "2001-08-01" "2001-09-01" "2001-10-01" "2001-11-01" "2001-12-01"
[373] "2002-01-01" "2002-02-01" "2002-03-01" "2002-04-01" "2002-05-01" "2002-06-01"
[379] "2002-07-01" "2002-08-01" "2002-09-01" "2002-10-01" "2002-11-01" "2002-12-01"
[385] "2003-01-01" "2003-02-01" "2003-03-01" "2003-04-01" "2003-05-01" "2003-06-01"
[391] "2003-07-01" "2003-08-01" "2003-09-01" "2003-10-01" "2003-11-01" "2003-12-01"
[397] "2004-01-01" "2004-02-01" "2004-03-01" "2004-04-01" "2004-05-01" "2004-06-01"
[403] "2004-07-01" "2004-08-01" "2004-09-01" "2004-10-01" "2004-11-01" "2004-12-01"
[409] "2005-01-01" "2005-02-01" "2005-03-01" "2005-04-01" "2005-05-01" "2005-06-01"
[415] "2005-07-01" "2005-08-01" "2005-09-01" "2005-10-01" "2005-11-01" "2005-12-01"
[421] "2006-01-01" "2006-02-01" "2006-03-01" "2006-04-01" "2006-05-01" "2006-06-01"
[427] "2006-07-01" "2006-08-01" "2006-09-01" "2006-10-01" "2006-11-01" "2006-12-01"
[433] "2007-01-01" "2007-02-01" "2007-03-01" "2007-04-01" "2007-05-01" "2007-06-01"
[439] "2007-07-01" "2007-08-01" "2007-09-01" "2007-10-01" "2007-11-01" "2007-12-01"
[445] "2008-01-01" "2008-02-01" "2008-03-01" "2008-04-01" "2008-05-01" "2008-06-01"
[451] "2008-07-01" "2008-08-01" "2008-09-01" "2008-10-01" "2008-11-01" "2008-12-01"
[457] "2009-01-01" "2009-02-01" "2009-03-01" "2009-04-01" "2009-05-01" "2009-06-01"
[463] "2009-07-01" "2009-08-01" "2009-09-01" "2009-10-01" "2009-11-01" "2009-12-01"
[469] "2010-01-01" "2010-02-01" "2010-03-01" "2010-04-01" "2010-05-01" "2010-06-01"
[475] "2010-07-01" "2010-08-01" "2010-09-01" "2010-10-01" "2010-11-01" "2010-12-01"

###
monthly_discharge$date <- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))


ggplot(monthly_discharge, aes(date, value)) + geom_line() 

#+ facet_grid(pts~., scales = "free_y") 

monthly_discharge<- as.data.frame(monthly_discharge)
###
#Identified 2 point locations to compare discharge (upstream and downstream of Blantyre)

#Discharge at chiromo station (Downstream of Blantyre)
## Chiromo data station -16.55, 35.13 (data from GRDC location)

pt<- rbind(as.numeric(c('35.125', '-16.375')))

pt<- as.data.frame(pt)
colnames(pt)<- c('x','y')

# extract data for the point(s)
point_extract <- ncdf2raster(pth = sprintf("%s/%s", pth, fn),
                 transpose = TRUE, time = c(1, 12*40 ), spatial = pt)

point_extract$date <- as.Date(point_extract$time, origin = "1901-01-01")

#x$pts <- rep(c("outlet", "gauge2"), nrow(x) / 2)
# visualize
ggplot(point_extract, aes(date, value)) + geom_line() 

#+ facet_grid(pts~., scales = "free_y")  
#+ 
#+ #### Compare with data from CWatM output (both at point location and also running average CWatM basin wide discharge)
monthly_discharge_output<- as.data.frame(monthly_discharge$value)
monthly_discharge_output<- monthly_discharge_output[1:480,]
monthly_discharge_output<- as.data.frame(monthly_discharge_output[1:480])
monthly_discharge_output$group<- "CwatM average basin"
colnames(monthly_discharge_output)<- c("Value", "group")

measured_discharge_output<- read.csv("H:/MyDocuments/Validation data/GRDC Discharge data/1992900/monthly measured discharge output.csv")
measured_discharge_output<- as.data.frame(measured_discharge_output[1:480,])
measured_discharge_output$group<- "measured chiromo station"
colnames(measured_discharge_output)<- c("Value", "group")

comparisondf<- rbind(monthly_discharge_output,measured_discharge_output)

####
#comparisondf<- comparisondf[1:36,]

comparisondf$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))
comparisondf$Value<- as.numeric(comparisondf$Value)

Warning: NAs introduced by coercion

####

point_data<- as.numeric(point_extract$value)
point_data<- as.data.frame(point_data)
point_data$group<- "CWatM chiromo station"
point_data$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))
colnames(point_data)<- c("Value", "group", "date")

comparisondf<- rbind(point_data, comparisondf)

ggplot(comparisondf, aes(x= date, y= Value, colour=group)) + geom_line() 

######### 
#Compare with discharge station upstream of Blantrye
#Discharge at Liwonde station

#pt_liwonde<- rbind(as.numeric(c('35.2', '-15.07')))
pt_liwonde<- rbind(as.numeric(c('35.29167', '-15.04167')))

pt_liwonde<- as.data.frame(pt_liwonde)
colnames(pt_liwonde)<- c('x','y')

# extract data for the point(s)
point_extract_liwonde <- ncdf2raster(pth = sprintf("%s/%s", pth, fn),
                 transpose = TRUE, time = c(1, 12*40 ), spatial = pt_liwonde)

point_extract_liwonde$date <- as.Date(point_extract_liwonde$time, origin = "1901-01-01")

#x$pts <- rep(c("outlet", "gauge2"), nrow(x) / 2)
# visualize
ggplot(point_extract_liwonde, aes(date, value)) + geom_line() 

#+ facet_grid(pts~., scales = "free_y")  
#+ #############################
measured_discharge_output_liwonde<- read.csv("H:/MyDocuments/Validation data/GRDC Discharge data/1992700/Monthly discharge summary.csv")
measured_discharge_output_liwonde<- as.data.frame(measured_discharge_output_liwonde[1:480,2])
measured_discharge_output_liwonde$group<- "measured_liwonde"
colnames(measured_discharge_output_liwonde)<- c("Value", "group")
####

point_extract_liwonde<- as.numeric(point_extract_liwonde$value)
point_extract_liwonde<- as.data.frame(point_extract_liwonde)
point_extract_liwonde$group<- "cwatm_liwonde"
colnames(point_extract_liwonde)<- c("Value", "group")

comparisondf<- rbind(point_extract_liwonde, measured_discharge_output_liwonde)

start<- mondate("1971-01-01")
end<- mondate("2010-12-31")
comparisondf$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))
colnames(comparisondf)<- c("Value", "group", "date")

ggplot(comparisondf, aes(x= date, y= Value, colour=group)) + geom_line() 

###### Metrics for comparison
KGE <- function(o,s) {
  nans <- is.na(o)
  o <- o[!nans]
  s <- s[!nans]
  B <- mean(s, na.rm = TRUE) / mean(o, na.rm = TRUE)
  y <- (sd(s, na.rm = TRUE) / mean(s, na.rm = TRUE)) / (sd(o, na.rm = TRUE) / mean(o, na.rm = TRUE))
  r <- cor(s, o, method = "pearson")
  
  kge <- 1 -((r - 1) ^ 2 + (B - 1) ^ 2 + (y - 1) ^ 2)^0.5
  
  return(kge)
  
}
NSE <- function(o,s, normalize = TRUE) {
  nse <- 1 - sum((o - s) ^ 2, na.rm = TRUE) / sum((o - mean(o, na.rm = TRUE)) ^ 2, na.rm = TRUE)
  if(normalize) {
    nse <- 1 / (2 - nse)
  }
  return(nse)
}
nRMSE = function(o, s){
  nrmse <- ((sum((o - s) ^ 2) / length(o)) ^ 0.5/ sd(o))
  return(nrmse)
}
nRMSE_f = function(o, s){
  nanrm <- !is.na(o) * ! is.na(s)
  o <- o[nanrm]
  s <- s[nanrm]
  DF_o <- o - mean(o, na.rm = TRUE)
  DF_s <- s - mean(s, na.rm = TRUE)
  
  nrmse <- (100 / sd(DF_o)) * (sum((DF_o - DF_s) ^ 2) / length(DF_o)) ^ 0.5
  return(nrmse)
}
################
# Run metrics for stations 
KGE_chiromo<- KGE(as.numeric(measured_discharge_output$Value),as.numeric(point_data$Value))

Warning: NAs introduced by coercion

NSE_chiromo<- NSE(as.numeric(measured_discharge_output$Value),as.numeric(point_data$Value))

Warning: NAs introduced by coercion

nRMSE_chiromo<- nRMSE(as.numeric(measured_discharge_output$Value),as.numeric(point_data$Value))

Warning: NAs introduced by coercion

nRMSE_f_chiromo<- nRMSE_f(as.numeric(measured_discharge_output$Value),as.numeric(point_data$Value))

Warning: NAs introduced by coercion

###
KGE_liwonde<- KGE(as.numeric(measured_discharge_output_liwonde$Value),as.numeric(point_extract_liwonde$Value))

NSE_liwonde<- NSE(as.numeric(measured_discharge_output_liwonde$Value),as.numeric(point_extract_liwonde$Value))

nRMSE_liwonde<- nRMSE(as.numeric(measured_discharge_output_liwonde$Value),as.numeric(point_extract_liwonde$Value))

nRMSE_f_liwonde<- nRMSE_f(as.numeric(measured_discharge_output_liwonde$Value),as.numeric(point_extract_liwonde$Value))

##Print metrics values for stations
print (c("Chiromo Station metrics", "KGE_chiromo", KGE_chiromo, "NSE_chiromo", NSE_chiromo, "nRMSE_chiromo", nRMSE_chiromo, "nRMSE_f_chiromo", nRMSE_f_chiromo ))

[1] "Chiromo Station metrics" "KGE_chiromo"             "-2.90500309340093"      
[4] "NSE_chiromo"             "0.0117622627186949"      "nRMSE_chiromo"          
[7] NA                        "nRMSE_f_chiromo"         "143.112938375522"       

####
print (c("liwonde Station metrics", "KGE_liwonde", KGE_liwonde, "NSE_liwonde", NSE_liwonde, "nRMSE_liwonde", nRMSE_liwonde, "nRMSE_f_liwonde", nRMSE_f_liwonde ))

[1] "liwonde Station metrics" "KGE_liwonde"             "-0.120170562108403"     
[4] "NSE_liwonde"             "0.175589755533874"       "nRMSE_liwonde"          
[7] NA                        "nRMSE_f_liwonde"         "96.9130044495098"       

Evapotranspiration comparisons.

Run comparisons of daily evapotranspiration generated in CWatM with various metrics of evapotranspiration in generated in Google Earth Engine. Align any overlap between the dates available for the CWatM run and the available evapotranspiration estimates.

library(lubridate)

## Import CWatM evapotranspiration summary
ET_summary_CwatM<- "CWatM_outputs_evapotranspiration_original.csv"
CWatM_average_cell_ET<- read.csv(sprintf("%s/%s", pth, ET_summary_CwatM))

##Calculate ET for average cell by dividing total by number of cells
CWatM_average_cell_ET_average<- (CWatM_average_cell_ET$totalET)*1000/(as.numeric(CWatM_average_cell_ET[1,1]))

CWatM_average_cell_ET<- cbind(CWatM_average_cell_ET, CWatM_average_cell_ET_average)
CWatM_average_cell_ET<- CWatM_average_cell_ET[,2:4]
colnames(CWatM_average_cell_ET)<- c("Date", "Total Evap", "Average cell Evap")

## Convert date from character to date format
CWatM_average_cell_ET$Date <-as.Date(CWatM_average_cell_ET$Date,
        format =  "%d/%m/%Y")

#Plot average ET 
ggplot(CWatM_average_cell_ET, aes(x= Date, y= `Average cell Evap`)) + geom_line()

##### Comapring CWATM ET estimations with Google Earth Engine calculations

# GLDAS Evapotranspiration
ET_GLDAS<- read.csv(sprintf("%s/%s", pth, "ET_GLDAS.csv"))
ET_GLDAS$system.time_start<- mdy(ET_GLDAS$system.time_start)
colnames(ET_GLDAS)<- c("Date", "Value")
ET_GLDAS_CWatM<- merge(CWatM_average_cell_ET,ET_GLDAS, by= "Date")
colnames(ET_GLDAS_CWatM)<- c("Date", "CWatM total", "CWatM Average", "GLDAS Value")
ggplot(ET_GLDAS_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `GLDAS Value`, colour="GLDAS Evapotranspiration")) 

# Terraclimate Evapotranspiration
ET_Terraclimate<- read.csv(sprintf("%s/%s", pth, "ET_Terraclimate.csv"))
ET_Terraclimate$system.time_start<- mdy(ET_Terraclimate$system.time_start)
colnames(ET_Terraclimate)<- c("Date", "Value")
ET_Terraclimate_CWatM<- merge(CWatM_average_cell_ET,ET_Terraclimate, by= "Date")
colnames(ET_Terraclimate_CWatM)<- c("Date", "CWatM total", "CWatM Average", "Terraclimate Value")
ggplot(ET_Terraclimate_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `Terraclimate Value`, colour="Terraclimate Evapotranspiration")) 

# MODIS_500m Evapotranspiration
ET_MODIS_500m<- read.csv(sprintf("%s/%s", pth, "ET_MODIS_500m.csv"))
ET_MODIS_500m$system.time_start<- mdy(ET_MODIS_500m$system.time_start)
colnames(ET_MODIS_500m)<- c("Date", "Value")
ET_MODIS_500m_CWatM<- merge(CWatM_average_cell_ET,ET_MODIS_500m, by= "Date")
colnames(ET_MODIS_500m_CWatM)<- c("Date", "CWatM total", "CWatM Average", "MODIS_500m Value")
ggplot(ET_MODIS_500m_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `MODIS_500m Value`, colour="MODIS_500m Evapotranspiration")) 

# MODIS_1km Evapotranspiration
ET_MODIS_1km<- read.csv(sprintf("%s/%s", pth, "ET_MODIS_1km.csv"))
ET_MODIS_1km$system.time_start<- mdy(ET_MODIS_1km$system.time_start)
colnames(ET_MODIS_1km)<- c("Date", "Value")
ET_MODIS_1km_CWatM<- merge(CWatM_average_cell_ET,ET_MODIS_1km, by= "Date")
colnames(ET_MODIS_1km_CWatM)<- c("Date", "CWatM total", "CWatM Average", "MODIS_1km Value")
ggplot(ET_MODIS_1km_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `MODIS_1km Value`, colour="MODIS_1km Evapotranspiration")) 

# SMAP Evapotranspiration
# No overlapping dates for CWATM and SMAP as SMAP does not begin until after 2010 but enables comparisons for future runs
ET_SMAP<- read.csv(sprintf("%s/%s", pth, "ET_SMAP.csv"))
ET_SMAP$system.time_start<- mdy(ET_SMAP$system.time_start)
colnames(ET_SMAP)<- c("Date", "Value")
ET_SMAP_CWatM<- merge(CWatM_average_cell_ET,ET_SMAP, by= "Date")
colnames(ET_SMAP_CWatM)<- c("Date", "CWatM total", "CWatM Average", "SMAP Value")
ggplot(ET_SMAP_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `SMAP Value`, colour="SMAP Evapotranspiration")) 

# ETPML Evapotranspiration
ETPML<- read.csv(sprintf("%s/%s", pth, "ETPML.csv"))
ETPML$system.time_start<- mdy(ETPML$system.time_start)
colnames(ETPML)<- c("Date", "Ec", "Ei", "Es")
ETPML$Total<- (ETPML[,2]+ETPML[,3]+ETPML[,4])
ETPML_CWatM<- merge(CWatM_average_cell_ET,ETPML, by= "Date")
colnames(ETPML_CWatM)<- c("Date", "CWatM total", "CWatM Average", "ETPML Ec", "ETPML Ei", "ETPML Es", "ETPML Total")

###### Plot CWatM output alongside differnt components of ETPML Evapotranspiration estimate
ggplot(ETPML_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `ETPML Total`, colour="ETPML Total Evapotranspiration"))  + geom_line(aes(x= Date, y= `ETPML Ec`, colour="ETPML Ec")) + geom_line(aes(x= Date, y= `ETPML Ei`, colour="ETPML Ei"))  + geom_line(aes(x= Date, y= `ETPML Es`, colour="ETPML Es"))

#### Plot CWatM with only sum of different components of ETPML
ggplot(ETPML_CWatM) + geom_line(aes(x= Date, y= `CWatM Average`, colour="CWatM")) + geom_line(aes(x= Date, y= `ETPML Total`, colour="ETPML Total Evapotranspiration"))

With regards to evapotranspiration, the CWatM model appears to be well calibrated.

There is a question of whether the Google Earth Engine estimates account for lake evapotranspiration which is likely to constitute a significant component of the evapotranspiration total simulated by CWatM