---
title: "Anaerobic duration script, lake Arendsee"
author: "Richard LaBrie"
date: "2023-02-08"
output:
  pdf_document: default
  html_document: default
---

This script is used to produce figures 1 and 2 (only Arendsee) and figures S2 and S4 to S9 from the article title "Time span in anaerobic condition predicts biogeochemical consequences of oxygen depletion in lakes".

Note that this script was written using the "Tomorrow Night Bright" editor theme and "Lucida Console" editor font. These settings change the way the script looks like in your R Studio. If you want to use them, go to: Tools > Global Options > Appearance

If you have already worked with R markdown scripts, you can skip to the first chunk of code.
R Markdown scripts (.Rmd) have a few differences with normal .R script.
First, the path to the .Rmd file is stored and you don't need to set another working directory. The relative path can be used as usual with "./", same as going one folder up with "../". The path is only stored when executing lines from the script. If you want to run lines of code from the console, you need to set a working directory.
Second, all chunk of code can be run at once using the "play" symbol on the write of the chunk, and the entire script can be run with the "run" command at the top of script.
Third, plots that you generate from the script won't show in the plot section of R Studio, but under your current chunk of code. If you run the plot command from the console, it will show in the plot screen of R Studio.
Finally, if you click on the little triangle left to a chunk, it will collapse the chunk to a single line. It is much easier to navigate through the script using this feature.


#Loading libraries
```{R Libraries}
library(RColorBrewer)
library(dplyr)
library(segmented)
library(minpack.lm)
library(tidyr)
library(ggplot2)
library(plotrix)

#This step is only to avoid rerunning the script every time (not included with the raw data)
#load("./AnaerobicDuration.RData")
```

#Create sub-folders
To run smoothly, the script will export exploratory and manuscript figures in 2 different folders. This chunk here is simply to create them.
```{R Folder creation}
dir.create("Exploration output")
dir.create("Output")
```

#The three functions below are used later in the script.
The first one "Matrix_Inter.R" is used to interpolate missing days of oxygen profiles.
The second one "O2_DecayRate.R" returns a decay rate for each 1m step based on linear, log-linear or exponential plateau equations. 
The third one "O2_Model.R" is based on "O2_DecayRate.R" but returns the anaerobic duration values. The function works well for lake Arendsee but was too rigid to be used with other lakes. I recommend using "O2_DecayRate.R".
```{R Internal functions}
source("./Matrix_Inter.R", encoding = 'UTF-8')
source("./O2_DecayRate.R") 
source("./O2_Model.R") 
```


#Load data and create objects that will be reused throughout the script
```{R Create Master objects}
Bats <- read.csv("./Arendsee.alpha.csv") #Read Arendsee bathymetry file
Ar.full <- read.csv("./Arendsee.DO-logger.csv")

ar.depths = c(30,35,40,45,47) #Used later in for loops
ar.year = c(2017,2018,2019,2020) #Used later in for loops

#Load all file names
YSI.list = list.files()
#Keep only csv files
YSI.list = YSI.list[grep(pattern = "DO-YSI",list.files())]

#Oxygen threshold for supplementary figures S6 to S9 (2) and for the main story (0.5)
O2.threshold = c(2, 0.5)
```

#These are the metrics used to compare models. I put them here to make it easier to follow the code later
RMSE = sqrt(mean(summary(Model)$residuals^2)) #Used for all three model types
AIC = AIC(Model) #Used for all three model types

R2 for log and exponential plateau (nlsLM function doesn't calculate R2, so we do it manually)
pred <- predict(Model, alpha)
rss <- sum((pred - Observations) ^ 2)
tss <- sum((Observations - mean(Observations)) ^ 2)
rsq <- 1 - rss/tss


#The following chunk of code will calculate Jz for DO loggers.
I use two different ways, a linear regression by selecting dates beforehand, and a segmented regression for a more statistical approach. Comparing the two produces Figure S2. Note that figure S2 is not produced in this chunk of code, but later as it also contains the comparison using the automated YSI profiler. Also note that figures are all produced with Base and requires tens to hundreds of line of code. For this reason, they will always be in the own chunk of code and clear flagged as figures.
The objects created in think chunk will be reused later in the script to compare oxygen sampling methods
This section of the script is heavily annotated as the other sections below using other sampling methods are built identically.
```{R Jz - Loggers}
#Create an output object (years x depths) for linear regression
Jz.mat.day = as.data.frame(matrix(nrow=length(ar.depths), ncol = length(ar.year)))
colnames(Jz.mat.day) = ar.year #Column names are years
rownames(Jz.mat.day) = ar.depths #Row names are depths
Jz.mat.day.err <- Jz.mat.day #Create an "error" matrix 

#Select time stamps 
#These time were selected after looking at full time series (not included in this script), and were aggregated to daily values (division of hourly data by 24)
ar.from.day <- round(matrix(c(2000,1919,2000,1500,
                              2000,707,2000,1500,
                              2000,2000,2000,1500,
                              2000,1000,2000,1500,
                              2000,2000,2000,2000), nrow=5, ncol=4, byrow = T)/24) 
ar.to.day <- round(matrix(c(7500,7400,7300,7300,
                            7200,6000,7000,7200,
                            7000,7000,6800,6000,
                            5000,4800,5000,4600,
                            4500,4400,4800,4400), nrow=5, ncol=4, byrow = T)/24)

#Create a temporary matrix to store segmented regression slopes (Jz)
Ar.seg.dates.day = matrix(nrow=20,ncol=4)
colnames(Ar.seg.dates.day) = c("First date", "First DO", "Second date", "Second DO")
rownames(Ar.seg.dates.day) = as.vector(outer(ar.year, ar.depths, paste, sep=" "))
rownames(Ar.seg.dates.day) = paste0(rownames(Ar.seg.dates.day), "m")

#Create AIC table for comparison
AIC.table.day = matrix(nrow=20, ncol=2)
colnames(AIC.table.day) = c("AIC.lin", "AIC.2")
rownames(AIC.table.day) = rownames(Ar.seg.dates.day)

#Temporary lists to store information in the loop
Jz.seg.list.day = list() #This saves the entire segmented regression output
Jz.seg.coef.day = list() #This saves the coefficients of the segmented regression

{
  counter=1 #Used to store information in the lists
  
  #I store temporary graphs in large pdf.
  #Creating the graphs is not mandatory, but this loop is to calculate Jz
  pdf("./Exploration output/Ar.loggers.SegLm.day.100-200.TEST.pdf", width=8, height=3) #Save exploration analysis in a different folder that figures for the article
  par(mfrow=c(1,4))
  
  for(i in 1:length(ar.depths)){
    #Keep only DO
    Ar.DO <- filter(Ar.full, depth == ar.depths[i]) %>% filter(parameter == "DO_mgL")
    #Create Day of Year variable
    Ar.DO$DOY = as.numeric(strftime(Ar.DO$CET, format = "%j"))
    
    #Aggregate hourly values to daily values
    Ar.DO.day = Ar.DO %>% group_by(year, DOY) %>% summarize(value = mean(value))
    
    #Another for loop to go year by year
    for(j in 1:length(ar.year)){
      #Keep a single year
      Ar.DO.yr.day <- filter(Ar.DO.day, year == ar.year[j])
      #Keep select dates from spring to winter
      Ar.DO.yr.day <- Ar.DO.yr.day[which(Ar.DO.yr.day$DOY>=ar.from.day[i,j] & Ar.DO.yr.day$DOY<=ar.to.day[i,j]),]
      #Change tibble to dataframe
      Ar.DO.yr.day = as.data.frame(Ar.DO.yr.day)
      
      #Linear model
      Dec.Day = Ar.DO.yr.day$DOY
      value = Ar.DO.yr.day$value
      
      lin.mod = lm(value ~ Dec.Day)
      
      #Plot results (not mandatory)
      plot(value~ Dec.Day,
           xlab = "Time",
           ylab = "DO (mg/L)",
           main = paste0(ar.year[j]," ", ar.depths[i], "m"))
      abline(lin.mod, lwd = 2, col = "cyan")
      
      #Store AIC value in the matrix
      AIC.table.day[counter,1] = round(AIC(lin.mod))
      
      #############################################################################
      #This is the critical part to calculate Jz. For your own work, the next few lines are the one you need i.e., the segmented regression using 2 breakpoints. The rest of this for loop was written to automatically calculate Jz for many years and many depths.
      #############################################################################
      #Here I add dummy values at the beginning and end of "day of year" and "DO values". This way I force the segmented regression to find the most relevant ecological breakpoints when DO starts going down and when it stops
      #The other sampling scheme (automated profiler and 2 profiles) use the same logic
      Dec.Day.segLm = c(seq(Dec.Day[1]-400,Dec.Day[1]-1,1),Dec.Day, seq(Dec.Day[length(Dec.Day)]+1,Dec.Day[length(Dec.Day)]+400,1))
      value.segLm = c(rep(max(value),400),value, rep(min(value),400))
      lin.mod.segLm = lm(value.segLm ~ Dec.Day.segLm)
      
      #Piecewise regression (find the inflection points)
      #With 2 breakpoints (the most important one). The others are only for comparison
      psi = c(150,170) #else psi = c(30)
      segmented.mod <- segmented(lin.mod.segLm, psi=psi)
      plot.segmented(segmented.mod, add=T, conf.level=.95, rug=F, col = "chartreuse3", lwd=2)
      # Store AIC value for 2 breakpoints
      AIC.table.day[counter,2] = round(AIC(segmented.mod))
      
      #Store Jz and associated error
      Jz.mat.day[i,j] = summary(lin.mod)$coefficients[2,1]*-1 #Multiply by -1 to have a positive consumption rate per day
      Jz.mat.day.err[i,j] = summary(lin.mod)$coefficients[2,2]
      
      #Save regressions information in list and matrices
      Jz.seg.list.day[[counter]] = segmented.mod
      Ar.seg.dates.day[counter,2] = Ar.DO.yr.day[round(summary(Jz.seg.list.day[[counter]])$psi[1,2]), "value"]
      Ar.seg.dates.day[counter,4] = Ar.DO.yr.day[round(summary(Jz.seg.list.day[[counter]])$psi[2,2]), "value"]
      Jz.seg.coef.day[[counter]] = segmented.mod$coefficients
      legend("topright", legend = c(Ar.seg.dates.day[counter,1],Ar.seg.dates.day[counter,3]))

      counter = counter+1
    }
  }
  dev.off()
}

#Compare AIC. Segmented regression is always much better
AIC.table.day

###Calculate the linear model###
#The list of Jz contains all the information provided by the SEGMENTED REGRESSION
#Some of this information is not useful, hence the names "remove" and "remove1"
Ar.Jz.mat.day = matrix(unlist(Jz.seg.coef.day),nrow=6)
rownames(Ar.Jz.mat.day) = c("Intercept", "First slope", "Second slope", "Third slope", "Remove", "Remove1")
colnames(Ar.Jz.mat.day) = rownames(Ar.seg.dates.day)
Ar.Jz.mat.day <- Ar.Jz.mat.day[-c(5,6),]

#We need the first and second slopes from the segmented package to calculate Jz. This creates a vector with all values (depths and years)
Ar.Jz.day = vector(length=20)
for(i in 1:length(Ar.Jz.day))
{
  Ar.Jz.day[i] = Ar.Jz.mat.day[2,i] + Ar.Jz.mat.day[3,i]
}

#Separate Jz for each year, and multiply by -1 to have positive oxygen consumption rates
Ar.Jz.day.2017 = Ar.Jz.day[c(1,5,9,13,17)]*-1
Ar.Jz.day.2018 = Ar.Jz.day[c(1,5,9,13,17)+1]*-1
Ar.Jz.day.2019 = Ar.Jz.day[c(1,5,9,13,17)+2]*-1
Ar.Jz.day.2020 = Ar.Jz.day[c(1,5,9,13,17)+3]*-1

#With the linear Jz calculation
Ar.Jz.day.2017.lin = Jz.mat.day[,1]
Ar.Jz.day.2018.lin = Jz.mat.day[,2]
Ar.Jz.day.2019.lin = Jz.mat.day[,3]
Ar.Jz.day.2020.lin = Jz.mat.day[,4]

#Jz.mat.day contains Jz from the LINEAR REGRESSION
#Store the bathymetry information in the matrix for ease of use
Jz.mat.day[,5] = Bats %>% slice(c(30,35,40,45,47)) %>% dplyr::select(alpha)


#Below are the linear, log-linear and exponential plateau equations to model Jz as a function of alpha(z)

#Jz ~ Alpha(z), linear regression with linear Jz calculation
#Note. Living comes from the Living(stone) model
Ar.Living.2017.day = summary(lm(Ar.Jz.day.2017.lin ~ Jz.mat.day[,5]))
#R2 = 0.88; AIC = -25.3; RMSE = 0.0106
Ar.Living.2018.day = summary(lm(Ar.Jz.day.2018.lin ~ Jz.mat.day[,5]))
#R2 = 0.97; AIC = -26.12632; RMSE = 0.0097
Ar.Living.2019.day = summary(lm(Ar.Jz.day.2019.lin ~ Jz.mat.day[,5]))
#R2 = 0.8856; AIC = -27.42955; RMSE = 0.008549407
Ar.Living.2020.day = summary(lm(Ar.Jz.day.2020.lin ~ Jz.mat.day[,5]))
#R2 = 0.96; AIC = -28.97831; RMSE = 0.007322751

#Jz ~ Alpha(z), log linear equation with linear Jz calculation
alpha.temp = Jz.mat.day[,5]
Ar.Living.2017.day.lin.log = nlsLM(Ar.Jz.day.2017.lin ~ b*log10(k*alpha.temp),
                                   start = list(b = 0.08, k = 200))
#R2 = 0.97; AIC = -32.2; RMSE = 0.005280211
Ar.Living.2018.day.lin.log = nlsLM(Ar.Jz.day.2018.lin ~ b*log10(k*alpha.temp),
                                   start = list(b = 0.08, k = 200))
#R2 = 0.98; AIC = -28.32138; RMSE = 0.007819951
Ar.Living.2019.day.lin.log = nlsLM(Ar.Jz.day.2019.lin ~ b*log10(k*alpha.temp),
                                   start = list(b = 0.08, k = 200))
#R2 = 0.99; AIC = -40.3215; RMSE = 0.002355297
Ar.Living.2020.day.lin.log = nlsLM(Ar.Jz.day.2020.lin ~ b*log10(k*alpha.temp),
                                   start = list(b = 0.08, k = 200))
#R2 = 0.99; AIC = -38.08294; RMSE = 0.002946219

#Jz ~ Alpha, exponential plateau with linear Jz calculation
Ar.Living.2017.day.lin.exp = nlsLM(Ar.Jz.day.2017.lin ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                   start = list(j.m = max(Ar.Jz.day.2017.lin), b = 0.2, k = 1))
#R2 = 0.986; AIC = -33.92481; RMSE = 0.003655875
Ar.Living.2018.day.lin.exp = nlsLM(Ar.Jz.day.2018.lin ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                   start = list(j.m = max(Ar.Jz.day.2018.lin), b = 0.2, k = 1))
#R2 = 0.997; AIC = -35.16316; RMSE = 0.003230059
Ar.Living.2019.day.lin.exp = nlsLM(Ar.Jz.day.2019.lin ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                   start = list(j.m = max(Ar.Jz.day.2019.lin), b = 0.2, k = 1))
#R2 = 0.999; AIC = -48.7609; RMSE = .000829217
Ar.Living.2020.day.lin.exp = nlsLM(Ar.Jz.day.2020.lin ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                   start = list(j.m = max(Ar.Jz.day.2020.lin), b = 0.2, k = 1))
#R2 = 0.996; AIC = -39.03195; RMSE = 0.00219377


#Repeat the same steps as above, but with Jz calculated from the segmented regression
#Jz ~ Alpha(z), linear regression with linear Jz calculation
Ar.Living.2017.day.SegLm = summary(lm(Ar.Jz.day.2017 ~ Jz.mat.day[,5]))
#R2 = 0.94; AIC = -26.4; RMSE = 0.0095
Ar.Living.2018.day.SegLm = summary(lm(Ar.Jz.day.2018 ~ Jz.mat.day[,5]))
#R2 = 0.91; AIC = -17.27; RMSE = 0.0236
Ar.Living.2019.day.SegLm = summary(lm(Ar.Jz.day.2019 ~ Jz.mat.day[,5]))
#R2 = 0.8829; AIC = -26.45; RMSE = 0.0094
Ar.Living.2020.day.SegLm = summary(lm(Ar.Jz.day.2020 ~ Jz.mat.day[,5]))
#R2 = 0.90; AIC = -21.7; RMSE = 0.0151

#Jz ~ Alpha, log linear, with segmented regression O2 consumption
Ar.Living.2017.day.SegLm.log = nlsLM(Ar.Jz.day.2017 ~ b*log10(k*alpha.temp),
                                     start = list(b = 0.08, k = 100))
Ar.Living.2018.day.SegLm.log = nlsLM(Ar.Jz.day.2018 ~ b*log10(k*alpha.temp),
                                     start = list(b = 0.08, k = 100))
Ar.Living.2019.day.SegLm.log = nlsLM(Ar.Jz.day.2019 ~ b*log10(k*alpha.temp),
                                     start = list(b = 0.08, k = 100))
Ar.Living.2020.day.SegLm.log = nlsLM(Ar.Jz.day.2020 ~ b*log10(k*alpha.temp),
                                     start = list(b = 0.08, k = 100))


#Jz ~ Alpha, exponential plateau, with segmented regression O2 consumption
Ar.Living.2017.day.SegLm.exp = nlsLM(Ar.Jz.day.2017 ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                     start = list(j.m = max(Ar.Jz.day.2017), b = 0.2, k = 10))

Ar.Living.2018.day.SegLm.exp = nlsLM(Ar.Jz.day.2018 ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                     start = list(j.m = max(Ar.Jz.day.2018), b = min(Ar.Jz.day.2018), k = 10),
                                     lower = c(1.4*max(Ar.Jz.day.2018),0,0),
                                     upper = c(2*max(Ar.Jz.day.2018), max(Ar.Jz.day.2018),15)) #R2 = 0.963

Ar.Living.2019.day.SegLm.exp = nlsLM(Ar.Jz.day.2019 ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                     start = list(j.m = max(Ar.Jz.day.2019), b = 0.2, k = 10))
Ar.Living.2020.day.SegLm.exp = nlsLM(Ar.Jz.day.2020 ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                                     start = list(j.m = max(Ar.Jz.day.2020), b = 0.2, k = 10))

#Aggregate Jz (Segemented regression) into a vector for later use
Jz.Logger.SegLm = c(Ar.Jz.day.2017,Ar.Jz.day.2018,Ar.Jz.day.2019,Ar.Jz.day.2020)
```

#We will now repeat most of the previous steps but using automated YSI profilers instead of DO loggers
```{R Jz - YSI, warning=FALSE}
YSI.yrs = c(2017,2018,2019,2020,2021) #Here we don't use ar.years because we also have data from 2021

#Create output list for segmented regression
Jz.YSI.list = list()
Jz.YSI.list.err = list()
#And for linear regression
Jz.YSI.lin.list = list()
Jz.YSI.lin.list.err = list()

#These cutoff were chosen a posteriori using the breakpoint thresholds, then manually adjusted
#Similar procedure as with loggers
YSI.17.from = c(70, 73, 79,79,80,79,79,78,78,79,77,77,79,84,85,86,87,90,100)
YSI.17.to = c(315,314,308,306,306,308,309,306,299,307,307,305,297,279,270,262,250,229,202)

YSI.18.from = c(91,93,94,94,95,96,96,97,96,96,95,105,110,110,110,110,110,110)
YSI.18.to = c(330,328,339,338,336,339,329,326,326,327,328,300,275,260,250,239,223,210)

YSI.19.from = c(100,102,101,102,101,101,101,102,102,101,101,101,99,98,95,95,94,93)
YSI.19.to = c(314,301,308,305,304,300,311,305,300,280,278,270,266,259,250,239,230,210)

YSI.20.from = c(105,107,108,106,106,106,107,107,107,107,107,107,107,108,107,107,105,105)
YSI.20.to = c(253,251,251,253,252,250,250,245,240,233,230,219,212,206,200,195,189,181)

YSI.21.from = c(118,120,123,123,123,120,120,120,120,117,119,119,119,119,120,120,125,125)
YSI.21.to = c(255,250,255,255,255,244,241,230,225,222,222,218,214,205,199,188,183,175)


for(i in 1:length(YSI.yrs)){
  #Read the file
  Ar.YSI = read.csv(paste0("./",YSI.list[i]))
  
  #Remove NA at bottom of profiles with niminum value of the profile
  for(j in 1:ncol(Ar.YSI))
  {
    if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
  }
  
  #Transform to long format
  Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
  #Split Date into Year, months and day
  Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")
  
  Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))
  
  Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                               mm = substring(Date, 6, 7),
                               dd = substring(Date, 9, 10))
  
  #Keep only values below 30 meters to match other methods
  Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]
  
  #Order the dataframe
  Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]
  
  
  #Find Jz for each depth
  Jz.YSI = vector(length = length(unique(Ar.long.deep$Depth_m)))
  Jz.YSI.err = Jz.YSI
  Jz.YSI.lin = Jz.YSI
  Jz.YSI.lin.err = Jz.YSI.lin
  #Note. After looking at the pdf once, the second slope in the piecewise regression is always right
  #Note. The next lines of codes reproduce the steps from the logger section to calculate Jz
  # pdf(paste0("./Exploration output/Ar.YSI.",YSI.yrs[i],".pdf"))
  for(k in 1:length(Jz.YSI)){
    Ar.temp <- Ar.long.deep[Ar.long.deep$Depth_m == unique(Ar.long.deep$Depth_m)[k],]
    
    #All these if() are used to make the for loop run smoothly. each of them constrains data for a single year using the time brackets set above the first loop
    if(YSI.yrs[i] == 2017) {
      Ar.temp = Ar.temp[Ar.temp$DOY>YSI.17.from[k] & Ar.temp$DOY < YSI.17.to[k],]
      psi = c(100,200)}
    if(YSI.yrs[i] == 2018) {
      Ar.temp = Ar.temp[Ar.temp$DOY>YSI.18.from[k] & Ar.temp$DOY <  YSI.18.to[k],] 
      psi = c(110,200)}
    if(YSI.yrs[i] == 2019) {
      Ar.temp = Ar.temp[Ar.temp$DOY>YSI.19.from[k] & Ar.temp$DOY < YSI.19.to[k],]
      psi = c(100,200)}
    if(YSI.yrs[i] == 2020) {
      Ar.temp = Ar.temp[Ar.temp$DOY>YSI.20.from[k] & Ar.temp$DOY <  YSI.20.to[k],]
      psi = c(110,200)}
    if(YSI.yrs[i] == 2021) {
      Ar.temp = Ar.temp[Ar.temp$DOY>YSI.21.from[k] & Ar.temp$DOY < YSI.21.to[k],]
      psi = c(110,200)}
    
    #Store values in objects for ease of use
    Dec.Day = Ar.temp$DOY
    value = Ar.temp$DO_mgL
    
    #Remove any NA
    if(any(is.na(Dec.Day))) {
      Dec.Day.na = which(is.na(Dec.Day))
      Dec.Day = Dec.Day[-Dec.Day.na]
      value = value[-Dec.Day.na]
    }
    if(any(is.na(value))) {
      value.na = which(is.na(value))
      Dec.Day = Dec.Day[-value.na]
      value = value[-value.na]
    }
    
    #Linear regression. Same as with loggers
    lin.mod = lm(value[value != min(value)] ~ Dec.Day[value != min(value)])
    print(c(length(value[value != min(value)]),unique(Ar.long.deep$Depth_m)[k]))
    
    #Adding dummy values at the beginning and end, as with loggers
    Dec.Day.segLm = c(seq(Dec.Day[1]-200,Dec.Day[1]-1,1),Dec.Day, seq(Dec.Day[length(Dec.Day)]+1,Dec.Day[length(Dec.Day)]+300,1))
    value.dummy = max(value)
    value.segLm = c(rep(value.dummy,200),value, rep(min(value),300))
    lin.mod.segLm = lm(value.segLm ~ Dec.Day.segLm)
    
    plot(value ~ Dec.Day, las = 1,
         main = paste(unique(Ar.long.deep$Depth_m)[k], "m", sep=" "),
         ylab = expression(DO~mg~L^-1),
         xlab = "Day of year",
         ylim = c(0,13))
    abline(lin.mod, lwd = 2, col = "cyan")
    
    #Piecewise model
    segmented.mod <- segmented(lin.mod.segLm, psi=psi)
    plot.segmented(segmented.mod, add=T, conf.level=.95, rug=F, col = "chartreuse3", lwd=2)
    legend("topright", legend = c(paste("First breakpoint is", round(segmented.mod$psi[1,2]), "day", sep = " "),
                                  paste("Second breakpoint is", round(segmented.mod$psi[2,2]), "day", sep = " ")))
    
    Jz.YSI[k] = (segmented.mod$coefficients[2] + segmented.mod$coefficients[3])*-1
    Jz.YSI.err[k] = summary(segmented.mod)$coef[2,2] + summary(segmented.mod)$coef[3,2]
    Jz.YSI.lin[k] = summary(lin.mod)$coef[2,1] * -1
    Jz.YSI.lin.err[k] = summary(lin.mod)$coef[2,2]
  }
  dev.off()
  Jz.YSI.list[[i]] = Jz.YSI
  Jz.YSI.list.err[[i]] = Jz.YSI.err
  Jz.YSI.lin.list[[i]] = Jz.YSI.lin
  Jz.YSI.lin.list.err[[i]] = Jz.YSI.lin.err
  
  #Select the right alphas
  #Here again, if() are used to make the for loop run smoothly with the differences in raw data
  #2017: 33:48
  #2018-2020: 30:47
  if(i == 1){
    alpha.YSI.list=list()
    alpha.YSI.list[[1]] = Bats %>% slice(c(33:48)) %>% dplyr::select(alpha)
  }
  if(i == 2) alpha.YSI.list[[2]] = Bats %>% slice(c(30:47)) %>% dplyr::select(alpha)
}



#Jz~alpha(z) relationships with different equations. Same as with loggers
#Note. Sometimes there are warnings "NaNs produced". This is normal and does not impact the results

#Linear regressions
Ar.Jz.alpha.YSI.2017 = lm(Jz.YSI.list[[1]][-c(1:3)] ~ alpha.YSI.list[[1]][,1]) 
#R2 = 0.993; AIC = -146.5; RMSE = 0.002059
Ar.Jz.alpha.YSI.2018 = lm(Jz.YSI.list[[2]] ~ alpha.YSI.list[[2]][,1]) 
#R2 = 0.966; AIC = -119.55; RMSE = 0.007398443
Ar.Jz.alpha.YSI.2019 = lm(Jz.YSI.list[[3]] ~ alpha.YSI.list[[2]][,1]) 
#R2 = 0.90; AIC = -133.856; RMSE = 0.004972443
Ar.Jz.alpha.YSI.2020 = lm(Jz.YSI.list[[4]] ~ alpha.YSI.list[[2]][,1]) 
#R2 = 0.86; AIC = -109.9361; RMSE = 0.009663632
Ar.Jz.alpha.YSI.2021 = lm(Jz.YSI.list[[5]] ~ alpha.YSI.list[[2]][,1]) 
#R2 = 0.98; AIC = -105.2641; RMSE = 0.01100277

#Log-linear and exponential plateau regressions
#2017
temp = Jz.YSI.list[[1]][-c(1:3)]
temp.max = max(temp)
alpha = alpha.YSI.list[[1]][,1]
Ar.Jz.alpha.YSI.2017.log = nlsLM(temp ~ b*log10(k*alpha),
                                 start = list(b=0.08, k=200))
#R2 = 0.876; AIC = -100.1125; RMSE = 0.008782861
Ar.Jz.alpha.YSI.2017.exp = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*alpha),
                                 start = list(t.m = temp.max, b=0.05, k=5))
#R2 = 0.9926; AIC = -144.405; RMSE = 0.002067128


#2018
temp = Jz.YSI.list[[2]]
temp.max = max(temp)
alpha = alpha.YSI.list[[2]][,1]
Ar.Jz.alpha.YSI.2018.log = nlsLM(temp ~ b*log10(k*alpha),
                                 start = list(b=0.08, k=200)) 
#R2 = 0.9528; AIC = -113.2695; RMSE = 0.008808999
Ar.Jz.alpha.YSI.2018.exp = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*alpha),
                                 start = list(t.m = temp.max, b=0.05, k=5))
#R2 = 0.9930; AIC = 145.7226; RMSE = 0.003382937



#2019
temp = Jz.YSI.list[[3]]
temp.max = max(temp)
Ar.Jz.alpha.YSI.2019.log = nlsLM(temp ~ b*log10(k*alpha),
                                 start = list(b=0.08, k=200))
#R2 = 0.999; AIC = -221.2; RMSE = 0.00043
Ar.Jz.alpha.YSI.2019.exp = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*alpha),
                                 start = list(t.m = temp.max, b=0.05, k=5))
#R2 = 0.993; AIC = -180.7806; RMSE = 0.00127



#2020
temp = Jz.YSI.list[[4]]
temp.max = max(temp)
Ar.Jz.alpha.YSI.2020.log = nlsLM(temp ~ b*log10(k*alpha),
                                 start = list(b=0.08, k=200))
#R2 = 0.9918; AIC = ; RMSE = 
Ar.Jz.alpha.YSI.2020.exp = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*alpha),
                                 start = list(t.m = temp.max, b=0.05, k=5))
#R2 = 0.997; AIC = -179.1 ; RMSE = 0.0013

#Test with shorter dataset, from 5 to 15 depths in hypolimnion. 
Ar.Jz.alpha.YSI.2020.log.short = list()
for(i in 5:15){
  temp = Jz.YSI.list[[4]][1:i]
  temp.max = max(temp)
  alpha.short = alpha[1:i]
  Ar.Jz.alpha.YSI.2020.log.short[[i-4]] = nlsLM(temp ~ b*log10(k*alpha.short),
                                                start = list(b=0.08, k=200))
}

#Repeat, but with exponential plateau model
Ar.Jz.alpha.YSI.2020.exp.short = list()
for(i in 5:15){
  temp = Jz.YSI.list[[4]][1:i]
  temp.max = max(temp)
  temp.min = min(temp)
  alpha.short = alpha[1:i]
  Ar.Jz.alpha.YSI.2020.exp.short[[i-4]] = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*alpha.short),
                                                start = list(t.m = temp.max, b=temp.min, k=5),
                                                lower = c(1.4*temp.min, 0,0),
                                                upper = c(2*temp.max, temp.max,15))
}


#2021
temp = Jz.YSI.list[[5]]
temp.max = max(temp)
Ar.Jz.alpha.YSI.2021.log = nlsLM(temp ~ b*log10(k*alpha), 
                                 start = list(b=0.08, k=200))
#R2 = 0.937, AIC = -81.67; RMSE = 0.0212
Ar.Jz.alpha.YSI.2021.exp = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*alpha),
                                 start = list(t.m = temp.max, b=0.05, k=5))
#R2 = 0.994, AIC = -122.2; RMSE = 0.0065

rm(temp)
rm(temp.max)
rm(alpha)

```

#We will now print two figures; each in their own chunk of code
#Figure S2
```{R Figure S2}
  # png("./Output/Fig.S2.png", width = 8.36, height = 4.4, units = "in", res=300)
  #Remove comment from previous line to store the graph in a png file.
  #Make sure to remove the comment mark (#) below before dev.off()
  par(mfrow=c(1,2))
  plot(Jz.YSI.lin.list[[1]] ~ Jz.YSI.list[[1]],
       ylab = "Jz (linear fit)", 
       xlab = "Jz (segmented fit)",
       ylim = c(0.03, 0.3),
       xlim = c(0.03, 0.45),
       las = 1, pch = 16, col = "blue",
       main = "YSI profiler")
  
  #Add points from the different years
  for(i in 2:5){
    points(Jz.YSI.lin.list[[i]] ~ Jz.YSI.list[[i]],
           pch = 16, col = c("chartreuse3","red","pink","#73BAE6")[i-1])  }
  #Add regression line
  abline(lm(unlist(Jz.YSI.lin.list) ~ unlist(Jz.YSI.list)))
  #Add 1:1 line
  abline(0,1, lty = 2)
  
  legend("topleft", legend = c("2017","2018","2019","2020","2021"),
         text.col = c("blue", "chartreuse3","red","pink","#73BAE6"))
  
  legend("bottomright", lty = c(1,2), legend = c("regression fit", "1:1 fit"))
  
  #Repeat steps, but with loggers objects
  plot(c(as.matrix(Jz.mat.day[,-5])) ~ Jz.Logger.SegLm,
       xlim = c(0.025, 0.25),
       ylim = c(0.02, 0.2),
       ylab = "Jz (linear fit)", 
       xlab = "Jz (segmented fit)",
       las = 1, pch = 16,
       col = c("blue","blue","blue","blue","blue",
                     "chartreuse3","chartreuse3","chartreuse3","chartreuse3","chartreuse3",
                     "red","red","red","red","red",
                     "pink","pink","pink","pink","pink"),
                     main = "Loggers")
  abline(lm(c(as.matrix(Jz.mat.day[,-5])) ~Jz.Logger.SegLm))
  abline(0,1, lty=2)
  # dev.off()

summary(lm(c(as.matrix(Jz.mat.day[,-5]))~Jz.Logger.SegLm)) #R2 = 0.9606


```

#Figure S4
```{R Figure S4}
#Comparison of different thickness of oxic hypolimnion
#This simulates the second profile being anoxic
#png("./Output/Fig S4.png", res = 300, units = "in", height = 4, width = 8)
par(mfrow=c(1,2))
par(mar=c(4,5,3,1)+.1)
plot(Jz.YSI.list[[4]] ~ alpha.YSI.list[[2]][,1],
     ylim = c(0.05,0.2),
     xlim = c(0, 0.25),
     xlab = expression(alpha(z)),
     ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
     main = "Log-linear",
     las = 1)
#Add the main relationship curve
curve(coef(Ar.Jz.alpha.YSI.2020.log)[1]*log10(x*coef(Ar.Jz.alpha.YSI.2020.log)[2]),
      add=T, col ="blue", lwd = 2, to = 0.235)
#Add the partial anoxic simulations
for(i in 1:11){
  curve(coef(Ar.Jz.alpha.YSI.2020.log.short[[i]])[1]*log10(x*coef(Ar.Jz.alpha.YSI.2020.log.short[[i]])[2]),
        add=T, col ="black", lwd = 1, to = 0.235)
}
legend("topleft", lty = 1, col= c("blue", "black"), legend = c("Full profile", "Short profiles"), cex = 0.6)
mtext(text = "a", side = 3, at = -0.01)

#Repeat with exponential platerau equation
plot(Jz.YSI.list[[4]] ~ alpha.YSI.list[[2]][,1],
     ylim = c(0.05,0.2),
     xlim = c(0, 0.25),
     xlab = expression(alpha(z)),
     ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
     main = "Exponential plateau",
     las = 1)
curve(coef(Ar.Jz.alpha.YSI.2020.exp)[1] - (coef(Ar.Jz.alpha.YSI.2020.exp)[1] - 
                                             coef(Ar.Jz.alpha.YSI.2020.exp)[2])*
        exp(-1*coef(Ar.Jz.alpha.YSI.2020.exp)[3]*x), 
      add=T, col = "Blue", lwd = 2, to = 0.235)
for(i in 1:11){
  curve(coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[1] - (coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[1] - 
                                                          coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[2])*
          exp(-1*coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[3]*x), 
        add=T, col = "Black", to = 0.235)
  # text(x = c(rep(c(0.245,0.255),5),0.255)[i], y =  coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[1] - (coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[1] - 
  #                                                                             coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[2])*
  #                             exp(-1*coef(Ar.Jz.alpha.YSI.2020.exp.short[[i]])[3]*0.25), 
  #      labels = seq(5,15,1)[i], cex = 0.5)

# legend("topleft", legend=("n: number of depths considered"), cex = 0.5)
}
  mtext(text = "b", side = 3, at = -0.01)

#dev.off()

```

#In the following chunk of code, we will calculate Jz for every pair of oxic days taken 28 days apart. We did this to emulate a more normal sampling scheme
```{R Jz - 2 oxic profiles, warning=FALSE}

#These are start and end day of year to simulate using only 2 oxic profiles
#Dates were chosen as for loggers (first scan to make sure we excluded anoxic profiles)
DOY.ini = c(94,122,108,114,115)
DOY.end = c(200,182,201,182, 169)

#Create R2 matrix
R2.mat = matrix(nrow=5, ncol = 3)
rownames(R2.mat) = c("2017", "2018", "2019", "2020", "2021")
colnames(R2.mat) = c("Lin", "Log-lin", "Exp")
AIC.mat <- RMSE.mat <- R2.mat

#Create output lists, as we did before
Jz.lowres.slice.list = list()
ar.alpha.list = list()
Ar.Jz.lm.best.list = list()
Jz.lowres = list()

#The structure of this for loop is very similar to the ones we did before with loggers and YSI. For this one, I also automated the process of storing RMSE, AIC and R2 into matrices
for(i in 1:length(YSI.list)){
  Ar.YSI = read.csv(paste0("./",YSI.list[i]))
  
  #Remove NA at bottom of profiles with niminum value of the profile
  for(j in 1:ncol(Ar.YSI))
  {
    if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
  }
  
  #Transform to long format
  Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
  #Split Date into Year, months and day
  Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")
  
  Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))
  
  Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                               mm = substring(Date, 6, 7),
                               dd = substring(Date, 9, 10))
  
  #Keep only values below 30 meters to match other methods
  Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]
  
  #Order the dataframe
  Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]
  
  #2 oxic profiles - 28 days apart
  #Keep only profiles within 2 dates
  Ar.long.deep = Ar.long.deep[Ar.long.deep$DOY >= DOY.ini[i] & Ar.long.deep$DOY <= DOY.end[i],]
  
  #Keep only values higher than 2mgO2/L
  Ar.long.deep = Ar.long.deep[Ar.long.deep$DO_mgL >=2,]
  
  Depths = unique(Ar.long.deep$Depth_m)
  Jz.lowres = matrix(nrow=length(Depths), ncol = 2)
  colnames(Jz.lowres) = c("Slope", "Error") 
  rownames(Jz.lowres) = Depths
  for(j in 1:length(Depths)){
    Ar.long.deep.temp = Ar.long.deep[Ar.long.deep$Depth_m == Depths[j],]
    Ar.long.deep.temp = Ar.long.deep.temp %>% group_by(DOY) %>% summarize(DO_mgL = mean(DO_mgL))
    
    #Fill in missing days (linear interpolation)
    Ar.long.filled = t(Ar.long.deep.temp)
    colnames(Ar.long.filled) = Ar.long.filled[1,]
    Ar.long.filled = t(as.matrix(Ar.long.filled[-1,]))
    Ar.long.filled = interpolation_fct(Ar.long.filled)
    
    #Calculate Jz between each pair of points
    Slope.temp = vector(length=length(Ar.long.filled)-28)
    for(k in 1:(length(Ar.long.filled)-28)){
      Slope.temp[k] = (Ar.long.filled[,k]-Ar.long.filled[,k+28])/28
    }
    Jz.lowres[j,1] = mean(Slope.temp)
    Jz.lowres[j,2] = sd(Slope.temp)/sqrt(length(Slope.temp))
    # print(c(length(Slope.temp), Depths[j]))
  }
  slice.beg = c(30,30,30,30,30)
  slice.end = c(47,47,47,47,47)
  
  # Jz.lowres.end = c(8,3,3,3,4)
  
  
  ar.alpha <- Bats%>% slice(slice.beg[i]:slice.end[i]) %>% dplyr::select(alpha)
  # if(i == 1) ar.alpha = as.matrix(ar.alpha[-11,])
  ar.alpha.list[[i]] = ar.alpha
  Jz.lowres.slice = Jz.lowres #[-c(1:Jz.lowres.end[i])]# Remove the first few  values because they are identical (likely because sediments to volume is too low)
  if(i == 1) Jz.lowres.slice = Jz.lowres.slice[-19,]
  Jz.lowres.slice.list[[i]] = Jz.lowres.slice
  
  plot(Jz.lowres.slice[,1] ~ ar.alpha[,1])
  Ar.Jz.lowres = lm(Jz.lowres.slice[,1] ~ ar.alpha[,1])
  # Ar.Jz.lowres.log = lm(Jz.lowres.slice ~ log10(ar.alpha[,1]))
  
  #Save relevant comparison metrics
  R2.mat[i,1] = summary(Ar.Jz.lowres)$r.squared
  AIC.mat[i,1] = AIC(Ar.Jz.lowres)
  RMSE.mat[i,1] = sqrt(mean(Ar.Jz.lowres$residuals^2))
  
  temp = Jz.lowres.slice[,1]
  temp.max = max(Jz.lowres.slice[,1])
  temp.min = min(Jz.lowres.slice[,1])
  Ar.Jz.lowres.log = nlsLM(temp ~ b*log10(k*ar.alpha[,1]),
                           start = list(b=0.08, k=200))
  
  #Calculate R2
  pred <- predict(Ar.Jz.lowres.log, ar.alpha[,1])
  rss <- sum((pred - temp) ^ 2)
  tss <- sum((temp - mean(temp)) ^ 2)
  rsq <- 1 - rss/tss
  
  #Save relevant comparison metrics
  R2.mat[i,2] = rsq
  AIC.mat[i,2] = AIC(Ar.Jz.lowres.log)
  RMSE.mat[i,2] = sqrt(mean(summary(Ar.Jz.lowres.log)$residuals^2))
  
  Ar.Jz.lowres.exp = nlsLM(temp ~ t.m - (t.m - b) * exp(-k*ar.alpha[,1]),
                           start = list(t.m = 1.6*temp.max, b=temp.min, k=5))
  pred <- predict(Ar.Jz.lowres.exp, ar.alpha[,1])
  rss <- sum((pred - temp) ^ 2)
  tss <- sum((temp - mean(temp)) ^ 2)
  rsq <- 1 - rss/tss
  
  #Save relevant comparison metrics
  R2.mat[i,3] = rsq
  AIC.mat[i,3] = AIC(Ar.Jz.lowres.exp)
  RMSE.mat[i,3] = sqrt(mean(summary(Ar.Jz.lowres.exp)$residuals^2))
  
  Ar.Jz.lm.best.list[[i]] = list(Ar.Jz.lowres, Ar.Jz.lowres.log, Ar.Jz.lowres.exp)
}
```

#The following chunk is very similar to the previous one. We will again calculate Jz using profiles taken 28 days apart, but we will do it only for 2018. This will be used for figure S5, bottom row. I chose to rewrite the code above and modify it slightly rather than adding conditions in the previous loop to store more information only for year 2018
```{R Jz - 2 oxic profiles 2018}
DOY.ini.2018 = 122
DOY.end.2018 = 182

Jz.lowres.slice.list.2018 = list()
ar.alpha.list.2018 = list()

Ar.Jz.lm.best.list.2018 = list()
Jz.lowres.2018 = list()
Ar.YSI = read.csv(paste0("./",YSI.list[2]))

#Remove NA at bottom of profiles with niminum value of the profile
for(j in 1:ncol(Ar.YSI))
{
  if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
}

#Transform to long format
Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
#Split Date into Year, months and day
Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")

Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))

Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                             mm = substring(Date, 6, 7),
                             dd = substring(Date, 9, 10))

#Keep only values below 30 meters to match other methods
Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]

#Order the dataframe
Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]

#2 oxic profiles - 28 days apart
#Keep only profiles within 2 dates
Ar.long.deep = Ar.long.deep[Ar.long.deep$DOY >= DOY.ini.2018 & Ar.long.deep$DOY <= DOY.end.2018,]

#Keep only values higher than 2mgO2/L
Ar.long.deep = Ar.long.deep[Ar.long.deep$DO_mgL >=2,]

Depths = unique(Ar.long.deep$Depth_m)
Jz.lowres.2018 = matrix(nrow = length(Depths), ncol = 33)
rownames(Jz.lowres.2018) = Depths
for(j in 1:length(Depths)){
  Ar.long.deep.temp = Ar.long.deep[Ar.long.deep$Depth_m == Depths[j],]
  Ar.long.deep.temp = Ar.long.deep.temp %>% group_by(DOY) %>% summarize(DO_mgL = mean(DO_mgL))
  
  #Fill in missing days (linear interpolation)
  Ar.long.filled = t(Ar.long.deep.temp)
  colnames(Ar.long.filled) = Ar.long.filled[1,]
  Ar.long.filled = t(as.matrix(Ar.long.filled[-1,]))
  Ar.long.filled = interpolation_fct(Ar.long.filled)
  
  #Calculate Jz between each pair of points
  Slope.temp = vector(length=length(Ar.long.filled)-28)
  for(k in 1:(length(Ar.long.filled)-28)){
    Slope.temp[k] = (Ar.long.filled[,k]-Ar.long.filled[,k+28])/28
  }
  
  Jz.lowres.2018[j,] = Slope.temp
  # print(c(length(Slope.temp), Depths[j]))
}

```

#Figure S5
```{R Figure S5, warning=FALSE}
ar.alpha <- Bats%>% slice(30:47) %>% dplyr::select(alpha)
# This is used to manually permutate loggers Jz to use each triplets
matrix.perm = matrix(c(1,2,3,1,2,4,1,2,5,2,3,4,2,3,5,3,4,5), ncol = 3, byrow = T)

#pdf("./Output/Fig.S5.pdf", width = 4.5, height = 3)
  par(mfrow=c(2,3))
  par(mar=c(4,5,1,1)+.1)
  #Three loggers, linear relationship
  plot(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5], las = 1,
       xlab = "",
       ylab = expression(Jz~(mg~O[2]~L^-1~d^-1)),
       ylim = c(0.04, 0.2), xaxt = "n",
       cex = 1.4, cex.axis=0.8, cex.lab = 0.9)
  abline(lm(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5]), col = "Blue", lwd = 2)
  for(i in 1:nrow(matrix.perm)){
    abline(lm(Ar.Jz.day.2018.lin[matrix.perm[i,]] ~ Jz.mat.day[matrix.perm[i,], 5]), lty = 1, col = alpha("black", 0.4))
  }
  
  #3 loggers, log-linear relationship
  plot(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5], las = 1,
       xlab = "",
       ylab = "", yaxt = "n",xaxt = "n",
       ylim = c(0.04, 0.2),
       cex = 1.4, cex.axis=0.8, cex.lab = 0.9)
  Figure.temp <- Jz.mat.day[,5]
  main = nlsLM(Ar.Jz.day.2018.lin ~ b*log10(k*Figure.temp),
               start = list(b = 0.08, k = 200))
  curve(coef(main)[1]*log10(x*coef(main)[2]), add=T, col = "Blue", lwd = 2)
  for(i in 1:nrow(matrix.perm)){
    alpha.temp = Jz.mat.day[matrix.perm[i,],5]
    Jz.temp = Ar.Jz.day.2018.lin[matrix.perm[i,]]
    temp = nlsLM(Jz.temp ~ b*log10(k*alpha.temp),
                 start = list(b = 0.08, k = 200))
    curve(coef(temp)[1]*log10(x*coef(temp)[2]), add=T, lty = 1, col = alpha("black", 0.4))
  }
  
  #3 loggers, exponential plateau
  plot(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5], las = 1,
       xlab = "",
       ylab = "", yaxt = "n",xaxt = "n",
       ylim = c(0.04, 0.2),
       cex = 1.4, cex.axis=0.8, cex.lab = 0.9)
  
  main = nlsLM(Ar.Jz.day.2018.lin ~ j.m - (j.m-b)*exp(-k*Figure.temp),
               start = list(j.m = max(Ar.Jz.day.2018.lin), b = 0.2, k = 1))
  curve(coef(main)[1] - (coef(main)[1] - coef(main)[2]) * exp(-coef(main)[3]*x),
        add=T, col = "Blue", lwd = 2)
  
  for(i in 1:nrow(matrix.perm)){
    alpha.temp = Jz.mat.day[matrix.perm[i,],5]
    Jz.temp = Ar.Jz.day.2018.lin[matrix.perm[i,]]
    temp = nlsLM(Jz.temp ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                 start = list(j.m = max(Jz.temp), b = 0.2, k = .9))
    curve(coef(temp)[1] - (coef(temp)[1] - coef(temp)[2]) * exp(-coef(temp)[3]*x),
          add=T, lty = 1, col = alpha("black", 0.4))
  }
  
  #2 oxic profiles, linear relationship
  plot(Jz.lowres.slice.list[[2]][,1] ~ ar.alpha$alpha, las = 1,
       xlab = expression(alpha(z)),
       ylab = expression(Jz~(mg~O[2]~L^-1~d^-1)),
       ylim = c(0.04, 0.2),
       cex = 1.4, cex.axis=0.8, cex.lab = 0.9)
  abline(Ar.Jz.lm.best.list[[2]][[1]], col = "Blue", lwd = 2)
  
  for(i in 1:ncol(Jz.lowres.2018)){
    abline(lm(Jz.lowres.2018[,i] ~ ar.alpha$alpha), lty = 1, col = alpha("black", 0.4))
  }
  
  #2 oxic profiles, log-linear
  plot(Jz.lowres.slice.list[[2]][,1] ~ ar.alpha$alpha, las = 1,
       xlab = expression(alpha(z)),
       ylab = "", yaxt = "n",
       ylim = c(0.04, 0.2),
       cex = 1.4, cex.axis=0.8, cex.lab = 0.9)
  Jz.temp = Jz.lowres.slice.list[[2]][,1]
  alpha.temp = ar.alpha$alpha
  main = nlsLM(Jz.temp ~ b*log10(k*alpha.temp),
               start = list(b = 0.08, k = 200))
  curve(coef(main)[1]*log10(x*coef(main)[2]), add=T, col = "Blue", lwd = 2)
  
  for(i in 1:ncol(Jz.lowres.2018)){
    Jz.temp = Jz.lowres.2018[,i]
    temp = nlsLM(Jz.temp ~ b*log10(k*alpha.temp),
                 start = list(b = 0.08, k = 200))
    curve(coef(temp)[1]*log10(x*coef(temp)[2]), add=T, lty = 1, col = alpha("black", 0.4))
  }
  
  #2 oxic profiles, exponential plateau
  plot(Jz.lowres.slice.list[[2]][,1] ~ ar.alpha$alpha, las = 1,
       xlab = expression(alpha(z)),
       ylab = "", yaxt = "n",
       ylim = c(0.04, 0.2),
       cex = 1.4, cex.axis=0.8, cex.lab = 0.9)
  Jz.temp = Jz.lowres.slice.list[[2]][,1]
  
  main = nlsLM(Jz.temp ~ j.m - (j.m-b)*exp(-k*alpha.temp),
               start = list(j.m = max(Ar.Jz.day.2018.lin), b = 0.2, k = 1))
  curve(coef(main)[1] - (coef(main)[1] - coef(main)[2]) * exp(-coef(main)[3]*x),
        add=T, col = "Blue", lwd = 2)
  
  for(i in 1:ncol(Jz.lowres.2018)){
    Jz.temp = Jz.lowres.2018[,i]
    temp = nlsLM(Jz.temp ~ j.m - (j.m-b)*exp(-k*alpha.temp),
                 start = list(j.m = max(Jz.temp), b = 0.3, k = 1))
    curve(coef(temp)[1] - (coef(temp)[1] - coef(temp)[2]) * exp(-coef(temp)[3]*x),
          add=T, lty = 1, col = alpha("black", 0.4))
  }
  
  #dev.off()
```


#Notice. The following chunk is particularly long. It includes three different steps, calculating AD for all data source - equations - years triplets, comparing them to observations and producing figures S6 to S9. I normally write different chunk for figures, but as they were produced using a different oxygen threshold (2mg/L), I print them here in a for loop to avoid duplicating the two previous sections. All will be annonated clearly. If you want to run the code manually, set "z" in the foor for at 1 for "2mg/L" and at 2 for "0.5mg/L".
#We have now calculated Jz using all the sampling scenario. We will proceed in the next chunk to calculate the Anaerobic duration (AD). We will reuse the objects created in the previous chunks. We will then compare the modelled anaerobic duration to observation
#We will now produce S6 to S9 (figure 1B will be produced later, but uses the same lines of code). 
Note. The figures are also reworked in Inkscape to visually improve them
We will use the function "O2_Model". It is a wrapper function written specifically for lake Arendsee and this specific dataset. We have to calculate AD for Loggers, automated profilers and 2 oxic profiles scenario, for multiple years and multiple equations (linear, log-linear and exponential plateau), hence a specific function to avoid copy-pasting the same lines of code 42 times, with slight differences (e.g., which object contains Jz).
Warning. This chunk can take some time depending on your computer RAM.
#I will briefly describe here what the different inputs of the wrapper function mean
b, k and JzMax are the coefficients of the linear, log-linear and exponential plateau equations
Model is either "Livingstone" (linear relationship), "log"(-linear) or "exp"(opential plateau)
alpha is alpha(z)
data.source is always "YSI". We use the YSI profiles as truth and compare the modelled daily profiles to one YSI profile
year is used to select the right YSI year for comparison
```{R AD calculation, warning=FALSE}
ar.alpha <- Bats%>% slice(30:50)
#Anaerobic duration for Loggers - 2017
#Note. Object names are called ".AA" because originally we termed Anaerobic duration "Anoxic Age". I kept the original name in object names to keep a bit of History.
#Each output is a list containing 3 objects. The first is the most important, a matrix of anaerobic duration for each day and each depth. The second is which modelled profile was closest to the YSI profile. The last one is the associated RMSE

#WARNING. Make sure that "ar.alpha" is a matrix in which the second column are alpha values. Otherwise this function won't work properly
for(z in 1:2){
  Loggers.linear.AA.2017 = O2_Model(b = coef(Ar.Living.2017.day.SegLm)[1,1],
                                    k = coef(Ar.Living.2017.day.SegLm)[2,1],
                                    Model = "Livingstone",
                                    alpha = ar.alpha,
                                    data.source = "YSI",
                                    year = 2017,
                                    O2_threshold = O2.threshold[z])
  
  Loggers.log.AA.2017 = O2_Model(b = coef(Ar.Living.2017.day.lin.log)[1],
                                 k = coef(Ar.Living.2017.day.lin.log)[2],
                                 Model = "Log",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2017,
                                 O2_threshold = O2.threshold[z])
  
  Loggers.exp.AA.2017 = O2_Model(b = coef(Ar.Living.2017.day.lin.exp)[2],
                                 k = coef(Ar.Living.2017.day.lin.exp)[3],
                                 JzMax = coef(Ar.Living.2017.day.lin.exp)[1],
                                 Model = "Exp",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2017,
                                 O2_threshold = O2.threshold[z])

#Anaerobic duration for Loggers - 2018
Loggers.linear.AA.2018 = O2_Model(b = coef(Ar.Living.2018.day.SegLm)[1,1],
                                   k = coef(Ar.Living.2018.day.SegLm)[2,1],
                                   Model = "Livingstone",
                                   alpha = ar.alpha,
                                   data.source = "YSI",
                                   year = 2018,
                                   O2_threshold = O2.threshold[z])
  
  Loggers.log.AA.2018 = O2_Model(b = coef(Ar.Living.2018.day.lin.log)[1],
                                 k = coef(Ar.Living.2018.day.lin.log)[2],
                                 Model = "Log",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2018,
                                 O2_threshold = O2.threshold[z])
  
  Loggers.exp.AA.2018 = O2_Model(b = coef(Ar.Living.2018.day.lin.exp)[2],
                                 k = coef(Ar.Living.2018.day.lin.exp)[3],
                                 JzMax = coef(Ar.Living.2018.day.lin.exp)[1],
                                 Model = "Exp",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2018,
                                 O2_threshold = O2.threshold[z])

#Anaerobic duration for Loggers - 2019
Loggers.linear.AA.2019 = O2_Model(b = coef(Ar.Living.2019.day.SegLm)[1,1],
                                   k = coef(Ar.Living.2019.day.SegLm)[2,1],
                                   Model = "Livingstone",
                                   alpha = ar.alpha,
                                   data.source = "YSI",
                                   year = 2019,
                                   O2_threshold = O2.threshold[z])
  
  Loggers.log.AA.2019 = O2_Model(b = coef(Ar.Living.2019.day.lin.log)[1],
                                 k = coef(Ar.Living.2019.day.lin.log)[2],
                                 Model = "Log",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2019,
                                 O2_threshold = O2.threshold[z])
  
  Loggers.exp.AA.2019 = O2_Model(b = coef(Ar.Living.2019.day.lin.exp)[2],
                                 k = coef(Ar.Living.2019.day.lin.exp)[3],
                                 JzMax = coef(Ar.Living.2019.day.lin.exp)[1],
                                 Model = "Exp",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2019,
                                 O2_threshold = O2.threshold[z])

#Anaerobic duration for Loggers - 2020
Loggers.linear.AA.2020 = O2_Model(b = coef(Ar.Living.2020.day.SegLm)[1,1],
                                   k = coef(Ar.Living.2020.day.SegLm)[2,1],
                                   Model = "Livingstone",
                                   alpha = ar.alpha,
                                   data.source = "YSI",
                                   year = 2020,
                                   O2_threshold = O2.threshold[z])
  
  Loggers.log.AA.2020 = O2_Model(b = coef(Ar.Living.2020.day.lin.log)[1],
                                 k = coef(Ar.Living.2020.day.lin.log)[2],
                                 Model = "Log",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2020,
                                 O2_threshold = O2.threshold[z])
  
  Loggers.exp.AA.2020 = O2_Model(b = coef(Ar.Living.2020.day.lin.exp)[2],
                                 k = coef(Ar.Living.2020.day.lin.exp)[3],
                                 JzMax = coef(Ar.Living.2020.day.lin.exp)[1],
                                 Model = "Exp",
                                 alpha = ar.alpha,
                                 data.source = "YSI",
                                 year = 2020,
                                 O2_threshold = O2.threshold[z])



#Anaerobic duration for High temporal Res (YSI) - 2017
YSI.linear.AA.2017 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2017)[1],
                               k = coef(Ar.Jz.alpha.YSI.2017)[2],
                               Model = "Livingstone",
                               alpha = ar.alpha,
                               data.source = "YSI",
                               year = 2017,
                               O2_threshold = O2.threshold[z])
  
  YSI.log.AA.2017 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2017.log)[1],
                             k = coef(Ar.Jz.alpha.YSI.2017.log)[2],
                             Model = "Log",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2017,
                             O2_threshold = O2.threshold[z])
  
  YSI.exp.AA.2017 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2017.exp)[2],
                             k = coef(Ar.Jz.alpha.YSI.2017.exp)[3],
                             JzMax = coef(Ar.Jz.alpha.YSI.2017.exp)[1],
                             Model = "Exp",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2017,
                             O2_threshold = O2.threshold[z])

#Anaerobic duration for High temporal Res (YSI) - 2018
YSI.linear.AA.2018 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2018)[1],
                               k = coef(Ar.Jz.alpha.YSI.2018)[2],
                               Model = "Livingstone",
                               alpha = ar.alpha,
                               data.source = "YSI",
                               year = 2018,
                               O2_threshold = O2.threshold[z])
  
  YSI.log.AA.2018 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2018.log)[1],
                             k = coef(Ar.Jz.alpha.YSI.2018.log)[2],
                             Model = "Log",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2018,
                             O2_threshold = O2.threshold[z])
  
  YSI.exp.AA.2018 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2018.exp)[2],
                             k = coef(Ar.Jz.alpha.YSI.2018.exp)[3],
                             JzMax = coef(Ar.Jz.alpha.YSI.2018.exp)[1],
                             Model = "Exp",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2018,
                             O2_threshold = O2.threshold[z])

#Anaerobic duration for High temporal Res (YSI) - 2019
YSI.linear.AA.2019 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2019)[1],
                               k = coef(Ar.Jz.alpha.YSI.2019)[2],
                               Model = "Livingstone",
                               alpha = ar.alpha,
                               data.source = "YSI",
                               year = 2019,
                               O2_threshold = O2.threshold[z])
  
  YSI.log.AA.2019 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2019.log)[1],
                             k = coef(Ar.Jz.alpha.YSI.2019.log)[2],
                             Model = "Log",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2019,
                             O2_threshold = O2.threshold[z])
  
  YSI.exp.AA.2019 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2019.exp)[2],
                             k = coef(Ar.Jz.alpha.YSI.2019.exp)[3],
                             JzMax = coef(Ar.Jz.alpha.YSI.2019.exp)[1],
                             Model = "Exp",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2019,
                             O2_threshold = O2.threshold[z])

#Anaerobic duration for High temporal Res (YSI) - 2020
YSI.linear.AA.2020 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2020)[1],
                               k = coef(Ar.Jz.alpha.YSI.2020)[2],
                               Model = "Livingstone",
                               alpha = ar.alpha,
                               data.source = "YSI",
                               year = 2020,
                               O2_threshold = O2.threshold[z])
  
  YSI.log.AA.2020 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2020.log)[1],
                             k = coef(Ar.Jz.alpha.YSI.2020.log)[2],
                             Model = "Log",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2020,
                             O2_threshold = O2.threshold[z])
  
  YSI.exp.AA.2020 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2020.exp)[2],
                             k = coef(Ar.Jz.alpha.YSI.2020.exp)[3],
                             JzMax = coef(Ar.Jz.alpha.YSI.2020.exp)[1],
                             Model = "Exp",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2020,
                             O2_threshold = O2.threshold[z])
  
#Anaerobic duration for High temporal Res (YSI) - 2021
YSI.linear.AA.2021 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2021)[1],
                               k = coef(Ar.Jz.alpha.YSI.2021)[2],
                               Model = "Livingstone",
                               alpha = ar.alpha,
                               data.source = "YSI",
                               year = 2021,
                               O2_threshold = O2.threshold[z])
  
  YSI.log.AA.2021 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2021.log)[1],
                             k = coef(Ar.Jz.alpha.YSI.2021.log)[2],
                             Model = "Log",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2021,
                             O2_threshold = O2.threshold[z])
  
  YSI.exp.AA.2021 = O2_Model(b = coef(Ar.Jz.alpha.YSI.2021.exp)[2],
                             k = coef(Ar.Jz.alpha.YSI.2021.exp)[3],
                             JzMax = coef(Ar.Jz.alpha.YSI.2021.exp)[1],
                             Model = "Exp",
                             alpha = ar.alpha,
                             data.source = "YSI",
                             year = 2021,
                             O2_threshold = O2.threshold[z])


#Anaerobic duration for Best case (YSI) - 2017
  #Note. I used "Best" here because I ran out of imagination to name my objects and wanted some positivity deep in the script
Best.linear.AA.2017 = O2_Model(b = coef(Ar.Jz.lm.best.list[[1]][[1]])[1],
                                k = coef(Ar.Jz.lm.best.list[[1]][[1]])[2],
                                Model = "Livingstone",
                                alpha = ar.alpha,
                                data.source = "YSI",
                                year = 2017,
                                O2_threshold = O2.threshold[z])
  
  Best.log.AA.2017 = O2_Model(b = coef(Ar.Jz.lm.best.list[[1]][[2]])[1],
                              k = coef(Ar.Jz.lm.best.list[[1]][[2]])[2],
                              Model = "Log",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2017,
                              O2_threshold = O2.threshold[z])
  
  Best.exp.AA.2017 = O2_Model(b = coef(Ar.Jz.lm.best.list[[1]][[3]])[2],
                              k = coef(Ar.Jz.lm.best.list[[1]][[3]])[3],
                              JzMax = coef(Ar.Jz.lm.best.list[[1]][[3]])[1],
                              Model = "Exp",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2017,
                              O2_threshold = O2.threshold[z])

#Anaerobic duration for Best case (YSI) - 2018
Best.linear.AA.2018 = O2_Model(b = coef(Ar.Jz.lm.best.list[[2]][[1]])[1],
                                k = coef(Ar.Jz.lm.best.list[[2]][[1]])[2],
                                Model = "Livingstone",
                                alpha = ar.alpha,
                                data.source = "YSI",
                                year = 2018,
                                O2_threshold = O2.threshold[z])
  
  Best.log.AA.2018 = O2_Model(b = coef(Ar.Jz.lm.best.list[[2]][[2]])[1],
                              k = coef(Ar.Jz.lm.best.list[[2]][[2]])[2],
                              Model = "Log",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2018,
                              O2_threshold = O2.threshold[z])
  
  Best.exp.AA.2018 = O2_Model(b = coef(Ar.Jz.lm.best.list[[2]][[3]])[2],
                              k = coef(Ar.Jz.lm.best.list[[2]][[3]])[3],
                              JzMax = coef(Ar.Jz.lm.best.list[[2]][[3]])[1],
                              Model = "Exp",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2018,
                              O2_threshold = O2.threshold[z])

#Anaerobic duration for Best case (YSI) - 2019
Best.linear.AA.2019 = O2_Model(b = coef(Ar.Jz.lm.best.list[[3]][[1]])[1],
                                k = coef(Ar.Jz.lm.best.list[[3]][[1]])[2],
                                Model = "Livingstone",
                                alpha = ar.alpha,
                                data.source = "YSI",
                                year = 2019,
                               O2_threshold = O2.threshold[z])
  
  Best.log.AA.2019 = O2_Model(b = coef(Ar.Jz.lm.best.list[[3]][[2]])[1],
                              k = coef(Ar.Jz.lm.best.list[[3]][[2]])[2],
                              Model = "Log",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2019,
                              O2_threshold = O2.threshold[z])
  
  Best.exp.AA.2019 = O2_Model(b = coef(Ar.Jz.lm.best.list[[3]][[3]])[2],
                              k = coef(Ar.Jz.lm.best.list[[3]][[3]])[3],
                              JzMax = coef(Ar.Jz.lm.best.list[[3]][[3]])[1],
                              Model = "Exp",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2019,
                              O2_threshold = O2.threshold[z])

#Anaerobic duration for Best case (YSI) - 2020
Best.linear.AA.2020 = O2_Model(b = coef(Ar.Jz.lm.best.list[[4]][[1]])[1],
                                k = coef(Ar.Jz.lm.best.list[[4]][[1]])[2],
                                Model = "Livingstone",
                                alpha = ar.alpha,
                                data.source = "YSI",
                                year = 2020,
                                O2_threshold = O2.threshold[z])
  
  Best.log.AA.2020 = O2_Model(b = coef(Ar.Jz.lm.best.list[[4]][[2]])[1],
                              k = coef(Ar.Jz.lm.best.list[[4]][[2]])[2],
                              Model = "Log",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2020,
                              O2_threshold = O2.threshold[z])
  
  Best.exp.AA.2020 = O2_Model(b = coef(Ar.Jz.lm.best.list[[4]][[3]])[2],
                              k = coef(Ar.Jz.lm.best.list[[4]][[3]])[3],
                              JzMax = coef(Ar.Jz.lm.best.list[[4]][[3]])[1],
                              Model = "Exp",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2020,
                              O2_threshold = O2.threshold[z])

#Anaerobic duration for Best case (YSI) - 2021
Best.linear.AA.2021 = O2_Model(b = coef(Ar.Jz.lm.best.list[[5]][[1]])[1],
                                k = coef(Ar.Jz.lm.best.list[[5]][[1]])[2],
                                Model = "Livingstone",
                                alpha = ar.alpha,
                                data.source = "YSI",
                                year = 2021,
                                O2_threshold = O2.threshold[z])
  
  Best.log.AA.2021 = O2_Model(b = coef(Ar.Jz.lm.best.list[[5]][[2]])[1],
                              k = coef(Ar.Jz.lm.best.list[[5]][[2]])[2],
                              Model = "Log",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2021,
                              O2_threshold = O2.threshold[z])
  
  Best.exp.AA.2021 = O2_Model(b = coef(Ar.Jz.lm.best.list[[5]][[3]])[2],
                              k = coef(Ar.Jz.lm.best.list[[5]][[3]])[3],
                              JzMax = coef(Ar.Jz.lm.best.list[[5]][[3]])[1],
                              Model = "Exp",
                              alpha = ar.alpha,
                              data.source = "YSI",
                              year = 2021,
                              O2_threshold = O2.threshold[z])

############################SECOND SECTION##############################
#We first need to calculate anaerobic duration based on observations only. Each year will be stored as a different item in the list
output=list()

for(k in 1:length(YSI.yrs)){
  #Read the file
  Ar.YSI = read.csv(paste0("./",YSI.list[k]))
  
  #Remove NA at bottom of profiles with minimum value of the profile
  for(j in 1:ncol(Ar.YSI))
  {
    if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
  }
  
  #Transform to long format
  Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
  #Split Date into Year, months and day
  Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")
  
  Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))
  
  Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                               mm = substring(Date, 6, 7),
                               dd = substring(Date, 9, 10))
  
  #Keep only values below 30 meters to match other methods
  Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]
  
  #Order the dataframe
  Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]
  
  temp.output = matrix(nrow = length(unique(Ar.long.deep$Depth_m)), ncol = length(unique(Ar.long.deep$DOY)))
  colnames(temp.output) = unique(Ar.long.deep$DOY)
  rownames(temp.output) = unique(Ar.long.deep$Depth_m)
  
  for(i in 1:length(unique(Ar.long.deep$Depth_m))){
    temp = vector(length = dim(Ar.long.deep[Ar.long.deep$Depth_m == unique(Ar.long.deep$Depth_m)[i],"DOY"])[1])
    temp = Ar.long.deep[Ar.long.deep$Depth_m == unique(Ar.long.deep$Depth_m)[i],] %>% group_by(DOY) %>% summarise(DO = mean(DO_mgL, na.rm=T))
    temp = as.data.frame(temp)
    temp.output[i,] = temp[,2]
  }
  output[[k]] = temp.output
}

#Home function to interpolate data on missing days. This function is not super flexible, but works great with Arendsee data. I do not recommend using it as is but rather to look at the code behind it and make sure it works with other datasets
interpolated.2017 = interpolation_fct(output[[1]])
interpolated.2018 = interpolation_fct(output[[2]][,-283])
interpolated.2019 = interpolation_fct(output[[3]])
interpolated.2020 = interpolation_fct(output[[4]])
interpolated.2021 = interpolation_fct(output[[5]])

#Here we select our oxygen threshold for anaerobic conditions
O2_threshold =  O2.threshold[z]

#Transform oxygen values below the oxygen threshold into "1" (yes/no values)
interpolated.2017.bin <- ifelse(interpolated.2017 <= O2_threshold, 1, 0)
interpolated.2018.bin <- ifelse(interpolated.2018 <= O2_threshold, 1, 0)
interpolated.2019.bin <- ifelse(interpolated.2019 <= O2_threshold, 1, 0)
interpolated.2020.bin <- ifelse(interpolated.2020 <= O2_threshold, 1, 0)
interpolated.2021.bin <- ifelse(interpolated.2021 <= O2_threshold, 1, 0)

#We transform 1 and 0 into TURE/FALSE
YSI.aa.2017 <- t(apply(interpolated.2017.bin, 1, function(X) X==1))
YSI.aa.2018 <- t(apply(interpolated.2018.bin, 1, function(X) X==1))
YSI.aa.2019 <- t(apply(interpolated.2019.bin, 1, function(X) X==1))
YSI.aa.2020 <- t(apply(interpolated.2020.bin, 1, function(X) X==1))
YSI.aa.2021 <- t(apply(interpolated.2021.bin, 1, function(X) X==1))

#Create a short function to find first day of continuous anaerobic conditions for each year and each depth
First.aa.doy.fct <- function(x){
  output = vector(length=nrow(x))
  for(i in 1:nrow(x)){
    output[i] = min(which(x[i,]==TRUE)) + as.numeric(colnames(x)[1])
  }
  return(output)
}

#Observed profiles
#Note. I removed the first few months to make sure I did not capture winter anaerobic conditions
YSI.aa.2017.cut = First.aa.doy.fct(YSI.aa.2017[,-c(1:120)])
YSI.aa.2018.cut = First.aa.doy.fct(YSI.aa.2018[,-c(1:120)])
YSI.aa.2019.cut = First.aa.doy.fct(YSI.aa.2019[,-c(1:120)])
YSI.aa.2020.cut = First.aa.doy.fct(YSI.aa.2020[,-c(1:120)])
YSI.aa.2021.cut = First.aa.doy.fct(YSI.aa.2021[,-c(1:120)])

#Repeat the previous steps with all modelled profiles
#Loggers-2017
{Loggers.linear.AA.2017.bin = t(apply(Loggers.linear.AA.2017[[1]][,-1], 1, function(X) X==1))
  Loggers.log.AA.2017.bin = t(apply(Loggers.log.AA.2017[[1]][,-1], 1, function(X) X==1))
  Loggers.exp.AA.2017.bin = t(apply(Loggers.exp.AA.2017[[1]][,-1], 1, function(X) X==1))
  
  Loggers.linear.AA.2017.FirstaaDOY = First.aa.doy.fct(Loggers.linear.AA.2017.bin)
  Loggers.log.AA.2017.FirstaaDOY = First.aa.doy.fct(Loggers.log.AA.2017.bin)
  Loggers.exp.AA.2017.FirstaaDOY = First.aa.doy.fct(Loggers.exp.AA.2017.bin)}

#Loggers-2018
{Loggers.linear.AA.2018.bin = t(apply(Loggers.linear.AA.2018[[1]][,-1], 1, function(X) X==1))
  Loggers.log.AA.2018.bin = t(apply(Loggers.log.AA.2018[[1]][,-1], 1, function(X) X==1))
  Loggers.exp.AA.2018.bin = t(apply(Loggers.exp.AA.2018[[1]][,-1], 1, function(X) X==1))
  
  Loggers.linear.AA.2018.FirstaaDOY = First.aa.doy.fct(Loggers.linear.AA.2018.bin)
  Loggers.log.AA.2018.FirstaaDOY = First.aa.doy.fct(Loggers.log.AA.2018.bin)
  Loggers.exp.AA.2018.FirstaaDOY = First.aa.doy.fct(Loggers.exp.AA.2018.bin)}

#Loggers-2019
{Loggers.linear.AA.2019.bin = t(apply(Loggers.linear.AA.2019[[1]][,-1], 1, function(X) X==1))
  Loggers.log.AA.2019.bin = t(apply(Loggers.log.AA.2019[[1]][,-1], 1, function(X) X==1))
  Loggers.exp.AA.2019.bin = t(apply(Loggers.exp.AA.2019[[1]][,-1], 1, function(X) X==1))
  
  Loggers.linear.AA.2019.FirstaaDOY = First.aa.doy.fct(Loggers.linear.AA.2019.bin)
  Loggers.log.AA.2019.FirstaaDOY = First.aa.doy.fct(Loggers.log.AA.2019.bin)
  Loggers.exp.AA.2019.FirstaaDOY = First.aa.doy.fct(Loggers.exp.AA.2019.bin)}

#Loggers-2020
{Loggers.linear.AA.2020.bin = t(apply(Loggers.linear.AA.2020[[1]][,-1], 1, function(X) X==1))
  Loggers.log.AA.2020.bin = t(apply(Loggers.log.AA.2020[[1]][,-1], 1, function(X) X==1))
  Loggers.exp.AA.2020.bin = t(apply(Loggers.exp.AA.2020[[1]][,-1], 1, function(X) X==1))
  
  Loggers.linear.AA.2020.FirstaaDOY = First.aa.doy.fct(Loggers.linear.AA.2020.bin)
  Loggers.log.AA.2020.FirstaaDOY = First.aa.doy.fct(Loggers.log.AA.2020.bin)
  Loggers.exp.AA.2020.FirstaaDOY = First.aa.doy.fct(Loggers.exp.AA.2020.bin)}

#YSI - 2017
{YSI.linear.AA.2017.bin = t(apply(YSI.linear.AA.2017[[1]][,-1], 1, function(X) X==1))
  YSI.log.AA.2017.bin = t(apply(YSI.log.AA.2017[[1]][,-1], 1, function(X) X==1))
  YSI.exp.AA.2017.bin = t(apply(YSI.exp.AA.2017[[1]][,-1], 1, function(X) X==1))
  
  YSI.linear.AA.2017.FirstaaDOY = First.aa.doy.fct(YSI.linear.AA.2017.bin)
  YSI.log.AA.2017.FirstaaDOY = First.aa.doy.fct(YSI.log.AA.2017.bin)
  YSI.exp.AA.2017.FirstaaDOY = First.aa.doy.fct(YSI.exp.AA.2017.bin)}

#YSI - 2018
{YSI.linear.AA.2018.bin = t(apply(YSI.linear.AA.2018[[1]][,-1], 1, function(X) X==1))
  YSI.log.AA.2018.bin = t(apply(YSI.log.AA.2018[[1]][,-1], 1, function(X) X==1))
  YSI.exp.AA.2018.bin = t(apply(YSI.exp.AA.2018[[1]][,-1], 1, function(X) X==1))
  
  YSI.linear.AA.2018.FirstaaDOY = First.aa.doy.fct(YSI.linear.AA.2018.bin)
  YSI.log.AA.2018.FirstaaDOY = First.aa.doy.fct(YSI.log.AA.2018.bin)
  YSI.exp.AA.2018.FirstaaDOY = First.aa.doy.fct(YSI.exp.AA.2018.bin)}

#YSI - 2019
{YSI.linear.AA.2019.bin = t(apply(YSI.linear.AA.2019[[1]][,-1], 1, function(X) X==1))
  YSI.log.AA.2019.bin = t(apply(YSI.log.AA.2019[[1]][,-1], 1, function(X) X==1))
  YSI.exp.AA.2019.bin = t(apply(YSI.exp.AA.2019[[1]][,-1], 1, function(X) X==1))
  
  YSI.linear.AA.2019.FirstaaDOY = First.aa.doy.fct(YSI.linear.AA.2019.bin)
  YSI.log.AA.2019.FirstaaDOY = First.aa.doy.fct(YSI.log.AA.2019.bin)
  YSI.exp.AA.2019.FirstaaDOY = First.aa.doy.fct(YSI.exp.AA.2019.bin)}

#YSI - 2020
{YSI.linear.AA.2020.bin = t(apply(YSI.linear.AA.2020[[1]][,-1], 1, function(X) X==1))
  YSI.log.AA.2020.bin = t(apply(YSI.log.AA.2020[[1]][,-1], 1, function(X) X==1))
  YSI.exp.AA.2020.bin = t(apply(YSI.exp.AA.2020[[1]][,-1], 1, function(X) X==1))
  
  YSI.linear.AA.2020.FirstaaDOY = First.aa.doy.fct(YSI.linear.AA.2020.bin)
  YSI.log.AA.2020.FirstaaDOY = First.aa.doy.fct(YSI.log.AA.2020.bin)
  YSI.exp.AA.2020.FirstaaDOY = First.aa.doy.fct(YSI.exp.AA.2020.bin)}

#YSI - 2021
{YSI.linear.AA.2021.bin = t(apply(YSI.linear.AA.2021[[1]][,-1], 1, function(X) X==1))
  YSI.log.AA.2021.bin = t(apply(YSI.log.AA.2021[[1]][,-1], 1, function(X) X==1))
  YSI.exp.AA.2021.bin = t(apply(YSI.exp.AA.2021[[1]][,-1], 1, function(X) X==1))
  
  YSI.linear.AA.2021.FirstaaDOY = First.aa.doy.fct(YSI.linear.AA.2021.bin)
  YSI.log.AA.2021.FirstaaDOY = First.aa.doy.fct(YSI.log.AA.2021.bin)
  YSI.exp.AA.2021.FirstaaDOY = First.aa.doy.fct(YSI.exp.AA.2021.bin)}

#Best-2017
{Best.linear.AA.2017.bin = t(apply(Best.linear.AA.2017[[1]][,-1], 1, function(X) X==1))
  Best.log.AA.2017.bin = t(apply(Best.log.AA.2017[[1]][,-1], 1, function(X) X==1))
  Best.exp.AA.2017.bin = t(apply(Best.exp.AA.2017[[1]][,-1], 1, function(X) X==1))
  
  Best.linear.AA.2017.FirstaaDOY = First.aa.doy.fct(Best.linear.AA.2017.bin)
  Best.log.AA.2017.FirstaaDOY = First.aa.doy.fct(Best.log.AA.2017.bin)
  Best.exp.AA.2017.FirstaaDOY = First.aa.doy.fct(Best.exp.AA.2017.bin)}

#Best-2018
{Best.linear.AA.2018.bin = t(apply(Best.linear.AA.2018[[1]][,-1], 1, function(X) X==1))
  Best.log.AA.2018.bin = t(apply(Best.log.AA.2018[[1]][,-1], 1, function(X) X==1))
  Best.exp.AA.2018.bin = t(apply(Best.exp.AA.2018[[1]][,-1], 1, function(X) X==1))
  
  Best.linear.AA.2018.FirstaaDOY = First.aa.doy.fct(Best.linear.AA.2018.bin)
  Best.log.AA.2018.FirstaaDOY = First.aa.doy.fct(Best.log.AA.2018.bin)
  Best.exp.AA.2018.FirstaaDOY = First.aa.doy.fct(Best.exp.AA.2018.bin)}

#Best-2019
{Best.linear.AA.2019.bin = t(apply(Best.linear.AA.2019[[1]][,-1], 1, function(X) X==1))
  Best.log.AA.2019.bin = t(apply(Best.log.AA.2019[[1]][,-1], 1, function(X) X==1))
  Best.exp.AA.2019.bin = t(apply(Best.exp.AA.2019[[1]][,-1], 1, function(X) X==1))
  
  Best.linear.AA.2019.FirstaaDOY = First.aa.doy.fct(Best.linear.AA.2019.bin)
  Best.log.AA.2019.FirstaaDOY = First.aa.doy.fct(Best.log.AA.2019.bin)
  Best.exp.AA.2019.FirstaaDOY = First.aa.doy.fct(Best.exp.AA.2019.bin)}

#Best-2020
{Best.linear.AA.2020.bin = t(apply(Best.linear.AA.2020[[1]][,-1], 1, function(X) X==1))
  Best.log.AA.2020.bin = t(apply(Best.log.AA.2020[[1]][,-1], 1, function(X) X==1))
  Best.exp.AA.2020.bin = t(apply(Best.exp.AA.2020[[1]][,-1], 1, function(X) X==1))
  
  Best.linear.AA.2020.FirstaaDOY = First.aa.doy.fct(Best.linear.AA.2020.bin)
  Best.log.AA.2020.FirstaaDOY = First.aa.doy.fct(Best.log.AA.2020.bin)
  Best.exp.AA.2020.FirstaaDOY = First.aa.doy.fct(Best.exp.AA.2020.bin)}

#Best-2021
{Best.linear.AA.2021.bin = t(apply(Best.linear.AA.2021[[1]][,-1], 1, function(X) X==1))
  Best.log.AA.2021.bin = t(apply(Best.log.AA.2021[[1]][,-1], 1, function(X) X==1))
  Best.exp.AA.2021.bin = t(apply(Best.exp.AA.2021[[1]][,-1], 1, function(X) X==1))
  
  Best.linear.AA.2021.FirstaaDOY = First.aa.doy.fct(Best.linear.AA.2021.bin)
  Best.log.AA.2021.FirstaaDOY = First.aa.doy.fct(Best.log.AA.2021.bin)
  Best.exp.AA.2021.FirstaaDOY = First.aa.doy.fct(Best.exp.AA.2021.bin)
}


#Combine into multiple and 1 big table
{First.2017.all = rbind(Observed2017 = YSI.aa.2017.cut[-19],
                        Loggers.linear.AA.2017.FirstaaDOY[-c(19:21)],
                        Loggers.log.AA.2017.FirstaaDOY[-c(19:21)],
                        Loggers.exp.AA.2017.FirstaaDOY[-c(19:21)],
                        YSI.linear.AA.2017.FirstaaDOY[-c(19:21)],
                        YSI.log.AA.2017.FirstaaDOY[-c(19:21)],
                        YSI.exp.AA.2017.FirstaaDOY[-c(19:21)],
                        Best.linear.AA.2017.FirstaaDOY[-c(19:21)],
                        Best.log.AA.2017.FirstaaDOY[-c(19:21)],
                        Best.exp.AA.2017.FirstaaDOY[-c(19:21)])
  First.2018.all = rbind(Observed2018 = YSI.aa.2018.cut,
                         Loggers.linear.AA.2018.FirstaaDOY[-c(19:21)],
                         Loggers.log.AA.2018.FirstaaDOY[-c(19:21)],
                         Loggers.exp.AA.2018.FirstaaDOY[-c(19:21)],
                         YSI.linear.AA.2018.FirstaaDOY[-c(19:21)],
                         YSI.log.AA.2018.FirstaaDOY[-c(19:21)],
                         YSI.exp.AA.2018.FirstaaDOY[-c(19:21)],
                         Best.linear.AA.2018.FirstaaDOY[-c(19:21)],
                         Best.log.AA.2018.FirstaaDOY[-c(19:21)],
                         Best.exp.AA.2018.FirstaaDOY[-c(19:21)])
  First.2019.all = rbind(Observed2019 = YSI.aa.2019.cut,
                         Loggers.linear.AA.2019.FirstaaDOY[-c(19:21)],
                         Loggers.log.AA.2019.FirstaaDOY[-c(19:21)],
                         Loggers.exp.AA.2019.FirstaaDOY[-c(19:21)],
                         YSI.linear.AA.2019.FirstaaDOY[-c(19:21)],
                         YSI.log.AA.2019.FirstaaDOY[-c(19:21)],
                         YSI.exp.AA.2019.FirstaaDOY[-c(19:21)],
                         Best.linear.AA.2019.FirstaaDOY[-c(19:21)],
                         Best.log.AA.2019.FirstaaDOY[-c(19:21)],
                         Best.exp.AA.2019.FirstaaDOY[-c(19:21)])
  First.2020.all = rbind(Observed2020 = YSI.aa.2020.cut,
                         Loggers.linear.AA.2020.FirstaaDOY[-c(19:21)],
                         Loggers.log.AA.2020.FirstaaDOY[-c(19:21)],
                         Loggers.exp.AA.2020.FirstaaDOY[-c(19:21)],
                         YSI.linear.AA.2020.FirstaaDOY[-c(19:21)],
                         YSI.log.AA.2020.FirstaaDOY[-c(19:21)],
                         YSI.exp.AA.2020.FirstaaDOY[-c(19:21)],
                         Best.linear.AA.2020.FirstaaDOY[-c(19:21)],
                         Best.log.AA.2020.FirstaaDOY[-c(19:21)],
                         Best.exp.AA.2020.FirstaaDOY[-c(19:21)])
  First.2021.all = rbind(Observed2021 = YSI.aa.2021.cut[-c(19:21)],
                         YSI.linear.AA.2021.FirstaaDOY[-c(19:21)],
                         YSI.log.AA.2021.FirstaaDOY[-c(19:21)],
                         YSI.exp.AA.2021.FirstaaDOY[-c(19:21)],
                         Best.linear.AA.2021.FirstaaDOY[-c(19:21)],
                         Best.log.AA.2021.FirstaaDOY[-c(19:21)],
                         Best.exp.AA.2021.FirstaaDOY[-c(19:21)])
  
  
  First.DOY.all = rbind(First.2017.all, First.2018.all, First.2019.all, First.2020.all,First.2021.all)
}
#Export a csv file if wanted
# write.csv(file = "./Output/First.DOY.all.csv", x = First.DOY.all, fileEncoding = "UTF-8")

###################################PRINTING FIGURES ONLY WHEN z=1##############################################
if(z==1){
#2017
{
  pdf("./Output/FigS6.pdf", height=8, width=8)
  par(mfrow=c(3,3))
  par(mar=c(4,5,1,0.5)+.1)
  for(i in 2:nrow(First.2017.all)){
    if(i == 2){
      plot(First.2017.all[1,] ~ First.2017.all[i,],
           xlab = "",
           xlim = c(180,350),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1, cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Linear model"))
      legend("topleft", lty = 1, legend = ("1:1 line"))
      abline(0,1)
    }
    if(i == 3 | i ==4){
      plot(First.2017.all[1,] ~ First.2017.all[i,],
           xlab = "",
           xlim = c(180,350),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Log-linear model", "Exponential plateau model")[i-2])
      abline(0,1)
    }
    if(i == 5){
      plot(First.2017.all[1,] ~ First.2017.all[i,],
           xlab = "",
           xlim = c(180,350),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 6 | i== 7){
      plot(First.2017.all[1,] ~ First.2017.all[i,],
           xlab = "",
           xlim = c(180,350),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 8){
      plot(First.2017.all[1,] ~ First.2017.all[i,],
           xlim = c(180,350),
           xlab = "First day of anoxia (modelled)",
           ylab = "First day of anoxia (observed)",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    if(i == 9| i == 10){
      plot(First.2017.all[1,] ~ First.2017.all[i,],
           xlim = c(180,350),
           xlab = "First day of anoxia (modelled)",
           ylab = "",
           yaxt = "n",
           las = 1,cex = 1.8, cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
  }
  dev.off()
}

#2018
{
  pdf("./Output/FigS7.pdf", height=8, width=8)
  par(mfrow=c(3,3))
  par(mar=c(4,5,1,0.5)+.1)
  for(i in 2:nrow(First.2018.all)){
    if(i == 2){
      plot(First.2018.all[1,] ~ First.2018.all[i,],
           xlab = "",
           xlim = c(190,350),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1, cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Linear model"))
      legend("topleft", lty = 1, legend = ("1:1 line"))
      abline(0,1)
    }
    if(i == 3 | i ==4){
      plot(First.2018.all[1,] ~ First.2018.all[i,],
           xlab = "",
           xlim = c(190,350),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Log-linear model", "Exponential plateau model")[i-2])
      abline(0,1)
    }
    if(i == 5){
      plot(First.2018.all[1,] ~ First.2018.all[i,],
           xlab = "",
           xlim = c(190,350),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 6 | i== 7){
      plot(First.2018.all[1,] ~ First.2018.all[i,],
           xlab = "",
           xlim = c(190,350),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 8){
      plot(First.2018.all[1,] ~ First.2018.all[i,],
           xlim = c(190,350),
           xlab = "First day of anoxia (modelled)",
           ylab = "First day of anoxia (observed)",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    if(i == 9| i == 10){
      plot(First.2018.all[1,] ~ First.2018.all[i,],
           xlim = c(190,350),
           xlab = "First day of anoxia (modelled)",
           ylab = "",
           yaxt = "n",
           las = 1,cex = 1.8, cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
  }
  dev.off()
}


#2020
{
  pdf("./Output/FigS8.pdf", height=8, width=8)
  par(mfrow=c(3,3))
  par(mar=c(4,5,1,0.5)+.1)
  for(i in 2:nrow(First.2020.all)){
    if(i == 2){
      plot(First.2020.all[1,] ~ First.2020.all[i,],
           xlab = "",
           xlim = c(170,330),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1, cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Linear model"))
      legend("topleft", lty = 1, legend = ("1:1 line"))
      abline(0,1)
    }
    if(i == 3 | i ==4){
      plot(First.2020.all[1,] ~ First.2020.all[i,],
           xlab = "",
           xlim = c(170,330),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Log-linear model", "Exponential plateau model")[i-2])
      abline(0,1)
    }
    if(i == 5){
      plot(First.2020.all[1,] ~ First.2020.all[i,],
           xlab = "",
           xlim = c(170,330),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 6 | i== 7){
      plot(First.2020.all[1,] ~ First.2020.all[i,],
           xlab = "",
           xlim = c(170,330),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 8){
      plot(First.2020.all[1,] ~ First.2020.all[i,],
           xlim = c(170,330),
           xlab = "First day of anoxia (modelled)",
           ylab = "First day of anoxia (observed)",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    if(i == 9| i == 10){
      plot(First.2020.all[1,] ~ First.2020.all[i,],
           xlim = c(170,330),
           xlab = "First day of anoxia (modelled)",
           ylab = "",
           yaxt = "n",
           las = 1,cex = 1.8, cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
  }
  dev.off()
}

#2021
{
  pdf("./Output/FigS9.pdf", height=5.6, width=8)
  par(mfrow=c(2,3))
  par(mar=c(4,5,1,0.5)+.1)
  for(i in 2:nrow(First.2021.all)){
    if(i == 2){
      plot(First.2021.all[1,] ~ First.2021.all[i,],
           xlab = "",
           xlim = c(160,330),
           ylab = "First day of anoxia (observed)",
           xaxt = "n",
           las = 1, cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Linear model"))
      legend("topleft", lty = 1, legend = ("1:1 line"))
      abline(0,1)
    }
    if(i == 3 | i ==4){
      plot(First.2021.all[1,] ~ First.2021.all[i,],
           xlab = "",
           xlim = c(160,330),
           ylab = "",
           xaxt="n",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = c("Log-linear model", "Exponential plateau model")[i-2])
      abline(0,1)
    }
    if(i == 5){
      plot(First.2021.all[1,] ~ First.2021.all[i,],
           xlab = "First day of anoxia (modelled)",
           xlim = c(160,330),
           ylab = "First day of anoxia (observed)",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
    if(i == 6 | i== 7){
      plot(First.2021.all[1,] ~ First.2021.all[i,],
           xlab = "First day of anoxia (modelled)",
           xlim = c(160,330),
           ylab = "",
           yaxt="n",
           las = 1,cex = 1.8,cex.axis=1.3, cex.lab = 1.4,
           main = "")
      abline(0,1)
    }
    
  }
  dev.off()
}
}
}
```

#Now is the time to produce figure 1. This chunk of code is 300 lines, used to produced a multipanel (18 panels) figure
```{R Figure 1}
{
  pdf("./Output/Fig.1.pdf", width = 9, height = 4.5, useDingbats = F)
  par(mfrow=c(3,6))
  par(mar=c(4,5,1,0.5)+.1)
  ############################################Loggers########################################
  plot(Ar.Jz.day.2017.lin ~ Jz.mat.day[, 5], las =1,
       # xlab = expression(alpha(z)),
       xlab = "",
       xaxt = "n",
       xlim = c(0,0.25),
       ylab = "",#expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       ylim = c(0.04, 0.25),cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       main = "Linear", col ="blue",pch=16)
  text(x = 0.03, y = 0.24, labels = "A1")
  # abline(lm(Ar.Jz.2017 ~ Jz.mat[, 5]), col ="blue", lty =2) 
  #NOTE. These are greyed out because they did not respect our conditions to draw the regression (R2 > 0.97)
  
  points(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5], col ="chartreuse3",pch=16,cex = 1.4)
  abline(lm(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5]), col ="chartreuse3", lty = 2)
  points(Ar.Jz.day.2019.lin ~ Jz.mat.day[, 5], col = "red",pch=16,cex = 1.4)
  # abline(lm(Ar.Jz.2019 ~ Jz.mat[, 5]), col ="red")
  
  points(Ar.Jz.day.2020.lin ~ Jz.mat.day[,5], col="pink", pch = 16,cex = 1.4)
  # abline(lm(Ar.Jz.2020 ~ Jz.mat[, 5]), col ="pink")
  
  #Log-linear
  plot(Ar.Jz.day.2017.lin ~ Jz.mat.day[, 5], las =1,
       # xlab = expression(alpha(z)),
       # ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       xlab = "",
       ylab = "",
       xaxt = "n",
       yaxt = "n",
       xlim = c(0,0.25),
       ylim = c(0.04, 0.25),cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       main = "Log-linear", col ="blue",pch=16)
  text(x = 0.03, y = 0.24, labels = "A2")
  curve(coef(Ar.Living.2017.day.lin.log )[1]*log10(x*coef(Ar.Living.2017.day.lin.log)[2]),
        add=T, col ="blue", lty = 2)
  points(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5], col ="chartreuse3",pch=16,cex = 1.4)
  curve(coef(Ar.Living.2018.day.lin.log )[1]*log10(x*coef(Ar.Living.2018.day.lin.log )[2]),
        add=T, col ="chartreuse3", lty = 2)
  points(Ar.Jz.day.2019.lin ~ Jz.mat.day[, 5], col = "red",pch=16,cex = 1.4)
  curve(coef(Ar.Living.2019.day.lin.log)[1]*log10(x*coef(Ar.Living.2019.day.lin.log)[2]),
        add=T, col ="red")
  
  points(Ar.Jz.day.2020.lin ~ Jz.mat.day[, 5], col="pink", pch = 16,cex = 1.4)
  curve(coef(Ar.Living.2020.day.lin.log)[1]*log10(x*coef(Ar.Living.2020.day.lin.log)[2]),
        add=T, col ="pink")
  
  #Exponential plateau
  plot(Ar.Jz.day.2017.lin ~ Jz.mat.day[, 5], las =1,
       # xlab = expression(alpha(z)),
       # ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       xlab = "",
       ylab = "",
       xaxt = "n",
       yaxt = "n",
       xlim = c(0,0.25),
       ylim = c(0.04, 0.25),cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       main = "Exponential plateau", col ="blue",pch=16)
  text(x = 0.03, y = 0.24, labels = "A3")
  curve(coef(Ar.Living.2017.day.lin.exp )[1] - (coef(Ar.Living.2017.day.lin.exp )[1] - coef(Ar.Living.2017.day.lin.exp )[2])*
          exp(-coef(Ar.Living.2017.day.lin.exp )[3]*x),
        add=T, col ="blue")
  points(Ar.Jz.day.2018.lin ~ Jz.mat.day[, 5], col ="chartreuse3",pch=16,cex = 1.4)
  curve(coef(Ar.Living.2018.day.lin.exp)[1] - (coef(Ar.Living.2018.day.lin.exp)[1] - coef(Ar.Living.2018.day.lin.exp)[2])*
          exp(-coef(Ar.Living.2018.day.lin.exp)[3]*x),
        add=T, col ="chartreuse3", lty = 2)
  points(Ar.Jz.day.2019.lin ~ Jz.mat.day[, 5], col = "red",pch=16,cex = 1.4)
  curve(coef(Ar.Living.2019.day.lin.exp)[1] - (coef(Ar.Living.2019.day.lin.exp)[1] - coef(Ar.Living.2019.day.lin.exp)[2])*
          exp(-coef(Ar.Living.2019.day.lin.exp)[3]*x),
        add=T, col ="red")
  
  points(Ar.Jz.day.2020.lin ~ Jz.mat.day[, 5], col="pink", pch = 16,cex = 1.4)
  curve(coef(Ar.Living.2020.day.lin.exp)[1] - (coef(Ar.Living.2020.day.lin.exp)[1] - coef(Ar.Living.2020.day.lin.exp)[2])*
          exp(-coef(Ar.Living.2020.day.lin.exp)[3]*x),
        add=T, col ="pink")
  
  #First day of anaerobic conditions, loggers
  plot(First.2019.all[1,] ~ First.2019.all[2,],
       xlab = "",
       xlim = c(190,350),
       ylab = "",
       xaxt = "n",
       las = 1, cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = c("Linear model"),
       col = "red")
  text(x=210, y=320, labels="B1")
  # legend("topleft", lty = 1, legend = ("1:1 line"))
  abline(0,1)
  
  plot(First.2019.all[1,] ~ First.2019.all[3,],
       xlab = "",
       xlim = c(190,350),
       ylab = "",
       xaxt="n",
       yaxt="n",
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "Log-linear model",
       col = "red")
  text(x=210, y=320, labels="B2")
  abline(0,1)
  
  plot(First.2019.all[1,] ~ First.2019.all[4,],
       xlab = "",
       xlim = c(190,350),
       ylab = "",
       xaxt="n",
       yaxt="n",
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "Exponential plateau model",
       col = "red")
  text(x=210, y=320, labels="B3")
  abline(0,1)
  
  #Automated profiler
  plot(Jz.YSI.list[[2]] ~ alpha.YSI.list[[2]][,1],
       las = 1,
       xlim = c(0,0.25),
       ylim = c(0.05, 0.40),
       xaxt = "n",
       xlab = "",
       # xlab = expression(alpha(z)),
       ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       main = "",cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       pch = 16, col = "chartreuse3")
  text(x=0.03, y=0.38, labels="A4")
  abline(Ar.Jz.alpha.YSI.2018, col = "chartreuse3", lty = 2)
  points(Jz.YSI.list[[1]][-c(1:3)] ~ alpha.YSI.list[[1]][,1],
         pch = 16, col = "blue",cex = 1.4)
  abline(Ar.Jz.alpha.YSI.2017, col = "blue")
  points(Jz.YSI.list[[3]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "red",cex = 1.4)
  # abline(Ar.Jz.alpha.YSI.2019, col = "red")
  points(Jz.YSI.list[[4]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "pink",cex = 1.4)
  # abline(Ar.Jz.alpha.YSI.2020, col = "pink")
  points(Jz.YSI.list[[5]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "#73BAE6",cex = 1.4)
  abline(Ar.Jz.alpha.YSI.2021, col = "#73BAE6", lty = 2)
  # legend("topleft", legend = c("2017","2018","2019","2020", "2021"),
  #        text.col = c("blue","chartreuse3","red","pink", "#73BAE6"))
  
  
  plot(Jz.YSI.list[[2]] ~ alpha.YSI.list[[2]][,1],
       las = 1,
       xlim = c(0,0.25),
       ylim = c(0.05, 0.40),
       xaxt = "n",
       yaxt = "n",
       ylab = "",
       xlab = "",
       # xlab = expression(alpha(z)),
       # ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       main = "",cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       pch = 16, col = "chartreuse3")
  text(x=0.03, y=0.38, labels="A5")
  # curve(coef(Ar.Jz.alpha.YSI.2018.log)[1]*log10(x*coef(Ar.Jz.alpha.YSI.2018.log)[2]),
  #       add=T, col ="chartreuse3")
  points(Jz.YSI.list[[1]][-c(1:3)] ~ alpha.YSI.list[[1]][,1],
         pch = 16, col = "blue",cex = 1.4)
  # curve(coef(Ar.Jz.alpha.YSI.2017.log)[1]*log10(x*coef(Ar.Jz.alpha.YSI.2017.log)[2]),
  #       add=T, col ="blue")
  points(Jz.YSI.list[[3]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "red",cex = 1.4)
  curve(coef(Ar.Jz.alpha.YSI.2019.log)[1]*log10(x*coef(Ar.Jz.alpha.YSI.2019.log)[2]),
        add=T, col ="red")
  points(Jz.YSI.list[[4]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "pink",cex = 1.4)
  curve(coef(Ar.Jz.alpha.YSI.2020.log)[1]*log10(x*coef(Ar.Jz.alpha.YSI.2020.log)[2]),
        add=T, col ="pink")
  points(Jz.YSI.list[[5]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "#73BAE6",cex = 1.4)
  # curve(coef(Ar.Jz.alpha.YSI.2021.log)[1]*log10(x*coef(Ar.Jz.alpha.YSI.2021.log)[2]),
  #       add=T, col ="#73BAE6")
  
  plot(Jz.YSI.list[[2]] ~ alpha.YSI.list[[2]][,1],
       las = 1,
       xlim = c(0,0.25),
       ylim = c(0.05, 0.40),
       xaxt = "n",
       yaxt = "n",
       ylab = "",
       xlab = "",
       # xlab = expression(alpha(z)),
       # ylab = expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       main = "",cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       pch = 16, col = "chartreuse3")
  text(x=0.03, y=0.38, labels="A6")
  curve(coef(Ar.Jz.alpha.YSI.2018.exp)[1] - (coef(Ar.Jz.alpha.YSI.2018.exp)[1] - 
                                               coef(Ar.Jz.alpha.YSI.2018.exp)[2])*
          exp(-1*coef(Ar.Jz.alpha.YSI.2018.exp)[3]*x), 
        add=T, col = "chartreuse3")
  
  points(Jz.YSI.list[[1]][-c(1:3)] ~ alpha.YSI.list[[1]][,1],
         pch = 16, col = "blue",cex = 1.4)
  curve(coef(Ar.Jz.alpha.YSI.2017.exp)[1] - (coef(Ar.Jz.alpha.YSI.2017.exp)[1] - 
                                               coef(Ar.Jz.alpha.YSI.2017.exp)[2])*
          exp(-1*coef(Ar.Jz.alpha.YSI.2017.exp)[3]*x), 
        add=T, col = "blue")
  
  points(Jz.YSI.list[[3]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "red",cex = 1.4)
  curve(coef(Ar.Jz.alpha.YSI.2019.exp)[1] - (coef(Ar.Jz.alpha.YSI.2019.exp)[1] - 
                                               coef(Ar.Jz.alpha.YSI.2019.exp)[2])*
          exp(-1*coef(Ar.Jz.alpha.YSI.2019.exp)[3]*x), 
        add=T, col = "red")
  
  points(Jz.YSI.list[[4]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "pink",cex = 1.4)
  curve(coef(Ar.Jz.alpha.YSI.2020.exp)[1] - (coef(Ar.Jz.alpha.YSI.2020.exp)[1] - 
                                               coef(Ar.Jz.alpha.YSI.2020.exp)[2])*
          exp(-1*coef(Ar.Jz.alpha.YSI.2020.exp)[3]*x), 
        add=T, col = "pink")
  
  points(Jz.YSI.list[[5]] ~ alpha.YSI.list[[2]][,1],
         pch = 16, col = "#73BAE6",cex = 1.4)
  curve(coef(Ar.Jz.alpha.YSI.2021.exp)[1] - (coef(Ar.Jz.alpha.YSI.2021.exp)[1] - 
                                               coef(Ar.Jz.alpha.YSI.2021.exp)[2])*
          exp(-1*coef(Ar.Jz.alpha.YSI.2021.exp)[3]*x), 
        add=T, col = "#73BAE6")
  
  #First day of anaerobic conditions, automated profiler
  plot(First.2019.all[1,] ~ First.2019.all[5,],
       xlab = "",
       xlim = c(190,350),
       ylab = "First day of anoxia (observed)",
       xaxt = "n",
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "",
       col = "red")
  text(x=210, y=320, labels="B4")
  abline(0,1)
  
  plot(First.2019.all[1,] ~ First.2019.all[6,],
       xlab = "",
       xlim = c(190,350),
       ylab = "",
       xaxt="n",
       yaxt="n",
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "",
       col = "red")
  text(x=210, y=320, labels="B5")
  abline(0,1)
  
  plot(First.2019.all[1,] ~ First.2019.all[7,],
       xlab = "",
       xlim = c(190,350),
       ylab = "",
       xaxt="n",
       yaxt="n",
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "",
       col = "red")
  text(x=210, y=320, labels="B6")
  abline(0,1)
  
  
  #2 oxic profiles
  plot(Jz.lowres.slice.list[[1]][,1] ~ ar.alpha.list[[1]][,1],
       las = 1,
       xlim = c(0,0.25),
       ylim = c(0.05, 0.2),
       xlab = "",
       ylab = "",#expression(J[z]~(mg~O[2]~L^-1~d^-1)),
       main = "",cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       pch = 16, col = "blue")
  text(x=0.03, y=0.19, labels="A7")
  # abline(best.2017.short.lm, col = "blue")
  points(Jz.lowres.slice.list[[2]][,1] ~ ar.alpha.list[[2]][,1],
         pch = 16, cex = 1.4,col = "chartreuse3")
  abline(Ar.Jz.lm.best.list[[2]][[1]], col = "chartreuse3")
  points(Jz.lowres.slice.list[[3]][,1] ~ ar.alpha.list[[3]][,1],
         pch = 16, cex = 1.4,col = "red")
  # abline(Ar.Jz.lm.best.list[[3]][[1]], col = "red")
  points(Jz.lowres.slice.list[[4]][,1] ~ ar.alpha.list[[4]][,1],
         pch = 16, cex = 1.4,col = "pink")
  # abline(Ar.Jz.lm.best.list[[4]][[1]], col = "pink")
  points(Jz.lowres.slice.list[[5]][,1] ~ ar.alpha.list[[5]][,1],
         pch = 16, cex = 1.4,col = "#00B0F0")
  abline(Ar.Jz.lm.best.list[[5]][[1]], col = "#00B0F0", lty = 2)
  
  plot(Jz.lowres.slice.list[[1]][,1] ~ ar.alpha.list[[1]][,1],
       las = 1,
       xlim = c(0,0.25),
       ylim = c(0.05, 0.2),
       xlab = expression(alpha(z)),
       ylab = "",
       yaxt = "n",cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       pch = 16, col = "blue")
  text(x=0.03, y=0.19, labels="A8")
  curve(coef(Ar.Jz.lm.best.list[[1]][[2]])[1]*log10(x*coef(Ar.Jz.lm.best.list[[1]][[2]])[2]),
        add=T, col ="blue", lty = 2)
  points(Jz.lowres.slice.list[[2]][,1] ~ ar.alpha.list[[2]][,1],
         pch = 16, cex = 1.4,col = "chartreuse3")
  # curve(coef(Ar.Jz.lm.best.list[[2]][[2]])[1]*log10(x*coef(Ar.Jz.lm.best.list[[2]][[2]])[2]),
  #       add=T, col ="chartreuse3", lty = 2)
  points(Jz.lowres.slice.list[[3]][,1] ~ ar.alpha.list[[2]][,1],
         pch = 16, cex = 1.4,col = "red")
  curve(coef(Ar.Jz.lm.best.list[[3]][[2]])[1]*log10(x*coef(Ar.Jz.lm.best.list[[3]][[2]])[2]),
        add=T, col ="red", lty = 2)
  points(Jz.lowres.slice.list[[4]][,1] ~ ar.alpha.list[[2]][,1],
         pch = 16, cex = 1.4,col = "pink")
  curve(coef(Ar.Jz.lm.best.list[[4]][[2]])[1]*log10(x*coef(Ar.Jz.lm.best.list[[4]][[2]])[2]),
        add=T, col ="pink")
  points(Jz.lowres.slice.list[[5]][,1] ~ ar.alpha.list[[5]][,1],
         pch = 16, cex = 1.4,col = "#00B0F0")
  # curve(coef(Ar.Jz.lm.best.list[[5]][[2]])[1]*log10(x*coef(Ar.Jz.lm.best.list[[5]][[2]])[2]),
  #       add=T, col ="#00B0F0")
  
  plot(Jz.lowres.slice.list[[1]][,1] ~ ar.alpha.list[[1]][,1],
       las = 1,
       xlim = c(0,0.25),
       ylim = c(0.05, 0.2),
       xlab = "",
       ylab = "",
       yaxt = "n",cex = 1.4, cex.axis=0.9, cex.lab = 0.8,cex.main = 1.2,
       pch = 16, col = "blue")
  text(x=0.03, y=0.19, labels="A9")
  curve(coef(Ar.Jz.lm.best.list[[1]][[3]])[1] - (coef(Ar.Jz.lm.best.list[[1]][[3]])[1] - 
                                                   coef(Ar.Jz.lm.best.list[[1]][[3]])[2])*
          exp(-1*coef(Ar.Jz.lm.best.list[[1]][[3]])[3]*x), 
        add=T, col = "blue")
  
  points(Jz.lowres.slice.list[[2]][,1] ~ ar.alpha.list[[2]][,1], pch = 16, cex = 1.4,col = "chartreuse3")
  curve(coef(Ar.Jz.lm.best.list[[2]][[3]])[1] - (coef(Ar.Jz.lm.best.list[[2]][[3]])[1] - 
                                                   coef(Ar.Jz.lm.best.list[[2]][[3]])[2])*
          exp(-1*coef(Ar.Jz.lm.best.list[[2]][[3]])[3]*x), 
        add=T, col = "chartreuse3")
  
  points(Jz.lowres.slice.list[[3]][,1] ~ ar.alpha.list[[3]][,1], pch = 16, cex = 1.4,col = "red")
  curve(coef(Ar.Jz.lm.best.list[[3]][[3]])[1] - (coef(Ar.Jz.lm.best.list[[3]][[3]])[1] - 
                                                   coef(Ar.Jz.lm.best.list[[3]][[3]])[2])*
          exp(-1*coef(Ar.Jz.lm.best.list[[3]][[3]])[3]*x), 
        add=T, col = "red")
  
  points(Jz.lowres.slice.list[[4]][,1] ~ ar.alpha.list[[4]][,1], pch = 16, cex = 1.4,col = "pink")
  curve(coef(Ar.Jz.lm.best.list[[4]][[3]])[1] - (coef(Ar.Jz.lm.best.list[[4]][[3]])[1] -
                                                   coef(Ar.Jz.lm.best.list[[4]][[3]])[2])*
          exp(-1*coef(Ar.Jz.lm.best.list[[4]][[3]])[3]*x),
        add=T, col = "pink")
  
  points(Jz.lowres.slice.list[[5]][,1] ~ ar.alpha.list[[5]][,1], pch = 16, cex = 1.4,col = "#00B0F0")
  curve(coef(Ar.Jz.lm.best.list[[5]][[3]])[1] - (coef(Ar.Jz.lm.best.list[[5]][[3]])[1] -
                                                   coef(Ar.Jz.lm.best.list[[5]][[3]])[2])*
          exp(-1*coef(Ar.Jz.lm.best.list[[5]][[3]])[3]*x),
        add=T, col = "#00B0F0")
  
  #First day of anaerobic conditions, 2 oxic profiles
  plot(First.2019.all[1,] ~ First.2019.all[8,],
       xlim = c(190,350),
       xlab = "",
       ylab = "",
       pch = 1,
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "",
       col = "red")
  text(x=210, y=320, labels="B7")
  abline(0,1)
  
  plot(First.2019.all[1,] ~ First.2019.all[9,],
       xlim = c(190,350),
       xlab = "First day of anoxia (modelled)",
       ylab = "",
       yaxt = "n",
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "",
       col = "red")
  text(x=210, y=320, labels="B8")
  abline(0,1)
  
  plot(First.2019.all[1,] ~ First.2019.all[10,],
       xlim = c(190,350),
       xlab = "",
       ylab = "",
       yaxt = "n",
       pch = 16,
       las = 1,cex = 1.4, cex.axis=0.8, cex.lab = 0.9,
       main = "",
       col = "red")
  text(x=210, y=320, labels="B9")
  abline(0,1)
  
  dev.off()
}
```

#We will now proceed to match nutrients with anaerobic distribution.
```{R AD - nutrients}
#Load nutrient data
NH4_SRP.data = read.csv("./SRP-NH4.Arendsee.2017-2019.csv")
#Create a Day of Year variable
NH4_SRP.data$DOY = as.numeric(strftime(NH4_SRP.data$Date, format = "%j"))

#Add column for anaerobic duration
NH4_SRP.data$AD <- NA

#Reselect alpha(z) values
ar.alpha <- Bats %>% dplyr::slice(30:50)

#We start with 2017 and will do it for the other years below. I did not create a for loop here because each year has specificity that would be too hard to include in a for loop. Hence, there are hundreds on lines of codes because I repeat the same large chunk of code multiple times.
#Note. This section is basically the code behind O2_Model.R, but slightly modified. This is the section you probably need for your own Anaerobic Duration calculation
#In the first example below, the object "Ar.Jz.alpha.YSI.2017.exp" is the nlsLM object for the exponential plateau fit between Jz and alpha(z) in lake Arendsee, 2017. We use the three coefficients of that equation, the model type and a vector of alpha(z) values. The output object is a vector containing oxygen consumption rates for each depth (in mg O2 L-1 d-1), which we will use a few lines afters to calculate daily oxygen profiles, assuming that the entire hypolimnion starts fully oxygenated.

#2017
#Create new oxygen consumption rate vector with higher resolution (some nutrients were taken at 0.5m interval)
  O2.consumption.rate <- O2.decay(b = coef(Ar.Jz.alpha.YSI.2017.exp)[2], 
                                  k = coef(Ar.Jz.alpha.YSI.2017.exp)[3],
                                  JzMax = coef(Ar.Jz.alpha.YSI.2017.exp)[1],
                                  Model = "Exp",
                                  alpha = ar.alpha[,2])
  nday = 271 #Numbers of days to simulate
  dbeg = 135 #Arbitrary number
  
  #Create an empty matrix with the good dimensions 
  O2.modelled.matrix = as.data.frame(matrix(nrow=length(O2.consumption.rate), ncol = nday))
  #Set first profile at 12.3 mg O2/L (100% saturation in Arendsee)
  O2.modelled.matrix[,1] = 12.3
  
  #Model daily profiles based on consumption rate
  for(i in 2:ncol(O2.modelled.matrix)){
    O2.modelled.matrix[,i] = 12.3 - O2.consumption.rate*(i-1)
  }
  #Change negative values to 0
  O2.modelled.matrix[O2.modelled.matrix<0] = 0
  
  #Anaerobic duration threshold, here 0.5
  O2.threshold = 0.5
  
  #Change values below threshold to 1, others to 0
  O2.modelled.matrix.bin <- ifelse(O2.modelled.matrix <= O2.threshold, 1, 0)
  
  #Create new object to calculate Anaerobic duration
  O2.modelled.AnaerobicDuration <- O2.modelled.matrix.bin
  for(j in 2:ncol(O2.modelled.AnaerobicDuration)){
    temp = O2.modelled.AnaerobicDuration[,j-1]+O2.modelled.AnaerobicDuration[,j]
    temp = ifelse(temp > O2.modelled.AnaerobicDuration[,j-1], temp, 0)
    O2.modelled.AnaerobicDuration[,j] = temp
  }
  
  #Add a column for depths and colnames
  O2.modelled.AnaerobicDuration = cbind(ar.alpha[,1], O2.modelled.AnaerobicDuration)
  colnames(O2.modelled.AnaerobicDuration)[1] = "Depth"
  colnames(O2.modelled.AnaerobicDuration)[2:(nday+1)] = paste(seq(dbeg,dbeg+nday-1,1))
  
  #Transform into dataframe
  O2.modelled.AnaerobicDuration = as.data.frame(O2.modelled.AnaerobicDuration)
  
  #Transform to long format
  O2.modelled.AnaerobicDuration.long = pivot_longer(O2.modelled.AnaerobicDuration,
                                            cols = 2:ncol(O2.modelled.AnaerobicDuration),
                                            names_to = "DOY", values_to = "AD")
  O2.modelled.AnaerobicDuration.long = as.data.frame(O2.modelled.AnaerobicDuration.long)
  
  #See results in ugly graph
  #NOTE. The object input is really long because ggplot doesn't handle well uneven Y axis (there ar1 and 0.5m steps). To avoid rendering issues, I removed the problematic depths just for visualisation   
  print(ggplot(O2.modelled.AnaerobicDuration.long[O2.modelled.AnaerobicDuration.long$Depth!=48.5 & O2.modelled.AnaerobicDuration.long$Depth!=47.5,])+
          geom_raster(aes(x=as.numeric(DOY), y=Depth, fill=AD))+
          scale_y_reverse()+
          scale_fill_gradient(low="white", high="red")+
          ggtitle(2017, "- Progression of anoxia, excluding lake mixing")+
          xlab("DOY")+
          labs(fill = bquote(AD~(.(O2.threshold)~mg~L^-1))))
  
  #Align the modeled oxygen profile with the real data to associate Anaerobic duration with 
  #reduced compounds concentration using the lowest sum of squared difference
  year.i = grep(2017, YSI.list)
  Ar.YSI = read.csv(paste0("./",YSI.list[year.i]))
  #Remove NA at bottom of profiles with niminum value of the profile
  for(j in 1:ncol(Ar.YSI))
  {
    if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
  }
  
  #Transform to long format
  Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
  #Split Date into Year, months and day
  Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")
  
  Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))
  
  Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                               mm = substring(Date, 6, 7),
                               dd = substring(Date, 9, 10))
  
  #Keep only values below 30 meters to match other methods
  Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]
  
  #Order the dataframe
  Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]
  
  #Select the date matching chemistry observations
  #In 2020, Arendsee was sampled on 2020-05-26 (DOY 147) and 2020-08-25
  #Until now, the matrix of daily modelled oxygen profiles starts at 100% saturation and gradually declines using the empirical oxygen consumption rate, but no date (or day of year) are associated. This next step is to set dates for the modelled daily oxygen profiles using observations. We use the lowest RMSE to "match" modelled profiles with observations. In the manuscript, this step is described in the paragraph just after equation 2.
  
  Ar.first.profile = Ar.long.deep[,-2] %>% filter(mm == "05", dd == "26") %>% 
    group_by(Depth_m) %>% summarize(DO_mgL = mean(DO_mgL))
  
  Match.temp.1 = O2.modelled.matrix[-c(19,21),]
  for(i in 1:ncol(O2.modelled.matrix)){
    Match.temp.1[,i] = O2.modelled.matrix[,i] - Ar.first.profile[,2]
  }
  
  Match.temp.1.abs = (Match.temp.1)^2
  Match.day.1 = which(colSums(Match.temp.1.abs) == min(colSums(Match.temp.1.abs)))
  Match.day.1.RMSE = sqrt(sum(Match.temp.1.abs[Match.day.1])/length(Match.temp.1.abs[Match.day.1]))
  #Stratified early but makes sense with the data
  
  O2.modelled.AnaerobicDuration.2017 = O2.modelled.AnaerobicDuration
  #Make sure column 54 is DOY 147
  #Note. DOY 147 is the day when the YSI profile was taken
  #Note. Column 54 is the value of Match.day.1 +1 for the first column (depths) of this matrix
  colnames(O2.modelled.AnaerobicDuration.2017)[c(2:ncol(O2.modelled.AnaerobicDuration.2017))] = seq(95,ncol(O2.modelled.AnaerobicDuration.2017)+93,1)
  
  #Transform to long format
  O2.modelled.AnaerobicDuration.2017.long = pivot_longer(O2.modelled.AnaerobicDuration.2017,
                                                 cols = 2:ncol(O2.modelled.AnaerobicDuration.2017),
                                                 names_to = "DOY",
                                                 values_to = "AD")
  
  #We use these values ONLY to fill gaps from observed Anaerobic duration
  #Change TRUE/FALSE with 1/0
  YSI.aa.2017[YSI.aa.2017==FALSE] = 0
  YSI.aa.2017[YSI.aa.2017==TRUE] = 1
  
  #Compile Anaerobic duration
  YSI.aa.2017.cumul = YSI.aa.2017[,-c(1:120)]
  for(j in 2:ncol(YSI.aa.2017.cumul)){
    temp = YSI.aa.2017.cumul[,j-1]+YSI.aa.2017.cumul[,j]
    temp = ifelse(temp > YSI.aa.2017.cumul[,j-1], temp, 0)
    YSI.aa.2017.cumul[,j] = temp
  }
  
  YSI.aa.2017.cumul = cbind(Depth=c(30:48),YSI.aa.2017.cumul)
  
  YSI.aa.2017.cumul.long = pivot_longer(as.data.frame(YSI.aa.2017.cumul),
                                        cols = 2:ncol(YSI.aa.2017.cumul), 
                                        names_to = "DOY",
                                        values_to = "AD")
  
  #Join modelled 48.5m to observed Anaerobic duration
  Combined.long = rbind(YSI.aa.2017.cumul.long,O2.modelled.AnaerobicDuration.2017.long[O2.modelled.AnaerobicDuration.2017.long$Depth==48.5,])
  
  #See if it makes sense
  Combined.aa.2017 = pivot_wider(Combined.long, names_from = DOY, values_from = AD)
  #It does, continue with the long format
  
  #Subset nutrient data by years
  NH4_SRP.data.2017 = NH4_SRP.data[grep("2017", NH4_SRP.data$Date),]
  NH4_SRP.data.2017[NH4_SRP.data.2017$Depth_m=="48.5", "Depth_m"] = 48
  
  #Match calculated Anaerobic duration with nutrient data
  for(i in 1:nrow(NH4_SRP.data.2017)){
    if(NH4_SRP.data.2017[i,"DOY"] < 147){
      NH4_SRP.data.2017[i,"AD"] = 0
      next
    }
    NH4_SRP.data.2017[i,"AD"] = Combined.long [Combined.long $Depth == NH4_SRP.data.2017[i, "Depth_m"] & Combined.long $DOY==NH4_SRP.data.2017[i,"DOY"],"AD"]
  }
  
  #Remove early dates (DOY < 94). This is the first column above where column 54 is DOY 147 (see above)
  NH4_SRP.data.2017 = NH4_SRP.data.2017[-which(NH4_SRP.data.2017$DOY<=94),]
  #Include oxic duration
  Oxic.df.2017 = data.frame(DOY = c(95:365), OxicDuration = c(1:271))
  NH4_SRP.data.2017$OxicDuration <- Oxic.df.2017[match(NH4_SRP.data.2017$DOY, Oxic.df.2017$DOY),"OxicDuration"]
  
######################Repeat all the steps before with 2019 and 2020#####################
  
#With 2019
#Create new anaerobic duration vector with higher resolution
  O2.consumption.rate <- O2.decay(b = coef(Ar.Jz.alpha.YSI.2019.exp)[2], 
                                  k = coef(Ar.Jz.alpha.YSI.2019.exp)[3],
                                  JzMax = coef(Ar.Jz.alpha.YSI.2019.exp)[1],
                                  Model = "Exp",
                                  alpha = ar.alpha[,2])
  # nday = 271
  # dbeg = 135 
  
  #Create empty matrix with the good dimensions
  O2.modelled.matrix = as.data.frame(matrix(nrow=length(O2.consumption.rate), ncol = nday))
  #Set first profile at 12.3 mg O2/L
  O2.modelled.matrix[,1] = 12.3
  
  #Model daily profiles based on consumption rate
  for(i in 2:ncol(O2.modelled.matrix)){
    O2.modelled.matrix[,i] = 12.3 - O2.consumption.rate*(i-1)
  }
  #Change negative values to 0
  O2.modelled.matrix[O2.modelled.matrix<0] = 0
  
  #anaerobic duration threshold
  O2.threshold = 0.5
  
  #Change values below threshold to 1, others to 0
  O2.modelled.matrix.bin <- ifelse(O2.modelled.matrix <= O2.threshold, 1, 0)
  
  #Create new object to calculate anaerobic duration
  O2.modelled.AnaerobicDuration <- O2.modelled.matrix.bin
  for(j in 2:ncol(O2.modelled.AnaerobicDuration)){
    temp = O2.modelled.AnaerobicDuration[,j-1]+O2.modelled.AnaerobicDuration[,j]
    temp = ifelse(temp > O2.modelled.AnaerobicDuration[,j-1], temp, 0)
    O2.modelled.AnaerobicDuration[,j] = temp
  }
  
  #Add a column for depths and colnames
  O2.modelled.AnaerobicDuration = cbind(ar.alpha[,1], O2.modelled.AnaerobicDuration)
  colnames(O2.modelled.AnaerobicDuration)[1] = "Depth"
  colnames(O2.modelled.AnaerobicDuration)[2:(nday+1)] = paste(seq(dbeg,dbeg+nday-1,1))
  
  #Transform into dataframe
  O2.modelled.AnaerobicDuration = as.data.frame(O2.modelled.AnaerobicDuration)
  
  #Transform to long format
  O2.modelled.AnaerobicDuration.long = pivot_longer(O2.modelled.AnaerobicDuration,
                                            cols = 2:ncol(O2.modelled.AnaerobicDuration),
                                            names_to = "DOY", values_to = "AD")
  O2.modelled.AnaerobicDuration.long = as.data.frame(O2.modelled.AnaerobicDuration.long)
  
  #See results in ugly graph
  print(ggplot(O2.modelled.AnaerobicDuration.long[O2.modelled.AnaerobicDuration.long$Depth!=48.5 & O2.modelled.AnaerobicDuration.long$Depth!=47.5,])+
          geom_raster(aes(x=as.numeric(DOY), y=Depth, fill=AD))+
          scale_y_reverse()+
          scale_fill_gradient(low="white", high="red")+
          ggtitle(2019, "- Progression of anoxia, excluding lake mixing")+
          xlab("DOY")+
          labs(fill = bquote(AD~(.(O2.threshold)~mg~L^-1))))
  
  #Align the modeled oxygen profile with the real data to associate anaerobic duration with 
  #reduced compounds concentration using the lowest sum of squared difference
  year.i = grep(2019, YSI.list)
  Ar.YSI = read.csv(paste0("./",YSI.list[year.i]))
  #Remove NA at bottom of profiles with niminum value of the profile
  for(j in 1:ncol(Ar.YSI))
  {
    if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
  }
  
  #Transform to long format
  Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
  #Split Date into Year, months and day
  Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")
  
  Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))
  
  Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                               mm = substring(Date, 6, 7),
                               dd = substring(Date, 9, 10))
  
  #Keep only values below 30 meters to match other methods
  Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]
  
  #Order the dataframe
  Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]
  
  #Select the date matching chemistry observations
  #In 2020, Arendsee was sampled on 2020-05-26 (DOY 147) and 2020-08-25	
  Ar.first.profile = Ar.long.deep[,-2] %>% filter(mm == "05", dd == "26") %>% 
    group_by(Depth_m) %>% summarize(DO_mgL = mean(DO_mgL))
  
  Match.temp.1 = O2.modelled.matrix[-c(19:21),]
  for(i in 1:ncol(O2.modelled.matrix)){
    Match.temp.1[,i] = O2.modelled.matrix[,i] - Ar.first.profile[,2]
  }
  
  Match.temp.1.abs = (Match.temp.1)^2
  Match.day.1 = which(colSums(Match.temp.1.abs) == min(colSums(Match.temp.1.abs)))
  Match.day.1.RMSE = sqrt(sum(Match.temp.1.abs[Match.day.1])/length(Match.temp.1.abs[Match.day.1]))
  
  
  O2.modelled.AnaerobicDuration.2019 = O2.modelled.AnaerobicDuration
  #Make sure column 41 is DOY 147
  colnames(O2.modelled.AnaerobicDuration.2019)[c(2:ncol(O2.modelled.AnaerobicDuration.2019))] = seq(108,ncol(O2.modelled.AnaerobicDuration.2019)+106,1)
  
  
  #Change a few values at 47.5 and 48m deep based on independent profiles 
  O2.modelled.AnaerobicDuration.2019[19,c(93:272)] = seq(1,180,1)
  O2.modelled.AnaerobicDuration.2019[20,c(92:272)] = seq(1,181,1)
  O2.modelled.AnaerobicDuration.2019[21,c(91:272)] = seq(1,182,1)
  
  
  #Transform to long format
  O2.modelled.AnaerobicDuration.2019.long = pivot_longer(O2.modelled.AnaerobicDuration.2019, cols = 2:ncol(O2.modelled.AnaerobicDuration.2019), names_to = "DOY", values_to = "AD")
  
  
  #We use these values ONLY to fill gaps from observed anaerobic duration
  #Change TRUE/FALSE with 1/0
  YSI.aa.2019[YSI.aa.2019==FALSE] = 0
  YSI.aa.2019[YSI.aa.2019==TRUE] = 1
  
  #Compile anaerobic duration
  YSI.aa.2019.cumul = YSI.aa.2019[,-c(1:120)]
  for(j in 2:ncol(YSI.aa.2019.cumul)){
    temp = YSI.aa.2019.cumul[,j-1]+YSI.aa.2019.cumul[,j]
    temp = ifelse(temp > YSI.aa.2019.cumul[,j-1], temp, 0)
    YSI.aa.2019.cumul[,j] = temp
  }
  
  YSI.aa.2019.cumul = cbind(Depth=c(30:47),YSI.aa.2019.cumul)
  
  YSI.aa.2019.cumul.long = pivot_longer(as.data.frame(YSI.aa.2019.cumul),
                                        cols = 2:ncol(YSI.aa.2019.cumul), 
                                        names_to = "DOY",
                                        values_to = "AD")
  
  #Join modelled 48.5m to observed anaerobic duration
  Combined.long = rbind(YSI.aa.2019.cumul.long,
                        O2.modelled.AnaerobicDuration.2019.long[O2.modelled.AnaerobicDuration.2019.long$Depth==47.5,],
                        O2.modelled.AnaerobicDuration.2019.long[O2.modelled.AnaerobicDuration.2019.long$Depth==48,])
  
  #See if it makes sense
  Combined.aa.2019 = pivot_wider(Combined.long, names_from = DOY, values_from = AD)
  #Some depths are weird, but let's see
  
  
  #Subset nutrient data by years
  NH4_SRP.data.2019 = NH4_SRP.data[grep("2019", NH4_SRP.data$Date),]
  NH4_SRP.data.2019[NH4_SRP.data.2019$Depth_m==47.8,"Depth_m"] = 48
  NH4_SRP.data.2019[NH4_SRP.data.2019$Depth_m==47.6,"Depth_m"] = 48
  NH4_SRP.data.2019[NH4_SRP.data.2019$Depth_m==47.5,"Depth_m"] = 47
  
  #Match calculated anaerobic duration with nutrient data
  for(i in 1:nrow(NH4_SRP.data.2019)){
    if(NH4_SRP.data.2019[i,"DOY"] < 147){
      NH4_SRP.data.2019[i,"AD"] = 0
      next
    }
    NH4_SRP.data.2019[i,"AD"] = Combined.long[Combined.long $Depth == NH4_SRP.data.2019[i,"Depth_m"] & Combined.long $DOY==NH4_SRP.data.2019[i,"DOY"],"AD"]
    }
  
    #Remove early dates (DOY<120)
  NH4_SRP.data.2019 = NH4_SRP.data.2019[-which(NH4_SRP.data.2019$DOY<108),]

  #Remove 1 line that looks like an outlier (lowest values for no apparent reason)
  NH4_SRP.data.2019 =  NH4_SRP.data.2019[-which(NH4_SRP.data.2019$SRP_mgL==min(NH4_SRP.data.2019$SRP_mgL)),]
  
  #Add Oxic Duration
  Oxic.df.2019 = data.frame(DOY = c(108:368), OxicDuration = c(1:261))
  NH4_SRP.data.2019$OxicDuration <- Oxic.df.2019[match(NH4_SRP.data.2019$DOY, Oxic.df.2019$DOY),"OxicDuration"]

#With 2020
Chemistry <- read.csv("./SRP-NH4.Arendsee.2020.csv")

#Create new anaerobic duration vector with higher resolution
  ar.alpha.2020.NH4_SRP = ar.alpha
  
  O2.consumption.rate <- O2.decay(b = coef(Ar.Jz.alpha.YSI.2020.exp)[2], 
                                  k = coef(Ar.Jz.alpha.YSI.2020.exp)[3],
                                  JzMax = coef(Ar.Jz.alpha.YSI.2020.exp)[1],
                                  Model = "Exp",
                                  alpha = ar.alpha.2020.NH4_SRP[,2])
  # nday = 271 #Numbers of days to simulate
  # dbeg = 135 #First day of simulated stratification
  
  #Create empty matrix with the good dimensions
  O2.modelled.matrix = as.data.frame(matrix(nrow=length(O2.consumption.rate), ncol = nday))
  #Set first profile at 12.3 mg O2/L
  O2.modelled.matrix[,1] = 12.3
  
  #Model daily profiles based on consumption rate
  for(i in 2:ncol(O2.modelled.matrix)){
    O2.modelled.matrix[,i] = 12.3 - O2.consumption.rate*(i-1)
  }
  #Change negative values to 0
  O2.modelled.matrix[O2.modelled.matrix<0] = 0
  
  #anaerobic duration threshold
  O2.threshold = 0.5
  
  #Change values below threshold to 1, others to 0
  O2.modelled.matrix.bin <- ifelse(O2.modelled.matrix <= O2.threshold, 1, 0)
  
  #Create new object to calculate anaerobic duration
  O2.modelled.AnaerobicDuration <- O2.modelled.matrix.bin
  for(j in 2:ncol(O2.modelled.AnaerobicDuration)){
    temp = O2.modelled.AnaerobicDuration[,j-1]+O2.modelled.AnaerobicDuration[,j]
    temp = ifelse(temp > O2.modelled.AnaerobicDuration[,j-1], temp, 0)
    O2.modelled.AnaerobicDuration[,j] = temp
  }
  
  #Add a column for depths and colnames
  O2.modelled.AnaerobicDuration = cbind(ar.alpha.2020.NH4_SRP[,1], O2.modelled.AnaerobicDuration)
  colnames(O2.modelled.AnaerobicDuration)[1] = "Depth"
  colnames(O2.modelled.AnaerobicDuration)[2:(nday+1)] = paste(seq(dbeg,dbeg+nday-1,1))
  
  #Transform into dataframe
  O2.modelled.AnaerobicDuration = as.data.frame(O2.modelled.AnaerobicDuration)
  
  #Transform to long format
  O2.modelled.AnaerobicDuration.long = pivot_longer(O2.modelled.AnaerobicDuration,
                                            cols = 2:ncol(O2.modelled.AnaerobicDuration),
                                            names_to = "DOY", values_to = "AD")
  O2.modelled.AnaerobicDuration.long = as.data.frame(O2.modelled.AnaerobicDuration.long)
  
  #See results in ugly graph
    print(ggplot(O2.modelled.AnaerobicDuration.long[O2.modelled.AnaerobicDuration.long$Depth!=48.5 & O2.modelled.AnaerobicDuration.long$Depth!=47.5,])+
          geom_raster(aes(x=as.numeric(DOY), y=Depth, fill=AD))+
          scale_y_reverse()+
          scale_fill_gradient(low="white", high="red")+
          ggtitle(2020, "- Progression of anoxia, excluding lake mixing")+
          xlab("DOY")+
          labs(fill = bquote(AD~(.(O2.threshold)~mg~L^-1))))

  
  #Align the modeled oxygen profile with the real data to associate anaerobic duration with 
  #reduced compounds concentration using the lowest sum of squared difference
  year.i = grep(2020, YSI.list)
  Ar.YSI = read.csv(paste0("./",YSI.list[year.i]))
  #Remove NA at bottom of profiles with niminum value of the profile
  for(j in 1:ncol(Ar.YSI))
  {
    if(any(is.na(Ar.YSI[,j]))) Ar.YSI[is.na(Ar.YSI[,j]),j] = min(Ar.YSI[,j], na.rm=T)
  }
  
  #Transform to long format
  Ar.long = pivot_longer(Ar.YSI, cols = 2:ncol(Ar.YSI), names_to = "Date", values_to = "DO_mgL")
  #Split Date into Year, months and day
  Ar.long$Date = substring(Ar.long$Date, 2,20) %>% strptime("%Y.%m.%d.%H.%M.%S")
  
  Ar.long$DOY = as.numeric(strftime(Ar.long$Date, format = "%j"))
  
  Ar.long = Ar.long %>% mutate(Year = substring(Date, 1, 4),
                               mm = substring(Date, 6, 7),
                               dd = substring(Date, 9, 10))
  
  #Keep only values below 30 meters to match other methods
  Ar.long.deep = Ar.long[Ar.long$Depth_m>=30,]
  
  #Order the dataframe
  Ar.long.deep = Ar.long.deep[order(Ar.long.deep$Depth_m ,Ar.long.deep$DOY),]
  
  #Select the date matching chemistry observations
  #In 2020, Arendsee was sampled on 2020-05-26 (DOY 147) and 2020-08-25	
  Ar.first.profile = Ar.long.deep[,-2] %>% filter(mm == "05", dd == "26") %>% 
    group_by(Depth_m) %>% summarize(DO_mgL = mean(DO_mgL))
  
  Match.temp.1 = O2.modelled.matrix[-c(14,20:22),]
  for(i in 1:ncol(O2.modelled.matrix)){
    Match.temp.1[,i] = O2.modelled.matrix[,i] - Ar.first.profile[,2]
  }
  
  Match.temp.1.abs = (Match.temp.1)^2
  Match.day.1 = which(colSums(Match.temp.1.abs) == min(colSums(Match.temp.1.abs)))
  Match.day.1.RMSE = sqrt(sum(Match.temp.1.abs[Match.day.1])/length(Match.temp.1.abs[Match.day.1]))
  
  O2.modelled.AnaerobicDuration.2020 = O2.modelled.AnaerobicDuration
  #Make sure column 37 is DOY 147
  colnames(O2.modelled.AnaerobicDuration.2020)[c(2:ncol(O2.modelled.AnaerobicDuration.2020))] = seq(112,ncol(O2.modelled.AnaerobicDuration.2020)+110,1)
  
  
  #Transform to long format
  O2.modelled.AnaerobicDuration.2020.long = pivot_longer(O2.modelled.AnaerobicDuration.2020,  cols = 2:ncol(O2.modelled.AnaerobicDuration.2020),  names_to = "DOY",  values_to = "AD")
  
  
  #We use these values ONLY to fill gaps from observed anaerobic duration
  #Change TRUE/FALSE with 1/0
  YSI.aa.2020[YSI.aa.2020==FALSE] = 0
  YSI.aa.2020[YSI.aa.2020==TRUE] = 1
  
  #Compile anaerobic duration
  
  YSI.aa.2020.cumul = YSI.aa.2020[,-c(1:111)]
  
  for(j in 2:ncol(YSI.aa.2020.cumul)){
    temp = YSI.aa.2020.cumul[,j-1]+YSI.aa.2020.cumul[,j]
    temp = ifelse(temp > YSI.aa.2020.cumul[,j-1], temp, 0)
    YSI.aa.2020.cumul[,j] = temp
  }
  
  YSI.aa.2020.cumul = cbind(Depth=c(30:47),YSI.aa.2020.cumul)
  
  YSI.aa.2020.cumul.long = pivot_longer(as.data.frame(YSI.aa.2020.cumul),
                                        cols = 2:ncol(YSI.aa.2020.cumul), 
                                        names_to = "DOY",
                                        values_to = "AD")
  
  
  #Join modelled 48.5m to observed anaerobic duration
  Combined.long = rbind(YSI.aa.2020.cumul.long,  O2.modelled.AnaerobicDuration.2020.long[O2.modelled.AnaerobicDuration.2020.long$Depth==48,])
  
  
  #See if it makes sense (it does)
  Combined.aa.2020 = pivot_wider(Combined.long, names_from = DOY, values_from = AD)
  
  Combined.aa.2020 <- as.data.frame(Combined.aa.2020) 
  Combined.aa.2020[19,] <- c(48,rep(0,57),seq(1,214,1))
  
  Combined.long <- pivot_longer(data = Combined.aa.2020,
                                cols = 2:ncol(Combined.aa.2020),
                                names_to = "DOY",
                                values_to = "AD")
  
  
  
  #See if it makes sense
  Combined.aa.2020 = pivot_wider(Combined.long, names_from = DOY, values_from = AD)
  
  
  
  #Select NH4 and SRP from other campaign dataset
  NH4_SRP.data.2020 <- filter(Chemistry, depth >=30)

  #Add DOY column
  NH4_SRP.data.2020$DOY = 0
  #Set DOY according to dates
  NH4_SRP.data.2020[grep("05-26", NH4_SRP.data.2020$date),"DOY"] = 147
  NH4_SRP.data.2020[grep("08-25", NH4_SRP.data.2020$date),"DOY"] = 238
  
  #Change 42.5 by 43m because it makes no difference with anaerobic duration
  NH4_SRP.data.2020[NH4_SRP.data.2020$depth==42.5,"depth"] = 42
  
  #Add anaerobic duration column to NH$-SRP table
  NH4_SRP.data.2020$AD = 0
  
  #Match calculated anaerobic duration with nutrient data
  for(i in 1:nrow(NH4_SRP.data.2020)){
    if(NH4_SRP.data.2020[i,"DOY"] < 147){
      NH4_SRP.data.2020[i,"AD"] = 0
      next
    }
    NH4_SRP.data.2020[i,"AD"] = Combined.long[Combined.long $Depth==NH4_SRP.data.2020[i,"depth"] & Combined.long $DOY==NH4_SRP.data.2020[i,"DOY"],"AD"]
  }
  
  #Change colnames to match other data source
  colnames(NH4_SRP.data.2020)[c(1,2,3,4)] = c("Date","Depth_m","SRP_mgL","NH4.N_mgL")
  
  #Change SRP into mgL (was in ugL)
  NH4_SRP.data.2020[,3] = NH4_SRP.data.2020[,3]/1000
  
  
  #Add Oxic Duration
  Oxic.df.2020 = data.frame(DOY = c(112:368),
                            OxicDuration = c(1:257))
  NH4_SRP.data.2020$OxicDuration <- Oxic.df.2020[match(NH4_SRP.data.2020$DOY, Oxic.df.2020$DOY),"OxicDuration"]
  
  
#Prepare to make Figure 2
  
#Add a column for pch values and years. This makes it easier to produce the graphs
NH4_SRP.data.2017$year = 1
NH4_SRP.data.2017$cex = 1
NH4_SRP.data.2019$year = 2
NH4_SRP.data.2019$cex = 2
NH4_SRP.data.2020$year = 3
NH4_SRP.data.2020$cex = 1


NH4_SRP.data.all = rbind(NH4_SRP.data.2017,
                         NH4_SRP.data.2019,
                         NH4_SRP.data.2020)

#Some relationships

#NH4 relationship with stagnant age
NH4.value = NH4_SRP.data.all$NH4.N_mgL[NH4_SRP.data.all$AD==0]
Oxic.value = NH4_SRP.data.all$OxicDuration[NH4_SRP.data.all$AD==0]
NH4.nlm.Oxic <- nlsLM(NH4.value ~ k/(Oxic.value + b),
                      start = list(k=1, b = 3))

NH4.value.rsq <- NH4.value[!is.na(Oxic.value)]
Oxic.value.rsq <- Oxic.value[!is.na(Oxic.value)]
Oxic.value.rsq <- Oxic.value.rsq[!is.na(NH4.value.rsq)]
NH4.value.rsq <- NH4.value.rsq[!is.na(NH4.value.rsq)]

NH4.nlm.Oxic.exp <- nlsLM(NH4.value.rsq ~ b*exp(-k*Oxic.value.rsq),
                          start = list(k=0.5,b=1))

pred <- predict(NH4.nlm.Oxic.exp, Oxic.value.rsq)
rss <- sum((pred - NH4.value.rsq) ^ 2)
tss <- sum((NH4.value.rsq - mean(NH4.value.rsq)) ^ 2)
rsq <- 1 - rss/tss
```

#This will produce Figure 2, middle row (Lake Arendseed ata). The other lakes were analysed in other scripts that are very similar to this one, then imported to construct the final figure. The figure is printed in the Output folder as a jpeg file.
```{R Figure 2 - Arendsee}

 {
  jpeg("./Output/Fig2.jpg", width = 8, height = 4, units = "in", res=300)
  #Oxic and anaerobic durations, SRP
  par(mfrow=c(1,2))
  par(mar=c(4,5,1,1.2)+.1)
  #Lake Arendsee
  plot(NH4_SRP.data.all$SRP_mgL[NH4_SRP.data.all$AD==0] ~ NH4_SRP.data.all$OxicDuration[NH4_SRP.data.all$AD==0],
       las = 1,
       ylim = c(0.1,0.36),
       ylab = expression(SRP~(mg~L^-1)),
       xlab = "",
       xlim = c(0, 400),
       xaxt = "n",
       pch = c(15,16,17)[NH4_SRP.data.all$year[NH4_SRP.data.all$AD==0]],
       col = c("#73BAE6", "#FFD700","#C7144C")[NH4_SRP.data.all$year[NH4_SRP.data.all$AD==0]],
       cex = 1.5,
       cex.lab = 1.5,
       cex.axis = 1.4)
  
  #Add "parenthesis" to bottom 4 points
  draw.arc(x=115, y=0.12, radius = 62, angle1 = 0.9*pi, angle2 = 1.08*pi)
  draw.arc(x=10.2, y=0.12, radius = 62, angle1 = 0 , angle2 = 0.1*pi)
  draw.arc(x=10.2, y=0.12, radius = 62, angle1 = 1.92*pi , angle2 = 2*pi)
  
  #Add axis values and labels
  axis(side = 1, at = c(0,100,200), labels = c("0", "100", "200"), cex.axis =1.4)
  mtext(expression(paste("Oxic duration")), side = 1, at = 126,line = 2, cex = 0.85, font=1)
  mtext(expression(paste("(days)")), side = 1, at = 126,line = 3, cex = 0.85, font=1)
  #Add panel letters and legend
  mtext("e", side = 3, at = -10)
  mtext("f", side = 3, at = 246)
  legend("topleft",
         pch = c(22,21,24),
         pt.cex = c(1.5,1.5,1.5),
         legend = c("2017", "2019","2020"),
         pt.bg = c("#73BAE6", "#FFD700","#C7144C"))
  
  #Cut the outside box to make it look like 2 panels
  clip(x1 = 252.2, x2 = 259.8, y1=0, y2=1) #Set boundaries where the plot can be modified
  abline(h = 0.0892, col = "white", lwd=2) #Add a thick white line on bottom axis
  abline(h = 0.37, col = "white", lwd=2) #Add a thick white line on top
  #Add vertical lines to make "two panels"
  clip(0,400,0,1)
  abline(v = 252)
  abline(v = 260)
  
  clip(-1,230,0.1,0.3) #Plot boundaries for thw regression line
  abline(lm(NH4_SRP.data.all$SRP_mgL[NH4_SRP.data.all$AD==0][-c(9:12)] ~ NH4_SRP.data.all$OxicDuration[NH4_SRP.data.all$AD==0][-c(9:12)]),
         lwd = 2)
  
  #Reset boundaries to full plot
  clip(235,410,0.05,0.4)
  #To make it look like two separate plots, I add a fixed value to all AD values to make them appear after the points of the oxic duration part of the plot
  Temp.srp = NH4_SRP.data.all$AD[NH4_SRP.data.all$AD!=0] + 270
  points(NH4_SRP.data.all$SRP_mgL[NH4_SRP.data.all$AD!=0] ~ Temp.srp,
         las = 1,
         ylim = c(0.1,0.36),
         xlim = c(0,400),
         pch = c(15,16,17)[NH4_SRP.data.all$year[NH4_SRP.data.all$AD!=0]],
         col = c("#73BAE6", "#FFD700","#C7144C")[NH4_SRP.data.all$year[NH4_SRP.data.all$AD!=0]],
         cex = 1.5 )
  
  mtext("Anaerobic duration", side = 1, at = 340.5,line = 2, cex = 0.85, font=1)
  mtext("(days)", side = 1, at = 340.5,line = 3, cex = 0.85, font=1)
  clip(271,410,0.2,0.35)
  abline(lm(NH4_SRP.data.all$SRP_mgL[NH4_SRP.data.all$AD!=0] ~ Temp.srp),
         lwd = 2) #slope = 0.88?g/L R8=0.53
  
  axis(side = 1, at = c(270,340, 410), labels = c("0", "70", "140"), cex.axis = 1.4)
  
  
  plot(NH4_SRP.data.all$NH4.N_mgL[NH4_SRP.data.all$AD==0] ~ NH4_SRP.data.all$OxicDuration[NH4_SRP.data.all$AD==0],
       las = 1,
       ylim = c(0,1.1),
       xlim = c(0,400),
       ylab = expression(NH[4]^"+"~(mg~L^-1)),
       xlab = "",
       xaxt = "n",
       pch = c(15,16,17)[NH4_SRP.data.all$year[NH4_SRP.data.all$AD==0]],
       col = c("#73BAE6", "#FFD700","#C7144C")[NH4_SRP.data.all$year[NH4_SRP.data.all$AD==0]],
       cex = 1.5,
       cex.lab = 1.5,
       cex.axis = 1.4)#NH4_SRP.data.all$cex[NH4_SRP.data.all$AD==0])
  axis(side = 1, at = c(0,100,200), labels = c("0", "100", "200"), cex.axis = 1.4)
  mtext("Oxic duration", side = 1, at = 117.5,line = 2, cex = 0.85, font=1)
  mtext("(days)", side = 1, at = 117.5,line = 3, cex = 0.85, font=1)
  mtext("g", side = 3, at = -10)
  mtext("h", side = 3, at = 246)
  curve(coef(NH4.nlm.Oxic.exp)[2]*exp(-coef(NH4.nlm.Oxic.exp)[1]*x), 
        add=T,
        from = 0.01,
        to = 250,
        lwd = 2)
  clip(x1 = 254.2, x2 = 261.8, y1=-1, y2=2)
  abline(h = -0.043, col = "white", lwd=4)
  abline(h = 1.143, col = "white", lwd=4)
  clip(0,400,-1,2)
  abline(v = 254)
  abline(v = 262)
  
  
  clip(235,410,-0.1,1.2)
  points(NH4_SRP.data.all$NH4.N_mgL[NH4_SRP.data.all$AD!=0] ~ Temp.srp,
         las = 1,
         ylim = c(0.1,0.36),
         xlim = c(0,400),
         pch = c(15,16,17)[NH4_SRP.data.all$year[NH4_SRP.data.all$AD!=0]],
         col = c("#73BAE6", "#FFD700","#C7144C")[NH4_SRP.data.all$year[NH4_SRP.data.all$AD!=0]],
         cex = 1.5)
  
  axis(side = 1, at = c(270,340, 410), labels = c("0", "70", "140"), cex.axis = 1.4)
  mtext("Anaerobic duration", side = 1, at = 340.5,line = 2, cex = 0.85, font=1)
  mtext("(days)", side = 1, at = 340.5,line = 3, cex = 0.85, font=1)
  clip(271,410,0,1.0)
  abline(lm(NH4_SRP.data.all$NH4.N_mgL[NH4_SRP.data.all$AD!=0] ~ Temp.srp),
         lwd = 2) #slope = 7.2?g/L R2 = 0.50
    mtext("Lake Arendsee", side = 4, line = 0.2)
  dev.off()
}

```

