1 Introduction

Over the ocean, direct measurements of turbulent fluxes are notoriously difficult because of platform motion, flow distortion, and sea spray. Air–sea exchanges of momentum, heat and moisture are thus usually derived using indirect methods. Flux–gradient relationships are required in multiple indirect methods such as the bulk aerodynamic, profile, and inertial dissipation methods that calculate the fluxes from the mean, profile, and high-frequency spectral measurements, respectively (Stull 1988). These flux–gradient relationships are also extensively used in numerical weather and climate numerical models to provide the lower boundary conditions for the atmospheric flow. Consequently, accurate parametrization of flux–gradient relationships is of crucial importance for capturing and simulating air–sea interaction.

The most commonly used flux–gradient relationships are based on the Monin–Obukhov similarity theory (hereafter MOST, see Monin and Obukhov 1954). MOST assumes that the non-dimensional vertical gradients of velocity, temperature, and specific humidity are universal functions of the atmospheric stability parameter \(\zeta \) (Stull 1988; Foken 2008). In this study, we focus on the scalar flux–gradient relationships

$$\begin{aligned} \left( {\frac{\kappa z}{\theta _*}} \right) \frac{\partial \overline{\theta }}{\partial z}&= \phi _T (\zeta ) \end{aligned}$$
(1)
$$\begin{aligned} \left( {\frac{\kappa z}{q_*}} \right) \frac{\partial \overline{q} }{\partial z}&= \phi _q (\zeta ) \end{aligned}$$
(2)

where the von Karman constant \(\kappa = 0.4, z\) is the height above the surface, \(\theta \) is the potential temperature and \(q\) is the specific humidity; \(\theta _*\) and \(q_*\) are the surface scaling parameters for potential temperature and specific humidity, respectively, with \(\theta _*=-{\overline{{w}^{\prime }{\theta }^{\prime }} }/{u_*}\) and \(q_*=-{\overline{{w}^{\prime }{q}^{\prime }}}/{u_*}\). Here, \(u_*=\left( {\overline{{u}^{\prime }{w}^{\prime }} ^{2}+\overline{{v}^{\prime }{w}^{\prime }} ^{2}} \right) ^{1/4}\) is the friction velocity, where \(u\) is the streamwise velocity component, \(v\) is the cross-stream velocity component and \(w\) is the vertical velocity component. Also, \(\zeta = z/L\) where \(L\) is the Obukhov length,

$$\begin{aligned} L=\frac{\overline{\theta _v}}{g\kappa }\frac{u_*^2}{\left( {\theta _*+0.61\overline{\theta _{vs}} q_*} \right) }, \end{aligned}$$
(3)

where \(g\) is the acceleration of gravity, \(\theta _{v}\) is the virtual potential temperature and \(\theta _{vs}\) is the virtual potential temperature at the surface. The overbar indicates the Reynolds-averaged values and the primes denote the turbulent perturbations. While these flux–gradient relationships are assumed to be universal, the exact functional forms of these relationships cannot be predicted by MOST (Stull 1988; Foken 2008). There have been new studies exploring the functional forms of these flux–gradient relationships from a theoretical perspective (Katul et al. 2011, 2013; Li et al. 2012b), which still rely significantly on experimental observations. The functional forms of these flux–gradient relationships have been the topic of many field campaigns, including the famous Kansas experiment (Dyer and Hicks 1970; Businger et al. 1971) over land and the Coupled Ocean-Atmosphere Response Experiment (COARE) (Fairall et al. 1996, 2003) over ocean. Numerous other experimental studies (Dyer 1974; Dyer and Bradley 1982; Högström 1988, 1996; Kader and Yaglom 1990; Brutsaert 1992; Rutgersson et al. 2001; Edson et al. 2004; Larsen et al. 2004; Smedman et al. 2007) have generated and assessed a range of empirical and semi-empirical functions (see Li et al. 2012b for a review), among which the most widely used are the Businger–Dyer relationships from the Kansas experiment (Businger 1988). The Businger–Dyer relationships for scalars under unstable conditions are:

$$\begin{aligned} \phi _{Kx} =\left( {1-\gamma _{kx} \zeta } \right) ^{-1/2} \end{aligned}$$
(4)

where \(\gamma _{kx}\) is an empirical constant (\(x =T, q\) for temperature and humidity, respectively). The universal function for moisture is usually assumed to be equal to that for temperature (i.e., \(\gamma _{kq}=\gamma _{kT}\)).

The Businger–Dyer relationships were only validated for a limited range of atmospheric stabilities (\(0.1< - \zeta < 2\), see Grachev et al. 2000). In the free convection limit (\(\zeta \rightarrow -\infty \)), due to the fact that \(u_*\) is no longer dynamically important (Priestley 1955), the universal functions are required to take the convective form (Carl et al. 1973; Wilson 2001; Li et al. 2012b),

$$\begin{aligned} \phi _{Cx} =\left( {1-\alpha _{Cx} \zeta } \right) ^{-1/3} \end{aligned}$$
(5)

where \(\alpha _{cx}\) is an empirical constant (\(x =T, q\) for temperature and humidity, respectively). Similarly, \(\alpha _{cq}\) is usually assumed to be identical to \(\alpha _{cT}\). In order to retain a function that covers the whole range of instability, Fairall et al. (1996) used an interpolation between the integrals of \(\phi _{kx}\) and \(\phi _{cx}\), which are expressed as

$$\begin{aligned} \psi _{kx,cx} (\xi )=\int _0^\xi {\frac{1-\phi _{kx,cx} (\xi )}{\xi }\text{ d }\xi }, \end{aligned}$$
(6)

and the interpolated function is

$$\begin{aligned} \psi _x (\zeta ) =\frac{\psi _{Kx} (\zeta ) +\zeta ^{2}\psi _{Cx} (\zeta )}{1+\zeta ^{2}}. \end{aligned}$$
(7)

This interpolated function \(\psi _{x}\) forms the basis of the COARE bulk algorithm. The newest COARE 3.0 algorithm adopts \(\gamma _{kx}= 15\) and \(\alpha _{cx}= 34.15\) in the \(\phi _{kx}\) and \(\phi _{cx}\), respectively (see Fairall et al. 2003). Edson et al. (2004) suggests that \(\gamma _{kx}= 13.4\) and \(\alpha _{cx} = 30\) should be used based on their measurements.

Grachev et al. (2000) report that only specific choices of the \(\alpha _{cx}\) coefficients in the COARE algorithm can result in smooth and monotonically decreasing functions for the interpolated specific choices of the \(\alpha _{cx}\) coefficients (Grachev et al. 2000). To relax the restrictions on the choices of \(\alpha _{cx}\), Akylas and Tombrou (2005, hereafter AT 2005) proposed another function that is based on a direct interpolation between \(\phi _{kx}\) and \(\phi _{cx}\), viz.

$$\begin{aligned} \phi _x (\zeta )=\frac{c^{2}\phi _{Kx} (\zeta )+\zeta ^{2}\phi _{Cx} (\zeta )}{c^{2}+\zeta ^{2}} \end{aligned}$$
(8)

where \(c\) is a constant (= 1 as suggested in AT2005). In this study, the AT2005 approach is adapted in order to parametrize the scalar flux–gradient relationships under a wide range of instabilities. The two different interpolations between the standard Businger–Dyer relation and the free convection form (i.e., the COARE bulk algorithm and the AT2005 scheme), together with two other parametrizations for the scalar flux–gradient relationships, are also inter-compared.

The similarity assumption that temperature and specific humidity share the same flux–gradient relationship has not received much attention within the context of air–sea interactions, despite an increasing understanding of the dissimilarity between the turbulent transport of heat and water vapour (see e.g., Detto et al. 2008; Cava et al. 2008; Katul et al. 2008; Cancelli et al. 2012; Li et al. 2012a). In this study, we aim to investigate the dissimilarity between the flux–gradient relationships of temperature and specific humidity in the marine surface layer, particularly under unstable conditions (i.e., \(\zeta <0\)) that favour the occurrence of scalar dissimilarity (Dias and Brutsaert 1996; Park et al. 2009; Cancelli et al. 2012; Li et al. 2012a). The paper is organized as follows: Sect. 2 introduces the experimental site and the data processing procedures. Section 3 presents the methodology, followed by results in Sect. 4, while Sect. 5 discusses the implications of the study.

2 Experimental Set-Up and Data Processing

2.1 Site and Instruments

The Bohe marine meteorological platform was constructed in a coastal zone of the South China Sea, as shown in Fig. 1. The platform is located at \(21.44^{\circ }\)N, \(111.39^{\circ }\)E, about 1.1 km east of the small island of Dazhuzhou. The top of the platform is approximately 36 m above the water surface, as shown in Fig. 2, with the mean water depth about 14 m. A wave monitor system \(\text{ WaMoS }^\circledR \) is located on the shore and its coverage is shown in Fig. 1.

Fig. 1
figure 1

A map of the coastal marine platform. The sector denotes the coverage of the wave monitor system \(\text{ WaMoS }^{\circledR }\) on the shore

Fig. 2
figure 2

A schematic view of the Bohe marine meteorological platform

Figure 2 also shows an overview of the instruments available from the platform: sensors for obtaining the profiles of mean wind speed, temperature and relative humidity are installed at five levels above the mean sea surface: at 13.4, 16.4, 20, 23.4 and 31.3 m. Two eddy-covariance systems, each including a three-dimensional sonic anemometer (Campbell Scientific CSAT3 or Gill WindMaster Pro) and an optical water vapour/\(\text{ CO }_{2}\) sensor (LICOR-7500), are installed at heights of 27.3 and 35.1 m. The sonic anemometer boom points to the east (i.e., \(90^{\circ }\) from the north). The atmospheric pressure is also measured at three levels: 13.4, 27.3 and 35.1 m, and atmospheric pressure at other heights is calculated by linear interpolation using the measured pressure at these three levels. All instruments are mounted on booms of length about 2 m. Detailed information about the instruments is presented in Table 1.

For the following analyses, the potential temperature is calculated using the air temperature and air pressure, and the specific humidity is calculated from the relative humidity, air temperature, and the air pressure (Brutsaert 2005).

Table 1 Basic information about the measurements

2.2 Turbulent Fluxes

Turbulent fluxes are calculated based on the eddy-covariance method and are denoted positive when directed upwards. In our study, 30 minutes is chosen to be the averaging interval for the following two reasons: first, it is the most commonly used averaging time interval; second, the Ogive test (see e.g., Oncley et al. 1996) also suggests that this time scale contains contributions of eddies of almost all sizes (not shown). In order to obtain high quality turbulence data, the time series are scrutinized and some basic quality controls are implemented. Further, statistical tests were performed to ensure stationarity in the time series. In addition, the fluxes at two levels are compared to validate the assumption of constant flux in the surface layer and to exclude the influence of horizontal advection. Finally, the footprints of the measurements are also analyzed.

2.2.1 Basic Tests

For the basis analyses of turbulent time series, three tests suggested by Vickers and Mahrt (1997) were employed: the absolute limit test, spike test, and higher-moment test. The absolute limit test is used to exclude unrealistically large or small values. The ranges of the absolute values of different variables are: [0, 30 m \(\text{ s }^{-1}\)] for the horizontal wind speed, and [0, 5 m \(\text{ s }^{-1}\)] for the vertical wind speed, [0, \(50^{\circ }\)C] for the air temperature, and [0, 50 g \(\text{ m }^{-3}\)] for the water vapour density. Any 30-min segment that has more than 1 % of data points outside the pre-defined range is excluded from the analyses.

In the spike test, the standard deviation of the differences between consecutive data points of a 30-min record is calculated. A point is considered as a “spike” if its corresponding time-differential values lie beyond the threshold of six times the standard deviation, and the discarded spike point is then linearly interpolated with its neighbourhood points. When more than four points continue to lie outside the threshold, they are no longer identified as spikes. If the number of spikes detected exceeds 1 % of the total number of data points, the segment is also excluded from the analysis. Note that the spike test was performed for a non-overlapping window.

The ranges of some second- to fourth-order moments are also pre-defined. The data segments are excluded if the skewness is outside of the range [\(-2\), 2] or the kurtosis is outside of the range [1, 8]. The record is also discarded when the standard deviations of the following variables are outside of the specified ranges: [0.01, 4] for the horizontal wind speed, [0.01, 3] for the vertical wind speed, [0.01, 0.5] for the air temperature, and [0.01, 0.5] for the water vapour density. These limits are empirically chosen based on e.g. Vickers and Mahrt (1997) and Park et al. (2009). The analyses show that the data points discarded through these higher order moment tests are \(<\)4.6 % of the total data.

2.2.2 Statistical Tests

Since the flux–gradient relationships are based on MOST, which requires that the turbulence flow is well developed and in steady state, two statistical metrics [non-stationarity test (NST) and integral turbulence characteristic (ITC)] are calculated (Foken and Wichura 1996).

The NST is assessed by comparing the 30-min scalar flux to the average of six consecutive 5-min scalar fluxes:

$$\begin{aligned} \textit{NS}T_x =\left| {\frac{\overline{{w}^{\prime }{x}^{\prime }_{30}} -\frac{1}{6}\sum _{i=1}^6 {\left( {\overline{{w}^{\prime }{x}^{\prime }_5} } \right) _i} }{\overline{{w}^{\prime }{x}^{\prime }_{30}} }} \right| , \end{aligned}$$
(9)

where \(x = u, T\), and \(q\). The ITC test is also employed to qualify the development of turbulence using the flux–variance method (Foken and Wichura 1996). The ITC is calculated as:

$$\begin{aligned} {\textit{ITC}}_y =\left| {\frac{\left( {\sigma _y/Y_*} \right) _{\mathrm{modeled}} -\left( {\sigma _y/Y_*} \right) _\mathrm{measured}}{\left( {\sigma _y/Y_*} \right) _{\mathrm{modeled}}}}\right| , \end{aligned}$$
(10)

where \(\sigma _{y}\) is the standard deviation, \(Y_{*}\) denotes the surface scaling parameter and \(y = u\) (i.e., the ITC test is applied only to the momentum flux). Records with the NST or ITC \(>\) 30 % are excluded. However, the validity of those functions previously determined (Lee et al. 2004, pp. 191–193) is not well validated under extreme instabilities, and runs with \(\zeta < -2\) are thus free of the ITC test.

2.2.3 Flux Corrections

The 10-Hz eddy-covariance data are corrected by the following procedures prior to flux calculation. Following the recommendations of Lee et al. (2004, pp. 33–66), the coordinate system is aligned with the mean flow streamlines using the planar fit method (Wilczak et al. 2001). The rotation is performed using measurements from October 27 to November 26, 2009, since the positions of the two sonic anemometers do not change significantly during the analysis period. More than 1,400 segments are used to obtain the regression plane and the resulting pitch and roll angles are about 2\(^{\circ }\) and 3\(^{\circ }\), respectively.

To obtain sensible heat fluxes, the sonic temperature is converted into the actual temperature (Kaimal and Gaynor 1991); a density correction (Webb et al. 1980) is applied to water vapour fluxes. Frequency response corrections are also applied to raw fluxes accounting for low-pass (sensor separation, dynamic frequency sensor response, scalar and vector path averaging, frequency response mismatch) and high-pass (blocking averaging) following (Moore 1986). The above processes are carried out using EdiRe developed by the University of Edinburgh (http://www.geos.ed.ac.uk/abs/research/micromet/EdiRe/). After these corrections, the friction velocity decreases 0.0076 m s\(^{-1}\) on average, ranging from \(-\)0.069 to 0.049 m s\(^{-1}\). The sensible heat fluxes decrease by 6.7 W m\(^{-2}\) on average with the changes ranging from \(-\)20 to 10 W m\(^{-2}\). The latent heat fluxes typically increase, with the increases 8 W m\(^{-2}\) on average and 23 W m\(^{-2}\) under most under unstable conditions.

2.2.4 Comparison Between Turbulent Fluxes at Two Levels

In order to exclude the effects of horizontal advection, the constant flux assumption is tested by comparing the fluxes measured at the two levels. Data segments with differences between fluxes at the two levels larger than 10 % are excluded. Among the 653 segments that pass the basic tests and quality control (see Sect. 2.2.1), there are 589 and 264 segments for sensible heat fluxes and latent heat fluxes, respectively, that satisfy the constant flux assumption, as shown in Fig. 3.

Fig. 3
figure 3

A comparison between latent (LE) and sensible heat fluxes (\(H\)) at \(z = 35.1\) and 27.3 m. Here, \(LE= \rho \lambda w^{{\prime }}q^{{\prime }}, H=\rho c p w^{{\prime }} \theta ^{{\prime }}\), \(\rho \) is the air density, \(\lambda \) is the latent heat of vaporization, \(c_{p}\) is the specific heat at constant pressure

2.2.5 Footprint Analysis

The concepts of fetch and the internal boundary layer are utilized in this study to examine the primary source area that affects the turbulent flux measurements (Garratt 1990). At this site, the overwater fetch is about 10 km, and nearly unlimited for wind directions from the north-east to south-west (\(045^{\circ }\)\(250^{\circ }\)). In other directions, the fetch is about 6–10 km and there are three small islands. In the unstable offshore flow regime, with the approximated 1:30 height-to-fetch ratio (Garratt 1990; Vickers and Mahrt 1997), the estimated internal boundary-layer heights vary from 200 to 330 m, so the whole tower is within the internal boundary layer and the influence of the land should be minimal. However, as illustrated in the map and mentioned above, there are three small islands near the tower that might affect the measurements. To quantify this effect, the flux footprint is estimated using an analytical footprint model (Kormann and Meixner 2001; Neftel et al. 2008). Since northerly winds are prevalent as can be seen from Fig. 5, Fig. 4 shows two representative footprints that are calculated with wind directions \(340^{\circ }\) and \(020^{\circ }\) from the north, respectively. The footprints are mainly composed of water surfaces. For wind directions in the range from \(340^{\circ }\) to \(020^{\circ }\) (clockwise), no less than 89 % of the footprint area is within the rectangle shown on Fig. 4, of which the contribution from the small island is no more than 1 %. It is thus concluded that the influence of small islands can also be neglected for our purpose.

Fig. 4
figure 4

A schematic of the footprint under unstable conditions. The two ellipses represent the footprints (\(f > 0.01 f_{\max }\), where \(f\) is the footprint function, as in Kormann and Meixner 2001) for wind directions \(340^{\circ }\) and \(020^{\circ }\) from the north, respectively. The rectangle is of size 2.7 km \(\times \) 3.8 km

Fig. 5
figure 5

Time series of (a) pressure, (b) temperature, (c) relative humidity, (d) wind speed, and (e) wind direction at 13.4 m. The sea-surface temperature is also shown in (b) for comparison

2.2.6 Data Selection

In addition to the selection criteria described above, small fluxes are excluded using the following threshold values: 10 and 15 W m\(^{-2}\) for sensible heat flux and latent heat flux, respectively. Disturbed turbulence conditions downwind of the tower (wind direction \(270^{\circ } \pm 30^{\circ }\) from the north) are also excluded. After the above selection, 194 and 234 segments are retained for the sensible heat and latent heat fluxes, respectively.

2.3 Vertical Gradients

In addition to turbulent momentum, sensible and latent heat fluxes, the vertical gradients of temperature and specific humidity are another important component in the flux–gradient relationships. Many methods have been used to calculate the vertical gradients in the literature and it has been shown that the results are not insensitive to the methods (e.g., see Park et al. 2009). In this study, two methods are used to calculate the vertical gradients of temperature and specific humidity. First, the vertical gradients are obtained by fitting the following log-square function to the profiles:

$$\begin{aligned} X(z)=a\ln ^{2}(z)+b\ln z+c, \end{aligned}$$
(11)

where \(X(z)\) represents either the potential temperature \(\theta \) or the specific humidity \(q\) at measurement height \(z\). The gradients are then determined by evaluating the differentiation of the function \(X(z)\) with respect to \(z\). The profiles of temperature and specific humidity are scrutinized in order to ensure the quality of the resulting vertical gradients. The specific humidity measurements at 23.4 m are constantly biased and are excluded in the calculation of vertical gradients. The temperature gradients are calculated at all five levels, i.e., 13.4, 16.4, 20, 23.4 and 31.3 m. Second, the vertical gradients are simply computed based on the temperature and specific humidity at the lowest level (13.4 m) and the highest level (31.3 m) using the first-order approximation. This simple approach is chosen to maximize the differences in temperature and specific humidity in calculating vertical gradients and thus to reduce the influence of measurement imprecision. The log-square fitting method and the simpler first-order approximation method are termed “method 1” and “method 2” hereafter. Given that method 2 only considers the differences in temperature and specific humidity between the lowest and highest level, it is expected that method 1 should then be more representative and yield a more realistic form of the profiles.

3 Methodology

As mentioned earlier, the function used in the COARE 3.0 algorithm is based on an interpolation between the integrals of \(\phi _{kx}\) and \(\phi _{cx}\) (Fairall et al. 1996; Grachev et al. 2000; Edson et al. 2004). The resulting \(\phi _{x}\) from Eq. 8 is:

$$\begin{aligned} \phi _x (\zeta )=\frac{\phi _{Kx} (\zeta ) +\zeta ^{2}\phi _{Cx} (\zeta )}{1+\zeta ^{2}}+2\zeta ^{2}\frac{\Psi _{Kx} (\zeta ) -\Psi _{Cx} (\zeta )}{({1+\zeta ^{2}})^{2}}, \end{aligned}$$
(12)

where \(\phi _{kx}\) and \(\phi _{cx}\) are given by

$$\begin{aligned} \phi _{Kx}&= ({1-15\zeta })^{-1/2},\end{aligned}$$
(13)
$$\begin{aligned} \phi _{Cx}&= \left( {1-34.15\zeta } \right) ^{-1/3}, \end{aligned}$$
(14)

and their integrals \(-\psi _{kx}\) and \(-\psi _{cx}\) are

$$\begin{aligned} \psi _{Kx}&= 2\ln \left( {\frac{1+\phi _{Kx}}{2}} \right) ,\end{aligned}$$
(15)
$$\begin{aligned} \psi _{Cx}&= 1.5\ln \left( {\frac{1+\phi _{Cx} +\phi _{Cx}^2}{3}} \right) -\sqrt{3}\tan ^{-1}\left( {\frac{1+2\phi _{Cx}}{\sqrt{3}}} \right) +\frac{\pi }{\sqrt{3}}. \end{aligned}$$
(16)

The other scheme proposed by Akylas and Tombrou (2005) is a direct interpolation between \(\phi _{kx}\) and \(\phi _{cx}\) (see Eq. 8). In this study, the coefficients \(\gamma _{kx}\) and \(\alpha _{cx}\) that are used in the AT2005 scheme are not specified; instead, they are to be determined using experimental data.

Two other parametrizations are also examined for comparison: the medium range forecast (MRF) scheme proposed by Hong and Pan (1996, hereafter referred to as MRF) and that of Park et al. (2009, hereafter P2009). The MRF scheme uses the following function for both temperature and water vapour under unstable conditions:

$$\begin{aligned} \phi _x (\zeta ) =\left( {1-16\zeta } \right) ^{-1/2}. \end{aligned}$$
(17)

In P2009, the functions for temperature and water vapour under unstable conditions are

$$\begin{aligned} \phi _T (\zeta )&= ({1-13.3\zeta })^{-1/2},\end{aligned}$$
(18)
$$\begin{aligned} \phi _q (\zeta )&= 1.21 ({1-13.1\zeta })^{-1/2}. \end{aligned}$$
(19)

Following Salesky et al. (2012) and Salesky and Chamecki (2012), the filtering method is applied to the time series of the instantaneous fluxes (e.g., \(u^{\prime }w^{\prime }\)), and the standard statistical methods are applied to the mean temperature and humidity profiles. After estimating the random errors in the turbulent quantities (e.g., \(u_{*}\)) and in the gradients \((\frac{\partial \bar{\theta }}{\partial z}, \frac{\partial \bar{q}}{\partial z})\), the errors in \(\zeta , \phi _{T}\) and \(\phi _{q}\) can be derived based on error propagation:

$$\begin{aligned} \varepsilon _{\phi _x} =\left( {\varepsilon _{u_*}^2 +\varepsilon _{\overline{w{\prime }x{\prime }}}^2 +\varepsilon _{\partial \overline{x}/\partial z}^2} \right) ^{1/2}=\left( {\frac{\varDelta _{u_*}^2}{u_*^2}+\frac{\varDelta _{\overline{w{\prime }x} {\prime }}^2}{\overline{w{\prime }x{\prime }} ^{2}}+\frac{\varDelta _{\partial \overline{x}/\partial z}^2}{\left( {\partial \overline{x}/\partial z} \right) ^{2}}} \right) ^{1/2} \end{aligned}$$
(20)

where \(\varepsilon \) denotes the relative error and \(\varDelta \) is the random error.

4 Results

4.1 Surface Meteorological Conditions and Air–Sea Fluxes

Two periods are selected in our study: one extends from November 13 to November 21, 2009 and the other from December 16 to December 20, 2009. The study periods include cold-air outbreak events during which cold, dry air masses move over warm coastal waters and relatively large air–sea temperature and specific humidity differences persist, with large sensible and latent heat fluxes that further create intense boundary-layer mixing and moist/forced convection (Bane and Osgood 1989; Grossman and Betts 1990; Dorman et al. 2004).

The mean meteorological conditions at the site during one of the study periods (from November 13 to November 21, 2009) are shown in Fig. 5. From November 16 to November 18, the pressure and wind speed increase while the air temperature and specific humidity decrease, implying a cold-air outbreak event. The mean flow during the study periods is predominantly from the north and the mean wind speeds range from 5 to 15 m \(\text{ s }^{-1}\), except during intermittent periods when the wind speeds are less than 2 m \(\text{ s }^{-1}\). The mean air temperature, sea surface temperature (SST) and relative humidity show no clear diurnal cycles, due to the strong influence of the cold-air outbreak. Air temperature ranges from 9 to \(20^{\circ }\)C and SST ranges from 10 to \(21 ^{\circ }\)C. After the onset of the cold-air outbreak, the SST is consistently higher than the air temperature. The difference between the surface and air temperatures reaches its maximum around November 18.

Figure 6 shows the turbulent fluxes (the friction velocity, sensible and heat fluxes) and the stability parameter \(\zeta \); the typical sensible heat fluxes are small at the air–sea interface (Mahrt et al. 1998). Nonetheless, during the cold-air outbreak event, the sensible heat flux reaches its maximum because of the large difference between the surface and air temperatures. The stability parameter \(\zeta \) is shown in Fig. 6d, with the negative stability parameter suggesting that unstable conditions prevail at the measurement site. Due to the large sensible and latent heat fluxes, as well as the variations in the friction velocity, the marine surface layer sometimes approaches free-convection conditions, which provides the basis for evaluating the scalar flux–gradient relationship under a wide range of instability conditions. Note that some segments (mainly during November 15–November 17) show countergradient latent heat fluxes under unstable conditions at both heights (27.3, 35.1 m), and such runs are discarded from the following analysis of the flux–gradient relationship.

Fig. 6
figure 6

Time series of a latent heat flux, b sensible heat flux, c friction velocity and d the stability parameter

4.2 Wave Conditions

The structure and dynamics of the atmospheric boundary layer is influenced by the sea waves, and that part of the atmospheric boundary layer that is directly influenced by waves is referred to as the wave boundary layer (Chalikov 1995). The wave boundary-layer height is usually related to the significant wave height. Figure 7 shows the hourly measurements of the significant wave height (\(H_\mathrm{sig}\)) and the mean wave period measured by the \(\text{ WaMoS }^\circledR \) (the location of \(\text{ WaMoS }^\circledR \) is shown on Fig. 1). The \(H_\mathrm{sig}\) ranges from 1.0 to 2.9 m, with an average value of 1.8 m and the maximum value reaching 3.6 m; the mean wave period ranges from 6.4 to 7.7 s and wind-driven waves (e.g., Holthuijsen 2007) dominate.

Fig. 7
figure 7

Time series of a significant wave height and b mean wave period

Based on the Chalikov’s 1995 parametrization relation (\(H_\mathrm{wbl} =3.7 H_\mathrm{sig}\), where \(H_\mathrm{wbl}\) is the wave boundary-layer height), the height of the wave boundary layer (WBL) rarely exceeds 13 m. However, it should be noted that the height of the WBL is also strongly dependent on the wave age, which is not considered here. As a result, the measurements of temperature and humidity at 13.4 m and above are not significantly affected by the wave conditions. The correlation coefficient between the significant wave height and the temperature/humidity time series at 13.4 m is 0.07/0.16, again implying the influence of sea-surface waves on temperature and moisture measurements is insignificant. In addition, studies have shown that the waves primarily influence the momentum flux (Chalikov and Rainchik 2011), while the influence on the scalar fluxes is controversial (Veron et al. 2011). In the following analyses, the influence of the sea-surface waves is not considered.

4.3 Flux–Gradient Relationships for Temperature (\(\phi _{T})\) and Humidity (\(\phi _{q})\)

Figure 8 shows the calculated \(\phi _{T}\) from the measurements and several parametrizations for \(\phi _{T}\) as a function of the stability parameter. The fitted AT2005 function based on the medians in 11 bins spaced approximately evenly on a logarithmic scale is

$$\begin{aligned} \phi _T =\frac{\left( {1-14.9\zeta } \right) ^{-1/2}+\zeta ^{2}\left( {1-319.6\zeta } \right) ^{-1/3}}{1+\zeta ^{2}}. \end{aligned}$$
(21)

As can be seen from Fig. 8, although the scatter is large under less unstable conditions, the bin-averaged values suggest that the two scaling laws in the AT2005 scheme are robust. For example, when \(\zeta < -2\), the “\(-\)1/3” law follows the observations fairly well; when \(\zeta > -2\), the “\(-\)1/2” scaling law also agrees with the measurements. The MRF and P2009 schemes that only include the “\(-\)1/2” scaling law do not capture the variations in \(\phi _{T}\) for the whole range of instability, which implies that the flux–gradient relationship for temperature does display different scaling laws under the two regimes (i.e., the slightly to moderately unstable regime and the convective regime). Recent phenomenological studies suggest that the anisotropy of turbulent eddies (Katul et al. 2011) and the interactions between turbulent eddies and the vertical temperature profile (Li et al. 2012b) are the causes of this “\(-\)1/2” scaling law under slightly to moderately unstable conditions. The “\(-\)1/3” scaling law is primarily due to the fact that cancels out in Eq. 1 and is also rooted in Kolmogorov’s “\(-\)1/3” scaling law in the inertial subrange (Katul et al. 2011; Li et al. 2012b). In order to parametrize the flux–gradient relationship of temperature over the whole range of instability, the two different scaling laws must be taken into consideration as in the COARE algorithm and the AT2005 scheme.

Fig. 8
figure 8

The flux–gradient relationship for temperature, \(\phi _{T}\), as a function of stability variable \(\zeta \). The error bar indicates the standard deviation within each bin

The fitted function Eq. 21 based on the AT2005 scheme is significantly lower in values than the COARE algorithm, particularly under convective conditions. The differences between the proposed function and the MRF and P2009 schemes are more or less within the uncertainty range. The scatter in the raw data and the differences between the proposed function and the other schemes imply that the temperature profiles are not well resolved. Similar features are also observed in Edson et al. (2004). However, the two different methods for computing the temperature gradients yield almost identical results, which provides some confidence in the parametrization (Eq. 21) proposed in our study.

Following Salesky and Chamecki (2012), the impact of random errors in turbulence measurements on \(\phi _{T}\) is assessed in our study. For each 30-min segment, the errors in \(\phi _{T}\) are calculated based on Eq. 20. As can be seen from Fig. 9, the relative errors in and the vertical gradient of potential temperature increase significantly with instability, which are primarily caused by the decreasing values of and the differences between the potential temperature at different heights, respectively. The errors in the sensible heat flux remain almost at the same level for all instability conditions. These features are in qualitative agreement with Salesky and Chamecki (2012). Figure 10 shows the error bars on \(\phi _{T}\) (derived from method 1) induced by the random errors in turbulence measurements. It is clear that the random errors in turbulence measurements do cause large errors in \(\phi _{T}\), particularly under close-to-neutral conditions. It is also clear that when the errors in \(\phi _{T}\) are considered, the predictions from the COARE 3.0 scheme still deviate from the measurements under unstable conditions. Among the 194 segments, there are 141, 47 and 6 segments in the stability range of \(\zeta \ge -2\), \(-2 > \zeta \ge -10\) and \(\zeta <-10\), respectively. In these three ranges, the \(\phi _{T}\) values predicted by the AT2005 scheme that are within the error bars of \(\phi _{T}\) are 113 (80 %), 45 (96 %) and 6 (100 %), respectively. The good performance of the AT2005 scheme arises from the fact that it is fitted using the same dataset. For the case of the COARE 3.0 scheme, the \(\phi _{T}\) values predicted that are within the error bars of \(\phi _{T}\) are 95 (67 %), 10 (21 %) and 2 (33 %), respectively. For the MRF scheme, the three values are 105 (74 %), 35 (74 %) and 6 (100 %); for the case of P2009, the three values are 96 (68 %), 28 (60 %) and 6 (100 %). These values are also listed in Table 2. As a result, when the random errors of turbulence measurements are taken into account, it can be seen that more than 65 % of the measurements of \(\phi _{T}\) are still in agreement with predictions by the MRF scheme and the P2009 scheme for the three instability ranges. This implies that the MRF scheme and the P2009 scheme still perform reasonably well for our site. This result is different from that given by Salesky and Chamecki (2012) who observed that the vertical error bars on \(\phi _{T}\) are outside of the empirical curves for many points under unstable conditions and concluded that additional non-dimensional parameters beyond the Monin–Obukhov similarity framework need to be included. In our study, the different schemes seem to capture a majority of data points when the random errors are considered, indicating that the dominant transport mechanism is adequately described by the Monin–Obukhov similarity theory.

Fig. 9
figure 9

Bin-averaged errors of the friction velocity, sensible heat flux and the gradient of potential temperature

Fig. 10
figure 10

The errors of \(\phi _{T}\) induced by random errors in turbulence measurements. The error bar indicates the standard deviation for each segment

Table 2 Number of segments that agree with predictions after considering the random errors

Similar to Fig. 8, Fig. 11 shows the flux–gradient relationship for specific humidity; the fitted function for \(\phi _{q}\) based on the AT2005 scheme is as follows:

$$\begin{aligned} \phi _q =\frac{\left( {1-16\zeta } \right) ^{-1/2}+\zeta ^{2}\left( {1-40.1\zeta } \right) ^{-1/3}}{1+\zeta ^{2}}. \end{aligned}$$
(22)

The fitted \(\phi _{q}\) in this study is very close to that used by the COARE 3.0 algorithm. In addition, the agreement between the measurements and the MRF and P2009 parametrizations is better for \(\phi _{q}\) than for \(\phi _{T}\), although the scatter in \(\phi _{q}\) is larger than that in \(\phi _{T}\). This is consistent with previous studies that also evaluated the flux–gradient relationships for temperature and specific humidity over open ocean (Edson et al. 2004).

Fig. 11
figure 11

The flux–gradient relationship for specific humidity, \(\phi _{q_B}\) The error bar indicates the standard deviation within each bin

The “\(-\)1/2” and “\(-\)1/3” scaling laws are not clearly distinguished by examining the raw data. Nonetheless, by examining the bin-averaged values, it is interesting to notice that the methods of calculating vertical gradients alter the results when the data points are sparse (i.e., when \(\zeta < -2\)). The method 1 seems to suggest a “\(-\)1/2” law and the method 2 favours the “\(-\)1/3” scaling law under convective conditions. These features are cautiously noted here and we expect that more experimental datasets will be needed to resolve these differences. The non-negligible sensitivity of humidity gradients to the fitting methods is also observed in P2009.

In addition, the errors of \(\phi _{q}\) derived from method 1 are plotted in Fig. 12, with comparison to the four empirical parametrizations. Among the selected 234 segments, there are 165, 56 and 13 segments in the stability range of \(\zeta \ge -2, -2 > \zeta \ge -10\) and \(\zeta <-10\), respectively, as shown in Table 2. In these three ranges, the values of \(\phi _{q}\) predictions by the AT2005 scheme that are within the error bars of \(\phi _{q}\) are 76 (46 %), 27 (48 %) and 8 (62 %), respectively; the values of \(\phi _{q}\) predictions by the COARE 3.0 scheme that are within the error bars of \(\phi _{q}\) are 77 (47 %), 26 (46 %) and 8 (62 %), respectively. As such, the COARE 3.0 scheme performs as well as the AT2005 scheme for \(\phi _{q}\) at this site. For the MRF scheme, the three values are 83 (50 %), 26 (46 %) and 9 (69 %); for the case of P2009, the three values are 51 (31 %), 29 (52 %) and 9 (69 %). These relatively low ratios for \(\phi _{q}\) as compared to those for \(\phi _{T}\) indicate that many measurements are truly outside of the predictions of the four schemes, even when the random errors in the turbulence measurements are taken into account. This is in agreement with Fig. 11, which displays larger variability in \(\phi _{q}\) as compared to \(\phi _{T}\) in Fig. 8.

Fig. 12
figure 12

The errors of \(\phi _{q}\) induced by random errors in turbulence measurements. The error bar indicates the standard deviation for each segment

4.4 Comparison Between \(\phi _{T}\) and \(\phi _{q}\)

The ratios of \(\phi -_{q}\) to \(\phi _{T}\) are shown in Fig. 13 in order to examine the similarity between the flux–gradient relationships for temperature and specific humidity. Most of the ratios are larger than unity under unstable conditions, which is in agreement with the study of P2009. However, unlike P2009 who suggest that the ratio of \(\phi _{q}\) to \(\phi _{T}\) is roughly a constant (= 1.2) as can be seen from Eqs. 18 and 19, the ratio of the fitted two functions (Eqs. 21 and 22) in our study increases with instability. Changing the method of calculating the vertical gradient does not alter this trend qualitatively, particularly under slightly/moderately unstable conditions. The ratio of \(\phi _{q}\) to \(\phi _{T}\) ranges from 1 under slightly unstable conditions to 2 under convective conditions. From Eq. 23 it can be seen that the ratio of the turbulent diffusivities of temperature and water vapour, which is equal to the ratio of \(\phi _{q}\) to \(\phi _{T}\), is larger than unity and increases with instability. The finding that heat is transported more efficiently than water vapour is consistent with many previous studies (conducted both over land and over water surfaces), which suggested that it is the active role of temperature that causes a larger transport efficiency as compared to water vapour (Warhaft 1976; Katul and Parlange 1994; Katul and Hsieh 1999; Li et al. 2012a) As instability increases, the buoyancy production of the sensible heat flux (\(g \overline{T^{\prime } T^{\prime }}/\overline{T}\)) becomes more efficient than that of the latent heat flux (\(g \overline{q^{\prime } T^{\prime }}/\overline{T}\)). In particular, the active role of temperature is theoretically studied within the context of flux–variance relationships (Katul and Hsieh 1999) and flux-structure parameter relationships (Li et al. 2012a),

$$\begin{aligned} \frac{K_T}{K_q}=\frac{\overline{{w}^{\prime }{\theta }^{\prime }}/\left( {{\text{ d }\overline{\theta }}/{\text{ d }z}} \right) }{\overline{{w}^{\prime }{q}^{\prime }}/\left( {{\text{ d }\overline{q}}/{\text{ d }z}} \right) }=\frac{\phi _q}{\phi _T}. \end{aligned}$$
(23)

Another measure of the dissimilarity, the turbulent transport efficiency (see e.g., Li and Bou-Zeid 2011), is also quantified. The ratio of the turbulent transport efficiencies of heat and water vapour is defined as

$$\begin{aligned} \frac{R_wT}{R_wq}=\frac{\overline{{w}^{\prime }{\theta }^{\prime }}/\sigma _w \sigma _\theta }{\overline{{w}^{\prime }{q}^{\prime }}/\sigma _w \sigma _q}. \end{aligned}$$
(24)
Fig. 13
figure 13

The ratios of \(\phi _{q}\) to \(\phi _{T}\). The black line is the ratio of the two proposed relations (Eqs. 21 and  22)

Figure 14 shows the ratios of \(R_{wT}\) to \(R_{wq}\) as a function of the stability parameter \(\zeta \), where it can be seen that the bin-averaged values of \(R_{wT}/R_{wq}\) show a slight increasing trend with instability (\(-\zeta \)). We only report the bin-averaged values within the range of \(\zeta \) from \(-\)20 to \(-\)0.01 since the number of data points in each bin is limited at high instabilities. This result underlines the previous finding that heat is transported more efficiently than water vapour.

Fig. 14
figure 14

The ratios of \(R_{wT}\) to \(R_{wq}\) (Eq. 24). The large black dots are bin-averaged values

5 Summary and Discussion

Using simultaneous measurements of eddy-covariance flux data and mean vertical profiles of temperature and specific humidity, the flux–gradient relationships for temperature and specific humidity are evaluated over a wide range of instability (\(0.1< -\zeta <50\)). Both \(\phi _{T}\) and \(\phi _{q}\) show departures from the predictions of Kansas-type functions and follow the free convection “\(-\)1/3” scaling as \(- \zeta \rightarrow \infty \). The AT2005 scheme captures the variations in the measured dimensionless temperature and humidity gradients well. For \(\phi _{q}\), the agreement between the COARE 3.0 algorithm predictions and the observations is also good. It should be noted that small flux values have been removed in this study, which might result in an underrepresentation of the near-neutral conditions.

The errors in the flux–gradient relationships caused by the random errors in the turbulence measurements are assessed in this study. The errors in the flux–gradient relationships mainly arise from errors in and the vertical gradient of temperature/specific humidity, both of which increase with instability. The error bars on \(\phi _{q}\) are larger than \(\phi _{T}\) and thus more data points do not agree with predictions for \(\phi _{q}\) than for \(\phi _{T}\), especially under close-to-neutral conditions. When the error bars on \(\phi _{T}\) are considered, the majority of data points agree with the AT2005 scheme, the MRF scheme and the P2009 scheme, implying that the dominant transport mechanism is adequately described by the Monin–Obukhov similarity theory at this site and the random errors do explain some of the differences among various schemes.

The results also indicate that the ratios of \(\phi _{q}\) to \(\phi _{T}\) range from 1 to 2 for \(0.1< -\zeta <10\). As instability (\(-\zeta \)) increases, the ratio approaches a constant value of 2, indicating that the turbulent transport of heat is more efficient than that of water vapour under unstable conditions. The ratio of turbulent transport efficiencies of heat and water vapour (\(R_{wT}/R_{wq}\)) underscores this finding.