Introduction

The halocline is a region with a strong vertical salinity gradient, between an overlying low-salinity layer and an underlying high-salinity layer.Halocline is a dominant contributor to the density stratification and thus coincides with the pycnocline.There are no agreed-upon methodsnfor inferring halocline depth,but a reasonable method might involve locating the depth at which \(dS/dp\) ,where S is salinity and p is pressure (kelley 2014 chapter 5).Calculating the derivative using e.g. diff(S)/diff(p) can be problematic because of sensitivity to noise, especially for data that have not been bin-averaged. Achieving smoothness with conventional filtering has problems at the end-points, which is particularly troublesome for a near-surface halocline (see the next blog entry). A possible solution to such problems is to calculate the derivative with a smoothing spline.(reference) kelley

Method

Pasted below is test code that does this with the ctd dataset in the oce package.we also apply this code to another CTD file. The function returns the pressure at which the smoothing spline has highest salinity derivative, and it can also plot the results (which is recommended). The parameter named deltap is used to set the value of df (degrees of freedom) for the spline. One might think of deltap as the thickness (in dbar) of the smoothing interval for each of the sub-components of the spline. kelley

library(oce) #load the oce package

findHalocline <- function(ctd, deltap=5, plot=TRUE)
{
    S <- ctd[['salinity']]
    p <- ctd[['pressure']]
    n <- length(p)
    ## trim df to be no larger than n/2 and no smaller than 3.
    N <- deltap / median(diff(p))
    df <- min(n/2, max(3, n / N))
    spline <- smooth.spline(S~p, df=df)
    SS <- predict(spline, p)
    dSSdp <- predict(spline, p, deriv=1)
    H <- p[which.max(dSSdp$y)]
    if (plot) {
        par(mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
        plotProfile(ctd, xtype="salinity")
        lines(SS$y, SS$x, col='red')
        abline(h=H, col='blue')
        mtext(sprintf("%.2f m", H), side=4, at=H, cex=3/4, col='blue')
        mtext(sprintf(" deltap: %.0f, N: %.0f, df: %.0f", deltap, N, df),
              side=1, line=-1, adj=0, cex=3/4)
    }
    return(H)
}
  
# Plot two panels to see influence of deltap.
par(mfrow=c(1, 2))
data(ctd)

findHalocline(ctd)      
## [1] 4.292
findHalocline(ctd, 1)        

## [1] 4.721

Now,we are going to apply this formula to another CTD data to see the result.

library(oce)#load the oce package
stn0860<- read.oce("sta0860.cnv") #load the data


find_halocline <- function(stn0860,deltap=5,plot=TRUE){
  s<- stn0860[["salinity"]]
  p<- stn0860[["pressure"]]
  n<- length(p)
  
  ## trim df to be no larger than n/2 and no smaller than 3.
  N<- deltap/median(diff(p))
  df<- min(n/2,max(3,n/N))
  spline<- smooth.spline(s~p,df=df)
  SS<- predict(spline,p)
  dssdp<- predict(spline,p,deriv=1)
  H<- p[which.max(dssdp$y)]
  
  if (plot){
    par(mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
    plotProfile(stn0860, xtype="salinity")
    lines(SS$y, SS$x, col='red')
    abline(h=H, col='blue')
    mtext(sprintf("%.2f m", H), side=4, at=H, cex=3/4, col='blue')
    mtext(sprintf(" deltap: %.0f, N: %.0f, df: %.0f", deltap, N, df),
          side=1, line=-1, adj=0, cex=3/4)
  }
  return(H)
}
# Plot two panels to see influence of deltap.
par(mfrow=c(1, 2))
data(stn0860)
find_halocline(stn0860)
## [1] 85
find_halocline(stn0860, 1)

## [1] 85

reference
http://dankelley.github.io/r/2014/01/11/inferring-halocline-depth.html