A Scheme for Estimating Time ‐ Varying Wind Stress Drag Coefficient in the Ekman Model with Adjoint Assimilation

: In this study, the time ‐ varying wind stress drag coefficient in the Ekman model was in ‐ verted by the cubic spline interpolation scheme based on the adjoint method. Twin experiments were carried out to investigate the influences of several factors on inversion results, and the conclu ‐ sions were (1) the inverted distributions with the cubic spline interpolation scheme were in good agreement with the prescribed distributions of the wind stress drag coefficients, and the cubic spline interpolation scheme was superior to direct inversion by the model scheme and Cressman interpo ‐ lation scheme; (2) the cubic spline interpolation scheme was more advantageous than the Cressman interpolation scheme even if there is moderate noise in the observations. The cubic spline interpo ‐ lation scheme was further validated in practical experiments where Ekman currents and wind speed derived from mooring data of ocean station Papa were assimilated. The results demonstrated that the variation of the time ‐ varying wind stress drag coefficient with time was similar to that of wind speed with time, and a more accurate inversion result could be obtained by the cubic spline inter ‐ polation scheme employing appropriate independent points. Overall, this study provides a poten ‐ tial way for efficient estimation of time ‐ varying wind stress drag coefficient.


Introduction
Wind stress over the sea surface is a primary driving force of upper ocean circulation and ocean surface waves and is an important mechanism by which the atmosphere conveys momentum to the ocean. Therefore, the exact estimation of wind stress is important for the simulation and forecasting of atmospheric and oceanic processes [1]. Meanwhile, wind stress is a key parameter for understanding physical processes in both the atmosphere and ocean, whose determination is important for coupled ocean-atmosphere modeling [2].
Usually, wind stress is parameterized by the wind stress drag coefficient (WSDC). It has long been recognized that the drag coefficient of the sea surface does not depend on wind speed alone, while it has been highlighted that the wind drag coefficient strongly depends on the sea roughness [3,4]. In previous studies, the WSDC was usually assumed to be a constant [5], considered as a linear function with respect to wind speed [6], or further believed to be a nonlinear function with respect to wind speed [7,8]. At present, it is widely accepted that wind stress coefficients vary nonlinearly with wind speed (and thus with time)-especially at high wind speeds [9,10]. Some scholars have found that spatially and temporally varying WSDC fields are needed in a variety of climate and air-sea interaction studies [11,12]. In particular, with the improvement of the time resolution of the model, it is meaningful to study the WSDC with daily variations. The determination of the WSDC is a difficult problem in the field of marine science and atmospheric science because of the interaction between air and sea; this difficulty stems from the problem of turbulence. The turbulence problem is difficult to solve because of its variability in time and space. Ocean observations at near-inertial frequencies show that the turbulent closure model is better at predicting wind speed than the plate model [13]. There are a variety of options for turbulent closure models, such as models depending on Richardson number, k-ε model, M-Y model, Prandtl mixed-length model, KPP model, etc., which play important roles in different regions and different events [14][15][16][17][18]. However, as the ocean environment often changes with space and time, and the observation data are insufficient, the applications of various turbulence models are easily limited. Specifically, the parameters calculated by different models vary greatly, and few models can be accurately applied to all sea areas or problems in the world.
During recent decades, data assimilation has played an important role in the determination of unknown parameters. The process that the parameters are estimated byassimilating observations with a data assimilation method-is known as parameter estimation. The adjoint method is a powerful tool for parameter estimation, which has been widely used by many scholars. Zhang et al. [19] successfully explored how eddy viscosity is developed with varying wind fields, but did not consider the WSDC. Cao et al. [20] proposed an estimation scheme to evaluate the eddy viscosity profile (EVP) at the bottom of the Ekman boundary layer based on the adjoint method. Zhang et al. [21] successfully inverted the time-varying vertical eddy viscosity coefficients in the Ekman model and also found that the inversion results were much more sensitive to the observations in the upper layers.
In order to improve the computational efficiency and accuracy of parameter estimation, independent point (IP) schemes were used by some scholars in the adjoint method [22][23][24]. The so-called IP scheme is where the WSDCs at some selected IPs are taken as control variables and those at other points are calculated by interpolating values at IPs. There are many interpolation schemes, such as the Cressman interpolation based on linear interpolation. However, the Cressman that was interpolation traditionally used in the IP scheme can result in unsmooth curves. Meanwhile, we found that the cubic spline interpolation (CSI) scheme was often used in practice because of its favorable smoothness and continuity [25][26][27][28]. Therefore, we try to use the CSI instead of linear interpolation to invert the time-varying WSDC in an Ekman model. We aim to improve the simulation capability of the model and provide a reference for the selection of WSDCs in other ocean models. In order to verify the feasibility and accuracy of the CSI scheme, twin experiments and practical experiments are conducted.
The organization of this paper is as follows: In Section 2, the numerical model and CSI scheme are introduced. In Section 3, a series of twin experiments are carried out to verify the feasibility of the CSI scheme. In Section 4, practical experiments are carried out to verify the inverted effect of the CSI scheme. Finally, conclusions are given in Section 5.

Governing Equations (Forward Model)
In this paper, the Ekman model [29] with time-varying WSDC on the sea surface is used, and its governing equations are: where  and  are Lagrange multipliers for u and v , respectively. The problem of minimizing the cost Function (5) is equivalent to finding the stationary point with respect to u , v , d C ,  and  under the condition that the gradients of the Lagrange functional are 0. Then, the following equation can be obtained: 0 0 ( , , , , , ) ( , , , , , ) , .
Equation (7) recover the original governing equations, while the Equation (8) result in the adjoint equations, which are formulated as follows: The corresponding boundary conditions and initial conditions are: The detailed finite-difference scheme for the above equations is given in Appendix A.

Independent Points Scheme and Cubic Spline Interpolation
The basic idea of the IP scheme is simple. For example, the WSDCs are assumed to be varying in time and I = {x1, x2,…, x10} are the time index. First, the indexes of the IPs L = {x1, x4, x7, x10} are a subset of I. For convenience, these IPs are selected uniformly. Second, the WSDCs at all IPs can be inverted by the adjoint method. Finally, the WSDCs at the other points (i.e., x2, x3, x5, x6, x8, and x9) can be interpolated from those at the IPs.
The CSI scheme is used for interpolations in this study and the cubic spline function is employed in this scheme.
Given a series of IPs within the interval [a b], where , and the corresponding function values are 0 1 , , , n y y y  , a cubic polynomial ( ) S x can be constructed in any interval which is expressed as: where i a , i b , i c , and i d are the spline interpolation coefficients, and the above Equation (15) must satisfy the condition of smooth continuity: Meanwhile, Equation (15) must satisfy the boundary condition. In this study, the natural boundary condition is taken, which can be formulated as: The selection and calculation of spline functions are presented in Appendix B. The calculation process of spline interpolation coefficients of ( ) S x can be found in Pan et al. [31].

Estimation of WSDCs with Different Distributions
As mentioned in the introduction, the WSDC and wind speed are time-varying. In addition, the variation of wind speed over time above the ocean is complex, and the specific distribution of the WSDC over time is unknown. Thus, we designed a group of twin experiments to verify the performance of the CSI scheme in the inversion of different prescribed WSDC distributions. Five cases of temporal WSDCs and wind speed are as follows: where T is the same as in the inverse model.
It can be inferred that the prescribed WSDC is increasingly more complex from Case 1 to Case 5. Note that the WSDC of Case 1 is set as 0.00024 in order to keep its variation range the same as that of Cases 2-5, and wind speed only exists in the east-west direction (i.e., 0 = a v ). In order to further highlight the effectiveness of the CSI scheme, experiments with direct inversion by the model (DIM, in which the WSDC at each point is independently inverted) and with Cressman interpolation (CI) schemes are also designed for comparison. The calculation formulas of the three schemes are shown in Appendix B. The number of selected IPs with different inversion schemes for each case is shown in Table 1. Model parameters are set in Table 2. The main purpose of this is to estimate time-varying WSDCs in this study. As mentioned above, the time-varying eddy viscosity has been studied in previous studies, such as Zhang et al. [21]. Therefore, we set the eddy viscosity as a variable that only varies with depth. The classical distribution of the eddy viscosity coefficient  .
where h is depth), which has been widely used in ocean models [32][33][34]. The initial guess value of WSDC is set as 1.2 × 10 −3 in all cases. Considering that the observation data are only available at some depths (5 m and 35 m), two layers (i.e., the first and the seventh layers) of observation data are assimilated in this section. The abbreviations of some variables are shown in Table 3.  The experimental steps are as follows: (1) Run the forward model with prescribed WSDCs, and the calculated velocities are regarded as pseudo "observations".
(2) Run the forward model with the initial guess value of WSDC (Cd0 = 1.2 × 10 −3 ) to obtain the simulated values of velocity.
(3) Run the inverse model with the difference between the "observations" and simulations by the model. The gradients of the cost function with respect to WSDCs can be calculated according to Equation (14). By employing an optimization algorithm, the WSDC can be optimized.
(4) By repeating steps (2) and (3) until the inversion results meet a convergence criterion, the final inversion values of WSDC are obtained. Since the cost functions are almost no longer declining after 500 iterations, we finally selected 500-step iteration for all experiments.
The prescribed and inverted WSDCs are shown in Figure 1. It can be seen that the prescribed distributions of WSDCs (Case 1-5) can be successfully inverted by the CI scheme and the CSI scheme, while they are partially inverted by the DIM scheme. As the distributions become increasingly more complex from Case 3 to Case 5, the consistencies between the inverted distributions obtained by the CI scheme and the prescribed distributions are getting weaker, while the inversions obtained by the CSI scheme are still in good agreement with the prescribed distributions, even at the wave peak and trough. The mean absolute errors (MAEs) and mean relative errors (MREs) between prescribed and inverted WSDCs after assimilation are shown in Table 4. The magnitude orders of MAEs for the DIM scheme were all 10 −4 , while those for the CI scheme and CSI scheme were about 10 −5 -10 −6 . In Case 1, the MRE for the DIM scheme was slightly larger than 25%, and furthermore, the MREs of Case 2-5 were more than 70%, indicating that the results of the DIM scheme were poor. When the distributions were relatively simple (Case 1-3), the MREs of the CI scheme and CSI scheme were less than 6%. Thus, the reliability of the inversion results was high. In addition, the MREs of the CI scheme increased to 17.73% and 20.00% for the two complex distributions (Case 4-5), but the MREs using the CSI scheme remained at a low level (4.48% and 6.34%), which shows that the inversion results of the CSI scheme were more reliable. We further analyzed the curves of the normalized cost function, which are shown in Figure 2. It is shown that the cost functions for the DIM scheme were only decreased by 1-2 orders of magnitude. On the contrary, those with the CI scheme and CSI were scheme generally decreased by 2-4 orders of magnitude. The results indicate that the latter two schemes can achieve better simulation results, demonstrating better inverted WSDCs. Although the cost functions of the CSI scheme decreased more slowly than those of the CI scheme during the several early iterations, they were lower after a certain amount of iterations, indicating that the inversion effect of the CSI scheme was more advantageous. The above results indicate that the CSI scheme showed good inversion ability for WSDCs of different complexities. Here, the simulation ability of current velocity was also assessed. The variations of the MAEs of the horizontal current velocity components (u and v) with the iterations are shown in Figure 3. It can be seen that the MAEs of the CI scheme and the CSI scheme in Cases 1-5 were obviously smaller than those of the DIM scheme. In cases 1-3, the MAEs of the CI scheme decreased more rapidly with increased iteration steps, but the final values of the MAEs were at the same level as those of the CSI scheme. However, for the relatively complex distributions (cases 4-5), the MAEs of the CSI scheme were still smaller than those of the CI scheme. In conclusion, the prescribed distributions of WSDCs can successfully be inverted by the adjoint method with CI or CSI schemes even if only two layers of data are assimilated. As the WSDCs from Case 1 to Case 5 were increasingly more complex, the CSI scheme obtained the best inversion result, and the CI scheme followed. At the same time, the DIM scheme obtained the worst inversion, again indicating that the IP scheme is an effective tool in parameter estimation based on the adjoint method.

Estimation of WSDCs under Different Observation Noise
In all the above cases, the perfect observations without any noise were directly produced by the forward model. However, there was noise in the actual observations of current velocity and wind speed, which might have affected the inversion results. Zhang and Lu [35] indicated that the average difference between the inversion results and true values would increase abnormally, so long as observation errors increased to a certain level. Therefore, we designed another group of experiments to investigate the potential effect of different data noise on WSDC estimation. Considering the distribution complexity, the WSDC in Case 5 was adopted in this group. Gaussian white noise was added to observa-  Table 5). Since group 1 showed that the CI and the CSI schemes were significantly better than the DIM scheme, the former two schemes were used in this group.

Group 2
Case 5-1 3% Case 5-2 5% Case 5-3 10% Case 5-4 15% The prescribed and inverted WSDCs for different observation noise are shown in Figure 4, and the correlation coefficients between them are shown in Table 6. It is clear that the prescribed distributions could be inverted successfully by the CI and CSI schemes in Case 5--1 and Case 5-2. However, in Case 5-3 and Case 5-4, the curves inverted with the CI scheme oscillated more violently than those with the CSI scheme, indicating that the inversion accuracy of the CSI scheme was higher than that of the CI scheme. Combining these results with Figure 4 and Table 6, it could be found that the CSI scheme was more advantageous, especially in Case 5-3 and Case 5-4. Moreover, the distortions of both schemes began to increase significantly with increasing noise level.  Additionally, it can be found in Figure 5 that the normalized cost functions decline greatly with the iterations. Specifically, at most iterations, the cost functions of the CSI scheme are smaller than those of the CI scheme, indicating that the current velocities simulated by the CSI scheme are closer to the observation. This further reveals the superiority of the CSI scheme. Additionally, the MAEs for the horizontal velocity are shown in Figure 6. It can be seen that, after 200 iterations, the MAEs of the CSI scheme in Case 5-1 were notably smaller than those of the CSI scheme, suggesting that the CSI scheme was more advantageous than the CI scheme. To sum up, the prescribed distributions of the WSDCs could still be well inverted by the CSI scheme even there were moderate noises in the observations. With the increase in observation errors, the performance of the CSI scheme became worse. Still, the CSI scheme showed a more stable tolerance of observation noise than the CI scheme.

Practical Experiments and Results
The twin experiments proved that WSDC can be inverted successfully by the adjoint method and that the CSI scheme is an advantageous scheme compared with other two schemes. However, its performance still needs to be testified with real observation data.
To that end, this section conducted practical experiments to invert the WSDCs with the CSI scheme.

Data and Experimental Design
Current velocity and wind speed data were obtained from the ocean mooring station (Papa, 50.1° N,144.9° W) in the mooring system of the Ocean Climate Station Project (OCS) of NOAA, which has been instrumented with upper ocean and surface sensors since June 2007. The current velocity data were measured by an acoustic Doppler current profiler at 5 m and 35 m below the sea surface. The wind speed at an altitude of 4 m measured by the anemorumbometer located at the buoy tower of the sea surface was the average value per hour.
The current velocity and wind speed data of 10 days from 0:00 on 3 July 2007 to 0:00 on 13 July 2007 were selected. As the wind speed was measured at 4 m above the sea surface, it needed to be transformed to the value at 10 m (U10 for short) by a bulk formula [36,37]. The time resolution was transformed from 1 h to 0.5 h by interpolation to meet the demands of the model. The original wind speed (U4 for short) and U10 are shown in Figure  7a. The time series of current velocity are given in Figure 7b,c. Considering the significant influence of the number of IPs on inversion results, the IP numbers were set at 21, 41, and 97 (i.e., the intervals of 12 h, 6 h, and 2.5 h, respectively), and the corresponding numerical experiments were named as Cases 6-8, respectively, as shown in Table 7. The other parameters were set the same as in the twin experiments. Although there was a gap between the zero initial conditions and the real scenarios, the model simulations were gradually adjusted towards the real scenarios by assimilating the observations made.

Results Analysis
The normalized cost function is shown in Figure 8. After assimilation, the cost function in Case 6 decreased by around 30%, and cost functions in Cases 7-8 decreased significantly by about 50%, which indicates that the CSI scheme is effective in reducing the difference between observed and modeled Ekman velocities. In addition, the parameter optimization effect could be improved by increasing the number of IPs, but not significantly if the number of IPs is further increased. The MAEs between the inverted and observed velocities could also demonstrate the effectiveness of the CSI scheme, as shown in Figure 9. It could be found that the MAEs in the three cases all showed a declining trend. The MAEs of U (V) in Cases 7-8 were decreased significantly, from 0.065 (0.07) ms −1 to about 0.05 (0.05) ms −1 , indicating that the distributions of the WSDCs could be inverted by the CSI scheme competently. The time-series of WSDCs and standard wind speed (U10) are shown in Figure 10. It was shown that the temporal variation of the WSDCs was closely related to the variations of wind speed. This consistency becomes more obvious with the increasing IPs. The WSDCs had a peak during 7-8 July 2007, which agreed well with the wind speed. After 11 July, the time-varying curve of wind speed had two relatively deep troughs and one sharp peak, and there were similar variations in the distributions of the WSDCs. But there was a certain phase difference between them. This may have been due to defects in the dynamic process of the model and in the observation data, that was available only at 5 m and 35 m depths. We also compared the time series of observed and inverted current velocities at 5 m and 35 m depth, as shown in Figure 11. For the inverted current velocity at 5 m depth, both the amplitude and phase coincided with observations well, especially in Case 7 and Case 8. However, there were some discrepancies between the modeled results and the observations at 35 m depth. The comparison was also consistent with previous conclusions that the model is more sensitive to upper current velocity. In summary, the above results had a significant downward trend, indicating that the inverted results were stable. As the IP number increased, the temporal variation of inverted WSDC was more similar to that of wind speed U10. This phenomenon is consistent with physical oceanography in the light of there being a positive correlation between WSDC and wind speed.

Conclusions
Based on the adjoint assimilation technique, the CSI scheme was used to invert the time-varying WSDC in the Ekman model. Two groups of twin experiments and one group of practical experiments were carried out. The effectiveness of the CSI scheme was successfully verified by the above experiments. The main conclusions are as follows: (1) The first group of experiments showed that the prescribed distributions of WSDCs were successfully inverted by the CI scheme and the CSI scheme, even if only the observations in the first and seventh layers were assimilated-while the prescribed distributions were partially inverted by the DIM scheme. The CSI scheme was more advantageous than the CI and the DIM schemes as the complexity of WSDC increased. (2) In the second group, even if there was some noise in the observations, the CSI scheme could invert the distribution of WSDCs successfully. With the observation error increasing, the stability of the CSI scheme was higher than that of the CI scheme, and the results of the CSI scheme were smoother. (3) In the practical experiment, the temporal distribution of WSDC inverted by the CSI scheme had a similar variation trend to that of wind speed U10, and a more accurate inverted result could be obtained by the CSI scheme with appropriate IPs.
Overall, this work extends our understanding of the Ekman model and provides an effective method to estimate the time-series of WSDC. However, in the inversion process, some influencing factors are simplified (e.g., the temporal eddy viscosity). Further work, in which the time-varying WSDCs and eddy viscosity are simultaneously estimated, will be the focus.

The DIM scheme
Based on Equation (14), the WSDCs can be optimized as: where  is step-lengths which are used to adjust the parameters preferably. e is the number of iterations.

Cressman interpolation scheme
The scheme selects several IPs evenly in the time interval. The interpolation scheme is as follows: where k d C is WSDC at an unknown time point; j d C represents WSDC at a given time point, which can be calculated through the adjoint method; n is the number of known grid points; and j W is a weighting function, whose general form is: where R is the interpolation radius, and k r is the distance from the unknown time point to the known time point.

Computation of Cubic Spline Interpolation Weights
The cubic spline ( ) S x in the interval [xi, xi+1] between two IPs was constructed as: Natural boundary conditions are selected in this study, namely We can get from Equations (A19) and (A20) that: So, the matrix P is: where m is the weight coefficient.