--- title: "Getting Started with DMSTAr" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting Started with DMSTAr} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Introduction **DMSTAr** is an R implementation of core concepts from the Dynamic Model for Stormwater Treatment Areas (DMSTA), focused on daily phosphorus transport and treatment. The package is designed to: - simulate single treatment cells or networks of cells ("cases") - reproduce DMSTA-style daily mass balance behavior - support transparent, reproducible workflows in R This vignette provides a practical introduction to the main ideas and a minimal working example. ## Core Concepts Before running a model, it hepls to understand a few core terms/concepts used throghout `DMSTAr` taken from DMSTA. For more information on DMSTA specifically, please visit the [DMSTA2](http://wwwalker.net/dmsta/) webpage. A **case** represents a single modeled unit (e.g., an STA cell, reservoir, or reach). Each case has: - inflows and outflows - a water volume - phosphorus mass and concentration Cases can be connected into a **network**, where outflow from one case becomes inflow to another. Routing is defined explicitly using fractions and (optionally) lags. DMSTAr operates on a daily timestep and conserves mass by tracking: - water volume - phosphorus load - concentration derived from mass and volume ## Minimal Examples Below are a couple of examples of how use the functions in this R-package to do simple simulations. The three main modeling functions in `DMSTAr` is `dmsta_flow_series(...)` (single cell hydrology only), `dmsta_flowP_series(...)` (single cell hydrology and P) and `dmsta_flowP_case(...)` (networked hydrology and P). Each function requires input parameters and input (i.e. forcing) data to simulate outflow conditions. For input data purposes we will use the internal data series in this package: ```{r} library(DMSTAr) data(series) ``` To see what this is use `?DMSTAr::series` ### Hydrology only For purposes of this demonstration we will use only the first year and half of the example input data frame. ```{r} # 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) ``` ```{r,eval = FALSE} str(out_hydro) ``` Each function provides a list of data.frames including results (`out$results`), daily water budget (`out$budgets$water`) and metadata (`out$meta`). ```{r 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) ``` ### Hydrology and P ```{r} ## 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) ``` ```{r,eval = FALSE} str(hydroP_out) ``` Similarly to the hydrology outputs `out` contains a list of data frames, the only difference in this function also include a P mass balance budget, to view see `out$budgets$mass`. ```{r 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) ``` ### Hydrology and P Case `dmsta_flowP_case` allows for the simulation of cells in series where for example, Cell 1 flows into Cell 2. It is part wrapper function but also links more complex interactions of seepage and recycled discharge/load between cells. ```{r} # 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 ``` ```{r,eval = FALSE,include = FALSE} str(hydroP_net_rslt) ``` Much like the other functions this function stores results, water and mass budgets and meta data. The difference being this function stores the data at the *"case"* level (inflow to cell 1 and outflow of cell 2) as well as cell specific information. ```{r 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 Network Simulation Here is example code on how to perform a case network simulation. Due to the involved nature of this process, this code is here just as an example but given processes above it can be easily replicated given the necessary data. ```{r 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 ) ) ``` First, you need to build the cases. It is simply a list with case name, input data and cell specific information. ```{r 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) ``` Next is defining the case network where you can specify the where the different outflow can be routed to. For instance, above networks `STA1_DW` outflow to `STA1W` then `STA1W` is discharged to `OUTLET 1`. This can be as complicated as you'd like. For instance here is the Eastern and Central flow-path network simulation for the Everglades STAs as part of the Restoration Strategies modeling efforts. ```{r 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) ``` ```{r case_net3, eval = FALSE} # 3) Run network out <- run_network_of_cases( cases = cases, routes = routes, Nsteps = 4L, return_cell_series = TRUE ) ``` In the third step, the actual simulation occurs similar to the `dmsta_flowP_case` example above but for each case within the estimated network. Outputs are similar but also includes some additional outputs.