ODNR makes budgets for P reduction in potential wetlands using the following method:
Put the lat/long of the outlet of the wetland into USGS Stream Stats to get the drainage area of the wetland (convert to acres)
Based on the size of the wetland correct for the amount of the drainage the wetland interacts with
# Drainage Area from Stream Stats
DA = 48
# Treated drainage area (corrected for drainage size)
tDA = if(DA <= 99){DA*0.95} else if(DA >= 100 & DA <= 999){DA*0.75} else if(DA >= 1000 & DA <= 9999){DA*0.25} else if(DA >= 10000 & DA <= 99999){DA*0.1} else if(DA >= 100000 & DA <= 999999){DA*0.01} else if(DA > 100000000){DA*0.005}
# units = acres
print(tDA)
## [1] 45.6
# P Loading Rates for Ohio from USGS SPARROW:
# HUC-8 Drainage basin ID
HUC8 = c(4100001,4100008,4100009,4100010,4100011,4100003,4100005,5120101,5040006,5090202,4120200,4120101,4100007,4100006,5040001,5080001,5060001,4100014,4110001,4110002,4110003,4110004)
# Total phosphorus yield for each basin (kg/km^2/year)
TPyieldkgkm=c(75.2,131,129,83.9,140,114,106,293,206,216,0,123,152,114,134,201,209,149,179,249,245,109)
loading=data.frame(HUC8,TPyieldkgkm)
# correct to lbs/acre/year from kg/km^2/year:
corr=112.0850759
loading$TPyield = loading$TPyieldkgkm/corr
# select the HUC-8 basin that the wetland is a part of
loadSJREHUC = loading$TPyield[loading$HUC8 == 4100003]
# multiply by the treated drainage area from step 2
loadSJRE = loadSJREHUC * tDA
# units = lbs / year
print(loadSJRE)
## [1] 46.37906
# P Removal Rates:
low = loadSJRE*0.18
hybrid = loadSJRE*0.64
high = loadSJRE*0.78
# units = lbs/year
print(low)
## [1] 8.34823
print(hybrid)
## [1] 29.6826
print(high)
## [1] 36.17566
4b. For floodplain wetlands estimates are based on constants from Noe et al. 2019: 3.2 pounds P/ acre / year removal
Not that many, but big assumptions
1. Drainage area is accurate
(it’s not)
2. SPARROW gives an accurate measurement of P loading
3. Literature estimates of loading are accurate
* Honestly this
probably isn’t the worst assumption if only because their estimates are
so wide ranging.
4. Seasonal variation does not matter to whoever is
looking for this budget.
5. Annual variation does not play a major
role in phosphorus loading, or retention, or hydrologic
dynamics.
Calculate drainage area and treated drainage area using digital elevation models:
# Drainage Area calculated by Bishwodeep using Digital Elevation Model (DEM)
tDA = 103.32 # units = acres
# Tile drainage very rough estimates by Melanie Coulter
tDANE = 16.8
tDAN = 25.9
tDANW = 38.5
# convert tDA to m^2
tDAm2NE=tDANE*4046.86
tDAm2N=tDAN*4046.86
tDAm2NW=tDANW*4046.86
tDAm2=tDA*4046.86
# units = m^2
print(tDAm2)
## [1] 418121.6
Then we import data on precipitation from the sensor network and from the nearby weather station at Ft. Wayne (from WeatherUnderground) which we will use to fill in our rainfall data.
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/SJRE/Data")
raindata<-read.csv("SJREsensors.csv")
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata_wayne<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%Y-%m-%d %H:%M:%S"))
dayrain_wayne = mutate(raindata_wayne, datetime = as.POSIXct(datetime, format = "%m/%d/%Y")) # rainfall data is in inches
# Here we will aggregate our rainfall values (just for the sensor data, the WeatherUnderground data is already the daily sum) so we have a sum of rainfall for each day
DF2 <- transform(raindata, Date = as.Date(datetime))
dayrain=aggregate(precipaccum ~ Date, DF2, sum) # rainfall data is in mm
We will calculate the runoff that flows through the wetland during rain events based on some information about soil properties.
\[Q=\begin{cases} 0 \text{ for } P \le I_a\\ {(P-I_a)^2 \over P-I_a+S} \text{ for } P > I_a \end{cases}\]
where:
Q = runoff (inches)
P = rainfall (inches)
S =
potential maximum soil moisture retention after runoff begins
(inches)
I_a = initial abstraction (inches; this is the amount of
rainfall needed for runoff to occur)
The runoff curve (CN) is used to calculate S \[S = {1000 \over CN} -10\] We will assume that: \[ I_a = 0.2S\] as a standard correction. (See here https://onlinelibrary.wiley.com/doi/10.1111/j.1752-1688.2006.tb04481.x for discussion of whether that is a good correction but we’re going to power through with it)
CN will depend on the soil of the site. In this case conventional
practice would place it around 80, fairly poorly drained (we know that
the surrounding fields flood often). However we are going to us 75 to be
further in line with NE Ohio tile drained fields Wessel et al. (2022)
# Start out by calculating S (inches)
CN=75
S = (1000/CN)-10
# Calculate I (inches)
Ia = 0.2*S
# convert rain to inches/day
dayrain$precip_in = dayrain$precipaccum/25.4
# caculate runoff (inches) for each day in the dataset
dayrain$Q = ifelse(dayrain$precip_in<=Ia, 0, ((dayrain$precip_in-Ia)^2)/(dayrain$precip_in-Ia+S))
# convert back to meters
dayrain$Qm = dayrain$Q * 0.0254
# scale up to the entire drainage area
dayrain$Qmfull=tDAm2*dayrain$Qm
# convert discharge to liters
dayrain$QL = dayrain$Qmfull * 1000
# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$Date,format="%Y-%m-%d")
# Summarize the dataset by month and year
monthrain=dayrain %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(QL=sum(QL,na.rm=T))
Now we add in the monthly measurements of Phosphorus concentrations to calculate P loading. We’ll start with pool B.
# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/SJRE/Data")
Pdata<-read.csv("SJRE P Budget Full inletsandoutlets.csv")
# subset out event sampling
Pdata<- subset(Pdata, Sampling == "baseflow")
# These are the first 6 rows to give you an idea what the datatable looks like
head(Pdata)
## Date Year Month Day month Sampling X1_Inlet_N.TP X2_Inlet_NW.TP
## 1 8/13/2021 2021 Aug 13 8 baseflow 0.206 0.230
## 2 10/11/2021 2021 Oct 11 10 baseflow NA NA
## 3 3/8/2022 2022 Mar 8 3 baseflow 0.063 0.073
## 4 6/8/2022 2022 Jun 8 6 baseflow NA NA
## 5 7/20/2022 2022 July 20 7 baseflow 0.037 0.123
## 6 8/2/2022 2022 Aug 2 8 baseflow 0.131 0.418
## X3_Inlet_NE.TP X18_Inlet_ditchExt.TP.OR.upstream17 X5_inlet_poolB
## 1 0.352 NA 0.0680
## 2 NA 0.232 NA
## 3 0.054 NA 0.0320
## 4 NA NA 0.2745
## 5 0.089 NA 0.4010
## 6 0.105 NA 0.4120
## X5_Outlet_poolB X15_Outlet_swale X17_Outlet_ditch.TP X19_River X1_Inlet_N.DRP
## 1 0.193 NA NA NA 0.070349937
## 2 0.077 NA 0.292 NA <NA>
## 3 0.046 0.094 NA NA 0.028
## 4 0.192 NA NA NA <NA>
## 5 0.082 0.216 0.124 NA bdl
## 6 0.062 NA 0.046 NA 0.014013
## X2_Inlet_NW.DRP X3_Inlet_NE.DRP X18_Inlet_ditchExt.DRP.OR.upstream17
## 1 0.1107251 0.2021195 <NA>
## 2 NA NA 0.1558
## 3 0.0571000 0.0225000 <NA>
## 4 NA NA <NA>
## 5 0.0493000 0.0495000 <NA>
## 6 0.1011440 0.0196180 <NA>
## X5_inlet_poolBDRP X5_Outlet_poolBDRP X15_Outlet_swaleDRP X17_Outlet_ditch.DRP
## 1 0.0173631 0.01563599 <NA> <NA>
## 2 <NA> 0.0269 <NA> 0.2369
## 3 bdl bdl 0.0146 <NA>
## 4 bdl bdl <NA> <NA>
## 5 0.0955 0.0121 bdl bdl
## 6 bdl bdl <NA> 0.043377
## X19_River.DRP X1_Inlet_N.TN X2_Inlet_NW.TN X3_Inlet_NE.TN
## 1 <NA> NA NA NA
## 2 <NA> NA NA NA
## 3 <NA> NA NA NA
## 4 <NA> NA NA NA
## 5 <NA> NA NA NA
## 6 <NA> NA NA NA
## X18_Inlet_ditchExt.TN.OR.upstream17 X15_Outlet_swale_TN X17_Outlet_ditch.TN
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
# merge the datasets
total <- merge(monthrain,Pdata,by=c("Month","Year"),all=TRUE)
# averaging the various inflows
total$inflowTP <- rowMeans(total[c('X1_Inlet_N.TP', 'X2_Inlet_NW.TP','X3_Inlet_NE.TP','X18_Inlet_ditchExt.TP.OR.upstream17')], na.rm=TRUE)
# Now calculate the average outflow TP
total$outflowTP <- rowMeans(total[c('X15_Outlet_swale', 'X17_Outlet_ditch.TP')], na.rm=TRUE)
# Measure the mass of P entering the wetland (units = mg P/year)
total$Pin=total$inflowTP*total$QL
# calculate the total mass of P exiting the wetland (units = mg P/year)
total$Pout=total$outflowTP*total$QL
# Convert to lbs P/year
total$Pinlb=total$Pin/453592
total$Poutlb=total$Pout/453592
# calculate the amount of P removed
total$Premoved=total$Pinlb-total$Poutlb
Removal = sum (total$Premoved,na.rm=T)
# This is the removal of P over the course of the 3 months included in this data
print(Removal)
## [1] 6.001022
# We can make a rough annual estimate by just multiplying by 4:
AnnRemoval=Removal*4
print(AnnRemoval)
## [1] 24.00409
# Calculate percent removal
total$percentP=total$Premoved/total$Pinlb
percentRemoval = mean (total$percentP,na.rm=T)
print(percentRemoval)
## [1] 0.9588235
Now we can start to look at how to deal with the months with rainfall data, but no phosphorus data
# make a dataset for only 2023
total23 <- subset(total, Year == "2023")
# clean up any repeat months
total23 <- total23 %>% group_by(Month) %>% summarise(across(c(Year, QL, inflowTP, outflowTP,Pin, Pout, Pinlb, Poutlb,Premoved,percentP), mean))
print(total23)
## # A tibble: 12 × 11
## Month Year QL inflowTP outflowTP Pin Pout Pinlb Poutlb
## <ord> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Jan 2023 NA 0.257 0.05 NA NA NA NA
## 2 Feb 2023 NA 0.167 0.31 NA NA NA NA
## 3 Mar 2023 NA 0.0733 NaN NA NA NA NA
## 4 Apr 2023 NA 0.05 0.31 NA NA NA NA
## 5 May 2023 NA 0.453 0.3 NA NA NA NA
## 6 Jun 2023 1723. 0.85 0.02 1464. 34.5 0.00323 7.60e-5
## 7 Jul 2023 2429094. 1.19 0.07 2890622. 170037. 6.37 3.75e-1
## 8 Aug 2023 0 1.73 0.04 0 0 0 0
## 9 Sep 2023 268703. NaN NaN NaN NaN NaN NaN
## 10 Oct 2023 529448. NaN NaN NaN NaN NaN NaN
## 11 Nov 2023 0 NaN NaN NaN NaN NaN NaN
## 12 Dec 2023 0 NaN NaN NaN NaN NaN NaN
## # ℹ 2 more variables: Premoved <dbl>, percentP <dbl>
Now we have a dataset with each month of the year where some months have phosphorus data but no discharge data, and some months have discharge data but no phoshorus data. Lets do the easy part first and fill in missing discharge data with rainfall data from Ft. Wayne:
# Now we come back to our previously loaded Fort Wayne Dataset
# don't need to convert rain to inches/day, but lets name it so our old code still works
# calculate runoff (inches) for each day in the dataset
dayrain_wayne$Q = ifelse(dayrain_wayne$precipaccum_in<=Ia, 0, ((dayrain_wayne$precipaccum_in-Ia)^2)/(dayrain_wayne$precipaccum_in-Ia+S))
# convert back to meters
dayrain_wayne$Qm = dayrain_wayne$Q * 0.0254
# scale up to the entire drainage area
dayrain_wayne$Qmfull=tDAm2*dayrain_wayne$Qm
# convert discharge to liters
dayrain_wayne$QL = dayrain_wayne$Qmfull * 1000
# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain_wayne$Date=as.POSIXct(dayrain_wayne$datetime,format="%Y-%m-%d")
# Summarize the dataset by month and year
monthrainft=dayrain_wayne %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(QL=sum(QL,na.rm=T))
Now we just need to slot our values from Ft. Wayne into SJRE. The data from Ft. Wayne is a good approximation of rainfall, but you miss important events, for example in November Ft. Wayne did not get enough rainfall to have runoff, whereas SJRE had a MAJOR rainfall event with more than 12 inches of rain in a single day.
# first we need to cut out the rows from the Ft. Wayne dataset where we already have rainfall data from the sensor network:
# As a note to future me: I am not including the June data from Ft. Wayne because EVEN THOUGH we only have 3 days of rainfall from the sensor network for June, those 3 days had MORE DISCHARGE THAN WAS AT FT. WAYNE FOR THE WHOLE MONTH. So from that we can assume the rest of June was pretty dry and our estimate even with only 3 days is better than Ft. Wayne for June. It's not as if it's less than the 3 days worth.
monthrainft = subset(monthrainft, Month == "Jan" | Month == "Feb" | Month == "Mar" | Month == "Apr" | Month == "May")
# now we need to slot in the Ft. Wayne data into our main dataset
totalplus23 <- merge(total23,monthrainft,by=c("Month","Year","QL"),all=TRUE)
total23final <- totalplus23 %>% group_by(Month) %>% summarise(across(c(Year, QL, inflowTP, outflowTP,Pin, Pout, Pinlb, Poutlb,Premoved,percentP), mean, na.rm=TRUE))
Now we have discharge estimates for all months of 2023 (some sensor data, some data from the Ft. Wayne Weather Station). We have phosphorus concentrations when water was available to be collected during monthly sampling. Our next steps are to 1. calculate loads and retention based on discharge from the Ft. Wayne data, 2. Estimate load and retention based on high and low values for months with missing concentration data.
# first make all NaNs into NAs
total23final <- total23final %>% mutate_all(~ifelse(is.nan(.), NA, .))
# percent phosphorus retention for each month
total23final$percentP=(total23final$inflowTP-total23final$outflowTP)/total23final$inflowTP
# We have one date where we have input data, but not output data. We're going to assume that it's the average retention of the previous and the prior month:
# the previous month is -0.86% retention, and the next month is -5.2% retention. The inflow concentration is 0.07333333
mean(c(-0.86,-5.2))
## [1] -3.03
-3.03*0.073
## [1] -0.22119
# replace the month in question (March 2023)
total23final$outflowTP=ifelse(test = total23final$Month == "3" & total23final$Year == "2023", yes = 0.22, no = total23final$outflowTP)
# Measure the mass of P entering the wetland (units = mg P/year)
total23final$Pin=total23final$inflowTP*total23final$QL
# calculate the total23final mass of P exiting the wetland (units = mg P/year)
total23final$Pout=total23final$outflowTP*total23final$QL
# Convert to lbs P/year
total23final$Pinlb=total23final$Pin/453592
total23final$Poutlb=total23final$Pout/453592
# calculate the amount of P removed
total23final$Premoved=total23final$Pinlb-total23final$Poutlb
# calculate the % P removed
total23final$percentP=total23final$Premoved/total23final$Pinlb
# calculate the average percent removal
mean(total23final$percentP,na.rm=TRUE)
## [1] -0.9673529
# calculate a percent removal by the LOAD of the wetland (this highlights the importance of load in calculating percent removal)
sum(total23final$Premoved,na.rm=TRUE)/sum(total23final$Pinlb,na.rm=TRUE)
## [1] 0.4561208
We can estimate the remaining months by estimating a low and a high
value for P loading and P retention.
We have discharge data for all
months, so we don’t want to just estimate loading, when we can estimate
concentration instead:
# our low estimate is the lowest concentration we saw:
lowPcon = min(total23final$inflowTP,na.rm = TRUE)
# our high estimate is the highest concentration we saw:
highPcon= max(total23final$inflowTP,na.rm = TRUE)
# then we can estimate the low end of % retention
lowPper = min(total23final$percentP,na.rm = TRUE)
# and then the high end of % retention
highPper = max(total23final$percentP,na.rm = TRUE)
# now calculate P load and retention using these concentrations to estimate the high and low end of P removal by the wetland:
# first for the low estimate
total23final$Premovedlow=ifelse(test = is.na(total23final$Premoved) , yes = total23final$QL*lowPcon/453592*lowPper, no = total23final$Premoved)
# then for the high estimate
total23final$Premovedhigh=ifelse(test = is.na(total23final$Premoved) , yes = total23final$QL*highPcon/453592*highPper, no = total23final$Premoved)
# Now we can sum up our full months:
lowRemoval = sum (total23final$Premovedlow,na.rm=T)
highRemoval = sum (total23final$Premovedhigh,na.rm=T)
print(lowRemoval)
## [1] 4.897512
print(highRemoval)
## [1] 8.327536
Now we can add in the sub-surface drainage that enters the wetland through the 3 tile drains: First we calculate the rainfall across the three different fields that drain into the wetland.
# reload our rainfall data:
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/SJRE/Data")
raindata<-read.csv("SJREsensors.csv")
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata_wayne<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%Y-%m-%d %H:%M:%S"))
dayrain_wayne = mutate(raindata_wayne, datetime = as.POSIXct(datetime, format = "%m/%d/%Y")) # rainfall data is in inches
# Here we will aggregate our rainfall values (just for the sensor data, the WeatherUnderground data is already the daily sum) so we have a sum of rainfall for each day
DF2 <- transform(raindata, Date = as.Date(datetime))
dayrain=aggregate(precipaccum ~ Date, DF2, sum) # rainfall data is in mm
# convert sensor data to meters
dayrain$precipaccum_m = dayrain$precipaccum/1000
dayrain_wayne$precipaccum_m = dayrain_wayne$precipaccum_in * 0.0254
# scale up to the drainage area for each of the three tiles in both of our datasets (this gets us meters cubed of rainfall)
dayrain_wayne$fullprecipN=tDAm2N*dayrain_wayne$precipaccum_m
dayrain_wayne$fullprecipNW=tDAm2NW*dayrain_wayne$precipaccum_m
dayrain_wayne$fullprecipNE=tDAm2NE*dayrain_wayne$precipaccum_m
dayrain$fullprecipN=tDAm2N*dayrain$precipaccum_m
dayrain$fullprecipNW=tDAm2NW*dayrain$precipaccum_m
dayrain$fullprecipNE=tDAm2NE*dayrain$precipaccum_m
# convert rainfall to liters
dayrain$precipLN = dayrain$fullprecipN * 1000
dayrain$precipLNW = dayrain$fullprecipNW * 1000
dayrain$precipLNE = dayrain$fullprecipNE * 1000
dayrain_wayne$precipLN = dayrain_wayne$fullprecipN * 1000
dayrain_wayne$precipLNW = dayrain_wayne$fullprecipNW * 1000
dayrain_wayne$precipLNE = dayrain_wayne$fullprecipNE * 1000
# Summarize the dataset by month and year
library(lubridate)
library(tidyverse)
dayrain_wayne$Date=as.POSIXct(dayrain_wayne$datetime,format="%Y-%m-%d")
dayrain$Date=as.POSIXct(dayrain$Date,format="%Y-%m-%d")
# Merge everything together. There's probably a more elegant way to code this, but I don't know how to do it and merging in R is super complicated.
# Ft Wayne data first
monthrainftN=dayrain_wayne %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(precipLN=sum(precipLN,na.rm=T))
monthrainftNW=dayrain_wayne %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(precipLNW=sum(precipLN,na.rm=T))
monthrainftNE=dayrain_wayne %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(precipLNE=sum(precipLN,na.rm=T))
monthrainft1 = merge(monthrainftN,monthrainftNW,by=c("Month","Year"),all=TRUE)
monthrainft = merge(monthrainftNE,monthrainft1,by=c("Month","Year"),all=TRUE)
# Now sensor data
monthrainN=dayrain %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(precipLN=sum(precipLN,na.rm=T))
monthrainNW=dayrain %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(precipLNW=sum(precipLN,na.rm=T))
monthrainNE=dayrain %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
group_by(Year,Month) %>% summarise(precipLNE=sum(precipLN,na.rm=T))
monthrain1 = merge(monthrainN,monthrainNW,by=c("Month","Year"),all=TRUE)
monthrain = merge(monthrainNE,monthrain1,by=c("Month","Year"),all=TRUE)
##
# Now we can cut out the months that we already have data for in the sensor data from the Ft. Wayne dataset
monthrainft = subset(monthrainft, Month == "Jan" | Month == "Feb" | Month == "Mar" | Month == "Apr" | Month == "May")
# now we slot in the Ft. Wayne data into our main dataset and have precipitation for the whole of 2023 for each of the 3 tile drains
monthrainfull <- merge(monthrain,monthrainft,by=c("Month","Year","precipLN","precipLNE","precipLNW"),all=TRUE)
# We expect that subsurface drainage will be ~ 31% of total runoff (Pease et al. 2018)
monthrainfull$QLN = monthrainfull$precipLN * 0.31
monthrainfull$QLNE = monthrainfull$precipLNE * 0.31
monthrainfull$QLNW = monthrainfull$precipLNW * 0.31
Now that we have rainfall and sub-surface discharge for each tile drain in 2023 we can add in phosphorus data to predict the full retention at the site.
# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/SJRE/Data")
Pdata<-read.csv("SJRE P Budget Full inletsandoutlets.csv")
# subset out event sampling
Pdata<- subset(Pdata, Sampling == "baseflow")
# subset out 2023 data
Pdata2023<- subset(Pdata, Year == 2023)
# combine rows where there are two measurements in a single month
Pdata2023 = Pdata2023 %>% group_by(Month) %>% summarise(across(c(X1_Inlet_N.TP, X2_Inlet_NW.TP,X3_Inlet_NE.TP,X18_Inlet_ditchExt.TP.OR.upstream17,X15_Outlet_swale, X17_Outlet_ditch.TP), mean))
# merge the datasets
total23 <- merge(monthrainfull,Pdata2023,by=c("Month"),all=TRUE)
# calculate the phosphorus load entering the system for each tile drain each month:
total23$PinN = total23$QLN * total23$X1_Inlet_N.TP
total23$PinNW = total23$QLNW * total23$X2_Inlet_NW.TP
total23$PinNE = total23$QLNE * total23$X3_Inlet_NE.TP
# now we can fill in the gaps with a high and low for each individual tile drain:
lowPconN = min(total23$X1_Inlet_N.TP,na.rm = TRUE)
lowPconNW = min(total23$X2_Inlet_NW.TP,na.rm = TRUE)
lowPconNE = min(total23$X3_Inlet_NE.TP,na.rm = TRUE)
highPconN= max(total23$X1_Inlet_N.TP,na.rm = TRUE)
highPconNW= max(total23$X2_Inlet_NW.TP,na.rm = TRUE)
highPconNE= max(total23$X3_Inlet_NE.TP,na.rm = TRUE)
# now we will fill in those months with a high and low estimate of P load and then convert to lbs:
total23$PinNlow = ifelse(test = is.na(total23$PinN) , yes = total23$QLN * lowPconN / 453592, no = total23$PinN / 453592)
total23$PinNWlow = ifelse(test = is.na(total23$PinNW) , yes = total23$QLNW * lowPconNW / 453592, no = total23$PinNW / 453592)
total23$PinNElow = ifelse(test = is.na(total23$PinNE) , yes = total23$QLNE * lowPconNE / 453592, no = total23$PinNE / 453592)
total23$PinNhigh = ifelse(test = is.na(total23$PinN) , yes = total23$QLN * highPconN / 453592, no = total23$PinN / 453592)
total23$PinNWhigh = ifelse(test = is.na(total23$PinNW) , yes = total23$QLNW * highPconNW / 453592, no = total23$PinNW / 453592)
total23$PinNEhigh = ifelse(test = is.na(total23$PinNE) , yes = total23$QLNE * highPconNE / 453592, no = total23$PinNE / 453592)
# now we combine all those for a total load of P into the system
total23$Pinlbslow = total23$PinNlow + total23$PinNElow + total23$PinNWlow
total23$Pinlbshigh = total23$PinNhigh + total23$PinNEhigh + total23$PinNWhigh
# calculate the average outflow TP
total23$outflowTP <- rowMeans(total23[c('X15_Outlet_swale', 'X17_Outlet_ditch.TP')], na.rm=TRUE)
# calculate the load of P leaving the system (assuming that all water in goes out)
total23$QLtotal = total23$QLN + total23$QLNE + total23$QLNW
total23$Poutlbs = total23$QLtotal * total23$outflowTP / 459592
# first make all NaNs into NAs
total23 <- total23 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# calculate % retention for each month (I'm using PinlbsHigh because it doesn't matter, we're only calculating for months with outflow data and we only have a range of values when there is no data collected)
total23$percentP <- (total23$Pinlbshigh - total23$Poutlbs) / total23$Pinlbshigh
# now we will apply a high and low retention value to each month where we're missing data in order to predict the P retention in those months:
# estimate the low end of % retention
lowPper = min(total23$percentP,na.rm = TRUE)
# and then the high end of % retention
highPper = max(total23$percentP,na.rm = TRUE)
lowPper
## [1] -5.119059
highPper
## [1] 0.969158
# Now removal of P by the wetland
total23$Premovedhighsub=ifelse(test = is.na(total23$Poutlbs) , yes = (highPper/100)*total23$Pinlbshigh, no = total23$Pinlbshigh - total23$Poutlbs)
total23$Premovedlowsub=ifelse(test = is.na(total23$Poutlbs) , yes = (lowPper/100)*total23$Pinlbslow, no = total23$Pinlbslow - total23$Poutlbs)
Now we want to compare our surface flow and sub-surface flow contributions to P retention:
# We pare down our dataframes so they're only the high and low estimates of retention (in lbs)
total23subsurface = total23[c("Month","Year","Premovedhighsub","Premovedlowsub")]
total23surface = total23final[c("Month","Year","Premovedhigh","Premovedlow")]
totalplot = merge(total23subsurface,total23surface,by=c("Month","Year"),all=TRUE)
# Make a combined variable:
totalplot$Premovedhightotal = totalplot$Premovedhighsub + totalplot$Premovedhigh
totalplot$Premovedlowtotal = totalplot$Premovedlowsub + totalplot$Premovedlow
# Now we can look at the high and low estimates of retention by the site:
sum(totalplot$Premovedhightotal)
## [1] 40.49657
sum(totalplot$Premovedlowtotal)
## [1] 19.66171
Now lets plot it out real quick to see the annual variation:
# lets put our dataframe into the tall format so we can plot it more easily:
dfplots = totalplot %>%
gather(highlow,Pretained, Premovedlowtotal, Premovedhightotal)
# now make the plot
ggplot(data=dfplots, aes(x=Month, y=Pretained, group=highlow)) +
geom_line(aes(linetype=highlow),size=1)+
geom_point(size=2)+
scale_x_continuous(name="Month", breaks = 1:12, limits=c(0, 12))+
theme_classic()