Phosphorus Budgets for FORB

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 = 63
# 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] 59.85

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 
loadFORBHUC = loading$TPyield[loading$HUC8 == 4100005]
# multiply by the treated drainage area from step 2
loadFORB = loadFORBHUC * tDA
# units = lbs / year
print(loadFORB)

## [1] 56.60076

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 = loadFORB*0.18
hybrid = loadFORB*0.64
high = loadFORB*0.78

# units = lbs/year
print(low)

## [1] 10.18814

print(hybrid)

## [1] 36.22448

print(high)

## [1] 44.14859

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 boy are they big assumptions
1. Drainage area is accurate (it’s not)
2. SPARROW gives an accurate measurement of P loading for wetlands
3. Literature estimates of retention 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: Rainfall at Ft. Wayne Budget, Drainage area from high resolution DEM

We’ll use data from Ft. Wayne Indiana (31 miles away) downloaded from Weather Underground.

We’ll calculating runoff using the curve number of the soils: 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 it’s probably around 80, fairly poorly drained, but we may want to come back and think abit more about what our curve number should be.

The drainage area we are using here accounts mainly for the drainage entering through the culvert which we expect the majority of to run through the wetland and out the stream into the Maumee river. The majority of this area is the field to the southeast of the FORB land parcel.

# Drainage Area calculated by DEM (164.41 acres, 66.54 hectares)
tDA = 164.41
# convert tDA to m^2
tDAm2=tDA*4046.86
# units = m^2
print(tDAm2)

## [1] 665344.3

Then we import data on precipitation from weather underground at the Ft. Wayne weather station (31 miles away)

setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))
# Rename so our old code keeps working 
dayrain=raindata

Then we calculate runoff as above.

# Start out by calculating S (inches)
CN=85
S = (1000/CN)-10
# Calculate I (inches)
Ia = 0.2*S
# don't need to convert rain to inches/day, but lets name it so our old code still works
dayrain$precip = dayrain$precipaccum_in
# calculate runoff (inches) for each day in the dataset
dayrain$Q = ifelse(dayrain$precip<=Ia, 0, ((dayrain$precip-Ia)^2)/(dayrain$precip-Ia+S))
# scale up to the entire drainage area
dayrain$Q=tDAm2*dayrain$Q
dayrain$fullprecip=tDAm2*dayrain$precipaccum_in
# convert back to meters 
dayrain$Qm = dayrain$Q * 0.0254
dayrain$precipm = dayrain$fullprecip * 0.0254
# convert discharge to liters
dayrain$QL = dayrain$Qm * 1000
dayrain$precipL = dayrain$precipm * 1000
# Sum up annual precipitation and discharge
sum(dayrain$precipL)

## [1] 570028366

sum(dayrain$QL)

## [1] 71443877

# Sum up data to get monthly values
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.2

## Warning: package 'forcats' was built under R version 4.3.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4     ✔ tibble  3.2.1
## ✔ purrr   1.0.2     ✔ tidyr   1.3.0
## ✔ readr   2.1.4

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

dayrain$Date=as.POSIXct(dayrain$datetime,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))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

Now we add in the monthly measurements of Phosphorus concentrations to calculate P loading.

# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB P Budget Full inletsandoutlets.csv")
# subset out event sampling 
Pdata<- subset(Pdata, Sampling == "baseflow")
# take the average outflow between streams 1 and 2 (for samples at the end of the stream)
Pdata$outflowTP <- rowMeans(Pdata[c('Outlet_stream1_FORB_17_END', 'Outlet_stream2_FORB_18')], na.rm=TRUE)
# Average the inflow between the main inlet tile and the start of complex 4 (the main inlet tile is not sampled as often as the start of complex 4)
Pdata$inflowTP <- rowMeans(Pdata[c('Inlet_wetland_FORB_9', 'Inlet_tile_FORB_15')], na.rm=TRUE)
# merge the datasets
total <- merge(Pdata,monthrain,by=c("Month","Year"),all=TRUE)
# Measure the mass of P entering the wetland (units = mg P/month)
total$Pin=total$inflowTP*total$QL
# calculate the total mass of P exiting the wetland (units = mg P/month)
total$Pout=total$outflowTP*total$QL
# Convert to lbs P/month
total$Pinlb=total$Pin/453592
total$Poutlb=total$Pout/453592
# calculate the amount of P removed
total$Premoved=total$Pinlb-total$Poutlb

Now we need to look at where we have data and attempt to extrapolate where we do not

# To start lets look at just 2023
total23<- subset(total, Year == "2023")
total23<- subset(total23, select=c(Month,Year,Pinlb,Poutlb,Premoved,QL,inflowTP,outflowTP))
print(total23)

##    Month Year      Pinlb    Poutlb   Premoved         QL inflowTP outflowTP
## 1    Apr 2023  0.8929178 0.1488196  0.7440982  1350068.0   0.3000     0.050
## 3    Aug 2023         NA        NA         NA  5348405.2       NA        NA
## 4    Dec 2023         NA        NA         NA   103484.1       NA        NA
## 5    Feb 2023  9.0789473 4.0714101  5.0075372 13479993.0   0.3055     0.137
## 6    Jan 2023  0.7081859 0.4902826  0.2179034  1235490.3   0.2600     0.180
## 8    Jul 2023 19.7221139       NaN        NaN 14860121.4   0.6020       NaN
## 10   Jun 2023        NaN       NaN        NaN   635357.9      NaN       NaN
## 12   Mar 2023  6.7371064 1.9248875  4.8122188 14551893.1   0.2100     0.060
## 13   May 2023 17.0749616 5.3359255 11.7390361 16135554.1   0.4800     0.150
## 15   Nov 2023        NaN       NaN        NaN        0.0      NaN       NaN
## 18   Oct 2023         NA        NA         NA  1695046.7       NA        NA
## 20   Sep 2023        NaN       NaN        NaN  2048463.7      NaN       NaN

We can assume there is no P retained when there is no discharge:

library(dplyr)
# first make all NaNs into NAs
total23 <- total23 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# now replace NAs where discharge = 0 as also 0
total23$Premoved=ifelse(test = total23$QL == 0, yes = "0", no = total23$Premoved)

We have one event with loading data but no retention data. Here we can estimate retention based on the average retention of the other months.

# tell r that Premoved is numeric: 
total23$Premoved=as.numeric(total23$Premoved)
# Calculate percent removal
total23$percentP=total23$Premoved/total23$Pinlb
percentRemoval = mean (total23$percentP,na.rm=T)
print(percentRemoval)

## [1] 0.6188732

There’s a pretty good amount of variation here, but lets power through and use the average.

# replace the month in question (July 2023)
total23$Premoved=ifelse(test = total23$Month == "Jul" & total23$Year == "2023", yes = total23$Pinlb*percentRemoval, no = total23$Premoved)

Now we can sum up the total P removed by the wetland over the course of the year:

Removal = sum (total23$Premoved,na.rm=T)
# This is the removal of P by pool by over the course of the 8 months that we have data for
print(Removal)

## [1] 34.72628

We can estimate the other 4 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(total23$inflowTP,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon= max(total23$inflowTP,na.rm = TRUE)

# then we can 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)


# 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
total23$Premovedlow=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*lowPcon/453592*lowPper, no = total23$Premoved)
# then for the high estimate
total23$Premovedhigh=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*highPcon/453592*highPper, no = total23$Premoved)

# Now we can sum up our full months: 
lowRemoval = sum (total23$Premovedlow,na.rm=T)
highRemoval = sum (total23$Premovedhigh,na.rm=T)
print(lowRemoval)

## [1] 36.1267

print(highRemoval)

## [1] 45.59897

ASSUMPTIONS IN THIS BUDGET:

Drainage area is accurate (it accounts for only topographic drainage area)
Curve number is a good choice for this site on average (we chose a higher number to better account for sub-surface drainage)
All inflowing water leaves the wetland as surfacewater outflow
Groundwater effects on water and phosphorus budgets are negligible
Tile drains are accounted for (they’re not)
Nearby rainfall data is accurate for the site
Storm events do not change phosphorus retention dynamics

H2: Rainfall from Ft. Wayne, expert estimate of drainage area

# Drainage Area is an expert estimate by the District Administrator of the Paulding County Soil and Water Conservation District which accounts for both surface and sub-surface drainage
tDA = 16.5
# convert tDA to m^2
tDAm2=tDA*4046.86
# units = m^2
print(tDAm2)

## [1] 66773.19

Then we import data on precipitation from weather underground at the Ft. Wayne weather station (31 miles away)

setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))
# Rename so our old code keeps working 
dayrain=raindata

Then we calculate runoff as above.

# Start out by calculating S (inches)
CN=85
S = (1000/CN)-10
# Calculate I (inches)
Ia = 0.2*S
# don't need to convert rain to inches/day, but lets name it so our old code still works
dayrain$precip = dayrain$precipaccum_in
# caculate runoff (inches) for each day in the dataset
dayrain$Q = ifelse(dayrain$precip<=Ia, 0, ((dayrain$precip-Ia)^2)/(dayrain$precip-Ia+S))
# scale up to the entire drainage area
dayrain$Q=tDAm2*dayrain$Q
dayrain$fullprecip=tDAm2*dayrain$precipaccum_in
# convert back to meters 
dayrain$Qm = dayrain$Q * 0.0254
dayrain$precipm = dayrain$fullprecip * 0.0254
# convert discharge to liters
dayrain$QL = dayrain$Qm * 1000
dayrain$precipL = dayrain$precipm * 1000
# Sum up annual precipitation and discharge
sum(dayrain$precipL)

## [1] 57207396

sum(dayrain$QL)

## [1] 7170026

# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$datetime,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))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

Now we add in the monthly measurements of Phosphorus concentrations to calculate P loading as above.

# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB P Budget Full inletsandoutlets.csv")
# subset out event sampling 
Pdata<- subset(Pdata, Sampling == "baseflow")
# take the average outflow between streams 1 and 2 (for samples at the end of the stream)
Pdata$outflowTP <- rowMeans(Pdata[c('Outlet_stream1_FORB_17_END', 'Outlet_stream2_FORB_18')], na.rm=TRUE)
# Average the inflow between the main inlet tile and the start of complex 4 (the main inlet tile is not sampled as often as the start of complex 4)
Pdata$inflowTP <- rowMeans(Pdata[c('Inlet_wetland_FORB_9', 'Inlet_tile_FORB_15')], na.rm=TRUE)
# merge the datasets
total <- merge(Pdata,monthrain,by=c("Month","Year"),all=TRUE)
# Measure the mass of P entering the wetland (units = mg P/month)
total$Pin=total$inflowTP*total$QL
# calculate the total mass of P exiting the wetland (units = mg P/month)
total$Pout=total$outflowTP*total$QL
# Convert to lbs P/month
total$Pinlb=total$Pin/453592
total$Poutlb=total$Pout/453592
# calculate the amount of P removed
total$Premoved=total$Pinlb-total$Poutlb

# Extrapolate where we don't have data
total23<- subset(total, Year == "2023")
total23<- subset(total23, select=c(Month,Year,Pinlb,Poutlb,Premoved,QL,inflowTP,outflowTP))
print(total23)

##    Month Year      Pinlb     Poutlb   Premoved         QL inflowTP outflowTP
## 1    Apr 2023 0.08961222 0.01493537 0.07467685  135491.28   0.3000     0.050
## 3    Aug 2023         NA         NA         NA  536759.84       NA        NA
## 4    Dec 2023         NA         NA         NA   10385.55       NA        NA
## 5    Feb 2023 0.91115279 0.40860207 0.50255072 1352836.72   0.3055     0.137
## 6    Jan 2023 0.07107273 0.04920420 0.02186853  123992.39   0.2600     0.180
## 8    Jul 2023 1.97928885        NaN        NaN 1491344.83   0.6020       NaN
## 10   Jun 2023        NaN        NaN        NaN   63763.79      NaN       NaN
## 12   Mar 2023 0.67612831 0.19317952 0.48294879 1460411.39   0.2100     0.060
## 13   May 2023 1.71362366 0.53550739 1.17811626 1619345.80   0.4800     0.150
## 15   Nov 2023        NaN        NaN        NaN       0.00      NaN       NaN
## 18   Oct 2023         NA         NA         NA  170112.95       NA        NA
## 20   Sep 2023        NaN        NaN        NaN  205581.48      NaN       NaN

library(dplyr)
# first make all NaNs into NAs
total23 <- total23 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# now replace NAs where discharge = 0 as also 0
total23$Premoved=ifelse(test = total23$QL == 0, yes = "0", no = total23$Premoved)
# We have one event with loading data but no retention data. Here we can estimate retention based on the average retention of the other months. 
# tell r that Premoved is numeric: 
total23$Premoved=as.numeric(total23$Premoved)
# Calculate percent removal
total23$percentP=total23$Premoved/total23$Pinlb
percentRemoval = mean (total23$percentP,na.rm=T)
print(percentRemoval)

## [1] 0.6188732

# replace the month in question (July 2023)
total23$Premoved=ifelse(test = total23$Month == "Jul" & total23$Year == "2023", yes = total23$Pinlb*percentRemoval, no = total23$Premoved)

Removal = sum (total23$Premoved,na.rm=T)
# This is the removal of P by pool by over the course of the 8 months that we have data for
print(Removal)

## [1] 3.48509

# our low estimate is the lowest concentration we saw: 
lowPcon = min(total23$inflowTP,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon= max(total23$inflowTP,na.rm = TRUE)

# then we can 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)


# 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
total23$Premovedlow=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*lowPcon/453592*lowPper, no = total23$Premoved)
# then for the high estimate
total23$Premovedhigh=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*highPcon/453592*highPper, no = total23$Premoved)

# Now we can sum up our full months: 
lowRemoval = sum (total23$Premovedlow,na.rm=T)
highRemoval = sum (total23$Premovedhigh,na.rm=T)
print(lowRemoval)

## [1] 3.625634

print(highRemoval)

## [1] 4.57626

ASSUMPTIONS IN THIS BUDGET:

Drainage area is accurate (it may not account for the full topographic drainage area)
Curve number is a good choice for this site on average (we chose a higher number to better account for sub-surface drainage)
All inflowing water leaves the wetland as surfacewater outflow
Groundwater effects on water and phosphorus budgets are negligible
Tile drains are accounted for (they’re not)
Nearby rainfall data is accurate for the site
Storm events do not change phosphorus retention dynamics

H3: Rainfall from Ft. Wayne, expert estimate of drainage area, estimates of runoff from Pease et al. 2018

# Drainage Area delineated by expert opinion
tDA = 16.3
tDAsub = 7.3
# convert tDA to m^2
tDAm2=tDA*4046.86
tDAm2s=tDAsub*4046.86
# units = m^2

Then we import data on precipitation from weather underground at the Ft. Wayne weather station (31 miles away)

setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))
# Rename so our old code keeps working 
dayrain=raindata

Then we estimate runoff according to Pease et al. 2018 with ~7% of rainfall as surface runoff, and 31% of rainfall as subsurface runoff for a total of 38% of rainfall as runoff but it gets a bit tricky where there are different drainage areas for surface vs. subsurface water.

# don't need to convert rain to inches/day, but lets name it so our old code still works
dayrain$precip = dayrain$precipaccum_in
# scale up to the entire drainage area
dayrain$fullprecip=tDAm2*dayrain$precipaccum_in
dayrain$fullprecipsub=tDAm2s*dayrain$precipaccum_in
totalprecipL = sum(dayrain$fullprecip)*0.0254*1000
# We expect that subsurface drainage will be ~ 31% of total rainfall and surface runoff will be ~ 7% of total rainfall (Pease et al. 2018)
dayrain$Qsur = dayrain$fullprecip * 0.07
dayrain$Qsub = dayrain$fullprecipsub * 0.31
# now add those up:
dayrain$Q = dayrain$Qsur+dayrain$Qsub
# convert back to meters 
dayrain$Qm = dayrain$Q * 0.0254
# convert discharge to liters
dayrain$QL = dayrain$Qm * 1000
# Sum up annual precipitation and discharge
sum(dayrain$pQL)

## [1] 0

# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$datetime,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))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

sum(monthrain$QL)

## [1] 11802059

Now we add in the monthly measurements of Phosphorus concentrations to calculate P loading.

# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB P Budget Full inletsandoutlets.csv")
# subset out event sampling 
Pdata<- subset(Pdata, Sampling == "baseflow")
# take the average outflow between streams 1 and 2 (for samples at the end of the stream)
Pdata$outflowTP <- rowMeans(Pdata[c('Outlet_stream1_FORB_17_END', 'Outlet_stream2_FORB_18')], na.rm=TRUE)
# Average the inflow between the main inlet tile and the start of complex 4 (the main inlet tile is not sampled as often as the start of complex 4)
Pdata$inflowTP <- rowMeans(Pdata[c('Inlet_wetland_FORB_9', 'Inlet_tile_FORB_15')], na.rm=TRUE)
# merge the datasets
total <- merge(Pdata,monthrain,by=c("Month","Year"),all=TRUE)
# Measure the mass of P entering the wetland (units = mg P/month)
total$Pin=total$inflowTP*total$QL
# calculate the total mass of P exiting the wetland (units = mg P/month)
total$Pout=total$outflowTP*total$QL
# Convert to lbs P/month
total$Pinlb=total$Pin/453592
total$Poutlb=total$Pout/453592
# calculate the amount of P removed
total$Premoved=total$Pinlb-total$Poutlb

# Extrapolate where we don't have data
total23<- subset(total, Year == "2023")
total23<- subset(total23, select=c(Month,Year,Pinlb,Poutlb,Premoved,QL,inflowTP,outflowTP))
print(total23)

##    Month Year     Pinlb     Poutlb  Premoved        QL inflowTP outflowTP
## 1    Apr 2023 0.4859780 0.08099633 0.4049817  734785.8   0.3000     0.050
## 3    Aug 2023        NA         NA        NA  867747.0       NA        NA
## 4    Dec 2023        NA         NA        NA  573832.7       NA        NA
## 5    Feb 2023 0.9803488 0.43963267 0.5407161 1455575.6   0.3055     0.137
## 6    Jan 2023 0.5595690 0.38739389 0.1721751  976215.4   0.2600     0.180
## 8    Jul 2023 2.4612086        NaN       NaN 1854459.4   0.6020       NaN
## 10   Jun 2023       NaN        NaN       NaN  479360.2      NaN       NaN
## 12   Mar 2023 0.6657899 0.19022568 0.4755642 1438080.7   0.2100     0.060
## 13   May 2023 1.5069946 0.47093583 1.0360588 1424084.8   0.4800     0.150
## 15   Nov 2023       NaN        NaN       NaN  227433.7      NaN       NaN
## 18   Oct 2023        NA         NA        NA 1200150.1       NA        NA
## 20   Sep 2023       NaN        NaN       NaN  570333.7      NaN       NaN

#We can assume there is no P retained when there is no discharge: 
library(dplyr)
# first make all NaNs into NAs
total23 <- total23 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# now replace NAs where discharge = 0 as also 0
total23$Premoved=ifelse(test = total23$QL == 0, yes = "0", no = total23$Premoved)
#We have one event with loading data but no retention data. Here we can estimate retention based on the average retention of the other months. 
# tell r that Premoved is numeric: 
total23$Premoved=as.numeric(total23$Premoved)
# Calculate percent removal
total23$percentP=total23$Premoved/total23$Pinlb
percentRemoval = mean (total23$percentP,na.rm=T)
print(percentRemoval)

## [1] 0.6188732

# replace the month in question (July 2023)
total23$Premoved=ifelse(test = total23$Month == "Jul" & total23$Year == "2023", yes = total23$Pinlb*percentRemoval, no = total23$Premoved)
Removal = sum (total23$Premoved,na.rm=T)
# This is the removal of P by pool by over the course of the 8 months that we have data for
print(Removal)

## [1] 4.152672

# our low estimate is the lowest concentration we saw: 
lowPcon = min(total23$inflowTP,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon= max(total23$inflowTP,na.rm = TRUE)

# then we can 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)

# 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
total23$Premovedlow=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*lowPcon/453592*lowPper, no = total23$Premoved)
# then for the high estimate
total23$Premovedhigh=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*highPcon/453592*highPper, no = total23$Premoved)

# Now we can sum up our full months: 
lowRemoval = sum (total23$Premovedlow,na.rm=T)
highRemoval = sum (total23$Premovedhigh,na.rm=T)
print(lowRemoval)

## [1] 4.710924

print(highRemoval)

## [1] 8.486876

ASSUMPTIONS IN THIS BUDGET:

Drainage area is accurate (it’s focused mostly on sub-surface drainage)
Curve number is a good choice for this site on average (we chose a higher number to better account for sub-surface drainage)
All inflowing water leaves the wetland as surfacewater outflow
Groundwater effects on water and phosphorus budgets are negligible
Tile drains are accounted for (they’re not)
Nearby rainfall data is accurate for the site
Storm events do not change phosphorus retention dynamics

H4: Rainfall from Ft. Wayne, high res DEM for surface drainage area, expert estimate of subsurface drainage area

# Drainage Area calculated by DEM (164.41 acres, 66.54 hectares)
tDAsur = 164.41
tDAsub = 7.3
# convert tDA to m^2
tDAm2sur=tDAsur*4046.86
tDAm2sub=tDAsub*4046.86

Then we import data on precipitation from weather underground at the Ft. Wayne weather station (31 miles away)

setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))
# Rename so our old code keeps working 
dayrain=raindata

Then we calculate runoff as above.

# Start out by calculating S (inches)
CN=80
S = (1000/CN)-10
# Calculate I (inches)
Ia = 0.2*S
# don't need to convert rain to inches/day, but lets name it so our old code still works
dayrain$precip = dayrain$precipaccum_in
# caculate runoff (inches) for each day in the dataset
dayrain$Q = ifelse(dayrain$precip<=Ia, 0, ((dayrain$precip-Ia)^2)/(dayrain$precip-Ia+S))
# scale up to the entire drainage area
dayrain$Q=tDAm2sur*dayrain$Q
dayrain$fullprecip=tDAm2sur*dayrain$precipaccum_in
# convert back to meters 
dayrain$Qm = dayrain$Q * 0.0254
dayrain$precipm = dayrain$fullprecip * 0.0254
# convert discharge to liters
dayrain$QL = dayrain$Qm * 1000
dayrain$precipL = dayrain$precipm * 1000
# Sum up annual precipitation and discharge
sum(dayrain$precipL)

## [1] 570028366

sum(dayrain$QL)

## [1] 40783055

# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$datetime,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))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

Now we add in the monthly measurements of Phosphorus concentrations to calculate P loading.

# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB P Budget Full inletsandoutlets.csv")
# subset out event sampling 
Pdata<- subset(Pdata, Sampling == "baseflow")
# take the average outflow between streams 1 and 2 (for samples at the end of the stream)
Pdata$outflowTP <- rowMeans(Pdata[c('Outlet_stream1_FORB_17_END', 'Outlet_stream2_FORB_18')], na.rm=TRUE)
# Average the inflow between the main inlet tile and the start of complex 4 (the main inlet tile is not sampled as often as the start of complex 4)
Pdata$inflowTP <- rowMeans(Pdata[c('Inlet_wetland_FORB_9', 'Inlet_tile_FORB_15')], na.rm=TRUE)
# merge the datasets
total <- merge(Pdata,monthrain,by=c("Month","Year"),all=TRUE)
# Measure the mass of P entering the wetland (units = mg P/month)
total$Pin=total$inflowTP*total$QL
# calculate the total mass of P exiting the wetland (units = mg P/month)
total$Pout=total$outflowTP*total$QL
# Convert to lbs P/month
total$Pinlb=total$Pin/453592
total$Poutlb=total$Pout/453592
# calculate the amount of P removed
total$Premoved=total$Pinlb-total$Poutlb

# To start lets look at just 2023
total23<- subset(total, Year == "2023")
total23<- subset(total23, select=c(Month,Year,Pinlb,Poutlb,Premoved,QL,inflowTP,outflowTP))
print(total23)

##    Month Year      Pinlb     Poutlb   Premoved         QL inflowTP outflowTP
## 1    Apr 2023  0.2941601 0.04902668 0.24513341   444762.2   0.3000     0.050
## 3    Aug 2023         NA         NA         NA  2994261.4       NA        NA
## 4    Dec 2023         NA         NA         NA        0.0       NA        NA
## 5    Feb 2023  5.2684047 2.36259064 2.90581404  7822278.9   0.3055     0.137
## 6    Jan 2023  0.1154885 0.07995355 0.03553491   201479.4   0.2600     0.180
## 8    Jul 2023 11.4346318        NaN        NaN  8615710.2   0.6020       NaN
## 10   Jun 2023        NaN        NaN        NaN   108595.3      NaN       NaN
## 12   Mar 2023  4.2681708 1.21947738 3.04869344  9219086.4   0.2100     0.060
## 13   May 2023 10.9287229 3.41522591 7.51349701 10327461.0   0.4800     0.150
## 15   Nov 2023        NaN        NaN        NaN        0.0      NaN       NaN
## 18   Oct 2023         NA         NA         NA   202085.5       NA        NA
## 20   Sep 2023        NaN        NaN        NaN   847334.4      NaN       NaN

#We can assume there is no P retained when there is no discharge: 
library(dplyr)
# first make all NaNs into NAs
total23 <- total23 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# now replace NAs where discharge = 0 as also 0
total23$Premoved=ifelse(test = total23$QL == 0, yes = "0", no = total23$Premoved)

# tell r that Premoved is numeric: 
total23$Premoved=as.numeric(total23$Premoved)
# Calculate percent removal
total23$percentP=total23$Premoved/total23$Pinlb
percentRemoval = mean (total23$percentP,na.rm=T)
print(percentRemoval)

## [1] 0.6188732

# replace the month in question (July 2023)
total23$Premoved=ifelse(test = total23$Month == "Jul" & total23$Year == "2023", yes = total23$Pinlb*percentRemoval, no = total23$Premoved)


# our low estimate is the lowest concentration we saw: 
lowPcon = min(total23$inflowTP,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon= max(total23$inflowTP,na.rm = TRUE)

# then we can 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)


# 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
total23$Premovedlow=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*lowPcon/453592*lowPper, no = total23$Premoved)
# then for the high estimate
total23$Premovedhigh=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*highPcon/453592*highPper, no = total23$Premoved)

# Now we can sum up our full months: 
lowRemoval = sum (total23$Premovedlow,na.rm=T)
highRemoval = sum (total23$Premovedhigh,na.rm=T)
print(lowRemoval)

## [1] 21.41676

print(highRemoval)

## [1] 25.41762

Now we repeat for subsurface runoff!

# Reload our rainfall dataset
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
raindata<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
raindata = mutate(raindata, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))
# Rename so our old code keeps working 
dayrain=raindata
# scale up to the entire drainage area
dayrain$fullprecip=tDAm2sub*dayrain$precipaccum_in
# convert back to meters cubed
dayrain$precipm = dayrain$fullprecip * 0.0254
# convert discharge to liters
dayrain$precipL = dayrain$precipm * 1000
# Sum up annual precipitation and discharge
sum(dayrain$precipL)

## [1] 25309939

# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$datetime,format="%Y-%m-%d")
# Summarize the rainfall by month and year but just for rainfall data
monthrain=dayrain %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
  group_by(Year,Month) %>% summarise(precipL=sum(precipL,na.rm=T))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

# We expect that subsurface drainage will be ~ 31% of total runoff (Pease et al. 2018)
monthrain$QL = monthrain$precipL * 0.31
sum(monthrain$QL)

## [1] 7846081

Now we add in the monthly measurements of Phosphorus concentrations to calculate P loading for subsurface.

# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB P Budget Full inletsandoutlets.csv")
# subset out event sampling 
Pdata<- subset(Pdata, Sampling == "baseflow")
# take the average outflow between streams 1 and 2 (for samples at the end of the stream)
Pdata$outflowTP <- rowMeans(Pdata[c('Outlet_stream1_FORB_17_END', 'Outlet_stream2_FORB_18')], na.rm=TRUE)
# Average the inflow between the main inlet tile and the start of complex 4 (the main inlet tile is not sampled as often as the start of complex 4)
Pdata$inflowTP <- rowMeans(Pdata[c('Inlet_wetland_FORB_9', 'Inlet_tile_FORB_15')], na.rm=TRUE)
# merge the datasets
total <- merge(Pdata,monthrain,by=c("Month","Year"),all=TRUE)
# Measure the mass of P entering the wetland (units = mg P/month)
total$Pin=total$inflowTP*total$QL
# calculate the total mass of P exiting the wetland (units = mg P/month)
total$Pout=total$outflowTP*total$QL
# Convert to lbs P/month
total$Pinlb=total$Pin/453592
total$Poutlb=total$Pout/453592
# calculate the amount of P removed
total$Premoved=total$Pinlb-total$Poutlb

# To start lets look at just 2023
total23<- subset(total, Year == "2023")
total23<- subset(total23, select=c(Month,Year,Pinlb,Poutlb,Premoved,QL,inflowTP,outflowTP))
print(total23)

##    Month Year     Pinlb     Poutlb  Premoved        QL inflowTP outflowTP
## 1    Apr 2023 0.3230811 0.05384686 0.2692343  488490.1   0.3000     0.050
## 3    Aug 2023        NA         NA        NA  576883.5       NA        NA
## 4    Dec 2023        NA         NA        NA  381487.5       NA        NA
## 5    Feb 2023 0.6517419 0.29227049 0.3594714  967675.6   0.3055     0.137
## 6    Jan 2023 0.3720049 0.25754183 0.1144630  648994.0   0.2600     0.180
## 8    Jul 2023 1.6362265        NaN       NaN 1232855.9   0.6020       NaN
## 10   Jun 2023       NaN        NaN       NaN  318681.6      NaN       NaN
## 12   Mar 2023 0.4426212 0.12646319 0.3161580  956044.9   0.2100     0.060
## 13   May 2023 1.0018592 0.31308101 0.6887782  946740.3   0.4800     0.150
## 15   Nov 2023       NaN        NaN       NaN  151199.3      NaN       NaN
## 18   Oct 2023        NA         NA        NA  797867.1       NA        NA
## 20   Sep 2023       NaN        NaN       NaN  379161.3      NaN       NaN

library(dplyr)
# first make all NaNs into NAs
total23 <- total23 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# now replace NAs where discharge = 0 as also 0
total23$Premoved=ifelse(test = total23$QL == 0, yes = "0", no = total23$Premoved)

# tell r that Premoved is numeric: 
total23$Premoved=as.numeric(total23$Premoved)
# Calculate percent removal
total23$percentP=total23$Premoved/total23$Pinlb
percentRemoval = mean (total23$percentP,na.rm=T)
print(percentRemoval)

## [1] 0.6188732

# replace the month in question (July 2023)
total23$Premoved=ifelse(test = total23$Month == "Jul" & total23$Year == "2023", yes = total23$Pinlb*percentRemoval, no = total23$Premoved)

# our low estimate is the lowest concentration we saw: 
lowPcon = min(total23$inflowTP,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon= max(total23$inflowTP,na.rm = TRUE)

# then we can 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)


# 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
total23$Premovedlow=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*lowPcon/453592*lowPper, no = total23$Premoved)
# then for the high estimate
total23$Premovedhigh=ifelse(test = is.na(total23$Premoved) , yes = total23$QL*highPcon/453592*highPper, no = total23$Premoved)

# Now we can sum up our full months: 
lowRemoval = sum (total23$Premovedlow,na.rm=T)
highRemoval = sum (total23$Premovedhigh,na.rm=T)
print(lowRemoval)

## [1] 3.131851

print(highRemoval)

## [1] 5.642127

ASSUMPTIONS IN THIS BUDGET:

Drainage area is accurate (maybe?)
Curve number is a good choice for this site on average (we chose a higher number to better account for sub-surface drainage)
All inflowing water leaves the wetland as surfacewater outflow
Groundwater effects on water and phosphorus budgets are negligible
Tile drains are accounted for (they’re not)
Nearby rainfall data is accurate for the site
Storm events do not change phosphorus retention dynamics

Hybrid Modelling Budget

This is not related to the series of hypotheses, but an attempt to look at how using real % retention in the ODNR budget improves it in comparison to estimating rainfall data for the site.

Calculate phosphorus loading as above:

# units = lbs / year
print(loadFORB)

## [1] 56.60076

Now we’ll calculate the P retention at a couple locations across the wetland.

We’ll start with Wetland Complex 4:

# Add in the P dataset (units here are mg/L for TP)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB P Budget Full inletsandoutlets.csv")
# subset out event sampling 
Pdata<- subset(Pdata, Sampling == "baseflow")
# Calculate percent retention in Wetland Complex 4
Pdata$retainedPcomp4=Pdata$Inlet_wetland_FORB_9-Pdata$Outlet_wetland_FORB_12
Pdata$percentPcomp4 = Pdata$retainedPcomp4 / Pdata$Inlet_wetland_FORB_9 * 100
percentRemoval = mean (Pdata$percentPcomp4,na.rm=T)
print(percentRemoval)

## [1] 8.721188

We get an average percent removal of 8.7% but see here how variable that is:

# make it realize it's a date time value
Pdata$Date <- as.Date(Pdata$Date, format="%m/%d/%Y")
# make the plot
library(ggplot2)
p1<-ggplot(data=Pdata, aes(x=Date, y=percentPcomp4)) +
  geom_bar(stat="identity")+scale_x_date(date_breaks = "1 month", date_labels =  "%b %Y")+theme_classic() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
print(p1)

## Warning: Removed 6 rows containing missing values (`position_stack()`).

Now let’s try out the whole flowpath, starting at the tile draining into Wetland Complex 4 and ending in the two streams that deposit into the Maumee:

# take the average outflow between streams 1 and 2 (for samples at the end of the stream)
Pdata$outflowTP <- rowMeans(Pdata[c('Outlet_stream1_FORB_17_END', 'Outlet_stream2_FORB_18')], na.rm=TRUE)
# inflow remains largely the same for the main flowpath because the inlet tile is not sampled often
Pdata$inflowTP <- rowMeans(Pdata[c('Inlet_wetland_FORB_9', 'Inlet_tile_FORB_15')], na.rm=TRUE)

# Calculate percent retention in whole wetland
Pdata$retainedP=Pdata$inflowTP-Pdata$outflowTP
Pdata$percentP = Pdata$retainedP / Pdata$inflowTP * 100
percentRemoval = mean (Pdata$percentP,na.rm=T)
print(percentRemoval)

## [1] 53.27277

We want to still look at the variability of this measurement:

# make it realize it's a date time value
Pdata$Date <- as.Date(Pdata$Date, format="%m/%d/%Y")
# make the plot
library(ggplot2)
p2<-ggplot(data=Pdata, aes(x=Date, y=percentP)) +
  geom_bar(stat="identity")+scale_x_date(date_breaks = "1 month", date_labels =  "%b %Y")+theme_classic() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
print(p2)

## Warning: Removed 8 rows containing missing values (`position_stack()`).

We see that the full wetland much more consistently retains P than just wetland complex 4.

P retention at the wetland ranges from 19 to 83 % retention, but that doesn’t account for seasonal differences so lets try to take annual averages:

# 2021
Pdata2021<- subset(Pdata, Year == "2021")
percentRemoval2021 = mean (Pdata2021$percentP,na.rm=T)
print(percentRemoval2021)

## [1] 45.16129

# 2022
Pdata2022<- subset(Pdata, Year == "2022")
percentRemoval2022 = mean (Pdata2022$percentP,na.rm=T)
print(percentRemoval2022)

## [1] 41.61901

# 2023
Pdata2023<- subset(Pdata, Year == "2023")
percentRemoval2023 = mean (Pdata2023$percentP,na.rm=T)
print(percentRemoval2023)

## [1] 61.88732

This gives us a range from 42% retention to 62% retention (but our data is not evenly distributed across different years)

We can now apply these measurements of P retention to the estimated phosphorus loading from SPARROW modelling.

# units = lbs / year
print(loadFORB)

## [1] 56.60076

# P Removal Rates:
low = loadFORB*0.42
high = loadFORB*0.62

# units = lbs/year
print(low)

## [1] 23.77232

print(high)

## [1] 35.09247

By this budgetting we can estimate that Forder Bridge retains between 24 and 35 lbs P / year

Lets try to account for the isolated pools:

We’re going to try and use depth data to estimate the mass of phosphorus in the water of each pool and assume that decreasing mass is a result of uptake (either into the sediment or via losses through groundwater).

For this we’re going to need estimates of water volume and phosphorus concentrations at various points in time.

We will use point measurements of depth to estimate water volume until we get some actual depth sensors in the pools.

Point measurements of depth are NOT IDEAL. They are not measured in the same place and they are not consistent across the pool. WE CAN DO BETTER. But we need more/better data to do better and this is just a rough estimate.

# inputting area for each pool as estimated by Bob Midden using Acrobat's measure tool (units are in m2):
aforb1=911.46
aforb2=665.74
aforb3=773.41
aforb5=440.08
aforb6=460.15
aforb7=305.56
aforb8=661.46

We will compare this method to an estimate of inflow based on rainfall and assuming that all P that enters the pools stays in the pools because they are isolated wetlands.

library(dplyr)
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB_Samp master 12.19.csv")
# change all -999 values to NA
Pdata$tp_mgL=na_if(Pdata$tp_mgL, -999)
Pdata$samp_wdepth_cm=na_if(Pdata$samp_wdepth_cm, -999)
site1<- subset(Pdata, location == "FORB_01")
# calculate the volume in m^3
site1$volume = (site1$samp_wdepth_cm/100)*aforb1
# Convert to liters (*1000)
site1$volumeL=site1$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site1$Pmass = site1$volumeL * site1$tp_mgL

Lets plot this out for the first site to see what’s happening:

library(ggplot2)
#convert datetime so R knows it's a date
site1$datetime=as.POSIXct(site1$datetime,format="%m/%d/%Y")
site1$datetime <- as.Date(site1$datetime, format="%Y-%m-%d")

# now plot it out
p1<-ggplot(data=site1, aes(x=datetime, y=Pmass)) +
  geom_bar(stat="identity",width=20)+scale_x_date(date_breaks = "1 month", date_labels =  "%b %Y")+theme_classic()+ylab("P mass in water (mg)")
print(p1)

## Warning: Removed 64 rows containing missing values (`position_stack()`).

## Warning: `position_stack()` requires non-overlapping x intervals

So we’re going to break this down into a minimum and maximum estimate of P being retained by these pools: Our highest value was January 2023: with over 40,000 mg of P in the pool. We know that at least this much P was in the pool at a given time and are assuming that there is not an outflow because these are isolated pools and that most P is retained.

So even if the P is being agitated and recycled and isn’t NEW there was at least this much P in the pool at any one point during the year. That is our minimum estimate of how much P is being retained by this pool.

# take the maximum value and convert to lbs to get the low end estimate of P retention
lowPpool1 = max(site1$Pmass,na.rm=TRUE)*2.20462e-6
lowPpool1

## [1] 0.0964523

Our maximum estimate assumes that EACH P measurement corresponds to a rain event and that ALL P from that event is retained and that the pool dries out in between events. This is likely not the case, but it’s our maximum where we are assuming that all the mass of P in our pool is being retained

# take the sum of all P mass measurements to get the high end P retention in pool 1
highPpool1 = sum(site1$Pmass,na.rm=TRUE)*2.20462e-6
highPpool1

## [1] 0.3997747

Now we can do the rest of the pools:

### Pool 2 ###
site6<- subset(Pdata, location == "FORB_02")
# calculate the volume in m^3
site6$volume = (site6$samp_wdepth_cm/100)*aforb2
# Convert to liters (*1000)
site6$volumeL=site6$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site6$Pmass = site6$volumeL * site6$tp_mgL
# calculate low end estimate
lowPpool2 = max(site6$Pmass,na.rm=TRUE)*2.20462e-6
# calculate high end estimate
highPpool2 = sum(site6$Pmass,na.rm=TRUE)*2.20462e-6

### Pool 3 ###
site3<- subset(Pdata, location == "FORB_03")
# calculate the volume in m^3
site3$volume = (site3$samp_wdepth_cm/100)*aforb3
# Convert to liters (*1000)
site3$volumeL=site3$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site3$Pmass = site3$volumeL * site3$tp_mgL
# calculate low end estimate
lowPpool3 = max(site3$Pmass,na.rm=TRUE)*2.20462e-6
# calculate high end estimate
highPpool3 = sum(site3$Pmass,na.rm=TRUE)*2.20462e-6

### Pool 5 ###
site5<- subset(Pdata, location == "FORB_05")
# calculate the volume in m^3
site5$volume = (site5$samp_wdepth_cm/100)*aforb5
# Convert to liters (*1000)
site5$volumeL=site5$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site5$Pmass = site5$volumeL * site5$tp_mgL
# calculate low end estimate
lowPpool5 = max(site5$Pmass,na.rm=TRUE)*2.20462e-6
# calculate high end estimate
highPpool5 = sum(site5$Pmass,na.rm=TRUE)*2.20462e-6

### Pool 6 ###
site6<- subset(Pdata, location == "FORB_06")
# calculate the volume in m^3
site6$volume = (site6$samp_wdepth_cm/100)*aforb6
# Convert to liters (*1000)
site6$volumeL=site6$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site6$Pmass = site6$volumeL * site6$tp_mgL
# calculate low end estimate
lowPpool6 = max(site6$Pmass,na.rm=TRUE)*2.20462e-6
# calculate high end estimate
highPpool6 = sum(site6$Pmass,na.rm=TRUE)*2.20462e-6

### Pool 7 ###
site7<- subset(Pdata, location == "FORB_07")
# calculate the volume in m^3
site7$volume = (site7$samp_wdepth_cm/100)*aforb7
# Convert to liters (*1000)
site7$volumeL=site7$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site7$Pmass = site7$volumeL * site7$tp_mgL
# calculate low end estimate
lowPpool7 = max(site7$Pmass,na.rm=TRUE)*2.20462e-6
# calculate high end estimate
highPpool7 = sum(site7$Pmass,na.rm=TRUE)*2.20462e-6

### Pool 8 ###
site8<- subset(Pdata, location == "FORB_08")
# calculate the volume in m^3
site8$volume = (site8$samp_wdepth_cm/100)*aforb8
# Convert to liters (*1000)
site8$volumeL=site8$volume * 1000
# calculate mass by multiplying concentration (mg/L) * liters = mg P
site8$Pmass = site8$volumeL * site8$tp_mgL
# calculate low end estimate (and covert to lbs at the same time)
lowPpool8 = max(site8$Pmass,na.rm=TRUE)*2.20462e-6
# calculate high end estimate (and covert to lbs at the same time)
highPpool8 = sum(site8$Pmass,na.rm=TRUE)*2.20462e-6

Now we can sum up all of those values to get an estimate of the low and high end retention for the isolated pools at FORB (excluding Wetland complex 2 which has never been observed to have water in it)

# Wetland complex 1 
Comp1Plow = lowPpool1+lowPpool2+lowPpool3
Comp1Phigh = highPpool1+highPpool2+highPpool3
Comp1Plow

## [1] 0.3227833

Comp1Phigh

## [1] 1.053466

# Wetland Complex 3
Comp3Plow = lowPpool5+lowPpool6+lowPpool7+lowPpool8
Comp3Phigh = highPpool5+highPpool6+highPpool7+highPpool8
Comp3Plow

## [1] 0.3471769

Comp3Phigh

## [1] 1.223496

# Both complexes: 
IsoPlow = Comp1Plow + Comp3Plow
IsoPhigh = Comp1Phigh + Comp3Phigh
IsoPlow

## [1] 0.6699602

IsoPhigh

## [1] 2.276963

ASSUMPTIONS IN THIS BUDGET:

Volume estimates are accurate (they’re not; measurements at staff guages or depth sensors will improve this)
Groundwater dynamics are captured
Measuring this at the monthly scale is sufficient to capture P dynamics at these sites.

Lets try a different method to account for these pools for comparison!

In this new method our goal is to estimate the amount of water flowing into each wetland complex and then use average concentrations in the complex to estimate the amount of P entering the complex during the month that we estimate discharge into the pool. We think this estimate is much more dicey than the previous one but wanted to compare and look for ways that it could be useful because if we could improve this method it would solve the issue with the last one that we are relying on monthly data to capture phosphorus dynamics.

We are using the full complexes as opposed to individual pools because there was a lot of overlap in terms of drainage area between the different pools

# load in the rainfall data
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
dayrain<-read.csv("Ftwaynerainfall.csv")
# make it recognize that it's a date-time value
dayrain = mutate(dayrain, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))

First we need to input the drainage area for each pool as measured by Bishwodeep with DEMs

DAc1 = 10.00 # this is a placeholder number
DAc3 = 10.00 # this is a placeholder number
DAc2 = 10.00 # this is a placeholder number
# convert tDA to m^2
DAc1m2=DAc1*4046.86
DAc3m2=DAc3*4046.86
DAc2m2=DAc2*4046.86

We’ll do the same thing as before and use pool 1 as an example before replicating for the whole batch. Basically the method for calculating runoff is the same as the budget using rainfall to estimate the main flow path (but with a much smaller drainage area).

# Start out by calculating S (inches)
CN=85
S = (1000/CN)-10
# Calculate I (inches)
Ia = 0.2*S
# don't need to convert rain to inches/day, but lets name it so our old code still works
dayrain$precip = dayrain$precipaccum_in
# calculate runoff (inches) for each day in the dataset
dayrain$Q = ifelse(dayrain$precip<=Ia, 0, ((dayrain$precip-Ia)^2)/(dayrain$precip-Ia+S))
# scale up to the entire drainage area
dayrain$Q=DAc1m2*dayrain$Q
# convert back to meters 
dayrain$Qm = dayrain$Q * 0.0254
# convert discharge to liters
dayrain$QL = dayrain$Qm * 1000
# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$datetime,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))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

Now we’re going to use the concentration in the pool at each month to estimate the amount of P entering the pool. This probably isn’t a great estimate because the concentration in the pool is likely already somewhat diluted, OR is being loaded by what’s already in the pool that is released when water enters it, but it might give us a decent estimate. Also not as great when comparing monthly concentrations to continuous rainfall data, but that problem likely isn’t going to be solved unless we’re able to find a relationship between P concentration and turbidity or something like that and then get sensors for turbidity in the pools.

Regardless of whether or not this is a good metric of concentration it should give us a decent estimate of the amount of runoff (where the majority of P should be coming from) that is entering the pool.

# bring in our P data
setwd("C:/Users/kande120/OneDrive - Kent State University/Phosphorus Budgets/FORB/Data")
Pdata<-read.csv("FORB_Samp master 12.19.csv")
# change all -999 values to NA
Pdata$tp_mgL=na_if(Pdata$tp_mgL, -999)
# make the date value usable
Pdata = mutate(Pdata, datetime = as.POSIXct(datetime, format = "%m/%d/%Y"))
Pdata$Date=as.POSIXct(Pdata$datetime,format="%Y-%m-%d")
# subset out complex 1 
site1<- subset(Pdata, location == "FORB_01" |location == "FORB_02" | location == "FORB_03" )
# Summarize the dataset by month and year (this takes the average of all values collected for each month between the different pools)
monthsite1=site1 %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
  group_by(Year,Month) %>% summarise(tp_mgL=mean(tp_mgL,na.rm=T))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

# merge the datasets
totalc1 <- merge(monthsite1,monthrain,by=c("Month","Year"),all=TRUE)
# convert NaNs to NAs
totalc1 <- totalc1 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# calculate the amount of P entering the pool for each month! (units are mg) 
totalc1$massP <- totalc1$tp_mgL*totalc1$QL

## make high and low estimates for months with discharge but no concentration data
# our low estimate is the lowest concentration we saw: 
lowPcon1 = min(totalc1$tp_mgL,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon1= max(totalc1$tp_mgL,na.rm = TRUE)
# first for the low estimate
totalc1$Pmasslow1=ifelse(test = is.na(totalc1$massP) , yes = totalc1$QL*lowPcon1, no = totalc1$massP)
# then for the high estimate
totalc1$Pmasshigh1=ifelse(test = is.na(totalc1$massP) , yes = totalc1$QL*highPcon1, no = totalc1$massP)

# sum up the mass of P for the total retained by the pool (currently just 2023 because of the data included)
sum(totalc1$Pmasslow1,na.rm=T)*2.20462e-6

## [1] 1.1431

sum(totalc1$Pmasshigh1,na.rm=T)*2.20462e-6

## [1] 7.890139

Notes to self for refining this method: add in multiple years precip data, take average of multiple years for estimating pool retention

Now we need to do all of this for wetland complex 3! Wetland complex 2 is always dry so we’re not really sure what we could use as P concentrations for that one.

# calculate runoff (inches) again for each day in the dataset
dayrain$Q = ifelse(dayrain$precip<=Ia, 0, ((dayrain$precip-Ia)^2)/(dayrain$precip-Ia+S))
# scale up to the entire drainage area
dayrain$Q=DAc2m2*dayrain$Q
# convert back to meters 
dayrain$Qm = dayrain$Q * 0.0254
# convert discharge to liters
dayrain$QL = dayrain$Qm * 1000
# Sum up data to get monthly values
library(lubridate)
library(tidyverse)
dayrain$Date=as.POSIXct(dayrain$datetime,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))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

# subset out complex 2
site2<- subset(Pdata, location == "FORB_05" |location == "FORB_06" | location == "FORB_07"|location == "FORB_08" )
# Summarize the dataset by month and year (this takes the average of all values collected for each month between the different pools)
monthsite2=site2 %>% mutate(Month=month(Date,label = T),Year=year(Date)) %>%
  group_by(Year,Month) %>% summarise(tp_mgL=mean(tp_mgL,na.rm=T))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

# merge the datasets
totalc2 <- merge(monthsite2,monthrain,by=c("Month","Year"),all=TRUE)
# convert NaNs to NAs
totalc2 <- totalc2 %>% mutate_all(~ifelse(is.nan(.), NA, .))
# calculate the amount of P entering the pool for each month! (units are mg) 
totalc2$massP <- totalc2$tp_mgL*totalc2$QL

## make high and low estimates for months with discharge but no concentration data
# our low estimate is the lowest concentration we saw: 
lowPcon2 = min(totalc2$tp_mgL,na.rm = TRUE)
# our high estimate is the highest concentration we saw: 
highPcon2= max(totalc2$tp_mgL,na.rm = TRUE)
# first for the low estimate
totalc2$Pmasslow2=ifelse(test = is.na(totalc2$massP) , yes = totalc2$QL*lowPcon2, no = totalc2$massP)
# then for the high estimate
totalc2$Pmasshigh2=ifelse(test = is.na(totalc2$massP) , yes = totalc2$QL*highPcon2, no = totalc2$massP)

# sum up the mass of P for the total retained by the pool (currently just 2023 because of the data included)
sum(totalc2$Pmasslow2,na.rm=T)*2.20462e-6

## [1] 3.81576

sum(totalc2$Pmasshigh2,na.rm=T)*2.20462e-6

## [1] 10.22504

Phosphorus Budgets for FORB

Kenneth Anderson

2023-10-20

H0: ODNR Budgets:

ASSUMPTIONS IN ODNR BUDGETS:

H1: Rainfall at Ft. Wayne Budget, Drainage area from high resolution DEM

ASSUMPTIONS IN THIS BUDGET:

H2: Rainfall from Ft. Wayne, expert estimate of drainage area

ASSUMPTIONS IN THIS BUDGET:

H3: Rainfall from Ft. Wayne, expert estimate of drainage area, estimates of runoff from Pease et al. 2018

ASSUMPTIONS IN THIS BUDGET:

H4: Rainfall from Ft. Wayne, high res DEM for surface drainage area, expert estimate of subsurface drainage area

ASSUMPTIONS IN THIS BUDGET:

Hybrid Modelling Budget

Lets try to account for the isolated pools:

ASSUMPTIONS IN THIS BUDGET:

Lets try a different method to account for these pools for comparison!