diff --git a/08-hydrology1_precip.Rmd b/08-hydrology1_precip.Rmd index 93a1b72..9cf8652 100644 --- a/08-hydrology1_precip.Rmd +++ b/08-hydrology1_precip.Rmd @@ -2,6 +2,10 @@ All of the earlier chapters of this book dealt with the behavior of water in different hydraulic systems, such as canals or pipes. Now we consider the bigger picture of where the water originates, and ultimately how we can estimate how much water is available for different uses, and how much excess (flood) water systems will need to be designed and built to accommodate. +The _hydromisc_ package will need to be installed to access some of the data used below. If it is not installed, do so following the instructions on [the github site for the package](https://github.com/EdM44/hydromisc). + +## The hydrologic cycle + A fundamental concept is the hydrologic cycle, depicted in Figure \@ref(fig:hydr-cycle). (ref:fighc1) The hydrologic cycle, from the [USGS](https://labs.waterdata.usgs.gov/visualizations/water-cycle) @@ -10,9 +14,131 @@ A fundamental concept is the hydrologic cycle, depicted in Figure \@ref(fig:hydr knitr::include_graphics('images/USGS_WaterCycle_anthropogenic.png') ``` -The primary variable in the hydrologic cycle from an engineering perspective is precipitation, since that is the source of the water used and managed in engineered systems. +The primary variable in the hydrologic cycle from an engineering perspective is precipitation, since that is the source of the water used and managed in engineered systems. When precipitation hits the land surface it can take many paths: principally it may infiltrate into the soil or form runoff. Water collected in ponds or depressions may evaporate, and additional water may evaporate from the soil surface or be drawn from the soil by vegetation to be transpired. In a simplified way, the water budget may be expressed as Equation \@ref(eq:hc-1). +\begin{equation} + Q = P - E + \Delta S + (\#eq:hc-1) +\end{equation} +where _Q_ is outflow (as in streamflow volume), _P_ is precipitation, _E_ is evapotranspiration, and _ΔS_ is the change in storage (in snowpack, the soil or groundwater). Over long time periods, _ΔS_ becomes very small compared to the other terms, so for what is called ‘climatological’ periods (often 30-year averages) the water budget can simplified to Equation \@ref(eq:hc-2) +\begin{equation} + Q = P - E + (\#eq:hc-2) +\end{equation} +All variables must have the same units, so typically _Q_ is expressed as a depth of water (or runoff) over the contributing drainage area. Example \@ref(exm:ex-hcx1) shows a demonstration of a simple water budget over a basin. -The _hydromisc_ package will need to be installed to access some of the data used below. If it is not installed, do so following the instructions on [the github site for the package](https://github.com/EdM44/hydromisc). +::: {.example #ex-hcx1} +Using only the generally available data P, T (air temperature), and Q, assess the water budget for the watershed contributing flow to the USGS stream gauge 11169500, Saratoga Creek, CA. +::: +First considering Q, streamflow for USGS stream gauges can be directly downloaded from the USGS either from [their website](https://dashboard.waterdata.usgs.gov/app/nwd/en/) or using the R package [dataRetrieval](https://cran.r-project.org/package=dataRetrieval). However, for many sites those streamflow observations will include effects of upstream impoundments (dams) or diversions, so will not represent the response of the watershed to precipitation. Another option is to use estimates of "naturalized" streamflows, which have been calculated by the USGS as part of its [National Hydrographic Dataset (NHD)](https://www.usgs.gov/national-hydrography/about-national-hydrography-products) products. The R package [nhdplusTools](https://cran.r-project.org/package=nhdplusTools) provides one way to access this. The packages [sf](https://cran.r-project.org/package=sf) and [terra](https://cran.r-project.org/package=terra) will be useful for manipulating some spatial objects created along the way. + +First the gauge and contributing basin are downloaded using [nhdplusTools](https://cran.r-project.org/package=nhdplusTools), and the river reach flowlines are identified. +```{r message=FALSE} +library(nhdplusTools) + +nwissite <- list(featureSource = "nwissite", featureID = "USGS-11169500") + +# Retrieve gauge location and basin boundary objects +gauge <- get_nldi_feature(nwissite) +basin <- get_nldi_basin(nwissite) +nldi_feature <- list(featureSource = "comid", featureID = as.integer(gauge$comid)[1]) +flowline_nldi <- navigate_nldi(nldi_feature, + mode = "upstreamTributaries", + distance_km = 1000) +``` + +With the basin and reaches defined, the nhdplus data for the basin contributing to the outlet gauge can be downloaded as a geopackage (.gpkg file). Here the current working directory is set as where the file will be written. +```{r message=FALSE, eval=FALSE} +# Download the geopackage to the working directory +output_file_download <- file.path(getwd(), "subset_download.gpkg") + +nhd_subset <-subset_nhdplus(comids = as.integer(flowline_nldi$UT$nhdplus_comid), + output_file = output_file_download, + flowline_only = FALSE, + nhdplus_data = "download", + return_data = TRUE, overwrite = TRUE) +flowline <- nhd_subset$NHDFlowline_Network +``` + +```{r message=FALSE, echo=FALSE} +ndhsubsetfile <- system.file("extdata", "subset_download.gpkg", package="hydromisc") +flowline <- sf::st_read(ndhsubsetfile, "NHDFlowline_Network", quiet = TRUE) +``` + +In the example above, _nhd_subset_ has four layers of _sf_ objects, one of which is the _NHDFlowline_ layer that is used in what follows. What is displayed in Figure \@ref(fig:wbfig) is a static image produced using the [mapview](https://cran.r-project.org/package=mapview) function _mapshot_. +```{r, message=FALSE, eval=FALSE} +mapview::mapview(basin, legend=FALSE) + + mapview::mapview(flowline, color = "cyan", col.regions = "white", lwd = 3, legend=FALSE) + + mapview::mapview(gauge, color = "red", col.regions = "white", cex = 3, legend=FALSE) +``` + +```{r wbfig, message=FALSE, echo=FALSE, fig.align = 'center', out.width = "70%", fig.cap = "Watershed boundary and nhdplus stream network for USGS gauge 11169500."} +knitr::include_graphics('images/nhdplus_basin.png') +``` + +In this example the flows with the prefix *qc* (representing the QCMA flows described in the [NHDPlus User's Guide](https://pubs.usgs.gov/publication/ofr20191096)) are extracted to represent naturalized flows. These are also converted to _mm_ using the total basin drainage area and the number of days per month. +```{r message=FALSE} +# Extract the data for only the outlet reach +outlet <- subset(flowline, flowline$comid == gauge$comid) +drainage_area_km2 <- outlet$totdasqkm +# 12 monthly flows for outlet - QC is natural flows adjusted using unimpaired gauges +df.flows <- dplyr::select(sf::st_drop_geometry(outlet), starts_with("qc_"), -qc_ma) + +df.flows <- df.flows |> tidyr::pivot_longer(cols = 1:12, names_to = "Month", values_to = "flow_cfs") |> + dplyr::mutate(Month = month.abb) +# Add a column that converts monthly flow in cfs into a depth in mm over the watershed +dpm <- lubridate::days_in_month(1:12) +df.flows$runoff <- df.flows$flow_cfs * 86400 * dpm * 304.8^3 / (drainage_area_km2 * 10^12) +``` + +To obtain P and T there are many options. A simple one is to use an online tool such as [Model My Watershed](https://modelmywatershed.org/) and read in the csv file of P and T. This example will obtain 30-year normals for P and T using the [prism R package](https://cran.r-project.org/package=prism). This is extracted from the [PRISM](https://prism.oregonstate.edu/normals/) archive, a widely used product for P and T in the U.S. The [terra](https://cran.r-project.org/package=terra) package is used to extract the portion of the layer within the basin and find the mean. +```{r message=FALSE, eval=FALSE} +library(prism) + +prism_set_dl_dir(tempdir()) +# Download the monthly 30-year averages at 4km resolution +get_prism_normals(type="tmean", resolution = "4km", mon = 1:12, keepZip = FALSE) +flist <- prism_archive_subset("tmean", "monthly normals", mon = 1:12, resolution = "4km") +flist_full <- pd_to_file(flist) +rtmp <- terra::rast(flist_full) +df.t_prism <- terra::extract(rtmp, basin, fun = mean) |> + tidyr::pivot_longer(cols = c(-ID), names_to = "Month", values_to = "tavg_C") |> + dplyr::mutate(Month = month.abb) |> dplyr::select(-c("ID")) + +get_prism_normals(type="ppt", resolution = "4km", mon = 1:12, keepZip = FALSE) +flist <- prism_archive_subset("ppt", "monthly normals", mon = 1:12, resolution = "4km") +flist_full <- pd_to_file(flist) +rppt <- terra::rast(flist_full) +df.p_prism <- terra::extract(rppt, basin, fun = mean) |> + tidyr::pivot_longer(cols = c(-ID), names_to = "Month", values_to = "precip") |> + dplyr::mutate(Month = month.abb) |> dplyr::select(-c("ID")) +``` + +```{r message=FALSE, echo=FALSE} +df.t_prism <- hydromisc::df.t_prism +df.p_prism <- hydromisc::df.p_prism +``` + +Finally, combine the runoff, P and T into a single data frame and plot them. The potential evapotranspiration (PET) is added using a function that only requires temperature and latitude (for day length) from the [SPEI](https://cran.r-project.org/package=SPEI) package. +```{r wbprism, message=FALSE, fig.align = 'center', out.width = "70%", fig.cap = "Partial water budget for the basin contributing to USGS gauge 11169500."} +df.wb <- list(df.flows, df.t_prism, df.p_prism) |> purrr::reduce(dplyr::left_join, by = "Month") +# Find basin centroid +pt1 <- sf::st_coordinates(sf::st_centroid(basin)) +df.wb$PET <- SPEI::thornthwaite(df.wb$tavg_C, lat = pt1[2], verbose=FALSE) + +df.plot <- df.wb |> dplyr::select(Month, runoff, PET, precip) |> + tidyr::pivot_longer(cols = c(runoff, PET), names_to = "Component", values_to = "Flux") + +# Convert months to integers for plotting +suppressPackageStartupMessages(library(ggplot2)) +ggplot(data = df.plot, aes(x = match(Month,month.abb))) + + geom_col(aes(y = Flux, fill = Component)) + + geom_line(aes(y = precip, lty = "precip"), linewidth = 1.25) + + scale_x_discrete(limits = month.abb) + + labs(x="",y="Precip, Runoff or PET, mm", linetype = NULL) + + theme_bw() +``` + +While incomplete, this water budget illustrates months in which the basin is water limited (runoff + PET > precipitation) or energy limited (runoff + PET < precipitation). Water is more available to recharge the soil moisture and groundwater during energy limited periods since actual evaporation, which cannot exceed PET, will leave excess water, following Equation \@ref(eq:hc-1). As another alternative, a simple water budget can be estimated from P and T and an estimate of the soil moisture holding capacity using the [bioclim](https://cran.r-project.org/package=bioclim) R package. ## Precipitation observations @@ -93,7 +219,6 @@ knitr::kable(monavg, digits = 2) |> The winter of 2016-2017 (water year 2017) was a record wet year for much of California. Figure \@ref(fig:prcp-2) shows a **hyetograph** the daily values for that year. ```{r prcp-2, message=FALSE, warning = FALSE, fig.align = 'center', out.width = "75%", fig.cap="Daily Precipitation for San Jose, CA for water year 2017"} -library(ggplot2) ghcn_prcp2 <- data.frame(date = ghcn_data$DATE, wy = ghcn_data$wy, prcp = ghcn_data$PRCP ) ggplot(subset(ghcn_prcp2, wy==2017), aes(x=date, y=prcp)) + geom_bar(stat="identity",color="red") + @@ -111,7 +236,7 @@ avgrainrate <- ghcn_prcp2[ghcn_prcp2$prcp > 0.01,] |> group_by(wy) |> nraindays <- ghcn_prcp2[ghcn_prcp2$prcp > 0.01,] |> group_by(wy) |> summarise(nraindays = length(prcp)) #Find length of consecutive dry and wet spells for the record -days.dry.wet <- seas::interarrival(ghcn_prcp, var = "prcp", p.cut = 0.01, inv = FALSE) +days.dry.wet <- seas::interarrival(ghcn_prcp2, var = "prcp", p.cut = 0.01, inv = FALSE) #add a water year column to the result days.dry.wet$wy <- sapply(days.dry.wet$date, wateryr) res <- days.dry.wet |> group_by(wy) |> summarise(cdd = median(dry, na.rm=TRUE), cwd = median(wet, na.rm=TRUE)) @@ -199,7 +324,6 @@ annual_sum <- data.frame(year = annual_data$Year, sum_X = cumsum(annual_data$Station_X)) #create scatterplot with a label on every point -library(ggplot2) library(ggrepel) ggplot(annual_sum, aes(sum_X,sum_A, label = year)) + geom_point() + diff --git a/09-hydrology2_processes.Rmd b/09-hydrology2_processes.Rmd index 7ae54ed..1b574ab 100644 --- a/09-hydrology2_processes.Rmd +++ b/09-hydrology2_processes.Rmd @@ -207,7 +207,7 @@ Here $E_p$ is in mm/d, $\Delta$ and $\gamma$ are in $kPa~K^{-1}$, $R_n$ is in $M Open water evaporation can also be calculated using a modified version of the Penman-Monteith equation \@ref(eq:pm1). In this latter case, vegetation coefficients are not needed, so Equation \@ref(eq:pm1) can be used with $r_s=0$ and $r_a=251/(1+0.536~u_2)$, following [Thom & Oliver, 1977](https://doi.org/10.1002/qj.49710343610). -The R package [Evaporation](https://cran.r-project.org/package=Evapotranspiration) has functions to calculate $ET_0$ using this and many other functions. This is especially useful when calculating PET over many points or through a long time series. +The R package [Evaporation](https://cran.r-project.org/package=Evapotranspiration) has functions to calculate $ET_0$ using this and many other functions. This is especially useful when calculating PET over many points or through a long time series. Satellite-observed ET has been estimated globally as part of the [OpenET project](https://etdata.org/); the data may be accessed using the [openet R package](https://github.com/codeswitching/openet). ## Snow diff --git a/docs/designing-for-floods-flood-hydrology.html b/docs/designing-for-floods-flood-hydrology.html index 6be2d40..fb8e8ce 100644 --- a/docs/designing-for-floods-flood-hydrology.html +++ b/docs/designing-for-floods-flood-hydrology.html @@ -25,7 +25,7 @@ - + @@ -256,10 +256,11 @@
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation

    Start with some assumed values and run the snow model.

    -
    Tmax_snow <- 3  
-Tmin_rain <- 2
-kd <- 1
-Tbase <- 0
-snow_estim <- hydromisc::snow.sim(DATA = snow_data_sel, 
-                                  Tmax = Tmax_snow, 
-                                  Tmin = Tmin_rain,
-                                  kd = kd,
-                                  Tmelt = Tbase)
    +
    Tmax_snow <- 3  
+Tmin_rain <- 2
+kd <- 1
+Tbase <- 0
+snow_estim <- hydromisc::snow.sim(DATA = snow_data_sel, 
+                                  Tmax = Tmax_snow, 
+                                  Tmin = Tmin_rain,
+                                  kd = kd,
+                                  Tmelt = Tbase)

    Now the simulated values can be compared to the observations. If not installed already, install the hydroGOF package, which has some useful functions for evaluating how well modeled output fits observations. In the plot that follows we specify three measures of goodness-of-fit:

    These are discussed in detail in other references, but the aim is to calibrate (change the input parameters) until these values are low.

    -
    obs <- snow_data_sel$snow_water_equivalent
-sim <- snow_estim$swe_simulated
-hydroGOF::ggof(sim, obs, na.rm = TRUE, dates=snow_data_sel$date,
-               gofs=c("MAE", "RMSE", "PBIAS"),
-               xlab = "", ylab="SWE, mm", 
-               tick.tstep="months", cex=c(0,0),lwd=c(2,2))
    +
    obs <- snow_data_sel$snow_water_equivalent
+sim <- snow_estim$swe_simulated
+hydroGOF::ggof(sim, obs, na.rm = TRUE, dates=snow_data_sel$date,
+               gofs=c("MAE", "RMSE", "PBIAS"),
+               xlab = "", ylab="SWE, mm", 
+               tick.tstep="months", cex=c(0,0),lwd=c(2,2))
    Simulated and Observed SWE at SNOTEL site 1050 for Winter 2008-2009.

    @@ -655,48 +656,48 @@

    9.4.2 Basic snowmelt theory and s

    9.4.3 Snow model calibration

    While manual model calibration can improve the fit, a more complete calibration involves optimization methods that search the parameter space for the optimal combination of parameter values. The R function optim, part of the stats package installed with base R, is efficient at finding a local minimum or maximum to a function. For most hydrologic models finding the global optimum means searching the entire parameter space, since many local minima and maxima usually exist. Some version of Monte Carlo methods are used to do this.

    There are many R packages to perform different variations of Monte Carlo simulations, but a simple implementation can be done manually using the foreach package. In this example the four main model parameters are each varied over 10 evenly spaced values, producing \(10^4\) unique parameters sets that need to be run. For each run, the function will return a value for RMSE.

    -
    # package needed for efficiently looping through all simulations
-library(foreach)
-
-# Create the sets of values of each parameter within a defined range
-p_Tmax <- seq(-1, 3, length.out = 10)
-p_Tmin <- seq(-1, 2, length.out = 10)
-p_kd <- seq(0.5, 5, length.out = 10)
-p_Tbase <- seq(-2, 2, length.out = 10)
-
-# Assemble them into a data frame with all parameter combinations
-p_all <- expand.grid(Tmax = p_Tmax, Tmin = p_Tmin, kd = p_kd, Tmelt = p_Tbase)
-
-# Define the function to return a single test statistic for each simulation
-fcn <- function(datain, Tmax, Tmin, kd, Tmelt){
-  snow_estim <- hydromisc::snow.sim(DATA=datain, Tmax=Tmax, Tmin=Tmin, kd=kd, Tmelt=Tmelt)
-  calib.stats <- suppressWarnings(hydroGOF::gof(snow_estim$swe_simulated,obs,na.rm=TRUE))
-  return(as.numeric(calib.stats['RMSE',]))
-}
    +
    # package needed for efficiently looping through all simulations
+library(foreach)
+
+# Create the sets of values of each parameter within a defined range
+p_Tmax <- seq(-1, 3, length.out = 10)
+p_Tmin <- seq(-1, 2, length.out = 10)
+p_kd <- seq(0.5, 5, length.out = 10)
+p_Tbase <- seq(-2, 2, length.out = 10)
+
+# Assemble them into a data frame with all parameter combinations
+p_all <- expand.grid(Tmax = p_Tmax, Tmin = p_Tmin, kd = p_kd, Tmelt = p_Tbase)
+
+# Define the function to return a single test statistic for each simulation
+fcn <- function(datain, Tmax, Tmin, kd, Tmelt){
+  snow_estim <- hydromisc::snow.sim(DATA=datain, Tmax=Tmax, Tmin=Tmin, kd=kd, Tmelt=Tmelt)
+  calib.stats <- suppressWarnings(hydroGOF::gof(snow_estim$swe_simulated,obs,na.rm=TRUE))
+  return(as.numeric(calib.stats['RMSE',]))
+}

    Now the function needs to be used to run the snow model for every parameters combination. What is shown below is how to do this using the foreach package in a brute force manner, running all simulations on a single CPU. This takes a while.

    -
    out <- foreach(i = 1:nrow(p_all), .combine='c') %do% {
-  fcn(snow_data_sel, p_all[i, "Tmax"], p_all[i, "Tmin"], p_all[i, "kd"], p_all[i, "Tmelt"])
-}
    +
    out <- foreach(i = 1:nrow(p_all), .combine='c') %do% {
+  fcn(snow_data_sel, p_all[i, "Tmax"], p_all[i, "Tmin"], p_all[i, "kd"], p_all[i, "Tmelt"])
+}

    Because the number of simulations can become very large, this command could be adapted to run the simulations in parallel using multiple cores by using the R package doParallel, which loads the parallel package (part of the base R installation). By running the foreach simulations with multiple processes, the run time using 7 of 8 available cores dropped from 5.9 minutes for a single CPU to 1.1 minutes. Aside from setting up the cluster before the command and stopping it after (essential for freeing up resources), the only change is the use of %dopar% instead of %do%.

    -
    library(doParallel)  
-no_cores <- detectCores() - 1  
-registerDoParallel(cores=no_cores)  
-cl <- makeCluster(no_cores, type="PSOCK")  # On mac or linux you can change type to "FORK"
-out <- foreach(i = 1:nrow(p_all), .combine='c') %dopar% {
-  fcn(snow_data_sel, p_all[i, "Tmax"], p_all[i, "Tmin"], p_all[i, "kd"], p_all[i, "Tmelt"])
-}
-stopCluster(cl)
    +
    library(doParallel)  
+no_cores <- detectCores() - 1  
+registerDoParallel(cores=no_cores)  
+cl <- makeCluster(no_cores, type="PSOCK")  # On mac or linux you can change type to "FORK"
+out <- foreach(i = 1:nrow(p_all), .combine='c') %dopar% {
+  fcn(snow_data_sel, p_all[i, "Tmax"], p_all[i, "Tmin"], p_all[i, "kd"], p_all[i, "Tmelt"])
+}
+stopCluster(cl)

    It can be helpful to look at the four parameters for the simulations that had the lowest RMSE. This presumes the tidyverse package has been installed.

    -
    # Combine results with parameters, sort by RMSE (ascending) 
-out_all <- cbind(p_all, RMSE=out) |>  dplyr::arrange(RMSE)
-
-library(ggplot2)
-# Plot the parameters for the 10 best runs
-dplyr::slice(out_all, 1:10) |> 
-  tidyr::pivot_longer(-RMSE, names_to = 'param', values_to = 'value') |>
-  ggplot(aes(x = RMSE, y = value)) + 
-  geom_point() +
-  facet_wrap(vars(param), scales = "free")
    +
    # Combine results with parameters, sort by RMSE (ascending) 
+out_all <- cbind(p_all, RMSE=out) |>  dplyr::arrange(RMSE)
+
+library(ggplot2)
+# Plot the parameters for the 10 best runs
+dplyr::slice(out_all, 1:10) |> 
+  tidyr::pivot_longer(-RMSE, names_to = 'param', values_to = 'value') |>
+  ggplot(aes(x = RMSE, y = value)) + 
+  geom_point() +
+  facet_wrap(vars(param), scales = "free")
    Parameter values for the 10 simulations with the lowest RMSE.

    @@ -704,18 +705,18 @@

    9.4.3 Snow model calibration

    There are several occurrences of equal RMSE values, so some points are superimposed on one another. Some parameter values are consistent for all of the 10 best simulations: kd is 2.5 and Tmelt is 2\(^\circ\)C, but optimal values can be achieved with a variety of different combinations of Tmax and Tmin. The results using the optimal parameters can be plotted to visualize the simulation.

    -
    Tmax_opt <- max(out_all$Tmax[1],out_all$Tmin[1])
-Tmin_opt <- min(out_all$Tmax[1],out_all$Tmin[1])
-kd_opt <- out_all$kd[1]
-Tmelt_opt <- out_all$Tmelt[1]
-
-snow_estim_opt <- hydromisc::snow.sim(DATA=snow_data_sel,Tmax=Tmax_opt, 
-                                      Tmin=Tmin_opt, kd=kd_opt, Tmelt=Tmelt_opt)
-sim_opt <- snow_estim_opt$swe_simulated
-hydroGOF::ggof(sim_opt, obs, na.rm = TRUE, dates=snow_data_sel$date,
-               gofs=c("MAE", "RMSE", "PBIAS"),
-               xlab = "", ylab="SWE, mm", 
-               tick.tstep="months", cex=c(0,0),lwd=c(2,2))
    +
    Tmax_opt <- max(out_all$Tmax[1],out_all$Tmin[1])
+Tmin_opt <- min(out_all$Tmax[1],out_all$Tmin[1])
+kd_opt <- out_all$kd[1]
+Tmelt_opt <- out_all$Tmelt[1]
+
+snow_estim_opt <- hydromisc::snow.sim(DATA=snow_data_sel,Tmax=Tmax_opt, 
+                                      Tmin=Tmin_opt, kd=kd_opt, Tmelt=Tmelt_opt)
+sim_opt <- snow_estim_opt$swe_simulated
+hydroGOF::ggof(sim_opt, obs, na.rm = TRUE, dates=snow_data_sel$date,
+               gofs=c("MAE", "RMSE", "PBIAS"),
+               xlab = "", ylab="SWE, mm", 
+               tick.tstep="months", cex=c(0,0),lwd=c(2,2))
    Simulated and Observed SWE at SNOTEL site 1050 for Winter 2008-2009 using the parameter set producing the lowest RMSE.

    @@ -728,19 +729,19 @@

    9.4.3 Snow model calibration

    9.4.4 Estimating climate change impacts on snow

    Once a reasonable calibration is obtained, the effect of increasing temperatures on SWE can be simulated by including the deltaT argument in the hydromisc::snow.sim command. Here a 3\(^\circ\)C uniform temperature increase is imposed on the optimal parameterization above.

    -
    dT <- 3.0 
-snow_plus3 <- hydromisc::snow.sim(DATA=snow_data_sel,
-                      Tmax=Tmax_opt, Tmin=Tmin_opt, kd=kd_opt, Tmelt=Tmelt_opt,
-                      deltaT = dT)
-simplusdT <- snow_plus3$swe_simulated
-# plot the results
-dTlegend  <- expression("Simulated"*+3~degree*C)
-plot(as.Date(snow_data_sel$date),obs,type = "l",xlab = "", ylab = "SWE (mm)")
-lines(as.Date(snow_estim$date),sim_opt,lty=2,col="blue")
-lines(as.Date(snow_estim$date),simplusdT,lty=3,col="red")
-legend("topright", legend = c("Observed", "Simulated",dTlegend),
-       lty = c(1,2,3), col=c("black","blue","red"))
-grid()
    +
    dT <- 3.0 
+snow_plus3 <- hydromisc::snow.sim(DATA=snow_data_sel,
+                      Tmax=Tmax_opt, Tmin=Tmin_opt, kd=kd_opt, Tmelt=Tmelt_opt,
+                      deltaT = dT)
+simplusdT <- snow_plus3$swe_simulated
+# plot the results
+dTlegend  <- expression("Simulated"*+3~degree*C)
+plot(as.Date(snow_data_sel$date),obs,type = "l",xlab = "", ylab = "SWE (mm)")
+lines(as.Date(snow_estim$date),sim_opt,lty=2,col="blue")
+lines(as.Date(snow_estim$date),simplusdT,lty=3,col="red")
+legend("topright", legend = c("Observed", "Simulated",dTlegend),
+       lty = c(1,2,3), col=c("black","blue","red"))
+grid()
    Observed SWE and simulated with observed meteorology and increased temperatures.

    diff --git a/docs/flow-in-open-channels.html b/docs/flow-in-open-channels.html index 49f7cf8..3b3b951 100644 --- a/docs/flow-in-open-channels.html +++ b/docs/flow-in-open-channels.html @@ -25,7 +25,7 @@ - + @@ -256,10 +256,11 @@

  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation
  • 8 The hydrologic cycle and precipitation
  • 9 Fate of precipitation