Optimized Estimation of Surface Layer Characteristics from Profiling Measurements

New sampling techniques such as tethered-balloon-based measurements or small unmanned aerial vehicles are capable of providing multiple profiles of the Marine Atmospheric Surface Layer (MASL) in a short time period. It is desirable to obtain surface fluxes from these measurements, especially when direct flux measurements are difficult to obtain. The profiling data is different from the traditional mean profiles obtained at two or more fixed levels in the surface layer from which surface fluxes of momentum, sensible heat, and latent heat are derived based on Monin-Obukhov Similarity Theory (MOST). This research develops an improved method to derive surface fluxes and the corresponding MASL mean profiles of wind, temperature, and humidity with a least-squares optimization method using the profiling measurements. This approach allows the use of all available independent data. We use a weighted cost function based on the framework of MOST with the cost being optimized using a quasi-Newton method. This approach was applied to seven sets of data collected from the Monterey Bay. The derived fluxes and mean profiles show reasonable results. An empirical bias analysis is conducted using 1000 synthetic datasets to evaluate the robustness of the method.


Introduction
Although the Marine Atmospheric Surface Layer (MASL) is normally less than 100 m in depth, the physical processes occurring in this layer determine the amount of momentum, heat, and water vapor exchange between the ocean and the atmosphere as part of the air-sea interaction processes, measured as momentum, heat, and water vapor fluxes, respectively.These three fluxes describe air-sea interaction processes at the interface and are some of the major factors in determining forecast accuracy at all time scales from the mesoscale to the climate time scales.Since 70% of the Earth's surface is covered by the ocean, methods to better quantify the surface fluxes are critical for weather forecast model improvement and validation and for understanding the climate and climate change [1].
In addition to obtaining surface fluxes, accurate depictions of the temperature, humidity, and wind profiles in the MASL has important applications such as wind energy forecasting, pollution, or hazardous gas dispersion through the atmosphere.It is also known that the vertical profiles of temperature and humidity in MASL are critical in quantifying the effects of the atmosphere on electromagnetic (EM) wave propagation, which affects all ranges of the EM spectrum including radar and communication signals.When traveling in the surface layer, EM waves can result in very complex interference and bending patterns [2].The signals may even be trapped in a wave-guide and focused along a path well beyond its normal range, a propagation phenomenon usually referred to as ducting.
The impact of the atmosphere on EM wave propagation is, in fact, the primary focus of many previous research efforts [3][4][5].
Surface fluxes over the ocean can be obtained in multiple ways.For direct flux measurements using the eddy correlation method, surface buoys or research vessels are normally used [6].This is, however, not always feasible given the complexity of the sensor system and the amount of high-rate sampled data for storage.Very often, ship-or buoy-based measurements only provide data at low sampling rates, resulting in mean wind, temperature, and humidity at a single level and measured sea surface temperature.These mean quantities are normally referred to as atmospheric and oceanographic bulk variables.Surface fluxes are then calculated using such measurements based on the relationships of the measured variables to the fluxes.The Coupled Ocean-Atmosphere Response Experiment (COARE) bulk algorithm is one of the bulk surface flux parameterization schemes for this purpose [7].The COARE scheme is the state-of-the-art surface flux scheme and the most frequently used algorithm in the air-sea interaction community.
Over land surfaces, profile measurements of the wind, temperature, and humidity in the atmospheric surface layer are often made on instrumented towers at more than two levels.When combined with turbulence measurements, the dataset is used to derive or validate the empirical surface layer flux-profile relationship [8].The profile measurements at more than two levels are also used to derive surface fluxes using the least square optimization method and the MOST theory [9,10].This method was referred to as the profile method.
As part of the Coupled Air-Sea Processes and EM ducting Research (CASPER) project, we made profile measurements from a small boat using a tethered balloon.These measurements provide multiple profiling measurements in the lowest 60 m of the MASL in a relatively undisturbed environment.Such a dataset is inherently different from those traditionally used to derive surface fluxes where measurements were at fixed level(s) that are capable of providing statistically significant mean measurements at each level.The profiling measurements we work with here may not be sufficient to provide a statistically meaningful mean value at a giving altitude or altitude bin.However, independent measurements are available at many vertical levels with high vertical resolution.The goal of this paper is to develop and test a new method for the estimation of the surface layer turbulence fluxes and the mean profiles using the profiling dataset.The optimized estimation method takes advantage of the large number of measurements from many vertical levels to maximize the use of observational data and therefore achieve high statistical significance.Here, we use a least-square optimization method.Although the basic concept of this method is similar to that of Nieuwstadt [9], a few significant new developments in the methodology of our approach is described to suit the different type of dataset used in this study.This method was tested on several sets of observational data.The results indicate that the method is accurate and reliable.Furthermore, the estimation process is empirically verified to be unbiased.

The Basic Theory
The Monin-Obukhov Similarity Theory (MOST) is the theoretical framework that governs the mean wind speed, potential temperature, and specific humidity profiles and their relationship to surface fluxes in the atmospheric surface layer.In a gradient relationship form, the MOST can be expressed in the following differential equations [8]: where z represents the height from the surface, and u, θ, and q are the mean wind speed, potential temperature, and specific humidity, respectively.In Equation ( 1), u ˚is the frictional velocity, θ ˚the temperature perturbation scale, and q ˚the water vapor perturbation scale.These scale quantities define the turbulence momentum, sensible heat, and latent heat fluxes.The Obukhov length, denoted by L, is defined as: where κ, the von Karman constant, is 0.4 following Fairall et al. [11], and g is gravity.The functions ϕ m pz{Lq, ϕ h pz{Lq, and ϕ q pz{Lq are the empirical non-dimensional gradients for momentum and scalars, respectively.Following Panofsky [12], the integrated form of Equation ( 1) can be written as Equation (3) for the mean wind, potential temperature, and specific humidity (u, θ, and q).Here, u 1 , θ 1 , and q 1 represent the mean wind speed, potential temperature, and specific humidity at their corresponding reference levels of z u1 , z θ 1 , and z q1 .
Ψ m , Ψ q , and Ψ h in Equation ( 3) are universal functions related to ϕ m pz{Lq, ϕ h pz{Lq, and ϕ q pz{Lq by: where ξ " z{L, and ξ 1 represents the ξ value corresponding to z u1 , z θ1 , and z q1 heights for wind, temperature, and humidity, respectively.Equation (3) uses the mean wind, potential temperature, and specific humidity as input variables instead of their gradients shown in Equation (1).Given the accuracy of the measurements of the bulk quantities, the use of Equation ( 3) is considered more accurate for flux estimates compared to Equation (1), as the gradients of the mean wind and scalars have much larger sampling error than the bulk quantities themselves [12].
When bulk measurements are available at a second level in addition to z u1 , z θ 1 , and z q1 , Equation (3) can be solved for u ˚, θ ˚and q ˚.This method is frequently used in obtaining surface fluxes.When only one level of measurements is available, surface layer fluxes can also be obtained by using the properties at the surface roughness heights for momentum, temperature, and humidity (z 0u , z 0θ , and z 0q ).Although surface roughness heights over land surfaces are prescribed, those over the ocean can be expressed as a function of wind speed or momentum flux.
In this study, many of the formulations and constants are taken from the COARE surface flux scheme.The basic structure of COARE is built on surface layer MOST described by Liu et al. [13] with modifications to include various formulations of the roughness/stress relations, the gustiness factor that modifies surface fluxes, and the extension of the surface layer scaling in convective limits in addition to effects of precipitation.Fairall et al. [11] gives a full description of the various physical processes affecting surface flux exchanges over the ocean and describes the various components of version 2.5 of the COARE.An improved version, version 3.0 of COARE, was released in 2003 and was described by Fairall et al. [7].Here, surface roughness lengths for momentum and scalars follow those in COARE 3.0, where the momentum roughness length, z 0u , is a combined formulation of the Charnock relationship for rough sea [14] and that of Liu et al. [13] for smooth sea.
The empirical momentum and scalar stability functions can take various forms from different studies.In this study, again, we followed COARE 3.0 on the universal stability functions, ϕ m,h pz{Lq, or its integrated form Ψ m,h pξq, where the unstable profile stability functions are a blend of those from the Kansas experiment [8] and those from the experiment performed by Grachev et al. [15] for the convective limit.A detailed description for the rationale of the choice of stability functions and constants are given in Fairall et al. [7,11].
Although the majority of applications of the MOST theory in obtaining surface fluxes use mean measurements or model output from one or two levels, there are cases when measurements from more than two levels are available.The additional data allows the use of optimization methods to provide flux estimates that have better statistical significance and are consistent with the constant flux surface layer.Two studies of this kind were mentioned in the introduction; Nieuwstadt [9] and Kramm [10] both belong to this category.Nieuwstadt [9] presented the basic framework of this method using mean measurements at 5-7 levels on a mast over an area with known surface roughness.Compared to the single-or two-level methods, the profile method utilizes a larger number of data in the flux estimation to ensure the statistical significance of the results and to provide a true "constant flux" layer solution.A variation of this method was presented in Kramm [10], where the multi-level measurements were used to derive surface roughness and displacement height in addition to fluxes.

The Testing Data Sets
In the marine environment, making measurements at multiple fixed levels in undisturbed air away from large ships or buoys is very difficult.However, information about the marine surface layer is much needed particularly for determining the effects of the marine surface layer on EM propagation.As part of the CASPER effort, data was collected from a Rigid Hull Inflatable Boat (RHIB) using radiosonde on a tethered balloon, where temperature and relative humidity profiles in the lowest 50 m were made repeatedly (Figure 1).Meanwhile, wind speed was measured at 2 m above the surface using a Vaisala WXT520 Weather Transmitter on a small mast attached to the RHIB.Seven independent sets of data collected on 2 May 2015 off the shore of Monterey Bay, CA are used here as test datasets.These seven sets of profiles were collected at different locations and times of the day within 13 km offshore.Measurements for each set were completed within ~30 min.
The empirical momentum and scalar stability functions can take various forms from different studies.In this study, again, we followed COARE 3.0 on the universal stability functions,  ,ℎ (   ⁄ ), or its integrated form Ψ ,ℎ (), where the unstable profile stability functions are a blend of those from the Kansas experiment [8] and those from the experiment performed by Grachev et al. [15] for the convective limit.A detailed description for the rationale of the choice of stability functions and constants are given in Fairall et al. [7,11].
Although the majority of applications of the MOST theory in obtaining surface fluxes use mean measurements or model output from one or two levels, there are cases when measurements from more than two levels are available.The additional data allows the use of optimization methods to provide flux estimates that have better statistical significance and are consistent with the constant flux surface layer.Two studies of this kind were mentioned in the introduction; Nieuwstadt [9] and Kramm [10] both belong to this category.Nieuwstadt [9] presented the basic framework of this method using mean measurements at 5-7 levels on a mast over an area with known surface roughness.Compared to the single-or two-level methods, the profile method utilizes a larger number of data in the flux estimation to ensure the statistical significance of the results and to provide a true "constant flux" layer solution.A variation of this method was presented in Kramm [10], where the multi-level measurements were used to derive surface roughness and displacement height in addition to fluxes.

The Testing Data Sets
In the marine environment, making measurements at multiple fixed levels in undisturbed air away from large ships or buoys is very difficult.However, information about the marine surface layer is much needed particularly for determining the effects of the marine surface layer on EM propagation.As part of the CASPER effort, data was collected from a Rigid Hull Inflatable Boat (RHIB) using radiosonde on a tethered balloon, where temperature and relative humidity profiles in the lowest 50 m were made repeatedly (Figure 1).Meanwhile, wind speed was measured at 2 m above the surface using a Vaisala WXT520 Weather Transmitter on a small mast attached to the RHIB.Seven independent sets of data collected on 2 May 2015 off the shore of Monterey Bay, CA are used here as test datasets.These seven sets of profiles were collected at different locations and times of the day within 13 km offshore.Measurements for each set were completed within ~30 min.An example of the surface profiles of temperature and humidity, and the derived modified index of refraction, M, is shown in Figure 2. The formulation and the physical meaning of M can be found in Turton [2].The scattered data points represent measurements from the tethered radiosonde.Concurrent with this profiling measurement, we also made wind, temperature, humidity, and pressure measurements from a small mast on the RHIB at an approximate height of 2 m together with water temperature measurements from a thermistor floating close to the surface.The water temperature measurements are used in place of sea surface temperature (SST).Figure 2 clearly shows the surface layer gradients in potential temperature and specific humidity, particularly in the lowest few meters.The resulting evaporative duct height at about 4 m is also clearly seen.Assuming the samples were taken from an ensemble of surface layer profiles where the mean follows the MOST theory, one would want to retrieve the MOST-following mean profiles and the associated surface fluxes.A simple approach is to take the mean values from a single altitude bin and the SST to obtain surface fluxes using the COARE algorithm.The two solid lines were generated MOST profiles with the COARE fluxes using measurements from the bin-averaged temperature and a humidity between 9 m and 11 m (blue line) and between 1 m and 3 m (cyan line).These two profiles represent different surface fluxes and would result in very different ducting characteristics.The differences in the results are likely caused by inadequate sampling at each level, resulting in the estimated mean deviating from the ensemble mean.In this example, there are limited numbers of data points, around 10-15, at each level.However, the measurements provide many data points at other levels that were not used by the traditional flux and profile reduction method, which uses two levels of measurements only.The approach derived from this research intends to improve the two-level data input methods to maximize the use of available observations in order to achieve higher statistical significance.
Atmosphere 2016, 7, 14 5 of 15 An example of the surface profiles of temperature and humidity, and the derived modified index of refraction, M, is shown in Figure 2. The formulation and the physical meaning of M can be found in Turton [2].The scattered data points represent measurements from the tethered radiosonde.Concurrent with this profiling measurement, we also made wind, temperature, humidity, and pressure measurements from a small mast on the RHIB at an approximate height of 2 m together with water temperature measurements from a thermistor floating close to the surface.The water temperature measurements are used in place of sea surface temperature (SST).Figure 2 clearly shows the surface layer gradients in potential temperature and specific humidity, particularly in the lowest few meters.The resulting evaporative duct height at about 4 m is also clearly seen.Assuming the samples were taken from an ensemble of surface layer profiles where the mean follows the MOST theory, one would want to retrieve the MOST-following mean profiles and the associated surface fluxes.A simple approach is to take the mean values from a single altitude bin and the SST to obtain surface fluxes using the COARE algorithm.The two solid lines were generated MOST profiles with the COARE fluxes using measurements from the bin-averaged temperature and a humidity between 9 m and 11 m (blue line) and between 1 m and 3 m (cyan line).These two profiles represent different surface fluxes and would result in very different ducting characteristics.The differences in the results are likely caused by inadequate sampling at each level, resulting in the estimated mean deviating from the ensemble mean.In this example, there are limited numbers of data points, around 10-15, at each level.However, the measurements provide many data points at other levels that were not used by the traditional flux and profile reduction method, which uses two levels of measurements only.The approach derived from this research intends to improve the two-level data input methods to maximize the use of available observations in order to achieve higher statistical significance.

Optimal Estimation Based on the Weighted Nonlinear Least-Squares
This section outlines a new estimation method and an optimization process to derive the mean profiles and the surface fluxes based on the profile measurements in the MASL.To formulate for the most general case, we assumed that all three variables were sampled at different altitudes and sampling rates.The following notations are used to represent the measurements of , , and  at a range of different heights:

Optimal Estimation Based on the Weighted Nonlinear Least-Squares
This section outlines a new estimation method and an optimization process to derive the mean profiles and the surface fluxes based on the profile measurements in the MASL.To formulate for the most general case, we assumed that all three variables were sampled at different altitudes and sampling rates.The following notations are used to represent the measurements of u, θ, and q at a range of different heights: z u " z u1 , z u2 , ¨¨¨, z un u ; u " u 1 , u 2 , ¨¨¨u n u z θ " z θ1 , z θ2 , ¨¨¨, z θn θ ; θ " θ 1 , θ 2 , ¨¨¨θ n θ z q " z q1 , z q2 , ¨¨¨, z qn q ; q " q 1 , q 2 , ¨¨¨q n q (6) where n u , n θ , n q are the numbers of measurements for the corresponding variables.The estimation process, which is based on least-squares and incorporates the functions defined in Equation (3), uses the following quadratic cost function: where, The weights are given as: Since the wind measurement was only made at one level, the reference level for wind in Equation ( 3), z u1 , can be set at the surface roughness height, z 0u , given in Equation ( 5), at which the mean wind approaches zero.For z θ 1 , and z q1 , we used the lowest heights at which θ and q were measured.
In Equation ( 8), the functions of u, θ, and q are defined in Equation (3).The term, Var, represents the variance of the corresponding variable.The datasets in Equation ( 6) have different sizes for u, θ, and q.In addition, the measured u, θ, and q have different standard deviations due to sensor accuracy and environmental variabilities.The weights in Equation (9) ensure that the cost function treats all variables equally despite different sizes.The weights are inversely proportional to the variances to avoid biased weighting on data with large variances.In the particular dataset used in this study, u is measured at a single height, n u " 1.The wind speed effects should be less significant than those measured at multiple altitudes.We hence made the W u half of the weight it would be if n u ą 1.In addition, each term in Equation ( 8) is also weighted based upon the height from the surface.This design of the weighting function reflects the fact that the gradients of all three variables are inversely proportional to height near the surface (near-neutral conditions), resulting in significant variation of the mean profiles with height near the surface, which is the most critical region in characterizing EM ducting conditions.With all the above considerations in the choices of weights in the cost function, an estimate of (u ˚, θ ˚and q ˚) is obtained by solving the problem of weighted nonlinear least-squares:

The Numerical Algorithm
Finding analytical solutions to the problem of least-squares in Equation ( 10) is impossible due to the nonlinear nature of the functions and integrals in Equations ( 7)- (9).Therefore, the problem was solved numerically in a MATLAB environment to find a local solution.We used a quasi-Newton method of minimization [16].Newton's method of minimization is a searching process from an initial guess of a minimum point.The cost function, J `u˚, θ ˚, q ˚, θ 1 , q 1 ˘, is approximated using its second-order Taylor expansion around the initial guess, and the approximated cost function is a quadratic polynomial.The resulting minimum point of the quadratic polynomial is assumed to be the new initial guess, and the process is repeated until the estimated minimum point is close to the true solution.In a quasi-Newton method, the Hessian in the Taylor expansion is approximated using a method that avoids the computation of the second order derivatives.The quasi-Newton method is available in the MATLAB toolbox.

Results and Discussions
The new method described in Section 2.2 for estimating turbulence fluxes and mean profiles in the atmospheric surface layer is applied to the RHIB-based measurements.In this section, we will first examine the results for surface fluxes and surface layer profiles for seven profiling measurements.The sensitivity of the results to the choice of stability functions will be discussed.We will also demonstrate the advantage of the height-dependent weighting function in this particular application.The final section includes an empirical verification of the unbiasedness of the estimation method.

Estimated Surface Layer Profiles and Fluxes
The optimal estimation method was applied to seven sets of profiling data collected on 2 May 2015 off the shore of Monterey Bay, CA.In each dataset, all measurements were made within 30 min and at altitudes less than 50 m above the surface.The largest dataset has over 1000 data points for each sampled variable.The optimization method was applied to each profiling dataset, and the surface layer turbulence characteristics, i.e., the frictional velocity, and perturbation scales of temperature and water vapor, were estimated.Table 1 gives the results for all profiling sets.The corresponding cost functions for each profiling dataset are also shown in Table 1 together with the resulting turbulence momentum, sensible, and latent heat fluxes.The surface layer fluxes are calculated as: Momentum flux " ´ρu 2 ˚, Heat flux " ´ρu ˚θ˚, Latent heat flux " ´ρu ˚q˚( 11) where the density of the air is given by ρ " 1.29 kg¨m ´3, the heat capacity under constant pressure C p " 1005 J¨K ´1¨kg ´1, and the latent heat of water vapor is L " 2.5 ˆ10 6 J¨kg ´1.
Table 1.Results of MASL scaling parameters, surface fluxes, and magnitude of the least-squares cost function.The momentum flux is shown in N¨m ´2, while the sensible and latent heat fluxes are shown in W¨m ´2.The units for u ˚, θ ˚, q ˚are m¨s ´1, K, and g¨kg ´1, respectively.Based on the estimated set of `u˚, θ ˚, q ˚, θ 1 , q 1 ˘, the profiles of mean potential temperature and specific humidity can be calculated using Equation (3).The resulted profiles for all seven cases are shown in Figure 3 together with the originally sampled data used in the optimization calculation.Since the true values of the mean profiles are unknown, the profiles derived from the optimization method are expected to fit the data well if the estimated fluxes are accurate.This comparison between the derived mean profile and the sample data points is used as a measure of the performance of the optimization method.The values of the cost function (7) in the least-squares can also be used to indicate the accuracy of the estimation, although this measure should be used in a comparative sense.The estimated profiles in Figure 3 are considered the best fit for the sampled data and are consistent with the MOST.By visual inspection, the profiles fit the general trend of the sampled data well, especially for specific humidity.Less fidelity is seen in the fit for potential temperature particularly for the Datasets 1-4, with Dataset 2 being the worst in which the potential temperature is clearly underestimated for the majority of the measurements.We further noticed that the first three The estimated profiles in Figure 3 are considered the best fit for the sampled data and are consistent with the MOST.By visual inspection, the profiles fit the general trend of the sampled data well, especially for specific humidity.Less fidelity is seen in the fit for potential temperature particularly for the Datasets 1-4, with Dataset 2 being the worst in which the potential temperature is clearly underestimated for the majority of the measurements.We further noticed that the first three datasets are all in low wind conditions where the mean wind speed is less than 2.5 m¨s ´1.Dataset 2 has the lowest mean wind speed of ~1 m¨s ´1.The low wind condition may be responsible for the less ideal performance of the method.The weak wind not only results in larger relative measurement error in wind speed, but is also known to present large uncertainty in the empirical stability functions.On the contrary, the best overall fits are seen in Datasets 5, 6, and 7, which had moderate mean speed of 4-6 m¨s ´1.We also found that while the temperature fit seems to be somewhat biased for some cases, the corresponding cost functions are not significantly different from others.This can be explained by the fact that the cost function represents the overall fit for all three variables (wind speed, potential temperature, and specific humidity).The wind speed profile is not shown here, as wind speed was only measured at a single level on the mast.

--
The flux data in Table 1 is consistent with the wind speed, where Dataset 2 gives the smallest momentum flux corresponding to its low wind speed.The largest sensible and latent heat fluxes are seen in Dataset 1, which is consistent with the largest air-sea temperature difference observed in this dataset among all sampled cases (not shown).
Our initial calculation used the original stability function as per Businger [8].The major difference between the COARE and the Businger stability functions is the extension of the surface layer scaling in convective limits.This difference affects the flux calculation in the low wind and unstable cases the most, which is the case for Dataset 1-3.Indeed, flux calculations using the Businger stability function gave unreasonably large sensible and latent heat fluxes for the three low-wind cases (not shown), while the results for Datasets 4-7 did not change significantly.Overall, the results from the COARE stability functions are more reasonable.

Choice of Weights and the Cost Function
In optimal estimation, the weights in a cost function are important to the accuracy and reliability of the method.Here, several factors were taken into consideration in determining the weights for the cost function defined in Equations ( 7)- (9).In meteorological systems, the variability of different variables may differ in magnitude, as they each represent a different aspect of the physical processes.A common approach, as seen in Van de Geer [17], is to use a weight that is inversely proportional to the variance of the data, which is the case in Equation (9).Another factor is the data size, as various sensors may be used to collect data at different sampling rates.This issue is solved with the averaging method in Equations ( 7)- (9).
The cost function in Equation ( 8) has a factor that is inversely proportional to the height of measurements for each sampled data point.This was found to be necessary given that the largest gradient in the atmospheric surface layer is close to the surface.This gradient layer is critical for identifying an evaporative duct for estimating the atmospheric effects on EM propagation.Meanwhile, a large portion of the data was measured at heights higher than 20 m.It is thus necessary to put more weight on samples from the lower levels.The 1/z factor in Equation ( 8) serves this purpose.Figure 4 shows an example of the mean profile calculation for Dataset 1 without the 1/z factor in the cost function.It is clear that the potential temperature is significantly over-estimated at nearly all altitudes.The cost function, on the other hand, is not high compared to the result using the cost function in Equation (8).Similar results were found for other profiling sets.
In the dataset for this testing, wind speed and direction was only measured at a single low altitude of about 2 m.This makes the wind data less effective for estimation than other variables.The cost function in Equation ( 8), with the 1/z factor, further helps to increase the weight of the wind measurements.The results are significantly improved with this refined cost function (Figure 3, Dataset 1).factor.The units for θ and  are the same as in Figure 3.

The Characteristics of the Estimator-An Empirical Analysis
Two important statistical characteristics of estimators of parameters are the bias and sampling variability.Bias refers to whether an estimator tends to either over-or underestimate the parameter.Sampling variability refers to how much the estimate varies from sample to sample.Nonlinear least-squares, in general, result in biased estimates.In this study, numerical experiments were carried out to analyze the bias of the newly-developed method.An example of this experiment is shown in this section.
We first assume the following properties of the surface layer: A sequence of heights that consists of  = 100 equally spaced data samples from  = 0.2 m to  = 50 m was generated.The values of  ̅,  ̅ , and  ̅ profiles at each sampling height were generated using the surface layer MOST theory outlined in Equations ( 1)-( 4) and the surface roughness in Equation (5).These profiles were treated as the true mean profiles in this numerical experiment.We then generated random noise at each height and added them to the true values of  ̅,  ̅ , and  ̅ at the corresponding heights.This procedure results in the "measured" data.The noise was generated using a Gaussian distribution with a variance comparable to those from observations (() = 0.2 m 2 •s −2 , () = 0.02 K 2 , () = 0.025 g 2 •kg −2 ).The optimal estimation program was then applied to the synthetic data to generate an estimation of ( * ,  * ,  * ,  ̅ 1 ,  ̅ 1 ).This process was repeated N times, where a family of  = 1000 datasets were generated.Each dataset consisted of  = 100 samples for each variable (wind velocity, temperature, and water vapor).The statistics of the estimated  * ,  * ,  * ,  ̅ 1 ,  ̅ 1 ) are summarized in Table 2.The histograms of each variable are shown in Figure 5.  factor.The units for θ and q are the same as in Figure 3.

The Characteristics of the Estimator-An Empirical Analysis
Two important statistical characteristics of estimators of parameters are the bias and sampling variability.Bias refers to whether an estimator tends to either over-or underestimate the parameter.Sampling variability refers to how much the estimate varies from sample to sample.Nonlinear least-squares, in general, result in biased estimates.In this study, numerical experiments were carried out to analyze the bias of the newly-developed method.An example of this experiment is shown in this section.
A sequence of heights that consists of n " 100 equally spaced data samples from z " 0.2 m to z " 50 m was generated.The values of u, θ, and q profiles at each sampling height were generated using the surface layer MOST theory outlined in Equations ( 1)-( 4) and the surface roughness in Equation (5).These profiles were treated as the true mean profiles in this numerical experiment.We then generated random noise at each height and added them to the true values of u, θ, and q at the corresponding heights.This procedure results in the "measured" data.The noise was generated using a Gaussian distribution with a variance comparable to those from observations (Var puq " 0.2 m 2 ¨s´2 , ar pθq " 0.02 K 2 , Var pqq " 0.025 g 2 ¨kg ´2).The optimal estimation program was then applied to the synthetic data to generate an estimation of (u ˚, θ ˚, q ˚, θ 1 , q 1 ).This process was repeated N times, where a family of N " 1000 datasets were generated.Each dataset consisted of n " 100 samples for each variable (wind velocity, temperature, and water vapor).The statistics of the estimated (u ˚, θ ˚, q ˚, θ 1 , q 1 ) are summarized in Table 2.The histograms of each variable are shown in Figure 5.  Due to the nonlinear nature of the problem, the estimators do not follow a Gaussian distribution.However, Figure 5 shows that all histograms are nearly bell-shaped, which is similar to a Gaussian distribution.From Table 2, both the mean and median of the estimation were almost identical to the true value.This implies that the bias of the estimation process is negligible.Table 2 also   Due to the nonlinear nature of the problem, the estimators do not follow a Gaussian distribution.However, Figure 5 shows that all histograms are nearly bell-shaped, which is similar to a Gaussian distribution.From Table 2, both the mean and median of the estimation were almost identical to the true value.This implies that the bias of the estimation process is negligible.Table 2 also shows the standard deviations, IQRs, and the bounds of estimates.The standard deviation and IQR of u ˚, θ 1 , q 1 are very small.The standard deviation and IQR of θ ˚, q ˚are not small, but satisfactory for meteorology applications.The box-and-whisker plots are shown in Figure 6.It is clear that the majority of the estimates of all variables are close to the true value.Due to the nonlinear nature of the problem, the estimators do not follow a Gaussian distribution.However, Figure 5 shows that all histograms are nearly bell-shaped, which is similar to a Gaussian distribution.From Table 2, both the mean and median of the estimation were almost identical to the true value.This implies that the bias of the estimation process is negligible.Table 2   The fluxes of momentum, heat, and latent heat were also estimated based on results of the N " 1000 synthetic datasets and the results are shown in Table 3.In Equation ( 11), the fluxes are nonlinear functions of (u ˚, θ ˚, q ˚), a fact that had some impact on the statistics of the estimation.For example, the estimation of heat flux has a small bias in its mean, which is about 5% of the true value.The estimate of the momentum flux shows a negligible bias.The bias of the latent heat flux is relatively small.The standard deviation and the IQR of the momentum fluxes are very small, implying that the estimates are reliable with little variability.The standard deviation and the IQR of the sensible and latent heat fluxes are relatively large.All histograms (Figure 7) have bell-shaped distributions.Relative to the true values, the variability shown in Figure 8 is small for the momentum flux, but not so small for the heat and the latent heat fluxes.Reducing the variability of the estimation is a part of our future research.This method can be used with any stability functions and formulations for surface roughness.In this study, we used the same stability function and roughness parameterization as in COARE version 3.0 to be consistent with many of the flux calculations in the literature.The method was applied to seven profiling datasets obtained during the pilot experiment of the Coupled Air-Sea Processes and EM ducting Research (CASPER) project off the shore of Monterey Bay.The results on the profiles and the momentum, heat, and latent heat fluxes indicate that the estimation fit the samples of the datasets reasonably well in moderate wind regime and less ideally in very weak winds.Testing of the method using the original Businger stability functions [8] gives very unreasonable fluxes in the weak wind cases and only small changes in the moderate wind cases.This suggests the importance of the stability function in this method as in all other MOST based flux calculation schemes.
To further understand the limitation of the method, an empirical bias analysis was made and revealed negligible or reasonably small bias in the estimation based on 1000 synthetic datasets with random noise.The estimation of the profile and the momentum flux also yielded small IQRs, indicating the reliability of the method.However, the heat and the latent heat fluxes had relatively large standard deviation and IQR values.Further research is needed to reduce the variability of the fluxes.
The advantage of the new method is its ability to use measurements from multiple levels at different sampling rates, maximizing the amount of sensor data that can be used.The optimization applies directly to the sample data, not their mean at a given altitude, which is different from previous studies using similar approaches [9,10].The method here is applicable to large datasets from profiling measurements in which data is collected independently for each variable and at many different altitudes, allowing us to obtain statistically significant profiles and fluxes.Hence, the improved method addressed in this paper facilitates the use of a much larger data source in deriving MASL flux and profile properties, which is not possible with the previous studies using the mean measurements from two or more fixed levels of measurements.
This new method opens up the possibility of obtaining surface fluxes and near surface profiles using new observational capabilities in various applications such as EM ducting research and air-sea interaction in hurricane and severe weather environments.With the development of the Global Positioning System (GPS) dropwind sonde [18], high-resolution wind and thermodynamic profile measurements down to the ocean surface have been made in hurricane conditions in the past 20 years.Powell et al. [19] used dropsonde measurements to obtain surface layer roughness length and hence the momentum drag coefficient from a least-square fit of the composite profiles of many soundings.However, the Powell et al. approach was limited to the assumption of neutral stability.The method presented here provides a viable alternative for flux calculation from similar measurements with the flexibility of different stability functions in an environment where in situ flux measurements near the surface are unlikely.Furthermore, improved sensor technology makes it practical to use systems of sensors that collect samples of data at a wide range of altitudes from small Unmanned Aerial System (UAS) or tethered balloons.These newly emerged environmental sensing technologies provide large amount of data in many different levels (or profiles) but not a sufficiently large amount of data in a given altitude.The current method is applicable to this type of surface layer profiling measurements.
The current study can be improved in several ways.First, the gustiness factor for low wind [20] was not included in the calculation.The addition of this process in the algorithm is expected to improve our results in a low-wind regime.Turbulent fluxes in the low wind conditions can also be affected by other factors such as surface waves [21], swells [22], and even the boundary layer height [23].Furthermore, our cases are limited to unstable thermal stratification; measurements in stable boundary layers should be tested.Such measurements are available from the recently completed CASPER-East field campaign and will be used in future efforts.Finally, a more comprehensive evaluation of the method should involve comparisons with in situ flux measurements.CASPER-East made concurrent flux measurements from nearby small buoys and ships, which can be used to compare with those derived from this optimization method.

Figure 1 .
Figure 1.RHIB-based surface layer profiling measurements offshore in Monterey Bay on 2 May 2015.The sensors and the corresponding measured variables are indicated in the image.

Figure 1 .
Figure 1.RHIB-based surface layer profiling measurements offshore in Monterey Bay on 2 May 2015.The sensors and the corresponding measured variables are indicated in the image.

Figure 2 .
Figure 2. The results of the MOST model using measurements at two levels.The blue and cyan lines are the corresponding results of input from mean values at 9 m and 11 m and at 1 m and 3 m, respectively.The sampled data are in red.The modified refraction, M, has no unit.

Figure 2 .
Figure 2. The results of the MOST model using measurements at two levels.The blue and cyan lines are the corresponding results of input from mean values at 9 m and 11 m and at 1 m and 3 m, respectively.The sampled data are in red.The modified refraction, M, has no unit.

Figure 3 .
Figure 3.The mean potential temperature and specific humidity profiles estimated with the optimization method (red solid lines) and the original sampled data (blue dots).The vertical axes are height from the surface in meters.The units for θ and  are Kelvin and g•kg −1 , respectively.

Figure 3 .
Figure3.The mean potential temperature and specific humidity profiles estimated with the optimization method (red solid lines) and the original sampled data (blue dots).The vertical axes are height from the surface in meters.The units for θ and q are Kelvin and g¨kg ´1, respectively.

Figure 4 .
Figure 4. Mean surface layer profiles obtained from Dataset 1 using a cost function without the 1  ⁄

Figure 4 .
Figure 4. Mean surface layer profiles obtained from Dataset 1 using a cost function without the 1/z factor.The units for θ and q are the same as in Figure3.

Figure 6 .
Figure 6.Box-and-whisker plots for (a)  * ; (b)  * ; (c)  * ; (d)  ̅ 1 ; and (e)  ̅ 1 showing the results based on 1000 sample datasets.Each plot shows the first and third quartile of the data, while the line in the center denotes the median.The two lines extending on either side of the box cover the data points within 2 × IQR and centered at the median.All remaining data points are shown as outliers.The unit for each variable is the same as in Figure 5.

Figure 6 .
Figure 6.Box-and-whisker plots for (a)  * ; (b)  * ; (c)  * ; (d)  ̅ 1 ; and (e)  ̅ 1 showing the results based on 1000 sample datasets.Each plot shows the first and third quartile of the data, while the line in the center denotes the median.The two lines extending on either side of the box cover the data points within 2 × IQR and centered at the median.All remaining data points are shown as outliers.The unit for each variable is the same as in Figure 5.

Figure 6 .
Figure 6.Box-and-whisker plots for (a) u ˚; (b) θ ˚; (c) q ˚; (d) θ 1 ; and (e) q 1 showing the results based on 1000 sample datasets.Each plot shows the first and third quartile of the data, while the line in the center denotes the median.The two lines extending on either side of the box cover the data points within 2 ˆIQR and centered at the median.All remaining data points are shown as outliers.The unit for each variable is the same as in Figure 5.

Table 2 .
Statistical characteristics of the estimator.The units are the same as in Table1.

Table 2 .
Statistical characteristics of the estimator.The units are the same as in Table1.

Table 3 .
Statistical characteristics of the fluxes.The units are the same as in Table1.