BloomSail_surface.Rmd

---
title: "MLD"
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(seacarb)
library(oce)
library(marelac)
library(patchwork)
library(metR)
```


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

# Approach

In order to test how (and how well) the depth-integrated C~T~ estimates can be reproduced if only surface CO2 data were available, the BloomSail observations were restricted to those made in surface water and two reconstruction approaches were tested:

1. MLD: Integration of surface observation across the MLD, assuming homogenious vertical patterns
2. Temperature: Vertical reconstruction of incremental C~T~ changes based on profiles of incremental changes in temperature

# Sensor data

1m gridded, downcast profiles, from which CO~2~ data other than 3.5 m water depth were removed, were used.

```{r read_Tina-V_Sensor_data}

ts_profiles_ID <-
  read_csv(here::here("Data/_merged_data_files", "ts_profiles_ID_long_cum_MLD.csv"))

ts_profiles_ID <- ts_profiles_ID %>% 
  select(ID, date_time_ID, dep, CT = value, CT_diff = value_diff, CT_cum = value_cum, sign, rho_lim, MLD)

ts_profiles_ID <- ts_profiles_ID %>% 
  mutate(CT = if_else(dep == 3.5, CT, NaN),
         rho_lim = as.factor(rho_lim))

```


# MLD approach

MLD calculation was previously described in the [CT dynamics](https://jens-daniel-mueller.github.io/BloomSail/CT_dynamics.html) chapter.

## Timeseries

For comparison to MLD, the effective penetration depth of NCP, z~eff~, was calculated as the ratio of the inremental, depth-integrated change of C~T~, divided by the change in surface C~T~, for all cases where the change in surface C~T~ was negative.

```{r NCP_penetration_depth}

ts_profiles_ID_surface <- ts_profiles_ID %>% 
  filter(sign=="neg",
         rho_lim == "0.5",
         dep == 3.5) %>% 
  group_by(ID, date_time_ID) %>% 
  summarise(CT_diff_surf = mean(CT_diff, na.rm = TRUE)) %>% 
  ungroup()

ts_profiles_ID_i <- ts_profiles_ID %>% 
  filter(sign=="neg",
         rho_lim == "0.5") %>% 
  group_by(ID, date_time_ID) %>% 
  summarise(CT_diff_i = sum(CT_diff, na.rm = TRUE))

zeff <- inner_join(ts_profiles_ID_surface, ts_profiles_ID_i)
rm(ts_profiles_ID_surface, ts_profiles_ID_i)

zeff <- zeff %>% 
  mutate(zeff = CT_diff_i / CT_diff_surf)

```



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

ts_profiles_ID %>% 
  ggplot()+
  geom_hline(yintercept = 0)+
  geom_line(data = zeff, aes(date_time_ID, zeff, linetype="zeff"))+
  geom_line(aes(date_time_ID, MLD, col=rho_lim, linetype="MLD"))+
  scale_y_reverse()+
  scale_color_viridis_d(name="Rho limit")+
  scale_linetype(name="Estimate")+
  labs(y="Depth (m)")+
  theme(axis.title.x = element_blank())

rm(zeff)

```

## iC~T~ calculation

Integrated C~T~ depletion was calculated as the product of observed incremental C~T~ changes in surface waters and the respective mixed layer depth. 

```{r NCP_calculation}

ts_profiles_ID_surface <- ts_profiles_ID %>% 
  filter(dep==3.5) %>% 
  group_by(rho_lim) %>% 
  arrange(date_time_ID) %>% 
  mutate(iCT_diff = CT_diff * MLD / 1000,
         iCT_cum = cumsum(replace_na(iCT_diff, 0))) %>% 
  ungroup()
  
```


## Incremental and cumulative timeseries

Total incremental and cumulative C~T~ changes inbetween cruise dates were calculated.


```{r plot_integrated_NCP_timeseries, fig.asp=1.5}

p_iCT <- ts_profiles_ID_surface %>% 
  ggplot(aes(date_time_ID, iCT_diff, fill= rho_lim))+
  geom_hline(yintercept = 0)+
  geom_col(col="black", position = "dodge")+
  scale_y_continuous(breaks = seq(-100, 100, 0.2))+
  scale_fill_viridis_d()+
  labs(y="integrated CT changes [mol/m2]")+
  theme(axis.title.x = element_blank())

p_iCT_cum <- ts_profiles_ID_surface %>% 
  ggplot(aes(date_time_ID,  iCT_cum, 
             col=rho_lim))+
  geom_line()+
  geom_hline(yintercept = 0)+
  scale_color_viridis_d()+
  scale_y_continuous(breaks = seq(-100, 100, 0.2))+
  theme(strip.background = element_blank(),
        strip.text = element_blank())+
  labs(y="integrated, cumulative CT changes [mol/m2]", x="date")


(p_iCT / p_iCT_cum)+
  plot_layout(guides = 'collect')

rm(p_iCT, p_iCT_cum)

rm(ts_profiles_ID, ts_profiles_ID_surface)

```


# Temperature approach


```{r read_ts_data}

ts_profiles_ID <-
  read_csv(here::here("Data/_merged_data_files", "ts_profiles_ID.csv"))

ts_profiles_ID <- ts_profiles_ID %>% 
  mutate(CT = if_else(dep == 3.5, CT, NaN),
         ID = as.factor(ID))

```


## C~T~ vs temperature

As primary production (negative changes in C~T~) and increase in seawater temperature have a common driver (light), the relation between both changes was investigated.

```{r CT_and_temperature_changes_in_surface}

ts_profiles_ID_diff <- ts_profiles_ID %>% 
  drop_na() %>% 
  arrange(date_time_ID) %>%
  mutate(CT_diff = CT - lag(CT, default = first(CT)),
         tem_diff = tem - lag(tem, default = first(tem)),
         factor = CT_diff / tem_diff,
         factor = if_else(is.na(factor), 0, factor))

ts_profiles_ID_diff %>% 
  ggplot(aes(tem_diff, CT_diff))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  geom_path()+
  geom_point()

```


## Reconstruction of C~T~ dynamics

The ratio of the incremental change of C~T~ with temperature at the seasurface was applied to calculate the C~T~ in other water depth based on the known change in temperature.

```{r reconctruction_CT_incremental_changes_profiles}

ts_profiles_ID_diff <- ts_profiles_ID_diff %>% 
  select(ID, factor)

ts_profiles_ID <- full_join(ts_profiles_ID, ts_profiles_ID_diff)
rm(ts_profiles_ID_diff)

ts_profiles_ID <- ts_profiles_ID %>%
  arrange(date_time_ID) %>%
  group_by(dep) %>% 
  mutate(tem_diff = tem - lag(tem, default = first(tem))) %>% 
  ungroup()

ts_profiles_ID <- ts_profiles_ID %>% 
  mutate(CT_diff = tem_diff * factor) %>% 
  select(-factor)

```

The reconstructed incremental changes are added up to derive cummulative C~T~ changes throughout the water column.

```{r reconctruction_CT_cumulative_changes_profiles}

ts_profiles_ID_long <- ts_profiles_ID %>% 
  select(ID, date_time_ID, dep, tem = tem_diff, CT = CT_diff) %>% 
  pivot_longer(4:5, names_to = "var", values_to = "value_diff") %>% 
  group_by(dep, var) %>% 
  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_daily = value_diff / date_time_ID_diff,
         value_cum = cumsum(value_diff)) %>% 
  ungroup()

```

## Profiles of incremental changes

Changes of seawater parameters at each depth were reconstructed from one cruise day to the next and divided by the number of days inbetween.

```{r incremental_changes_profiles}

ts_profiles_ID_long %>% 
  arrange(dep) %>% 
  ggplot(aes(value_diff, 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 parameters were calculated at each depth relative to the first cruise day on July 5.

```{r cumulative_changes_profiles}

ts_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()+
  labs(x="Cumulative change of value")+
  facet_wrap(~var, scales = "free_x")

```

## Hovmoeller plots

Hoevmoeller plots were generated for the reconstructed daily and cumulative changes in C~T~. Absolute values are not reproducible with this approach. Furthermore, it meets our expectations 

### Daily changes

```{r plot_hovmoeller_daily, fig.cap="Hovmoeller plots of daily reconstructed changes in C~T~ and temperature. Note: Daily changes are currently plotted against the day when they were observed compared to the previous transect, although plotting against the mean date would be more plausible.", fig.asp=0.5, eval=FALSE}

bin_CT <- 2.5

ts_profiles_ID_long %>%
  filter(var == "CT") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID_ref, y=dep, z=value_diff_daily),
                    breaks = MakeBreaks(bin_CT),
                    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_CT),
                       guide = "colorstrip",
                       name="CT (µmol/kg)")+
  scale_y_reverse()+
  theme_bw()+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

rm(bin_CT)

```

### Cumulative changes

```{r plot_hovmoeller_cum, fig.cap="Hovmoeller plots of daily reconstructed changes in C~T~ and temperature.", fig.asp=0.5, eval=FALSE}

bin_CT <- 20

ts_profiles_ID_long %>%
  filter(var == "CT") %>% 
  ggplot()+
  geom_contour_fill(aes(x=date_time_ID, y=dep, z=value_cum),
                    breaks = MakeBreaks(bin_CT),
                    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_CT),
                       guide = "colorstrip",
                       name="CT (µmol/kg)")+
  scale_y_reverse()+
  theme_bw()+
  labs(y="Depth (m)")+
  coord_cartesian(expand = 0)+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())


rm(bin_CT)

```


## Heat penetration depth

As an alternative approach to the integration over the MLD or the reconstruction of C~T~ profiles, we can estimate the mean penetration depth of the warming signal, which was defined the surface change in temperature devided by the integrated change in seawater temperature across depth.

```{r heat_penetration_depth}

tem_diff_surface <- ts_profiles_ID %>% 
  filter(dep == 3.5) %>% 
  select(ID, tem_diff_surface=tem_diff)

ts_profiles_ID <- full_join(ts_profiles_ID, tem_diff_surface) 
rm(tem_diff_surface)

tem_depth <- ts_profiles_ID %>%
  filter(dep < 18) %>% 
  group_by(ID, date_time_ID) %>% 
  summarise(tem_diff_int = sum(tem_diff),
            tem_diff_surface = mean(tem_diff_surface),
            tem_depth = tem_diff_int/tem_diff_surface) %>% 
  ungroup()

tem_depth %>% 
  ggplot(aes(date_time_ID, tem_depth))+
  geom_hline(yintercept = 0)+
  geom_line()+
  geom_point()+
  scale_y_reverse()

```




## iC~T~ time series

Total incremental and cumulative C~T~ changes inbetween cruise dates were calculated for the upper 10 m of the water body.

```{r integrated_NCP}

iCT_10 <- ts_profiles_ID_long %>%
  filter(dep < 10, 
         var == "CT") %>% 
  select(ID, date_time_ID, date_time_ID_ref, CT_diff=value_diff, CT_cum=value_cum) %>% 
  group_by(ID, date_time_ID, date_time_ID_ref) %>% 
  summarise(CT_i_diff = sum(CT_diff)/1000,
            CT_i_cum = sum(CT_cum)/1000) %>% 
  ungroup()

iCT_10 %>% 
  ggplot()+
  #geom_point(data = cruise_dates, aes(date_time_ID, 0), shape=21)+
  geom_col(aes(date_time_ID_ref, CT_i_diff),
           position = "dodge", alpha=0.3)+
  geom_line(aes(date_time_ID, CT_i_cum))+
  scale_color_viridis_d(name="Depth limit (m)")+
  scale_fill_viridis_d(name="Depth limit (m)")+
  labs(y="iCT [mol/m2]", x="")+
  theme_bw()


```