A geeky post. If you haven’t a clue what I’m talking about, then apologies. If you’d like to know more, let me know and I’ll start from scratch some time.
A few months back I published a paper in the Open Access journal, Ocean Science, which described how to calculate the air-sea transfer velocity for any gas. It’s long and pretty tedious but brings together a whole bunch of stuff to enable other people to use the scheme to get the numbers they need without having to spend months trawling the literature for solubility data, diffusion coefficients, the equation of state for seawater etc. So all in all I hope it’s a useful contribution to science.
Along with the paper was an implementation of the scheme in R which is free for anybody to use. Since then I’ve put a copy up on github under an open source license. AS well as the core scheme from the paper there is a set of R functions which allow the reproduction of all the figures in the paper. I will also add extensions to the scheme into the github repository as I deem them useable by the community and/or publish them in future papers. Coming soon: chemical enhancement in the liquid and gas phases. Very exciting.
A few people have recently taken an interest in using the program to calculate air-sea trace gas fluxes for their gases of interest so I thought it would be useful to document how you might do that, starting with a spreadsheet containing temperature, salinity, wind speed and concentrations in both phases of the gas in question on each row, when one row equates to one flux calculation.
Here’s an example of some input data. It’s not very realistic a) because all the data is generated by picking random numbers from a gaussian distribution – most environmental data is skewed, especially concentration data and b) I don’t have time to check what a realistic concentration range for CH3I in atmosphere or ocean is. If someone would like to populate that spreadsheet with more realistic data in the future that’d be grand.
Anyway, step one is to generate a comma spaced value (.csv) file from your spreadsheet. In the case of the above you can get a csv version here, or in the case of your own data in (god forbid) Excel or some other spreadsheet, there’ll be a save as csv option.
So you’ve got your CSV file. Here I’ll assume that you have a working install of R pointed at your working directory which already contains the files from the github repository which make up the transfer velocity scheme. To test if the latter is the case, type (at the R console)
>source("K_calcs_Johnson_OS.R")
If nothing happens (i.e. R doesn’t complain that it has no idea what you’re talking about) then you’re good to go. Put your CSV file in the same directory (or somewhere else of your choosing if your happy to manipulate paths etc).
Now, to load that file in:
>my_data<-read.csv("Exampleflluxcalcinputs.csv")
read.csv assumes that the first row is headings for each of the columns so this has worked out of the box. If your data is formatted differently you might need to use read.table, or specify header=FALSE in the command. For more info type
>help("read.csv")
Anyway, the data is in and all is good at my end so we’ll press on.
In an R data frame you can get to each column of data in a number of ways. e.g.
> my_data[,3]
[1] 35.76926 35.36087 35.61887 35.35569 35.78936 35.47348 35.67157 35.45667
[9] 35.38169 35.77596 35.13184 35.57471 35.74122 35.60888 35.26568 34.96617
[17] 35.24797 35.19071 35.32195 35.34880
> my_data[,"Salinity"]
[1] 35.76926 35.36087 35.61887 35.35569 35.78936 35.47348 35.67157 35.45667
[9] 35.38169 35.77596 35.13184 35.57471 35.74122 35.60888 35.26568 34.96617
[17] 35.24797 35.19071 35.32195 35.34880
for convenience, let’s define some variables which are easy to remember:
> sal<-my_data[,"Salinity"]
> temp<-my_data[,"Temperature"]
> wind<-my_data[,"Windspeed"]
> sw_concn<-my_data[,6]
> atm_concn<-my_data[,5]
Simples. Next we need to use the K_calcs… program to calculate the transfer velocity for each compound.
make sure it’s loaded:
>source("K_calcs_Johnson_OS.R")
The function Kw calculates total transfer velocity for a gas with respect to a water phase concentration gradient.
You can see what a function does in R by calling it with no parentheses:
> Kw
function(compound,T,u,S){
#calculate total transfer velocity (Kw) after Liss and Slater
1/(rg(compound,T,u,S)+rl(compound,T,u,S))
}
You can see from the above that the transfer velocity scheme is a series of nested function with the key parameters [compound, temperature, wind speed etc] all passed around between them. you can trace these through by looking at each function in turn in the same way as we did above for Kw. In fact you can call any function directly so if, e.g. you need the diffusivity of ammonia in air at 30 degrees you can call:
> D_air("NH3",30)
[1] 0.2111959
(diffusivities are in cm2/s)
Back to Kw. To calculate a single Kw value (e.g. for CH3I at 20 celcius, 5m/s windspeed and S=35) you’d do this:
> Kw("CH3I",20,5,35)
[1] 1.550430e-05
Which is great, but we want to calulate Kw for all the rows in our data at once. This is where R is particularly powerful. We can define a new column in my_data and use values in the other columns to calculate the value for each row of this new column, all at once!:
> my_data[,"Kw"]<-Kw("CH3I",temp,wind,sal)
my_data[,"Kw"]
[1] 7.708174e-06 5.370894e-05 8.074108e-05 2.681450e-05 1.489535e-05
[6] 3.510663e-05 3.820262e-05 5.879458e-05 4.813323e-05 1.059725e-04
[11] 4.354995e-05 2.856864e-05 1.078323e-05 1.036024e-05 4.048401e-06
[16] 3.977783e-05 6.687053e-05 2.779482e-05 3.426026e-05 3.584965e-07
Sweet huh?
So next let’s add a column which is the dimensionless henry’s law constant for CH3I at the temperature in each row:
> my_data[,"H"]<-KH("CH3I",temp,sal)
> my_data[,"H"]
[1] 0.1660876 0.1644287 0.1858428 0.1836445 0.1487224 0.1621694 0.1746626
[8] 0.2077991 0.1671034 0.1882404 0.2030620 0.1269990 0.1507726 0.1790700
[15] 0.1742638 0.1406145 0.1847995 0.1601485 0.2361064 0.1669096
Finally we can add a column which is the flux calculation, according to the flux equation:
F = Kw(Ca/H – Cw)
where F is the flux, Ca is the air concentration, Cw the water concentration, H the Henry’s law constant. [make sure all your units are correct -easy to go wrong here. The way I’m laying it out here is one correct way, but there are many others. Feel free to get in touch if you get stuck…]
In the case below I’ve multiplied by minus one to make a positive value of F represent a flux from sea to air – the way I like it ;-):
> my_data[,"F"]<-(my_data[,"Kw"]*((atm_concn/my_data[,"H"])-sw_concn))*(-1)
>my_data[,"F"]
[1] 3.415365e-05 2.715205e-04 2.514285e-04 2.717888e-04 7.330746e-05
[6] 4.749550e-05 3.759730e-06 5.432819e-04 2.955133e-04 8.373667e-04
[11] 1.502132e-04 9.293719e-05 7.623237e-05 5.816068e-05 2.731569e-05
[16] 2.445745e-04 2.239382e-04 -1.252216e-05 3.216200e-04 6.304133e-07
* disclaimer: this has been done in a rush in the half hour before going away on holiday for a week – please check carefully before you rely on this. I’ll check it in detail when I get back.
You can then tell R to write you a nice CSV of all the data:
>write.csv(my_data,"flux_calcs.csv",row.names=FALSE)
Download it here
Sorry about the ridiculous number of ‘significant’ figures throughout. R can be tamed in that respect but I haven’t had time. Example here.
Bob’s your uncle as they say.