Non‐turbulent motion identified from properties of transport and its influence on the calculation of turbulent flux

Turbulent fluxes play a critical role in atmospheric science and are usually calculated from the eddy covariance system. However, the presence of motions of larger scale and the bias from observational instruments often affects the procurement of turbulent fluxes. Based on the Hilbert–Huang transform, the properties of transport can be defined and used to distinguish non‐turbulent motions from the observational data, and a new way of decomposing the variables is thereby put forward to reconstruct the turbulent flux series. To quantify the influence of non‐turbulent motions on the calculation of turbulent flux, the non‐turbulent motions extracted from the five‐level turbulence data of the Tianjin 255‐m meteorological tower from July 1 to August 31, 2017 are examined. The results reveal that the presence of non‐turbulent motion can lead to a universal overestimation of turbulent flux, and the degree of overestimation increases with the complexity and the intensity of non‐turbulent motion. According to the consequences above, an empirical relationship, together with the corresponding coefficients, is given to provide a guidance on the correction of turbulent fluxes in practical use, such as the simulation of atmospheric turbulence and the parameterization in meteorological and climate models.


INTRODUCTION
Turbulent transport is the most important way of exchanging heat, moisture, and momentum between the Earth's surface and the atmosphere (Caughey et al., 1979;Garratt, 1994;Stull, 1988), and is achieved through the action of turbulent motions (Chan et al., 2021;Marusic et al., 2010;Robinson, 1991).The measurement of turbulent flux is important for a wide range of applications, and the most universal method is the eddy covariance (EC) system.
The EC system is based on the principle of measuring the fluctuations of vertical wind velocity and a second quantity of interest, and using these measurements to estimate the turbulent flux (Aubinet et al., 2012;Baldocchi et al., 1988;Burba, 2013;Verma, 1990).Although the EC system is widely used in studies over past decades (Aubinet et al., 1999;Baldocchi, 2003;Crawford et al., 2017;Webb et al., 1980;Wilson et al., 2002), there are some limitations that can affect the accuracy and reliability of the measurements, such as the representativeness, sensitivity, and the standardization in installation and quality control (Berger et al., 2001;Moore, 1986;Vickers & Mahrt, 2003).Regardless of the technical problems, the signal detected by the EC system crosses a wide range of frequencies, and motions with larger scales can also be detected.These motions can vary from synoptic processes to ones with timescales that are slightly greater than turbulence (Mahrt, 2007(Mahrt, , 2010b)), the latter of which are considered as submesoscale motions (Anfossi et al., 2005;Conangla et al., 2008;Mahrt, 2008;Sun et al., 2004).As the traditional time-averaging method fails to distinguish large-scale motions from turbulent motions, the calculated turbulent flux is often contaminated, making the applications of similarity relationships and the simulations become difficult or even invalid (Li et al., 2016;Ren et al., 2019), especially in stable stratifications, where submesoscale motions make a great difference as intermittency happens frequently (Mahrt, 1998(Mahrt, , 2010a(Mahrt, , 2014)).
To lessen the influence of motions of larger scale, the calculation of turbulent flux from the EC system usually goes through a detrending or high-pass filtering process on the raw time series, and the signal with larger timescales is referred to as "low-frequency bias" (Aubinet et al., 1999).The method of this filtering process includes linear detrending, block averaging, autoregressive filtering, and other detrending schemes (such as quadratic detrending), and the benefits and defects of these methods have been discussed in detail (Culf, 2000;Moncrieff et al., 2005;Rannik & Vesala, 1999).However, these methods are not well-defined in either mathematics or physics (Kaimal & Finnigan, 1994).On the one hand, the arbitrary usage of detrending may enhance the effect of low-frequency loss (Aubinet et al., 2012), which merely makes the systematic error shift to inexplicable fluctuations, rather than disappear.On the other hand, the signal extracted from the original data through detrending is by no means guaranteed to be the submesoscale motions or motions of larger scales.Therefore, whether or not to apply detrending depends on the researcher's aim, to a large extent.
In this situation, how to effectively distinguish, separate, and eliminate the motions of larger scales detected by the EC system has become a major problem in the calculation of turbulent flux.Based on studies into the turbulent intermittency, Vickers and Mahrt (2003) divided the turbulent and mesoscale flux by a cospectral gap; Muschinski et al. (2004) found a boundary between large-and small-scale motions, which is identified by a "plateau" indicating a maximum point or a range of large values in the spectrum.These efforts can both be of guidance to distinguishing the large-scale motions, yet to effectively separate the motions in the calculation of turbulent flux, a practical method including the algorithm is needed.Huang et al. (1998) put forward a methodology of analyzing nonlinear and non-stationary time series based on the Hilbert spectrum, which was later called the method of Hilbert-Huang transform.The arbitrary-order Hilbert spectral analysis derived from this method has been proven effective in studies concerning turbulence (Huang et al., 2008(Huang et al., , 2011;;Wei et al., 2016).Based on the method of Hilbert-Huang transform, Wei et al. (2017) separated fine-scale turbulence from large-scale motion by a spectral gap, and the original data can be reconstructed as pure turbulence time series.Ren et al. (2019) developed an automated algorithm to objectively find the frequency location of the spectral gap, and analyzed the effects of submesoscale motions on the turbulent transport in stable boundary layers.Urbancic et al. (2021) developed a general theory of finding submesoscale motions based on the negative minimum of the Eulerian autocorrelation function.Li et al. (2023) evaluated the effectiveness of the wavelet analysis method for turbulent flux calculation under non-stationary conditions, which suggests a promising future in quantifying turbulent fluxes in non-stationary conditions.Recent studies tend to focus on stable stratifications, as the similarity relationships are shown to be violated, and the turbulent intermittence also makes a huge difference in investigations on turbulent transport and the identification of submesoscale motions.However, non-turbulent motions are also ubiquitous under unstable stratifications (Lan et al., 2022;Ren et al., 2023).Represented by local circulations, large roll vortices, and entrainment-lead motions, non-turbulent motions sometimes share the same scale with atmospheric turbulence, and are thus difficult to distinguish under traditional methods of spectroscopy.Nevertheless, as essentially different kinds of motions, there still exists potential for identifying non-turbulent motions by the difference of their dynamics (Borrell & Jiménez, 2016;Chauhan et al., 2014;da Silva & dos Reis, 2011), and the emphasis for the investigation of turbulent flux will be on the properties of transport.
In this study, a new kind of algorithm is put forward to identify non-turbulent motions from the observational data, based on the properties of transport, and to effectively calculate the turbulent fluxes without them being disturbed by non-turbulent motions, a new criterion of decomposing time series with fewer assumptions is also introduced to substitute the Reynolds decomposition.The five-level turbulence data from the Tianjin 255-m meteorological tower are examined to verify the newly introduced algorithm and quantify the influence of non-turbulent motions on the calculated turbulent flux.In addition, the detected and separated non-turbulent motions are classified according to different criteria, to reveal how their properties (such as complexity, intensity, and monotonicity) affect the calculation of turbulent flux.Finally, a referential criterion is given to help researchers diagnose to what degree the data collected from the EC system are "contaminated" by non-turbulent motions, and how can one lessen their influence and acquire fluxes with pure turbulent transport.

Observational site
To test the consistency of the turbulent data and thereby justify the representativeness of the calculated turbulent flux, the data from a multiple-level turbulence observational station are required.In this study, the observational data from the Tianjin 255-m meteorological tower (39.08  (Wei et al., 2018;Ye et al., 2015), as the surrounding underlying surface consists of mid-rise buildings no more than 25 m high, grassland, roads, and vegetation.The meteorological tower can provide overall turbulence information from the surface layer to part of the mixed layer, as the height of 255 m is situated at the lower part of the mixed layer in Tianjin (Fernando, 2010;Guo et al., 2016;Li, 2019;Oke, 1988;Roth, 2000).
The sonic anemometers are supported by 5-m horizontal booms on the tower oriented to the East, and the sampling frequency is 10 Hz.They provide wind components and air temperature data, and there is also a net radiometer (CNR4, Kipp & Zonen, Delft, Netherlands) installed at 40 m with a sampling period of 1 min, to help verify the energy enclosure of turbulent flux.The turbulence observational data of the meteorological tower are from July 1 to August 31, 2017 (a period typically dominated by unstable stratifications), for all five sonic anemometers except the 200 m one, as thois one was under maintenance from July 1 to July 31.

Data pre-processing
To effectively detect the existence of non-turbulent motions, the size of data samples has to be properly chosen.Based on pre-experiments and prior experiences, the non-turbulent motions are weak with a short time-space scale in stable stratifications, and the usual choice of 30-min data samples is sufficient to distinguish them.However, in unstable conditions, the non-turbulent motions can be intense with a maximum timespan longer than 30 min, which means a short averaging time may not be able to completely reflect the properties of non-turbulent motions.Therefore, the averaging time is set to be one hour as a balanced choice for distinguishing ability and the volume of data, and if the proportion of spikes in a given 1-hr sample is more than 20%, it would be discarded.The raw turbulence observational data are pre-processed in two different ways: 1. Use the EddyPro software (Advanced 7.0.9,LI-COR Biosciences) to process the data, including error flags, despiking (Vickers and Mahrt, 1997), double coordinate rotations, Reynolds decomposition, and detrending.This is the universal method of processing the original turbulence data, which is used by most of the related research.The calculated turbulent flux based on this way of pre-processing will be marked as "original" flux in this study.2. Remove outliers using EddyPro software and subtract the mean values of each variable, but manually apply horizontal coordinate rotation.The original w-wind is kept, and there is no Reynolds decomposition or detrending process.Substitutive algorithms will be applied in subsequent processes.
In the first way of pre-processing, the second coordinate rotation (w-u rotation) is carried out to eliminate the influence of terrain, instrument inclination, and the swing of sonic anemometers due to the strong wind.However, this kind of coordinate rotation can also blur the existence of vertical turbulence motion and the possible vertical airflow.In this study, rigorous tests have been performed on the instruments, the measured w-wind is thus ensured to be the exact vertical wind speed, and most of the wind speed is no more than 5 m⋅s −1 .As this study casts enormous concern on the w-wind, w-u coordinate rotation is not applied in the second way of pre-processing.After applying horizontal rotation once, the values of u, v, and non-zero w-wind can be directly extracted from the original data.

Method of Hilbert-Huang transformation
Hilbert-Huang transformation is considered capable of handling nonlinear and non-stationary signals with multiple components (Huang et al., 1996(Huang et al., , 1998(Huang et al., , 2008)), such as atmosphere turbulence data, and the basic step of this method are: 1. Conducting Empirical Mode Decomposition (EMD): EMD is an algorithm to decompose a real signal X(t) into a sum of intrinsic mode functions C i (t) and a residual function r n (t).For each intrinsic mode function (IMF), it satisfies two criteria: (a) the difference between the number of local extrema and the number of zero-crossings must be zero or one; (b) the running mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.
2. Performing the Hilbert transform for each IMF according to the definition: where S i (t) is the analytical signal of the ith IMF, is the instantaneous amplitude, and ] is the instantaneous phase.The instantaneous frequency can be defined as The original signal can then be written as: 3. Calculating the Hilbert marginal spectrum: With the time series of  i (t) and A i (t), the Hilbert spectrum, H(, t), can be defined as a frequency-time distribution of the amplitude, and the Hilbert marginal spectrum h() can be defined as: From the time series of  i (t) and A i (t) for each IMF, the joint probability density functions (PDF) p(, A) can be calculated.According to joint PDF, the arbitrary-order Hilbert marginal spectrum can be defined in an amplitude-frequency space by a marginal integration of the joint PDF as where  and A are the instantaneous frequency and instantaneous amplitude at time t, respectively.
4. Determining IMF's contribution to flux: As the flux is dependent with two elements (e.g., momentum flux depends on u-and w-wind), the contribution of IMFs follows a two-dimensional grid distribution of the products of IMFs of each element.In this way, the contribution can be directly defined as the flux generated by certain pairs of IMFs normalized by the overall flux: where F ij is the flux calculated by the ith and jth IMFs of the two elements of the flux; for example, if the sensible-heat flux is considered, then F 45 refers to the flux that is calculated by the fourth IMF of w-wind and the fifth IMF of potential temperature.Based on the distribution of contribution of IMFs, the main source of flux can be understood, which can be used to quantify the influence of frequency ranges for turbulent transport and identify the dominant ones.In addition, the influence of non-turbulent motions can also be judged.

Different types of turbulent transport
Theoretically, the turbulent transport reflects the gradient along the vertical direction, and the time-averaged turbulent flux can describe the average gradient.As atmospheric turbulence varies greatly in scale, the generation of turbulent transport can also be different, depending on the stratification state.From the time series of covariances and the amplitude of analytical turbulent signals, two main types can be intuitively defined: the steady-motion-dominant transport, and the burst-motion-dominant transport.For the former type, the instant flux always fluctuates around its mean value, and for the latter type, the instant flux is lower than the mean value for most of the time, but can increase to a quite high value within a short time.
For a mathematical definition of these two types, the probability distribution P k (A) of the instant amplitude A for each analytical IMF signal (where k denotes the specific element of the flux), together with the relative contributions of each IMF, and Q k (A) are similar to each other, the most frequently appearing amplitude contributes most to the turbulent flux, which means the turbulent motion is steady, and the flux is steady-motion-dominant.If there is a significant difference between P k (A) and Q k (A), it means a recently occurring amplitude can make a huge contribution to the turbulent flux, which reveals a burst-motion-dominant case.Note that as Q k (A) contains the information from two elements of the flux, a "steady" event usually means the motions of both two elements make steady contributions to the flux, while a "burst" event means either one of the elements can make a huge contribution within a short time.
To quantify the difference between P(A) and Q(A), the relative entropy is practical.The Kullback-Leibler Divergence (KLD, Kullback & Leibler, 1951) is a widely accepted expression of relative entropy in probability theory and information theory, and its definition is: where p and q are the probability density and the relative contribution of a certain amplitude a i .As the value of KLD is asymmetric for p and q, it is often used as an intermediate quantity for calculating the Jensen-Shannon Divergence (JSD, cf.Endres and Schindelin, 2003;Nielsen, 2019): which is symmetric with a boundary of [0,1].If D JS → 0, the distribution of P and Q can be considered the same; if D JS → 1, P and Q are totally different.With the value of JSD, the two types of turbulent transport can be effectively distinguished.Given the fact that the non-turbulent motion, or the usually mentioned submesoscale motion, usually corresponds to abnormal values of fluxes (Vickers & Mahrt, 2006), the value of JSD can essentially provide the foundation for identifying the non-turbulent motions by quantifying the properties of turbulent transport.

Identifying non-turbulent motions through properties of transport
In practice, the bias from the synoptic motions, the systematic error from the observational instruments, and the so-called submesoscale motions can be categorized as non-turbulent motions.Apart from several detrend schemes, systematic methods of identifying the submesoscale motions in stable stratifications by means of finding the spectral gap have been put forward (Ren et al., 2019;Wei et al., 2017), and the IMFs with frequencies lower than the spectral gap frequency are considered as submesoscale components.Although the presence of a spectral gap has been widely confirmed in stable stratifications, and proven to be useful in identifying the submesoscale motions, the hybridization of turbulent and submesoscale motions brings some trouble in practice, especially in unstable conditions (Mahrt, 2014;Wei et al., 2017Wei et al., , 2021)).Moreover, different stratification states can vary hugely in turbulent scale, making it difficult for the identification algorithms in past research to self-adapt to the corresponding scales and appropriately find the submesoscale motions.Therefore, a new algorithm is put forward here, to distinguish turbulent and non-turbulent motions, based on the properties of transport revealed from the Hilbert-Huang transform.
From the contribution grid distribution of IMFs, Q ij , the values at the margin (i.e., the covariance of IMFs with a huge difference in number) reflects the interaction between large and small eddies, and a mutation of either one's amplitude could result in a significant exception value, which can be attributed to the disturbance from non-turbulent motion in the calculation of flux.As different kinds of motion, the properties of transport tend to be distinct, and a boundary analogous to the spectral gap in the last section would be helpful in identifying the non-turbulent motion.To investigate the difference in properties of transport, the gird distribution of JSD, namely D JS,k (i, ) is calculated and analyzed.There are overall three indexes of D JS , where i and j are the numbers of the IMFs of two corresponding elements, and k denotes which element is being investigated for its non-turbulent motion.Figure 1 reveals that as the frequency of IMF increases, the value of JSD would increase, which means the flux tends to shift from steady-dominant motion to burst-dominant motion as the scale of motion decreases.This is consistent with traditional turbulence theories (e.g., K41 theory), because the breaking of large eddies corresponds with smaller eddies with more complicated motions, which makes the turbulent transport more unsteady.And if the motions are non-ideal, the energy lost and the existence of turbulent intermittence may further disturb the steadiness of turbulent transport.The right two subfigures are the schematics of two kinds of motions with different JSD, and it can be judged that a higher D JS is usually accompanied by an inconsistent pair of PDF and contribution to flux, which marks the burst-dominant motion.For different kinds of stratifications of the boundary layer, the value of D JS can vary from 0.05 to 0.75, yet its overall changing pattern with scale remains invariant.However, it often appears that for a certain IMF, the D JS is abnormally high, implying a kind of motion other than the atmosphere's turbulence.Therefore, identifying the D JS that violates the correlation above is helpful to distinguish the turbulent and non-turbulent motions.For any presence of exception value of D JS , the Mann-Whitney U-test is conducted to tell whether the corresponding IMF is statistically significantly different.The main principle of the Mann-Whitney U-test can be summarized as (taking the w-wind in momentum flux as an example): 1.For a specific IMF of w with the first index of i, collect D JS,w (i, ) for all of the second indexes as D ′ 1 , and the remaining D JS (with the first indexes other than i) are collected as D ′ 2 ; 2. Mix the data from D ′ 1 and D ′ 2 , and sort them with rank numbers; 3. Calculate the sum of rank numbers T 1 and T 2 for D ′ 1 and D ′ 2 ; , where For n 1 , n 2 > 10, z is proved to obey a normal distribution;

For the original hypothesis
96, the original hypothesis is declined (for a significance level of 95%).
When the original hypothesis is rejected, although the ith IMF of w-wind yields an exception value of JSD, it still follows the overall negative correlation.However, for an exceptional IMF with the original hypothesis accepted, one can judge that the overall trend (the decrease of D JS as frequency decreases) revealed in Figure 1 is broken from the ith IMF of w-wind, and thus the IMFs with higher numbers should not be considered as turbulent motions.This process is repeated for u-wind to identify the non-turbulent motions of the other element of momentum flux.For a better understanding of the process of identification, a flow chart is given in Figure 2 to visually demonstrate this algorithm.
In fact, the boundary IMF that separates the non-turbulent from the turbulent motions that are revealed by the properties of transport can be seen as the motion near the spectral gap, as the non-turbulent motion extracted is proven to be the same as the submesoscale motion in stable stratifications.However, in other stratification states where the spectral gap may not exist (or the algorithm fails) and the sampling duration has to be lengthened, the identification can still make a difference, as long as the difference in properties of transport exists.Another advantage of this method is that no detrend process is needed, because the elimination of residuals and IMFs can play the same role, and their physical significance is clear.
Figure 3 shows the comparison between the detrended original motion and the distinguished turbulent and non-turbulent motions.It can be judged from panel (a) that the non-turbulent motion can significantly change the waveform of the time series, and further disturb the calculation of turbulent flux.Further, the reconstructed time series is already trend-free, which means that the process of eliminating the non-turbulent motions is compatible with the process of detrending.The comparison of energy spectrA in panel (b) reveals that the energy of frequencies lower than 10 −2 Hz is greatly eliminated.This means that the traditional way of detrending fails to distinguish turbulent motions and non-turbulent motions, which also confirms the effectiveness of the method of identifying non-turbulent motions through properties of transport.and latent heat flux E, is usually acquired by the eddy covariance method and defined as follows:

Definitions of turbulent fluxes from different components of motion
where  is the potential temperature,  and C p are the air density and specific heat capacity at constant pressure; c is the concentration of carbon dioxide,  is the latent heat, and q is the concentration of water vapor, respectively, and the Reynolds decomposition defines the turbulent fluctuations.However, under the algorithm introduced in Section 3.1, the main assumptions of Reynolds decomposition, namely the average of the fluctuating variable V ′ = 0, tend to fail, and Equations (9a)-(9d) are thus no longer effective.In this situation, to effectively calculate the turbulent flux, a new way of decomposing a certain physical quantity V is put forward: where Ṽ is the non-turbulent component of variable V, and V • is the turbulent component.Note that for w and V, the original time series can be rewritten in the sum of analytical signals, and suppose that for both w and V, the boundary IMF between turbulent and non-turbulent motions is marked by numbers m w (m V ), and there are: S w,i (t)e i w,i (t) + r w (t), (12a) where n w and n V are the overall number of IMFs for w and V, and therefore the four components in Equation ( 11) can be expressed by: (13d) where Re denotes the real part of the complex number.The specific forms of w and Ṽ are not expanded to avoid confusion.
Based on this new way of decomposing the original data, the components of overall flux are reconstructed.Only the last term, w • V • , can describe the turbulent flux of V, and the first term, w Ṽ, plays a similar role as advective flux, F V , yet it also contains the effect of non-turbulent motions.The second and third terms cannot be eliminated based on stationarity assumption, and thus are included as the interaction terms for the following analysis.
Figure 4 shows the examples of four kinds of newly defined fluxes and the comparison of two kinds of turbulent fluxes.The main difference between the newly defined turbulent flux and the original turbulent flux exists in the maximum values, and there are fewer outliers in the former.The time series of four kinds of fluxes reflect significantly different orders of magnitude, which implies a great influence from the non-turbulent motions.

Classification of non-turbulent motion
In practice, the calculated turbulent flux tends to be affected by non-turbulent motions (Mahrt, 2010b;Vickers & Mahrt, 2006), and how to identify and eliminate the influence of non-turbulent motions is critical in quantifying the turbulent transport for the further utility of climate and synoptic models.As mentioned above, the frequency of non-turbulent motion varies from 0 to the frequency of the spectral gap, which can significantly affect the calculation of fluxes, as Equations ( 13a)-(13d) show.Theoretically, these equations are sufficient to describe the influence of non-turbulent motion, but in practice, it is necessary to divide non-turbulent motions into several categories, in order to intuitively judge the over-/underestimate of turbulent flux.Also, the classification of non-turbulent motions will depend greatly on the new way of decomposing and the identification of non-turbulent motion.
For a certain IMF, its amplitude S V, can vary smoothly with time, but the frequency is limited around a certain value (namely the mean frequency of that IMF).Therefore, an assumption of constant frequency for each IMF is rational.Paying attention to the specific characteristics of exponential integration, one can integrate Equations ( 13a)-(13d) by parts: where an assumption is made that dt , and there is also  w,i ≡ d w,i (t) dt .Note that the first term with the time span T in the denominator can be ignored, and thus Equations ( 13a)-(13d) can be rewritten as . It can be judged that the expressions in Equations ( 14b)-(14d) can further be integrated by parts based on the constant-frequency assumption, yet it might cause the accumulation of system error, and thus is not applied.
Although non-turbulent motion plays no role in the calculation of turbulent flux w • V • , it does affect the calculation of w ′ V ′ .This means the difference between w • V • and w ′ V ′ is controlled in Equations ( 14a)-( 14c).The systematic factors from non-turbulent motions are the time derivatives of time series, namely their monotonicity.Under this circumstance, the zero points of the first derivative of every non-turbulent time series are identified, and the time series sample is cut into several segments by these zero points, making every segment either monotonically increasing or decreasing.Based on the technical processes above, the length of time span, the first-order time derivative of non-turbulent time series, and the number of monotonic segments of non-turbulent motions are of concern.Thus, the classification of non-turbulent motion can be defined in different criteria as follows: 1.By the monotonicity of two flux elements (e.g., w and T when calculating sensible-heat flux) in a single segment of time series, and there are four cases: "+, +", "−, +", "+, −", and "−, −".Note that the zero points of first derivative from two flux elements are put together to cut the time series, in order to ensure that for each segment of a time series, the two elements are both monotonic.From this kind of definition, the effects of monotonicity of non-turbulent motion and the interaction of two elements can be distinguished.2. By the complexity of the non-turbulent motion.Technically, the ratio of the total number of monotonic segments within the data sample to the length of the data sample can be collected and defined as the complexity property: where N i, is the number of monotonic segments of two elements of turbulent flux, and L is the length of the data sample.Both N i, and R i, serve as the index of complexity, and the latter is independent of the length of data.In this study, the number of N i, from 1 to 300 are identified as different cases, and also a set of two cases is identified: no less than 100 time series segments, and its contrary.If there are differences between relatively stable and fluctuating non-turbulent motions, they will be revealed under this kind of classification.3.By the intensity of the non-turbulent motion.The average amplitudes of non-turbulent motions of two elements within a data segment or the whole sample can be similarly calculated as the square root of their normalized sum: where A i and A  are two elements of turbulent flux, and they are normalized by their maximum values within the overall time span (e.g., two months in this study).This classification criterion can indicate to which degree the intensity of non-turbulent motions can affect the calculation of turbulent flux.
In addition, the original fluxes (those calculated based on Reynolds average) are also divided into two cases according to their signs.In this way, the changes in numerical value and in absolute value can be distinguished.The examples of calculating R u,w (Example Script 1) and I u,w (Example Script 2) are provided, based on the non-turbulent motion of u-and w-wind, ũ and w, extracted from the method in Section 3.1:

Comparison of flux values
The effects of reconstruction of turbulent flux can be quantified by the absolute/relative change in the reconstructed turbulent flux compared to the original turbulent flux, for a certain sample of data.Likewise, given the classification criteria mentioned above, whether there is a systematic over-/underestimate of turbulent flux can also be tested for different segments of time series and different cases.
For classification criterion (1), different segments of time series may be classified into the same group, and they are thereby comparable.In this study, the absolute changes of the reconstructed turbulent flux compared to Example Script 1. Calculating R u,w.

Example Script 2. Calculating I u,w
1. Calculate the instantaneous amplitude of ũ and w by Equations ( 1) and ( 2): A ũ′ = hilbert ( ũ); A w′ = hilbert ( w); 2. Normalize A ũ′ and A w′ by their maxima: 2 ) ; 4. I u,w can be cut into data segments according to the monotonicity-changing points in Example scripts 1, if needed.
the original turbulent flux are calculated: and the mean value and standard deviation of ΔF V are given for each case.The mean value and standard deviation are calculated according to the weight of length of time series segments, to avoid the influence of differences in data length.The relative changes are not considered for the possible abnormal values and near-zero original fluxes.A box plot is also given to reveal the distribution of over-/underestimation values within a certain case defined in Section 4.1.
For classification criterion (2), the absolute change of turbulent flux for every data sample, ΔF V | sample , is calculated, and all samples are sorted out into certain case groups based on the complexity indexes N i, or R i, .The mean value and 95% confidence intervals are also calculated to give a better understanding of the degree and level of significance of the influence from non-turbulent motions with different complexity.
For classification criterion (3), the relationship between ΔF V | sample and I non−turb will be shown in a diagram, which will be of guidance in the correction of turbulent flux under non-turbulent motions of different intensity.
In a word, from these three kinds of classification criteria, the reconstructed turbulent flux can be assessed by the value of the over-/underestimate, and the influence of the monotonicity, complexity and intensity of non-turbulent motions can be evaluated.

RESULTS AND DISCUSSIONS
The last sections discussed the methods of performing Hilbert-Huang transformation, identifying the non-turbulent motions from the properties of transport, and the reconstruction of turbulent flux.In this section, the comparison of reconstructed turbulent fluxes and original fluxes (which are calculated by the first method of pre-processing in Section 2.2) will be analyzed to find how the non-turbulent motions affect the calculation of fluxes, and whether there are some critical parameters that can be used to conveniently correct the turbulent flux for practical use.

The correction of turbulent flux within a monotonic segment
Figure 5 shows the correction of turbulent momentum flux under different cases from criterion (1), distinguished by different colors of box plots, and from criteria (2) and (3) distinguished by groups on the x-axis.Here the turbulent flux is composed of the prevailing and crossing wind, yet the cases are divided only by u-and w-wind, as the prevailing wind contributes most to the turbulent flux.As the values of momentum flux are negative in general, the overall pattern is clear that in most cases the flux tends to be overestimated, compared to the original flux that is calculated under the usual method (the first way of pre-processing in Section 2.2, to be exact), while there are also some cases of underestimation.Usually, the contribution of non-turbulent motions to the momentum flux tends to be in the same direction as the contribution of turbulent From the comparison of different cases in each group of the box plot (adjacent different colors, marked in the legends), it can be judged that different combinations of monotonicity of non-turbulent motion are irrelevant to the value of flux correction.Moreover, a comparison of the different cases among the four groups of each box plot (marked on the x-axis) reveals that non-turbulent motions of larger average amplitude and a larger number of monotonic segments can result in a higher overestimation of momentum flux.This indicates that the intensity and complexity of non-turbulent motions can be relevant for the correction of turbulent flux.
Figure 6 shows the correction of turbulent sensible-heat fluxes under different cases.The overall pattern is the same as for the momentum flux, namely no matter whether the original flux is positive or negative, the absolute value of flux tends to be overestimated under traditional methods; the correction also increases from 40 to 160 m, and decreases at 200 m; the non-turbulent motions with higher complexity also affect the calculation of flux to a higher degree (the influence of intensity is also similar, and thus omitted).However, when the non-turbulent component of potential temperature is monotonically increasing (denoted by red and blue boxes), the degree of overestimation tends to be higher.Different from the u-wind, an increasing or decreasing trend often exists in non-turbulent motions, especially in near-neutral conditions, and thus separating them into different cases is necessary.Moreover, the overestimation is much more serious under unstable stratifications (denoted by first and second groups), which reveals that the reconstruction of turbulent flux in unstable boundary layers is as important as in stable ones.
Figure 7 shows correction of turbulent CO 2 and water vapor flux under different cases at 80 m (as the measurements at 40 m are lacking).The overall pattern is also similar, but the influence from non-turbulent motions is only sensitive when N mono < 100, which means an increase of the complexity of non-turbulent motions might be able to blur the influence of concrete forms of motion, yet in practical use, the difference among these cases can still be neglected.
From the comparison of the correction of turbulent fluxes under different cases distinguished by criterion (1), a systematic image can be acquired: the contribution of non-turbulent motions to the momentum flux tends to be in the same direction as that of the turbulent motions,

F I G U R E 7
The correction of turbulent latent heat/CO 2 flux under different cases.The different styles, the groups on the x-axis, and the components of the box plots, please refer to Figure 5. and once the former is not eliminated in the calculation of turbulent flux, the result tends to be overestimated in absolute value.However, there are also some opposite circumstances leading to underestimation.Complex kinds of non-turbulent motion tend to overestimate mke the turbulent flux to a higher degree, and the overall influence increases with height within the near-surface layer.

The correction of turbulent flux within a data sample
The last section can provide a straightforward image on how the different combination of monotonic segments, the complexity, and the intensity of non-turbulent motions can affect the calculation of turbulent fluxes.However, based on criterion (1), the results can only be applied to a relatively short time span, which is ineffective in practical use.In order to provide information on how much the turbulent flux should be corrected, the criteria ( 2) and ( 3) based on the whole data sample can be beneficial.Figure 8 shows the correction of the turbulent momentum flux at different heights under different complexity and intensity of non-turbulent motions.From the overall values, similar conclusions can be drawn that the non-turbulent motions tend to make the turbulent flux overestimated, and there are also non-significant cases (where the confidence interval crosses zero), implicating occasional underestimation.The comparison of values among different heights is also compatible with the results from criterion (1).
From the left panel of Figure 8, the degree of overestimation resulting from non-turbulent motions is kept at a low level, and increases as the complexity increases from R u,w = 0∼3 × 10 −3 ( N u,w = 0∼100 ) , and reaches a fluctuating state with a steady average value from 0.15 to 0.2 kg ⋅ m −1 ⋅ s −2 for larger R u,w .It can be explained by the results in Figure 4 and E quations (14a)-(14d), because the influence from non-turbulent motion within a monotonic segment is generally overestimation, and complex kinds of motions are able to overcome these effects together, making the overestimation for the whole data sample more serious.However, when the value of complexity is too high, the uncertainty and the cross-terms defined by Equations ( 14b) and (14c) prevent further overestimation.
From the right panel of Figure 8, the overestimation evidently increases as the intensity I u,w increases, which is a relatively intuitive result: it is harder to eliminate more intense kinds of non-turbulent motion by traditional calculation of turbulent flux, and thusthis can make the overestimating effect more serious.However, the uncertainty of the increasing pattern tends to be larger for large intensity, making it harder to give a unified standard for flux correction.
Figure 9 shows the correction of the turbulent sensible-heat flux at different height under different complexity and intensity of non-turbulent motions.The values and overall patterns are quite similar to those for momentum flux, yet the overestimation problem is less serious for w ′  ′ < 0, and its variation with complexity and intensity of non-turbulent motion is also less significant.Moreover, the uncertainty of correction with intensity also tends to increase.
Figure 10 shows the correction of the turbulent latent heat/CO 2 flux at 80 m under different complexity and intensity of non-turbulent motions.It can be judged that the pattern of increasing correction with complexity and intensity of non-turbulent motions can also be applied to latent heat and CO 2 flux, yet there are still subtle differences, such as that the error is larger for w ′ q ′ < 0, but also for w ′ c ′ ≥ 0. The correction uncertainty of these two kinds of fluxes also increases with intensity.
For practical use, the corrections with complexity and intensity revealed by

F I G U R E 10
The correction of turbulent latent heat (upper panels) and CO 2 (lower panels) flux with respect to the properties of non-turbulent motions.For the line types, please refer to Figures 8 and 9.
be calculated as functions of the original fluxes and the empirical relationships of R i, and I i, : The coefficients for different kinds of turbulent fluxes can be acquired from the results of lowest layer of observations (for the constant-flux layer in the usual sense; 40 m for momentum and sensible heat flux, and 80 m for latent heat and CO 2 flux), and are listed in Table . 1:

Instruction of correction in practice
The last two sections have proved that the complexity and intensity of non-turbulent motions can affect the calculation of turbulent flux.The values of correction are all defined by absolute change, as small values of original flux can disturb the relative change to a large extent, making it hard to summarize a general pattern.In practice, a standardized correction concerning the properties of non-turbulent motions can be put forward, based on the relationship between reconstructed and the original flux.can be linked with usual statistics, the correction of turbulent flux can be acquired conveniently.
Figure 11 reveals a significant positive correlation between the standard deviation of w-wind and I u,w , and a weak positive correlation between I u,w and R u,w is also proven.This makes it possible to correct the turbulent flux with the help of standard deviation by the method in Section 3 without having to perform a Hilbert-Huang transform.One can apply the empirical relationships in Section 5.2 with the help of the approximated I u,w and R u,w from the standard deviation of w-wind (which is calculated by usual methods, for example, the first way of pre-processing in Section 2.2), and then the corrected turbulent flux can be acquired.
From the comparison of two kinds of correction with the original turbulent momentum flux, it can be judged that the theoretical correction provides a precise turbulent based on fine calculations, and the degree of correction can vary significantly.In contrast, the empirical correction based on the standard deviation of the w-wind and fitted coefficients provides a more conservative correction, which means the degree of correction tends to be lower than the theoretical one.In practice, there are three choices: 1. Strict conduction of theoretical correction as put forward in Section 3. In this way, the corrected flux is the most reliable, but the cost of calculation can also be the highest (about 500 times that of the other choices); 2. Correction based on the fitted lines of the theoretical correction.One can directly multiply the original flux by the slope of the fitted line, and the calculated flux can be theoretically accurate, but with larger uncertainty; 3. Empirical correction based on the standard deviation and coefficients.This would turn out to be a balanced choice of theoretical significance the uncertainty.
For turbulent sensible heat, CO 2 and latent heat fluxes, the empirical corrections are also similar to theoretical corrections under the verification of cross-validations of other datasets, and the results are thus omitted.However, the empirical correction is not always reliable, as the correction of CO 2 and latent heat fluxes may lose accuracy for small magnitudes of fluxes, which results from the failure of the cmpirical relationship to tell the difference in R w, at large values of I w, .Therefore, a thorough evaluation of the correction methods should be conducted before making the final choice.

CONCLUSIONS
In this study, to deal with the problem that the calculation of turbulent flux from the eddy covariance system is usually affected by system bias resulting from non-turbulent motions, the method of Hilbert-Huang transformation is applied to analyze the data from EC systems, and a systematic algorithm is put forward to distinguish, separate, and eliminate the non-turbulent motion in the observational data.Specifically, the IMFs for both of the two elements of a certain kind of flux are taken into consideration, the PDF of each IMF's amplitude and the contribution distribution of its amplitude are compared through the JSD.
The IMFs with an exceptional value of JSD imply the presence of non-turbulent motions, and those IMFs with lower frequencies are all extracted from the original data, and considered as non-turbulent motion.The remaining IMFs are accordingly taken as turbulent time series, and are used to reconstruct the turbulent flux.
To investigate to which degree non-turbulent motions can affect the calculation of turbulent flux, the five-level turbulence data from the Tianjin 255-m meteorological tower from July 1 to August 31, 2017 are examined.The non-turbulent motions for all the data samples are extracted and classified by three kinds of criteria according to the monotonicity, complexity, and intensity.The comparison of reconstructed flux and original flux shows that non-turbulent motions tend to overestimate turbulent flux, despite there also being a few cases of underestimation.Different combinations of monotonicity of non-turbulent motions are irrelevant for the degree of overestimation, but the higher the complexity and the intensity are, the more serious is the effect of overestimation.
Based on these properties, an instruction on how to correct the turbulent flux is put forward.For a set of original data, an algorithm for extracting non-turbulent motions can be applied, and the empirical function concerning properties of the non-turbulent motion is given.With the coefficients in these functions given in advance, the reconstructed turbulent flux can be conveniently calculated, which is opportunate for both the simulation of atmospheric turbulence and the parameterization in meteorological and climate models.

F
I G U R E 1 The statistical relationship between frequency and Jensen-Shannon Divergence of the sensible-heat flux at 40 m (left), and two specific cases of D JS ∼0.1 (steady-dominant motion) and D JS ∼0.5 (burst-dominant motion).PDF, probability density functions.
and choose the smaller one as the statistic quantity, U; 5. Construct the statistics z = U−E(U) √ D(U)

F
The flow chart of identifying the turbulent and non-turbulent motions from properties of transport, taking w-wind in calculating momentum flux as an example.EMD, empirical mode decomposition; IMF, intrinsic mode function; JS, Jensen-Shannon; PDF, probability density functions.
In many studies, the turbulent flux, including sensible-heat flux H, momentum flux , CO 2 flux F c , F I G U R E 3 Comparison of the time series (left) and energy spectrum (right) of original and reconstructed u-wind data.The data used in this figure were obtained at 40 m, 01:00-02:00, on July 1, 2017.

F
Comparison of different kinds of newly defined flux at 40 m (taking the covariance of w-wind and potential temperature as an example) among different definitions.

F
I G U R E 5 The correction of turbulent momentum flux under different cases.Different combinations of monotonicity of u-and w-wind are distinguished by the styles of the box plots, and are marked in the legend; different cases of complexity and intensity of non-turbulent motions are distinguished by the groups of the box plot, and are marked on the x-axis.The middle line in the box plot denotes the median value, and the black dot denotes the timespan-weighted average.motions, and therefore once the non-turbulent motions are eliminated, the absolute values of reconstructed turbulent fluxes are lower.In addition, from the comparison by vertical level, the value of overestimation increases with height from 40 m to 160 m, and decreases at 200 m (beyond the constant-flux layer), which means the detected influence of non-turbulent motion is lower near ground height, which can be attributed to the friction effect from the underlying surface to the non-turbulent motions.

F
I G U R E 6The correction of turbulent sensible-heat flux under different cases.For the different styles, the groups on the x-axis, and the components of the box plots, please refer to Figure5.

F
I G U R E 8The correction of turbulent momentum flux with respect to the properties of non-turbulent motions (left panel for complexity and right for intensity).The different types of the lines denote the results at different height, the translucent patches denote the 95% significance intervals of the results at 40 m, and the dashed lines denote the fitted curves of empirical relationships.
Figures 8-10 can be fitted into empirical relationships, f ( R i, ) and g ( I i, ) , which are controlled by different sets of free coefficients A-D and E-G.Assuming that the influences of complexity and intensity bear equal significance, the reconstructed turbulent fluxes can F I G U R E 9 The correction of turbulent sensible-heat flux regarding the properties of non-turbulent motions.The lines above the dotted lines are the correction for H orig < 0, and the ones below the dotted lines for H orig ≥ 0. For the different line types, please refer to Figure 8. c The relationship among standard deviation of w-wind and I u,w (left top panel), I u,w and R u,w (left bottom panel), and corrected turbulent momentum flux and original flux (right panel).The hollow points and dashed fitted line in the right panel denote the method in Section 3, and the solid ones are for the empirical correction in this section.
The coefficients in empirical relationships of correction for different kinds of turbulent fluxes.
According to Table 1, the correction of turbulent flux can be determined by two properties of non-turbulent motions, R i, and I i, .It can be proposed that if these two parameters TA B L E 1