Sea-gas interface momentum flux calculation method coupled with sea wave data

Abstract

The invention discloses a sea-gas interface momentum flux calculation method coupled with sea wave data, and relates to the technical field of marine atmospheric boundary layer flux observation, and the method comprises the steps: carrying out the threshold screening and wild point removal processing of turbulence flux elements needed by the calculation based on a vortex correlation method, and carrying out the coordinate rotation of a three-dimensional wind speed; processing the obtained sea surface fluctuation data to obtain a sea wave spectrum and a wave parameter, and calculating a combined field peak parameter Qp based on the sea wave spectrum; on the basis of the Monin-Obhoff similarity theory, the Monin-Obhoff length, the atmospheric stability and the wind speed U10 of 10 meters on the sea are calculated; based on the data, the friction speed is split into a wind speed item and a wave mixing item for fitting, and finally the momentum flux is obtained. According to the method, key data such as the wave period and the wave direction of the sea waves are coupled into a calculation system, the Hendan peak parameter is introduced to represent the sea wave state, the influence of the sea waves on the momentum flux is comprehensively considered, the calculation precision is improved, and meanwhile the requirement for an observation instrument is lowered.

Description

Technical Field

[0001] This invention relates to the field of marine atmospheric boundary layer flux observation technology, and in particular to a method for calculating the momentum flux of the air-sea interface coupled with ocean wave data. Background Technology

[0002] In marine meteorological research, accurately calculating the momentum flux of the air-sea boundary layer is crucial for understanding the interaction between the ocean and the atmosphere. Accurate acquisition of this momentum flux contributes to a deeper understanding of the changing patterns of the marine climate system and has significant implications for many fields, including weather forecasting and marine environmental monitoring.

[0003] Currently, the commonly used methods for calculating air-sea boundary layer momentum flux include the eddy covariance method and various parameterization schemes. Among them, the eddy covariance method calculates momentum flux by directly measuring high-frequency wind speed fluctuations, making it the most accurate method for obtaining momentum flux data. However, the eddy covariance method has several limitations in practical applications. First, this method involves a large number of observations, requiring not only high-precision measurement of wind speed components but also accurate measurement of parameters such as water vapor density and carbon dioxide concentration. This places extremely high demands on the observation instruments, which are not only expensive but also require stringent accuracy and stability, increasing observation costs and maintenance difficulties. Second, in actual marine observation environments, especially in complex sea areas, there is a lack of stable setup space, and the method is often affected by marine environmental factors such as waves and currents, making it difficult to meet the observation conditions for the eddy covariance method, thus affecting the accuracy and reliability of the observation data.

[0004] Besides the eddy covariance method, the bulk parameter method is also a commonly used approach. COARE (Coupled Ocean-Atmosphere Response Experiment) 36 is a flux algorithm widely used in air-sea interaction research. Based on a series of physical parameterization schemes, it comprehensively considers various meteorological factors such as sea surface temperature, wind speed, and humidity to calculate momentum flux, heat flux, and water vapor flux at the air-sea interface. Compared to the eddy covariance method, its advantages lie in its relatively lower requirements for observation equipment and its lack of environmental limitations, enabling operational observations at more stations and providing a large amount of air-sea flux data for global ocean regions. However, COARE 36 also has some drawbacks. On the one hand, although it considers many meteorological factors, it does not adequately consider factors such as ocean waves when facing complex ocean environments. In reality, ocean waves have a significant impact on air-sea boundary layer momentum flux. Wave height, wave period, and wave direction parameters change sea surface roughness, thus affecting the magnitude of momentum flux, leading to discrepancies between the calculated results and actual conditions.

[0005] Furthermore, existing simple parameterization schemes typically only consider wind speed at a depth of 10 meters at sea, neglecting many other factors that significantly influence momentum flux. As research progresses, platform data verification reveals that factors such as actual ocean wave periods have a substantial impact on flux. Parameterization schemes relying solely on a single wind speed element cannot comprehensively and accurately reflect the true state of momentum flux in the air-sea boundary layer, leading to significant discrepancies between calculated and actual values, and thus failing to meet the needs of marine meteorological research and practical applications. Summary of the Invention

[0006] To overcome the aforementioned problems in the prior art, this invention proposes a method for calculating the momentum flux of the air-sea interface coupled with ocean wave data.

[0007] The technical solution adopted by this invention to solve its technical problem is: a method for calculating the momentum flux of the air-sea interface coupled with ocean wave data, comprising the following steps: Step 1: Threshold screening and outlier removal are performed on the turbulent flux elements required for calculation based on the eddy covariance method, and the coordinates of the three-dimensional wind speed are rotated. Step 2: Process the acquired sea surface undulation data to obtain the wave spectrum and wave parameters, and calculate the Aida peak parameter Q based on the wave spectrum. p ; Step 3: Calculate the Monin-Obukhov length, atmospheric stability, and 10-meter wind speed U at sea based on the Monin-Obukhov similarity theory. 10 ; Step 4: Based on the data obtained in Steps 1-3, the friction velocity is decomposed into a wind speed term and a wave mixture term for fitting, and finally the momentum flux is obtained. Friction speed in step 4 The specific calculation formula is as follows: ; ; ; Where T p For the period of the spectral peak, Let F be the angle between the wind and waves, and F be the wave term function. Q is the wave impact factor. p Here, A, B, C, and D are the peak parameters of the Aida peak; these are the fitted parameters. and The fitting parameters take different values ​​at different times.

[0008] The above-mentioned method for calculating the momentum flux at the air-sea interface coupled with ocean wave data, wherein step 1 specifically includes: Step 1.1: Collect data on ultrasonic false temperature diagnostic values, infrared gas analyzer diagnostic values, water vapor concentration, CO2 concentration, CO2 signal intensity, H2O signal intensity, three-dimensional wind speed, air temperature, and air pressure parameters, and perform threshold filtering and outlier removal on the collected parameter data. Step 1.2: Perform a first rotation of the three-dimensional wind speed in the xy-plane around the z-axis, so that the xz-plane is aligned with the average wind direction, and the average components... Satisfaction = 0; Step 1.3: Rotate the xz plane obtained in Step 1.2 around the y-axis to make the average components... =0, so that the average wind speed u satisfies .

[0009] The above-mentioned method for calculating the air-sea interface momentum flux coupled with ocean wave data, wherein the outlier removal principle in step 1.1 is as follows: Calculate the outlier for each parameter data type within a 30-minute time period; calculate the difference Δx between adjacent time series data, and use Δx to obtain the overall standard deviation σΔx of the differences between adjacent points; compare all data differences with the standard deviation one by one; if a data point satisfies the formula: Δx ≥ 6* σΔx, Then directly remove that point; compare all data differences with the standard deviation one by one, if a data point satisfies the formula: Δx ≥3 If *σΔx is calculated, the point is marked as an outlier. If 5 or more consecutive points are marked as outliers, they are marked as normal data. All outliers are evaluated, and data values ​​adjacent to an outlier that differ from the outlier's value by 30% are considered outliers. All outliers are removed and linear interpolation is performed to supplement them. If there are too many outliers within a 30-minute time period, the entire time period is removed.

[0010] The above-mentioned method for calculating the air-sea interface momentum flux coupled with ocean wave data, wherein step 2 involves the peak parameter Q. p The calculation method is as follows: ; in, It is the wave power spectral density function The zeroth moment (reflecting the total wave energy); The wave frequency.

[0011] The above-mentioned method for calculating the momentum flux at the air-sea interface coupled with ocean wave data, wherein step 3 is specifically calculated using the following formula: ; ; ; ; ; Where κ is the von Kármán constant. U is the friction speed. z Z represents the average wind speed at a height z above the sea surface, and z0 represents the sea surface roughness. For stability general functions, This is a dimensionless profile formula, where L is the Moning-Obukhov length. where is the average air temperature; g is the acceleration due to gravity; The horizontal line represents the Reynolds average, and the apostrophe indicates the turbulent fluctuation of the physical quantity.

[0012] The above-mentioned method for calculating the momentum flux of the air-sea interface coupled with ocean wave data, wherein... The specific calculation formula is as follows: ; Where a, b, c, and d are fitting parameters, reflecting the 10-meter wind speed U at sea. 10 The basic relationship of influence on friction speed and The fitting parameters take different values ​​at different times.

[0013] The beneficial effects of this invention are: (1) Improved calculation accuracy. Traditional parameterization schemes based solely on wind speed, due to neglecting key factors such as ocean waves, result in significant deviations between the calculated results and the actual momentum flux. In contrast, this invention, by coupling ocean wave data, comprehensively considers the influence of ocean wave elements on momentum flux and calculates friction velocity using different formulas segmented according to the Qp value, thus more accurately simulating the true momentum flux situation of the air-sea boundary layer. Verification using platform data shows that, compared to existing technologies, this scheme significantly improves the accuracy of momentum flux calculation results, providing a more reliable data foundation for marine meteorological research and contributing to a deeper understanding of the intrinsic mechanisms of air-sea interaction.

[0014] (2) Reduced observation difficulty and cost. Compared with the eddy covariance method, the number of observations required by this invention is significantly reduced. The eddy covariance method requires high-frequency and precise measurement of multiple parameters, which places extremely high demands on the observation instruments, making the equipment expensive and complex to maintain. In contrast, this invention only requires the acquisition of relatively conventional data such as 10m wind speed, wind direction, wave parameters, and wave spectrum over 30 minutes, greatly reducing the requirements for observation instruments. This enables the calculation of air-sea boundary layer momentum flux in more ocean areas, especially complex sea areas like this one, at a lower cost, improving the practical feasibility and universality of the scheme.

[0015] (3) Enhanced adaptability to complex marine environments. Marine environments are complex and variable, with significant differences in wave conditions. This invention utilizes Q... p As a segmentation criterion, it can flexibly handle momentum flux calculations under different wave conditions. Whether in relatively calm (Q) waves...p <2) The waves are still quite violent (Q) p In cases where ≥2), the friction velocity can be accurately calculated using the corresponding formulas, thereby obtaining reliable momentum flux results. This effectively overcomes the problem of inaccurate calculation results in some existing schemes when facing complex wave conditions. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the air-sea interface momentum flux of the present invention; Figure 2 This is the polynomial fitting result of the relationship between wind speed at 10m sea and friction speed in the embodiment of the present invention, where (a) is Q p <2 Fitting results; (b) is Q p ≥2 fitting results; Figure 3 This is the power function fitting result of the relationship between wind speed at 10m sea and friction speed in the embodiment of the present invention, where (a) is Q p <2 Fitting results; (b) is Q p ≥2 fitting results; Figure 4 This is the result of fitting an exponential function to the relationship between wind speed at 10m sea level and friction speed in the embodiment, where (a) is Q p <2 Fitting results; (b) is Q p ≥2 fitting results; Figure 5 The vertical velocity spectrum (blue line) during the period from January 29 to February 3, 2024, in this embodiment of the invention. ≥2, red line <2) Relationship diagram with wave spectrum (black line); Figure 6 The vertical velocity spectrum (blue line) during the period from February 17 to February 22, 2024, in this embodiment of the invention. ≥2, red line <2) Relationship diagram with wave spectrum (black line); Figure 7 This is the polynomial fitting result of function F in the embodiment of the present invention, where (a) is Q. p <2 Fitting results; (b) is Q p ≥2 fitting results; Figure 8 This is the fitting result of the power function F in the embodiment of the present invention, where (a) is Q. p <2 Fitting results; (b) is Q p ≥2 fitting results; Figure 9 This is the fitting result of the exponential function F in the embodiment, where (a) is Q. p <2 Fitting results; (b) is Qp ≥2 fitting results. Detailed Implementation

[0017] To enable those skilled in the art to better understand the technical solution of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0018] This embodiment discloses a method for calculating the momentum flux of the air-sea interface coupled with ocean wave data. A schematic diagram of the air-sea interface momentum flux is shown below. Figure 1 As shown, Q p It is a parameter used to classify sea states, when Q p When Q < 2, it corresponds to a wide frequency spectrum and a flat peak value, meaning that the waves are composed of multiple frequency components (such as waves caused by the coexistence of wind waves and swells, or waves after complex refraction / diffraction near the shore); when Q p When the frequency is ≥2, the waves exhibit a narrow spectrum and sharp peaks in the frequency domain, which means that the wave energy is concentrated near a certain dominant frequency (such as pure swells or single-peak strong waves caused by typhoons).

[0019] The calculation method in this embodiment specifically includes: Step 1: Threshold screening and outlier removal are performed on the turbulent flux elements required for calculation based on the eddy covariance method, and the coordinates of the three-dimensional wind speed are rotated.

[0020] This embodiment utilizes data acquired by the observation system aboard the Yangjiang platform, leveraging the eddy covariance observation system (including one CSAT3A three-dimensional ultrasonic anemometer, an EC150 carbon dioxide and water vapor analyzer, one PTB110 atmospheric pressure sensor, and an HMP155 thermometer) and the underwater wave observation system (AD2CP). The three-dimensional ultrasonic anemometer provides the 10Hz wind speed components (u, v, w) and ultrasonic virtual temperature Ts required for flux calculation; the carbon dioxide and water vapor analyzer observes the water vapor and carbon dioxide densities in the measured area; the atmospheric pressure sensor provides the atmospheric pressure value P in the measured area; and the thermometer monitors the atmospheric temperature T in the area. Both observation systems maintained their attitude during the experiment, acquiring meteorological and wave data for a total of 31 days.

[0021] For parameters such as ultrasonic false temperature diagnostic values, infrared gas analyzer diagnostic values, water vapor concentration, CO2 concentration, CO2 signal intensity, and H2O signal intensity, their values ​​are screened and evaluated to check for data anomalies and remove values ​​exceeding physically reasonable thresholds; the outlier removal principles adopted are as follows: 1. Perform calculations for each type of parameter data within a 30-minute time period; 2. Calculate the difference Δx between adjacent time series data, and use Δx to obtain the overall standard deviation (σΔx) of the difference between adjacent points.

[0022] 3. Compare all data differences with the standard deviation one by one. If a data point satisfies the formula: Δx≥6*σΔx, then remove the point directly. 4. Compare all the data differences with the standard deviation one by one again. If a data point satisfies the formula: Δx≥3*σΔx, then mark the point as an outlier. 5. If five or more consecutive points are marked as outliers, then mark them as normal data. 6. Identify all outliers. Data points adjacent to an outlier that differ from the outlier's value by 30% are considered outliers. 7. Remove all outliers and perform linear interpolation to fill in the gaps; 8. If there are too many wild points within a 30-minute time period, remove that time period directly.

[0023] Outlier removal is performed on three-dimensional wind speed (u, v, w), ultrasonic virtual temperature, water vapor concentration, CO2 concentration, air temperature, and air pressure. When the eddy covariance system is working, environmental factors such as rain, snow, and dust particles, or unstable power supply voltage and power outages can interfere with the sensors, resulting in large instantaneous noise, which are outliers. Outliers may significantly affect the results of data calculation of variance and covariance, thus affecting the flux calculation results.

[0024] The three-dimensional wind speeds (u, v, w) are rotated according to coordinates. Since the original turbulence data's wind speed components are measured by an ultrasonic instrument (such as CSAT3), the data corresponding to the three components are not in the natural coordinate system, but rather in the so-called 'ultrasonic coordinate system.' If the instrument is tilted, it will severely affect the measurement accuracy of the three components. This error is eliminated by a double rotation: 1. Achieve the first rotation of the xy-plane around the z-axis using a rotation matrix, aligning the xz-plane with the mean wind direction, and the average components... Satisfaction = 0: 2. The second rotation is a new rotation of the xz plane around the y-axis, which in turn affects the average components. =0: Ultimately, the average wind speed u is satisfied. .

[0025] Step 2: After acquiring sea surface undulation data for 1024 seconds per hour, the data is processed using Nortek's SignatureWaves post-processing software to obtain the wave spectrum and wave parameters (wave period, wave direction). Based on the wave spectrum, the Aida peak parameter Q is calculated. p .

[0026] Wave spectra contain a wealth of wave information, Q pAs the Goda peak parameter, it can effectively characterize the wave state. Due to the complexity and variability of wave states, a single parameter is insufficient to fully reflect its impact on momentum flux. Qp integrates key information from the wave spectrum and can serve as an effective segmentation criterion to distinguish the ways in which different wave conditions affect momentum flux.

[0027] ; in: It is the wave power spectral density function The zeroth moment (reflecting the total wave energy) is calculated using the following formula; The wave frequency.

[0028] ; Q p = 2 is the segmentation standard. In the following calculations of wind speed and wave terms, piecewise functions are used for calculation.

[0029] Step 3: Calculate the Monin-Obukhov length, atmospheric stability, and 10-meter wind speed U at sea based on the Monin-Obukhov similarity theory. 10 .

[0030] Based on the Monin–Obukhov (MO) similarity theory (MOST), the wind speed at the observation height is converted into a 10-meter wind speed U by applying a logarithmic wind profile. 10 This process requires atmospheric stability correction, and the conversion formula is as follows: ; ; In the formula, κ = 0.4 is the von Karman constant, and U z Let z be the average wind speed at a height z above the sea surface. Let z be the air friction velocity, and z0 be the sea surface roughness. This is a general stability function, which has a value of 0 under near-neutral conditions. Its specific expression is as follows: ; ; Where L is the Monin-Obukhov length (MO length), its expression is: ; In the formula: ρ is the average air temperature; g is the acceleration due to gravity, taken as 9.8 m / s². 2 ; The horizontal line represents the Reynolds average, and the apostrophe indicates the turbulent fluctuation of the physical quantity.

[0031] Quality control of the above data includes: 1. Removal of wind shadow area data Because the airflow from behind the ultrasonic anemometer (wind shadow area) is affected by the tower structure during observation, the uncertainty of the data increases. Based on the anemometer's installation orientation, observation data within the 225°–315° range are discarded. 2. Stability parameter selection Based on the Moning-Obukhov similarity theory, an excessively large stability parameter z / L (the ratio of observation height to Moning-Obukhov length) can lead to unreasonable wind speed conversion results (such as negative wind speeds at low wind speeds). Referring to existing research standards, data with z / L > 2 and z / L < -2 were excluded. 3. Low wind speed data removal When the wind speed is below 1 m / s, the uncertainty of wind direction observation increases significantly, and the relevant data cannot accurately reflect the wind stress characteristics, so U is directly discarded. 10 Low wind speed data <1m / s; 4. Rainfall data removal During rainfall, raindrops passing through the probe's measurement path may affect the ultrasonic propagation time, thus causing errors. Therefore, the rainfall rate P is directly discarded. rain Low wind speed data ≥0.5m / s; 5.95 confidence interval screening With friction speed With a wind speed of 10 meters U 10 Using the second-order polynomial fitting curve as a benchmark, a 95% confidence interval was defined, and outlier data points outside the interval were removed.

[0032] Step 4: Based on the data obtained in Steps 1-3, the friction velocity is decomposed into a wind speed term and a wave mixture term for fitting, and finally the momentum flux is obtained.

[0033] The momentum flux at the air-sea interface is often represented by wind stress, and the formula for calculating wind stress is: ; ; in, air density, This is the drag coefficient.

[0034] The calculation method in this embodiment uses friction speed. The momentum flux is calculated by separating the momentum flux into wind speed and wave mixture terms and fitting them separately.

[0035] 1. Wind speed calculation Numerous studies have found that the drag coefficient is primarily controlled by sea wind speed and is typically expressed as a linear function of wind speed. This embodiment first uses the least squares method to fit the relationship between the 10m sea wind speed and the friction velocity. The data used are ultrasonic wind speed data from an ultrasonic anemometer mounted on the Yangjiang platform and friction velocity calculated based on the eddy covariance method. Furthermore, based on MOST, the ultrasonic wind speed data is converted to a 10m wind speed U by applying a logarithmic wind profile. 10。 In the fitting stage, the Levenberg-Marquardt algorithm is used to iteratively find the optimal solution for the parameters. Fitting forms include polynomials, power functions, and exponential functions. The obtained parameter estimates satisfy the unbiasedness, efficiency, and consistency requirements of the Gauss-Markov theorem. Finally, the results are calculated based on the coefficient of determination (R²). 2 The optimal results are selected by using quantitative evaluation indicators such as root mean square error (RMSE), residual distribution law, and 95% confidence interval of parameters. Finally, the optimal result is obtained by polynomial fitting.

[0036] The results of different fitting forms are shown in Figure 1 and Table 1. To determine the optimal fitting relationship between variables, this study constructed three fitting models: polynomial, power function, and exponential function, and used the coefficient of determination (R²). 2 Ri and root mean square error (RMSE) are used as evaluation metrics for model fit. Among them, Ri... 2 The R-squared value represents the model's ability to explain data variation; a value closer to 1 indicates a better fit. RMSE reflects the average deviation between the fitted values ​​and the measured values; a smaller value indicates higher model fitting accuracy. The fitting evaluation indices for the three models are shown in Table 1. As can be seen from the table, the R-squared value of the polynomial fitting model... 2 In Q p When <2, it is 0.6584, Q p When the value is ≥2, the value is 0.4092, which is higher than that of the power function model and the exponential function model. This indicates that the model has the strongest explanatory power for the data variation pattern. Therefore, in this embodiment, the polynomial fitting model is selected as the representation model for the relationship between variables. ; Where Q p When Q is greater than or equal to 2, a = -0.002802, b = 0.03832, c = -0.118, d = 0.2023; when Q p When the value is less than 2, a = -0.0006677, b = 0.0088, c = 0.02284, d = 0.009237, to accommodate different Q values. p The influence of wind speed on friction velocity is determined. The initial friction velocity is calculated based on the above formula. These parameters were derived by fitting in-situ observation data from Yangjiang with theoretical research. The fitting process used the least squares method as the core optimization objective, and the Levenberg-Marquardt algorithm was employed to achieve optimal parameter estimation for the nonlinear parameter model. This satisfies the unbiased, efficient, and consistent requirements of the Gauss-Markov theorem, reflecting the 10-meter wind speed at sea. The basic relationship between friction speed and friction velocity.

[0037] Table 1

[0038] 2. Wave Term Calculation Introducing the peak factor Q p As a criterion, the sample was further subdivided into narrowband surges (Q). p ≥2 (blue curve) and broadband surge (Q p <2, red curve) two groups, the results are as follows Figure 5 and Figure 6 As shown.

[0039] First, the spectral peak frequency determines the vertical attenuation scale of wave disturbances, which is a geometric prerequisite for whether the wave signal can reach the observation height. (Comparison time period 1) Figure 5 ) and time period 2 ( Figure 6 The data shows that although both have high wave energy, the vertical velocity spectrum of period 2 follows the classic Kolmogorov-5 / 3 decay law and has no obvious spectral peak, indicating that the wave was dominated by high-frequency short waves that day and the wave boundary layer (WBL) thickness was insufficient to cover the sensor height; while the low-frequency swell (about 0.15 Hz) of period 1 decayed more slowly and successfully transmitted the wave signal to the observation layer.

[0040] Detailed spectral data from February 17 to February 22 ( Figure 6 In ), Q p The turbulent energy spectrum of ≥2 sample groups shows a clear peak at fp; in contrast, the Q spectrum of the same period shows a clear peak at fp. p <2 sample groups, although still dominated by swells, had their signals "flattened" in the vertical velocity spectrum, failing to form significant spikes and blending well with the background turbulence spectrum; during February 17 to February 22, even Q p For the sample group with ≥2, the peak intensity did not change significantly, which may be related to the lower absolute energy level or wider frequency band of the surge during that period.

[0041] Therefore, the wave peak period is related to Q. p As the core of this parameterization, the wave term is represented by the function F. The Levenberg-Marquardt algorithm is also used in the parameter solution stage to perform iterative optimal parameter calculation. The obtained parameter estimates satisfy the requirements of the Gauss-Markov theorem. Finally, based on the coefficient of determination (R²), the parameters are determined.2 The optimal parameter solution is obtained by using quantitative evaluation indicators such as root mean square error (RMSE), residual distribution law, and 95% confidence interval of parameters. Among these, when Q... p When Q is greater than or equal to 2, A = 0.0001176, B = -0.003764, C = 0.03181, D = -0.05076; when Q p When the result is less than 2, A = 0.001136, B = -0.0154, C = 0.04709, D = -0.01002 (fitting results are as follows). Figure 7-9 (as shown) ; ; Where T p For the period of the spectral peak, The angle between wind and waves, its value is between 0 and... 180°, if the wind direction is clockwise with the wave direction, it is positive; otherwise, it is negative. When the angle is small, the value of cos(angle_off) is close to 1, making the influence of the wave term relatively large. Combined with Q... p The influence of Q on the frictional velocity of the wave term reflects the effect of Q. p The modulation effect of the value on the wave.

[0042] 3. Calculation of final friction velocity ; The friction velocity obtained from the preliminary calculation Adding this to the wave term F yields the final frictional velocity. This result comprehensively reflects the combined effects of wind speed and ocean waves on momentum transport in the air-sea boundary layer, providing a more accurate picture of the actual air-sea boundary layer momentum flux.

[0043] In summary, by obtaining specific wind speed, wind direction, wave direction, wave period, and stability measurement data, the flux in the same area can be converted using the method of this embodiment.

[0044] To verify the effectiveness of the method in this embodiment, flux data calculated based on the eddy covariance method and flux data calculated based on COARE36 are also provided.

[0045] 1. Flux data calculated based on the eddy covariance method The eddy covariance method divides the measured quantity into average quantity and fluctuating quantity through Reynolds averaging. Then, it calculates the covariance between the concentration fluctuating quantity of different observed elements and the vertical wind speed component w fluctuating quantity to calculate the turbulent flux. Therefore, it is necessary to first calculate the covariance matrix by combining each component to obtain the covariance between each parameter and the vertical wind speed component w, as well as the mean of each component. (1) Momentum flux Tau (N / m 2 The calculation formula is: ; in, Friction speed (m / s) can be expressed as: ; in, Let be the covariance of the u-component and w-component of the three-dimensional ultrasonic wind. Let V be the covariance between the V and W components of the three-dimensional ultrasonic wind. ρ is the air density (kg / m³) 3 The following formula can be used to calculate the result: ; Where P is air pressure (Pa), T is air temperature (°C), and R is air pressure (Pa). d = 287 J / mol / K is the gas constant, e is the saturated water vapor pressure, calculated by the following formula: ; Where H2O is the water vapor density and Ts is the ultrasonic virtual temperature (°C). (2) Sensitive heat flux H (W / m 2 This can be represented as: ; In the formula C represents the covariance of the ultrasonic virtual temperature Ts and the three-dimensional ultrasonic wind w component. p The specific heat at constant pressure can be calculated using the following formula: ; In the formula C pd = 1004.67 (J / kg / K) is the specific heat at constant pressure of dry air, and q is the specific humidity, calculated by the following formula: .

[0046] (3) Ultrasonic virtual temperature correction: The initial sensible heat flux was calculated by subtracting the covariance between the ultrasonic imaginary temperature fluctuation and the vertical wind component. The ultrasound virtual temperature is obtained, but there is a certain deviation between the actual temperature and the virtual temperature, which needs to be verified. The final sensible heat flux H can be derived by reverse calculation: ; .

[0047] 2. Flux data calculated based on COARE36 COARE (Coupled Ocean-Atmosphere Response Experiment) 36 is a flux algorithm widely used in air-sea interaction research. Based on a series of physical parameterization schemes, it comprehensively considers various meteorological factors such as sea surface temperature, wind speed, and humidity, and uses complex formulas and algorithms to calculate momentum flux, heat flux, and water vapor flux at the air-sea interface. To use COARE 36 to calculate flux data, the following input data is required: sea surface temperature, wind speed, wind direction, atmospheric temperature and humidity, and atmospheric pressure P.

[0048] The flux calculated based on the eddy covariance method was used as the true value and compared with COARE36 and the proposed parameterization scheme. The results from January 10 to February 11, 2024 (a total of 33 days) were analyzed. The evaluation indicators for the calculation results included the coefficient of determination (R²), mean absolute error (MAE), and root mean square error (RMSE).

[0049] Coefficient of determination (R) 2 The value is used to measure the goodness of fit of the model to the observed data. The closer the value is to 1, the stronger the model's ability to interpret the data, that is, the higher the correlation between the calculated results and the true values.

[0050] ; Where n is the total number of data points. The observed value (here referring to the true flux value calculated based on the eddy covariance method). This is the predicted value (the flux value recalculated by this parameterization scheme or COARE36). That is, the average of the observed values. Table 2 shows the coefficient of determination (R²) of friction velocity calculated by COARE36 and the method of this embodiment with respect to the observed values. 2 As can be seen from Table 2, the coefficient of determination (R²) of the method in this embodiment is... 2 The values ​​are all higher than COARE36, which indicates that the structure calculated by the method in this embodiment is more correlated with the true value.

[0051] Table 2

[0052] The mean absolute error (MAE) reflects the magnitude of the average error between the predicted and actual values. The smaller the value, the closer the model's prediction is to the actual value.

[0053] ; Where n is the total number of data points. The observed value (here referring to the true flux value calculated based on the eddy covariance method). The values ​​are predicted (flux values ​​recalculated by this parameterization scheme or COARE36). Table 3 shows a comparison of the mean absolute error (MAE) between the friction velocity calculated by COARE36 and this scheme and the observed values. As can be seen from Table 3, the mean absolute error (MAE) of the method in this embodiment is smaller than that of COARE36, which indicates that the prediction results of the method in this embodiment are closer to the true values.

[0054] Table 3

[0055] The root mean square error (RMSE) not only considers the average size of the error, but also measures the degree of fluctuation of the error. The smaller the value, the higher the stability and accuracy of the model prediction.

[0056] ; Where n is the total number of data points. The observed value (here referring to the true flux value calculated based on the eddy covariance method). The values ​​are predicted values ​​(flux values ​​recalculated by the method in this embodiment or COARE36). Table 4 shows a comparison of the root mean square error (RMSE) between the friction velocity calculated by COARE36 and the method in this embodiment and the observed values. As can be seen from Table 4, the root mean square error (RMSE) of the method in this embodiment is smaller than that of COARE36, which indicates that the method in this embodiment has higher prediction stability and accuracy.

[0057] Table 4

[0058] The above embodiments are merely exemplary embodiments of the present invention and are not intended to limit the present invention. Those skilled in the art can make various modifications or equivalent substitutions to the present invention within its scope and spirit, and such modifications or equivalent substitutions should also be considered to fall within the scope of protection of the present invention.

Claims

1. A method of calculating momentum flux at the air-sea interface coupled with sea wave data, characterized by, The method comprises the following steps: Step 1, threshold screening and outlier processing are performed on the turbulent flux elements required for calculation based on the eddy correlation method, and coordinate rotation is performed on the three-dimensional wind speed; Step 2, process the acquired sea surface elevation data to obtain sea wave spectrum and wave parameters, and calculate the significant wave height parameter Q based on the sea wave spectrum p ; Step 3, calculation of Monin-Obukhov length, atmospheric stability, sea 10-meter wind speed U based on Monin-Obukhov similarity theory 10 ; Step 4, based on the data obtained in steps 1-3, the friction velocity is divided into a wind speed term and a wave mixing term for fitting, and finally the momentum flux is obtained; The friction velocity in step 4 The specific formula is: ； ； ； where T p is the spectral peak period, is the wind wave fetch angle, F is the wave term function, W p is the sea wave influence factor, Q p is the combined peak value parameter; A, B, C, and D are fitting parameters, and the fitting parameters have different values.

2. The method of claim 1, wherein, The step 1 specifically comprises: Step 1.1, collecting ultrasonic virtual temperature diagnostic values, infrared gas analyzer diagnostic values, water vapor concentration, CO2 concentration, CO2 signal intensity, H2O signal intensity, three-dimensional wind speed, air temperature, air pressure parameter data, and performing threshold screening and outlier processing on the collected parameter data; Step 1.2, First rotation of the three-dimensional wind speed in the x-y plane about the z axis to align the x-z plane with the mean wind direction, the mean component satisfies = 0; Step 1.

3. Rotate the x-z plane obtained in step 1.2 around the y axis such that the average component = 0, achieving that the average wind speed u satisfies .

3. The method of claim 2, wherein the method is characterized by, The outlier principle in step 1.1 is: for each parameter data type within a 30-minute period, calculate one by one; calculate the difference Δx of adjacent time series data, and obtain the overall standard deviation σΔx of the difference of adjacent points from Δx; compare all data differences with the standard deviation one by one, if a data point satisfies the formula: Δx≥6*σΔx, then directly exclude the point; compare all data differences with the standard deviation one by one, if a data point satisfies the formula: Δx≥3*σΔx, then mark the point as an outlier; if 5 or more consecutive points are marked as outliers, then mark them as normal data; judge all outliers, if the data value adjacent to the outlier differs by 30% from the outlier, it is considered as an outlier; remove all outliers and perform linear interpolation; if the number of outliers is too large within a 30-minute period, the period is directly excluded.

4. The method of claim 1, wherein, The step 2 combines the peak parameters Q p The calculation method is specifically: ； wherein is the zeroth moment of the wave power spectral density function reflecting the total amount of wave energy; is the wave frequency.

5. The method of claim 1, wherein, The specific calculation formula of step 3 is: ； ； ； ； ； where k is the von Karman constant, is the friction velocity, U z is the average wind speed at height z over the sea surface, z0 is the sea surface roughness, is the stability function, is the dimensionless profile formula, L is the Monin-Obukhov length, is the average air temperature; g is the gravitational acceleration; The overline represents the Reynolds average, and the prime indicates the turbulent fluctuation of the physical quantity.

6. The method of claim 1, wherein, The The specific calculation formula is: ； where a, b, c, d are fitting parameters, reflecting the sea 10-meter wind speed U 10 The basic influence relationship of friction speed, and The fitting parameters are different.

Images