sea gas exchange

The airsea package provides tools for calculating net community production (NCP) from a continuous oxygen record, if provided with appropriate companion physical data such as wind speed. Previous studies have approached this problem by observing the change in oxygen concentration between time steps ($dC$), while calculating the average rate of air-sea gas exchange over this period. The air-sea component of the change in oxygen is deducted from the observed change and the result considered NCP. This approach thus in effect assumes an instantaneous production of oxygen at the end of the time step.

Our implementation solves for NCP throughout the time step. That is to say we assume a constant rate of NCP and factor this into our air-sea gas transfer calculation. Further details and a case study can be found in @Hull2015.

In this vignette we briefly outline the implementation of the method in R and provide our validation code for O2NCP and associated functions.

library(airsea)

Simple test cases

Here we construct some fairly obvious test cases for the O2NCP function to work on.

Firstly lets try 0 change in oxygen, we should get 0 NCP if we switch off bubble supersaturation processes.

test2 = data.frame(timePeriod = 1000,
                   C0 = Csat(10, 35), C1 = Csat(10, 35) , # saturation conditions
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # stable wind (although this won't be a factor)
                   h0 = 1, h1 = 1, # shallow unchanging 1m mld
                   Pslp0 = 1015.25, Pslp1 = 1013.25 # standard  sea level pressure (again not used)
                   )
print(O2NCP.mean(test2, bubbles = F)) # expect 0

Next lets look at close to instant production at saturation conditions. We should expect a value very close, but slightly more than our $dC$. We won't exactly get our $dC$ due to outgassing of the produced oxygen during the time step.

test = data.frame(timePeriod = 1, # tiny timestep to neglect the calculation of out-gassing
                   C0 = Csat(10, 35), C1 = Csat(10, 35) + 5, # saturation conditions with 5 mmol m-3 increase
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 0.1, u1 = 0.1, # tiny wind speed to further neglect out-gassing
                   h0 = 1, h1 = 1, # unchanging 1m mld
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test, bubbles = F)) # very close, but slightly more than 5 mmol m^-3^ over time step

Now lets try negative NCP (net respiration)

test = data.frame(timePeriod = 1, # tiny time step to neglect the calculation of out-gassing
                   C0 = Csat(10, 35), C1 = Csat(10, 35) - 5, # saturation conditions with 5 mmol m-3 decrease
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 0.1, u1 = 0.1, # tiny wind speed to further neglect out-gassing
                   h0 = 1, h1 = 1, # shallow unchanging 1m mld
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test, bubbles = F)) # very close, but slightly less than -5 mmol m^-3^ over time step

Next lets turn on the bubbles again and see how they supersaturate the water, as $C1$ remains at saturation conditions, rather than bubble enhanced saturation our model predicts slightly negative NCP.

test = data.frame(timePeriod = 10*60, # 10 minute time step
                   C0 = Csat(10, 35), C1 = Csat(10, 35) , # saturation conditions
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # standard wind speed to induce some bubbles
                   h0 = 1, h1 = 1, # shallow unchanging 1m mld
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test)) # expect a bit less than 0

We can prove that this effect is due to bubbles by adjusting $C1$ to a bubble supersaturated value:

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) )
print(Csat(10, 35))
print(Cbubble)

test = data.frame(timePeriod = 10*60, # tiny timestep to neglect the calculation of out-gassing
                   C0 = Cbubble, C1 = Cbubble , # bubble correct saturation conditions
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # standard wind speed to induce some bubbles
                   h0 = 1, h1 = 1, # shallow unchanging 1m mld
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test)) # expect 0

Now lets check that the output scales correctly with large mixed layer depths. We should get a value of our $dC$ multiplied by our mixed layer depth ($h$).

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) )

test = data.frame(timePeriod = 1, # tiny timestep to neglect the calculation of out-gassing
                   C0 = Cbubble, C1 = Cbubble + 5 , # bubble correct saturation conditions with 5 mmol m-3 increase
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 1, u1 = 1, # tiny wind speed to further neglect out-gassing
                   h0 = 100, h1 = 100, # 100m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test)) # expect 500

If we express the NCP as concentration change rather than flux we should get out input dC.

print(O2NCP.mean(test, output_conc = T)) # expect 5

The amount out-gassed should increase if the conditions are further from equilibrium.

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) ) + 50 # start and end with supersturated waters

test = data.frame(timePeriod = 10*60, # 10 minute timestep
                   C0 = Cbubble, C1 = Cbubble + 5 , # bubble correct saturation conditions with 5 mmol m-3 increase
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # average wind speed
                   h0 = 100, h1 = 100, # 100m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test)) # expect a little more than 500

Indeed this means we should expect some positive NCP even with 0 change in $dC$ if we maintain supersaturation:

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) ) + 100 # start and end with supersturated waters

test = data.frame(timePeriod = 10*60, # 10 minute time step
                   C0 = Cbubble, C1 = Cbubble , # constant bubble correct saturation conditions
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # average wind speed
                   h0 = 100, h1 = 100, # 100m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test)) # expect more than 0

This should be equivilent to just the air-gas exchange:

kw("O2", 10, 7, 35, method = "NG00") * (Csat(10, 35) - Cbubble)*10*60

We should also expect net respiration if we maintain under-saturated conditions:

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) ) - 100 # start and end with under-saturated waters

test = data.frame(timePeriod = 10*60, # tiny timestep to neglect the calculation of out-gassing
                   C0 = Cbubble, C1 = Cbubble , # constant bubble correct saturation conditions
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # stable wind
                   h0 = 100, h1 = 100, # 100m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
print(O2NCP.mean(test)) # expect less than 0

Entrainment

Lets check that we get sensible measurements for a moving mixed layer.

  # conceptually:
  # 10 m-3 with concentration at saturation
v1 = 10
C1 = Csat(10, 35)
  # entrain 1 m-3 of water with concentration of 50
v2 = 1
C2 = 50
  # new concentration is total mass distributed over new volume
m1 = C1 * v1 # mass in first layer
m2 = C2 * v2 # mass in second layer
m3 = m1 + m2 # mass in final layer
v3 = v1 + v2 # volume of final layer
C3 = m3 / v3 # concentration in final layer
print(C3)

time = 1000

  # or dC is
dC = ((v3 - v1) * (C2 - C1))/v3
  # or dC per time is
(((v3 - v1) * (C2 - C1))/v3)/time
(((v3 - v1)/time) * (C2 - C1))/v3

  # final conc
C1 + ((v3 - v1) * (C2 - C1))/v3 == C3

test = data.frame(timePeriod = 1000,
                   C0 = C1, C1 = C3 , # as above
                   Cb0 = C2, Cb1 = C2 , # as above
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 0.1, u1 = 0.1, # stable wind (although this won't be a factor)
                   h0 = 10, h1 = 11, # as above, 1 meter deeper
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure (again not used)
                   )
kw("O2", 10, 0.1, 35) * (Csat(10, 35) - (C1 + C2)/2) * time
print(O2NCP.mean(test, bubbles = F, entrainment = T)) # expect 0

Numerical estimated tests

More complex tests, where starting and end conditions change, will require some numerical estimation of the true values to feed to the model.

Simple looping

Here we will solve the equation using a simple loop and very small time steps. We force the equation with a known NCP to calculate a new final oxygen concentration. This final concentration is then used by our O2NCP function, we should return our known NCP value.

loop_test <- function(test, NCP, accuracy = 1){
    # set starting conditions
  ti = 10 * accuracy # time scaling factor
  len = test$timePeriod * ti # do the calculation at scale
  C0 = test$C0 # set starting oxygen
  NCP = NCP / len # adjust NCP to rate per scale
    # calculate mean values
  k = (kw('O2', test$T0, test$u0, test$S0, method = "NG00") +
         kw('O2', test$T1, test$u1, test$S1, method = "NG00"))/2
  mld = test$h0
  Prs = ((test$Pslp0 + test$Pslp1)/2) / 1013.25  # surface pressure scaling
  B = (bubbleSat('O2', test$u0) + bubbleSat('O2', test$u1))/2
  Cs = (Csat(test$T0, test$S0) + Csat(test$T1, test$S1))/2
  dhdt = (test$h1 - test$h0)/len

  # entrainment
  if("Cb0" %in% colnames(test) & dhdt > 0){
    print("Cb found & dhdt > 0, entrainment calculated")
    Cb = (test$Cb0 + test$Cb1)/2
  }else{
      Cb = NA
    }

  Cs = Cs * Prs * (1 + B) # adjust Csat for bubbles and pressure

    # iterate
  for(i in 1:len){
      # F = k(C* - C)
    flux = k * (Cs - C0) * (1) # calculate air-sea flux
    flux = flux / ti
    if(!is.na(Cb)){
      # E = dhdt(Cb - C)
      E = dhdt * (Cb - C0)
    }else{ E = 0 }
      # dC/h = F + J
    mld = mld + dhdt
    C0 = C0 + ((NCP + flux + E)/mld) # adjust oxygen for fluxs
  }
  test$C1 = C0
  return(test)
}

lets test the model with known NCP forcings at saturation conditions:

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) )

test = data.frame(timePeriod = 60*60*1, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 10, # 10m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 5, accuracy = 5)
O2NCP.mean(test)

test = loop_test(test, -5, accuracy = 5)
O2NCP.mean(test)

test = loop_test(test, 0)
O2NCP.mean(test)

test = loop_test(test, 100)
O2NCP.mean(test)

Now lets keep NCP stable but change the enviroment:

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) )

test = data.frame(timePeriod = 60*60*1, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 12, # warming temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 10, # 10m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 5, accuracy = 5)
O2NCP.mean(test)

test = data.frame(timePeriod = 60*60*1, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 12, # warming temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 10, # 10m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 5)
O2NCP.mean(test)

test = data.frame(timePeriod = 60*60*1, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 12, T1 = 12, # warming temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 10, # 10m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 5)
O2NCP.mean(test)

test = data.frame(timePeriod = 60*60*1, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 15, T1 = 15, # warming temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 10, # 10m deep mixed layer
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 5)
O2NCP.mean(test) # expect 5

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 0.1) )

test = data.frame(timePeriod = 1, # 1 second
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 0.1, u1 = 0.1, # tiny wind speed to further reduce out-gassing
                   h0 = 10, h1 = 11, # as example
                   Cb0 = Cbubble - 50, Cb1 = Cbubble - 50 , # 
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

v1 = test$h0
C1 = test$C0
v2 = test$h1 - test$h0
C2 = test$Cb0
  # new concentration is total mass ditributed over new volume
m1 = C1 * v1
m2 = C2 * v2
m3 = m1 + m2
v3 = v1 + v2
m3 / v3 # C3 concentration
((v3 - v1) * (C2 - C1))/v3 # dC

test = loop_test(test, 0, accuracy = 1) # expect C1 to be m3 / v3
test$dC = test$C1-test$C0
print(test)
O2NCP.mean(test, entrainment = T, output_conc = T, debug = T) # expect close to 0

test = data.frame(timePeriod = 1*60*60, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 0.1, u1 = 0.1,
                   h0 = 10, h1 = 9, 
                   Cb0 = Cbubble - 10, Cb1 = Cbubble - 10 , # 
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 0)
O2NCP.mean(test, entrainment = T) # expect 0

test = data.frame(timePeriod = 1*60*60, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 0.1, u1 = 0.1,
                   h0 = 10, h1 = 15, 
                   Cb0 = Cbubble + 10, Cb1 = Cbubble + 10 , # 
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 0)
O2NCP.mean(test, entrainment = T) # expect 0

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) )

test = data.frame(timePeriod = 1*60*60, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # average winds
                   h0 = 10, h1 = 11, # as example
                   Cb0 = Cbubble - 10, Cb1 = Cbubble - 10 , # 
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )

test = loop_test(test, 10)
print(test)
O2NCP.mean(test, entrainment = T) # expect 10

Vectorisation

Lets check the function correctly deals with a data frame full of tests, which simulates a real dataset.

Cbubble = Csat(10, 35) * (1 + bubbleSat("O2", 7) )

test = data.frame(timePeriod = 1*60*60, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 11, # as example
                   Cb0 = Cbubble - 10, Cb1 = Cbubble - 10 , # 
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
test5 = data.frame(timePeriod = 1*60*60, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 11, # as example
                   Cb0 = NA, Cb1 = NA, # test NA
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
test6 = data.frame(timePeriod = 1*60*60, # 1 hour timestep
                   C0 = Cbubble, C1 = NA , # 
                   T0 = 10, T1 = 10, # stable temperature
                   S0 = 35, S1 = 35, # stable salinity
                   u0 = 7, u1 = 7, # tiny wind speed to further neglect out-gassing
                   h0 = 10, h1 = 11, # as example
                   Cb0 = 100, Cb1 = 100, # test NA
                   Pslp0 = 1013.25, Pslp1 = 1013.25 # standard  sea level pressure
                   )
test1 = loop_test(test, 10)
test2 = loop_test(test, -10)
test3 = loop_test(test, 5, accuracy = 10)
test4 = loop_test(test, 19)
test5 = loop_test(test5, 6)
test6 = loop_test(test6, -6)

test = rbind(test1, test2, test3, test4, test5, test6)
O2NCP.mean(test, entrainment = T) # expect 10, -10, 5, 19, 6, -6