Modeling for the performance of navigation, control and data post-processing of underwater gliders

Underwater gliders allow efficient monitoring in oceanography. In contrast to buoys, which log oceanographic data at individual depths at only one location, gliders can log data over a period of up to one year by following predetermined routes. In addition to the logged data from the available sensors, usually a conductivity-temperature-depth (CTD) sensor, the depth-average velocity can also be estimated using the horizontal glider velocity and the GPS update in a dead-reckoning algorithm. The horizontal velocity is also used for navigation or planning a long-term glider mission. This paper presents an investigation to determine the horizontal glider velocity as accurately as possible. For this, Slocum glider flight models used in practice will be presented and compared. A glider model for a steady-state gliding motion based on this analysis is described in detail. The approach for estimating the individual model parameters using nonlinear regression will be presented. In this context, a robust method to accurately detect the angle of attack is presented and the requirements of the logged vehicle data for statistically verified model parameters are discussed. The approaches are verified using logged data from glider missions in the Indian Ocean from 2016 to 2018. It is shown that a good match between the logged and the modeled data requires a time-varying model, where the model parameters change with respect to time. A reason for the changes is biofouling, where organisms settle and grow on the glider. The proposed method for deciphering an accurate horizontal glider velocity could serve to improve the dead-reckoning algorithm used by the glider for calculating depth-average velocity and for understanding its errors. The depthaverage velocity is used to compare ocean current models from CMEMS and HYCOM with the glider logged data.


Introduction
Today, underwater gliders are an inherent part of monitoring oceans. These platforms have proven their efficiency and robustness in the collection of oceanographic data in the last two decades [1][2][3]. The first operational underwater gliders, called "legacy gliders", were the Seaglider [4] built by the University of Washington, the Spray [5] built at the Scripps Institution of Oceanography, and the Slocum [6] developed by the Webb Research Corporation. This paper focuses on the Slocum glider. It should be noted that the equations and methods presented are also applicable to other glider types. The data used are from the Center for Ocean Observing Leadership (COOL) at Rutgers University [7]. The Rutgers glider team started with the first Slocum glider missions in 2003 [8]. Since that time, the team has conducted 505 missions, which have mapped ocean properties over 252,944 km during 13,563 days at sea. The team is also involved in the Challenger configuration [11], the fitted sensors [12] and biofouling [13] also have an important influence on the glider model and should be considered.
The angle of attack is of crucial importance in glider modeling. It will be used for the calculation of the horizontal and vertical glider velocity. The horizontal glider velocity is required for dead-reckoning navigation during a mission [14] and for estimating the depth-average velocity [15]. The depth-average velocity results from the difference between the glider velocity over ground and the horizontal glider velocity through water. The horizontal velocity is also applied in the planning for long-term glider missions to find an optimal/passable path from a defined start to a goal and to get information about the feasibility of the mission with regards to energy consumption and the estimated arrival time at the goal [16]. The difference between the vertical glider velocity and the depth rate was used in [17] to determine the vertical current velocity.
Section 2 presents the glider flight model with relations and dependencies of the angle of attack and the horizontal glider velocity, which are based on the analysis of Slocum glider models presented in Section 3. The order of the sections (first: glider model used, second: relevant glider models) relates to the fact that many issues presented in Section 3 are explained in Section 2. Approaches to estimate the model parameters and to determine the angle of attack will be described in Section 4. Section 5 presents the results of the parameter identification for various parameter sets. Requirements regarding the logged data in order to determine trustworthy model parameters will be discussed. Horizontal glider velocities from a long-term mission are used in Section 6 to calculate the depth-average velocities which are used for the comparison of ocean current models from CMEMS [18] and HYCOM [19] with the glider data.

Glider flight model
A generally used glider model for a steady-state gliding motion will be described in the following sections. Fig. 1 shows a schematic illustration of a glider with the defined reference frame, angles, velocities and forces.

Calculation of the horizontal glider velocity
It should be noted that the angle relations in Fig. 1 correspond to real glider conditions where a glider has a positive angle of attack α during the dive phase and a negative angle of attack during the climb phase.
The navigation of a glider, the planning of a long term glider mission or the in-situ estimation of the ocean current conditions during a mission require information about the correct horizontal glider velocity v x . This velocity is dependent on the vertical glider velocity v z and the glide path angle ξ, which is the result of the relation between the pitch angle θ and the angle of attack α = (1) The angle of attack α is defined as the angle between the projection of the total velocity vector of the glider = u v w V [ , , ] onto a vertical plane, formed by the body-fixed x b and z b axes (x z b b plane), and the body-fixed x b axis. The angle α can be defined with the body-fixed velocities u and w as = w u tan 1 (2) Assuming that the glider has zero roll and no yaw moment, this vertical plane can be used for a simplified glider flight model when the glider moves only in this plane. This assumption and the requirement for a symmetrical glider body will form the basis for the following modeling steps. This way, the resulting glide speed V can be calculated with the two body-fixed velocities u and w as which can also be described with the horizontal and vertical glider velocity through water v x , v z as and as function of the glide path angle ξ The horizontal glider velocity v x , which is of interest, can be computed as The vertical glider velocity (through water) v z results from the difference between the depth rate (vertical velocity over ground) z and the vertical current velocity v , zcurrent which is generally assumed to be zero In contrast to the variables vertical glider velocity v z and pitch angle θ, which can be directly derived from the logged glider data during the steady-state gliding, the angle of attack α has to be detected using additional glider parameters. The necessary steps and relationships for this are described below.

Force-velocity relations
The horizontal and vertical force equations for the glider in the vertical plane at equilibrium steady glides are where F D is the drag force, F L is the lift force, ξ is the glide path angle and F Bnet is the net buoyancy force of the glider given by . For a neutral buoyancy trimmed glider, m 0 corresponds to the variable ballast mass m b . For safety reasons, a glider is often trimmed slightly light, so that it floats when the buoyancy engine is set to . Therefore, a buoyancy trim offset Δm 0 has to be added to the resulting excess mass: Re-arranging Eqs. (8) and (9) into a separate description for F D and F L gives the buoyancy force components The hydrodynamic forces are modeled as where C D and C L are the non-dimensional drag and the lift coefficients which refer to the reference area A, ρ is the density of water, and V is the glide speed.

Angle of Attack relations
The coefficients C D and C L are functions of the angle of attack α and will be generally modeled as Using the horizontal force Eq. (8) and the hydrodynamic force Eqs. (13) and (14), the glide path angle ξ can be calculated as The necessary steps to solve α numerically will be presented in Section 4.2.
A characteristic value to describe the gliding flight is the lift-to-drag ratio L/D which can be described using Eqs. (13)-(16) as There is a maximum in L/D [20] at The smallest glide angle ξ min is at (L/D) max and can be calculated as where cot 1 is the arccotangent. Using Eq. (11) in (13) or (12) in (14), the glide speed V can be described as function of the related buoyancy force component to the hydrodynamic force as The vertical glider velocity v z can be calculated by substituting Eq.
Eq. (24) allows the calculation of the angle of attack for minimum vertical glider velocity v ,min z by maximizing the term C L (α) 3 /C D (α) 2 in the simplified equation, where the cosine term is ignored [21]. For typical, small glide angles ξ is the term cos(ξ) 3 close to one and constant. The solution of the extreme value problem In analogy to the approach described above it is possible to determine the glide path angle where the horizontal glider velocity v x has its maximum. Substitute Eq. (22) into (6) and maximize the term sin (ξ)cos (ξ) 2 in the simplified Eq. (28), where the drag coefficient C D (α) for small angles of attack is nearly constant. The extreme value problem for has one solution at The value x is valid for all types of gliders and independent of hydrodynamic coefficients [22]. Fig. 2 shows a graphical representation of the glider velocities v x and x z for dive. This form of representation corresponds to the glide polar of a sailplane [23] and can be used for graphical analysis and to determine the glider velocities v x and v z and angles α, ξ and θ for the operation point of interest. This plot is mirrored on the abscissa to the commonly used glide polar representation for a glider [22,24,25]. A glider should be flown within a glide angle range between ξ min and ξ x,max and between the operation points B and C. Here the slowest glide slope is defined at the maximum lift-to-drag ratio and results in the minimum specific energy consumption of the glider [22]. This slope is called stall glide slope in [25]. Flying at higher angles of attack will enable the lowest sink rate of a glider v z, min in operation point A. Because the linear relation of the lift curve C L (α) is not guaranteed for high angles of attack, this operation range should be avoided. The lift curve is dependent on the structural design of the glider and is mainly influenced by the shape of the wings. High angles of attack result in boundary layer separation and the wings stall [26].

Relevant work about slocum glider models
This section presents relevant work about the modeling of Slocum gliders. The interested reader will find information about model design, strategies to find the model parameters and relations between the individual models.

Relevant work
Vehicle Control Technologies, Inc. (VCT) [27] used computational fluid dynamics (CFD) computations to determine the coefficients C D and C L of the Slocum glider. These results are published in [24] and [26]. In this work, the coefficients C D and C L are normalized by the square length of the hull of 1.789 m. The reference area A is thus (1.789 m) 2 = 3.2 m 2 .
Graver [28] directly used the lift coefficient from [27] with a rescaling in his reference area, i.e. the frontal area of the vehicle for his glider flight model. With a glider diameter of 0.2127 m the reference area is A = 0.0355 m 2 . The drag coefficient was determined using the glider logging data and the given lift coefficient. An estimation of the buoyancy trim offset Δm 0 was necessary to avoid an asymmetrical drag coefficient curve.
Bhatta [29] used hydrodynamic force equations by analogy with Eqs. (37) and (38) so that the reference area definition could be omitted. Some parameters were estimated using glider logging data from sea trials in [28] and wind tunnel experiments in [26].
Williams et al. [30] used an iterative scheme to obtain estimates for the parameters C , D0 C D2 and C L1 for the lift and drag coefficients. The individual parameters were determined sequentially in a loop using the logged glider data: pitch angle θ, vertical glider velocity v z and the net buoyancy F B net . This loop is stopped when all parameters converge to a stable set. Here the reference area A is the frontal area of the glider. The values of the estimated parameters during dive are different from the values during climb. According to [30] this is due to the ballasting procedure. Williams's work also includes a detailed hydrodynamical analysis of a hull with a length-to-diameter ratio similar to a Slocum glider tested in a towing tank.
Merckelbach et al. [17] consider the hull and wings of a Slocum glider separately for the parameters of the lift and drag coefficients. The individual parameters for hull and wings are added together to give the total coefficients, so that they can be used in Eqs. (15) and (16). The buoyancy force calculation used considered the additional influences of water pressure P and water temperature T on the volume of the displaced fluid where V g is the glider volume at atmospheric pressure, ε is the compressibility of the hull, a T is the thermal expansion coefficient, and ΔV bp is the volume change resulting from the buoyancy engine. Therefore, the parameters V g and ε are estimated using the logged glider data. An additional estimated parameter is C D0 in Eq. (15). All other parameters in Eqs. (15) and (16) are the result of experiments or empirical formulas as described below. The determined parameters C D2 and C L1 from hull tow tests in Williams et al. [30] were used for the hull segment. Since these parameters refer to the frontal area A F , rescaling is applied. The reference area in this work is the wing surface area A W = 0.1 m 2 . The rescale factor is therefore A F /A W = 0.038/0.1 = 0.38. The parameter C L1 for the lift coefficient of the wings is the result of a semiempirical formula [31] for a lift-curve slope using the aspect ratio , and the wing sweep angle Λ. (A detailed analysis and comparison of possible formulas for a lift-curve slope for a Slocum glider is given in [32].) Prandtl's lifting line theory was used for modeling the drag coefficient of the wings. In this theory the drag coefficient C D has two terms, a parasitic drag C D0 as a result of the form drag or pressure drag, the skin friction drag and the interference drag and an inducted drag C Di calculated using Prandtl's lifting-line theory where (33) The aspect ratio is calculated as b 2 /A W , where b is the wing span and A W is the wing area. The span efficiency parameter e in Eq. (33) allows the consideration of a real lift distribution, which is usually disturbed through the addition of fuselage, engine nacelles or other parts [33]. It should be noted that some of the relations used to describe the hydrodynamic vehicle behavior are based on aerodynamic studies where there is a wealth of experience. Another name for this factor is the Oswald efficiency factor. (A wide range of approaches to calculate this factor are presented in [34].) Mahmoudian [35] also used Eq. (32) to describe the drag coefficient of the entire glider. The parameters used for the lift and drag coefficient in this model are from Bhatta [29].
Cooney [14] used exclusively empirical formulas from [33] to describe the glider behavior. Analogous to Merckelbach [17], a separate consideration of the hull and wings of a Slocum glider is conducted. The calculation of the drag coefficient is based on Eq. (32) where the parasitic drag C D0 is exclusively the result of the pressure drag. The calculation of the pressure drag is based on Hoerner [33] (3-12 equation (25)) where D is the hull diameter, L is the characteristic length and C (Re) f is the skin-friction drag coefficient for a flat plate in laminar flow which is a function of Reynolds number Re The Reynolds number Re represents the ratio of the dynamic forces relative to the friction forces [33] as where V is the glide speed, L is the characteristic length and ν the fluid kinematic viscosity, which is taken to be 1.35× 10 6 m 2 /s at 10 ∘ C. The speed-dependent behavior of the drag coefficient is a unique feature compared to the other models presented. Fig. 3 shows the pressure drag C D0 as a function of glide speed V. Thus, a lower glide speed V leads to a greater drag coefficient C D . The lift coefficient C L due to the wings is determined using the lifting-line theory of swept wings. The parameters of the empirical formula are the aspect ratio and the wing sweep angle Λ.

Comparison of the models
Although the models presented above describe the lift and drag coefficients for a Slocum glider, it is difficult to compare the individual parameters. The reasons for this are the different definition for the reference area A in Eqs. (13) and (14), different analytical equations to describe the drag and the lift coefficients, different types of Slocum glider and wing configurations. To compare the parameters used in the individual models, the following hydrodynamic force equations with the dimensioned lift and drag coefficients K L and K D will be used where p fw is the density of freshwater 1000 kg/m 3 . The relations of the dimensioned parameters K , L1 K D0 and K D2 to the non-dimensioned parameters C , L1 C D0 and C D2 in Eqs. (13) and (14) are located below the underbraces. This form of equation will also be used in this paper. Table 1  It should be noted that the lift coefficient in VCT [27] and Graver [28] is defined as which does not allow a direct comparison with the other models. In these cases, a linear regression can be used to estimate a linear coefficient C L1 for Eq. (16). Therefore, the data points for α will be created ranging from 0 ∘ to 10 ∘ , which form the elements of the regression matrix . The lift coefficient values for Eq. (39) can now be calculated using these data points. These calculated values form the output variable y in the regression. This results in a parameter K L1 equal to 219.93 ± 0.17 kg/m/rad for VCT and 219.43 ± 0.17 kg/m/rad for Graver.
To use a fixed drag coefficient for Cooney's work [14], two operating points, a glide speed = V 0.3 and = V 0.7 m/s, were chosen in Eq. (34). These give two K D0 values of 6.06 and 5 kg/m which allow a direct comparison with the other models. Table 1 shows a wide dispersion of the parameters and characteristic values. This also is evident in Figs. 4 and 5, which show the drag and lift coefficients as a function of α for the models. The curves of liftto-drag ratio in Fig. 6 as well as the curves of glide path angle in Fig. 7 for dive, where the angle of attack α is positive show four characteristic groups (The curves for negative angles of attack are mirrored in the diagonal quadrant.). The corresponding group number is shown to the right of the curves.
The first group includes the work from VCT [27], Bhatta [29], Merckelbach [17] and Mahmoudian [35] where the maximum in L/D lies between 5.7 and 7.2 for α between 9.4 ∘ and 12.2 ∘ . The second group contains the results from Cooney [14] where the maximum in L/D lies between 3.5 and 3.8 for α between 10.3 ∘ and 11.3 ∘ . A third group contains the results of Williams [30]. The L/D curves for the dive and climb are similar, but the maximum in L/D lies at 1.8, which equals approximately one-third of the value of the first group. A reason for this could be in the initial configuration for the iterative scheme, which includes parameters only of the glider hull detected in the towing tank tests. The α values are again similar to the first group. The last group includes Graver's model. The maximum in L/D is 2.35 at an α of 4.7 ∘ . There is also a significant difference in the drag coefficient curve in Fig. 4 compared to the other models. Possible explanations for the high drag coefficient used by Graver [28] are complex geometry through wing deformation and a CTD sensor in real or bad data logging conditions due to a sideslip angle caused by a static roll. Fig. 8 shows the determined angle of attack α and the resulting horizontal glider velocity v x using mean values for the pitch angle θ and the vertical glider velocity v z in Table 2 for all dives and climbs in the  Table 1 Parameters and characteristic values for the individual models.
VCT (2003)   period from 2018-09-05 to 2018-09-18 for the glider mission RU29-550 from Sri Lanka to Mauritius. As Fig. 8 shows, there is a maximal difference of 8% in the horizontal glider velocity in all models, except Williams. Although the models compared here have different parameters and L/D curves, the calculated horizontal vehicle velocities are close. The maximal difference in the horizontal glider velocity between Williams and the other models is around 20%.

Background
This section describes an approach to estimating the parameters for the lift and drag coefficients to determine α which is used in Eq. (6) for calculating the horizontal glider velocity v x . The idea is to minimize the difference between the logged vertical glider velocity v z , derived from the depth z, and the modeled vertical glider velocity v z using estimated parameters. The vertical force Eq. (9) is used for modeling the vertical glider velocity v z . Combining Eqs. (5), (37) and (38) into (9) gives Solving The estimated parameters are marked in bold and summarized in the parameter vector β In addition to the parameters for the lift and drag coefficients, the buoyancy trim offset Δm 0 is a further parameter to be estimated. This is necessary since the lift and drag coefficients are used for the modeled vertical glider velocity in both the dive and the climb. Without calibration of the buoyancy trim offset Δm 0 , the modeled lift and drag coefficients for the individual dive and climb lie above or below the ideal lines. The parameter vector β can be estimated using a search method to minimize the cost function C given by     Table 2.
which corresponds to the mean squared error (MSE) where n is the number of all logged dives and climbs within a defined time period. The values ρ, θ and m b used in Eq. (41) are average values computed for each dive and climb.

Determination of the angle of attack
The calculation of the lift and drag coefficient in Eq. (41) requires a known angle of attack α. This angle can be solved numerically by inserting Eqs. (37) and (38) into Eq. (17) to give A bisection or bracketing method can be used to solve this onedimensional optimization task. A bracketing method works without derivatives and finds the minimum through iterative decreasing of the interval until the desired tolerance ϵ is achieved, where the minimum lies. In this approach the bisection method in MATLAB fzero and the bracketing methods Golden section search [36], Fibonacci search [37] and Brent's method [38] were tested. Brent's method was used as it requires lesser cost function calls compared to the other methods. On average, 7-9 cost function calls are needed to determine the angle of attack. The MATLAB function name for this bracketing method is fminbnd. This method requires a cost function f(x) and a fixed interval [xl, xu] wherein the search parameter x lies. The cost function is defined as Good interval values can be found by using a quadratic approximation of Eq. (44) to determine α where the function tan(x) is approximated with x The resulting quadratic equation is therefore The approximated angle of attack α can be solved using the discriminant disc. An evaluation of the real roots gives In cases that The resulting interval is calculated using α from Eq. (48) or (49) by where Δα was chosen as 3 ∘ . This value worked well for all examined glider missions. The cost function C(β) in Eq. (43) thus includes an internal search procedure to detect the corresponding angle of attack α for every dive or climb.

Initial parameters setting
To guarantee a good convergence and to find the global minimum, initial parameters have to be defined. The parameters for the lift and drag coefficients K , L1 K , D0 and K D2 from Cooney [14] were used as initial parameters in this paper (The parameter K D0 corresponds to the mean value of this parameter in Table 1). The initial parameter for the buoyancy trim offset Δm 0 can be calculated by inserting Eq. (11) into Eq. (37) with the assumption of similar drag coefficients for all dives and climbs as follows where the unknown parameter angle of attack α in ξ was set to zero and mean values of all dives and climbs were used.

Results
Logged glider data from the Challenger Glider Mission [9,10] was used to evaluate the parameter identification approach presented in Section 4. The goal is the accurate determination of the horizontal glider velocity v x during a mission using the identified lift and drag coefficients K , L1 K D0 and K D2 . The glider data used in Eq. (41) is the mean values of the water density ρ, the pitch angle θ and the variable ballast mass m b computed for each dive and climb in the defined time period. The nonlinear regression function nlinfit from MATLAB was used to estimate the coefficients. This function uses the Levenberg-Marquardt nonlinear least squares algorithm and allows an easy determination of the confidence intervals for the coefficients using the MATLAB function nlparci after the minimization. The independent variable matrix X, the parameter vector β and the output vector y can be written as

Minimizations for different parameter sets
The first tests used data from the glider mission RU29-550 from Sri Lanka to Mauritius for a time period of two weeks from 2018-09-05 to 2018-09-18. The data of the glider state during the time period is shown in Fig. 9.
Five minimizations with different Parameter Sets PS* to be minimized were executed. The placeholder character * includes the positions of the estimated parameters in the parameter vector in Eq. (54). For example -PS24 means only parameter K D0 and Δm 0 will be estimated, parameter K L1 and K D2 have fixed values and will not be estimated. The estimated parameters and their 95% confidence intervals are shown as black text, while the unestimated parameters are shown in gray text in Table 3. An underlined parameter value represents a statistically insignificant value, where the corresponding confidence interval includes zero.
The unestimated parameters correspond to the parameters from Cooney [14]. The buoyancy trim offset Δm 0 was an estimated parameter in every minimization. The results for all minimizations is similar around −120 g (positively buoyant). This means that the glider is trimmed slightly light, which makes sense for safety reasons. The resulting drag and lift coefficients as a function of α for the individual minimizations are shown in Figs. 10 and 11.
All minimizations found solutions where the modeled vertical glider velocity v z agrees well with the logged vertical glider velocity v z , as shown in Fig. 12.
Likewise, the calculated drag and lift coefficients K D and K L , found using the logged glider data by substituting Eqs. (5) and (11) and by substituting Eqs. (5) and (12) into Eq. (38) Fig. 9. Glider data used for the minimization: pitch angle θ, vertical glider velocity v z , ballast mass m b , water density ρ.    Although not all possible parameters were estimated in the first four minimizations, the good correlation between the logged and modeled data in Fig. 12 is remarkable. In addition to a good agreement between the chosen unestimated parameters and the real system parameters is another reason for the good correlation the cluster-like distribution of the logged data. This can be seen in Figs (57). In the case that the parameters of only one coefficient need to be optimized, the parameters result from the intersection of the given coefficient curve and the curve to describe possible settings of α for the operating point. This can be seen for the first minimization PS14 in Fig. 15, where the intersection of the given drag coefficient curve and the curve to describe possible settings of α results in an angle of attack α of around 1.46 ∘ . The three minimizations PS24, PS34 and PS234 use a given lift coefficient and estimate one or both parameters of the drag coefficient. This corresponds to Graver's work, where the lift coefficient was used from VCT [27] and only the drag coefficient was determined. Although the values found for K D0 and K D2 in Table 3 are very different, the calculated angles of attack α and the horizontal glider velocities v x in Fig. 17 are similar. The reason for this is the cluster-like distribution of the logged data, which leads to an operating point of α = 3.4 ∘ and K D = 6.8 N(s/m) 2 using the intersection method for the lift coefficient in Fig. 16 described      Table 2. above. This operating point can be modeled by various quadratic function curves, which is clearly shown in Fig. 10, where this operating point corresponds to the intersection point of the curves PS24, PS34 and PS234 (red circles). It should be noted that the negative value for the estimated parameter K D2 obtained for minimization PS234 is a result of the cluster-like distribution of the logged data. A negative value for the parameter is impossible from the physical point of view, but it leads to the minimal cost for C(β) in Eq. (43) using the logged data.
The minimization PS1234 has no limitations regarding the lift or drag coefficient curve used. This means that the lift and drag coefficients used to model the operating point are the result of a defined angle of attack α on the possible setting curves in Figs. 15 and 16. Multiple combinations for K D and K L are thus possible. The parameter set found leads to an angle of attack of around 6.8 ∘ for dives. This value is twice as high as the results of Cooney's approach (see Fig. 8) or the other minimizations (see Fig. 17) and is not trustworthy when using data with a cluster-like distribution around an operating point.
This data distribution is shown in Figs. 18 and 19 where the data samples are greatly scattered around the ideal curves.
The shape of these data distributions can be described by the scale parameters interdecile range (IDR) between the 10th and 90th percentiles for the pitch angle θ and the residual standard deviation of a linear regression s y.x for the vertical velocity v z . A linear regression model = + v z 0 1 is admissible for a pitch angle between ± 15 ∘ and ± 35 ∘ where the ideal curves in Figs. 18 and 19 can be assumed as linear. The using of IDR instead of a standard deviation results from the fact that the pitch angle data is not normally distributed and multimodal in practice. The MATLAB call for IDR is IDR=diff(prctile(x,[10 90])). The residual standard deviation is defined as where n is the sample size, x i and y i are the individual data samples, x and ȳ are the sample means and s x and s y are the standard deviations.
The variables x and y correspond here to the pitch angle θ and the vertical glider velocity v z . A linear relation between the vertical glider velocity v z and the pitch angle θ leads to a correlation coefficient of -1. The calculated coefficients for the dives and climbs using the logged data are -0.77 and -0.802. A reason for these distributions could be an existing vertical current velocity v zcurrent in Eq. (7) caused by internal waves. This is also discussed in [15] as a possible influencing factor on the measurement accuracy of depth-average velocity.

Sensitivity and topographical analysis of the cost function
For a better understanding of the influence of the parameters β i on the cost function C(β) in Eq. (43), and thus on modeling the vertical glider velocity v z in Eq. (41), a global sensitivity analysis is carried out. The method used is based on cumulative distribution functions (CDFs) and is described in detail in [40]. The key idea is to analyze the differences between the conditional and unconditional CDFs of the output y using the Kolmogorov-Smirnov (KS) statistic as a measure of their distances. This distance correlates to the sensitivity of the input x i to the output y. The unconditional CDF is the result when all inputs vary simultaneously in defined ranges, whereas the conditional CDFs result when varying all inputs except x i , which is fixed at a defined value. The KS statistic provides a curve of sensitivity over all defined fixed values of x i . A new density-based sensitivity index, called PAWN was presented in [40], where a defined statistic such as the maximum, mean or median extracts a single value T i from the KS curve for every input x i . This value varies between 0 and 1. A low value of T i implies a smaller influence of x i on y.
In this paper, the inputs x i of the sensitivity analysis are the parameters β i and the output y is the cost function value C(β). The data used corresponds to that of the previous section. The Sensitivity Analysis for Everybody (SAFE) Toolbox for MATLAB [41] was used for the analysis. The curves for the KS statistic in Fig. 20 show a high sensitivity for the parameters K D0 and Δm 0 . The parameters K L1 and K D2 have low sensitivity within the whole range and lie under the critical level for a confidence level α of 0.05. This can explain the large confidence intervals for these parameters in Table 3 and their uncertain estimation in Section 5.1. An interesting point is the increase in the middle of the KS curves of the parameters K D0 and Δm 0 in Fig. 20 that correspond to the estimated parameters. These parameter areas show a higher sensitivity  to the cost function because they characterize the area of the cost function minimum. Fig. 21 shows the resulting PAWN sensitivity indices for the four parameters and their confidence intervals. It shows clearly that the most influential parameters of the cost function C(β) are K D0 and Δm 0 . The analysis of the cost function topography may also shed light on the influence of the parameters β i on the cost function C(β). To additionally examine the influence of the data distribution on the cost function topography, simulated data was generated. The parameter set of minimization PS24 in Table 3 was used to calculate the vertical glider velocities v z in Eq. (41) and define thus the minimum of the cost function. Two scenarios were examined. The first scenario includes one operating point for dive and climb where the pitch angles θ, the variable ballast masses m b and the water density ρ in Table 2 were used for the calculations in Eq. (41). The second scenario includes an additional operating point where the pitch angles θ are 7 ∘ smaller than in the first scenario. Fig. 22 shows the contour plot for the cost function C(β) with respect to the parameters K L1 and K D2 when the parameters K D0 and Δm 0 are fixed at their optimal values for the two operating points scenario. The minimum lies inside a long flat valley. Using the cost function for one operating point results in a similar valley without significant minimumarea. This can clearly be seen in Fig. 23, which shows the parameter K L1 and the minimum of the valley C min with   respect to the parameter K D2 for one and two operating points. The minimum of the valley C min for two operating points shows a significant minimumarea, whereas the minimum curve C min of one operating point is very flat, which makes it difficult or impossible to converge to the minimum. Multiple combinations of K L1 and K D2 are possible to fulfill the cost function requirement. This result is confirmed by the analysis of the minimization PS1234 in Section 5.1.

Requirements in parameter estimation
In order to detect the correct curve characteristics of the lift and the drag coefficient, it is necessary to use samples for dives and climbs which result in significantly variable angles of attack. This means using data with maximal informativeness about the system. This is a basic requirement in system identification or machine learning to create models which are valid for the whole operational area [42]. Tests with artificially generated noisy data using a parameter set assumed to be known in Eq. (41) and normally distributed pseudo-random numbers were carried out for minimization PS234 and PS1234.
These tests show that the minimization PS234 delivers significant parameters with a maximal 10% error to the given parameter K D2 and a maximal 1% error to K D0 and Δm 0 with an interdecile range IDR θ of at least 8 ∘ by two residual standard deviations s 2 v . z of 0.01 m/s corresponding to the logged data. In such a data distribution, the Pearson's correlation coefficient r is -0.99. Such a data distribution can be achieved by defining two operating points for the dives and climbs, where the pitch angles differ by about 7 ∘ . This is in line with [12] where three operating points θ = 16 ∘ , 19 ∘ , 27 ∘ for the dives and climbs were defined to detect the induced drag C D2 in Eq. (15). Such a requirement contradicts the energy optimal control of a glider during long-term missions. A compromise has to be reached between the energy optimality and the correct estimation of all glider model parameters during a mission. It would be sufficient when the glider is operated with non optimal operation point for two or three days, to collect enough samples.
Finding significant parameters for all four parameters in minimization PS1234 requires in addition to an interdecile range IDR θ of at least 8 ∘ a two residual standard deviations s 2 v . z of less than 0.001 m/s for the vertical velocity v z which corresponds to a tenth of the observed deviation in the logged data. Pearson's correlation coefficient r is here -1.0. This requirement is difficult to realize in practice (existing vertical current velocity v zcurrent and measurement and detection errors in θ and v z ), which makes the estimation of all four parameters impossible. Another alternative is to use known model parameters and estimate the unknown ones, which is described in Section 5.4.1.

Depth-average velocity analysis for RU29-492
The logged data from the glider mission RU29-492 from Perth, Australia to Sri Lanka [47] for a time period of seven months from 2016-11-17 to 2017-06-24 was used to analyze the depth-average velocity during a long-term mission which is described in Section 6. The horizontal glider velocity required for this could be calculated with time-invariant and time-varying model parameters. The time-varying model allows the assessment of biofouling during the mission. The individual steps to determine the horizontal glider velocity v x using a time-invariant and a time-varying model are shown in Fig. 24. The time-invariant model uses only one parameter set as result of the minimization of Eq. (43) using all logged glider data for the dives and climbs of the whole mission period, whereas the time-varying model estimates a parameter set for each dive and climb at time t(i) using only logged glider data which lie in a sliding window defined by a time interval between t i t ( ) 0.5 window and + t i t ( ) 0.5 window . The length of the sliding window is t window .
The data of the glider state during the time period is shown in Fig. 25. For most of this period, the glider was programmed to fly four yo profiles, starting from the surface and diving to 980 m, climbing to 100 m, diving again to 980 m, repeating this two additional times and returning to the surface. During the surface period of normally 10 min the glider sends the logged data via satellite to the Dockserver, the landor ship-based glider communication center, and receives new instructions for command and control for the next dive period.

Parameter estimation
Due to the existing distribution of the logged data a reliable parameter estimation for all parameters in Eq. (41) is not appropriate (see Section 5.3). Therefore, the quadratic parameter K D2 in drag coefficient Eq. (37) and the linear parameter K L1 in lift coefficient Eq. (38) are assumed to be known and were used from Cooney [14] (see Table 1). This strategy is equivalent to the minimization set PS24 in Section 5.1 and corresponds to the works of Graver [28] and Merckelbach [17] where only a part of the model parameters were estimated. A close look at factors that influence both parameters shows that a change in value as a result of biofouling is unlikely. The linear parameter K L1 only depends on constructive parameters/relations (aspect ratio , wing sweep angle Λ) where biofouling has no influence. Since in the liftingline theory, the quadratic parameter K D2 is calculated using the lift coefficient K L (see Eqs. (32) and (33)) and the aspect ratio , a change in value as a result of biofouling is therefore also unlikely.
The approach described in Section 4.1 was used to estimate the unknown parameters K D0 and Δm 0 . The first minimization estimates these two parameters using all logged data for the dives and climbs of the 7-month mission period in Eq. (43). This results in two single model parameters K D0 = 7.355 kg/m and Δm 0 = -0.085 g, which are timeinvariant during the whole mission period. The modeled vertical glider velocity v z for dive and climb using these parameters in Eq. (41) are shown in curve "Time-invariant Model" in Fig. 26. The comparison of the logged and the modeled vertical glider velocity between dive and climb shows differences. There is a good match of the curves for climb for the whole mission period, where the vertical glider velocity is approximately equal to 0.16 m/s during the whole mission period. The modeled curve for the dive shows a relatively constant progression of 0.14 m/s, whereas the logged vertical glider velocity has larger values at the first half of the mission and smaller values in the second part of the mission. The real progression can be approximated as a linear slope starting at 0.17 m/s and ending at 0.14 m/s. The consistent commanded pitch angle θ and the variable ballast mass m b (see Fig. 25) cannot explain the contradictory progression of the vertical glider velocity for the dive (time-varying progression) and the climb (time-invariant progression) during the mission.
To analyze the time-varying behavior of the vertical glider velocity for dive, an extended minimization process will be performed. In this process, the minimization will be executed for every mission day using the logged data over a sliding window of 3 weeks in Eq. (43). The resulting curves for the parameter K D0 and Δm 0 are shown in Fig. 27. The glider was trimmed slightly, which results in a buoyancy trim offset of Δm 0 = -65 g at the start of the mission. This was necessary to compensate for the lower density of surface waters starting at Sumatra and most likely near Sri Lanka. The line "Time-varying Model" in Fig. 26 shows the course of the vertical glider velocity v z using these parameters in Eq. (41). Now both modeled curves for dive and climb match very well with the real glider behavior. The reason for this good match is the trends of the parameters K D0 and Δm 0 during the mission. The parameter K , D0 which corresponds to the parasitic drag coefficient and the negative buoyancy trim offset Δm 0 increase with time. Biofouling could be responsible for this behavior. The biofouling grows during the mission, which increases the skin friction drag and generates an additional buoyancy. This is also described in [13], where the influence of biofouling on the glider behavior during long-term missions is explained in detail. The 0.4 kg/m 3 increase in water density during the mission (see Fig. 25) can also explain the increase of Δm 0 by about 30 g. (The used glider had an extended energy bay, which results in a glider volume of 72.2 l. This volume leads to an additional buoyancy of 28.88 g = 0.4 kg/m 3 · 0.0722 m 3 during the mission.) These two effects, a larger drag coefficient and a larger negative buoyancy trim offset, add up during dive and reduce the glide speed. A larger negative buoyancy leads to a larger glide speed during climb. This effect will be compensated for by the larger drag coefficient, which decreases the glide speed.
The result is an approximately equal glide speed for climb during the mission. Fig. 28 shows the horizontal glider velocity and an estimated linear trend for dive and climb using the parameters from the time varying model. The estimated linear trend equations, computed with regression,   • The horizontal glider velocity v x for dive and climb can be different.

Resulting velocities
This means to obtain a similar horizontal glider velocity for dive and climb during the mission the current buoyancy trim offset Δm 0 needs to be detected as accurately as possible. A good estimate of this parameter can be achieved using Eq. (51). Fig. 29 shows the curves estimated using this approach ("Calculated Parameter") and the results of the moving minimization ("Estimated Parameter") as described above. The sliding window for the logged data is 3 weeks in both approaches. Since the maximum error between the curves is only 3.9 g, the calculated buoyancy trim offset is a good guide to determine the right variable balance mass m b for climb and dive during the mission. A more accurate result can be achieved using an angle of attack α in Eq. (51). Fig. 30 shows the calculated angle of attack courses using the model parameters from Cooney, a time-invariant and a time-varying model and the logged pitch angle θ in the angle of attack determination presented in Section 4.2. The angle of attack α is approximately equal to ± 4 ∘ during the mission using the parameters from the time-invariant model. Using the time-varying model, the angle of attack will increase from ± 3.6 ∘ to ± 4.4 ∘ during the mission. Although Cooney's model does not include time-varying parameters a slight increase of the angle of attack α is observed for dives during the mission. The reason is the inclusion of the glide speed V in the pressure drag C D0 calculation in Eq. (34), which is used in the angle of attack determination. The glide speed V is decreased for dives and remains constant for climbs during the mission. Fig. 31 shows the calculated horizontal glider velocity v x using the angles of attack in Fig. 30, the vertical glider velocity v z and the pitch angle θ in Eq.

Evaluation of the depth-average velocity
The depth-average velocity is a unique calculated quantity in the    operating area of underwater gliders. It is used for glider navigation or as a score to evaluate the quality of ocean current models in and near the operating area of the glider. The evaluation of the ocean current models can help to choose a suitable current model for navigation during the mission. Another application for evaluation is the data postprocessing, where the depth-average velocities of the individual dive segments in a mission will be compared with ocean current models. This evaluation can help to determine the confidence in new data sources, to improve the ocean models using the spatial and temporal anomalies between modeled and observed data, and to modify an ocean model with additional data/information from a glider or calculation methods from ocean models with a better correlation to the logged data. The individual steps for an evaluation process are presented below.

Calculation of depth-average velocity
To evaluate the depth-average velocity of an ocean current model with that of a glider requires their previous calculation. Both calculation methods for a glider and an ocean model will be described in detail below.

Depth-average velocity of a glider
The depth-average velocity of a glider v cGlider is the difference between the velocity over ground v g and the horizontal glider velocity through water v h [15] = The velocity over ground v g is the result of the GPS fixes at the beginning x start and the end x end of a dive segment and the time required for it The horizontal glider speed through water v h can be described by its magnitude, the horizontal glider velocity v x , and its direction, the logged glider heading ψ As shown in Section 5.4.1, there are two different horizontal glider velocities for dive v xdive and climb v xclimb . To calculate an exclusive horizontal velocity v x used in Eq. (64) a temporal weighting where the velocities are active is necessary. The horizontal glider velocity for dive would be underestimated and the horizontal glider velocity for climb would be overestimated using a simple mean calculation, because the dive time period is larger than the climb time period. The vertical glider velocities v zdive and v zclimb will be used for temporal weighting

Depth-average velocity of an ocean model
The depth-average velocity can be understood as the mean ocean current in the operation area of the glider during the time period for a considered dive segment. Such a segment can consist of one or more yo profiles. To calculate this value, the ocean current conditions v c have to be detected at the position x start , at the time t start and at the position x end , at the time t end in several depth layers for the dive segment. Therefore, a defined number n of depth layers are equidistantly distributed between the surface and the dive-to depth which is shown in Fig. 32. The number n is taken to be 20. Finally, the mean value for all these ocean current components will be calculated according to the following equation Since the ocean current data, coming from the Ocean General Circulation Models (OGCM) as data files, will be provided only at discrete positions, depths and times with a nonlinear depth scale and a coarser time and length scale, so that the ocean current information cannot be taken directly from the files. Hence, a multi-dimensional interpolation scheme will be used to extract the desired ocean current data. Additional information about the interpolation scheme is presented in [16]. This approach is a simplification and works well for easy current conditions. In case of complex current conditions or to determine the depth-average velocity more accurately, a glider simulation should be used. Therefore, the exact yo movement of the glider in the current field of the ocean model using the commanded heading, pitch and dive profile for the dive segment will be simulated. The resulting surfacing position and the simulated horizontal glider velocity can then be used in Section 6.1.1 to determine the depth-average velocity.

Goodness-of-fit indicators
Two goodness-of-fit indicators were used for the evaluation. Both indicators use the distance information d between the calculated end positions using the depth-average velocity of the ocean model v cModel and the glider v cGlider after a defined time period Δt.
It should be noted that the horizontal glider velocity through water v h has no influence on Eq. (67) because this vector will be used in both calculations for glider and for ocean model, and thus will be offset. The first indicator is a non-dimensional skill score (ss) presented in [43]. This skill score was developed to compare simulated with observed drifter trajectories and has to be adapted for single dive segments. It has been applied in a variety of work [44][45][46] in the last decade. To overcome the difficulties with dive segments of different lengths a non-dimensional index s will be used to normalize the distance d with the length l of the dive segment where n is a tolerance threshold. A skill score of 1 implies a perfect match between model and glider. The second indicator is novel and uses a time period Δt of 24 hours in Eq. (67). The result of this indicator Da24h is given in km. This makes it possible to overcome the difficulties with dive segments of different time periods. Likewise, this parameter can be imagined easily, as it corresponds to the distance between the modeled and the real surfacing position of the glider after 24 hours. A small distance means a good match between model and glider.

Evaluation
Three ocean current models were chosen for the comparison of the modeled with the observed depth-average velocities along the mission route. The models used are GLBu008 (GOFS 3.0) and GLBv008 (GOFS 3.1) of HYCOM [19] and the 24 hourly global-analysis-forecast-phys-001-024 model of the Copernicus Marine Environment Monitoring Service (CMEMS) [18]. Fig. 33 shows a map with the mission route of RU29-492 from Perth, Australia to Sri Lanka with the weekly mean depth-average velocity vectors from the glider RU29 (blue) using the calculated horizontal glider velocity v x of the time-varying model and the 24 hourly CMEMS ocean model (red) for the several dive segments. This map includes additionally the smoothed skill score ss values for the several dive segments. A moving average algorithm with a sliding window of seven days was used for smoothing. There is a good match between the observed and modeled depth-average velocities for most of the mission period where the skill score ss is larger than 0.7. There are four regions during the mission with a bad evaluation:   34 shows the smoothed courses of the two goodness-of-fit indicators: skill score ss and distance after 24h Da24h for the three ocean current models using the calculated horizontal glider velocity v x of the time-varying model. The distance after 24h Da24h shows an inverse behavior to the skill score ss. The HYCOM models show a worse evaluation in comparison to the CMEMS model. As shown in Table 4, which includes the mean values of the indicators, the GLBv008 model has a slightly worse evaluation in comparison to the GLBu008 model. The bad evaluation of the HYCOM models in the area are noticeable: • 14.5 ∘ -15.5 ∘ S and 101 ∘ -103 ∘ E start of February 2017 (only GLBv008) • 0 ∘ S and 85.5 ∘ -87 ∘ E middle of June 2017. Table 4 shows that the mean values of the two indicators for the methods used to detect horizontal glider velocity v x are very similar. The higher model accuracy of the time-varying model does not lead to other significant results. The reason is the small differences in the calculated horizontal glider velocities in the individual methods (maximal 5%) (see Fig. 31). The CMEMS has a 20% better quality in ss and a 33% better quality in Da24h than the indicators of the HYCOM models. The GLBu008 shows a 6% better quality in ss and a 11% better quality in Da24h than the GLBv008 model.

Conclusions
This paper examines the need for an accurate glider model in navigation, control and data post-processing. A model which is generally used for a Slocum glider was presented in detail. The calculation of the horizontal glider velocity v x based on an accurate angle of attack detection was discussed in this context. A robust minimization approach to detect the model parameters and its limitations due to the available logged glider data were described. The correct parameter estimation requires a large data distribution in the whole operation range of a glider which contradicts the energy optimal control of a glider during   34. Skill score ss and distance after 24h Da24h curves for the three ocean models: 24 hourly global-analysis-forecast-phys-001-024 model of CMEMS (CMEMS 24HOURLY), GLBu008 (HYCOM GlBu008) and GLBv008 (HYCOM GLBv008) of HYCOM.

Table 4
Calculated mean values for skill score ss and distance after 24h Da24h for the three ocean models: 24 hourly global-analysis-forecast-phys-001-024 model of CMEMS, GLBu008 and GLBv008 of HYCOM and the three methods to detect the horizontal glider velocity v x . long-term missions. The analysis of RU29-492 from Perth, Australia to Sri Lanka shows an increase in the drag coefficient over the course of the mission caused by biofouling, which leads to an increase of the angle of attack of about 1 ∘ using similar pitch commands. In addition, an increase in the buoyancy trim offset as result of biofouling and changing salinity conditions during the mission could also be observed. This should be considered for the navigation and control of a glider through an online model parameter estimation during the mission. The decrease of the horizontal glider velocity from 0.3 to 0.26 m/s within 7 months should be considered in the global path planning. The comparison of the depth-average velocity of ocean current models from CMEMS and HYCOM with the observed velocity from the glider using chosen goodness-of-fit indicators allows a spatial and temporal model evaluation along the mission route. In this context, one question arises: How should the evaluation results be processed so that they can be used in ocean model design?