CT_dynamics.Rmd

---
title: "CT dynamics"
author: "Jens Daniel Müller"
date:  "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  workflowr::wflow_html:
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: false
editor_options:
  chunk_output_type: console
---


```{r global_options, include = FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
```

```{r packages}
library(tidyverse)
library(patchwork)
library(seacarb)
library(marelac)
library(metR)
library(scico)
library(lubridate)
library(zoo)
library(tibbletime)
library(sp) # check points in polygon

```

```{r ggplot_theme, include = FALSE}
theme_set(theme_bw())
```


# Sensor data

## Data preparation

```{r subsetting criteria}

dep_grid <- 1
max_dep <- 25
max_dep_gap <- 20
max_gap <- 3
surface_dep <- 6

stations_out <- c("PX1", "PX2", "P14", "P13", "P01")
dates_out <- c("180616","180820")
phases_in <- c("down", "up")

```


Profile data are prepared by:

- Ignoring those made on June 16 (pCO2 sensor not in operation)
- Removing HydroC Flush and Zeroing periods
- Selecting only continous downcast periods
- Gridding profiles to `r dep_grid`m depth intervals
- removing grids with pCO2 < 0 µatm (presumably RT correction artefact after zeroing)
- Discarding profiles with `r max_dep_gap` or more observation missing within upper `r max_dep`m
- assigning mean date_time_ID value to all profiles belonging to one cruise
- discarding "coastal" station P01, P13, P14
- Restricting profiles to upper `r max_dep`m

Please note that:

- The label ID representm the start date of the cruise ("YYMMDD"), not the exact mean sampling date


```{r read_prepare_sensor_data}

tm <-
 read_csv(here::here("data/_merged_data_files/response_time",
                      "tm_RT_all.csv"),
               col_types = cols(ID = col_character(),
                                pCO2_analog = col_double(),
                                pCO2_corr = col_double(),
                                Zero = col_character(),
                                Flush = col_character(),
                                mixing = col_character(),
                                Zero_counter = col_integer(),
                                deployment = col_integer(),
                                lon = col_double(),
                                lat = col_double(),
                                pCO2 = col_double()))


# Filter relevant rows and columns
tm_profiles <- tm %>% 
  filter(type == "P",
         Flush == "0",
         Zero == "0",
         !ID %in% dates_out,
         !(station %in% c("PX1", "PX2"))) %>% 
  select(date_time, ID, station, lat, lon, dep, sal, tem, pCO2_corr, pCO2, duration)

stations <- tm_profiles %>% 
  group_by(station) %>% 
  summarise(lat = mean(lat),
            lon = mean(lon)) %>% 
  ungroup()

tm_profiles <- tm_profiles %>% 
  filter(!(station %in% c("P14", "P13", "P01")))

# Assign meta information
tm_profiles <- tm_profiles %>% 
  group_by(ID, station) %>% 
  mutate(duration = as.numeric(date_time - min(date_time))) %>%
  arrange(date_time) %>% 
  ungroup()

meta <- read_csv(here::here("Data/TinaV/Sensor",
                            "Sensor_meta.csv"),
                 col_types = cols(ID = col_character()))

meta <- meta %>% 
    filter(!ID %in% dates_out,
           !(station %in% stations_out))

tm_profiles <- full_join(tm_profiles, meta)
rm(meta)


# creating descriptive variables
tm_profiles <- tm_profiles %>% 
  mutate(phase = "standby",
         phase = if_else(duration >= start & duration < down & !is.na(down) & !is.na(start),
                         "down", phase),
         phase = if_else(duration >= down  & duration < lift & !is.na(lift) & !is.na(down ),
                         "low",  phase),
         phase = if_else(duration >= lift  & duration < up   & !is.na(up  ) & !is.na(lift  ),
                         "mid",  phase),
         phase = if_else(duration >= up    & duration < end  & !is.na(end ) & !is.na(up   ),
                         "up",   phase))

tm_profiles <- tm_profiles %>% 
  select(-c(start, down, lift, up, end, comment, p_type, duration))


# select downcasst only
tm_profiles <- tm_profiles %>% 
  filter(phase %in% phases_in)

#tm_profiles_highres <- tm_profiles

# grid observation to 1m depth intervals
tm_profiles <- tm_profiles %>% 
  mutate(dep_grid = as.numeric(as.character(cut(dep, seq(0,40,1), seq(0.5,39.5,1))))) %>% 
  group_by(ID, station, dep_grid, phase) %>%
  summarise_all("mean", na.rm = TRUE) %>% 
  ungroup() %>% 
  select(-dep, dep=dep_grid)

# Remove zero pCO2 data
tm_profiles <- tm_profiles %>%
  filter(pCO2 >= 0)

# subset complete profiles
profiles_in <- tm_profiles %>% 
  filter(dep < max_dep_gap,
         phase == "down") %>% 
  group_by(ID, station) %>% 
  summarise(nr_na = max_dep_gap/dep_grid - n()) %>% 
  mutate(select = if_else(nr_na < max_gap,
                          "in", "out")) %>% 
  select(-nr_na) %>% 
  ungroup()

tm_profiles <- full_join(tm_profiles, profiles_in)
rm(profiles_in)

```

## pCO~2~ profile overview

```{r pCO2_profiles_overview, fig.cap="Overview pCO~2~ profiles at stations (P02-P12) and cruise dates (ID). y-axis restricted to displayed range.", fig.asp = 1.5}

tm_profiles %>% 
  arrange(date_time) %>% 
  ggplot(aes(pCO2, dep, col=select, linetype=phase))+
  geom_hline(yintercept = 25)+
  geom_path()+
  scale_y_reverse()+
  scale_x_continuous(breaks = c(0, 600), labels = c(0, 600))+
  scale_color_brewer(palette = "Set1", direction = -1)+
  coord_cartesian(xlim = c(0,600))+
  facet_grid(ID~station)

```

## Subset

```{r subset_profiles_downcast_in}

tm_profiles <- tm_profiles %>%
  filter(select == "in",
         phase == "down") %>%
  select(-c(select, phase)) %>% 
  filter(dep < max_dep)

```


```{r remove_subsetting criteria}

rm(dep_grid,max_dep_gap,max_gap,
   stations_out,dates_out,phases_in)

```

## Cruise dates

```{r mean_cruise_dates}
# assign mean date_time stamp

cruise_dates <- tm_profiles %>% 
  group_by(ID) %>% 
  summarise(date_time_ID = mean(date_time),
            date_ID = format(as.Date(date_time_ID), "%b %d")) %>% 
  ungroup()

# inner_join remove P14 observations lacking date_time_ID 
tm_profiles <- inner_join(cruise_dates, tm_profiles)

```



## Station map

```{r subsetting_criteria_fm_and_plot}

# borders for fm data in BloomSail area
poly_lon <- c(18.9, 19.5, 19.5, 18.9)
poly_lat <- c(57.4, 57.2, 57.45, 57.61)

# plot area
lat_lo <- 57.25
lat_hi <- 57.6

lon_lo <- 18.6
lon_hi <- 19.7

```

```{r read_subset_fm_data}

fm <-
  read_csv(here::here("Data/_summarized_data_files",
                       "fm.csv")) %>% 
  filter(lat <= lat_hi, lat >= lat_lo, lon >= lon_lo)

fm <- fm %>% 
  mutate(Area = point.in.polygon(point.x = lon,
                                 point.y = lat,
                                 pol.x = poly_lon,
                                 pol.y = poly_lat),
         Area = as.character(Area),
         Area = if_else(Area == "1", "utilized", "sampled"))

fm %>% 
  filter(Area == "utilized") %>% 
  select(-Area) %>% 
  write_csv(here::here("Data/_summarized_data_files",
                       "fm_bloomsail.csv"))

```


```{r map, fig.cap="Location of stations sampled between the east coast of Gotland and Gotland deep."}


map <- read_csv(here::here("data/Maps","Bathymetry_Gotland_east_small.csv")) %>% 
  filter(lat < lat_hi, lat > lat_lo,
         lon < lon_hi, lon > lon_lo)

map_low_res <- map %>% 
  mutate(lat = cut(lat,
                   breaks = seq(57,58,0.01),
                   labels = seq(57.005,57.995,0.01)),
         lon = cut(lon,
                   breaks = seq(18,22,0.01),
                   labels = seq(18.005,21.995,0.01))) %>% 
  group_by(lat,lon) %>% 
  summarise_all(mean, na.rm=TRUE) %>% 
  ungroup() %>% 
  mutate(lat = as.numeric(as.character(lat)),
         lon = as.numeric(as.character(lon)))

tm_track <- tm %>% 
  arrange(date_time) %>% 
  slice(which(row_number() %% 50 == 1))

ggplot()+
  geom_contour_fill(data = map_low_res,
                    aes(x=lon, y=lat, z=-elev),
                    na.fill = TRUE,
                    breaks = seq(0,300,30))+
  geom_raster(data=map %>% filter(is.na(elev)),
              aes(lon, lat), fill="darkgrey")+
  geom_path(data = tm_track, aes(lon, lat, group=ID, col="sampled"))+
  geom_path(data = fm, aes(lon, lat, group=ID, col=Area))+
  geom_label(data = stations %>% filter(!(station %in% c("P14", "P13", "P01"))),
             aes(lon, lat, label=station, col="utilized"))+
  geom_label(data = stations %>% filter(station %in% c("P14", "P13", "P01")),
             aes(lon, lat, label=station, col="sampled_station"))+
  coord_quickmap(expand = 0, ylim = c(lat_lo+0.01, lat_hi-0.01))+
  labs(x="Longitude (°E)", y="Latitude (°N)")+
  scale_fill_gradient(low = "lightsteelblue1", high = "dodgerblue4",
                   name="Depth (m)", breaks = seq(0,180,30),
                   guide = "colorstrip")+
  scale_color_manual(values = c("white", "darkgrey", "orangered"), guide=FALSE)

ggsave(here::here("output/Plots/Figures_publication/article", "station_map.pdf"),
       width = 220, height = 150, dpi = 300, units = "mm")

rm(map, map_low_res,
   lat_hi, lat_lo, lon_hi, lon_lo, 
   poly_lat, poly_lon,
   fm, tm_track)

rm(tm)

```


## Data coverage

```{r data_coverage, fig.cap="Spatio-temporal data coverage, indicated as station visits over time."}

cover <- tm_profiles %>% 
  group_by(ID, station) %>% 
  summarise(date = mean(date_time),
            date_time_ID = mean(date_time_ID)) %>% 
  ungroup()

cover %>% 
  ggplot(aes(date, station, fill=ID))+
  geom_vline(aes(xintercept = date_time_ID, col=ID))+
  geom_point(shape=21)+
  scale_fill_viridis_d(labels = cruise_dates$date_ID,
                       name = "Mean cruise date")+
  scale_color_viridis_d(labels = cruise_dates$date_ID,
                       name = "Mean cruise date")+
    scale_x_datetime(date_breaks = "week",
                   date_labels = "%b %d")+
  theme(axis.title.x = element_blank())

ggsave(here::here("output/Plots/Figures_publication/article", "data_coverage.pdf"),
       width = 130, height = 65, dpi = 300, units = "mm")

rm(cover)

```

# Bottle C~T~ and A~T~

At stations P07 and P10 discrete samples for lab measurmentm of C~T~ and A~T~ were collected. Please note that - in contrast to the pCO~2~ profiles - samples were taken on June 16, but removed here for harmonization of results.

```{r read_bottle_CO2}

tb <-
  read_csv(here::here("Data/_summarized_data_files", "tb.csv"),
           col_types = cols(ID = col_character()))

tb <- tb %>% 
  filter(station %in% c("P07", "P10"),
         dep <= max_dep) %>% 
  mutate(ID = if_else(ID == "180722", "180723", ID))

tb <- inner_join(tb, cruise_dates)

```

## Mean alkalinity

In order to derive C~T~ from measured pCO~2~ profiles, the alkalinity mean + sd in the upper `r max_dep`m and both stations was calculated as:

```{r bottle_AT}

AT_mean <- tb %>% 
  filter(dep <= max_dep) %>% 
  summarise(AT = mean(AT, na.rm = TRUE)) %>%
  pull()

AT_mean

AT_sd <- tb %>% 
  filter(dep <= 20) %>% 
  summarise(AT = sd(AT, na.rm = TRUE)) %>%
  pull()

AT_sd

```

Likewise, the mean salinity amounts to:

```{r bottle_sal}

sal_mean <- tb %>% 
  filter(dep <= 20) %>% 
  summarise(sal = mean(sal, na.rm = TRUE)) %>%
  pull()

sal_mean

```

```{r write_fixed_values}

bind_cols(start = min(tm_profiles$date_time), 
          end = max(tm_profiles$date_time),
          AT = AT_mean,
          AT_sd = AT_sd,
          sal = sal_mean) %>% 
  write_csv(here::here("Data/_summarized_data_files", "tb_fix.csv"))

rm(AT_sd)
```

## nC~T~ calculation

The alkalinity-normalized C~T~, nC~T~, was calculated.

```{r nCT_bottle_calculation}

tb <- tb %>% 
  mutate(nCT = CT/AT * AT_mean)

```

## Vertical profiles

```{r bottle_CO2_profiles}

tb_long <- tb %>% 
  pivot_longer(c(sal:AT, nCT), names_to = "var", values_to = "value")

tb_long %>% 
  ggplot(aes(value, dep))+
  geom_path(aes(col=ID))+
  geom_point(aes(fill=ID), shape=21)+
  scale_y_reverse()+
  scale_fill_viridis_d(labels = cruise_dates$date_ID)+
  scale_color_viridis_d(labels = cruise_dates$date_ID)+
  facet_grid(station~var, scales = "free_x")+
  theme(legend.position = "bottom",
        legend.title = element_blank())

```


## Surface time series

```{r plot_bottle_CO2_surface_timeseries, fig.cap="Time series of bottle data. Shown are mean values of samples collected at water depths < 10m (usually collected at 0 and 5 m).", fig.asp=1.1}

tb_surface <- tb_long %>% 
  filter(dep < surface_dep) %>% 
  group_by(ID, date_time_ID, var, station) %>% 
  summarise(value = mean(value, na.rm = TRUE)) %>% 
  ungroup()

tb_surface_station_mean <- tb_long %>% 
  filter(dep < surface_dep) %>% 
  group_by(ID, date_time_ID, var) %>% 
  summarise(value_mean = mean(value, na.rm = TRUE),
            value_sd = sd(value, na.rm = TRUE)) %>% 
  ungroup()

tb_long %>%
  filter(dep<11) %>% 
  ggplot()+
  geom_line(data = tb_surface, aes(date_time_ID, value, col="Individual"))+
  geom_line(data = tb_surface_station_mean, aes(date_time_ID, value_mean, col="Both (mean)"))+
  geom_point(aes(date_time_ID, value, fill=dep), shape=21)+
  scale_fill_scico(palette = "oslo", direction = -1, name="Depth (m)")+
  scale_color_brewer(palette = "Set1", name="Station surface mean")+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  facet_grid(var~station, scales = "free_y")+
  labs(x="Mean transect date")


rm(tb_long, tb_surface, tb)
```


Important notes:
- nC~T~ drop and temporal patterns agree well with those found in the nC~T~ time series derived from pCO~2~ measurements (below).


# nC~T~ profiles

## Calculation from pCO~2~

Alkalinity normalized C~T~ (nC~T~) profiles were calculated from sensor pCO~2~ and T profiles, and constant salinity and alkalinity values. Note that the impact of fixed vs. measured salinity has only a negligible impact on nC~T~ profiles.

```{r CT_calculation}

tm_profiles <- tm_profiles %>% 
  mutate(nCT = carb(24, var1=pCO2, var2=AT_mean*1e-6,
                   S=sal_mean, T=tem, P=dep/10, k1k2="m10", kf="dg", ks="d",
                   gas="insitu")[,16]*1e6)

rm(sal_mean)

tm_profiles %>% 
  write_csv(here::here("Data/_merged_data_files/CT_dynamics", "tm_profiles.csv"))

```

## Plot all profiles

```{r profiles_all_by_date}

tm_profiles <-  tm_profiles %>% 
  arrange(date_time_ID)

p_tem <- 
  tm_profiles %>% 
  ggplot(aes(tem, dep, col=ID, group = interaction(station, ID)))+
  geom_path()+
  scale_y_reverse(expand = c(0,0))+
  labs(x = "Temperature (\u00B0C)",
       y = "Depth (m)")+
  scale_color_viridis_d(guide = FALSE)

p_pCO2 <- 
  tm_profiles %>% 
  ggplot(aes(pCO2, dep, col=ID, group = interaction(station, ID)))+
  geom_path()+
  scale_y_reverse(expand = c(0,0))+
  labs(x = expression(pCO[2]~(µatm)))+
  scale_color_viridis_d(guide = FALSE)+
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank())

p_nCT <- 
  tm_profiles %>% 
  ggplot(aes(nCT, dep, col=ID, group = interaction(station, ID)))+
  geom_path()+
  scale_y_reverse(expand = c(0,0))+
  labs(x = expression(nC[T]~(µmol~kg^{-1})))+
  scale_color_viridis_d(labels = cruise_dates$date_ID)+
  theme(legend.title = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank())

p_tem + p_pCO2 + p_nCT

ggsave(here::here("output/Plots/Figures_publication/article", "profiles_all.pdf"),
       width = 180, height = 80, dpi = 300, units = "mm")

rm( p_tem, p_pCO2, p_nCT)

```



## Mean profiles

Mean vertical profiles were calculated for each cruise day (ID).

```{r calculate_profile_mean_sd}

tm_profiles_ID_mean <- tm_profiles %>% 
  select(-c(station,lat, lon, pCO2_corr, date_time)) %>% 
  group_by(ID, date_time_ID, dep) %>% 
  summarise_all(list(mean), na.rm=TRUE) %>% 
  ungroup()

tm_profiles_ID_sd <- tm_profiles %>% 
  select(-c(station,lat, lon, pCO2_corr, date_time)) %>% 
  group_by(ID, date_time_ID, dep) %>% 
  summarise_all(list(sd), na.rm=TRUE) %>% 
  ungroup()

tm_profiles_ID_sd_long <- tm_profiles_ID_sd %>% 
  pivot_longer(sal:nCT, names_to = "var", values_to = "sd")
  
tm_profiles_ID_mean_long <- tm_profiles_ID_mean %>% 
  pivot_longer(sal:nCT, names_to = "var", values_to = "value")
  
tm_profiles_ID_long <- inner_join(tm_profiles_ID_mean_long, tm_profiles_ID_sd_long)

tm_profiles_ID_mean %>% 
  write_csv(here::here("Data/_merged_data_files/CT_dynamics", "tm_profiles_ID.csv"))

rm(tm_profiles_ID_sd_long, tm_profiles_ID_sd, tm_profiles_ID_mean_long, tm_profiles_ID_mean)

```


```{r nCT_profiles_all, fig.cap="Mean vertical profiles per cruise day across all stations."}

tm_profiles_ID_long %>% 
  ggplot(aes(value, dep, col=ID))+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  scale_color_viridis_d()+
  facet_wrap(~var, scales = "free_x")
  
```

```{r nCT_profiles_individual, fig.cap="Mean vertical profiles per cruise day across all stations plotted indivdually. Ribbons indicate the standard deviation observed across all profiles at each depth and transect.", fig.asp=1.5}

all <- tm_profiles_ID_long %>% 
  filter(var %in% c("nCT", "tem")) %>% 
  rename(group = ID)

tm_profiles_ID_long %>% 
  filter(var %in% c("nCT", "tem")) %>% 
  ggplot()+
  geom_path(data=all, aes(value, dep, group=group))+
  geom_ribbon(aes(xmin = value-sd, xmax=value+sd, y=dep, fill=ID), alpha=0.5)+
  geom_path(aes(value, dep, col=ID))+
  scale_y_reverse()+
  scale_color_viridis_d()+
  scale_fill_viridis_d()+
  facet_grid(ID~var, scales = "free_x")

rm(all)

```

Important notes:

- the standard deviation of C~T~ in the upper 10m increases on June 30.

## Individual profiles

C~T~, pCO~2~, S, and T profiles were plotted individually 
[pdf here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plotm/CT_dynamics/tm_profiles_pCO2_tem_sal_CT.pdf){target="_blank"} 
and grouped by ID 
[pdf here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plotm/CT_dynamics/tm_profiles_ID_pCO2_tem_sal_CT.pdf){target="_blank"}. The later gives an idea of the differences between stations at one point in time.

```{r plot_all_discretized_profiles_individual, eval=FALSE}

# tm_profiles_highres <- tm_profiles_highres %>% 
#   filter(phase == "down")

pdf(file=here::here("output/Plots/CT_dynamics",
    "tm_profiles_pCO2_tem_sal_CT.pdf"), onefile = TRUE, width = 9, height = 5)

for(i_ID in unique(tm_profiles$ID)){
  for(i_station in unique(tm_profiles$station)){

    if (nrow(tm_profiles %>% filter(ID == i_ID, station == i_station)) > 0){
      
      
      # i_ID      <-      unique(tm_profiles$ID)[1]
      # i_station <- unique(tm_profiles$station)[1]

      p_pCO2 <- 
        tm_profiles %>%
        arrange(date_time) %>% 
        filter(ID == i_ID,
               station == i_station) %>%
        ggplot(aes(pCO2, dep, col="grid_RT"))+
        # geom_point(data = tm_profiles_highres %>% 
        #              arrange(date_time) %>% 
        #              filter(ID == i_ID, station == i_station),
        #            aes(pCO2_corr, dep, col="raw"))+
        # geom_point(data = tm_profiles_highres %>% 
        #              arrange(date_time) %>% 
        #              filter(ID == i_ID, station == i_station),
        #            aes(pCO2, dep, col="raw_RT"))+
        geom_point(aes(pCO2_corr, dep, col="grid"))+
        geom_point()+
        geom_path()+
        scale_y_reverse()+
        scale_color_brewer(palette = "Set1")+
        labs(y="Depth [m]", x="pCO2 [µatm]", title = str_c(i_ID," | ",i_station))+
        coord_cartesian(xlim = c(0,200), ylim = c(30,0))+
        theme_bw()+
        theme(legend.position = "left")
      
      p_tem <- 
        tm_profiles %>%
        arrange(date_time) %>% 
        filter(ID == i_ID,
               station == i_station) %>%
        ggplot(aes(tem, dep))+
        geom_point()+
        geom_path()+
        scale_y_reverse()+
        labs(y="Depth [m]", x="Tem [°C]")+
        coord_cartesian(xlim = c(14,26), ylim = c(30,0))+
        theme_bw()
      
      p_sal <- 
        tm_profiles %>%
        arrange(date_time) %>% 
        filter(ID == i_ID,
               station == i_station) %>%
        ggplot(aes(sal, dep))+
        geom_point()+
        geom_path()+
        scale_y_reverse()+
        labs(y="Depth [m]", x="Tem [°C]")+
        coord_cartesian(xlim = c(6.5,7.5), ylim = c(30,0))+
        theme_bw()
      
      p_nCT <- 
        tm_profiles %>%
        arrange(date_time) %>% 
        filter(ID == i_ID,
               station == i_station) %>%
        ggplot(aes(nCT, dep))+
        geom_point()+
        geom_path()+
        scale_y_reverse()+
        labs(y="Depth [m]", x="nCT* [µmol/kg]")+
        coord_cartesian(xlim = c(1400,1700), ylim = c(30,0))+
        theme_bw()
      

      print(
            p_pCO2 + p_tem + p_sal + p_nCT
            )
      
      rm(p_pCO2, p_sal, p_tem, p_nCT)

      
    }
  }
}

dev.off()

rm(i_ID, i_station, tm_profiles_highres)

```

```{r plot_all_discretized_profiles_by_ID, eval=FALSE}

tm_profiles_long <- tm_profiles %>%
  select(-c(lat, lon, pCO2_corr)) %>% 
  pivot_longer(sal:nCT, values_to = "value", names_to = "var")


pdf(file=here::here("output/Plots/CT_dynamics",
    "tm_profiles_ID_pCO2_tem_sal_CT.pdf"), onefile = TRUE, width = 9, height = 5)

for(i_ID in unique(tm_profiles$ID)){

  #i_ID <- unique(tm_profiles$ID)[1]
  
  sub_tm_profiles_long <- tm_profiles_long %>% 
        arrange(date_time) %>% 
        filter(ID == i_ID)
 
  print(
  
  sub_tm_profiles_long %>% 
    ggplot()+
    geom_path(data = tm_profiles_long,
              aes(value, dep, group=interaction(station, ID)), col="grey")+
    geom_path(aes(value, dep, col=station))+
    scale_y_reverse()+
    labs(y="Depth [m]", title = str_c("ID: ", i_ID))+
    theme_bw()+
    facet_wrap(~var, scales = "free_x")
   
  )  
  rm(sub_tm_profiles_long)
}

dev.off()

rm(i_ID, tm_profiles_long)

```



## Profiles of incremental changes

Changes of seawater vars at each depth are calculated from one cruise day to the next and divided by the number of days inbetween.

```{r calculate_profiles_incremental_cumulative_changes}

tm_profiles_ID_long <- tm_profiles_ID_long %>%
    group_by(var, dep) %>%
    arrange(date_time_ID) %>%
    mutate(date_time_ID_diff = as.numeric(date_time_ID - lag(date_time_ID)),
           date_time_ID_ref  = date_time_ID - (date_time_ID - lag(date_time_ID))/2,
           value_diff = value     - lag(value, default = first(value)),
           value_diff_daily = value_diff / date_time_ID_diff,
           value_cum = cumsum(value_diff)) %>% 
  ungroup()
```

```{r incremental_changes_profiles}

tm_profiles_ID_long %>% 
  arrange(dep) %>% 
  ggplot(aes(value_diff_daily, dep, col=ID))+
  geom_vline(xintercept = 0)+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  scale_color_viridis_d()+
  facet_wrap(~var, scales = "free_x")+
  labs(x="Change of value inbetween cruises per day")
  
```

## Profiles of cumulative changes

Cumulative changes of seawater vars were calculated at each depth relative to the first cruise day on July 5.

```{r cumulative_changes_profiles}

tm_profiles_ID_long %>% 
  arrange(dep) %>% 
  ggplot(aes(value_cum, dep, col=ID))+
  geom_vline(xintercept = 0)+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  scale_color_viridis_d()+
  facet_wrap(~var, scales = "free_x")+
  labs(x="Cumulative change of value")

```


Important notes:

- Salinity in the upper 10m decreases by >0.1 on June 30, and returns to average conditions already on Aug 02.

Cumulative positive and negative changes of seawater vars were calculated separately at each depth relative to the first cruise day on July 5.


```{r calculate_profiles_directional_changes}

tm_profiles_ID_long <- tm_profiles_ID_long %>% 
  mutate(sign = if_else(value_diff < 0, "neg", "pos")) %>% 
  group_by(var, dep, sign) %>%
  arrange(date_time_ID) %>%
  mutate(value_cum_sign = cumsum(value_diff)) %>% 
  ungroup()

```


```{r cumulative_directional_changes_profiles}

tm_profiles_ID_long %>% 
  arrange(dep) %>% 
  ggplot(aes(value_cum_sign, dep, col=ID))+
  geom_vline(xintercept = 0)+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  scale_color_viridis_d()+
  scale_fill_viridis_d()+
  facet_wrap(~interaction(sign, var), scales = "free_x", ncol=4)+
  labs(x="Cumulative directional change of value")


# tm_profiles_ID_long %>%
#   write_csv(here::here("Data/_merged_data_files/CT_dynamics", "tm_profiles_ID_long_cum.csv"))

```

# Timeseries

## Timeseries depth intervals

Mean seawater parameters were calculated for 5m depth intervals.

```{r observed_timeseries, fig.asp=1.3}

tm_profiles_ID_long_grid <- tm_profiles_ID_long %>% 
  mutate(dep = cut(dep, seq(0,30,5))) %>% 
  group_by(ID, date_time_ID, dep, var)  %>% 
  summarise_all(list(mean), na.rm=TRUE) 

tm_profiles_ID_long_grid %>% 
  ggplot(aes(date_time_ID, value, col=as.factor(dep)))+
  geom_path()+
  geom_point()+
  scale_color_viridis_d(name="Depth (m)")+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  facet_wrap(~var, scales = "free_y", ncol=1)+
  theme(axis.title.x = element_blank())

rm(tm_profiles_ID_long_grid)
```


## Hovmoeller plots

### Absolute values

```{r plot_hovmoeller_absolut, fig.cap="Hovmoeller plotm of absolute changes in C~T~ and temperature."}

bin_nCT <- 30

p_nCT_hov <- tm_profiles_ID_long %>% 
  filter(var == "nCT") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value),
                    breaks = MakeBreaks(bin_nCT),
                    col="black")+
  geom_point(aes(x=date_time_ID, y=c(24.5)), size=3, shape=24, fill="white")+
  scale_fill_scico(breaks = MakeBreaks(bin_nCT),
                   guide = "colorstrip",
                   name="nCT (µmol/kg)",
                   palette = "davos", direction = -1)+
  scale_y_reverse()+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  theme_bw()+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        legend.position = "left")

bin_Tem <- 2

p_tem_hov <- tm_profiles_ID_long %>% 
  filter(var == "tem") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value),
                    breaks = MakeBreaks(bin_Tem),
                    col="black")+
  geom_point(aes(x=date_time_ID, y=c(24.5)), size=3, shape=24, fill="white")+
  scale_fill_viridis_c(breaks = MakeBreaks(bin_Tem),
                       guide = "colorstrip",
                       name="Tem (°C)",
                       option = "inferno")+
  scale_y_reverse()+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        legend.position = "left")

p_nCT_hov / p_tem_hov

rm(p_nCT_hov, bin_nCT, p_tem_hov, bin_Tem)

```


```{r pub_fig_nCT_tem_hov_abs_profiles_cum, include=FALSE}

bin_nCT <- 30

p_nCT_hov <- tm_profiles_ID_long %>% 
  filter(var == "nCT") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value),
                    breaks = MakeBreaks(bin_nCT),
                    col="black")+
  geom_vline(aes(xintercept = date_time_ID), col="white", linetype = "1f")+
  scale_fill_scico(breaks = MakeBreaks(bin_nCT),
                   guide = "colorstrip",
                   name= expression(nC[T]~(µmol~kg^{-1})),
                   palette = "davos", direction = -1)+
  scale_y_reverse(sec.axis = dup_axis())+
  scale_x_datetime(breaks = "week",
                   date_labels = "%b %d")+
  theme_bw()+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        axis.title.y.right = element_blank(),
        axis.text.y.right = element_blank(),
        legend.position = "left")

bin_Tem <- 2

p_tem_hov <- tm_profiles_ID_long %>% 
  filter(var == "tem") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value),
                    breaks = MakeBreaks(bin_Tem),
                    col="black")+
  geom_vline(aes(xintercept = date_time_ID), col="white", linetype = "1f")+
  scale_fill_viridis_c(breaks = MakeBreaks(bin_Tem),
                       guide = "colorstrip",
                       name="Temperature (\u00B0C)",
                       option = "inferno")+
  scale_y_reverse(sec.axis = dup_axis())+
  scale_x_datetime(breaks = "week",
                   date_labels = "%b %d")+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        axis.title.y.right = element_blank(),
        axis.text.y.right = element_blank(),
        legend.position = "left")

  
p_nCT_ID_cum <- 
  tm_profiles_ID_long %>%
  filter(var == "nCT") %>% 
  ggplot(aes(value_cum, dep, col=ID, fill=ID))+
  geom_path()+
  geom_point(shape=21, col="black")+
  scale_y_reverse(expand = c(0,0),
                  position = "right",
                  sec.axis = dup_axis())+
  labs(x = expression(nC[T]~cumulative~changes~(µmol~kg^{-1})),
       y = "Depth (m)")+
  scale_color_viridis_d(labels = cruise_dates$date_ID)+
  scale_fill_viridis_d(labels = cruise_dates$date_ID)+
  theme(legend.position = "top",
        axis.title.y.left = element_blank(),
        axis.text.y.left = element_blank(),
        legend.title = element_blank())


p_tem_ID_cum <- 
  tm_profiles_ID_long %>%
  filter(var == "tem") %>% 
  ggplot(aes(value_cum, dep, col=ID, fill=ID))+
  geom_path()+
  geom_point(shape=21, col="black")+
  scale_y_reverse(expand = c(0,0),
                  position = "right",
                  sec.axis = dup_axis())+
  labs(x = "Temperature cumulative changes (\u00B0C)",
       y = "Depth (m)")+
  scale_color_viridis_d(labels = cruise_dates$date_ID, guide = FALSE)+
  scale_fill_viridis_d(labels = cruise_dates$date_ID, guide = FALSE)+
  theme(axis.title.y.left = element_blank(),
        axis.text.y.left = element_blank())


p_nCT_hov + p_nCT_ID_cum +
p_tem_hov + p_tem_ID_cum +
  plot_layout(widths = c(2, 1))

ggsave(here::here("output/Plots/Figures_publication/article", "Hov_abs_profiles_cum.pdf"),
       width = 300, height = 190, dpi = 300, units = "mm")

rm(p_nCT_hov, bin_nCT, p_tem_hov, bin_Tem, p_nCT_ID_cum, p_tem_ID_cum)

```


### Incremental changes

```{r plot_hovmoeller_daily, fig.cap="Hovmoeller plotm of daily changes in C~T~ and temperature. Note that calculated  value of change (in contrast to absolute and cumulative values) are referred to the mean dates inbetween cruise, and are not extrapolated to the full observational period."}

bin_nCT <- 2.5

nCT_hov <- tm_profiles_ID_long %>% 
  filter(var == "nCT") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID_ref, y=dep, z=value_diff_daily),
                    breaks = MakeBreaks(bin_nCT),
                    col="black")+
  geom_point(aes(x=date_time_ID, y=c(24.5)), size=3, shape=24, fill="white")+
  scale_fill_divergent(breaks = MakeBreaks(bin_nCT),
                       guide = "colorstrip",
                       name="nCT (µmol/kg)")+
  scale_y_reverse()+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  theme_bw()+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

bin_Tem <- 0.1

Tem_hov <- tm_profiles_ID_long %>% 
  filter(var == "tem") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID_ref, y=dep, z=value_diff_daily),
                    breaks = MakeBreaks(bin_Tem),
                    col="black")+
  geom_point(aes(x=date_time_ID, y=c(24.5)), size=3, shape=24, fill="white")+
  scale_fill_divergent(breaks = MakeBreaks(bin_Tem),
                       guide = "colorstrip",
                       name="Tem (°C)")+
  scale_y_reverse()+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  theme_bw()+
  labs(x="",y="Depth (m)")+
  coord_cartesian(expand = 0)

nCT_hov / Tem_hov

rm(nCT_hov, bin_nCT, Tem_hov, bin_Tem)

```


### Cumulative changes

```{r plot_hovmoeller_cumulative, fig.cap="Hovmoeller plotm of cumulative changes in C~T~ and temperature."}

bin_nCT <- 20

nCT_hov <- tm_profiles_ID_long %>% 
  filter(var == "nCT") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value_cum),
                    breaks = MakeBreaks(bin_nCT),
                    col="black")+
  geom_point(aes(x=date_time_ID, y=c(24.5)), size=3, shape=24, fill="white")+
  scale_fill_divergent(breaks = MakeBreaks(bin_nCT),
                       guide = "colorstrip",
                       name="nCT (µmol/kg)")+
  scale_y_reverse()+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  theme_bw()+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

bin_Tem <- 2

Tem_hov <- tm_profiles_ID_long %>% 
  filter(var == "tem") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value_cum),
                    breaks = MakeBreaks(bin_Tem),
                    col="black")+
  geom_point(aes(x=date_time_ID, y=c(24.5)), size=3, shape=24, fill="white")+
  scale_fill_divergent(breaks = MakeBreaks(bin_Tem),
                       guide = "colorstrip",
                       name="Tem (°C)")+
  scale_y_reverse()+
  scale_x_datetime(breaks = "week", date_labels = "%d %b")+
  theme_bw()+
  labs(x="",y="Depth (m)")+
  coord_cartesian(expand = 0)

nCT_hov / Tem_hov

rm(nCT_hov, bin_nCT, Tem_hov, bin_Tem)

```



# Depth-integration C~T~

A critical first step for the determination of net community production (NCP) is the integration of observed changes in C~T~ over depth to derive iC~T~. Two approaches were tested:

- Integration of changes in nCT over a predefined, fixed water depth
- Integration of changes in nCT over a mixed layer depth (MLD)

Both aproaches deliver depth-integrated, incremental changes of C~T~ inbetween cruise dates. Those were summed up to derive a trajectory of cummulative iC~T~ changes.

## Fixed depths approach

Incremental and cumulative C~T~ changes inbetween cruise dates were integrated across the water colums down to predefined depth limitm. This was done separately for observed positive/negative changes in C~T~, as well as for the total observed changes.

```{r integrated_inCT_fixed_depth, fig.asp=1.1}


inCT_grid_sign <- tm_profiles_ID_long %>% 
  select(ID, date_time_ID, date_time_ID_ref) %>% 
  unique() %>% 
  expand_grid(sign = c("pos", "neg"))

inCT_grid_total <- tm_profiles_ID_long %>% 
  select(ID, date_time_ID, date_time_ID_ref) %>% 
  unique() %>% 
  expand_grid(sign = c("total"))

# dep_i <- 10
#rm(inCT, dep_i)

for (dep_i in seq(9,13,1)) {


inCT_sign_temp <- tm_profiles_ID_long %>% 
  filter(var == "nCT", dep < dep_i) %>% 
  mutate(sign = if_else(ID == "180705" & dep == 0.5, "neg", sign)) %>% 
  group_by(ID, date_time_ID, date_time_ID_ref, sign) %>% 
  summarise(nCT_i_diff = sum(value_diff)/1000) %>% 
  ungroup()

inCT_sign_temp <- inCT_sign_temp %>% 
  group_by(sign) %>%
  arrange(date_time_ID) %>% 
  mutate(nCT_i_cum = cumsum(nCT_i_diff)) %>% 
  ungroup()

inCT_sign_temp <- full_join(inCT_sign_temp, inCT_grid_sign) %>% 
  arrange(sign, date_time_ID) %>% 
  fill(nCT_i_cum)


inCT_total_temp <- tm_profiles_ID_long %>% 
  filter(var == "nCT", dep < dep_i) %>% 
  group_by(ID, date_time_ID, date_time_ID_ref) %>% 
  summarise(nCT_i_diff = sum(value_diff)/1000) %>% 
  ungroup()

inCT_total_temp <- inCT_total_temp %>% 
  arrange(date_time_ID) %>% 
  mutate(nCT_i_cum = cumsum(nCT_i_diff)) %>% 
  ungroup() %>% 
  mutate(sign = "total")

inCT_total_temp <- full_join(inCT_total_temp, inCT_grid_total) %>% 
  arrange(sign, date_time_ID) %>% 
  fill(nCT_i_cum)

inCT_temp <- bind_rows(inCT_sign_temp, inCT_total_temp) %>% 
    mutate(dep_i = dep_i)


if (exists("inCT")) {
  inCT <- bind_rows(inCT, inCT_temp)
  } else {inCT <- inCT_temp}

rm(inCT_temp, inCT_sign_temp, inCT_total_temp)

}

rm(inCT_grid_sign, inCT_grid_total)

inCT <- inCT %>% 
  mutate(dep_i = as.factor(dep_i))

inCT %>% 
  ggplot()+
  geom_point(data = cruise_dates, aes(date_time_ID, 0), shape=21)+
  geom_col(aes(date_time_ID_ref, nCT_i_diff, fill=dep_i),
           position = "dodge", alpha=0.3)+
  geom_line(aes(date_time_ID, nCT_i_cum, col=dep_i))+
  scale_color_viridis_d(name="Depth limit (m)")+
  scale_fill_viridis_d(name="Depth limit (m)")+
  labs(y="inCT [mol/m2]", x="")+
  facet_grid(sign~., scales = "free_y", space = "free_y")+
  theme_bw()

inCT_fixed_dep <- inCT
rm(inCT, dep_i)
# inCT %>%
#   write_csv(here::here("Data/_merged_data_files", "inCT_dep_limitm.csv"))

```

## MLD approach

As an alternative to fixed depth levels, vertical integration as low as the mixed layer depth was tested.

### Density Calculation

Seawater density Rho was determined from S, T, and p according to TEOS-10.

```{r calulcate_seawater_density}

tm_profiles <- tm_profiles %>% 
  mutate(rho = swSigma(salinity = sal, temperature = tem, pressure = dep/10))

```

### Density profiles

```{r plot_hydrographical_profiles_all, fig.cap="Mean vertical profiles per cruise day across all stations."}

tm_profiles_ID_mean_hydro <- tm_profiles %>% 
  select(-c(station,lat, lon, pCO2_corr, pCO2, nCT, date_time)) %>% 
  group_by(ID, date_time_ID, date_ID, dep) %>% 
  summarise_all(list(mean), na.rm=TRUE) %>% 
  ungroup()

tm_profiles_ID_sd_hydro <- tm_profiles %>% 
  select(-c(station,lat, lon, pCO2_corr, pCO2, nCT, date_time)) %>% 
  group_by(ID, date_time_ID, date_ID, dep) %>% 
  summarise_all(list(sd), na.rm=TRUE) %>% 
  ungroup()


tm_profiles_ID_sd_hydro_long <- tm_profiles_ID_sd_hydro %>% 
  pivot_longer(sal:rho, names_to = "var", values_to = "sd")
  
tm_profiles_ID_mean_hydro_long <- tm_profiles_ID_mean_hydro %>% 
  pivot_longer(sal:rho, names_to = "var", values_to = "value")
  
tm_profiles_ID_hydro_long <- inner_join(tm_profiles_ID_mean_hydro_long, tm_profiles_ID_sd_hydro_long)
tm_profiles_ID_hydro <- tm_profiles_ID_mean_hydro

rm(tm_profiles_ID_mean_hydro_long,
   tm_profiles_ID_mean_hydro,
   tm_profiles_ID_sd_hydro_long,
   tm_profiles_ID_sd_hydro)

tm_profiles_ID_hydro_long %>% 
  ggplot(aes(value, dep, col=ID))+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  scale_color_viridis_d()+
  facet_wrap(~var, scales = "free_x")
  
```

### MLD calculation

Mixed layer depth (MLD) was determined based on the difference between density at the surface and at depth, for a range of density criteria

```{r calculate_MLD}
# density criterion

tm_profiles_ID_hydro <- expand_grid(tm_profiles_ID_hydro, rho_lim = c(0.1,0.2,0.5))

MLD <- tm_profiles_ID_hydro  %>% 
  arrange(dep) %>% 
  group_by(ID, date_time_ID, rho_lim) %>% 
  mutate(d_rho = rho - first(rho)) %>% 
  filter(d_rho > rho_lim) %>% 
  summarise(MLD = min(dep)) %>% 
  ungroup()

```

### Daily density profiles

```{r plot_density_profiles, fig.cap="Overview density profiles at stations (P01-P14) and cruise dates (ID). Horizontal lines indicate determined MLD", fig.asp=1.2}

tm_profiles_ID_hydro <- 
  full_join(tm_profiles_ID_hydro, MLD)

tm_profiles_ID_hydro %>% 
  arrange(dep) %>% 
  ggplot(aes(rho, dep))+
  geom_hline(aes(yintercept = MLD, col=as.factor(rho_lim)))+
  geom_path()+
  scale_y_reverse()+
  scale_color_brewer(palette = "Set1", name= "Rho limit")+
  facet_wrap(~ID)+
  theme_bw()

```

### MLD timeseries

```{r plot_MLD_time_series, fig.asp=0.4}

MLD %>% 
  ggplot(aes(date_time_ID, MLD, col=as.factor(rho_lim)))+
  geom_hline(yintercept = 0)+
  geom_point()+
  geom_path()+
  scale_color_brewer(palette = "Set1", name= "Rho limit")+
  scale_y_reverse()+
  labs(x="")

```


### inCT calculation

```{r integrated_inCT_MLD}

# inCT_grid_sign <- tm_profiles_ID_long %>% 
#   seleCT(ID, date_time_ID, date_time_ID_ref) %>% 
#   unique() %>% 
#   expand_grid(sign = c("pos", "neg"))
# 
# inCT_grid_total <- tm_profiles_ID_long %>% 
#   seleCT(ID, date_time_ID, date_time_ID_ref) %>% 
#   unique() %>% 
#   expand_grid(sign = c("total"))


inCT <- tm_profiles_ID_long %>% 
  filter(var == "nCT")

inCT <- full_join(inCT, MLD)

# inCT %>%
#     write_csv(here::here("Data/_merged_data_files", "tm_profiles_ID_long_cum_MLD.csv"))

inCT <- inCT %>% 
  filter(dep <= MLD)

inCT <- inCT %>% 
  group_by(ID, date_time_ID, date_time_ID_ref, rho_lim) %>% 
  summarise(nCT_i_diff = sum(value_diff)/1000) %>% 
  ungroup()

inCT <- inCT %>% 
  group_by(rho_lim) %>% 
  arrange(date_time_ID) %>% 
  mutate(nCT_i_cum = cumsum(nCT_i_diff)) %>% 
  ungroup()

# inCT_total_temp <- full_join(inCT_total_temp, inCT_grid_total) %>% 
#   arrange(sign, date_time_ID) %>% 
#   fill(nCT_i_cum)

inCT <- inCT %>% 
  mutate(rho_lim = as.factor(rho_lim))

inCT %>% 
  ggplot()+
  geom_point(data = cruise_dates, aes(date_time_ID, 0), shape=21)+
  geom_col(aes(date_time_ID_ref, nCT_i_diff, fill=rho_lim),
           position = "dodge", alpha=0.3)+
  geom_line(aes(date_time_ID, nCT_i_cum, col=rho_lim))+
  scale_color_viridis_d(name="Rho limit")+
  scale_fill_viridis_d(name="Rho limit")+
  labs(y="inCT [mol/m2]", x="")+
  #facet_grid(sign~., scales = "free_y", space = "free_y")+
  theme_bw()

inCT_MLD <- inCT
rm(inCT, MLD)

# inCT %>%
#   write_csv(here::here("Data/_merged_data_files", "inCT_dep_limitm.csv"))

```

## Comparison of approaches

In the following, all cummulative iC~T~ trajectories are displayed. Highlighted are those obtained for the fixed depth approach with 10 m limit, and the MLD approach with a high density threshold of 0.5 kg/m3.

```{r comparing_inCT_estimates}

inCT <- full_join(inCT_fixed_dep, inCT_MLD)

inCT <- inCT %>% 
  mutate(group = paste(as.character(sign), as.character(dep_i), as.character(rho_lim)))

inCT %>% 
  arrange(date_time_ID) %>% 
  ggplot()+
  geom_hline(yintercept = 0)+
  geom_point(data = cruise_dates, aes(date_time_ID, 0), shape=21)+
  geom_line(aes(date_time_ID, nCT_i_cum,
                group=group), col="grey")+
  geom_line(data = inCT_fixed_dep %>% filter(dep_i==12, sign=="total"),
            aes(date_time_ID, nCT_i_cum, col="12m - total"))+
  geom_line(data = inCT_MLD %>% filter(rho_lim == 0.1),
            aes(date_time_ID, nCT_i_cum, col="MLD - 0.1"))+
  scale_color_brewer(palette = "Set1", name="")+
  labs(y="inCT [mol/m2]", x="")

rm(inCT, inCT_MLD)

```

# NCP determination

In order to derive an estimate of the net community production NCP (which is equivalent to the formed organic matter that can be exported from the investigated surface layer), two steps are required:

- decision about the most appropiate iC~T~ trajectory
- correction of quantifyable CO~2~ fluxes in and out of the investigated water volume during the period of interest, this includes:
  - Air-sea CO~2~ fluxes
  - CO~2~ fluxes due to vertical mixing
  - CO~2~ fluxes due to lateral transport of water masses (not corrected here)

## Best iC~T~ estimate

To determine the optimum depth for the nC~T~ integration we investigated the vertical distribution of cumulative temperature and nC~T~ changes on the peak of the productivity signal on June 23:


```{r i_nCT_depth_distribution_180723}

tm_profiles_ID_long_180723 <- tm_profiles_ID_long %>% 
  filter(ID == 180723,
         var == "nCT")

p_tm_profiles_ID_long <- tm_profiles_ID_long_180723 %>% 
  arrange(dep) %>% 
  ggplot(aes(value_cum, dep))+
  geom_vline(xintercept = 0)+
  geom_hline(yintercept = 12, col="red")+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  labs(x="Cumulative change of nCT on July 23 (180723)")+
  theme(legend.position = "left")

tm_profiles_ID_long_180723_dep <- tm_profiles_ID_long_180723 %>% 
  select(dep, value_cum) %>% 
  filter(value_cum < 0) %>%
  arrange(dep) %>% 
  mutate(value_cum_i = sum(value_cum),
         value_cum_dep = cumsum(value_cum),
         value_cum_i_rel = value_cum_dep/value_cum_i*100)

p_tm_profiles_ID_long_rel <- tm_profiles_ID_long_180723_dep %>% 
  ggplot(aes(value_cum_i_rel, dep))+
  geom_hline(yintercept = 12, col="red")+
  geom_vline(xintercept = 90)+
  geom_point()+
  geom_line()+
  scale_y_reverse(limits = c(25,0))+
  scale_x_continuous(breaks = seq(0,100,10))+
  labs(y = "Depth (m)", x = "Relative contribution on July 23")+
  theme_bw()

p_tm_profiles_ID_long + p_tm_profiles_ID_long_rel

rm(tm_profiles_ID_long_180723,
   tm_profiles_ID_long_180723_dep,
   p_tm_profiles_ID_long,
   p_tm_profiles_ID_long_rel)

```


```{r i_tem_distribution_180723}

tm_profiles_ID_long_180723 <- tm_profiles_ID_long %>% 
  filter(ID == 180723,
         var == "tem")

p_tm_profiles_ID_long <- tm_profiles_ID_long_180723 %>% 
  arrange(dep) %>% 
  ggplot(aes(value_cum, dep))+
  geom_vline(xintercept = 0)+
  geom_hline(yintercept = 12, col="red")+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  labs(x="Cumulative change of Temp on July 23")+
  theme(legend.position = "left")

tm_profiles_ID_long_180723_dep <- tm_profiles_ID_long_180723 %>% 
  select(dep, value_cum) %>% 
  filter(value_cum > 0) %>%
  arrange(dep) %>% 
  mutate(value_cum_i = sum(value_cum),
         value_cum_dep = cumsum(value_cum),
         value_cum_i_rel = value_cum_dep/value_cum_i*100)

p_tm_profiles_ID_long_rel <- tm_profiles_ID_long_180723_dep %>% 
  ggplot(aes(value_cum_i_rel, dep))+
  geom_hline(yintercept = 12, col="red")+
  geom_vline(xintercept = 90)+
  geom_point()+
  geom_line()+
  scale_y_reverse(limits = c(25,0))+
  scale_x_continuous(breaks = seq(0,100,10))+
  labs(y = "Depth (m)", x = "Relative contribution on July 23")+
  theme_bw()

p_tm_profiles_ID_long + p_tm_profiles_ID_long_rel

rm(tm_profiles_ID_long_180723,
   tm_profiles_ID_long_180723_dep,
   p_tm_profiles_ID_long,
   p_tm_profiles_ID_long_rel)

```


The cummulative iC~T~ trajectory determined by integration of C~T~ to a fixed water depth of 12 m was used for NCP calculation for the following reasons:

- During the first productivity pulse that lasted until July 23:
  - no negative nC~T~ changes were detected below that depth
  - cumulative nC~T~ switch sign at that depth
  - 95% of the cumulative warming signal appears across that depth
- MLD were too shallow to cover all observed negative C~T~ changes
  
## Air-Sea CO~2~ flux

### Surface water data

The cruise mean pCO~2~ recorded in profiling-mode (stations only) and depths < 6m was used for gas exchange calcualtions.

```{r surface_water_time_series}

tm_profiles_surface_long <- tm_profiles %>% 
  filter(dep < surface_dep) %>% 
  select(date_time = date_time_ID, ID, tem, pCO2_water = pCO2, nCT) %>% 
  pivot_longer(tem:nCT, values_to = "value", names_to = "var")

tm_profiles_surface_long_ID <- tm_profiles_surface_long %>% 
  group_by(ID, date_time, var) %>% 
  summarise_all(list(~mean(.), ~sd(.), ~min(.), ~max(.))) %>% 
  ungroup()

# tb_surface_ID <- tb_surface %>% 
#   group_by(ID, var) %>% 
#   summarise(date_time_ID = mean(date_time_ID, na.rm = TRUE),
#             value_mean = mean(value, na.rm = TRUE),
#             value_sd = sd(value, na.rm = TRUE)) %>% 
#   ungroup()

p_pCO2_surf <- tm_profiles_surface_long_ID %>%
  filter(var == "pCO2_water") %>% 
  ggplot(aes(x=date_time))+
  geom_ribbon(aes(ymin=mean-sd, ymax=mean+sd), alpha=0.2)+
  geom_path(aes(y=mean))+
  geom_point(aes(y=mean))+
  scale_fill_discrete(guide=FALSE)+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = expression(atop(pCO[2], (mu*atm))),
       title = "Surface water observations")+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

p_tem_surf <- tm_profiles_surface_long_ID %>%
  filter(var == "tem") %>% 
  ggplot(aes(x=date_time))+
  geom_ribbon(aes(ymin=mean-sd, ymax=mean+sd), alpha=0.2)+
  geom_path(aes(y=mean))+
  geom_point(aes(y=mean))+
  scale_fill_discrete(guide=FALSE)+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = "temp \n (\u00B0C)")+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

p_nCT_surf <-
tm_profiles_surface_long_ID %>%
  filter(var == "nCT") %>% 
  ggplot()+
  geom_point(data = tb_surface_station_mean %>% 
               filter(var == "nCT"),
             aes(x = date_time_ID,
                 y = value_mean,
                 color = "discrete")) +
  geom_linerange(data = tb_surface_station_mean %>%
                   filter(var == "nCT"),
             aes(x = date_time_ID,
                 ymin = value_mean - value_sd,
                 ymax = value_mean + value_sd,
                 color = "discrete"))+
  geom_ribbon(aes(x=date_time, ymin=mean-sd, ymax=mean+sd), alpha=0.2)+
  geom_path(aes(x=date_time, y=mean))+
  geom_point(aes(x=date_time, y=mean))+
  scale_color_manual(values = "red")+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = expression(atop(nC[T],
                           (mu*mol~kg^{-1}))))+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.position = c(0.35,0.75),
        legend.title = element_blank(),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(colour = "black", fill = "white"))


p_pCO2_surf + p_tem_surf + p_nCT_surf + 
  plot_layout(ncol = 1)

#rm(tm_profiles_surface)

start <- min(tm_profiles_surface_long_ID$date_time)
end   <- max(tm_profiles_surface_long_ID$date_time)

```

### Wind and atm. pCO~2~

Metrological data were recorded on the flux tower located on Ostergarnsholm island.

```{r read_wind_Ostergarn_tower}

og <- read_csv(here::here("data/_summarized_data_files",
                          "og.csv"))

og <- og %>%
  filter(date_time > start,
         date_time < end)

rm(end, start)

```

Data sets for atmospheric and seawater observations were merged and interpolated to a common time stamp.

```{r join_interpolate_wind_water_data}

tm_profiles_surface_ID <- tm_profiles_surface_long_ID %>% 
  filter(var %in% c("pCO2_water", "tem")) %>% 
  select(date_time:mean) %>% 
  pivot_wider(names_from = "var", values_from = "mean")

tm_og <- full_join(og, tm_profiles_surface_ID) %>% 
  arrange(date_time)

tm_og <- tm_og %>% 
  mutate(pCO2_water = approxfun(date_time, pCO2_water)(date_time),
         tem = approxfun(date_time, tem)(date_time),
         wind = approxfun(date_time, wind)(date_time)) %>% 
  filter(!is.na(pCO2_atm))

rm(tm_profiles_surface_ID, og)

```

```{r calculate_running_mean, eval=FALSE}

rolling_mean   <- rollify(~mean(.x, na.rm = TRUE), window = 48)

tm_og <- tm_og %>% 
  mutate(wind_daily = rolling_mean(wind),
         pCO2_atm_daily = rolling_mean(pCO2_atm))

```



```{r Ostergarn_tower_time_series}

p_pCO2_atm <- tm_og %>%
  ggplot(aes(x=date_time))+
  geom_path(aes(y=pCO2_atm))+
  # geom_path(aes(y=pCO2_atm_daily, col="24h"))+
  # scale_color_manual(values = c("red", "black"),
  #                    guide = FALSE)+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = expression(atop(pCO["2,atm"], (mu*atm))),
       title = "Atmospheric observations")+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

p_wind <- tm_og %>%
  ggplot(aes(x=date_time))+
  geom_path(aes(y=wind))+
  # geom_path(aes(y=wind_daily, col="24h"))+
  # scale_color_brewer(palette = "Set1")+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = expression(atop(windspeed, (m~s^{-1}))))+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank())

p_pCO2_atm + p_wind +
  plot_layout(ncol = 1) + 
  plot_layout(guides = 'collect')

```

### Air-sea fluxes

F = k * dCO2

with

dCO2 = K0 * dpCO2 and 

k = coeff * U^2 * (660/Sc)^0.5


Unitm used here are:

- dpCO~2~: µatm
- K0: mol atm-1 kg-1
- dCO~2~: µmol kg-1

- wind speed U: m s-1

- coeff for k calculation (eg 0.251 in W14): cm hr-1 (m s-1)-2 
- gas transfer velocities k: cm hr-1 (= 60 x 60 x 100 m s-1)

- air sea CO~2~ flux F: mol m–2 d–1

- conversion between the unit of k * dCO2 and F requires a factor of 10-5 * 24

```{r gas_exchange_calc}

Sc_W14 <- function(tem) {
  2116.8 - 136.25 * tem + 4.7353 * tem^2 - 0.092307 * tem^3 + 0.0007555 * tem^4
}

# Sc_W14(20)

tm_og <- tm_og %>% 
  mutate(dpCO2 = pCO2_water - pCO2_atm,
         dCO2  = dpCO2 * K0(S=6.92, T=tem),
         # W92 = gas_transfer(t = tem, u10 = wind, species = "CO2",
         #                      method = "Wanninkhof1")* 60^2 * 100,
         #k_SM18 = 0.24 * wind^2 * ((1943-119.6*tem + 3.488*tem^2 - 0.0417*tem^3) / 660)^(-0.5),
         k = 0.251 * wind^2 * (Sc_W14(tem)/660)^(-0.5))
  # pivot_longer(9:10, names_to = "k_para", values_to = "k_value")

# calculate flux F [mol m–2 d–1]

tm_og <- tm_og %>% 
  mutate(flux = k*dCO2*1e-5*24)
#         flux_daily = rolling_mean(flux))

rm(Sc_W14)

```

```{r flux_air_sea_daily_CO2, fig.asp=0.3}

p_flux_daily <- tm_og %>%
  ggplot(aes(x=date_time))+
  geom_path(aes(y=flux))+
  # geom_path(aes(y=flux_daily, col="24h"))+
  # scale_color_brewer(palette = "Set1")+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = expression(atop(flux[daily], (mol~m^{-2}~d^{-1}))),
       title = "Air-sea fluxes")+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank())

```

```{r flux_air_sea_cumulative_CO2, fig.asp=0.3}
# scale flux to time interval

tm_og <- tm_og %>% 
  mutate(scale = 24*2) %>% 
  mutate(flux_scale = flux / scale) %>% 
  #group_by(freq, k_para) %>% 
  arrange(date_time) %>% 
  mutate(flux_cum = cumsum(flux_scale)) %>% 
  ungroup()

p_flux_cum <- tm_og %>%
  ggplot(aes(x=date_time))+
  geom_path(aes(y=flux_cum))+
  scale_fill_discrete(guide=FALSE)+
  scale_x_datetime(date_breaks = "week",
                   sec.axis = dup_axis())+
  labs(y = expression(atop(flux[cum],
                           (mol~m^{-2}))))+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

p_flux_daily + p_flux_cum +
  plot_layout(ncol = 1)

```

## iC~T~ correction

The cumulative integrated nCT (inCT) time series obtained through integration across the upper 12m of the water column was used for further calculations of NCP.

Correction of inCT for air-sea CO2 fluxes will be based on estimates derived from observation with 30min measurement interval and calculation according to Wanninkhof (2014).

To derive an integrated NCP estimated, the observed change in inCT must be corrected for the air-sea flux of CO2. inCT was determined for the upper 12m of the water column. The MLD was always shallower 12m, except for the last cruise day. Therefore:

- Cumulative air-sea fluxes can be added completely to inCT before Aug 7.
- Between Aug 7 and the last cruise on Aug 15 it was assumed, that the CO2 flux was homogenously mixed down to the deepend thermocline at 17m. The flux correction applied to the upper 12m can therefore be scaled with a factor 12/17.
    
During the last cruise, deeper mixing up to 17m water depth was observed, resulting in increased inCT at 0-12 m and a decrease of inCT in 12-17m. The loss of nCT in 12-17m can be assumed to be entirely cause by mixing with low-nCT surface water. However, some of the observed nCT loss is balanced through nCT input attributable to the air-sea flux. Therefore, the observed loss, corrected for 5/17 of the air-sea-flux, was added to the integrated nCT changes in 0-12m.

```{r inCT_correction_criteria}

i_dep_lim <- 12
i_dep_mix_lim <- 17

```



```{r inCT_correction, fig.asp=0.7}

# extract CT data for fixed depth approach, depth limit 10m
inCT_10 <- inCT_fixed_dep %>% 
  filter(dep_i == i_dep_lim, sign=="total") %>% 
  select(-c(sign, dep_i))

rm(inCT_fixed_dep)

inCT_10 <- inCT_10 %>%
  select(ID, date_time = date_time_ID, date_time_ID_ref, nCT_i_diff, nCT_i_cum)

# date of the second last cruise
date_180806 <- unique(inCT_10$date_time)[7]

# calculate cumulative air-sea fluxes affecting 0-10m
tm_og_flux <- tm_og %>% 
  mutate(flux_scale = if_else(date_time > date_180806,
                              i_dep_lim/i_dep_mix_lim * flux_scale,
                              flux_scale)) %>%
  arrange(date_time) %>% 
  mutate(flux_cum = cumsum(flux_scale)) %>% 
  select(date_time, flux_cum)

# calculate cumulative air-sea fluxes affecting 10-17m
tm_og_flux_dep <- tm_og %>% 
  filter(date_time > date_180806) %>% 
  mutate(flux_scale = (i_dep_mix_lim - i_dep_lim)/i_dep_mix_lim * flux_scale) %>%
  arrange(date_time) %>% 
  mutate(flux_cum = cumsum(flux_scale)) %>% 
  select(date_time, flux_cum)

inCT_10_flux <- full_join(inCT_10, tm_og_flux) %>% 
  arrange(date_time)

rm(tm_og_flux, inCT_10, tm_og)


# linear interpolation of cumulative changes to frequency of the flux estimates estimates 
inCT_10_flux <- inCT_10_flux %>% 
  mutate(nCT_i_cum = approxfun(date_time, nCT_i_cum)(date_time),
         flux_cum = approxfun(date_time, flux_cum)(date_time)) %>% 
    fill(flux_cum) %>% 
  mutate(nCT_i_flux_cum = nCT_i_cum + flux_cum)


# calculate cumulative fluxes inbetween cruises 
inCT_10_flux_diff <- inCT_10_flux %>% 
  filter(!is.na(date_time_ID_ref)) %>% 
  mutate(flux_diff = flux_cum - lag(flux_cum, default = 0)) %>% 
  select(date_time_ID_ref, observed=nCT_i_diff, flux=flux_diff) %>% 
  pivot_longer(cols = 2:3, names_to = "var", values_to = "value_diff")


# calculate mixing with deep waters, corrected for air sea fluxes
inCT_10_mix <- tm_profiles_ID_long %>% 
  filter(ID == "180815",
         var == "nCT",
         dep < i_dep_mix_lim,
         dep > i_dep_lim) %>% 
  group_by(ID, date_time_ID,
           date_time_ID_ref) %>% 
  summarise(value_diff = sum(value_diff)/1000 + min(tm_og_flux_dep$flux_cum)) %>% 
  ungroup()

rm(tm_og_flux_dep)

inCT_10_mix_diff <- inCT_10_mix %>% 
  select(date_time_ID_ref, value_diff) %>% 
  mutate(var="mixing")

inCT_10_flux_mix_diff <- 
  full_join(inCT_10_flux_diff, inCT_10_mix_diff)

inCT_10_flux_mix <-
  full_join(inCT_10_flux,
            inCT_10_mix %>% rename(mix_cum = value_diff))

rm(inCT_10_mix, inCT_10_mix_diff, inCT_10_flux, inCT_10_flux_diff, date_180806)

inCT_10_flux_mix <- inCT_10_flux_mix %>% 
  arrange(date_time) %>% 
  fill(ID) %>% 
  mutate(mix_cum = if_else(ID %in% c("180806", 180815), mix_cum, 0),
         mix_cum = na.approx(mix_cum),
         nCT_i_flux_mix_cum = nCT_i_flux_cum + mix_cum)

  
p_inCT <- inCT_10_flux_mix %>% 
  arrange(date_time) %>% 
  ggplot()+
  geom_col(data = inCT_10_flux_mix_diff,
           aes(date_time_ID_ref, value_diff, fill=var),
           position = position_dodge2(preserve = "single"),
           alpha=0.5)+
  geom_hline(yintercept = 0)+
  geom_point(data = cruise_dates, aes(date_time_ID, 0), shape=21)+
  geom_line(aes(date_time, nCT_i_cum, col="observed"))+
  geom_line(aes(date_time, nCT_i_flux_mix_cum, col="mixing + flux corrected"))+
  geom_line(aes(date_time, nCT_i_flux_cum, col="flux corrected"))+
  scale_x_datetime(date_breaks = "week",
                   date_labels = "%b %d",
                   sec.axis = dup_axis())+
  scale_fill_brewer(palette = "Set1", name="incremental changes")+
  scale_color_brewer(palette = "Set1", name="cumulative changes")+
  labs(y=expression(atop(integrated~nC[T], (mol~m^{-2}))),
       title = "Water column inventory changes")+
  guides(guide_colourbar(order = 1))+
  theme(axis.title.x = element_blank(),
        axis.text.x.top = element_blank(),
        legend.position = "bottom",
        legend.direction = "vertical")

p_inCT

inCT_10_flux_mix %>%
  write_csv(here::here("Data/_merged_data_files", "tm_NCP_cum.csv"))

inCT_10_flux_mix_diff %>%
  write_csv(here::here("Data/_merged_data_files", "tm_NCP_inc.csv"))

```



```{r plot_inCT_time_series, include=FALSE}

p_pCO2_surf + p_tem_surf + p_nCT_surf + 
  p_pCO2_atm + p_wind + p_flux_daily + p_flux_cum + 
  p_inCT +
  plot_layout(ncol = 1,
              heights = c(rep(1,7), 4))

ggsave(here::here("output/Plots/Figures_publication/article", "atm_water_timeseries.pdf"),
       width = 140, height = 270, dpi = 300, units = "mm")

```

# Open tasks / questions

- clean and harmonize chunk labeling (label: plot, 1 plot per chunk, etc)
- included removed stations in coverage plot
- Significance of changes in AT for calculated nCT changes
  - Calculate AT-S ratios, reconstruct AT profiles, calculate true nCT profiles, normalize nCT profiles to mean AT
- demostrate strong permanent thermocline at around 25 m