## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----------------------------------------------------------------------------- library(DMSTAr) data(series) ## ----------------------------------------------------------------------------- # for example, limit input file series <- series[1:540,]; # Data formatting series$Qi <- cfs_to_hm3d(series$Flow) # cfs to hm3/d series$Rain <- in_to_m(series$Rainfall) # inches to meters per day series$Et <- in_to_m(series$ET) series$Zcontrol <- cm_to_m(0) # meters; setting to zero to see what happens # If you have release series; otherwise set to 0 series$Qr0 <- 0 # constrained outflow (forced Q) if used series$Qr1 <- 0 # release 1 series$Qr2 <- 0 # release 2 series$Ci <- series$Conc # input parameters params <- list( A_cell = 2.19, # cell area; km2 Zmin = 2, # minimum depth; cm Vmin = 0, # minimum volume; hm3 Q_a = 1.0, # discharge coef Q_b = 4.0, # discharge exponent Zweir = 0, # depth offset for outflow computation; cm Q_zmin = 38, # minimum depth of discharge; cm Qomax = 0.0, # maximum discharge hm3/day Qimax = 0, # maximum inflow (hm3/day) Width = 1.55, # cell widthl km Bypass_elev = 0, # mean depth at which bypass begins, cm Seepout_Rate = 0.00789, # outflow seepage rate per unit head; cm/d/cm Seepout_Elev = 0.0, # elevation controlling outflow seepage rate; cm Seepin_Rate = 0.0, # seepage inflow rate; ; cm/d/cm Seepin_Elev = 0.0, # elevation controlling inflow seepage rate; cm ShutdownET = TRUE, force_Q_out = FALSE, wrap_interp = TRUE, Zinit = 40, # initial water column depth; cm Qin_Frac = 0.22, # fraction of basin flows going into this cell Zrelease = 0, # minimum depth for releases; cm RecycleQ = 0 ) # initial volume (based on input values) V_init <- (cm_to_m(params$Zinit) * params$A_cell) out_hydro <- dmsta_flow_series( series = series, params = params, Nsteps = 4) hydro_rslt <- out_hydro$results # convert water depth from meter to centimeters hydro_rslt$Z_end_cm <-m_to_cm(hydro_rslt$Z_end) ## ----eval = FALSE------------------------------------------------------------- # str(out_hydro) ## ----hydro_plot, echo=FALSE, fig.align = "center",fig.cap = "Simulated outflow discharge (top) and water level (bottom).",fig.width = 7,fig.height = 5,fig.align="center"---- hydro_rslt$Date <- as.Date(hydro_rslt$Date) oldpar <- par(no.readonly = TRUE) par(mar=c(2,3,1,0.5),oma = c(2,2,0.5,0.5)) layout(matrix(1:2,2,1)) plot(Q_treated~Date,hydro_rslt,type="l",las=1,ann=F, ylim=c(0,0.4), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines( Qin~Date,hydro_rslt,col="dodgerblue1") mtext(side=3,adj=0.01,"WY 1966") mtext(side=3,adj=0.99,"Cell 1") mtext(side=2,line=3,expression("Outflow Q (hm"^"3"~"d"^"-1"*")")) legend("topleft",inset = c(0.01,0.01),legend=c("Inflow","Outflow"), pch=c(NA),lty=c(1),lwd=c(1), col=c("dodgerblue1",'black'),pt.bg=NA, pt.cex=1.5,ncol=1,cex=0.75,bty="n",y.intersp=1.25,x.intersp=0.5,xpd=NA,xjust=0.5,yjust=1) plot(Z_end_cm~Date,hydro_rslt,type="l",las=1,ann=F, ylim=c(min(hydro_rslt$Depth,hydro_rslt$Z_end_cm), max(hydro_rslt$Depth,hydro_rslt$Z_end_cm)), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) mtext(side=2,line=3,"Depth (cm)") mtext(side=1,line=2.5,"Date") par(oldpar) ## ----------------------------------------------------------------------------- ## Phosphorous Modeling parameters pparams <- list( DutyCycle = 0.95, Cmax = 2000, C1000 = 22, Cstar = 3, Ks_per_yr = 16.8, Z1 = 40, Z2 = 100, Z3 = 200, Chalf = 300, K2Coef1 = 0, Ytrans = 0, Ysigma = 0, Czero = 0, C_rain = 10, DryDepo = 20, SeasonalFactor = 0, C1000_2 = NULL, ks_2 = 0, zh_2 = 0, k_depth_penalty = 1, seepage_c = 20, seepin_conc = 0, C_init_ppb = 30, Y_init_mgm2 = 3387.67297548954, n_tanks = 3, Nsteps = 4 ) # Add the hydrology parameters (above) to the P parameters params <- modifyList(params,pparams) ## choose model structure ttankS <- params$n_tanks # tanks in series (can be fractional) Nsteps <- params$Nsteps # RK substeps per day ## build tank geometry once tanks <- dmsta_build_tanks(params$A_cell, ttankS) ## build kinetics once: 3 modules STA/PSTA/RES ppar <- build_P_kin_slots( mods = c("STA", "PSTA", "RES"), pparams = params, Dpy = 365.25, DutyCycle = params$DutyCycle) # A function to validation/check input parameters validate_P_paramsK(ppar) ## constants constants <- list( Cmax = params$Cmax, C_rain = params$C_rain, DryDepo = params$DryDepo / 365.25, # convert to mg/m2-day seepin_conc = params$seepin_conc, seepout_conc_max = params$seepage_c, fseep_recycle = 0, fseep_out = 0 ) ## initial conditions Z_init_m <- cm_to_m(params$Zinit) V_init <- params$A_cell * Z_init_m P_state0 <- dmsta_p_init_state( tanks, Z_init_m = Z_init_m, C_init_ppb = params$C_init_ppb, Y_init_mgm2 = params$Y_init_mgm2 ) hydroP_out <- dmsta_flowP_series( series = series, params = params, ttankS = ttankS, Nsteps = Nsteps, tanks = tanks, ppar = ppar, constants = constants, V_init = V_init, init_P_state = P_state0, return_steps = FALSE ) hydroP_rslt <- hydroP_out$results hydroP_rslt$Z_end_cm <-m_to_cm(hydroP_rslt$Z_end) ## ----eval = FALSE------------------------------------------------------------- # str(hydroP_out) ## ----hydroP_plot, echo=FALSE, fig.align = "center",fig.cap = "Simulated outflow discharge (top left), water level (bottom left), TP load (top right) and TP concentration (bottom right).",fig.width = 7,fig.height = 5---- oldpar <- par(no.readonly = TRUE) par(mar=c(2,3,1,1.5),oma = c(2,2,0.5,0.5)) layout(matrix(1:4,2,2)) plot(Q_treated~Date,hydroP_rslt,type="l",las=1,ann=F, ylim=c(0,0.4), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines( Qin~Date,hydroP_rslt,col="dodgerblue1") mtext(side=3,adj=0.01,"WY 1966") mtext(side=2,line=3,expression("Outflow Q (hm"^"3"~"d"^"-1"*")")) legend("topleft",inset = c(0.01,0.01),legend=c("Inflow","Outflow"), pch=c(NA),lty=c(1),lwd=c(1), col=c("dodgerblue1",'black'),pt.bg=NA, pt.cex=1.5,ncol=1,cex=0.75,bty="n",y.intersp=1.25,x.intersp=0.5,xpd=NA,xjust=0.5,yjust=1) plot(Z_end_cm~Date,hydroP_rslt,type="l",las=1,ann=F, ylim=c(min(hydro_rslt$Z_end_cm), max(hydro_rslt$Z_end_cm)), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) mtext(side=2,line=3,"Depth (cm)") mtext(side=1,line=2.5,"Date") plot(L_out~Date,hydroP_rslt,type="l",las=1,ann=F, ylim=c(0, max(hydroP_rslt$Lin)*1.05), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines(Lin~Date,hydroP_rslt,col="dodgerblue1") mtext(side=2,line=2.75,expression("TP Load (kg"~"d"^"-1"*")")) mtext(side=3,adj=0.99,"Cell 1") plot(C_out~Date,hydroP_rslt,type="l",las=1,ann=F, ylim=c(0, max(hydroP_rslt$Cin)*1.05), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines(Cin~Date,hydroP_rslt,col="dodgerblue1") mtext(side=2,line=2.75,expression("TP FWM ("*mu*"g"~"L"^"-1"*")")) mtext(side=1,line=2.5,"Date") par(oldpar) ## ----------------------------------------------------------------------------- # input parameters # 1) Base hydrology params (shared structure) hydro_base <- list( A_cell = 2.19, # km2 # depths in cm Zmin = 2, # cm Zinit = 40, # cm Zweir = 0, # cm Q_zmin = 38, # cm Zrelease = 0, # cm Bypass_elev = 0, # hydraulics Q_a = 1.0, Q_b = 4.0, Width = 1.55, # km Qomax = 0.0, # hm3/day; 0 disables max cap in this implementation Qimax = 0.0, # hm3/day; 0 disables inflow cap # seepage (rates in m/day per m head; elevations in cm) Seepout_Rate = 0.0, Seepout_Elev = 0.0, # cm Seepin_Rate = 0.0, Seepin_Elev = 0.0, # cm ShutdownET = TRUE, force_Q_out = FALSE, DutyCycle = 0.95, Cmax = 2000 ) # 2) Base P params (shared structure) P_base <- list( # STA module C1000 = 22, Cstar = 3, Ks_per_yr = 16.8, Z1 = 40, Z2 = 100, Z3 = 200, Chalf = 300, K2Coef1 = 0, SeasonalFactor = 0, # keep 0 for base parity # PSTA (NEWS transition) Ytrans = 0, Ysigma = 0, Czero = 0, C1000_2 = NULL, ks_2 = 0, zh_2 = 0, # RES depth penalty k_depth_penalty = 1, # atmos + seepage water quality C_rain = 10, # ppb (ug/L) DryDepo = 20, # mg/m2-yr seepage_c = 20, # ppb cap for seep outflow seepin_conc = 0, # ppb # initial P state C_init_ppb = 30, Y_init_mgm2 = 1000 ) # 3) Cell-specific params params_cell1 <- modifyList(hydro_base, modifyList(P_base, list( Qin_Frac = 0.22, Seepout_Rate = 0.00789, Ks_per_yr = 16.8, Y_init_mgm2 = 3387.67297548954 ))) params_cell2 <- modifyList(hydro_base, modifyList(P_base, list( Qin_Frac = 0, Seepout_Rate = 0.00155, Ks_per_yr = 52.5, Y_init_mgm2 = 768.480041186681 ))) # 4) Build cells cells <- list( dmsta_make_cell( label = "CELL1", params = params_cell1, ttankS = 3.0, DownCell = 2L, Qin_Frac = params_cell1$Qin_Frac, RecycleIndex = 1L # self; can omit if your validator maps NA/0 -> self ), dmsta_make_cell( label = "CELL2", params = params_cell2, ttankS = 3.0, DownCell = 0L, Qin_Frac = params_cell2$Qin_Frac, RecycleIndex = 2L # self ) ) cells <- dmsta_validate_cells(cells) # Run case hydroP_net_rslt <- dmsta_flowP_case( series = series, cells = cells, Nsteps = 4L, max_iter = 1L, return_cell_series = TRUE, keep_Q17 = TRUE ) hydroP_case_rslt <- hydroP_net_rslt$results$case hydroP_cells_rslt <- hydroP_net_rslt$results$cells ## ----eval = FALSE,include = FALSE--------------------------------------------- # str(hydroP_net_rslt) ## ----hydroP_case_plot, echo=FALSE, fig.align = "center",fig.cap = "Simulated outflow discharge, water level, TP load and TP concentration for cell 1 (top) and cell 2 (bottom).",fig.width = 7,fig.height = 5,fig.align="center"---- oldpar <- par(no.readonly = TRUE) par(mar=c(3,3,1,1.5),oma = c(2,2,0.5,0.5)) layout(matrix(1:8,2,4,byrow=T)) for(i in 1:2){ tmp <- hydroP_cells_rslt[[i]] tmp$Z_end_cm <-m_to_cm(tmp$Z_end) tmp$C_in_total <- with(tmp,ifelse(Q_in_total==0,NA,C_in_total)) plot(Q_out_treated~Date,tmp,type="l",las=1,ann=F, ylim=c(0,0.4), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines( Q_in_total~Date,tmp,col="dodgerblue1") if(i==1){mtext(side=3,adj=0.01,"WY 1966")} mtext(side=2,line=3,expression("Outflow Q (hm"^"3"~"d"^"-1"*")")) legend("topleft",inset = c(0.01,0.01),legend=c("Inflow","Outflow"), pch=c(NA),lty=c(1),lwd=c(1), col=c("dodgerblue1",'black'),pt.bg=NA, pt.cex=1.5,ncol=1,cex=0.75,bty="n",y.intersp=1.25,x.intersp=0.5,xpd=NA,xjust=0.5,yjust=1) if(i ==2){mtext(side=1,line=2.5,"Date")} plot(Z_end_cm~Date,tmp,type="l",las=1,ann=F, ylim=c(min(hydro_rslt$Z_end_cm), max(hydro_rslt$Z_end_cm)), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) mtext(side=2,line=3,"Depth (cm)") if(i ==2){mtext(side=1,line=2.5,"Date")} plot(L_out_treated~Date,tmp,type="l",las=1,ann=F, ylim=c(0,max(tmp$L_in_tota)*1.05), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines(L_in_total~Date,tmp,col="dodgerblue1") mtext(side=2,line=2.75,expression("TP Load (kg"~"d"^"-1"*")")) if(i ==2){mtext(side=1,line=2.5,"Date")} plot(C_out_treated~Date,tmp,type="l",las=1,ann=F, ylim=c(0,max(tmp$C_in_total,na.rm=T)*1.05), xlim = c(as.Date("1965-05-01"),as.Date("1966-05-01"))) lines(C_in_total~Date,tmp,col="dodgerblue1") mtext(side=2,line=2.75,expression("TP FWM ("*mu*"g"~"L"^"-1"*")")) if(i ==2){mtext(side=1,line=2.5,"Date")} mtext(side=3,adj=0.99,paste("Cell",i)) } par(oldpar) ## ----case_net1, eval = FALSE-------------------------------------------------- # # 1) Build cases # cases <- list( # STA1_DW = list( # series_base = STA1_DW_input, # data.frame with Date, Qi, Ci, Rain, Et, Zcontrol # cells = STA1_DW_cell # list of dmsta_make_cell(...) objects # ), # STA1W = list( # series_base = STA1W_input, # cells = STA1W_cell # ) # ) ## ----case_net2, eval = FALSE-------------------------------------------------- # # 2) Build/parse routes # net <- data.frame( # CaseName = c("STA1_DW", "STA1W"), # Bypass_to = c("", ""), # Release1_to = c("", ""), # Release2_to = c("", ""), # Outflow_to = c("STA1W", "1"), # Seepage_to = c("", ""), # stringsAsFactors = FALSE # ) # routes <- build_routes_from_net_table(net, outlet_count = 1L) # ## ----net_example-------------------------------------------------------------- EC_network <- data.frame(CaseName = c("FEB55A_N", "FEB_S5A", "FEBS5A_OUT", "STA1_DW", "STA1W", "STA1E", "FEB_34", "FEB34_OUT", "STA2B", "STA34"), Bypass_to = c("FEB_S5A", "STA1_DW", "STA1_DW", "STA1E", "1", "2", "STA34", "STA34", "3", "4"), Release1_to = c("FEB_S5A", "STA1_DW", NA, NA, NA, NA, "STA34", NA, NA, NA), Release2_to = c("5", "5", NA, NA, NA, NA, "STA2B", NA, NA, NA), Outflow_to = c("FEB_S5A", "FEBS5A_OUT", "STA2B", "STA1W", "1", "2", "FEB34_OUT", "STA2B", "3", "4"), Seepage_to = c(NA, NA, NA, NA, NA, NA, "FEB34_OUT", NA, NA, NA)) EC_routes <- build_routes_from_net_table(EC_network, outlet_count = 5) t(EC_network) ## ----case_net3, eval = FALSE-------------------------------------------------- # # 3) Run network # out <- run_network_of_cases( # cases = cases, # routes = routes, # Nsteps = 4L, # return_cell_series = TRUE # )