A method for calculating momentum flux of sea-air interface coupled with sea wave data

Abstract

The application discloses a kind of sea-air interface momentum flux calculation methods of coupling sea wave data, it is related to marine atmospheric boundary layer flux observation technical field, including the threshold value screening of required turbulent flux elements based on the calculation of eddy correlation method, outlier processing is handled, and three-dimensional wind speed is rotated to coordinate;The sea surface undulation data obtained are processed to obtain sea wave spectrum and wave parameter, and based on sea wave spectrum calculation Q p ; Based on the theory of Monin-Obukhov similarity, Monin-Obukhov length, atmospheric stability and 10-meter sea wind U 10 ; Based on the above data, the friction velocity is divided into wind speed term and wave mixing term for fitting, and finally the momentum flux is obtained. The application couples the wave period, wave direction and other key data of sea wave into the calculation system, introduces the Q p peak parameter to represent the sea wave state, fully considers the influence of sea wave on momentum flux, improves the calculation accuracy, and reduces the requirement for observation instrument.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of marine atmospheric boundary layer flux observation, and particularly relates to a sea-air interface momentum flux calculation method coupled with sea wave data. BACKGROUND

[0002] In marine meteorology research, accurate calculation of sea-air boundary layer momentum flux is crucial for understanding the interaction between the ocean and the atmosphere. Accurate acquisition of sea-air boundary layer momentum flux helps to better understand the change rule of the marine climate system, and has important significance for meteorological forecasting, marine environment monitoring and many other fields.

[0003] At present, in the calculation of sea-air boundary layer momentum flux, the commonly used methods mainly include the eddy correlation method and various parameterization schemes. Among them, the eddy correlation method calculates the momentum flux by directly measuring the high-frequency fluctuation of wind speed, which is the most accurate method for obtaining momentum flux data at present. However, the eddy correlation method has many limitations in practical application. First of all, this method involves a large number of observations, not only the high-precision measurement of wind speed components, but also the accurate measurement of parameters such as water vapor density and carbon dioxide concentration, which has very high requirements for observation instruments. Not only is the instrument equipment expensive, but also the accuracy and stability of the instrument are demanding, increasing the observation cost and maintenance difficulty. Secondly, in the actual marine observation environment, especially in complex sea areas, there is no stable construction space, and it is affected by the marine environment such as sea waves and sea currents. The observation conditions of the eddy correlation method are often difficult to meet, which affects the accuracy and reliability of the observation data.

[0004] In addition to the eddy correlation method, the bulk parameter method is also a commonly used method. COARE (Coupled Ocean-Atmosphere Response Experiment) 36 is a widely used flux algorithm for studying the interaction between sea and air. It is based on a series of physical parameterization schemes, and considers various meteorological elements such as sea surface temperature, wind speed, humidity to calculate the momentum flux, heat flux and water vapor flux of the sea-air interface. Compared with the eddy correlation method, it has the advantages of relatively low requirements for observation equipment and is not limited by the environment, which can carry out business observation at more observation sites and provide a large amount of sea-air flux data for global marine areas. However, COARE36 also has some shortcomings. On the one hand, although it considers more meteorological elements, it still does not fully consider the factors such as sea waves when facing complex marine environments. The actual sea waves have a significant impact on the sea-air boundary layer momentum flux. The wave height, wave period, wave direction and other parameters of the sea waves will change the roughness of the sea surface, and then affect the size of the momentum flux, and the calculation result has deviation from the actual situation.

[0005] In addition, some existing simple parameterization schemes usually only consider the wind speed value at 10m above sea level, ignoring many other factors that have important influence on momentum flux. With the deepening of research, it is found through platform data verification that factors such as actual sea wave period have a non-negligible influence on flux. The parameterization scheme relying only on a single wind speed element cannot comprehensively and accurately reflect the real situation of the momentum flux of the sea-air boundary layer, resulting in a large deviation between the calculation result and the actual value, and it is difficult to meet the needs of marine meteorological research and practical application. SUMMARY

[0006] In order to overcome the above problems existing in the prior art, the present application provides a sea-air interface momentum flux calculation method coupled with sea wave data.

[0007] The technical scheme adopted by the present application to solve its technical problems is: a sea-air interface momentum flux calculation method coupled with sea wave data, comprising the following steps:

[0008] Step 1, threshold screening and outlier removal processing are performed on the required turbulent flux elements calculated based on the eddy correlation method, and coordinate rotation is performed on the three-dimensional wind speed;

[0009] Step 2, the obtained sea surface undulation data is processed to obtain sea wave spectrum and wave parameters, and the peak parameter Q of the combined field is calculated based on the sea wave spectrum p ;

[0010] Step 3, the Monin-Obukhov length, atmospheric stability and sea wind speed U at 10m are calculated based on the Monin-Obukhov similarity theory 10 ;

[0011] Step 4, based on the data obtained in steps 1-3, the friction velocity is split into a wind speed term and a wave mixing term for fitting, and finally the momentum flux is obtained;

[0012] The friction velocity in step 4 is The specific calculation formula is:

[0013] ;

[0014] ;

[0015] ;

[0016] Where T p is the peak wave period of the spectrum, is the wind wave angle, F is the wave term function, is the sea wave influence factor, Q p is the peak parameter of the combined field; A, B, C and D are fitting parameters, and the fitting parameters have different values.

[0017] The sea-air interface momentum flux calculation method coupled with sea wave data, wherein the step 1 specifically comprises:

[0018] 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 removal processing on the collected parameter data;

[0019] Step 1.2, performing a first rotation of the three-dimensional wind speed around the z-axis in the x-y plane to make the x-z plane consistent with the average wind direction, and the average component satisfies = 0;

[0020] Step 1.3, rotating the x-z plane obtained in step 1.2 around the y-axis to make the average component = 0, so that the average wind speed u satisfies .

[0021] The sea-air interface momentum flux calculation method coupled with sea wave data, wherein the outlier removal 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, and if the data values adjacent to the outliers differ by 30%, they are considered as outliers; remove all outliers and perform linear interpolation to supplement; if the number of outliers is too large within a 30-minute period, directly exclude the period.

[0022] The sea-air interface momentum flux calculation method coupled with sea wave data, wherein the step 2 peak value parameter Q p The calculation method specifically comprises:

[0023] ;

[0024] wherein, is the zeroth moment of the wave power spectral density function (reflecting the total wave energy); is the wave frequency.

[0025] The sea-air interface momentum flux calculation method coupled with sea wave data, wherein the step 3 specifically comprises the calculation formula:

[0026] ;

[0027] ;

[0028] ;

[0029] ;

[0030] ;

[0031] wherein, kappa is von Karman constant, is friction velocity, U z is the average wind speed at height z above sea surface, z0 is sea surface roughness. is stability function, is dimensionless profile formula, L is Monin-Obukhov length, is the average temperature of air; g is the acceleration of gravity; The upper horizontal line of represents the Reynolds average, and the apostrophe represents the turbulent fluctuation of the physical quantity.

[0032] The above-mentioned sea-air interface momentum flux calculation method coupled with sea wave data, The specific calculation formula is:

[0033] ;

[0034] wherein, a, b, c, d are fitting parameters, reflecting the sea 10m wind speed U 10 The basic influence relationship of friction velocity, and The fitting parameters are different.

[0035] The beneficial effects of the present application are that (1) the calculation precision is improved. The traditional parameterization scheme based on wind speed only ignores key factors such as sea waves, resulting in a large deviation between the calculation results and the actual momentum flux. And the present application fully considers the influence of sea wave elements on momentum flux by coupling sea wave data, and according to the Qp value, different formulas are used to calculate the friction velocity, which can more accurately simulate the real momentum flux of the sea-air boundary layer. Compared with the prior art, the platform data verification shows that the present application significantly improves the accuracy of the momentum flux calculation results, provides a more reliable data basis for marine meteorological research, and helps to better understand the internal mechanism of sea-air interaction.

[0036] (2) The difficulty and cost of observation are reduced. Compared with the eddy correlation method, the required observation amount of the present application is greatly reduced. The eddy correlation method requires high-frequency and accurate measurement of multiple parameters, and the observation instrument is required to be extremely high, the equipment is expensive and the maintenance is complex. The present application only needs to obtain 30 minutes of 10m wind speed, wind direction, wave parameters and sea wave spectrum and other relatively conventional data, greatly reducing the requirements for observation instruments, enabling the calculation of sea-air boundary layer momentum flux to be carried out at a lower cost in more marine areas, especially in complex sea areas like this, improving the practical application feasibility and universality of the scheme.

[0037] (3) The adaptability to complex marine environment is enhanced. The marine environment is complex and changeable, and the sea wave state differs significantly. The present application takes Q p as a segmented discrimination standard, which can flexibly cope with momentum flux calculation under different sea wave conditions. Whether in the case of relatively stable sea waves (Q p <2) or relatively severe sea waves (Q p ≥2), the friction velocity can be accurately calculated through the corresponding formula, and reliable momentum flux results can be obtained, effectively overcoming the problem that the calculation results of some existing schemes are inaccurate when facing complex sea wave conditions. BRIEF DESCRIPTION OF DRAWINGS

[0038] Figure 1 is a schematic diagram of the sea-air interface momentum flux of the present application;

[0039] Figure 2 is the polynomial fitting result of the relationship between the 10m wind speed and the friction velocity at sea in the embodiment of the present application, wherein (a) is the fitting result of Q p <2; (b) is the fitting result of Q p ≥2;

[0040] Figure 3 is the power function fitting result of the relationship between the 10m wind speed and the friction velocity at sea in the embodiment of the present application, wherein (a) is the fitting result of Q p <2; (b) is the fitting result of Q p ≥2;

[0041] Figure 4 is the exponential function fitting result of the relationship between the 10m wind speed and the friction velocity at sea in the embodiment, wherein (a) is the fitting result of Q p <2; (b) is the fitting result of Q p ≥2;

[0042] Figure 5 is a relationship diagram between the vertical velocity spectrum (blue line ≥2, red line <2) and the sea wave spectrum (black line) during January 29, 2024 to February 3 in the embodiment of the present application;

[0043] Figure 6 is the vertical velocity spectrum (blue line ≥ 2, red line < 2) and the sea wave spectrum (black line) during the period from February 17 to February 22, 2024 in an embodiment of the present application;

[0044] Figure 7 is the polynomial fitting result of the function F in an embodiment of the present application, where (a) is the fitting result of Q p < 2; (b) is the fitting result of Q p ≥ 2;

[0045] Figure 8 is the power function fitting result of the function F in an embodiment of the present application, where (a) is the fitting result of Q p < 2; (b) is the fitting result of Q p ≥ 2;

[0046] Figure 9 is the exponential function fitting result of the function F in an embodiment, where (a) is the fitting result of Q p < 2; (b) is the fitting result of Q p ≥ 2. DETAILED DESCRIPTION

[0047] In order for those skilled in the art to better understand the technical solutions of the present application, the present application will be described in detail below in combination with the drawings and specific embodiments.

[0048] The present embodiment discloses a sea-air interface momentum flux calculation method coupled with sea wave data, and a sea-air interface momentum flux schematic diagram is shown in Figure 1 Q p is a parameter for dividing sea conditions. When Q p < 2, it corresponds to a sea wave spectrum case with wide spectrum and flat peak, which means that the sea wave is superimposed by multiple frequency components at this time (such as coexistence of wind wave and swell, sea wave after complex refraction / diffraction near the coast); when Q p ≥ 2, the sea wave in the frequency domain shows narrow spectrum and sharp peak, which means that the sea wave energy is concentrated near a dominant frequency (such as pure swell, single-peak strong wave caused by typhoon).

[0049] The present embodiment calculation method specifically includes:

[0050] Step 1, threshold screening and outlier processing are performed on the required turbulent flux elements based on the calculation of the eddy correlation method, and coordinate rotation is performed on the three-dimensional wind speed.

[0051] This example is based on the data obtained by the observation system carried by Yangjiang platform, with the help of the eddy correlation observation system (including 1 CSAT3A three-dimensional ultrasonic anemometer, EC150 carbon dioxide water vapor analyzer, 1 PTB110 atmospheric pressure sensor and HMP155 air thermometer) and the underwater sea wave observation system (AD2CP) data carried by Yangjiang platform. The three-dimensional ultrasonic anemometer provides 10Hz wind speed three components (u, v, w) and ultrasonic virtual temperature Ts required for flux calculation, the carbon dioxide water vapor analyzer is used to observe the water vapor density and carbon dioxide density in the measured area, the atmospheric pressure sensor provides the atmospheric pressure value P in the measured area, and the air thermometer is used to monitor the atmospheric temperature T in the area. The two observation systems keep the same posture during the experiment, and meteorological and sea wave data of 31 days are obtained.

[0052] For ultrasonic virtual temperature diagnostic value, infrared gas analyzer diagnostic value, water vapor concentration, CO2 concentration, CO2 signal intensity, H2O signal intensity and other parameter data, the numerical values are screened and evaluated, and the data abnormality is checked and the values exceeding the physically reasonable threshold are removed. The principle of removing outliers is as follows:

[0053] 1. Calculate each parameter data type within 30 minutes;

[0054] 2. Calculate the difference Δx of adjacent time series data, and obtain the overall standard deviation (σΔx) of the difference of adjacent points by Δx.

[0055] 3. Compare all data differences with standard deviations one by one, if a data point satisfies the formula: Δx≥6*σΔx, then directly remove the point;

[0056] 4. Compare all data differences with standard deviations one by one again, if a data point satisfies the formula: Δx≥3*σΔx, then mark the point as an outlier;

[0057] 5. If 5 or more consecutive points are marked as outliers, they are marked as normal data;

[0058] 6. Judge all outliers, if the data value adjacent to the outlier is different from the outlier by 30%, it is considered as an outlier;

[0059] 7. Remove all outliers and perform linear interpolation;

[0060] 8. If the number of outliers is too large within 30 minutes, directly remove the period.

[0061] 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.

[0062] The three-dimensional wind speeds (u, v, w) are rotated in 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:

[0063] 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:

[0064] 2. The second rotation is a new rotation of the xz plane around the y-axis, which in turn affects the average components. =0:

[0065] Ultimately, the average wind speed u is satisfied. .

[0066] 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 .

[0067] Wave spectra contain a wealth of wave information, Q p As 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.

[0068] ;

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

[0070] ;

[0071] Qp = 2 is a piecewise standard, which is calculated in the form of a piecewise function in the following wind speed term and wave term calculations.

[0072] Step 3, calculate the Monin-Obukhov length, atmospheric stability, and 10-meter sea wind speed U 10 .

[0073] Based on the Monin-Obukhov Similarity Theory (MOST), the observed wind speed at the height is converted to the 10-meter wind speed U 10 by applying the logarithmic wind profile, which requires atmospheric stability correction. The conversion formula is as follows:

[0074] ;

[0075] ;

[0076] where κ = 0.4 is the von Karman constant, U z is the average wind speed at height z above the sea surface, is the air friction velocity, z0 is the sea surface roughness, is the stability universal function, which is 0 under near-neutral conditions, and its specific expression is as follows:

[0077] ;

[0078] ;

[0079] where L is the Monin-Obukhov length (MO length), and its expression is:

[0080] ;

[0081] where: is the average air temperature; g is the acceleration of gravity, taken as 9.8 m / s 2 ; The upper horizontal line represents Reynolds averaging, and the apostrophe indicates the turbulent fluctuation of the physical quantity.

[0082] Quality control is performed on the above data, including:

[0083] 1. Wind shadow zone data rejection

[0084] Due to the interference of the tower structure from the rear of the instrument (wind shadow zone) during observation by the ultrasonic anemometer carried by the observation system, the data uncertainty increases. According to the installation orientation of the anemometer, the observation data in the 225°-315° direction range is rejected;

[0085] 2. Stability parameter screening

[0086] Based on the Monin-Obukhov similarity theory, unreasonable results (e.g. negative wind speed at low wind speed) may occur when the stability parameter z / L (ratio of observation height to Monin-Obukhov length) is too large. According to the existing research standard, data with z / L>2 and z / L<-2 are removed.

[0087] 3. Low wind speed data screening

[0088] When the wind speed is lower than 1 m / s, the uncertainty of wind direction observation increases significantly, and the relevant data cannot accurately reflect the characteristics of wind stress. Therefore, low wind speed data with U 10 <1 m / s is directly removed.

[0089] 4. Rainfall data screening

[0090] During rainfall, raindrops may pass through the measurement path of the probe, affecting the ultrasonic propagation time and causing errors. Therefore, low wind speed data with P rain ≥0.5 m / s is directly removed.

[0091] 5. 95% confidence interval screening

[0092] The friction velocity is fitted with a second-order polynomial curve of 10-meter wind speed U 10 , and the abnormal data points outside the 95% confidence interval are removed.

[0093] Step 4: Based on the data obtained in steps 1-3, the friction velocity is divided into wind speed term and wave mixing term for fitting, and finally the momentum flux is obtained.

[0094] The momentum flux at the sea-air interface is commonly represented by wind stress, and the calculation formula of wind stress is:

[0095] ;

[0096] ;

[0097] wherein, is the air density, is the drag coefficient.

[0098] The calculation method of this embodiment divides the friction velocity into wind speed term and wave mixing term for fitting, to achieve more accurate and simple momentum flux calculation.

[0099] 1. Wind speed term calculation

[0100] A large number of studies have found that the drag coefficient is mainly controlled by the offshore wind speed, which is usually expressed as a linear function of the wind speed. In this embodiment, the relationship between the 10 m offshore wind speed and the friction velocity is first fitted using the least squares method. The data is selected from the ultrasonic wind speed data based on the ultrasonic anemometer mounted on the Yangjiang platform and the friction velocity calculated based on the eddy correlation method. The ultrasonic wind speed data is converted to 10 m wind speed U 10。 The Levenberg-Marquardt algorithm is used to complete the iterative optimal solution of the parameters in the fitting stage. The fitting form includes polynomial, power function, exponential function, etc. The obtained parameter estimation results meet the requirements of unbiasedness, effectiveness and consistency of the Gauss-Markov theorem. Finally, based on the quantitative evaluation indexes such as the determination coefficient (R 2 ), root mean square error (RMSE), residual distribution rule and 95% confidence interval of parameters, the optimal result-polynomial fitting is selected and the optimal parameter solution is obtained.

[0101] The results of different fitting forms are shown in the following figure and Table 1. In order to determine the optimal fitting relationship between variables, this study constructs three fitting models of polynomial, power function and exponential function, and uses the determination coefficient (R 2 ) and root mean square error (RMSE) as evaluation indexes of model fitting effect. Among them, R 2 characterizes the explanation ability of the model to data variation, and the value tends to 1, indicating that the fitting effect is better; RMSE reflects the average deviation degree of fitted value and measured value, and the smaller the value, the higher the model fitting precision. The fitting evaluation indexes of the three models are shown in Table 1. As shown in the table, the R 2 of the polynomial fitting model is 0.6584 when Q p < 2, and 0.4092 when Q p ≥ 2, which is higher than that of the power function model and the exponential function model, indicating that the model has the strongest explanation ability to the data variation rule. Therefore, this embodiment selects the polynomial fitting model as the representation model of the relationship between variables:

[0102] ;

[0103] When Q p is greater than or equal to 2, a = -0.002802, b = 0.03832, c = -0.118, d = 0.2023; when Q p is less than 2, a = -0.0006677, b = 0.0088, c = 0.02284, d = 0.009237, in order to adapt to the influence characteristics of wind speed on friction velocity under different Q p values. The preliminary friction velocity U The parameters are based on the data of the Yangjiang in-situ observation experiment and theoretical research. The least square method is used as the core optimization target in the fitting process. The Levenberg-Marquardt algorithm is used to achieve the optimal estimation of the parameters of the nonlinear parameter model, which meets the unbiased, effective and consistent requirements of the Gauss-Markov theorem, and reflects the 10-meter wind speed The basic influence relationship of friction velocity.

[0104] Table 1

[0105]

[0106] 2. Wave term calculation

[0107] Introducing the spectral peak factor Q p As a criterion, the samples are divided into two groups, narrow-band swell (Q p ≥ 2, blue curve) and wide-band swell (Q p < 2, red curve), and the results are shown in Figure 5 and Figure 6 .

[0108] First, the spectral peak frequency determines the vertical decay scale of sea wave disturbance, which is the geometric prerequisite for sea wave signal to reach the observation height. Comparing the data of period 1 ( Figure 5 ) and period 2 ( Figure 6 ), it can be seen that although both have high sea 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 on that day, high-frequency short waves are dominant, and the wave boundary layer (WBL) thickness is not enough to cover the sensor height; while the low-frequency swell (about 0.15 Hz) of period 1 decays slowly, successfully transmitting the wave signal to the observation layer.

[0109] In the sub-spectrum graph from February 17 to February 22 ( Figure 6 ), the turbulence energy spectrum of the Q p ≥ 2 sample group has a clear spectral peak at fp; on the contrary, the Q p < 2 sample group of the same period, although it is still dominated by swell, its signal is "smoothed" in the vertical velocity spectrum and cannot form a significant protrusion, and is well integrated with the background turbulence spectrum; during February 17 to February 22, even the Q p ≥ 2 sample group, the intensity of the peak does not change significantly, which may be related to the lower absolute energy level or wider frequency band of the swell during this period.

[0110] Therefore, the wave spectral peak period and Q pAs 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 value is less than 2, A = 0.001136, B = -0.0154, C = 0.04709, D = -0.01002 (fitting results are as follows). Figures 7-9 (as shown)

[0111] ;

[0112] ;

[0113] 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.

[0114] 3. Calculation of final friction velocity

[0115] ;

[0116] 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.

[0117] 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.

[0118] 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.

[0119] 1. Flux data calculated based on eddy correlation method

[0120] The eddy correlation method calculates the turbulent flux by Reynolds averaging, dividing the measured quantity into mean and fluctuation quantities, and then calculating the covariance of the concentration fluctuation of different observation elements and the vertical wind speed component w fluctuation. Therefore, the covariance matrix needs to be calculated by combining each component first, so as to obtain the covariance of each parameter and the vertical component of wind speed w and the mean value of each component.

[0121] (1) Momentum flux Tau (N / m 2 ) is calculated by the formula:

[0122] ;

[0123] Wherein, is the friction velocity (m / s), which can be expressed as:

[0124] ;

[0125] Wherein, is the covariance of the three-dimensional ultrasonic wind u component and the w component, is the covariance of the three-dimensional ultrasonic wind v component and the w component;

[0126] ρ is the air density (kg / m 3 ), which is calculated by the formula:

[0127] ;

[0128] Wherein, P is the air pressure (Pa), T is the air temperature (℃), R d = 287 J / mol / K is the gas constant, e is the saturated water vapor pressure, which is calculated by the formula:

[0129] ;

[0130] Wherein, H2O is the water vapor density, Ts is the ultrasonic virtual temperature (℃);

[0131] (2) Sensible heat flux H (W / m 2 ) can be expressed as:

[0132] ;

[0133] In the formula, is the covariance of the ultrasonic virtual temperature Ts and the three-dimensional ultrasonic wind w component, C p is the constant pressure specific heat, which can be calculated by the formula:

[0134] ;

[0135] In the formula, C pd= 1004.67 (j / kg / K) is the specific heat at constant pressure of dry air, q is the specific humidity, calculated by the following formula:

[0136] .

[0137] (3) Ultrasonic virtual temperature correction:

[0138] The initial sensible heat flux is obtained by calculating the covariance of ultrasonic virtual temperature fluctuation and vertical wind component , but there is a certain deviation between ultrasonic virtual temperature and actual temperature, which needs to be corrected by to obtain the final sensible heat flux H:

[0139] ;

[0140] .

[0141] 2. Flux data calculated based on COARE36

[0142] COARE (Coupled Ocean-Atmosphere Response Experiment) 36 is a widely used flux algorithm for studying the interaction between the ocean and the atmosphere. It is based on a series of physical parameterization schemes, which comprehensively considers various meteorological elements such as sea surface temperature, wind speed, humidity, etc. Through complex formulas and algorithms, it calculates the momentum flux, heat flux and water vapor flux at the sea surface. To calculate the flux data using COARE36, the following input data is required: sea surface temperature, wind speed, wind direction, atmospheric temperature and humidity, atmospheric pressure P, etc.

[0143] The flux calculated based on the eddy correlation method is taken as the true value, and compared with COARE36 and the parameterization scheme respectively. The results from January 10 to February 11, 2024 (33 days) in the experimental data are analyzed. The evaluation indicators of the calculation results include the determination coefficient (R2), the mean absolute error (MAE) and the root mean square error (RMSE).

[0144] The determination coefficient (R 2 ) 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 explain the data, i.e. the higher the correlation between the calculated results and the true value.

[0145] ;

[0146] where n is the total number of data, is the observed value (here referring to the true value of the flux calculated based on the eddy correlation method), is the predicted value (the flux value recalculated by the parameterization scheme or COARE36), The mean of the observed values. Table 2 shows the determination coefficient (R 2 ) of the friction velocity calculated by the COARE36 and the method of the present embodiment and the observed values, from which it can be seen that the determination coefficient (R 2 ) of the method of the present embodiment is higher than that of the COARE36, which shows that the correlation of the calculated structure of the method of the present embodiment with the true value is higher.

[0147] Table 2

[0148]

[0149] The mean absolute error (MAE) reflects the average error between the predicted value and the true value, and the smaller the value, the closer the model prediction result is to the true value.

[0150] ;

[0151] where n is the total number of data, is the observed value (here, the true value of the flux calculated based on the eddy correlation method), is the predicted value (the flux value recalculated by the parameterization scheme or the COARE36). Table 3 shows the mean absolute error (MAE) of the friction velocity calculated by the COARE36 and the method of the present embodiment and the observed values, from which it can be seen that the mean absolute error (MAE) of the method of the present embodiment is smaller than that of the COARE36, which shows that the prediction result of the method of the present embodiment is closer to the true value.

[0152] Table 3

[0153]

[0154] 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, and the smaller the value, the higher the stability and accuracy of the model prediction.

[0155] ;

[0156] where n is the total number of data, is the observed value (here, the true value of the flux calculated based on the eddy correlation method), is the predicted value (the flux value recalculated by the method of the present embodiment or the COARE36). Table 4 shows the root mean square error (RMSE) of the friction velocity calculated by the COARE36 and the method of the present embodiment and the observed values, from which it can be seen that the root mean square error (RMSE) of the method of the present embodiment is smaller than that of the COARE36, which shows that the prediction stability and accuracy of the method of the present embodiment are higher.

[0157] Table 4

[0158]

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

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. 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; 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.

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 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 upper horizontal line represents the Reynolds average, and the apostrophe indicates the turbulent fluctuation of the physical quantity.

Images