Uncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts
UUncertainty Quantification of Tra jectory ClusteringApplied to Ocean Ensemble Forecasts
Guilherme S. Vieira , ∗ , Irina I. Rypina , † , and Michael R. Allshouse , ‡ Department of Mechanical and Industrial Engineering, Northeastern University,Boston, MA 02115, USA Department of Physical Oceanography, Woods Hole Oceanographic Institution,Woods Hole, MA 02543, USASubmitted to
Fluids on August 21, 2020
Abstract
Partitioning ocean flows into regions dynamically distinct from their surroundings based on materialtransport can assist search-and-rescue planning by reducing the search domain. The spectral clusteringmethod partitions the domain by identifying fluid particle trajectories that are similar. The partitioningvalidity depends on the accuracy of the ocean forecasting, which is subject to several sources of uncer-tainty: model initialization, limited knowledge of the physical processes, boundary conditions, and forcingterms. Instead of a single model output, multiple realizations are produced spanning a range of potentialoutcomes, and trajectory clustering is used to identify robust features and quantify the uncertainty ofthe ensemble-averaged results. First, ensemble statistics are used to investigate the cluster sensitivity tothe spectral clustering method free-parameters and the forecast parameters for the analytic Bickley jet,a geostrophic flow model. Then, we analyze an operational coastal ocean ensemble forecast and com-pare the clustering results to drifter trajectories south of Martha’s Vineyard. This approach accuratelyidentifies regions of low uncertainty where drifters released within a cluster remain there throughoutthe window of analysis. Drifters released in regions of high uncertainty tend to either enter neighboringclusters or deviate from all predicted outcomes.
Fluid flows, even if unsteady and aperiodic, may admit persistent patterns generally referred to as coherentstructures that reveal flow characteristics related to the transport of fluid particles [McWilliams, 1984;Provenzale, 1999; Haller, 2015]. Coherent structures of the elliptic type [Haller and Beron-Vera, 2012;Froyland and Padberg-Gehle, 2015; Allshouse and Peacock, 2015] are portions of fluid that do not significantlymix with the rest of the domain. From a Lagrangian perspective, the perimeter delimiting material withinthese structures remains nearly constant as they move [Allshouse and Thiffeault, 2012; Froyland, 2015], andfluid can be transported over long distances while surrounded by more vigorous mixing [Provenzale, 1999;Beron-Vera et al., 2013; Haller and Beron-Vera, 2013]. Eulerian methods compute coherent structures directlyfrom instantaneous velocity fields and may not highlight features that have a persistent impact on transport.Additionally, these methods require full velocity fields, which may be unavailable or hard to reconstruct fromsparse sets of observed trajectories. Lagrangian methods, in turn, can identify coherent structures based onthe trajectories themselves and identify the dominant features over a given time interval [Provenzale, 1999;Boffetta et al., 2001; Peacock and Dabiri, 2010; Allshouse and Thiffeault, 2012; Haller and Beron-Vera, 2013].One approach for identifying these coherent structures uses cluster analysis. Clustering algorithms revealunderlying structures in data sets by partitioning the data so that similar elements are assigned to the samecluster, while dissimilar elements are assigned to different ones. These algorithms have been extensivelystudied and widely applied in image segmentation and anomaly detection, as well as biological and physicalprocesses [Jain et al., 1999; Everitt et al., 2011]. For fluid transport analysis, clustering techniques can ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ phy s i c s . a o - ph ] A ug . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 efficiently identify elliptic structures when applied to fluid particle trajectories [Froyland and Padberg-Gehle,2015; Hadjighasem et al., 2016], which we refer to as trajectory clustering.Here, we analyze particle trajectories using the spectral clustering algorithm [Shi and Malik, 2000; Nget al., 2002; von Luxburg, 2007; Filippone et al., 2008]. This method has been used to identify coherentstructures in analytic and simulated flows [Hadjighasem et al., 2016, 2017; Vieira and Allshouse, 2020]. Theclustering performs a systematic partitioning of the trajectories into coherent and incoherent sets, providing aconceptual simplification of the underlying dynamical system for a general flow. This method is also frame-invariant, and hence the identified structures are the same in all frames that translate or rotate relativeto each other. One drawback of the spectral clustering method, however, is that there are a number offree-parameters that the user has to select, which could impact the results [Hadjighasem et al., 2017].While ocean drifter data can be clustered a posteriori to identify coherent structures, an a priori analysisto predict coherent structures relies on trajectories obtained from a forecast model. The clustering accuracyis limited by that of the trajectories, which, in turn, is limited by the accuracy of the velocity model. Inocean modeling, however, several sources of uncertainties pose a challenge to the Lagrangian approach [Ler-musiaux, 2006; Lermusiaux et al., 2006]. To simplify ocean models and reduce computational expenses, thegoverning equations are only resolved on a restricted range of spatial and temporal scales, and the influenceof scales outside this window is either parameterized or neglected. Uncertainties also arise from the limitedknowledge of processes within the scale window, which require approximate representations or parameter-izations. Moreover, measurements used for model initialization and parameter estimation are limited incoverage and accuracy, leading to imprecise initial conditions and model parameters. Finally, models ofinteractions between the ocean and the atmosphere are approximate, and ocean boundary conditions areinexact. All of the above lead to differences between the actual values and the modeled values of physicalfields and properties.To account for model uncertainties intrinsic to the modeling process, an effective option is to performensemble statistics. Different sets of model parameter values generate an ensemble of possible outcomes,which are then processed to provide probabilistic information about the variability on the end results. Search-and-rescue planning, for instance, already considers ensemble statistics to produce probability-distributionmaps for the target’s location [Kratzke et al., 2010]. The vast uncertain parameter space together with thecontinuous motion of floating objects driven by unsteady flows, however, can lead to error accumulation in thepredicted trajectories [Serra et al., 2020]. Coherent structures have been shown to depend less on individualtrajectories and are more robust to model parameter variations and noise, highlighting the main structuresin the flow even in the event of imperfect or scarce trajectories [Haller, 2002; Lermusiaux et al., 2006]. Therobustness of the clustering algorithm to perturbations in the individual trajectories via an ensemble set ofrealizations has yet to be studied and could aid in quantifying uncertainty for the trajectory clusters.The detection of robust clusters could then be used in the deployment of drifters to observe such featuresin the ocean. Other Lagrangian methods have been used in the past in experiments to interpret the observedbehavior of drifters deployed in the ocean [Olascoaga et al., 2013; Jacobs et al., 2014; Beron-Vera et al., 2015;Williams et al., 2015; Rypina et al., 2014, 2016; Rypina and Pratt, 2017]. In these studies, drifters werereleased, tracked, and their trajectories were compared to results from Lagrangian analyses performed aposteriori using velocities from models, satellite altimetry or radar measurements. Few studies, however,have used Lagrangian methods to plan and execute field experiments [Haza et al., 2007, 2010; Serra et al.,2020]. The drifter-release experiment presented in [Filippi et al., 2020], whose drifter data is used in thispaper, is one of the first experiments targeting coherent structures based on an a priori trajectory clusteringobtained from forecast simulations.Our approach accounts for and quantifies uncertainty when clustering trajectories. We apply the spec-tral clustering algorithm varying the method free-parameters to understand how the clustering results aresensitive to the implementation, analyze ensemble simulations to understand how the model parametersimpact the resulting clusters, and finally apply the method to a forecast data set to compare the clusteringprediction to experimental drifter data. The method by Hadjighasem et al. [2016] is modified with a morebroadly used similarity function and soft clustering membership probabilities, allowing for a probabilisticview of the clusters for each individual model realization. Two different systems are analyzed with ourapproach: an analytic flow model and forecast simulations provided by a coastal ocean model. First, theBickley jet analytic flow [Rypina et al., 2007; Beron-Vera et al., 2010] with model parameter variations isused to mimic model uncertainty. Then, we analyze ensemble-forecasts generated by the MIT-MSEAS oceanmodel [Haley and Lermusiaux, 2010; Haley Jr et al., 2015] of the coastal region of the Martha’s Vineyard2 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts island. Ensemble statistics of the resulting clusters provide a probabilistic view of the coherent structuresidentified by the method. The forecast clustering results are used to identify coherent structures in a drifter-release experiment, and the drifter trajectories are compared to the forecast cluster behavior and associateduncertainty.The paper is organized as follows. Section 2 presents the spectral clustering method used in this work andhow clustering results are processed to provide ensemble statistics. Section 3 introduces the quasi-periodicBickley jet system, performs ensemble statistics for a single set of particle trajectories while varying themethod free-parameters to quantify the sensitivity of the clustering results to the implementation. Section 4varies Bickley jet system model parameters to asses the robustness of the identified coherent structures in ascenario of uncertainty in the velocity field. Section 5 presents the clustering analysis used to target coherentstructures in a drifter-release experiment south of Martha’s Vineyard, and compares the forecast results tothe observed drifter trajectories. Finally, conclusions are presented in Section 6. This section presents the spectral clustering method and how uncertainty is measured. We adapt the methodof Hadjighasem et al. [2016] to use soft clustering (fuzzy c -means) to assign a cluster membership probabilityto each particle trajectory. Section 2.1 introduces the method along with specifics of the similarity measureselected and the soft clustering approach. Section 2.2 presents the general procedure we use to calculatestatistics and quantify uncertainty when considering an ensemble of possible clustering results. To analyze flows from a Lagrangian perspective, we consider massless fluid particles that move with thevelocity field u ( t, x ). The trajectory x i ( t ) of fluid particle i departing from x i ( t ) at time t = t is x i ( t ) = x i ( t ) + (cid:90) tt u ( τ, x i ( τ )) dτ. (1)A total of N particles are initialized in a grid that uniformly covers the targeted domain at t , and theirindividual trajectories are numerically integrated and output at discrete time instances within the timeinterval [ t , t f ].We use the set of trajectories { x i ( t ) } ≤ i ≤ N to partition the spatial domain into clusters. The spectralclustering algorithm performs an eigenanalysis to project the trajectory set onto a subspace that may yieldclusters maximizing the within-cluster similarity and minimizing the between-clusters similarity [Nock et al.,2009]. Particles clustered together should move as a compact group, with limited mixing with particlesoutside of the cluster, while particles in different clusters should experience dissimilar trajectories [Froylandand Padberg-Gehle, 2015; Hadjighasem et al., 2016].The spectral analysis requires the construction of a positive, weighted, undirected graph known as thesimilarity graph. This graph quantifies the pairwise similarities between trajectories. Each graph noderepresents a particle trajectory, and the graph edges are weighted by the similarity between the nodes theyconnect [von Luxburg, 2007; Hadjighasem et al., 2016; Shi and Malik, 2000]. To compute these weights,one must first define a measure for similarity between trajectories. One possible metric for the dissimilaritybetween trajectories x i and x j is their time-averaged distance r ij = 1 t f − t (cid:90) t f t dist( x i ( τ ) , x j ( τ )) dτ, (2)where dist( · , · ) is the Euclidean distance. Then, a similarity measure for the edge weights w ij = w ( r ij ) mustbe chosen as a function of the time-averaged distance. The only restriction on this functional dependenceis that the weight should be a monotonically decreasing function of the distance. This choice of similarityfunction controls the graph edge weight distribution, and therefore can impact the clustering results [vonLuxburg, 2007].Hadjighasem et al. [2016] use a similarity function that is the inverse of the time-averaged distance, w ij = l x /r ij (a constant l x is included here to make the weight dimensionless and does not impact the3 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 results). A more widely used function for graph partitioning is the Gaussian function [Shi and Malik, 2000;Ng et al., 2002; von Luxburg, 2007] w ij = exp (cid:0) − r ij / σ (cid:1) , (3)which we will use. The similarity radius σ > w ii = 1 and w ij → r ij (cid:29) σ . Thechoice of the Gaussian measure has a number of advantages over the l x /r ij choice: first, it is bounded, with0 ≤ w ij ≤
1; second, no offset value is set in the diagonal entries [Hadjighasem et al., 2016]; and third,the sparsification step can be skipped (or at least has a negligible impact on the results, see SupplementaryMaterial A), as (3) goes to zero faster than the measure in [Hadjighasem et al., 2016] for large r ij . Thesimilarity radius σ in (3) determines the rate of decay of w ij to zero as r ij increases, which is analogous tothe graph sparsification in [Hadjighasem et al., 2016] where w ij = 0 for r ij above a defined threshold. Adirect comparison between these two similarity functions and the impact of the σ -choice are presented inSections 3.1 and 3.2, respectively.The similarity matrix W ∈ R N × N stores the w ij values, and the diagonal degree matrix D is computedsuch that d ii = (cid:80) Nj =1 w ij . The generalized eigenvectors q of the unnormalized graph Laplacian L = D − W are computed from the generalized eigenproblem Lq = λ Dq . (4)The normalized vectors q , q , . . . , q N , corresponding to the generalized eigenvalues 0 ≤ λ ≤ . . . ≤ λ N , dif-ferentiate properties in the graph and facilitate the clustering process [von Luxburg, 2007]. While all ofthe eigenvectors provide information, the dominant eigenvectors, associated with the smallest eigenvalues,reveal the most dynamically relevant characteristics. Each trajectory is characterized by a value within eacheigenvector, and these are the values ultimately used to cluster trajectories.The eigenvector matrix Q ∈ R N × M whose columns are the dominant eigenvectors q , . . . , q M , with M (cid:28) N , stores the characteristics used to cluster the trajectories. The number of retained eigenvectors M is specified depending on the system. Let y i ∈ R M be the characterization vector corresponding to the i -th row of Q , which contains the condensed differentiating information of trajectory i . The eigenanalysisprovides a suitable low dimensional representation of the data set. The characterization vectors { y i } ≤ i ≤ N are then partitioned into K clusters. The relationship between the number of clusters K and the number ofdominant eigenvectors M used for clustering depends on characteristics of the system, and is specified foreach of the case studies in the following sections.We cluster the characterization vectors using a fuzzy c -means algorithm [Froyland and Padberg-Gehle,2015; Filippone et al., 2008], instead of the conventional k -means often performed at this step [von Luxburg,2007; Hadjighasem et al., 2016]. The c -means algorithm assigns to each trajectory i probabilities p i,k ∈ [0 ,
1] of being a member of cluster k , with 1 ≤ k ≤ K , such that (cid:80) Kk =1 p i,k = 1. This step requires theprescription of a fuzziness parameter m > m →
1, the clustering result approaches the k -means result where p i,k ∈ { , } . For greater m , the “looser” distribution in membership probabilities allows the identification of particles that presentintermediate behaviors between clusters, as opposed to always assigning a trajectory as member of a singlecluster. The cluster centers are initialized based on the dominant eigenvectors before starting the c -meansiterative process detailed by Froyland and Padberg-Gehle [2015]. Finally, Hadjighasem et al. [2016] use aneigengap heuristic to determine M and K that relies on a large cluster differentiation to work properly,and in many practical cases an eigengap may be less pronounced or nonexistent [von Luxburg, 2007]. Weleave the number of dominant eigenvectors and clusters as method free-parameters and study the relateduncertainty associated to the clustering results for each K -choice. In Sections 4 and 5, we will be interested not in applying the spectral clustering algorithm to a singlerealization of the flow, but in collecting clustering results from different realizations and combining them toget statistical information about the variability of the clusters with the method free-parameters σ , m and K ,or with velocity model parameters. We therefore need a method to combine clustering results from differentrealizations. Regardless of the differences between realizations, we initialize the N particles on the same grid,so that trajectory i is uniquely labelled across realizations, based on x i ( t ). Each realization I ∈ { , . . . , R } ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts generates a full set of membership probabilities p ( I ) i,k , with i ∈ { , . . . , N } and k ∈ { , . . . , K } . The clusterlabels for different realizations are matched based on the similarity between cluster centers.To quantify the uncertainty of a trajectory i being a member of cluster k , the probabilities p ( I ) i,k are usedto compute the mean membership probabilities over all R realizations p i,k = 1 R (cid:32) R (cid:88) I =1 p ( I ) i,k (cid:33) , (5)and the corresponding sample standard deviations S i,k = (cid:118)(cid:117)(cid:117)(cid:116) R − R (cid:88) I =1 (cid:16) p ( I ) i,k − p i,k (cid:17) . (6)Both the mean and the standard deviation of the realizations are bounded values, with 0 ≤ p i,k ≤ ≤ S i,k ≤ .
5. We perform this calculation for each trajectory, for each cluster.
We analyze the analytic, quasi-periodic Bickley jet system to evaluate the spectral clustering method free-parameter and velocity field model sensitivity. This system is an idealized model of a meandering zonal jetunder geostrophic balance and has been extensively used to illustrate coherent structures in fluid flows [Ryp-ina et al., 2007; Beron-Vera et al., 2010; Schlueter-Kuck and Dabiri, 2017; Hadjighasem et al., 2017]. Themodel features a sheared zonal flow on which a superposition of Rossby-like waves propagate. The stream-function ψ prescribing the two-dimensional velocity field u = − ∂ y ψ e x + ∂ x ψ e y with a superposition of threewaves is ψ ( t, x, y ) = − U L tanh (cid:16) yL (cid:17) + U L sech (cid:16) yL (cid:17) (cid:88) n =1 A n cos [ k n ( x − c n t + φ n )] , (7)where U and L are a characteristic speed and length, respectively. The domain is periodic in the x -direction,with periodicity l x = 2 πR e cos(60 ◦ ), where R e = 6371 km is Earth’s radius. The rectangular domaincorresponds to x/l x ∈ [0 , y/l x ∈ [ − . , . A n , wave numbers k n = 2 πn/l x , phase speeds c n , and phases φ n , for n = 1 , ,
3. To model the self-consistent state obtainedby Rypina et al. [2007] for modes 2 and 3 on the periodic domain, we fix U = 62 .
74 ms − , L = 1767 km, c /U = 0 . c /U = 0 . c /U = 0 . A (cid:28) min( A , A ),the mode-1 wave acts as a small perturbation to the system’s periodicity. In this section, we fix the modeamplitudes to A = 0 . A = 0 .
15, and A = 0 .
30, values used by Hadjighasem et al. [2016], and allthree waves are in phase: φ = φ = φ = 0. Note, however, that while the system dynamics depend on thevalues of A n and φ n [Rypina et al., 2007], there is no physical basis for the stated choice of amplitudes andphases, and the impact of varying these parameters will be explored in Section 4.The present section discusses the application of the spectral clustering algorithm to the Bickley jet system,and studies the sensitivity of the results to user-defined method free-parameters. We present the impact ofthe choice of similarity measure in Section 3.1. Then, study the sensitivity to the method free-parameters:similarity radius σ , in Section 3.2, tightness of the cluster memberships m , in Section 3.3, and number ofclusters K , in Section 3.4. To cluster the Bickley jet system according to the method described in Section 2.1, particles are initializedin a uniform grid of 400 by 120 positions uniformly covering the domain, and advected from t = 0 to t f = 40 days, matching [Hadjighasem et al., 2016]. The distance function in (2) takes into consideration the x -periodicity of the domain, and the M = 6 dominant eigenvectors of (4) are used to partition the domaininto K = 7 clusters, to account for 6 materially coherent vortices, which are the coherent clusters, and anincoherent cluster, the chaotic sea [Hadjighasem et al., 2016]. The method uses a fuzziness parameter m = 2.5 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 Figure 1: Spectral clustering membership probabilities for clusters k ∈ { , . . . , } identifying materiallycoherent vortices, with fuzziness m = 2 and number of clusters K = 7 (incoherent cluster probabilitiesare omitted). The similarity functions w ij used are ( a ) the proposed Gaussian similarity measure (3) with σ/l x = 0 . b ) the inverse measure l x /r ij from [Hadjighasem et al., 2016] sparsifying values less than1/0.075. The six color maps presented here are also used in Figures 2, 4(a,b) and 6.The membership probabilities for the clusters identifying the 6 materially coherent vortices at time t are presented in Figure 1. Figure 1(a) presents our results, obtained using a Gaussian similarity function(3), with similarity radius σ/l x = 0 . t due to the x -periodicity. The membership probabilitiesof being part of the chaotic sea, the seventh cluster, are complementary to the ones plotted and are omittedthroughout. Note that the use of a soft membership assigns to particles located at the periphery of thevortices lower probabilities of being a member, which relates to lower similarity in the dynamics (some ofthem may, for instance, be trapped inside the vortices for just a fraction of the time window of analysis,then merge with the chaotic sea, or vice-versa).The results for the Gaussian similarity measure in Figure 1(a) are similar to those obtained with the l x /r ij similarity measure used by Hadjighasem et al. [2016], presented in Figure 1(b). For the latter, matrixentries w ij = 0 for r ij /l x > . σ for the similarity measure (3). Both similarity functions yield similar results for the selectedparameter values, with smoother transitions in membership probabilities (from 1 to 0) for the Gaussianmeasure.When using the l x /r ij similarity measure, the degree of sparsification is an additional parameter for theclustering method that can impact the results. This parameter can be eliminated for the Gaussian similaritymeasure as no sparsification was necessary for the result in Figure 1(a), but it is worthwhile to sparsifythe matrix to reduce computational costs and storage. We demonstrate in the Supplementary MaterialA that sparsifying entries of W that satisfy w ij < exp( − / ≈ · − , corresponding to r ij > σ , hasnegligible impact on the clustering membership probability results, and hence we sparsify according to thisrule hereafter. Note that this result is valid for the Bickley jet, but the impact of sparsification may vary foran arbitrary flow. We highlighted in the previous section how the similarity radius σ is closely related to the sparsificationthreshold used by Hadjighasem et al. [2016]. While the sparsification threshold selection was mentionedin [Hadjighasem et al., 2016], the clustering results are highly sensitive to this parameter [Filippi et al.,2020], and we now address this sensitivity.While there may be some intuition on the size of the structures of interest, this is not always helpful inprescribing σ or in understanding how the σ -choice impacts the resulting clusters. Hadjighasem et al. [2016]choose their threshold by defining which values of W to keep based on a fixed percent sparsification of the6 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts Figure 2: Membership probabilities for vortices k ∈ { , . . . , } , for different choices of the parameter σ ,plotted at t = t . Values correspond to ( a ) σ/l x = 0 . b ) 0 . c ) 0 .
030 and ( d ) 0 . σ in the Gaussian similarity function (3), we vary σ to demonstrating the clustering sensitivity to changes insparsification. First, we define an interval bounding all relevant σ -values to be tested. For σ/l x distributedon the interval [0 . , . m = 2 to identify K = 7 clusters. Figure 2 presents the membership probabilities for each of the six coherent clusters for σ/l x = 0 . σ/l x = 0 .
020 presented in Figure 1(a),smaller values of σ (Figures 2(a,b)) tend to shrink the coherent clusters to the vortex cores only, while largerchoices of σ (Figures 2(c,d)) assign higher membership probabilities to filaments that correspond to particlesthat do not belong to the coherent vortex from the start, but have a long residence time on the perimeter ofthe vortex.To determine how the membership probabilities for each trajectory depend on σ , we use the informa-tion from the different realizations to compute the membership probability means p i,k and sample standarddeviations S i,k for each trajectory and cluster, as described in Section 2.2. These statistical measures arepresented in Figure 3. In Figure 3(a), we present p i, for vortex k = 1, and in Figure 3(b) the correspond-ing standard deviation S i, . Figures 3(c,d) condense the information for all vortices k ∈ { , . . . , } . Thesuperimposed mean p i, ˆ k in Figure 3(c) and the superimposed standard deviation S i, ˆ k in Figure 3(d) aresuch that, for each trajectory i , ˆ k is the cluster that maximizes the mean membership probability, henceˆ k = argmax k p i,k .Figure 3(a,c) demonstrates that the vortex cores have the greatest mean membership, which relates tothe fact that those particles are identified with high membership probabilities in all realizations. Particlesfurther away from the cores have a lower mean, as a result of lower probabilities of being part of the respectivevortices, in particular for low σ values. We also notice sharp drops in the probability after a certain vortexsize. The uncertainty of particles being part of a vortex is highlighted by the standard deviations in Figure3(b,d), where we again see negligible standard deviation (low uncertainty) on the membership probabilities forthe cores. Because there are no realizations that identify the central jet and some portions of the chaotic sea7 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 Figure 3: Clustering statistics varying σ , plotted at t = t , for σ uniformly distributed in [0 . , . a ) mean for k = 1, ( b ) standard deviation for k = 1, ( c ) superimposed mean and( d ) superimposed standard deviation for k ∈ { , . . . , } . The two color maps presented here are also used inFigures 4(c,d), 5, 7 and 8(a).as part of the clusters, there is also a low standard deviation for those trajectories. The highest uncertaintyis obtained for particles between the core and the chaotic sea, and some filaments are also highlighted withhigher standard deviation. Those filaments correspond to particles consistently identified as part of a vortex,with high membership probabilities, for σ greater than a threshold, but not for lower values of σ . Next, we consider the sensitivity of the choice of the fuzziness parameter m , introduced by the c -meansstep that assigns cluster membership probabilities. In this analysis, σ/l x = 0 .
020 while m is varied on theinterval [1 . , .
00] with steps of 0.05. The results are presented in Figure 4. The membership probabilitiesfor the extreme values ( m = 1 .
05 and 3.00) are illustrated in Figure 4(a,b), and the superimposed mean andstandard deviation in Figure 4(c,d).While m = 1 .
05 (in Figure 4(a)) yields a bimodal probability distribution, with 99% of the p i,k > . .
95 (so tending to the corresponding k -means result), the use of m = 3 .
00 (Figure4(b)) produces smoother probability transitions from the vortex core to the perimeter, as well as overalllower membership probabilities, with only 10% of the p i,k > . .
95. The superimposedmean and standard deviations reveal a similar trend to the one observed for the dependence on σ , with twomajor differences: (i) while varying m introduces uncertainty in the membership probability of trajectoriesstarting at the vortex perimeters, the m -choice does not lead to the same significant capture of filaments asmember of the vortex for the current σ/l x = 0 .
020 value, and (ii) the magnitude of the standard deviationsrelated to m are about half of the ones associated to σ (see Figure 4(d) in comparison to Figure 3(d)). Finally, the number of clusters in the system is not necessarily known beforehand, and here we address howthe clustering results statistically vary for different choices of K . As presented in [Hadjighasem et al., 2016],for a suitable sparsification level of the similarity matrix, the eigengap heuristic can be used to infer the use of M = 6 dominant eigenvectors to choose K = 7 for the Bickley jet, and thus identify the 6 materially coherent8 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts Figure 4: Membership probabilities for k ∈ { , . . . , } , with ( a ) m = 1 .
05, and ( b ) m = 3 .
00. Clusteringstatistics plotted at t = t for m uniformly distributed in [1 . , . c ) superim-posed mean and ( d ) superimposed standard deviation for k ∈ { , . . . , } . Corresponding colorbars presentedin Figures 1 and 3.vortices, plus the chaotic sea. It is known, however, that the existence of an eigengap is not guaranteed forany system, and depends on, among other properties, the connectivity of the similarity graph [von Luxburg,2007]. The number of clusters for systems not as distinctly separated as the Bickley jet will, therefore, bemore challenging to determine. Without any prior knowledge about the number of clusters in the system,one could consider clustering based on different numbers of dominant eigenvectors of (4) into different K ,which might result in merging clusters, splitting clusters, missing clusters, or even identifying new ones. Aswe vary K , we fix the relationship between M and K to K = M + 1 to always include the chaotic sea as anindependent cluster [Hadjighasem et al., 2016].Because the free-parameter K cannot be varied continuously, and because varying K while fixing σ and m would only mean changing the number of eigenvectors to use in the c -means step, we adopt a differentstrategy to quantify the K -uncertainty. For each choice of K , we perform statistics using realizations in which σ and m are varied. The ( σ/l x , m ) pairs are drawn from uniform distributions over [0 . , . × [1 . , . K = 6, 7 and 8. One should now expectmissing or merging vortices for K <
7, while splitting vortices or identifying new structures for
K >
7. Forthe ensemble statistics, if vortex k is not identified in realization I , we set p ( I ) i,k = 0 for all i . For K = 8, anew coherent cluster corresponding to the jet is consistently identified, in all realizations. The jet is thereforeconsidered a seventh coherent cluster.The superimposed means and standard deviations for K = 6, 7 and 8 are presented in Figure 5. For K = 6, Figure 5(a) reveals that different ( σ/l x , m ) free-parameter choices can result in different vorticesnot being identified. Notably, vortices 3, 5 and 6 are missed in some realizations, which reduces theirmean probability and increases uncertainty for those clusters. Figure 5(b), for K = 7, presents six vorticesidentified with similar probability distributions to the ones obtained by varying σ alone in Section 3.1 (inFigure 3(c,d)). Such similarity between these cases highlights the dominance of the σ -sensitivity over thechoice of m , and the fact that these two parameters do not compound each other. Finally, Figure 5(c)presents the results for K = 8, highlighting the consistent identification of the jet as a coherent cluster in thesystem. The detection of the jet has little impact on the probability distributions for the six vortices. Boththe means and standard deviations for the vortices remain similar to the ones obtained for K = 7, with a9 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 Figure 5: Clustering statistics varying σ and m , for different numbers of clusters K , with ( σ/l x , m ) drawnfrom a uniform distribution over [0 . , . × [1 . , . p i, ˆ k (left) and the corresponding superimposed standard deviation S i, ˆ k (right), for ( a ) K = 6, ( b ) K = 7 and( c ) K = 8. Corresponding colorbars presented in Figure 3.slight standard deviation increase for the cores. The jet cluster has an overall higher sensitivity to σ and m ,with higher standard deviations than the vortex cores in Figure 5(b).While these ensemble statistics could be used as a basis for setting the method free-parameters, a thoroughinvestigation on this is beyond the focus of the parameter sensitivity study presented in this paper. Basedon the results for K = 6, 7, and 8, it would be reasonable to state that K = 7 could be chosen over the otheroptions, as it is the case leading to the smallest space-averaged uncertainty. At the same time, however, thejet is a structure that differentiates itself from the rest of the flow, by behaving in a distinct and coherentway compared to the remaining particles in the chaotic sea. It could be argued that the choice of K = 8 forthis system is an equally possible way of clustering the system, and provides extra information about thejet cluster, at the price of increasing the result sensitivity to other method free-parameters. Further work isnecessary to automate this selection for a general system. In the previous section, a single set of model parameters of the Bickley jet system, with fixed physicalparameters, was used to demonstrate how the choice of the method free-parameters impact the clustering10 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts results. Here, we focus on model parameters that influence the dynamics of the system by changing thevelocity field and the resulting trajectories. Multiple realizations of the system are used to determine theclustering sensitivity to model parameters, allowing for a characterization of the robustness of the identifiedclusters. Section 4.1 explains how system parameters are randomized to generate multiple dynamically differ-ent realizations. Section 4.2 presents the statistics over the described realizations, leading to an uncertaintyquantification of the clustering results to model parameters.
While system parameters such as U , L , c , and c are set by physical arguments (as discussed in Section 3),the wave amplitudes A n and phases φ n are not, despite exerting a major influence on the system dynam-ics [Rypina et al., 2007]. We use these values as unknown model parameters to introduce variability in thesystem dynamics. The amplitudes and phases are drawn from normal distributions centered around thevalues used in Section 3 (and [Hadjighasem et al., 2016]). The realization presented in Section 3 is hereafterreferred to as the central realization, and corresponds to amplitudes A n = A n , with A = 0 . A = 0 . A = 0 .
30, and phases φ n = 0. Therefore, each realization I in this section is generated by A ( I ) n / A n ∼ N (cid:32) , (cid:18) (cid:19) (cid:33) and φ ( I ) n / l x ∼ N (cid:32) , (cid:18) (cid:19) (cid:33) , (8)where N (cid:0) µ, Σ (cid:1) denotes a normal distribution of mean µ and standard deviation Σ. The standard deviationsfor A ( I ) n are scaled by the corresponding mean values, while the standard deviation for φ ( I ) n is chosen smallenough so that the vortex centers for each realization are likely to be inside of the area covered by the vorticesin the central realization.A total of R = 1000 realizations are generated from these distributions, and the spectral clusteringalgorithm described in Section 2 is applied with method parameters fixed to σ/l x = 0 . m = 2 and K = 7. While it is possible that individual realizations may require a different method free-parameterselection, our goal in presenting this section is to separate the effects of method free-parameters from thoseof model parameters. By sampling model parameters together with method free-parameters, the clusteruncertainty resulting from the model would be obfuscated. We thus fix the method free-parameters basedon the central realization only, while considering multiple realizations with varying model parameters.Figure 6 presents the membership probabilities for four realizations. While the realizations do modifythe position, shape, and dynamics of the vortices (see videos in Supplementary Material), their presenceand number, as well as the presence of the shear jet, mostly remain unchanged. Figure 6(a,b) presents howvariable the initial cluster sizes and shapes can be, as well as the effect of the wave phases on the nonuniformspacing between the identified vortices at t . While for most realizations all six vortices are identified, thereare cases where some of the expected vortices are not identified. Figure 6(c) presents a case in which the jetis identified and one of the vortices is missed. For that case, vortex 5 gets identified as part of the chaoticsea (not plotted), while a highly asymmetric jet is identified as another coherent cluster in the system,and trajectories are assigned membership probabilities p i,jet . Figure 6(d) presents a case for which a moresymmetric jet is identified as a cluster and two of the vortices (2 and 4) are merged into a single cluster.Other realizations, not illustrated here, result in cases where only a few (or even none) of the vortices areidentified. For those realizations, there is no clear Eulerian signature of the six vortices. In what follows,only clusters that identify one and only one vortex are considered for statistical purposes. The clustering results for all R = 1000 realizations are analyzed to measure their uncertainty statistics.For each one of the vortices, the mean and standard deviation over the ensemble of parameter value setsof the Bickley jet are computed according to (5) and (6). Figure 7 presents the superimposed membershipprobability mean, p i, ˆ k , and corresponding superimposed sample standard deviation, S i, ˆ k , for the ensemble.Even with a variety of behaviors observed in the ensemble, as illustrated in Figure 6, averaging the ensemblesmooths the vortices, resulting in cluster shapes and positions that are similar to the ones obtained for thecentral realization (Figure 1(a)). However, the mean probabilities in Figure 7(a) highlight how introducinguncertainty in the wave amplitudes and phases leads to smaller vortex cores with high membership probabil-ities. The mean probabilities now peak at 0.91 rather than 1.00, due to the phase shift between realizations.11 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 Figure 6: Examples of membership probabilities for the six identified clusters, plotted at t = t , for differentmodel parameters { A n } and { φ n } drawn from normal distributions. Cases correspond to parameters ( A , A , A , φ /l x , φ /l x , φ /l x ) equal to ( a ) (0.0087, 0.19, 0.20, 0.01, 0.06, -0.03), ( b ) (0.0069, 0.26, 0.25, 0.00,-0.08, 0.09), ( c ) (0.0102, 0.35, 0.25, 0.01, -0.01, 0.01), where the jet is identified and one of the vorticesmissed, and ( d ) (0.0077, 0.11, 0.32, 0.00, -0.01, 0.02), where the jet (grayscale) is identified and two vorticesare merged. Corresponding colorbars for the vortices presented in Figure 1. (See Supplementary Materialfor a video of the time evolution of the clustered trajectories.)The membership probability decay from the vortex cores to the perimeters is also more gradual than for thecentral realization, and the averaging clears out previously identified filaments that are realization-dependent.Figure 7(b) highlights a more spread out distribution and overall higher magnitude for the superimposedstandard deviation. Moreover, higher standard deviations are now observed at the vortex cores, which arisesfrom the phase parameter uncertainty. Also, notice that most of the low uncertainty regions associated to thejet in Figure 3(b) have higher uncertainty in Figure 7(b). With the varying model parameters and dynamics,the vortex positions, shapes and trajectories are more variable. All of these introduce uncertainties that arenot observed for the central set of parameter values. The ensemble analysis, therefore, highlights structuresensitivities (and robustness) that are not apparent from the central realization alone.To illustrate the ramifications of regions of high uncertainty, we demonstrate how different ensembletrajectories are when initialized at high and low uncertainty locations. At t = 0, the orange particles inFigure 8(a) are initially located at the core of cluster 1 and correspond to p i, = 0 .
883 and S i, = 0 . p j, = 0 .
469 and S j, = 0 . k = 1 and 4 for the central realization (for particles suchthat p i,k ≥ . t = 40 days, as opposed to only 45.0% for the higher uncertaintycase. The presence of the jet separating the top and bottom of the domain keeps most orange and redparticles from moving to the opposite half.Regions of higher uncertainty, not characterized when considering the central realization only, are revealedby the ensemble analysis, and correspond to flow regions for which particle cluster membership is mostunknown. Our analysis shows that the cores of the vortices are robust even if the model parameters arevaried, but the narrow filaments identified in the central realizations should be viewed as less robust as12 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts Figure 7: Clustering statistics, plotted at t = t , for an ensemble of R = 1000 realizations of the Bickleyjet system with variable amplitudes A n and phases φ n drawn from normal distributions around the centralrealization values. Membership probability ( a ) superimposed mean p i, ˆ k and ( b ) superimposed standarddeviation S i, ˆ k , for vortices k ∈ { , . . . , } . Corresponding colorbars presented in Figure 3.Figure 8: Evolution of trajectories from the 1000 realizations released from a low (orange, top) and high (red,bottom) uncertainty position. ( a ) Orange particles are initialized at the core of vortex 1, while red particlesare initialized at the perimeters of vortex 4. Positions are plotted on top of the superimposed standarddeviation S i, ˆ k computed using vortices 1 and 4, and plotted at time t = 0. The positions for the multiplerealizations and the central realization vortex boundaries are presented at ( b ) t = 5, ( c ) t = 20 and ( d ) t f = 40 days. Transparency is used to highlight high or low concentration of particles. The vortex boundariesfrom the central realization are plotted in gray and enclose particles with membership probabilities greaterthan 0.5, for clusters 1 and 4. Corresponding colorbar for ( a ) presented in Figure 3. (See SupplementaryMaterial for a video of the time evolution.)demonstrated by the ensemble analysis. Further, while both particles in Figure 8 are initialized and remaininside of their respective clusters for the central realization, the same happens for less than half of therealizations considered, once model parameter uncertainty is incorporated and accounted for.13 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020 Having demonstrated the clustering uncertainty quantification method on an analytic system, we applythe technique to a real ensemble forecast of a nested primitive equation ocean model, with intrinsic modeluncertainty. The model is used to forecast the three-dimensional velocity field for the coastal region nearMartha’s Vineyard, an island located south of Cape Cod, Massachusetts, USA. Trajectories from driftersdeployed [Filippi et al., 2020] for the corresponding day are then compared to the cluster behaviors. Sec-tion 5.1 presents characteristics of the Martha’s Vineyard region, the model used to forecast the velocityfield, and the results of the clustering uncertainty quantification analysis applied to the ensemble forecast.Section 5.2 details the drifter experiment and compares the drifter trajectories to the forecast trajectoryclustering results and the associated uncertainty.
The island of Martha’s Vineyard, with an area of almost 250 km , is the largest island in New England,and lies 11 km off the coast of Cape Cod. The prevailing currents in the coastal region south of the islandare associated with wind-driven coastal divergence, tidal forcing and a mean southward drift [Rypina et al.,2014]. During the summer months, the region experiences a mean westward surface current that reachesvelocities of 15 cm s − [Rypina et al., 2016]. The drifter deployment experiments targeting predicted co-herent structures [Filippi et al., 2020], presented in Section 5.2, took place around the 2 . uninhabitedisland of Nomans Land, south of Martha’s Vineyard. The channel between the two islands has a width ofapproximately 5 km and an average depth of 10 m.We used the MIT Multidisciplinary Simulation, Estimation, and Assimilation Systems (MIT-MSEAS)primitive equation ocean modeling system [Haley and Lermusiaux, 2010; Haley Jr et al., 2015] to computeocean surface velocity forecasts in the Martha’s Vineyard coastal region during August 2018. The modelingsystem provided forecasts of the ocean state variable fields (three-dimensional velocity, temperature, salinity,and sea-surface height) every hour, with a spatial resolution of 200 m. More details about the model forecastinitialization, tidal forcing, atmospheric flux forcing, and CTD data assimilation can be found in [Serra et al.,2020; Filippi et al., 2020]. The deterministic two-way nesting ocean forecast initialized from the estimatedocean state conditions at a particular time is referred to as the central forecast, and ensemble forecasts wereinitialized using Error Subspace Statistical Estimation procedures [Lermusiaux, 2002]. The forecasts withinthe ensemble were initialized from perturbed initial conditions of all state variables and forced by perturbedtidal forcing, atmospheric forcing fluxes and lateral boundary conditions. These perturbations were createdto represent the expected uncertainties in each of these quantities. Finally, parameter uncertainties (bottomdrag, mixing coefficients, etc.) were also modeled by perturbing the values of parameters for each forecast.A total of 71 forecasts were considered for the present study. Fluid particle trajectories confined tothe surface were analyzed over a 6 h-time-window, between t = 16:00 and t f = 22:00 UTC on August 7,2018. Such a short timing is critical in search-and-rescue operations, as after six hours the likelihood ofrescuing people alive drops significantly [Serra et al., 2020]. The forecast velocity fields used for trajectoryintegration were generated the night before the experiment. Synthetic trajectories were computed usingthe web-based gateway Trajectory Reconstruction and Analysis for Coherent Structure Evaluation [Ameliand Shadden, 2019, 2020]. At time t , particles are uniformly distributed on a 250-by-250 grid covering thedomain [70 . ◦ W , . ◦ W] × [41 . ◦ N , . ◦ N], from which portions corresponding to land are removed.This grid is approximately 21 km by 22 km.Figure 9(a) presents the initial particle distribution, superposed with the velocity field for the centralforecast. Trajectories are integrated using an adaptive 7 th -order Runge-Kutta-Fehlberg method and bicubicspline spatial interpolation, with free-slip boundary conditions applied near land. Particle positions areoutput every 5 minutes. Figure 9(b) presents the final positions of the particles for the central forecast.Darker regions correspond to particles collecting, which is mostly observed along the coast. The modelvelocity field captures the reversal of the tide, as can be observed from the flipping in the average flowdirection between t in Figure 9(a) and t f in Figure 9(b). Trajectories obtained in each of the 71 forecastsare presented in Figure 9(c) in contrast with the central forecast, for particles initialized at distinct positions,to demonstrate the degree of trajectory variability between forecasts.A similar study to the one in Section 3 was performed for the central forecast, to determine the spectralclustering sensitivity to method free-parameters. We select the following method free-parameters: similarityradius σ = 1 km, fuzziness m = 2 and number of clusters K = 4, which are used for all R = 71 forecasts.14 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts Figure 9: Martha’s Vineyard coastal area, central forecast particle distribution, and superimposed forecastmodel velocity field at ( a ) the initial time t = 16:00, ( b ) the final time t f = 22:00 UTC. ( c ) Trajectoriesobtained in each of the 71 forecasts for 4 different initial positions. Circles represent the initial positions andcrosses the final positions, after 6 hours. The darker trajectories correspond to the central forecast.While it is possible to break the domain into fewer clusters, four clusters minimize the space-averagedstandard deviation of the results. Because the graph is fully connected with this choice of σ , λ = 0 andthe components of q are constant [von Luxburg, 2007]. Therefore, the clustering (in the absence of achaotic sea) into K = 4 clusters is performed using the information from eigenvectors q , q and q only.The membership probabilities for the central forecast are presented in Figure 10(a) for trajectories with p i,k ≥ .
5. The domain is partitioned into 4 quadrants of similar size, and gaps between clusters correspondto particles with membership probabilities lower than 0 . p i,k and sample standard deviations S i,k for k ∈ { , . . . , } .The superimposed mean p i, ˆ k and standard deviation S i, ˆ k are presented in Figure 10(b,c). The parameter-ization used for the model produces only a modest variation in the trajectory outcomes over six hours asdemonstrated in Figure 9(c). Regions of highest uncertainty for this system correspond to identifying theedges of the clusters accurately, but this level of uncertainty is significantly lower than those observed in theBickley jet example. The most pronounced uncertainty regions, in Figure 10(c), correspond to the boundarybetween clusters 1 and 4, followed by that between clusters 1 and 2. The experiment targeting predicted coherent structures consisted of surface drifter releases that took placeon August 7, 2018, in the vicinity of Nomans Land [Filippi et al., 2020]. The CODE drifters used inthe experiment have technical specifications listed in [Serra et al., 2020; Filippi et al., 2020]. Drifters ofthe same design are routinely used by the U.S. Coast Guard in search-and-rescue operations, as well as inprevious field experiments in the coastal region near Martha’s Vineyard [Rypina et al., 2014, 2016]. Thedrifters were equipped with GPS transmitters that provided positioning fixes every 5 min, with an accuracyof a few meters. An elliptical route around Nomans Land was used for the drifter deployment, employingtwo WHOI vessels to minimize ship time, so that all drifters were in water by the start of the interval ofanalysis. Eighteen drifters were deployed in the water around the predicted locations of the targeted coherentstructures. The position of the drifters at the starting time of our analysis, t = 16:00 UTC, is presented inFigure 11(a).Figure 11(a-d) presents the mean cluster trajectories for each of the 4 clusters, on which we superposethe drifter positions from [Filippi et al., 2020]. To plot the evolution of p i,k , the entire domain is first splitinto 175 ×
175 bins. Then, at each time instance, an average membership probability for each bin, over the15 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020
Figure 10: ( a ) Central forecast membership probabilities at t = t , for clusters k ∈ { , . . . , } . The black boxencloses the domain where trajectories are initialized. Clustering statistics for ensemble of R = 71 forecasts:superimposed membership probability ( b ) mean p i, ˆ k and ( c ) standard deviation S i, ˆ k . The color mapspresented in ( a ) are also used in Figure 11.forecasts, is calculated. This average is based on particles inside each bin at time t , and weighted by theirmembership probabilities p i,k . The overall system dynamics are similar to what was observed for the centralforecast in Figure 9(a,b), with the clusters highlighting groups of different behavior. Using the labels fromFigure 10(a), while cluster 1 (purple) travels a longer distance to the northeast, ending north of Martha’sVineyard, cluster 2 (blue) moves a shorter distance northward, while clusters 3 and 4 (red and yellow) areless dynamic and remain mostly to the east of Nomans Land.Drifters are marker-coded according to the mean membership probabilities p i,k of their spatial positionat t . Triangles correspond to cluster 1, squares to cluster 2, circles to cluster 3, and no drifters wereinitialized in cluster 4. Drifters initialized in higher uncertainty regions (see Figure 10(c)) are plotted ascrosses, and correspond to regions where S i, ˆ k > .
1, or p i, ˆ k < .
7. The results demonstrate that the drifterspredominately remain inside of the forecast model clusters during the first four hours of the time interval,and exceptions to this behavior are associated with higher uncertainty. Consider, for example, the threedrifters represented as crosses, initially west of Nomans Land. Two of them move along the northern coastof the island, with one of them beaching and the other ending on cluster 2. The third drifter headed southfirst, then east, also toward cluster 2. This fact highlights the value of the uncertainty quantification analysisto understand different dynamical behaviors in the flow which are not captured by the central forecast alone.All drifters were eventually advected eastward after 20:00.Some discrepancies between the clustering predicted behavior and the drifter trajectories were observed.Wind gusts that occurred between 16:00 and 20:00 significantly affected the drifter trajectories. These gustshad not been predicted by the forecasts from the National Centers for Environmental Prediction used bythe MIT-MSEAS model for the atmospheric forcing [Filippi et al., 2020]. The study was therefore limitedby differences between observed flows and model flows, and the coherent structure analysis accuracy islimited to the accuracy of the velocity field used for processing. Nonetheless, drifters released in the clusterswith high membership probabilities tend to behave similarly to the clusters, even if they do not preciselymatch the predicted trajectory behavior. This fact highlights the coherent structure robustness, even inuncertain conditions. The clustering analysis partitions the domain into robust regions where a drifterreleased in the region will remain there or nearby over the period of analysis. The uncertainty quantificationanalysis helped identifying the key structures delimiting regions with different transport behaviors, furthershowing and expanding the applicability of trajectory clustering for studying oceanic flows, despite modelimperfections. 16 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts
Figure 11: Time evolution of the average clusters in a binned domain. The 18 drifters are marker-codeddepending on the cluster in which they started. Crosses represent the 4 drifters initialized at locations ofhigher uncertainty, with S i, ˆ k > .
1. The times correspond to ( a ) t = 16:00, ( b ) t = 18:00, ( c ) t = 20:00, and( d ) t f = 22:00 UTC. Triangles are initialized on the purple cluster (3 drifters), squares on the blue cluster(5 drifters) and circles on the orange cluster (6 drifters). Corresponding colorbars presented in Figure 10(a).(See Supplementary Material for a video of the time evolution.) Ensemble statistics of the trajectory clustering results provide a partitioning of the fluid domain that mayprovide critical information in emergency response situations, such as search-and-rescue operations, whenoperational decisions about optimal resource allocation need to be made quickly, accurately, and account formodel uncertainties. We presented a modified version of the spectral clustering method with soft membershipprobabilities, and applied it to fluid particle trajectories to identify coherent structures first in an analyticflow model, then in forecast simulations provided by a coastal ocean model. Uncertainty quantificationwas applied to assess both the result sensitivity to the clustering method free-parameters and the clustervariability with unknown parameters of the model data. The method sensitivity study, performed on theanalytic quasi-periodic Bickley jet system, identified the similarity radius as the free-parameter to which theclusters are most sensitive. To mimic model uncertainty, the Bickley jet parameters were varied to performa model sensitivity study that highlights the robustness of vortex cores compared to the more uncertainvortex perimeters. 17 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020
Finally, the method was applied to an ocean ensemble forecast of the coastal region of Martha’s Vineyard,and the clustering results were compared to drifter trajectories from a drifter release experiment targetingpredicted coherent structures. The forecast clusters from the ensemble analysis provided a good baseline forthe drifter behavior, as drifters deployed within each cluster performed similar motions to their correspondingclusters. Ocean transport predictions are challenging due to the complexity of the underlying dynamicsgoverning the flow, and while the Lagrangian approach ultimately depends on the accuracy of the availablevelocity fields and the quality of the model data, the results presented here demonstrate that the identifiedclusters are robust to uncertainties in the model and able to predict the main elements of the organizationof flow transport. The coupling of the clustering approach to uncertainty quantification can provide a morecomplete and informative description of flow transport and areas of higher and lower uncertainty withindifferent clusters. Despite not having been the case for the drifter release presented here, a clusteringforecast analysis could be used in planning a drifter deployment for future experiments.Further refinement of the trajectory clustering method is highly desirable, in particular aiming to reduceparameter sensitivity. The uncertainty quantification with respect to method free-parameters could beused to select parameters that minimize cluster variability, in order to identify clusters that are physicalstructures, and not a byproduct of the system and parameters chosen. On a more fundamental level,while the method presented here provides robust clusters, it may be possible to improve the method byincorporating fundamental changes to the similarity measure, rather than addressing the sensitivity to free-parameters only. To quantify trajectory similarity, one could not only consider the similarity between theparticle spatial coordinates in time, but also that of their velocity vectors, and propose a hybrid notion ofsimilarity. Finally, applying the method to other oceanic forecasts and using the forecast clustering resultsto plan and execute drifter release experiments can be a promising path to more effective experiments,by increasing the likelihood of released drifters capturing targeted coherent structures and ocean transportbarriers.
Author Contributions
Conceptualization, G.S.V. and M.R.A.; Formal analysis, G.S.V.; Funding acquisition, I.I.R.; Methodology,G.S.V. and M.R.A.; Software, G.S.V.; Supervision, M.R.A.; Visualization, G.S.V.; Writing – original draft,G.S.V.; Writing – review & editing, G.S.V., I.I.R. and M.R.A.
Funding
This research was funded by National Science Foundation Division of Atmospheric and Geospace Sciences,award number 1520825.
Acknowledgments
We thank Margaux Filippi and Thomas Peacock for sharing the drifter data set, Pierre Lermusiaux forproviding the forecast ensemble, and H. M. Aravind for his helpful suggestions.
Conflicts of Interest
The authors declare no conflict of interest.
References
Allshouse, M. R. and Peacock, T. (2015). Lagrangian based methods for coherent structure detection.
Chaos:An Interdisciplinary Journal of Nonlinear Science , 25(9):097617.Allshouse, M. R. and Thiffeault, J.-L. (2012). Detecting coherent structures using braids.
Physica D:Nonlinear Phenomena , 241(2):95–105.Ameli, S. and Shadden, S. C. (2019). A transport method for restoring incomplete ocean current measure-ments.
Journal of Geophysical Research: Oceans , 124(1):227–242.18 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts
Ameli, S. and Shadden, S. C. (2020). Trajectory reconstruction and analysis for coherent structure evaluation(trace). http://transport.me.berkeley.edu/trace/.Beron-Vera, F. J., Olascoaga, M. J., Brown, M. G., Ko¸cak, H., and Rypina, I. I. (2010). Invariant-tori-likelagrangian coherent structures in geophysical flows.
Chaos: An Interdisciplinary Journal of NonlinearScience , 20(1):017514.Beron-Vera, F. J., Olascoaga, M. J., Haller, G., Farazmand, M., Tri˜nanes, J., and Wang, Y. (2015). Dissipa-tive inertial transport patterns near coherent lagrangian eddies in the ocean.
Chaos: An InterdisciplinaryJournal of Nonlinear Science , 25(8):087412.Beron-Vera, F. J., Wang, Y., Olascoaga, M. J., Goni, G. J., and Haller, G. (2013). Objective detection ofoceanic eddies and the agulhas leakage.
Journal of Physical Oceanography , 43(7):1426–1438.Boffetta, G., Lacorata, G., Redaelli, G., and Vulpiani, A. (2001). Detecting barriers to transport: a reviewof different techniques.
Physica D: Nonlinear Phenomena , 159(1-2):58–70.Everitt, B., Landau, S., Leese, M., and Stahl, D. (2011).
Cluster Analysis . Wiley Series in Probability andStatistics. Wiley, 5 edition.Filippi, M., Rypina, I. I., Hadjighasem, A., and Peacock, T. (2020). A parameter-free spectral clusteringapproach to coherent structure detection in geophysical flows.
Submitted to Fluids .Filippone, M., Camastra, F., Masulli, F., and Rovetta, S. (2008). A survey of kernel and spectral methodsfor clustering.
Pattern recognition , 41(1):176–190.Froyland, G. (2015). Dynamic isoperimetry and the geometry of lagrangian coherent structures.
Nonlinearity ,28(10):3587.Froyland, G. and Padberg-Gehle, K. (2015). A rough-and-ready cluster-based approach for extracting finite-time coherent sets from sparse and incomplete trajectory data.
Chaos: An Interdisciplinary Journal ofNonlinear Science , 25(8):087406.Hadjighasem, A., Farazmand, M., Blazevski, D., Froyland, G., and Haller, G. (2017). A critical comparisonof lagrangian methods for coherent structure detection.
Chaos: An Interdisciplinary Journal of NonlinearScience , 27(5):053104.Hadjighasem, A., Karrasch, D., Teramoto, H., and Haller, G. (2016). Spectral-clustering approach to la-grangian vortex detection.
Physical Review E , 93(6):063107.Haley, P. J. and Lermusiaux, P. F. (2010). Multiscale two-way embedding schemes for free-surface primitiveequations in the “multidisciplinary simulation, estimation and assimilation system”.
Ocean dynamics ,60(6):1497–1537.Haley Jr, P. J., Agarwal, A., and Lermusiaux, P. F. (2015). Optimizing velocities and transports for complexcoastal regions and archipelagos.
Ocean Modelling , 89:1–28.Haller, G. (2002). Lagrangian coherent structures from approximate velocity data.
Phys. Fluids A , 14:1851–1861.Haller, G. (2015). Lagrangian coherent structures.
Annual Review of Fluid Mechanics , 47:137–162.Haller, G. and Beron-Vera, F. J. (2012). Geodesic theory of transport barriers in two-dimensional flows.
Physica D , 241:1680–1702.Haller, G. and Beron-Vera, F. J. (2013). Coherent lagrangian vortices: The black holes of turbulence.
Journalof Fluid Mechanics , 731.Haza, A., Griffa, A., Martin, P., Molcard, A., ¨Ozg¨okmen, T., Poje, A., Barbanti, R., Book, J., Poulain, P.,Rixen, M., et al. (2007). Model-based directed drifter launches in the adriatic sea: Results from the dartexperiment.
Geophysical Research Letters , 34(10). 19 . S. Vieira, I. I. Rypina, and M. R. Allshouse, 2020
Haza, A. C., ¨Ozg¨okmen, T. M., Griffa, A., Molcard, A., Poulain, P.-M., and Peggion, G. (2010). Transportproperties in small-scale coastal flows: relative dispersion from vhf radar measurements in the gulf of laspezia.
Ocean Dynamics , 60(4):861–882.Jacobs, G. A., Bartels, B. P., Bogucki, D. J., Beron-Vera, F. J., Chen, S. S., Coelho, E. F., Curcic, M.,Griffa, A., Gough, M., Haus, B. K., et al. (2014). Data assimilation considerations for improved oceanpredictability during the gulf of mexico grand lagrangian deployment (glad).
Ocean Modelling , 83:98–117.Jain, A. K., Murty, M. N., and Flynn, P. J. (1999). Data clustering: a review.
ACM computing surveys(CSUR) , 31(3):264–323.Kratzke, T. M., Stone, L. D., and Frost, J. R. (2010). Search and rescue optimal planning system. In , pages 1–8. IEEE.Lermusiaux, P. (2002). On the mapping of multivariate geophysical fields: Sensitivities to size, scales, anddynamics.
Journal of Atmospheric and Oceanic Technology , 19(10):1602–1637.Lermusiaux, P. F. (2006). Uncertainty estimation and prediction for interdisciplinary ocean dynamics.
Journal of Computational Physics , 217(1):176–199.Lermusiaux, P. F., Chiu, C.-S., Gawarkiewicz, G. G., Abbot, P., Robinson, A. R., Miller, R. N., Haley, P. J.,Leslie, W. G., Majumdar, S. J., Pang, A., et al. (2006). Quantifying uncertainties in ocean predictions.
Oceanography , 19:92–105.McWilliams, J. C. (1984). The emergence of isolated coherent vortices in turbulent flow.
Journal of FluidMechanics , 146:21–43.Ng, A. Y., Jordan, M. I., and Weiss, Y. (2002). On spectral clustering: Analysis and an algorithm. In
Advances in neural information processing systems , pages 849–856.Nock, R., Vaillant, P., Henry, C., and Nielsen, F. (2009). Soft memberships for spectral clustering, withapplication to permeable language distinction.
Pattern Recognition , 42(1):43–53.Olascoaga, M. J., Beron-Vera, F. J., Haller, G., Trinanes, J., Iskandarani, M., Coelho, E., Haus, B. K.,Huntley, H., Jacobs, G., Kirwan, A., et al. (2013). Drifter motion in the gulf of mexico constrained byaltimetric lagrangian coherent structures.
Geophysical Research Letters , 40(23):6171–6175.Peacock, T. and Dabiri, J. (2010). Introduction to focus issue: Lagrangian coherent structures.Provenzale, A. (1999). Transport by coherent barotropic vortices.
Annual review of fluid mechanics , 31(1):55–93.Rypina, I., Brown, M. G., Beron-Vera, F. J., Ko¸cak, H., Olascoaga, M. J., and Udovydchenkov, I. (2007). Onthe lagrangian dynamics of atmospheric zonal jets and the permeability of the stratospheric polar vortex.
Journal of the Atmospheric Sciences , 64(10):3595–3610.Rypina, I., Kirincich, A., Limeburner, R., and Udovydchenkov, I. (2014). Eulerian and lagrangian correspon-dence of high-frequency radar and surface drifter data: Effects of radar resolution and flow components.
Journal of Atmospheric and Oceanic Technology , 31(4):945–966.Rypina, I. I., Kirincich, A., Lentz, S., and Sundermeyer, M. (2016). Investigating the eddy diffusivity conceptin the coastal ocean.
Journal of Physical Oceanography , 46(7):2201–2218.Rypina, I. I. and Pratt, L. J. (2017). Trajectory encounter volume as a diagnostic of mixing potential influid flows.
Nonlinear Processes in Geophysics , 24(2).Schlueter-Kuck, K. L. and Dabiri, J. O. (2017). Coherent structure colouring: identification of coherentstructures from sparse data using graph theory.
Journal of Fluid Mechanics , 811:468–486.Serra, M., Sathe, P., Rypina, I., Kirincich, A., Ross, S. D., Lermusiaux, P., Allen, A., Peacock, T., andHaller, G. (2020). Search and rescue at sea aided by hidden flow structures.
Nature Communications ,11:2525. 20 ncertainty Quantification of Trajectory Clustering Applied to Ocean Ensemble Forecasts
Shi, J. and Malik, J. (2000). Normalized cuts and image segmentation.
IEEE Transactions on patternanalysis and machine intelligence , 22(8):888–905.Vieira, G. S. and Allshouse, M. R. (2020). Internal wave boluses as coherent structures in a continuouslystratified fluid.
Journal of Fluid Mechanics , 885:A35.von Luxburg, U. (2007). A tutorial on spectral clustering.
Statistics and computing , 17(4):395–416.Williams, M. O., Rypina, I. I., and Rowley, C. W. (2015). Identifying finite-time coherent sets from limitedquantities of lagrangian data.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 25(8):087408.21
Similarity Matrix Sparsification
As mentioned in Section 2.1, the use of the Gaussian similarity measure w ij = exp (cid:0) − r ij / σ (cid:1) (9)in the similarity matrix W greatly reduces the impact of the sparsification step compared to the l x /r ij measure used in Hadjighasem et al. [2016]. This result is due to w ij in (9) approaching zero faster for large r ij .It is mentioned in Section 3.1, however, that sparsifying the matrix is worthwhile to reduce computationalcosts and storage, and we therefore sparsify entries of W that satisfy w ij < exp( − / r ij > σ . This sparsification rule is used in Sections 3.2, 3.3, 3.4 and 4.1, under the assumption that suchsparsification has negligible impact on the clustering membership probability results.Here, we present a study for the Bickley jet system with model parameters as in Section 3 and methodfree-parameters σ/l x = 0 . m = 2, and K = 7, where the clustering is performed using different levels ofsparsification. We vary the sparsification levels by setting w ij = 0 for r ij /σ > α , with α ∈ { . , , , , } and α = ∞ representing the non-sparsified result matching Figure 1(a).Figure 12(a) presents a histogram with the distribution of r ij values. The vast majority of these valuesrepresent the time-averaged distance between particles that describe dissimilar trajectories, and thereforecorrespond to graph edges of negligible weights. Figure 12(b) presents the weight dependence on r ij /σ , andtherefore the maximum magnitude of matrix entries that are dropped when W is sparsified for a choice of α .Figure 12: ( a ) Histogram of the time-averaged distance values, made dimensionless by σ = 0 . /l x , witha bin width of 0.25. ( b ) Weight dependence on r ij /σ for the Gaussian measure (9). The w ij values for r ij /σ = 0 .
5, 1, 2, 4 and 8 are shown in blue.The impact of the sparsification parameter α on the number of retained (nonzero) entries in the similaritymatrix is shown in Table 1. For all finite choices of α presented here, the percent sparsification of W isgreater than 95%. Table 1: Dependence of W entries on the sparsification parameter α . Sparsification parameter α Number of nonzero entries in W Percent Sparsification (%) ∞ α are presented in Figure 13.We notice that sparsifying with α = 0 .
5, 1 or 2 clearly has an effect on the resulting membership probabilities,s clusters shrink as α is reduced. For α ≥
4, however, the probabilities no longer depend on α , and nodifference is observed between the method output for α = 4 or the non-sparsified case ( α = ∞ ).Figure 13: Membership probabilities for vortices k ∈ { , . . . , } , for different choices of the sparsificationparameter α , plotted at t = t . Values correspond to ( a ) α = 0 .
5, ( b ) 1, ( c ) 2, ( d ) 4, ( e ) 8, and ( f ) + ∞ .Corresponding colorbars presented in Figure 1Sparsifying with α = 4 also reduces the number of nonzero entries in W by about 60 times, and thereforefacilitates storage and reduces the computational time when solving the generalized eigenproblem. Note thatthis result is valid for the Bickley jet parameterized as in Section 3, but the impact of sparsification mayvary for an arbitrary flow. Video Description
The four supplementary videos, for which captions are available below, correspond to the time evolution ofthe results presented in Figures 2, 6, 8 and 11. (cid:63) video02.avi
Membership probabilities for vortices k ∈ { , . . . , } , for different choices of the parameter σ . Valuescorrespond to ( a ) σ/l x = 0 . b ) 0 . c ) 0 .
030 and ( d ) 0 . (cid:63) video04.avi Examples of membership probabilities for the six identified clusters, for different model parameters { A n } and { φ n } drawn from normal distributions. Cases correspond to parameters ( A , A , A , φ /l x , φ /l x , φ /l x ) equal to ( a ) (0.0087, 0.19, 0.20, 0.01, 0.06, -0.03), ( b ) (0.0069, 0.26, 0.25, 0.00, -0.08,0.09), ( c ) (0.0102, 0.35, 0.25, 0.01, -0.01, 0.01), where the jet is identified and one of the vorticesmissed, and ( d ) (0.0077, 0.11, 0.32, 0.00, -0.01, 0.02), where the jet (grayscale) is identified and twovortices are merged. Corresponding colorbars presented in Figures 1 and 4. (cid:63) video08.avi Evolution of trajectories from the 1000 realizations released from a low (orange, top) and high (red,bottom) uncertainty position. Orange particles are initialized at the core of vortex 1, while red particlesare initialized at the perimeters of vortex 4. Transparency is used to highlight high or low concentra-tion of particles. The vortex boundaries from the central realization are plotted in gray and encloseparticles with membership probabilities greater than 0.5, for clusters 1 and 4. (cid:63) video11.avi
Time evolution of the average clusters in a binned domain. The 18 drifters are marker-coded dependingon the cluster in which they started. Crosses represent the 4 drifters initialized at locations of higheruncertainty, with S i, ˆ k > ..