Phosphorus Budgets for SJRE

H0: ODNR Budgets:

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

0-99 acres: 95%
100-999 acres: 75%
1,000-9,999 acres: 25%
10,000-99,999 acres: 10%
100,000-999,999 acres: 1%
1,000,000 acres +: 0.5%

# 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

Estimate P loading rate using the USGS SPARROW model (units are pounds P / year) based on total drainage area treated (calculated in step 2)

# 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

Estimate P removal based on literature (as proposed by TetraTech; units: pounds P/ acre / year):

Low estimate: 18% removal of P loading
Hybrid estimate: 64% removal of P loading
High estimate: 78% of P loading

# 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

ASSUMPTIONS IN ODNR BUDGETS:

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.

H1: Main flow path budget calculating Runoff from Rainfall

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()

ASSUMPTIONS IN BUDGET:

Drainage area is accurate
All inflowing water leaves the wetland as surfacewater outflow
Groundwater effects on water and phosphorus budgets are negligible
Baseflow P reduction = storm flow P reduction
All water entering the wetland leaves the wetland
The average of inflow concentrations represents average P concentration entering the wetland (This likely is NOT true because the amounts of water flowing through each input is different)
Pool B represents the whole wetland

Phosphorus Budgets for SJRE

Kenneth Anderson

2023-10-20

H0: ODNR Budgets:

ASSUMPTIONS IN ODNR BUDGETS:

H1: Main flow path budget calculating Runoff from Rainfall

ASSUMPTIONS IN BUDGET: