The authors have made significant changes to the manuscript. I am satisfied with the authors‘ responses as well as the changes. I think the quality of the manuscript was significantly improved. In general, I only have one major comment regarding the F*yw and the age threshold of F*yw (please see below).
Major comments
As I understood from the manuscript. The authors use F*yw from 22 Swiss catchments from Freyberg et al. (2018). I am not sure if the authors and Freyberg et al. (2018) used the same code (method) for the calculation F*yw or not. In the attached data along with this review, I demonstrated that using different source codes (methods) could result in different values of Fyw (or F*yw) and the age thresholds of F*yw. If so, the differences in F*yw among catchments could be driven by using different methods (source code) rather than the physical factors (e.g., precipitation, elevation, LFD, Fbf,...). Therefore, I would suggest the authors to use the same source code (method) for calculating F*yw for all catchments and provide the source code for review/or make it publicly accessible.
When applying the sine wave fitting, the authors will get F*yw corresponding to a specific age threshold (τyw), depending on the shape factor (alpha) of the gamma distribution (Equations 13-14 from Kirchner 2016). The age threshold could be 1 to 3 months (Figure 10, Kirchner 2016) and even beyond this range. Therefore, the variation in F*yw might be because of different age thresholds rather than physical factors. In this sense, the F*yw or Fyw might not be a useful matrix for catchment intercomparison studies. Therefore, I suggest showing the alpha and age thresholds for each catchment in the results (these values were often not shown in previous studies). Discuss the consequences of the differences found with the age threshold in the limitation section of the paper (e.g., will the relations between F*yw and other factors (e.g., precipitation, elevation, LFD, Fbf,...) change if somehow we can calculate F*yw that corresponds to the same age threshold).
Minor comments
Figure 2c: Might be using a triangle (for P) and a circle (for Q) in combination with different colors (for rain, hybrid, and snow) for easier to differentiate between them
L237-238: Should be noted that this demonstration uses thought experiments (not from real catchments)
Figure 6: Why did the authors only show quaternary deposits for 7 catchments, I think either showing for all catchments or not showing the figure or not showing (Table 2 is sufficient)
L207-211: “Our data set includes NBPV, whose area is 42% glaciercovered and consequently shows a characteristic glacier-dominated streamflow regime with a monthly peak in late summer. Thus, NBPV may belong to a fourth category of glacier-dominated catchments, for which the effect of glacier-melt on F*yw cannot be neglected, and this was partially discussed by Ceperley et al. (2020)”: Why did the author still classify NBPV as “snow-dominated” catchment? Does glacier-melt account for when calculating F*yw?
Section 4.1. Why did the authors discuss the Winter Flow Index (WFI) - F*yw relation in this section “The role of Quaternary deposits”?
Section 4.2. To be consistent, I would suggest using another title, e.g., “the role of groundwater flow (baseflow) in F*yw”
L414: I think the author mentions the fraction of based flow, not the amount of stored water
#-------------------------------------
#Rscript for Fyw and alpha
#1. How to run this script
# Copy this code to a text file and save as "youngfraction.R"
# Then open this file in R
# All of the required functions to run this Rsript was attached at the end of this file (please run these function before running this code)
#2. Download isotope data from the below link and extract in the same folder with the file "youngfraction.R"
# https://zenodo.org/record/4057967#.Y00oMHZBxPY
# then you will have two folders "Precipitation" and "isotopes_streamflow"
#3. The IRLS_hess-2017-720.R file was downloaded from
#https://hess.copernicus.org/articles/22/3841/2018/hess-22-3841-2018-supplement.zip
#I made 3 modifications in this function (please search with keywords "original code", where modifications were made)
#-------------------------------------
# Load library
library(lubridate)
library(lhs)
# set working directory - PLEASE CHANGE TO YOUR WORKING DIRECTORY --
setwd("./isotope_data_sine_fitting")
#-------------------------------------------------------------------------------
# 1. Fit sine wave to isotope in precipitation and streamflow (not using IRLS method)
# 20k parameter sets were generated within a given range
# best parameter set was selected based on weighted root mean square error
# this procedure is repeated n times to have the uncertainty in Fyw and alpha
#-------------------------------------------------------------------------------
# Read isotope in precipitation and convert date to decimal
isotope_P <- read.csv("./Precipitation/RIA.isoPrcp", header = TRUE, sep = "", skip = 12)
isotope_P$date <- toDecimal(as.Date(isotope_P$date, format = "%Y-%m-%d"))
# Read isotope in streamflow and convert date to decimal
isotope_Q <- read.csv("./isotopes_streamflow/RIA.isoStrm", header = TRUE, sep = "", skip = 22)
isotope_Q$Sampling_date <- toDecimal(as.Date(isotope_Q$Sampling_date, format = "%Y-%m-%d"))
# I cannot find Q data, now calculate with weight is = 1, in other words, flow non-weighted Fyw
weight_P <- rep(1, nrow(isotope_P))
weight_Q <- rep(1, nrow(isotope_Q))
# Initialize result
Fyw <- c()
alpha <- c()
# Make n iterations to get the uncertainty range of Fyw and alpha
n = 5
for (iter in 1:n){
# Generate parameter set
nsim <- 20000
parameter <- randomLHSrange(data.frame(
A = c(0,15),
phi = c(0,2*pi),
k = c(-14,-5)),
nsim)
# Fit isotope in P to sinewave - Plot result - parameter(A, phi, k)
bestParameterP <- sinefitbest(parameter, isotope_P$d18O_prcp, isotope_P$date, weight_P)
# plot(isotope_P$d18O_prcp) lines(isotope_Q$d18O_mean, add = TRUE)
# lines(sinefit(bestParameterP[1], bestParameterP[2], bestParameterP[3], isotope_P$date), col = "blue", lwd = 3)
# Fit istope data in Q
bestParameterQ <- sinefitbest(parameter, isotope_Q$d18O_mean, isotope_Q$Sampling_date, weight_Q)
# plot(isotope_Q$d18O_mean)
# lines(sinefit(bestParameterQ[1], bestParameterQ[2], bestParameterQ[3], isotope_Q$Sampling_date), col = "blue", lwd = 3)
# Find young water fraction and alpha value of the gamma distribution
Fyw <- c(Fyw, bestParameterQ[1]/bestParameterP[1])
alpha <- c(alpha, findAlpha(bestParameterQ,bestParameterP))
}
# Get the range of the estimated Fyw and alpha
summary(Fyw)
summary(alpha)
# Write result to result.csv file
write.csv(data.frame(Simulation = c(1:length(Fyw)), Fyw = Fyw, Alpha = alpha),
file = "Fyw_alpha.csv", row.names = FALSE)
#-------------------------------------------------------------------------------
# 2. Fit sine wave to isotope in precipitation and streamflow (using IRLS method)
# https://hess.copernicus.org/articles/22/3841/2018/hess-22-3841-2018-supplement.zip
#-------------------------------------------------------------------------------
# Fit sine wave to isotope concentrations in Q
fitQ <- IRLS(isotope_Q$d18O_mean,
data.frame(cos = cos(2*pi*isotope_Q$Sampling_date),sin = sin(2*pi*isotope_Q$Sampling_date)),
weight_Q,
type="Cauchy")
AQ <- sqrt(sum(fitQ$fit$coefficients[2:3]^2))
phiQ <- atan(fitQ$fit$coefficients[3]/fitQ$fit$coefficients[2])
# # Fit sine wave to isotope concentrations in P
fitP <- IRLS(isotope_P$d18O_prcp,
data.frame(cos = cos(2*pi*isotope_P$date),sin = sin(2*pi*isotope_P$date)),
weight_P,
type="Cauchy")
AP <- sqrt(sum(fitP$fit$coefficients[2:3]^2))
phiP <- atan(fitP$fit$coefficients[3]/fitP$fit$coefficients[2])
# Find Fyw
Fyw = AQ/AP
Fyw
# The phase shift must be positive and iwthin the range [0, 2*pi]
phiP <- phiP%%(2*pi)
phiQ <- phiQ%%(2*pi)
alpha <- findAlpha(c(AQ, phiQ%%(2*pi)), c(AP, phiP%%(2*pi)))
alpha
# Please run all functions after this line before running the main code
##############################################
#iteratively reweighted least squares (IRLS) with optional point weights
# (in addition to residual weights that are adjusted to downweight data with unusually large residuals,
# as in conventional IRLS)
#
#Author: James Kirchner, ETH Zurich
#Public use allowed under GNU Public License v. 3 (see http://www.gnu.org/licenses/gpl-3.0.en.html)
#
##############################################
IRLS <- function(Y, X, ww=rep(1,length(Y)), wt=rep(1,length(Y)), type="Cauchy") {
# Y is a numeric vector representing a response variable
# X is a numeric vector or array representing one or more explanatory variables
# ww is an optional numeric vector of point weights (such as masses for mass-weighted regressions)
# the default weight vector gives all ww's the value of 1
# wt is an optional numeric vector of initial values for the residual weights for the IRLS
# the default weight vector gives all initial wt's the value of 1
# type is an optional string constant indicating the weight function to use.
# If type is specified, it must be "bisquare", "Welsch", or "Cauchy". Anything else generates an error.
# the default weight function is Cauchy
# Y, ww and X (or each vector comprising X, if X is two-dimensional) must be of the same length, otherwise an error results.
# Y and X must not be exactly collinear (that is, there must be some nonzero residuals). This is not checked.
# IRLS returns an object of class "lm", generated by a call lm(Y ~ X, weights=wt, na.action="na.omit") with the iteratively determined weights.
# This object can then be handled just like any other return from lm.
# However, the variable names in this "lm" object will be Y and X (as defined internally here) regardless of what names were passed to this function!
##############################################
#define bisquare weight function
##############################################
bisquare <- function(x, MAR) {
if (MAR==0) w <- ifelse(x==0, 1, 0) #if MAR is zero, only 0 or 1 weights are possible
else {
w <- (1-(x/(6*MAR))^2)^2 #bisquare function
w[x>6*MAR] <- 0 #replace with zero whenever X is more than 6 times the median absolute residual
}
return(w)
}
##############################################
#define Welsch weight function
##############################################
Welsch <- function(x, MAR) {
if (MAR==0) w <- ifelse(x==0, 1, 0) #if MAR is zero, only 0 or 1 weights are possible
else w <- exp(-(x/(4.4255*MAR))^2) #Welsch weight function
return(w)
}
##############################################
#define Cauchy weight function
##############################################
Cauchy <- function(x,MAR) {
if (MAR==0) w <- ifelse(x==0, 1, 0) #if MAR is zero, only 0 or 1 weights are possible
else w <- 1/(1+(x/(3.536*MAR))^2) #Cauchy weight function
return(w)
}
if (type!="bisquare") {
if (type!="Welsch") {
if (type!="Cauchy") stop("IRLS stopped: no valid weight type specified. Valid functions are 'bisquare', 'Welsch', and 'Cauchy'")
}
}
if (length(Y) != length(ww)) stop("IRLS stopped: supplied weight vector must be same length as Y")
wwwt <- ww*wt
# original code fit <- lm(Y ~ X, weights=wwwt, na.action="na.omit")
fit <- lm(Y ~ X$cos + X$sin, weights=wwwt, na.action="na.omit") #initial least-squares fit
wt_chg <- 999.0 #initialize weight change
iter <- 0 #initialize the iteration counter
##############################################
#Here's the iteration loop, which runs until the largest weight change for any point is less than 0.01, or iteration limit is exceeded
##############################################
while ( (max(wt_chg,na.rm=TRUE) > 0.01) & !all(iter>10, summary(fit)$r.squared>0.999) ){
if (iter>1000) stop("IRLS stopped: more than 1000 interations, sorry!")
# This error can arise when Y is perfectly collinear with X for more than half the points,
# and thus the MAR fluctuates near zero, with the weights never stabilizing.
# That should normally be handled by the r-squared criterion for loop exiting as defined above.
iter <- iter+1 #increment the iteration counter
old_wt <- wt #save the old vector of weights for comparison with the next one
const <- fit$coefficients[1] #this is the intercept
slope <- fit$coefficients[2:length(fit$coefficients)] #this is the vector of regression coefficients
#explicit calculation of residuals
# original code if (length(slope)==1) resid <- Y - const - X*slope
if (length(slope)==1) resid <- Y - const - as.matrix(X)*slope
else resid <- as.vector(Y - const - as.matrix(X) %*% slope)
#can't use residuals(fit) because missing values will mess up the assignment of weights in the steps that follow
#note %*% is matrix multiplication in R
abs_resid <- abs(resid)
abs_resid_nonzero <- ifelse(Y==0, NA, abs_resid) #exclude residuals corresponding to exact zeroes from median, to avoid blowup when there are many repeated zeroes in any flow decile
MAR <- median(abs_resid_nonzero, na.rm=TRUE)
if (MAR==0.0) stop("IRLS stopped. Solution has collapsed: median absolute residual is zero!")
if (type=="bisquare") wt <- bisquare(abs_resid,MAR) #use bisquare weights
else if (type=="Welsch") wt <- Welsch(abs_resid,MAR) #use Welsch weights
else if (type=="Cauchy") wt <- Cauchy(abs_resid,MAR) #use Cauchy weights
wwwt <- ww*wt
# original code fit <- lm(Y ~ X, weights=wwwt, na.action="na.omit")
fit <- lm(Y ~ X$cos + X$sin, weights=wwwt, na.action="na.omit") #run multiple regression
wt_chg <- abs(wt-old_wt) #calculate change in weights from previous iteration
} #end while
qq <- list("fit"=fit, "wt"=wt)
return(qq)
}
##############################################
#END OF iteratively reweighted least squares
##############################################
##############################################
# The below Rscript was written by the reviewer
# This script contains functions for fitting the sinewave isotope data in streamflow and precipitation
# There is also functions for finding youngwater fraction and alpha value
##############################################
# Create a sine function describing isotope in precipitation or streamflow
sinefit <- function(A, phi,k,t){
sinefit <- A * sin(2*pi *t - phi) + k
return(sinefit)
}
# Find best fit parameter to sinewave fitting given parameter sets and objective function is RMSE
sinefitbest <- function(parameter, isotopeData, inputDate, weight){
modelPerformance <- c()
for (i in 1:nrow(parameter)){
modelPerformance <- c(modelPerformance,
weighted_least_square(isotopeData, sinefit(parameter[i,1], parameter[i,2], parameter[i,3], inputDate), weight))
}
#return best par
return(parameter[which(modelPerformance == min(modelPerformance)),])
}
# Model performance
weighted_least_square <- function(obs, sim, weight){
missing <- which(is.na(obs))
if(length(missing)>0){
weighted_least_square <- sqrt(mean(((obs[-missing] - sim[-missing])*weight)^2))
} else {
weighted_least_square <- sqrt(mean(((obs - sim)*weight)^2))
}
return(weighted_least_square)
}
# Convert date to decimal
toDecimal <- function(inputDate){
toDecimal <- decimal_date(as.Date(inputDate, format = "%Y-%m-%d"))
toDecimal <- toDecimal - trunc(toDecimal)
return(toDecimal)
}
# Latin hypercube with input range
randomLHSrange <- function(inputRange, n){
randomLHSrange <- randomLHS(n,ncol(inputRange))
for(i in 1:ncol(inputRange)){
randomLHSrange[,i] <- inputRange[1,i] + (inputRange[2,i] - inputRange[1,i])*randomLHSrange[,i]
}
return(randomLHSrange)
}
# find alpha value
findAlpha <- function(parameterQ, parameterP){
error <- c()
alpha <- seq(0.01,20,0.01)
for (i in 1:length(alpha)){
error <- c(error, (parameterQ[2] - parameterP[2]) - alpha[i] * atan(sqrt((parameterQ[1]/parameterP[1])^(-2/alpha[i]) -1)))
}
error <- abs(error)
return(alpha[which(error == min(error))])
} |
