Introduction

This document outlines the methodology behind the manuscript titled Glacial meltwater increases coastal CO₂ uptake and sensitivity to biogeochemical change. Here, we create seawater–freshwater mixing models to explore carbonate system dynamics in Greenlandic fjords. Additionally, we perform a sensitivity analysis to assess how individual biogeochemical processes influence the marine CO₂ system under different freshening scenarios.

Constructing Mixing Models

To start, we build a seawater–freshwater mixing model using endmember chemical concentrations. These endmembers were measured in August 2023 from seawater in Young Sound, NE Greenland, and glacial meltwater from the Tyroler River. Because salinity, dissolved inorganic carbon (DIC), and total alkalinity (TA) mix conservatively, we applied a linear dilution between the seawater and freshwater endmembers. From the resulting linear gradients in DIC and TA, we calculated the remaining carbonate system parameters—including the concentrations of individual inorganic carbon species—which cannot be directly measured.

Set up data frame using Seawater and Glacial Melt endmembers

# Step 2: Create the dataframe of endmembers
df <- data.frame(
  Parcel = c("River", "SW", "MQ"),
  TA = c(194.85, 2198.0433, 97.06),     # μmolkg−1
  DIC = c(195.95, 2045.953, 30.84),     # μmolkg−1
  Salinity = c(0, 31.32, 0),            # Practical Salinity Scale (PSS-78)
  Temperature = c(0, 0, 0)              # Celsius
)



print(df)
##   Parcel       TA      DIC Salinity Temperature
## 1  River  194.850  195.950     0.00           0
## 2     SW 2198.043 2045.953    31.32           0
## 3     MQ   97.060   30.840     0.00           0
#The TA/DIC ratio of Seawater and Glacial melt are:
df$TA[1]/df$DIC[1] #river
## [1] 0.9943863
df$TA[2]/df$DIC[2] # SW
## [1] 1.074337

Read in Data of mesaured pCO₂ from a controlled seawater-freshwater mixing experiment

measured_mix <- read.csv("~/Library/CloudStorage/OneDrive-Aarhusuniversitet/Aarhus PhD/Young Sound 2023/Analysis/Data/Table of Mixed pCO2.csv", sep=";")

measured_mix <- measured_mix %>% rename(Salinity = salinity)

Now create a mixing model between seawater and glacial melt

# Function to calculate parameters of the resulting mix of River and SW based on conservative mixing of parameters
calculate_mix_parameters <- function(df, mixing_ratio) {
  mix_TA <- df$TA[1] * mixing_ratio + df$TA[2] * (1 - mixing_ratio)
  mix_DIC <- df$DIC[1] * mixing_ratio + df$DIC[2] * (1 - mixing_ratio)
  mix_Salinity <- df$Salinity[1] * mixing_ratio + df$Salinity[2] * (1 - mixing_ratio)
  mix_Temperature <- df$Temperature[1] * mixing_ratio + df$Temperature[2] * (1 - mixing_ratio)
  
  return(data.frame(
    Parcel = "Mix",
    TA = mix_TA,
    DIC = mix_DIC,
    Salinity = mix_Salinity,
    Temperature = mix_Temperature,
    x = mixing_ratio
  ))
}

# Example mixing ratios
mixing_ratios <- seq(0, 1, by = 0.001)

# Calculate mix parameters for each mixing ratio
result <- lapply(mixing_ratios, function(ratio) calculate_mix_parameters(df, ratio))

# Combine the results into a single dataframe
result_df <- do.call(rbind, result)

# Print the resulting dataframe
tibble(result_df)

Using the output from our mixing model, we can calculate pCO₂ from TA and DIC

result_df <- result_df %>% mutate(
TA = TA/1000000, #converting µmol/kg to mol/kg
DIC = DIC/1000000) #converting µmol/kg to mol/kg

carb.df <- carb(flag = 15, 
                var1 = result_df$TA, 
                var2 = result_df$DIC, 
                k1k2 = "l", 
                S=result_df$Salinity, 
                T = result_df$Temperature, 
                eos = "teos10", 
                b = "l10") #calculate pco2 from TA and DIC

carb.df <- carb.df %>% 
  dplyr::rename(TA = ALK,
                pCO2_calc = pCO2,
                pH_calc = pH,
                salinity = S,
                temperature = T,
                pressure = P)

#The ocean acidification community agrees upon using Leuker et al (2000) constants. However, Sulpis et al (2020) illustrate they do not perform the best at low temperatures. 

#Below we also calculate pCO2 using K1 and K2 constants from Sulpis et al (2020)

carb.df_sulpis <- carb(flag = 15, 
                var1 = result_df$TA, 
                var2 = result_df$DIC, 
                k1k2 = "s20", 
                S=result_df$Salinity, 
                T = result_df$Temperature, 
                eos = "teos10", 
                b = "l10") #calculate pco2 from TA and DIC

carb.df_sulpis <- carb.df_sulpis %>% 
  dplyr::rename(TA = ALK,
                pCO2_calc = pCO2,
                pH_calc = pH,
                salinity = S,
                temperature = T,
                pressure = P)




carb.dfsens <- carb.df # set up base dataframe for sensitivity analysis below

tibble(carb.df)

Mixing Plot

# Define fixed ranges for DIC, TA, and pCO2
DIC_range <- range(carb.df$DIC * 1e6, na.rm = TRUE)
TA_range <- range(carb.df$TA * 1e6, na.rm = TRUE)
pCO2_range <- c(min(carb.df$pCO2_calc, na.rm = TRUE), max(measured_mix$Tyro))
pCO2_range <- c(0,407)

# Rescale pCO2 for the plot to match the TA range
pCO2_rescaled <- (carb.df$pCO2_calc - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)
pCO2_rescaled_sulpis <- (carb.df_sulpis$pCO2_calc - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Rescale Tyro points using the same fixed transformation
Tyro_rescaled <- (measured_mix$Tyro - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Combine data into a single dataframe
plot_data <- carb.df %>%
  mutate(
    pCO2_rescaled = pCO2_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC and TA
  geom_line(aes(y = DIC * 1e6, color = "DIC"), size = 3) +
  geom_line(aes(y = TA * 1e6, color = "TA"), size = 3) +

  # Secondary axis: pCO2
  
  geom_line(aes(y = pCO2_rescaled_sulpis, color = "pCO2 (S20)"), linetype = "dashed", size = 2) +
  geom_line(aes(y = pCO2_rescaled, color = "pCO2"), linetype = "dashed", size = 2) +
  geom_point(data = measured_mix, aes(y = Tyro_rescaled, x = Salinity, color = "pCO2"), size = 5) + # Ensure points have the same color mapping

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("DIC & TA (µmol kg"^-1~")"), # Label for primary y-axis
    sec.axis = sec_axis(
      ~ (. - min(TA_range)) / diff(TA_range) * diff(pCO2_range) + min(pCO2_range),
      name = expression(pCO[2] ~ (µatm))
    )
  ) +

  # Color scale and labels
  scale_color_manual(
    values = c("pCO2" = "blue", "pCO2 (S20)" = "grey60", "DIC" = "hotpink1", "TA" = "forestgreen"),
    labels = c("DIC", expression(pCO[2]), expression(pCO[2]~"(s20)"),  "TA"),
    name = "Legend"
  ) +

  # Add labels and theme
  labs(x = "Salinity") +
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "blue"),
    axis.line.y.right = element_line(color = "blue"),
    axis.text.y.right = element_text(color = "blue"),
    legend.position = "top",
    legend.title = element_blank())

Here we find that the calculated pCO₂ is much lower than that of the directly measured pCO₂. Potential sources of this mismatch could be: atmospheric CO₂ intrusion between sample collection and mixing, warming of the water mixture in the lab (not reflected in the calculations), or internal uncertainty in the carbonate system at such cold temperatures and low salinities. Indeed many have reported discrepancies between measured and calculated pCO₂. For example, Sulpis et al. (2020) found that below 8 °C, calculated and measured pCO₂ often diverge. While Sulpis et al. developed their own K₁ and K₂ dissociation constants to improve low-temperature performance, their constants apply to salinities between ~31–38 psu—completely outside the range of our mixtures (0–31 psu).

Indeed, Raimondi et al. (2019) note that “the propagated uncertainty (i.e. combined standard uncertainty) of calculated parameters of the carbonate system is (a) always larger than the analytical precision of the measurements themselves; (b) strongly dependent on the choice of input parameters and uncertainties; (c) less dependent on choice of the specific set of constants.”

Plot of carbonate species for mismatched endmembers

# Rescale pH and other DIC forms for proper secondary y-axis alignment
CO2_rescaled <- scales::rescale(carb.df$CO2 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
CO3_rescaled <- scales::rescale(carb.df$CO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
HCO3_rescaled <- scales::rescale(carb.df$HCO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))

# Combine data into a single dataframe for ggplot
plot_data <- carb.df %>%
  mutate(
    CO2_rescaled = CO2_rescaled,
    CO3_rescaled = CO3_rescaled,
    HCO3_rescaled = HCO3_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC forms
  geom_line(aes(y = CO2_rescaled, color = "CO2"), linetype = "dashed", size = 2) +
  geom_line(aes(y = CO3_rescaled, color = "CO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = HCO3_rescaled, color = "HCO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = pH_calc, color = "pH"), linetype = "solid", size = 2) +

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("pH"),
    sec.axis = sec_axis(
      ~ scales::rescale(., to = range(carb.df$CO2*1e6, na.rm = TRUE)),
      name = expression("DIC form (µmol kg"^-1~")"))) +

  # Color scale and labels
  scale_color_manual(
    values = c("CO2" = "blue", "CO3" = "coral1", "HCO3" = "turquoise2", "pH" = "black"),
    labels = c(expression(CO[2]), expression(CO[3]^"-2"), expression(HCO[3]^"-1"), "pH"),
    name = "Legend"
  ) +

  # Add labels and title
  labs(x = "Salinity", 
       #title = expression(SW + "Tyroler River")
       ) +
  
  # Add theme and styling
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "black"),
    axis.text.y.right = element_text(color = "blue"),
    axis.line.y.right = element_line(color = "blue"),
    legend.position = "top",
    legend.title = element_blank()
  )

Now we can try to adjust DIC slightly to better match measured pCO₂.

If our goal is to understand the source of the pCO₂ nonlinearity observed in the lab mixing experiments, we first need a better match between measured and calculated pCO₂. Only four parameters (pH, pCO₂, TA, and DIC) can be directly measured to describe the marine carbonate system. The concentrations of inorganic carbon species (CO₂, HCO₃⁻, CO₃²⁻) must be calculated. Since pCO₂ is directly related to dissolved CO₂ concentration, any mismatch between measured and calculated pCO₂ suggests similar errors in the inferred concentrations of CO₂, HCO₃⁻, and CO₃²⁻.

To improve this match, we adjust the DIC endmember values by estimating DIC from measured pCO₂ and TA (only two carbonate parameters are required to fully describe the system). Then, using TA and the adjusted DIC values, we rebuild the linear mixing model. This better reflects the evolution of the carbonate system as observed in the controlled mixing experiment.

Set up like before now using calculated DIC for endmembers

#As a reminder, the TA/DIC ratio of Seawater and Glacial melt were:
df$TA[1]/df$DIC[1] #river
## [1] 0.9943863
df$TA[2]/df$DIC[2] # SW
## [1] 1.074337
sw_dic_sat <- carb(flag = 24, 
                   var1 = 402, # measured pCO2 of seawater
                   var2 = df$TA[2]/1000000, 
                   k1k2 = "l", 
                   S = 31.32, 
                   T = 0, 
                   eos = "teos10", 
                   b = "l10")$DIC*1000000 

riv_dic_sat <- carb(flag = 24, 
                    var1 = 397, # measured pCO2 of glacial melt 
                    var2 = df$TA[1]/1000000, 
                    k1k2 = "l", 
                    S = 0, 
                    T = 0, 
                    eos = "teos10", 
                    b = "l10")$DIC*1000000 

#So if we assume that at the measured TA, DIC was adjusted to better reflect measured  pCO2 we would get DICs of:
sw_dic_sat
## [1] 2103.142
riv_dic_sat
## [1] 225.1866
#The change in DIC would then be:
(df$DIC[1] - riv_dic_sat)
## [1] -29.23665
(df$DIC[2] - sw_dic_sat)
## [1] -57.18919
#With a percentage change of:
(riv_dic_sat - df$DIC[1])/df$DIC[1] * 100
## [1] 14.92046
(sw_dic_sat - df$DIC[2])/df$DIC[2] * 100
## [1] 2.795235
#For freshwater umol/kg = umol/L

# Step 2: Create the new dataframe of endmembers using adjusted DIC
df <- data.frame(
  Parcel = c("River", "SW", "MQ"),
  TA = c(194.85, 2198.0433, 97.06),            # μmolkg−1
  DIC = c(riv_dic_sat, sw_dic_sat, 30.84),     # μmolkg−1
  Salinity = c(0, 31.32, 0),                   # Practical Salinity Scale (PSS-78)
  Temperature = c(0, 0, 0)                     # Celsius
)


print(df)
##   Parcel       TA       DIC Salinity Temperature
## 1  River  194.850  225.1866     0.00           0
## 2     SW 2198.043 2103.1422    31.32           0
## 3     MQ   97.060   30.8400     0.00           0
#New TA/DIC ratios
df$TA[1]/df$DIC[1] #river
## [1] 0.8652822
df$TA[2]/df$DIC[2] # SW
## [1] 1.045123

Now calculate mixing ratios for SW + River

# Function to calculate parameters of the resulting mix of River and SW based on conservative mixing of parameters
calculate_mix_parameters <- function(df, mixing_ratio) {
  mix_TA <- df$TA[1] * mixing_ratio + df$TA[2] * (1 - mixing_ratio)
  mix_DIC <- df$DIC[1] * mixing_ratio + df$DIC[2] * (1 - mixing_ratio)
  mix_Salinity <- df$Salinity[1] * mixing_ratio + df$Salinity[2] * (1 - mixing_ratio)
  mix_Temperature <- df$Temperature[1] * mixing_ratio + df$Temperature[2] * (1 - mixing_ratio)
  
  return(data.frame(
    Parcel = "Mix",
    TA = mix_TA,
    DIC = mix_DIC,
    Salinity = mix_Salinity,
    Temperature = mix_Temperature,
    x = mixing_ratio
  ))
}

# Example mixing ratios
mixing_ratios <- seq(0, 1, by = 0.001)

# Calculate mix parameters for each mixing ratio
result <- lapply(mixing_ratios, function(ratio) calculate_mix_parameters(df, ratio))


# Combine the results into a single dataframe
result_df <- do.call(rbind, result)

# Print the resulting dataframe
tibble(result_df)

Calculating pCO₂ again from the mixing model DIC and TA (now adjusted DIC endmembers)

result_df <- result_df %>% mutate(
TA = TA/1000000, #converting µmol/kg to mol/kg
DIC = DIC/1000000) #converting µmol/L to mol/kg  using sigmaT as density conversion factor



carb.df <- carb(flag = 15, var1 = result_df$TA, var2 = result_df$DIC, k1k2 = "l", S=result_df$Salinity, T = result_df$Temperature, eos = "teos10", b = "l10") #calculate pco2 from TA and DIC



carb.df <- carb.df %>% dplyr::rename(TA = ALK,
                                    pCO2_calc = pCO2,
                                    pH_calc = pH,
                                    salinity = S,
                                    temperature = T,
                                    pressure = P) %>%
                       mutate(TADICratio = TA/DIC)

Making plot of mixing model

# Define fixed ranges for DIC, TA, and pCO2
DIC_range <- range(carb.df$DIC * 1e6, na.rm = TRUE)
TA_range <- range(carb.df$TA * 1e6, na.rm = TRUE)
pCO2_range <- range(carb.df$pCO2_calc, na.rm = TRUE)


# Rescale pCO2 for the plot to match the TA range
pCO2_rescaled <- (carb.df$pCO2_calc - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Rescale Tyro points using the same fixed transformation
Tyro_rescaled <- (measured_mix$Tyro - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Combine data into a single dataframe
plot_data <- carb.df %>%
  mutate(
    pCO2_rescaled = pCO2_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC and TA
  geom_line(aes(y = DIC * 1e6, color = "DIC"), size = 3) +
  geom_line(aes(y = TA * 1e6, color = "TA"), size = 3) +

  # Secondary axis: pCO2
  geom_line(aes(y = pCO2_rescaled, color = "pCO2"), linetype = "dashed", size = 2) +
  geom_point(data = measured_mix, aes(y = Tyro_rescaled, x = Salinity, color = "pCO2"), size = 5) + # Ensure points have the same color mapping

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("DIC & TA (µmol kg"^-1~")"), # Label for primary y-axis
    sec.axis = sec_axis(
      ~ (. - min(TA_range)) / diff(TA_range) * diff(pCO2_range) + min(pCO2_range),
      name = expression(pCO[2] ~ (µatm))
    )
  ) +

  # Color scale and labels
  scale_color_manual(
    values = c("pCO2" = "blue", "DIC" = "hotpink1", "TA" = "forestgreen"),
    labels = c("DIC", expression(pCO[2]), "TA"),
    name = "Legend"
  ) +

  # Add labels and theme
  labs(x = "Salinity") +
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "blue"),
    axis.line.y.right = element_line(color = "blue"),
    axis.text.y.right = element_text(color = "blue"),
    legend.position = "top",
    legend.title = element_blank())

When DIC in the glacial meltwater and seawater endmembers is slightly increased (by 29 and 57 µmol kg⁻¹, equivalent to 15% and 2.8%, respectively), the calculated pCO₂ matches the measured values across both endmembers and all intermediate dilution points.

This adjustment helps reveal the sources of nonlinearity we observe: first, DIC dilution by glacial melt reduces pCO₂; then, below ~9 psu, acidification shifts the dominant inorganic carbon species toward CO₂, which increases pCO₂ again.

# Rescale pH and other DIC forms for proper secondary y-axis alignment
CO2_rescaled <- scales::rescale(carb.df$CO2 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
CO3_rescaled <- scales::rescale(carb.df$CO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
HCO3_rescaled <- scales::rescale(carb.df$HCO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))

# Combine data into a single dataframe for ggplot
plot_data <- carb.df %>%
  mutate(
    CO2_rescaled = CO2_rescaled,
    CO3_rescaled = CO3_rescaled,
    HCO3_rescaled = HCO3_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC forms
  geom_line(aes(y = CO2_rescaled, color = "CO2"), linetype = "dashed", size = 2) +
  geom_line(aes(y = CO3_rescaled, color = "CO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = HCO3_rescaled, color = "HCO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = pH_calc, color = "pH"), linetype = "solid", size = 2) +

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("pH"),
    sec.axis = sec_axis(
      ~ scales::rescale(., to = range(carb.df$CO3*1e6, na.rm = TRUE)),
      name = expression("DIC form (µmol kg"^-1~")"))) +

  # Color scale and labels
  scale_color_manual(
    values = c("CO2" = "blue", "CO3" = "coral1", "HCO3" = "turquoise2", "pH" = "black"),
    labels = c(expression(CO[2]), expression(CO[3]^"-2"), expression(HCO[3]^"-1"), "pH"),
    name = "Legend"
  ) +

  # Add labels and title
  labs(x = "Salinity", 
       #title = expression(SW + "Tyroler River")
       ) +
  
  # Add theme and styling
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "black"),
    axis.text.y.right = element_text(color = "coral1"),
    axis.line.y.right = element_line(color = "coral1"),
    legend.position = "top",
    legend.title = element_blank()
  )

How does the glacial endmember affect the shape of pCO₂ nonlinearity?

We saw above that changing the TA:DIC ratio (by altering the endmember composition) impacted the nonlinear reduction of pCO₂ during mixing. To further examine how the glacial endmember affects this nonlinearity, we can rerun mixing models using different glacial endmembers. This allows us to explore how the carbonate system evolves under varying scenarios, representative of freshwater sources found in different catchments across the Arctic.

To do this, we keep the seawater endmember fixed at the values measured in Young Sound.

Set up data frame with a freshwater endmember having a higher TA:DIC ratio than seawater

# Step 2: Create the dataframe of endmembers
df <- data.frame(
  Parcel = c("River", "SW", "MQ"),
  TA = c(194.85, 2198.0433, 97.06),     # μmolkg−1
  DIC = c(178, 2045.953, 30.84),        # μmolkg−1
  Salinity = c(0, 31.32, 0),            # Practical Salinity Scale (PSS-78)
  Temperature = c(0, 0, 0)              # Celsius
)


df
df$TA[1]/df$DIC[1] #river
## [1] 1.094663
df$TA[2]/df$DIC[2] # SW
## [1] 1.074337
194/1.089474
## [1] 178.0676
carb(flag = 15, var1 = df$TA[2]/1000000, var2 = df$DIC[2]/1000000, k1k2 = "l", S = 31.32, T = 0 , eos = "teos10", b = "l10")$pCO2
## [1] 263.8324
carb(flag = 15, var1 = df$TA[1]/1000000, var2 = df$DIC[1]/1000000, k1k2 = "l", S = 0, T = 0 , eos = "teos10", b = "l10")$pCO2
## [1] 6.500036
#matches atm pco2

now construct mixing model

# Function to calculate parameters of the resulting mix of River and SW
calculate_mix_parameters <- function(df, mixing_ratio) {
  mix_TA <- df$TA[1] * mixing_ratio + df$TA[2] * (1 - mixing_ratio)
  mix_DIC <- df$DIC[1] * mixing_ratio + df$DIC[2] * (1 - mixing_ratio)
  mix_Salinity <- df$Salinity[1] * mixing_ratio + df$Salinity[2] * (1 - mixing_ratio)
  mix_Temperature <- df$Temperature[1] * mixing_ratio + df$Temperature[2] * (1 - mixing_ratio)
  
  return(data.frame(
    Parcel = "Mix",
    TA = mix_TA,
    DIC = mix_DIC,
    Salinity = mix_Salinity,
    Temperature = mix_Temperature,
    x = mixing_ratio
  ))
}

# Example mixing ratios
mixing_ratios <- seq(0, 1, by = 0.001)

# Calculate mix parameters for each mixing ratio
result <- lapply(mixing_ratios, function(ratio) calculate_mix_parameters(df, ratio))

# Combine the results into a single dataframe
result_df <- do.call(rbind, result)


# Print the resulting dataframe
tibble(result_df)

Calculating carbonate system

result_df <- result_df %>% mutate(
TA = TA/1000000, #converting µmol/kg to mol/kg
DIC = DIC/1000000) #converting µmol/L to mol/kg  


carb.df <- carb(flag = 15, 
                var1 = result_df$TA, 
                var2 = result_df$DIC, 
                k1k2 = "l", 
                S=result_df$Salinity, 
                T = result_df$Temperature, 
                eos = "teos10", 
                b = "l10") #calculate pco2 from TA and DIC




carb.df <- carb.df %>% dplyr::rename(TA = ALK,
                                    pCO2_calc = pCO2,
                                    pH_calc = pH,
                                    salinity = S,
                                    temperature = T,
                                    pressure = P)

Plot

# Define fixed ranges for DIC, TA, and pCO2
DIC_range <- range(carb.df$DIC * 1e6, na.rm = TRUE)
TA_range <- range(carb.df$TA * 1e6, na.rm = TRUE)
pCO2_range <- range(carb.df$pCO2_calc, na.rm = TRUE)
pCO2_range <- c(0, 407)

# Rescale pCO2 for the plot to match the TA range
pCO2_rescaled <- (carb.df$pCO2_calc - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Rescale Tyro points using the same fixed transformation
Tyro_rescaled <- (measured_mix$Tyro - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Combine data into a single dataframe
plot_data <- carb.df %>%
  mutate(
    pCO2_rescaled = pCO2_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC and TA
  geom_line(aes(y = DIC * 1e6, color = "DIC"), size = 3) +
  geom_line(aes(y = TA * 1e6, color = "TA"), size = 3) +

  # Secondary axis: pCO2
  geom_line(aes(y = pCO2_rescaled, color = "pCO2"), linetype = "dashed", size = 2) +

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("DIC & TA (µmol kg"^-1~")"), # Label for primary y-axis
    sec.axis = sec_axis(
      ~ (. - min(TA_range)) / diff(TA_range) * diff(pCO2_range) + min(pCO2_range),
      name = expression(pCO[2] ~ (µatm))
    )
  ) +

  # Color scale and labels
  scale_color_manual(
    values = c("pCO2" = "blue", "DIC" = "hotpink1", "TA" = "forestgreen"),
    labels = c("DIC", expression(pCO[2]), "TA"),
    name = "Legend"
  ) +

  # Add labels and theme
  labs(x = "Salinity") +
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "blue"),
    axis.line.y.right = element_line(color = "blue"),
    axis.text.y.right = element_text(color = "blue"),
    legend.position = "top",
    legend.title = element_blank())

# Rescale pH and other DIC forms for proper secondary y-axis alignment
CO2_rescaled <- scales::rescale(carb.df$CO2 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
CO3_rescaled <- scales::rescale(carb.df$CO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
HCO3_rescaled <- scales::rescale(carb.df$HCO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))

# Combine data into a single dataframe for ggplot
plot_data <- carb.df %>%
  mutate(
    CO2_rescaled = CO2_rescaled,
    CO3_rescaled = CO3_rescaled,
    HCO3_rescaled = HCO3_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC forms
  geom_line(aes(y = CO2_rescaled, color = "CO2"), linetype = "dashed", size = 2) +
  geom_line(aes(y = CO3_rescaled, color = "CO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = HCO3_rescaled, color = "HCO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = pH_calc, color = "pH"), linetype = "solid", size = 2) +

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("pH"),
    sec.axis = sec_axis(
      ~ scales::rescale(., to = range(carb.df$CO2*1e6, na.rm = TRUE)),
      name = expression("DIC form (µmol kg"^-1~")"))) +

  # Color scale and labels
  scale_color_manual(
    values = c("CO2" = "blue", "CO3" = "coral1", "HCO3" = "turquoise2", "pH" = "black"),
    labels = c(expression(CO[2]), expression(CO[3]^"-2"), expression(HCO[3]^"-1"), "pH"),
    name = "Legend"
  ) +

  # Add labels and title
  labs(x = "Salinity", 
       #title = expression(SW + "Tyroler River")
       ) +
  
  # Add theme and styling
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "black"),
    axis.text.y.right = element_text(color = "blue"),
    axis.line.y.right = element_line(color = "blue"),
    legend.position = "top",
    legend.title = element_blank()
  )

We can see that when freshwater TA:DIC exceeds that of seawater, the pH change is unidirectionally positive, and minimum pCO₂ occurs at a salinity of zero.

Set up data frame with a freshwater endmember with TA:DIC lower than seawater

# Step 2: Create the dataframe of endmembers
df <- data.frame(
  Parcel = c("River", "SW", "MQ"),
  #TA = c(194, 2198, 88),                # μmolkg−1
  #DIC = c(178, sw_dic_sat, 30.78),     # μmolkg−1
  TA = c(194.85, 2198.0433, 97.06),     # μmolkg−1
  DIC = c(225.2, 2045.953, 30.84),     # μmolkg−1
  Salinity = c(0, 31.32, 0),            # Practical Salinity Scale (PSS-78)
  Temperature = c(0, 0, 0)              # Celsius
)


df
df$TA[1]/df$DIC[1] #river
## [1] 0.8652309
df$TA[2]/df$DIC[2] # SW
## [1] 1.074337
194/1.089474
## [1] 178.0676
carb(flag = 15, var1 = df$TA[2]/1000000, var2 = df$DIC[2]/1000000, k1k2 = "l", S = 31.32, T = 0 , eos = "teos10", b = "l10")$pCO2
## [1] 263.8324
carb(flag = 15, var1 = df$TA[1]/1000000, var2 = df$DIC[1]/1000000, k1k2 = "l", S = 0, T = 0 , eos = "teos10", b = "l10")$pCO2
## [1] 397.1702
#matches atm pco2

now construct mixing model

# Function to calculate parameters of the resulting mix of River and SW
calculate_mix_parameters <- function(df, mixing_ratio) {
  mix_TA <- df$TA[1] * mixing_ratio + df$TA[2] * (1 - mixing_ratio)
  mix_DIC <- df$DIC[1] * mixing_ratio + df$DIC[2] * (1 - mixing_ratio)
  mix_Salinity <- df$Salinity[1] * mixing_ratio + df$Salinity[2] * (1 - mixing_ratio)
  mix_Temperature <- df$Temperature[1] * mixing_ratio + df$Temperature[2] * (1 - mixing_ratio)
  
  return(data.frame(
    Parcel = "Mix",
    TA = mix_TA,
    DIC = mix_DIC,
    Salinity = mix_Salinity,
    Temperature = mix_Temperature,
    x = mixing_ratio
  ))
}

# Example mixing ratios
mixing_ratios <- seq(0, 1, by = 0.001)

# Calculate mix parameters for each mixing ratio
result <- lapply(mixing_ratios, function(ratio) calculate_mix_parameters(df, ratio))

# Combine the results into a single dataframe
result_df <- do.call(rbind, result)


# Print the resulting dataframe
tibble(result_df)

Calculating carbonate system

result_df <- result_df %>% mutate(
TA = TA/1000000, #converting µmol/kg to mol/kg
DIC = DIC/1000000) #converting µmol/L to mol/kg  


carb.df <- carb(flag = 15, 
                var1 = result_df$TA, 
                var2 = result_df$DIC, 
                k1k2 = "l", 
                S=result_df$Salinity, 
                T = result_df$Temperature, 
                eos = "teos10", 
                b = "l10") #calculate pco2 from TA and DIC




carb.df <- carb.df %>% dplyr::rename(TA = ALK,
                                    pCO2_calc = pCO2,
                                    pH_calc = pH,
                                    salinity = S,
                                    temperature = T,
                                    pressure = P)

Plot

# Define fixed ranges for DIC, TA, and pCO2
DIC_range <- range(carb.df$DIC * 1e6, na.rm = TRUE)
TA_range <- range(carb.df$TA * 1e6, na.rm = TRUE)
pCO2_range <- range(carb.df$pCO2_calc, na.rm = TRUE)
pCO2_range <- c(0, 407)

# Rescale pCO2 for the plot to match the TA range
pCO2_rescaled <- (carb.df$pCO2_calc - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Rescale Tyro points using the same fixed transformation
Tyro_rescaled <- (measured_mix$Tyro - min(pCO2_range)) / diff(pCO2_range) * diff(TA_range) + min(TA_range)

# Combine data into a single dataframe
plot_data <- carb.df %>%
  mutate(
    pCO2_rescaled = pCO2_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC and TA
  geom_line(aes(y = DIC * 1e6, color = "DIC"), size = 3) +
  geom_line(aes(y = TA * 1e6, color = "TA"), size = 3) +

  # Secondary axis: pCO2
  geom_line(aes(y = pCO2_rescaled, color = "pCO2"), linetype = "dashed", size = 2) +

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("DIC & TA (µmol kg"^-1~")"), # Label for primary y-axis
    sec.axis = sec_axis(
      ~ (. - min(TA_range)) / diff(TA_range) * diff(pCO2_range) + min(pCO2_range),
      name = expression(pCO[2] ~ (µatm))
    )
  ) +

  # Color scale and labels
  scale_color_manual(
    values = c("pCO2" = "blue", "DIC" = "hotpink1", "TA" = "forestgreen"),
    labels = c("DIC", expression(pCO[2]), "TA"),
    name = "Legend"
  ) +

  # Add labels and theme
  labs(x = "Salinity") +
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "blue"),
    axis.line.y.right = element_line(color = "blue"),
    axis.text.y.right = element_text(color = "blue"),
    legend.position = "top",
    legend.title = element_blank())

We observe that when freshwater has a lower TA:DIC ratio than seawater (particularly when TA:DIC < 1), mixing produces a U-shaped pCO₂ curve associated with acidification at low salinities.

# Rescale pH and other DIC forms for proper secondary y-axis alignment
CO2_rescaled <- scales::rescale(carb.df$CO2 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
CO3_rescaled <- scales::rescale(carb.df$CO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))
HCO3_rescaled <- scales::rescale(carb.df$HCO3 * 1e6, to = range(carb.df$pH_calc, na.rm = TRUE))

# Combine data into a single dataframe for ggplot
plot_data <- carb.df %>%
  mutate(
    CO2_rescaled = CO2_rescaled,
    CO3_rescaled = CO3_rescaled,
    HCO3_rescaled = HCO3_rescaled
  )

# Create the ggplot
ggplot(plot_data, aes(x = salinity)) +
  # Primary axis: DIC forms
  geom_line(aes(y = CO2_rescaled, color = "CO2"), linetype = "dashed", size = 2) +
  geom_line(aes(y = CO3_rescaled, color = "CO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = HCO3_rescaled, color = "HCO3"), linetype = "dashed", size = 2) +
  geom_line(aes(y = pH_calc, color = "pH"), linetype = "solid", size = 2) +

  # Scale for primary and secondary y-axes
  scale_y_continuous(
    name = expression("pH"),
    sec.axis = sec_axis(
      ~ scales::rescale(., to = range(carb.df$HCO3*1e6, na.rm = TRUE)),
      name = expression("DIC form (µmol kg"^-1~")"))) +

  # Color scale and labels
  scale_color_manual(
    values = c("CO2" = "blue", "CO3" = "coral1", "HCO3" = "turquoise2", "pH" = "black"),
    labels = c(expression(CO[2]), expression(CO[3]^"-2"), expression(HCO[3]^"-1"), "pH"),
    name = "Legend"
  ) +

  # Add labels and title
  labs(x = "Salinity", 
       #title = expression(SW + "Tyroler River")
       ) +
  
  # Add theme and styling
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "black"),
    axis.text.y.right = element_text(color = "turquoise2"),
    axis.line.y.right = element_line(color = "turquoise2"),
    legend.position = "top",
    legend.title = element_blank()
  )

df
df$TA[1]/df$DIC[1] #river
## [1] 0.8652309
df$TA[2]/df$DIC[2] # SW
## [1] 1.074337

Conducting pCO₂ sensitivity Analysis

To examine the effects of glacial meltwater runoff on biogeochemical processes, we can examine these processes along a range of conditions. To do this we develop an approach that isolates the effects of individual processes—such as primary production, microbial respiration, and CaCO₃ dissolution or precipitation—along a salinity–alkalinity gradient.

We define pCO₂ sensitivity for a given process as: \[S_{pCO_2}^{\text{process}} = \frac{\Delta pCO_{2,\text{process}}}{\Delta X}\] where ΔX is the magnitude of the perturbation (e.g., +1 µmol kg⁻¹ DIC for respiration), and \[S_{pCO_2}^{\text{process}}\] has units of µatm per µmol kg⁻¹.

ΔpCO₂ is the difference between perturbed and reference states: \[\Delta pCO_{2,\text{process}} = pCO_{2,\text{perturbed}} - pCO_{2,\text{reference}}\]

This formulation captures the local (absolute) response of pCO₂ per unit biogeochemical change along a dilution gradient, accounting for changes in alkalinity and buffer capacity.

While conceptually related to the Revelle Factor—which describes relative pCO₂ sensitivity to DIC perturbations—our method differs by quantifying absolute sensitivity (µatm per µmol kg⁻¹) and incorporating changes in both DIC and TA. This enables a process-specific interpretation of carbonate system reactivity under dynamic freshwater–seawater mixing.

In practice this looks like the following:

First, we calculate the baseline DIC values required to maintain constant pCO₂ along a seawater-freshwater gradient. We use a single pCO₂ value (e.g., 402 µatm) and the linear TA variation between our seawater and freshwater endmembers to calculate the DIC variation needed to create a pCO₂ isoline in the Deffeyes diagram below. We utilize the DIC representing a single pCO₂ value in order to isolate the pCO₂ sensitivity signal from the nonlinear mixing effects that are created when freshwater enters seawater.

We then can simulate biogeochemical processes:

  1. Primary production and respiration are modeled by removing or adding 1 µmol kg⁻¹ DIC, with no change in TA
  2. CaCO₃ precipitation and dissolution are modeled by removing or adding DIC and TA in a 1:2 molar ratio

For each case, pCO₂ was recalculated using the perturbed DIC and TA, while holding salinity, temperature, and pressure constant. Finally, we can see the change in pCO₂ between the perturbed state and the reference isoline. This yields process-specific pCO₂ sensitivity values across the mixing gradient from pure seawater to pure freshwater.

Calculate DIC for a pCO₂ isoline

carb.dfsens %>% select(salinity, temperature, pressure, TA) -> carb.dfsens

carb.dfsens$DIC <- carb(flag = 24,
     var1 = 402, #if we choose this as atmospheric
     var2 = carb.dfsens$TA, 
     S = carb.dfsens$salinity, 
     T = carb.dfsens$temperature, 
     P = carb.df$pressure,  
     k1k2 = "l",  eos = "teos10")$DIC

carb.dfsens$pCO2_calc <- carb(flag = 15,
     var1 = carb.dfsens$TA, 
     var2 = carb.dfsens$DIC, 
     S = carb.dfsens$salinity, 
     T = carb.dfsens$temperature, 
     P = carb.dfsens$pressure,  
     k1k2 = "l",  eos = "teos10")$pCO2

tibble(carb.dfsens)
names(carb.dfsens)
## [1] "salinity"    "temperature" "pressure"    "TA"          "DIC"         "pCO2_calc"

Perturb system representing biological activity

carb.dfsens <- carb.dfsens %>%
  mutate(ppmodDIC = DIC - (1/1e6), # to represent primary production we reduce DIC by 1 µmol kg-1
         brmodDIC = DIC + (1/1e6)) # to represent bacterial respiration we increase DIC by 1 µmol kg-1



carb.dfsens$ppmodpco2 <-
  carb(flag = 15,
     var1 = carb.dfsens$TA, 
     var2 = carb.dfsens$ppmodDIC, 
     S = carb.dfsens$salinity, 
     T = carb.dfsens$temperature, 
     P = carb.dfsens$pressure,  
     k1k2 = "l",  eos = "teos10")$pCO2

carb.dfsens$brmodpco2 <-
  carb(flag = 15,
     var1 = carb.dfsens$TA, 
     var2 = carb.dfsens$brmodDIC, 
     S = carb.dfsens$salinity, 
     T = carb.dfsens$temperature, 
     P = carb.dfsens$pressure,  
     k1k2 = "l",  eos = "teos10")$pCO2

carb.dfsens <- carb.dfsens %>%
  mutate(ppchangepco2  =  ppmodpco2 - pCO2_calc,
         brchangepco2  =  brmodpco2 - pCO2_calc,
         TA = TA*1e6)



library(scales) # For rescale function

# Define min and max for Salinity and TA
sal_min <- min(carb.dfsens$salinity, na.rm = TRUE)
sal_max <- max(carb.dfsens$salinity, na.rm = TRUE)

ta_min <- min(carb.dfsens$TA, na.rm = TRUE)
ta_max <- max(carb.dfsens$TA, na.rm = TRUE)

(carb.dfsens %>%
  ggplot(aes(x = salinity)) +
  geom_ribbon(aes(ymin = ppchangepco2, ymax = brchangepco2), fill = "purple3", alpha = 0.3) +
  geom_line(aes(y = ppchangepco2), color = "purple3") +
  geom_line(aes(y = brchangepco2), color = "purple3") +
  labs(x = "Salinity",
       y = expression("S"[pCO[2]]^process~" (µatm per µmol kg"^-1~")")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(breaks = seq(-15, 15, by = 3)) +
    
  # Adding Secondary X-Axis for TA
  scale_x_continuous(
    sec.axis = sec_axis(
      trans = ~ rescale(., to = c(ta_min, ta_max), from = c(sal_min, sal_max)), 
      name = expression("TA (µmol kg"^-1~")")
    )
  ) +
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "black"),
    axis.title.x.top = element_text(color = "black"), # Top x-axis label color
    axis.text.x.top = element_text(color = "forestgreen"),
    axis.line.x.top = element_line(color = "forestgreen"),
    legend.position = "top",
    legend.title = element_blank())-> biosensitivityplot)

carb.dfsens <- carb.dfsens %>%
  mutate(TA = TA/1e6)

Perturb system representing mineral precipitation/dissolution

# Apply DIC and TA changes due to CaCO3 dissolution
carb.dfsens <- carb.dfsens %>%
  mutate(
    CaCO3dissDIC = DIC + (1/1e6),  
    CaCO3dissTA = TA + (2/1e6),  # Adding 1 µmol kg-1 DIC and 2 µmol -kg TA to represent CaCO3 dissolution
    CaCO3precipDIC = DIC - (1/1e6),  
    CaCO3precipTA = TA - (2/1e6)  # Subtracting 1 µmol kg-1 DIC and 2 µmol -kg TA to represent CaCO3 dissolution
  )

# Recalculate pCO2 after CaCO3 dissolution
carb.dfsens$CaCO3disspco2 <-
  carb(flag = 15,
       var1 = carb.dfsens$CaCO3dissTA, 
       var2 = carb.dfsens$CaCO3dissDIC, 
       S = carb.dfsens$salinity, 
       T = carb.dfsens$temperature, 
       P = carb.dfsens$pressure,  
       k1k2 = "l",  eos = "teos10")$pCO2

# Recalculate pCO2 after CaCO3 precip
carb.dfsens$CaCO3precippco2 <-
  carb(flag = 15,
       var1 = carb.dfsens$CaCO3precipTA, 
       var2 = carb.dfsens$CaCO3precipDIC, 
       S = carb.dfsens$salinity, 
       T = carb.dfsens$temperature, 
       P = carb.dfsens$pressure,  
       k1k2 = "l",  eos = "teos10")$pCO2

# Calculate pCO2 change due to CaCO3 dissolution
carb.dfsens <- carb.dfsens %>%
  mutate(disschangepco2  = CaCO3disspco2 - pCO2_calc,
         precipchangepco2  = CaCO3precippco2 - pCO2_calc)


carb.dfsens <- carb.dfsens %>%
  mutate(TA = TA*1e6) # to plot in better units

# Define min and max for Salinity and TA
sal_min <- min(carb.dfsens$salinity, na.rm = TRUE)
sal_max <- max(carb.dfsens$salinity, na.rm = TRUE)

ta_min <- min(carb.dfsens$TA, na.rm = TRUE)
ta_max <- max(carb.dfsens$TA, na.rm = TRUE)

# Plot
(carb.dfsens %>%
  ggplot(aes(x = salinity)) +
    
  geom_ribbon(aes(ymin = precipchangepco2, ymax = disschangepco2), fill = "navajowhite", alpha = 0.9) +
  geom_line(aes(y = disschangepco2), color = "saddlebrown") +
  geom_line(aes(y = precipchangepco2), color = "saddlebrown") +
  labs(x = "Salinity", 
        y = expression("S"[pCO[2]]^process~" (µatm per µmol kg"^-1~")")) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(breaks = seq(-15, 15, by = 3)) +
  
  # Adding Secondary X-Axis for TA
  scale_x_continuous(
    sec.axis = sec_axis(
      trans = ~ rescale(., to = c(ta_min, ta_max), from = c(sal_min, sal_max)), 
      name = expression("TA (µmol kg"^-1~")")
    )
  ) +
  theme_bw(base_size = 20) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # Black border
    panel.grid.major.x = element_blank(), # Removes major vertical grid lines
    panel.grid.minor.x = element_blank(), # Removes minor vertical grid lines
    panel.grid.major.y = element_blank(), # Keeps major horizontal grid lines
    panel.grid.minor.y = element_blank(), # Removes minor horizontal grid lines
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "black"),
    axis.title.x.top = element_text(color = "black"), # Top x-axis label color
    axis.text.x.top = element_text(color = "forestgreen"),
    axis.line.x.top = element_line(color = "forestgreen"),
    legend.position = "top",
    legend.title = element_blank()
  ) -> caco3sensitivityplot)

carb.dfsens <- carb.dfsens %>%
  mutate(TA = TA/1e6) # back to calculating units (mol kg-1)



#plot_grid(biosensitivityplot, caco3sensitivityplot, nrow = 2)

Deffeyes Diagram

Finally we illustrate a Deffeyes diagram where we can once again see that the impact of biogeochemical processes on pCO₂ is greater under freshened conditions (low TA and DIC) compared to under undiluted conditions.

# Load required libraries
library(seacarb)
library(tidyverse)
library(ggplot2)
library(reshape2)
library(metR) # For geom_text_contour()

df <- data.frame(
  Parcel = c("River", "SW", "MQ"),
  TA = c(194.85, 2198.0433, 97.06),            # μmolkg−1
  DIC = c(riv_dic_sat, sw_dic_sat, 30.84),     # μmolkg−1
  Salinity = c(0, 31.32, 0),                   # Practical Salinity Scale (PSS-78)
  Temperature = c(0, 0, 0)                     # Celsius
)# back to dataframe representing laboratory conditions

# Define constants
salinity <- 20          # Salinity in PSU
temperature <- 6       # Temperature in °C
pressure <- 0           # Pressure in dbar
k1k2 <- "l"             # Use Lueker et al. (2000) constants
eos <- "teos10"         # TEOS-10 standard

TA_values <- seq(1050, 2350, by = 5)   #  µmol/kg)
DIC_values <- seq(900, 2200, by = 5)  #  µmol/kg)




# Convert TA and DIC values to mol/kg for use with seacarb
TA_values <- TA_values / 1e6  # Convert to mol/kg
DIC_values <- DIC_values / 1e6  # Convert to mol/kg

# Create a grid of all TA and DIC combinations
grid <- expand.grid(TA = TA_values, DIC = DIC_values)

# Calculate pCO2 for each TA and DIC combination
grid$pCO2 <- apply(grid, 1, function(row) {
  result <- carb(
    flag = 15,                  # TA and DIC as input
    var1 = row["TA"],           # Total Alkalinity
    var2 = row["DIC"],          # Dissolved Inorganic Carbon
    S = salinity,               # Salinity
    T = temperature,            # Temperature
    P = pressure,               # Pressure
    k1k2 = k1k2,                # Dissociation constants
    eos = eos                   # Equation of state
  )
  result$pCO2  # Extract and return pCO2 (in µatm)
})

# Convert data to a format suitable for ggplot
grid_long <- melt(grid, id.vars = c("TA", "DIC"), measure.vars = "pCO2")

# Convert TA and DIC back to µmol/kg for plotting
grid_long <- grid_long %>%
  mutate(TA = TA * 1e6, DIC = DIC * 1e6)

grid_long <- grid_long %>% filter(value < 1010)

# Plot the Deffeyes diagram with color background and annotated isolines
ggplot(grid_long, aes(x = DIC, y = TA, z = value)) +
  # Background color gradient for pCO2
  geom_raster(aes(fill = value)) +
  scale_fill_viridis_c(name = expression(pCO[2] ~ (µatm)), option = "plasma") +
  #scale_fill_gradientn(colors = c("darkblue","goldenrod1"),  name = expression(pCO[2] ~ (µatm))) +
  # Contour lines for pCO2
  geom_contour(aes(z = value), color = "black") +
  # Annotate contour lines with pCO2 values
  geom_text_contour(aes(z = value), color = "grey90", size = 3, check_overlap = TRUE) +
  #now add vector for undiluted
  geom_segment(aes(x = df$DIC[2], y = df$TA[2], xend = df$DIC[2]+37.5, yend = df$TA[2]+75), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+ 
  geom_segment(aes(x = df$DIC[2], y = df$TA[2], xend = df$DIC[2]-75, yend = df$TA[2]+7.5), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
  geom_segment(aes(x = df$DIC[2], y = df$TA[2], xend = df$DIC[2]+75, yend = df$TA[2]-7.5), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
  #####now 25% diluted
  geom_segment(aes(x = (df$DIC[1]*0.25+df$DIC[2]*0.75), y = (df$TA[1]*0.25+df$TA[2]*0.75), 
                    xend = (df$DIC[1]*0.25+df$DIC[2]*0.75)+37.5, yend = (df$TA[1]*0.25+df$TA[2]*0.75)+75), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
  geom_segment(aes(x = (df$DIC[1]*0.25+df$DIC[2]*0.75), y = (df$TA[1]*0.25+df$TA[2]*0.75), 
                   xend = (df$DIC[1]*0.25+df$DIC[2]*0.75)-75, yend = (df$TA[1]*0.25+df$TA[2]*0.75)+7.5), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
   geom_segment(aes(x = (df$DIC[1]*0.25+df$DIC[2]*0.75), y = (df$TA[1]*0.25+df$TA[2]*0.75), 
                   xend = (df$DIC[1]*0.25+df$DIC[2]*0.75)+75, yend = (df$TA[1]*0.25+df$TA[2]*0.75)-7.5), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
  #####now 50% diluted
  geom_segment(aes(x = (df$DIC[1]*0.5+df$DIC[2]*0.5), y = (df$TA[1]*0.5+df$TA[2]*0.5), 
                    xend = (df$DIC[1]*0.5+df$DIC[2]*0.5)+37.5, yend = (df$TA[1]*0.5+df$TA[2]*0.5)+75), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
    geom_segment(aes(x = (df$DIC[1]*0.5+df$DIC[2]*0.5), y = (df$TA[1]*0.5+df$TA[2]*0.5), 
                   xend = (df$DIC[1]*0.5+df$DIC[2]*0.5)-75, yend = (df$TA[1]*0.5+df$TA[2]*0.5)+7.5), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
  geom_segment(aes(x = (df$DIC[1]*0.5+df$DIC[2]*0.5), y = (df$TA[1]*0.5+df$TA[2]*0.5), 
                   xend = (df$DIC[1]*0.5+df$DIC[2]*0.5)+75, yend = (df$TA[1]*0.5+df$TA[2]*0.5)-7.5), 
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"), color = "linen", linetype = 1)+
   #Add axis labels and title
  labs(
    x = "DIC (µmol kg"^-1~")",
    y = "TA (µmol kg"^-1~")"
  )+
  scale_x_continuous(breaks = seq(900, 2200, 200)) +
  scale_y_continuous(breaks = seq(1000, 2400, 200)) +
  # Customize plot theme
  theme_bw(base_size = 13) +
  theme(
    #panel.background = element_render(color = "grey"),
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    axis.title.y = element_text(color = "black"),
    legend.position = c(0.85, 0.2))

#for PP of 10
carb(flag = 15,
     var1 = df$TA[2]/ 1e6, 
     var2 = df$DIC[2]/ 1e6, 
     S = salinity, 
     T = temperature, 
     P = pressure,  
     k1k2 = k1k2,  eos = eos)$pCO2 - carb(flag = 15,
                                          var1 = df$TA[2]/ 1e6, 
                                          var2 = (df$DIC[2]-10)/ 1e6, 
                                          S = salinity, 
                                          T = temperature, 
                                          P = pressure,  
                                          k1k2 = k1k2,  
                                          eos = eos)$pCO2 # for pure seawater
## [1] 31.15015
carb(flag = 15,
     var1 = (df$TA[1]*0.5 + df$TA[2]*0.5)/ 1e6, 
     var2 = (df$DIC[1]*0.5+df$DIC[2]*0.5)/ 1e6,
     S = salinity, 
     T = temperature, 
     P = pressure,  
     k1k2 = k1k2,  eos = eos)$pCO2 - carb(flag = 15,
                                          var1 = (df$TA[1]*0.5 + df$TA[2]*0.5)/ 1e6, 
                                          var2 = ((df$DIC[1]*0.5+df$DIC[2]*0.5)-10)/ 1e6, 
                                          S = salinity, 
                                          T = temperature, 
                                          P = pressure,  
                                          k1k2 = k1k2,  
                                          eos = eos)$pCO2 # for 50/50 mix
## [1] 48.90052
#for CaCO3 dissolve 
carb(flag = 15,
     var1 = df$TA[2]/ 1e6, 
     var2 = df$DIC[2]/ 1e6, 
     S = salinity, 
     T = temperature, 
     P = pressure,  
     k1k2 = k1k2,  eos = eos)$pCO2 - carb(flag = 15,
                                          var1 = (df$TA[2]+10)/ 1e6, 
                                          var2 = (df$DIC[2]+5)/ 1e6, 
                                          S = salinity, 
                                          T = temperature, 
                                          P = pressure,  
                                          k1k2 = k1k2,  
                                          eos = eos)$pCO2
## [1] 13.33212
carb(flag = 15,
     var1 = (df$TA[1]*0.5 + df$TA[2]*0.5)/ 1e6, 
     var2 = (df$DIC[1]*0.5+df$DIC[2]*0.5)/ 1e6,
     S = salinity, 
     T = temperature, 
     P = pressure,  
     k1k2 = k1k2,  eos = eos)$pCO2 - carb(flag = 15,
                                          var1 = ((df$TA[1]*0.5 + df$TA[2]*0.5)+10)/ 1e6, 
                                          var2 = ((df$DIC[1]*0.5+df$DIC[2]*0.5)+5)/ 1e6, 
                                          S = salinity, 
                                          T = temperature, 
                                          P = pressure,  
                                          k1k2 = k1k2,  
                                          eos = eos)$pCO2 # for 50/50 mix
## [1] 22.34846