An Uncertainty-Aware Assessment Framework for Multi-Glider Cooperative Observation of Mesoscale Eddies

1. Introduction

Mesoscale eddies are ubiquitous rotating structures in the global ocean, typically spanning 50–300 km in diameter and persisting for weeks to months [1,2]. These coherent vortices play a critical role in the redistribution of oceanic energy, heat, salt, and biogeochemical tracers, significantly influencing large-scale ocean circulation, regional climate variability, and marine ecosystem dynamics [3,4,5,6]. Understanding their three-dimensional structure is essential for improving ocean prediction models and quantifying their contribution to global energy budgets [2,7]. Satellite remote sensing detects surface signatures but lacks subsurface information [1,8]. Subsurface profiling platforms, including Argo floats, provide depth-resolved measurements but lack the spatial resolution and maneuverability to resolve the fine-scale internal structure of individual eddies [9]. This observational gap motivates the deployment of coordinated mobile platforms capable of three-dimensional sampling within targeted eddy features.

Autonomous underwater gliders have emerged as a promising platform for in-situ eddy observation, owing to their extended endurance, depth-profiling capability, and programmable navigation [10,11,12]. Recent studies have demonstrated their effectiveness in capturing eddy internal structures [13,14]. Deploying multiple gliders in a coordinated network has been validated in field experiments [15], and enables systematic sampling of an eddy’s spatial structure through cross-sectional transects [16,17]. Several path planning strategies have been proposed, including butterfly patterns that traverse the eddy core along intersecting diameters [18], zigzag lawn-mower patterns for area coverage [19] and adaptive approaches such as heading control for eddy tracking [20] and real-time waypoint adjustment [21,22] . These methods optimize controllable parameters, such as glider speed, heading angle, centerline orientation, and intersection geometry under the assumption of a known, deterministic eddy environment [23]. However, real ocean eddies are inherently variable: their intensity, vertical structure, and ambient current conditions cannot be precisely known prior to deployment. Whether parameters optimized for one eddy condition—glider speed, in particular—remain effective across a range of eddy intensities is an open question.

Existing studies on glider-based eddy observation typically assume a known environment when evaluating path designs [18,24]. The impact of ocean uncertainty on glider performance has been investigated at the individual platform level. Ocean currents alter glider motion dynamics and trajectory accuracy [25,26], time-varying currents challenge path planning [27,28], and stochastic conditions affect energy consumption [29,30]. These studies address trajectory-level errors under specific uncertainty sources. However, mesoscale eddy observation is a mission-level task, whose outcome is measured by the quality of reconstructed oceanographic fields, spatial coverage, shape and center accuracies. How multi-source uncertainties propagate from individual glider trajectories to these mission-level metrics remains largely unexplored. Systematic frameworks integrating adaptive sampling and uncertainty quantification have been developed for ocean observation [31], but their application to multi-glider eddy observation under multi-source uncertainty has not been addressed. A systematic framework that integrates uncertainty propagation, multi-dimensional evaluation, and parameter sensitivity analysis for glider network observation is lacking.

To address these limitations, this paper proposes the SCOPE (Stochastic Cooperative Observation Performance Evaluation) framework. The principal contributions are:

1.

SCOPE framework. A closed-loop simulation framework couples a parameterized multi-scale eddy environment with a kinematic glider model, propagating multi-source uncertainties to a 27-metric evaluation system spanning seven quality categories. An Adaptive Core Metric Selection (ACMS) algorithm choose core metrcis them while preserving category coverage.

2.

Sensitivity analysis and ablation diagnosis. Correlation and interaction analyses of the 10 uncertain parameters yield a ranked importance hierarchy and identify three physically coupled parameter pairs. Derived variables constructed from these couplings reveal a non-monotonic, intensity-dependent effect of the dominant parameter. Structured ablation experiments validate these findings by decomposing marginal and combined effects of the identified design rules.

3.

Methodology transfer to real ocean eddies. The ACMS algorithm and sensitivity ranking methodology are validated on four HYCOM eddies of different intensities. The core metrics and parameter rankings obtained by the same SCOPE methodology are compared between the analytical model and real eddies. The velocity ratio trend observed in the analytical model is reproduced in real eddies, with a systematic shift in magnitude attributable to differences in environmental complexity.

The remainder of this paper is organized as follows. Section 2 presents the SCOPE framework. Section 3 reports the two-stage sensitivity analysis, dose-response characterization, and ablation validation. Section 4 examines the framework’s transferability to real ocean eddies. Section 5 discusses implications and concludes.

2. The SCOPE Framework

2.1. Framework Overview

The SCOPE (Stochastic Cooperative Observation Performance Evaluation) framework is proposed to quantify and diagnose multi-glider observation performance under multi-source uncertainty. Figure 1 shows its workflow.

First, a Parameterized Eddy Model (Section 2.2) provides gridded temperature and velocity fields under main uncertain parameters. Second, the Uncertainty-Integrated Cooperative Motion Model (Section 2.3) propagates these uncertainties together with controllable mission parameters and uncontrollable motion factors, outputting glider trajectories and temperature measurements via discrete-time integration. A Monte Carlo experiment is then conducted by Latin Hypercube Sampling (LHS) of the uncertain parameters. Each LHS-sampled mission is scored by 27 metrics across seven quality categories, and ACMS algorithm compresses them preserving ranking fidelity (Section 2.4). From the composite score, the framework supports global sensitivity analysis, parameter importance ranking, and characterization of key parameters (Section 3).

2.2. Parameterized Eddy Field

The eddy environment field is constructed by parameterized velocity field and temperature structure to enable controlled uncertainty propagation and analysis.

2.2.1. Temperature Field

A gridded three-dimensional temperature field (Figure 2a) serves as the observation ground truth. It features a subsurface-intensified warm-core structure consistent with South China Sea anticyclonic eddies, constructed from three superimposed components [32]: a hyperbolic-tangent thermocline background, a warm-core anomaly, and sub-mesoscale perturbations. The dominant warm-core anomaly is modeled as:

$$T_{\mathrm{eddy}}(r,z)=A_T·e^{-\beta r_n^2}·e^{-{(z-z_T)}^2/\left(2H_T^2\right)}·\left[1+\gamma e^{-{(r_n-1)}^2/\left(2{\sigma }_e^2\right)}\right]$$

(1)

where $r_n=r/r_{\mathrm{edge}}\left(\varphi \right)$ is the normalized radial coordinate, $z_T$ is the depth of maximum thermal anomaly, and $H_T$ is the vertical e-folding half-width. The final factor produces enhanced thermal fronts at the eddy boundary [33]. Sub-mesoscale perturbations are added as multi-frequency sinusoidal fluctuations (amplitudes 0.2–0.3°C) to represent filamentary structures and internal waves [34]. Gaussian observation noise ($\sigma =0.08$°C) is superimposed. The eddy boundary is defined as a perturbed ellipse:

$$r_{\mathrm{edge}}\left(\varphi \right)=r_{\mathrm{ellipse}}\left(\varphi \right)+A_{\mathrm{w}3}sin\left(3\varphi \right)+A_{\mathrm{w}8}sin(8\varphi +\pi /4)+{ϵ}_r\left(\varphi \right)$$

(2)

with semi-axes $a=60$ km (aligned with the major axis), $b=51$ km (minor axis), $A_{\mathrm{w}3}=5$ km, $A_{\mathrm{w}8}=2$ km, and ${ϵ}_r\sim \mathcal{N}(0,1.5^2)$ km.

2.2.2. Flow Field

The Flow field (Figure 2b) uses a Lamb–Oseen vortex model [35]: The tangential velocity at radial distance r from the eddy center is:

$$V_{\theta }\left(r\right)=V_{max}·{\displaystyle \frac{r}{R_{max}}}·exp\left[{\displaystyle \frac{1}{2}}\left(1-{\displaystyle \frac{r^2}{R_{max}^2}}\right)\right]$$

(3)

where $V_{max}$ is the maximum tangential velocity at radius $R_{max}$. Vertical decay follows $f\left(z\right)=exp[-{(z-z_m)}^2/\left(2{\sigma }_z^2\right)]$ with depth of maximum swirl $z_m$ and vertical e-folding scale ${\sigma }_z$. A spatially uniform background current $U_c$ with direction ${\theta }_c$ is superimposed to represent mesoscale advection:

$$\mathbf{u}(r,\varphi ,z)=V_{\theta }\left(r\right)f\left(z\right){\widehat{\mathbf{e}}}_{\varphi }+U_c(cos{\theta }_c,sin{\theta }_c)$$

(4)

Vertical velocity (negligible relative to horizontal components) is neglected to simplify computations .

This parameterized formulation allows the key environmental variables—$V_{max}$, $z_m$, $U_c$, and ${\theta }_c$—to be treated as controllable inputs for uncertainty propagation. By varying these parameters within their prescribed ranges (Table 1), the model generates an ensemble of eddy realizations that supports Latin Hypercube Sampling and Monte Carlo experiments.

2.3. Uncertainty-Integrated Cooperative Motion Model

The cooperative motion model simulates fleet trajectories under a single design principle: no adaptive control or real-time optimization is applied, so that observation quality reflects only the propagation of environmental and instrumental uncertainties. The model comprises two coupled layers: an upper mission-planning layer (Section 2.3.1) and a lower motion-integration layer (Section 2.3.2). Uncertain parameters from both layers are combined in an integrated kinematics model that propagates the fleet state through discrete-time stepping (Section 2.3.3).

2.3.1. Upper Layer: Mission Planning

The planning layer generates fleet configuration and sampling waypoints for the motion-integration layer. An edge–center architecture (Figure 3a) is chosen to support comprehensive multi-dimensional evaluation. Edge circumnavigation resolves boundary shape, which is essential for tracking eddy evolution [15,36,37,38,39]. Center transection captures radial cross-sections of temperature and thermocline structure [40,41]. Both modes have been independently validated in prior studies as effective strategies for eddy observation. Together, they cover horizontal and vertical reconstruction, thermocline characterization, and eddy center localization, ensuring that all quality dimensions are observable.

Based on the eddy dimensions from the parameterized environmental field (Section 2.2), an architecture with three edge gliders and two center gliders (Figure 3a) is designed under the principle of spatiotemporal synchronization to reduce spatiotemporal aliasing during observation [10] and to provide a baseline for evaluating synchronization quality under uncertainty.

The two center gliders follow diametrical transects parameterized by the centerline angle ${\alpha }_{\mathrm{cline}}$ and intersection angle ${\alpha }_{\mathrm{inter}}$. A total of $n_c$ sampling waypoints are distributed uniformly along the two transects. These three parameters (${\alpha }_{\mathrm{cline}}$, ${\alpha }_{\mathrm{inter}}$, $n_c$) are treated as controllable uncertain inputs in the LHS experiment, enabling analysis of how mission planning choices propagate to observation quality.

The three edge gliders equally partition the eddy perimeter and sample along the boundary. Their waypoints are determined by a kriging-variance criterion based on equal-angle sectors:

$$n_i=\mathrm{round}\left({\displaystyle \frac{L_i}{\Delta s}}·\left[1+{\lambda }_{\mathrm{K}}·{\displaystyle \frac{{\sigma }_{\mathrm{K},i}^2-{\overline{\sigma }}_{\mathrm{K}}^2}{s_{{\sigma }^2}}}\right]\right)$$

(5)

where $L_i$ is the arc length of sector i ($i=1,\dots ,12$), $\Delta s$ is the base profile spacing, ${\sigma }_{\mathrm{K},i}^2$ is the kriging variance of sector i, ${\overline{\sigma }}_{\mathrm{K}}^2$ is the mean variance across all sectors, $s_{{\sigma }^2}$ is the standard deviation of the sector variances, and ${\lambda }_{\mathrm{K}}=0.5$ controls the sensitivity to variance. Spacing is constrained to $[5,20]$ km.

2.3.2. Lower Layer: Motion Integration

In three-dimensional space, the glider moves at speed $v_g$ in a sawtooth pattern between the surface and maximum depth (Figure 3b). Neglecting the angle of attack, the glide angle is represented by the pitch angle $\theta$, and the horizontal heading by the heading angle $\psi$. The velocity components are:

$$\dot{x}=v_{g}cos\theta cos\psi ,\dot{y}=v_{g}cos\theta sin\psi ,\dot{z}=v_{g}sin\theta$$

(6)

The heading angle is updated at each surfacing event to follow the planned path. A new heading command ${\psi }_{\mathrm{cmd}}$ is computed from the current surfacing position to the next profile endpoint. A first-order lag model simulates the physical inertia of glider turning:

$$\dot{\psi }=\mathrm{sat}\left({\displaystyle \frac{{\psi }_{\mathrm{cmd}}-\psi }{{\tau }_{\psi }}},{\dot{\psi }}_{max}\right)$$

(7)

where ${\tau }_{\psi }$ is the time constant and ${\dot{\psi }}_{max}$ is the maximum turn rate.

2.3.3. Uncertainty Integration and Propagation

Environmental and instrumental uncertainties are superimposed on the nominal kinematics. The eddy current field $\mathbf{u}(\mathbf{r},z)$ from Section 2.2 adds to the horizontal velocity components:

$$\dot{x}=v_{g}cos\theta cos\psi +u(\mathbf{r},z),\dot{y}=v_{g}cos\theta sin\psi +v(\mathbf{r},z)$$

(8)

Instrumental biases ${\delta }_{\psi }$ and ${\delta }_{\theta }$ are applied to the heading and pitch angles at each surfacing event:

$$\psi \leftarrow \psi +{\delta }_{\psi },\theta \leftarrow \theta +{\delta }_{\theta }$$

(9)

The 10-dimensional uncertain input space consists of three planning parameters, three motion parameters, and four environmental parameters (Table 2). Feeding these parameters into the cooperative motion model yields the fleet sampling plan and glider trajectories, as illustrated in Figure 4 under mean parameter values.

2.4. Performance Assessment and Core Metric Selection

An LHS [42] sample matrix is generated with $N=5000$ draws in MATLAB. Each draw simulates a complete mission and computes 27 evaluation metrics, producing the response matrix $\mathbf{Y}\in {\mathbb{R}}^{5000\times 27}$. A composite score is first constructed as a scalar measure of overall observation quality. The ACMS algorithm then compresses the metric set to a compact core subset that preserves the score-based ranking fidelity.

2.4.1. Multi-Dimensional Metric System and Composite Score

Seven quality categories of metrics are designed,including trajectory error by uncertainty, eddy center localization accuracy, eddy boundary fidelity, horizontal temperature mapping quality, vertical section profiling quality, energy efficiency for sustainability assessment, and fleet synchronization to quantify operational consistency (see Supplementary Material, Appendix A for the complete derivation). Each metric is normalized to $[0,1]$ according to its directionality, with bounds fixed from the 5000-sample ensemble. The composite score is the equally weighted mean of the 27 normalized metrics, serving as the reference for metric selection. The robustness of equal weighting against alternative schemes (entropy-weight, category-weight, CV-weight, and random) is validated in Supplementary Material, Appendix B.

2.4.2. Adaptive Core Metric Selection

The ACMS algorithm is designed to reduce redundancies among metrics while maximally preserving the multi-dimensional evaluation information. Starting from all 27 metrics, ACMS iteratively removes the metric whose elimination causes the smallest drop in Spearman rank correlation ${\rho }_S$ with the reference composite score, subject to the constraint that each of the seven categories retains at least one representative:

$${\mathcal{M}}^{*}=arg\underset{\mathcal{M}}{min}\left|\mathcal{M}\right|\mathrm{s}.\mathrm{t}.{\rho }_S(\mathbf{S}\left(\mathcal{M}\right),{\mathbf{S}}_{\mathrm{ref}})\ge 0.95$$

(10)

ACMS operates in three phases (full pseudo-code in Supplementary Material, Appendix C). Phase 1 pre-filters by coefficient of variation (CV $>0.10$) and enforces category coverage. Phase 2 iteratively removes the metric whose elimination causes the smallest drop in Spearman correlation ${\rho }_S$ with the reference composite score, stopping when ${\rho }_S<0.95$ or further removal would empty a category. Phase 3 computes the final composite score from the retained metrics.

Eighteen elimination steps yield $|{\mathcal{M}}^{*}|=9$ core metrics with ${\rho }_S=0.967$ (Figure 5), a 67% reduction while preserving all seven categories. The nine core metrics are: net swing error, time CV, temperature-weighted center error, Procrustes rotation angle, thermocline RMSE, surface max error, vertical max error, vertical grid coverage, and energy per unit distance. Five-fold cross-validation confirms no overfitting (details in Supplementary Material, Appendix D). The 9-metric composite score is used for all subsequent analyses.

3. Sensitivity Analysis Under Uncertainty

The effects of the 10 uncertain parameters on observation quality are analyzed in this section. First, correlation and interaction analysis is performed on the original parameters to rank their importance and identify physical couplings. Second, derived variables are constructed to resolve these couplings, followed by dose-response characterization of the dominant parameters. All results are based on the 5000-sample LHS ensemble and the 9-metric composite score from Section 2.4.

3.1. Correlation and Interaction Analysis

Parameter importance is assessed using Spearman’s rank correlation ${\rho }_S$ and k-NN mutual information (MI, $k=5$) [43,44]. With $N=5000$ samples, even small effects achieve nominal significance. Practical significance is therefore assessed using Cohen’s d [45], with results reported as meaningful only when $p<0.05$ and $\left|d\right|\ge 0.10$ [46]. As shown in Figure 6, $V_{max}$ and $v_g$ are the top-ranked parameters, suggesting a strong coupling between them. This motivates the interaction analysis in Section 3.2.

3.2. Interaction Analysis and Derived Variables

To resolve the couplings suggested by the correlation analysis, all $\left({\displaystyle \genfrac{}{}{0pt}{}{10}{2}}\right)=45$ pairs of the 10 input parameters are assessed using factorial ANOVA with Bonferroni correction (${\alpha }^{*}=0.05/45=0.0011$). Three pairs exhibit notable interactions (Figure 7): $V_{max}\times v_g$ (${\eta }^2=0.035$), $U_c\times v_g$ (${\eta }^2=0.010$), and ${\theta }_c\times {\alpha }_{\mathrm{cline}}$ (${\eta }^2=0.007$). The $V_{max}\times v_g$ interaction is the strongest, confirming that glider speed and eddy intensity share a coupled physical mechanism. The ${\theta }_c\times {\alpha }_{\mathrm{cline}}$ interaction reveals that performance depends on the relative angle between background current and transect heading. The $U_c\times v_g$ interaction (${\eta }^2=0.010$) does not admit a simple derived variable and $U_c$ is rarely identifiable in real ocean environments, so this coupling is not further resolved. The remaining 93% of interactions have ${\eta }^2<0.005$, indicating that the system is approximately additive once the dominant couplings are resolved.These couplings motivate the construction of physically independent derived variables. The velocity ratio:

$${\eta }_v={\displaystyle \frac{v_g}{V_{max}}}$$

(11)

absorbs the $v_g\times V_{max}$ interaction into a single dimensionless parameter. The background current is decomposed relative to the transect heading:

$$U_c^{\perp }=U_{c}sin({\theta }_c-{\alpha }_{\mathrm{cline}}),U_c^{\Vert }=U_{c}cos({\theta }_c-{\alpha }_{\mathrm{cline}})$$

(12)

where $U_c^{\perp }$ measures lateral drift and $U_c^{\Vert }$ measures headwind/tailwind effects [43], resolving the ${\theta }_c\times {\alpha }_{\mathrm{cline}}$ coupling into physically interpretable components.

The sensitivity analysis is repeated with the derived variables using the same methods. ${\eta }_v$ ranks first by MI (0.375). its low $|{\rho }_S|=0.141$ relative to the high MI confirms a non-monotonic effect: both excessively low and excessively high ${\eta }_v$ degrade performance. The dose-response characteristics of ${\eta }_v$ are examined in Section 3.3 to characterize this relationship. $U_c^{\perp }$ ranks second among environmental factors ($|{\rho }_S|=0.163$, MI $=0.037$), indicating that lateral drift is the dominant pathway through which background currents affect observation quality, while the along-track component $U_c^{\Vert }$ shows negligible influence ($|{\rho }_S|=0.004$). Table 3 compares the two parameter groups, providing a ranked reference for mission design. Parameters with $p\ge 0.05$ are marked as “Not sig” (not statistically significant).

3.3. Dose-Response Characterization

The dominance of ${\eta }_v$ motivates a detailed examination of its effect on observation quality. The 5000-sample LHS ensemble is stratified by ${\eta }_v$, and the mean composite score is computed for each bin.

The dose-response curve reveals a non-monotonic relationship with an optimal plateau (Figure 8a). The peak occurs near ${\eta }_v\approx 0.5$, with the 95%-of-peak interval spanning $[0.40,0.76]$. This wide plateau indicates that moderate deviations from the optimal value incur limited performance penalties, providing practical tolerance for speed selection under imperfect environmental knowledge. Score distributions shown in Figure 8b are tightest near the optimal ${\eta }_v$ and widen toward both extremes, indicating that the optimal configuration is robust to small speed variations.

Since the dose-response curve exhibits attenuation at both low and high ${\eta }_v$, the LHS samples are further stratified by $V_{max}$ and $U_c$ to examine whether the optimal ${\eta }_v$ depends on eddy intensity. Figure 9a reveals a systematic shift: the optimal ${\eta }_v$ decreases monotonically with $V_{max}$ (from ∼0.58 for weak eddies to ∼0.47 for strong eddies, Pearson $r=-0.87$), and the optimal interval narrows. The overall pattern is that weak eddies favor higher ${\eta }_v$ to ensure sufficient spatial coverage, while strong eddies benefit from lower ${\eta }_v$, allowing the glider to move with the flow and accumulate more profiles within the high-velocity region. Figure 9b shows the optimal ${\eta }_v$ is largely insensitive to $U_c$, as the dose-response curves for different background current magnitudes nearly overlap. This confirms that $V_{max}$, not $U_c$, is the primary environmental driver of ${\eta }_v$ selection. This trend is further examined on real eddies in Section 4.2.

3.4. Ablation Validation

The sensitivity analysis above identifies two actionable design rules: (i) match ${\eta }_v$ to eddy intensity within the optimal plateau, and (ii) align ${\alpha }_{\mathrm{cline}}$ perpendicular to the background current direction. To test whether these rules causally improve observation quality rather than merely correlate with it, a structured ablation experiment is conducted under a fixed environment ($V_{max}=0.85$ m/s, $U_c=0.25$ m/s). Five groups are compared under $N=200$ paired uncertainty draws (Table 4): A (deterministic ceiling via grid search), representing the optimal design under certainty; B (literature baseline, $v_g=0.30$ m/s) [10], representing the degradation of A under uncertainty; C (${\eta }_v$ rule, $v_g$ sampled from ${\eta }_v\sim \mathcal{U}[0.40,0.60]\times V_{max}$, a suitable interval for the moderate eddy intensity used in this experiment, Section 3.3); D (${\alpha }_{\mathrm{cline}}={\theta }_c+{90}^{\circ }$), testing the alignment rule identified in Section 3.2; and E (C and D combined). All groups share identical secondary parameters from Group A. Paired differences are assessed using the Wilcoxon signed-rank test [47]. The Normalized Gap Closure $\mathrm{NGC}=({\overline{S}}_{\mathrm{group}}-{\overline{S}}_B)/({\overline{S}}_A-{\overline{S}}_B)\times 100\%$ quantifies the fraction of the A–B gap recovered by each strategy [48].

Figure 10a shows the score distributions. All group means lie below the deterministic ceiling, and the ceiling falls within the upper tail of the distributions, indicating that the grid search may not have fully captured the parameter space. Group B exhibits the largest variance, confirming the baseline strategy as the least robust. Groups C–E all show markedly reduced variance. Figure 10b,c show the paired differences and NGC values (Table 4).

The A→B gap (0.059) quantifies a 6.6% performance penalty for using a fixed literature-default speed. While the absolute score differences are modest on the $[0,1]$ scale, the NGC metric normalizes this gap relative to the achievable improvement, isolating the contribution of each design rule. Whether a 6.6% improvement translates to operationally meaningful gains depends on the specific mission’s tolerance for reconstruction error. The ${\alpha }_{\mathrm{cline}}$ rule alone recovers 45.8% of this gap ($d=0.52$, $p<{10}^{-11}$; effect sizes by Cohen’s d), outperforming the ${\eta }_v$ rule (25.0%, $d=0.27$). The combined strategy (37.8%) shows mild negative synergy, demonstrating the framework’s ability to detect suboptimal parameter interactions. All rule-based strategies reduce performance variability by 36–44%, providing downside protection independent of mean improvement.

The ablation yields three practical insights. First, transect orientation is a high-leverage design choice, recovering nearly half of the performance gap. Although perpendicular alignment maximizes lateral drift, the drift disperses sampling points over a wider area, improving spatial coverage and temperature reconstruction; the composite score registers a net gain because these coverage benefits outweigh the associated trajectory penalties. Second, the ${\eta }_v$ rule provides a smaller but consistent improvement, robust to the specific choice of interval within the optimal plateau. Third, the combined strategy does not add the two improvements, revealing a synergy structure otherwise invisible to single-factor analysis—an advantage of the framework’s ability to decompose joint effects.

4. Methodology Transfer to Real Ocean Eddies

The analyses in Section 3 are based on the idealized temperature model and Lamb–Oseen eddy model. These analytical simplifications—axisymmetric vortex geometry, uniform background current, and a parameterized temperature anomaly—enable controlled uncertainty propagation but may produce findings that differ quantitatively from real ocean conditions. To evaluate whether the SCOPE framework’s core methodology (ACMS metric selection, parameter ranking, and ${\eta }_v$ dose-response characterization) transfers to non-idealized environments, four real mesoscale eddies in the South China Sea (Figure 11) are extracted from the HYCOM GLBy0.08 reanalysis product [49,50] (1/12° resolution). Table 5 summarizes their properties.

The eddies span $V_{max}=0.28$–$0.98$ m/s, covering weak to strong regimes. Eddy boundaries are identified using a multi-criterion consensus method combining the Okubo–Weiss parameter, stream function closed contours, and velocity e-folding contour (see see Supplementary Material, Appendix E for computational details). The actual 3D velocity field from HYCOM replaces the analytical model in the glider motion simulator, while all other modules remain unchanged. The Lamb–Oseen radial profile correlates with the real eddy velocity fields at a mean Pearson $r=0.918$ across the four eddies (computed on azimuthally averaged radial profiles), and the radial temperature anomaly profile correlates at $r=0.981$, confirming that the analytical model captures the essential velocity and temperature structure. Three eddy centers are shown in Figure 11: the temperature-derived center (red star), the vorticity-derived center (white cross), and the surface cold-water center (green circle) of E4, which is a surface-cold anticyclonic eddy. The temperature-derived center is used as the reference for eddy center localization metrics (Section 2.4).

4.1. ACMS Framework Transferability and Parameter Ranking

The ACMS algorithm is independently applied to each real eddy to test whether the SCOPE assessment framework operates meaningfully on non-idealized data. For each eddy, 500 LHS simulations using the actual HYCOM 3D velocity field yield a response matrix $\mathbf{Y}\in {\mathbb{R}}^{500\times 27}$, with the same 27 candidate metrics used in the analytical model. All four eddies achieve ${\rho }_S>0.95$ (Table 6), confirming reliable convergence on non-idealized data. We note that the 500-sample ensemble per real eddy, while sufficient for Spearman correlation estimation (statistical power $>0.80$ for $|{\rho }_S|\ge 0.12$), may yield conservative mutual information estimates for weak parameter dependencies. The consistent ranking of ${\eta }_v$ across all four eddies and its agreement with the analytical 5000-sample result (MI rank = 1) mitigate this concern.

The numbers of retained core metrics range from 13 to 20, larger than the analytical model’s 9, reflecting the richer variability structure of real eddies: asymmetric temperature fields and complex flows activate additional quality dimensions. The variation in $|{\mathcal{M}}^{*}|$ across eddies reveals a systematic relationship with eddy intensity. E3 (moderate, $V_{max}=0.51$ m/s) retains the largest set (20 metrics) because because its velocity field lies within the analytical model’s calibrated range and its intermediate intensity yields sufficient performance variation to make many metrics discriminable. E1 (weak, $V_{max}=0.28$ m/s) retains only 13 metrics because the weaker velocity field reduces the dynamic range of trajectory and shape metrics, pushing more metrics below the CV threshold. E4 (strong) occupies an intermediate position (15 metrics) because its intense velocity field amplifies trajectory errors but suppresses the discriminability of energy metrics—the buoyancy subsystem dominates energy consumption regardless of parameter choices. This intensity-dependent metric activation is a feature of the framework rather than a shortcoming: it reveals which quality dimensions are practically monitorable for a given eddy regime, guiding mission designers toward a leaner, regime-appropriate evaluation set.

Among the 9 analytical core metrics, 5 are stably reproduced across ≥3 of 4 eddies: vertical max error and vertical coverage (4/4), Procrustes rotation, thermocline RMSE, and energy per km (3/4). E3 ($V_{max}=0.51$ m/s, closest to the analytical model’s range) achieves 89% overlap (8/9). The complete per-eddy core metric lists are provided in the Supplementary Material, Appendix F.

The rankings of the uncertain parameters in real eddies are then analyzed by the method proposed in Section 3.1. Cross-environment consistency is assessed using Kendall’s coefficient of concordance [51]:

$$W={\displaystyle \frac{12{\sum }_{j=1}^k{(R_j-\overline{R})}^2}{p^2(k^3-k)}}$$

(13)

where $R_j$ is the rank-sum for parameter j across $p=4$ judges (the four real eddies) and $k=6$ parameters. The result $W=0.593$ (${\chi }^2=11.86$, $p=0.037$) indicates statistically significant agreement among the four real eddies. We note that the observed concordance is moderate and the p-value is near the conventional threshold of 0.05. This reflects genuine variability in parameter importance across eddy regimes rather than a failure of the methodology. For instance, ${\alpha }_{\mathrm{cline}}$ ranks second overall but drops to fifth in E4—the strong, near-axisymmetric eddy where centerline orientation matters less because the velocity field is more radially symmetric. Caution is therefore warranted when extrapolating parameter rankings from one eddy regime to another; the framework identifies which parameters are consistently important across regimes (${\eta }_v$) and which are regime-dependent (${\alpha }_{\mathrm{cline}}$, $n_c$).

The global ranking is ${\eta }_v>{\alpha }_{\mathrm{cline}}>n_c>{\delta }_{\theta }>{\delta }_{\psi }>{\alpha }_{\mathrm{inter}}$, with ${\eta }_v$ ranking first in 2 of 4 eddies and within the top 3 in all cases. Including the analytical model as a fifth judge yields $W=0.515$ ($p=0.025$). The dominant parameter (${\eta }_v$) and the least influential parameters (${\alpha }_{\mathrm{inter}}$, ${\delta }_{\psi }$) are consistent across environments, while intermediate rankings vary with eddy intensity. The per-eddy sensitivity rankings are documented in the Supplementary Material, Appendix F.

4.2. ${\eta }_v$ Dose-Response Transferability

The intensity-dependent ${\eta }_v$ dose-response observed in the analytical model (Section 3.3) is examined in the real eddies to investigate the systematic shift of the optimal ${\eta }_v$ with $V_{max}$. The same binned method is applied to the 500-sample LHS ensemble of each real eddy.

Figure 12 summarizes the results. The optimal ${\eta }_v$—marked by the peak of each dose-response curve—decreases monotonically with eddy intensity in both datasets. For the analytical model, the peak shifts from ${\eta }_v\approx 0.58$ (weak eddies) to ${\eta }_v\approx 0.47$ (strong eddies), with Pearson $r=-0.87$ between ${\eta }_v^{\mathrm{opt}}$ and $V_{max}$. For the real eddies, the peak shifts from ${\eta }_v\approx 2.35$ (E1, $V_{max}=0.28$ m/s) to ${\eta }_v\approx 0.35$ (E4, $V_{max}=0.98$ m/s), with $r=-0.96$. The stars in Figure 12a–c mark the peak of each dose-response curve. The dashed lines in Figure 12d mark the LHS-optimal ${\eta }_v$ from the 500-sample ensemble, confirming that the binned peaks are consistent with the empirical optima.

The decreasing trend is qualitatively consistent between simulation and real eddies: stronger eddies consistently favor lower ${\eta }_v$. Quantitatively, however, the real-eddy optimal ${\eta }_v$ values are systematically higher than the analytical predictions, with the largest discrepancy for the weakest eddy (2.35 vs. 0.58 for E1). This offset reflects two differences between the idealized and real environments.

First, the irregular, non-axisymmetric velocity field of real eddies induces larger and less predictable trajectory deviations than the smooth Lamb–Oseen vortex. In weak eddies, a higher glider speed enables the glider to resist lateral drift and maintain course corrections, producing a wave-like trajectory around the nominal path that increases the number of sampling points and benefits eddy center localization and temperature reconstruction. In strong eddies, by contrast, the glider has limited ability to counteract the intense flow; a high speed amplifies energy consumption and trajectory inconsistency without improving coverage. A lower speed allows the glider to move with the flow along a spiral path, accumulating more profiles within the eddy interior at lower energetic cost.

Second, the real eddies’ non-axisymmetric geometry introduces localized advective surges absent in the analytical model, further increasing the effective speed demand in weak-eddy conditions. Despite the numerical offset, the monotonic decreasing trend is robust across both datasets, indicating that the intensity-dependent nature of ${\eta }_v$ selection is a genuine physical characteristic. Calibrating a corrected ${\eta }_v^{\mathrm{opt}}\left(V_{max}\right)$ function using a broader range of real eddy data is left for future work.

Despite the numerical offset, the monotonic decreasing trend is robust across both datasets, confirming that the intensity-dependent trend of optimal ${\eta }_v$ transfers from the analytical model to real ocean environments. Calibrating a corrected ${\eta }_v^{\mathrm{opt}}\left(V_{max}\right)$ function using a broader range of real eddy data is left for future work.

5. Discussion and Conclusions

5.1. Discussion

The SCOPE framework is proposed for evaluating and diagnosing multi-glider cooperative observation of mesoscale eddies under uncertainty. Its core methodology includes a multi-dimensional metric system (27 metrics compressed to 9 core indicators), ACMS adaptive core metric selection, two-stage sensitivity analysis, and ablation-based diagnosis. These methods are independent of any specific simulation model and can be applied to both pre-deployment assessment with analytical environments and post-mission diagnosis with observed data.

SCOPE differs from conventional Observing System Simulation Experiments (OSSEs) in two respects. First, whereas OSSEs evaluate a proposed observing system’s impact on forecast skill by assimilating synthetic observations from a single Nature Run, SCOPE does not involve data assimilation. Instead, it samples uncertain environmental parameters via LHS and evaluates how observation quality varies across an ensemble of eddy realizations, yielding a distribution of outcomes that quantifies robustness rather than a single impact score. Second, OSSEs typically assess performance through a single forecast error metric, whereas SCOPE employs a 27-metric system spanning seven quality dimensions, providing a multi-faceted characterization of observation quality. In addition, while Lermusiaux (2007) established a general framework for adaptive sampling and uncertainty quantification in ocean observation, their focus was on data assimilation and adaptive path planning; SCOPE complements this by providing a diagnostic layer that decomposes performance into metric-level contributions and identifies dominant uncertainty drivers.

In the analytical model, the two-stage sensitivity analysis identifies physical couplings among parameters and reveals the non-monotonic dominant effect of ${\eta }_v$ through derived variable construction. $U_c^{\perp }$ ranks second among environmental factors ($|{\rho }_S|=0.163$, MI = 0.037), indicating that lateral drift is the dominant pathway through which background currents affect observation quality. This result is specific to the analytical model; in real eddies, the background current is spatially varying and cannot be decomposed into a single direction, so $U_c^{\perp }$ is not separately evaluated. The majority of parameter interactions (93%) are negligible, indicating approximate additivity once key couplings are resolved. The dose-response analysis shows that the optimal ${\eta }_v$ decreases monotonically with $V_{max}$: weak eddies require higher ${\eta }_v$, and strong eddies require lower ${\eta }_v$.

The ablation experiment validates the framework’s diagnostic capability by decomposing marginal effects, detecting negative synergy, and distinguishing mean improvement from variance reduction. The ablation is performed exclusively on the analytical model because isolating causal effects requires a controlled environment with fixed environmental parameters and design rules as the only varying factors. For real eddies, the stratified dose-response analysis serves a parallel role: the natural variation in eddy intensity across the four HYCOM cases reproduces the monotonic ${\eta }_v^{\mathrm{opt}}$–$V_{max}$ trend ($r=-0.96$), confirming that the intensity dependence of optimal ${\eta }_v$ is observable in realistic ocean fields.

The transferability of the framework’s methodology is tested on four HYCOM eddies in the South China Sea. The ACMS algorithm converges on all four eddies (${\rho }_S>0.95$), and the parameter importance rankings are statistically consistent across environments (Kendall $W=0.59$, $p=0.037$). The decreasing trend of optimal ${\eta }_v$ with $V_{max}$ is reproduced in the real eddies (Pearson $r=-0.96$), matching the analytical model ($r=-0.87$). However, the real-eddy optimal ${\eta }_v$ values are systematically higher than the analytical predictions, particularly for weak eddies (2.35 vs. 0.58). This offset is attributed to two mechanisms. First, the irregular, non-axisymmetric velocity field of real eddies induces larger trajectory deviations than the smooth Lamb–Oseen vortex; in weak eddies, higher glider speed produces a wave-like trajectory that increases sampling density, while in strong eddies lower speed allows the glider to move with the flow along a spiral path, accumulating more profiles at lower energetic cost. Second, the real eddies’ non-axisymmetric geometry introduces localized advective surges that further increase the speed demand in weak-eddy conditions.

Four limitations should be noted. First, the analytical model uses an axisymmetric Lamb–Oseen vortex with a uniform background current, which does not capture the azimuthal asymmetry and complex flows of real eddies, leading to the systematic offset in optimal ${\eta }_v$. Second, a continuous function ${\eta }_v^{\mathrm{opt}}\left(V_{max}\right)$ has not been derived; only discrete intensity-group trends are available, and calibration across a denser eddy intensity spectrum is needed. Third, the validation uses four eddies from a single region (South China Sea); broader geographic applicability remains to be tested. Fourth, the 500-sample ensemble per real eddy, while adequate for rank-correlation estimation, may yield conservative MI estimates for weak dependencies; larger ensembles would improve the resolution of secondary parameter effects.

5.2. Conclusions

1.

SCOPE framework. A 27-metric evaluation system spanning seven quality dimensions is constructed. The ACMS algorithm compresses the 27 metrics to 9 core indicators (${\rho }_S=0.967$) while preserving category coverage and over 95% of ranking information. The framework supports both pre-deployment assessment and post-mission diagnosis.

2.

Sensitivity hierarchy and ablation diagnosis. Two-stage sensitivity analysis identifies ${\eta }_v$ as the dominant parameter (MI = 0.375) with a non-monotonic effect; $U_c^{\perp }$ emerges as the leading environmental factor. The majority of parameter interactions (93%) are negligible. The optimal ${\eta }_v$ decreases monotonically with $V_{max}$ in both the analytical model (0.58→0.47, $r=-0.87$) and real eddies (2.35→0.35, $r=-0.96$). The ablation experiment quantifies uncertainty costs, decomposes marginal effects, detects negative synergy, and characterizes variance reduction.

3.

Methodology transfer. The ACMS algorithm and parameter ranking methodology transfer to real eddy environments. ACMS converges on all four HYCOM eddies (${\rho }_S>0.95$), with vertical profile metrics universally retained. Parameter importance rankings are statistically consistent across environments (Kendall $W=0.59$, $p=0.037$), with ${\eta }_v$ consistently ranked first.