extfieldsofstudy:
Computer Science
Mathematics

text:
Divergence-conforming methods for transient doubly-diffusive flows: A priori and a posteriori error analysis

The analysis of the double-diffusion model and $\mathbf{H}(\mathrm{div})$-conforming method introduced in [B\"urger, M\'endez, Ruiz-Baier, SINUM (2019), 57:1318--1343] is extended to the time-dependent case. In addition, the efficiency and reliability analysis of residual-based {\it a posteriori} error estimators for the steady, semi-discrete, and fully discrete problems is established. The resulting methods are applied to simulate the sedimentation of small particles in salinity-driven flows. The method consists of Brezzi-Douglas-Marini approximations for velocity and compatible piecewise discontinuous pressures, whereas Lagrangian elements are used for concentration and salinity distribution. Numerical tests confirm the properties of the proposed family of schemes and of the adaptive strategy guided by the {\it a posteriori} error indicators.

Introduction and problem formulation
1.1. Scope. A number of physical problems of relevance in industrial applications involve coupled incompressible flow and double-diffusion transport. We are interested in numerical schemes for the approximation of a class of coupled equations that arise as models of sedimentation of small particles under the effect of salinity of the ambient fluid. The governing model can be stated as follows (cf., e.g., [19,33]): posed on a spatial domain Ω ⊂ R d , d = 2 or d = 3, where t ∈ (0, t end ] is time, u is the fluid velocity, ν is the concentration-dependent viscosity, ρ m is the mean density of the fluid, p is the fluid pressure, ρ is density, g is the gravity acceleration, s is the salinity concentration, c is the concentration of solid particles, Sc = ν ref /κ s is the Schmidt number, where κ s is the diffusivity of salinity, ν ref is a reference viscosity in the absence of solid particles, and τ = κ s /κ c is the inverse of the diffusivity ratio, where κ c is the diffusivity of solid particles, and e z is the upward-pointing unit vector. We relate the densities through a linearised equation of state ρ = ρ m (αs + βc).
Again as in [19,33], the solid particles are assumed mono-sized with radius r, and settle at dimensionless veloc- is the Stokes velocity (settling velocity of a single particle in an unbounded fluid). The coupling mechanisms between flow and transport are only due to advection for concentration and salinity (where the advecting velocity for concentration, u − v p e z , is also divergence-free), and through the concentration-dependent viscosity. Further details are provided in later parts of the paper.
To put the paper further into the proper perspective, we mention that there exists an abundant body of literature devoted to constructing accurate finite element and related schemes for doubly-diffusive flows. Some recent contributions include variational multiscale stabilised schemes, least-squares methods, divergence-conforming mixed methods, volume-averaging discretisations, spectral elements, vorticity-based finite element formulations, and similar methods applied to, e.g., flows with heat and mass transport [1], reactive Boussinesq flows [2], nonlinear advection-reaction-diffusion in the context of bioconvective flows [10,31], cross-diffusion and boundary layer effects in doubly-diffusive Navier-Stokes-Brinkman equations [18,20] and in Darcy-Brinkman equations [39], or phase change models [21,38,40,41]; where the list is far from exhaustive.
The solvability analysis for the continuous and discrete problems usually follows energy and fixed-point schemes as done for classical Boussinesq equations, and this is also the approach we follow here. The discretisation in space uses an interior penalty divergence-conforming method for the flow equations (in this case, Brezzi-Douglas-Marini (BDM) elements of degree k ≥ 1 for the velocity and discontinuous elements of degree k − 1 for the pressure following [11,29]), combined with Lagrangian elements for the diffusive quantities, and the development stands as a natural extension of the formulation in [18] to the transient case. As such, it also features exactly divergencefree velocity approximations ensuring local conservativity and energy stability, and the error estimates of velocity are pressure-robust. The chosen time discretisation is the backward differentiation formula of degree 2 (BDF2), which for k = 2 gives a method of order 2 in space and time. Existence of discrete solutions is established by the Brouwer fixed-point theory similarly as in [18], and the error analysis in the semi-discrete and fully-discrete settings is adapted from the theory of [5] for the Boussinesq equations.
In many applications where double-diffusion effects occur, complicated flow patterns exist in zones far from boundary layers and sufficiently refined meshes are needed essentially in the whole spatial domain. However, for salinity-driven settling of solid particles that result in mathematical models such as (1.1), many of the flow features are clustered near zones of high-gradients of concentration, which is where the typical plumes are observed [19,31]. This motivates the use of adaptive mesh refinement guided by a posteriori error indicators. For instance, in the context of phase change models there are some results based on error-related metric change [21,34] and on goal-oriented adaptivity [41]. Regarding the design and rigorous analysis of residual-based a posteriori error estimators for flow-transport couplings, the literature is predominantly focused on the stationary case (see, e.g., [3,6,7,9,13,22,37,42] and the references therein). Only a few results are available for the time-dependent regime, from which we mention the adaptive mixed method for Richards equation in porous media [15], the remeshing scheme based on goal-oriented adaptivity for solidification problems advanced in [14], the collection of adaptive schemes for reactive flow discussed in [16] and for heat transfer in [30]. However, none of these theoretical frameworks is directly applicable to (1.1) using divergence-conforming approximations.
The a posteriori error analysis we advance here is of residual type, and its analysis uses ideas from the abstract results related to spatial estimators for discontinuous Galerkin schemes applied to parabolic problems in [24]. The approach hinges on a decomposition of the discrete solution into a conforming and a non-conforming contribution, along with a reconstruction technique (see also [32]). This has also been exploited for the construction of a posteriori estimators of time-dependent Stokes and Navier-Stokes equations [12,43]. Our a posteriori error analysis is divided into three parts. In first part, we present the error estimator for the steady coupled problem. In second part, we extend the a posteriori error estimation to the semi-discrete method, and finally we present the a posteriori error estimator for the unsteady coupled problem. We restrict that part of the analysis to the simpler backward Euler method. To the best of our knowledge, the a posteriori error estimation advanced in this paper is the first comprehensive study targeted for transient doubly-diffusive flows.
The remainder of this paper is organised as follows. In what is left of this section we outline the weak formulation of (1.1) and state the stability the continuous problem. In Section 2 we introduce the divergenceconforming method in fully discrete form, show existence of discrete solutions using fixed-point arguments, and rigorously establish a priori error estimates. Section 3 is devoted to the construction and analysis of efficiency and reliability for a residual-based a posteriori error estimator tailored for the stationary problem. In turn, these upper and lower bounds are used to establish properties of a second family of estimators for the transient case, and addressed in Sections 4 and 5. In Section 6 we collect numerical tests that verify the theoretical convergence rates predicted by the a priori error analysis, confirm the robustness of the proposed a posteriori error estimators, and illustrate the advantages of adaptive methods in the simulation of doubly-diffusive flows.
A problem similar to (1.2) has been studied in [4]. Assuming that F belongs to L 2 (0, t end ; H −1 (Ω)), that the initial velocity u 0 belongs to L 2 (Ω) and the initial data for the coupled species (s, c in the context of our problem) belongs to L 2 (Ω), the authors showed existence of a solution by using the Galerkin method and applying the Cauchy-Lipschitz theorem, and proved its uniqueness in two dimensions. Such analysis can be applied to (1.2) by noting that F is a Lipschitz-continuous function, and assuming the initial data belongs to appropriate spaces. This is, however, not the focus of the paper.

Finite element discretisation and a priori error bounds
We discretise space by a family of regular partitions, denoted T h , of Ω ⊂ R d into simplices K (triangles in 2D or tetrahedra in 3D) of diameter h K . We label by K − and K + the elements adjacent to an edge, while h e stands for the maximum diameter of the edge. If v and w are smooth vector and scalar fields defined on T h , then (v ± , w ± ) denote the traces of (v, w) on e that are the extensions from the interior of K + and K − , respectively. Let n e denote the outward unit normal vector to e on K. The average {{·}} and jump · operators are defined as whereas for boundary jumps and averages we adopt the conventions {{v}} = v = v, and {{w}} = w = w. In addition, ∇ h will denote the broken gradient operator.
For k ≥ 1 and a mesh T h on Ω, let us consider the discrete spaces (see e.g. [17]) [29]). Here P k (K) denotes the local space spanned by polynomials of degree up to k and V h is the space of divergence-conforming BDM elements. We then state the following semi-discrete Galerkin formulation for problem (1.2) The discrete versions of the trilinear forms a h 1 (·; ·, ·) and c h 1 (·; ·, ·) are defined using a symmetric interior penalty and an upwind approach, respectively (see, e.g., [11,29]): and u e h is the trace of u h taken from within the exterior of K.
We partition the interval [0, t end ] into N subintervals [t n−1 , t n ] of length ∆t. We use the implicit BDF2 scheme where all first-order time derivatives are approximated using the centred operator (similarly for ∂ t c) and for the first time step a first-order backward Euler method is used from t 0 to t 1 , starting from the interpolates u 0 h , s 0 h of the initial data. In what follows, we define the difference operator Dy n+1 := 3y n+1 − 4y n + y n−1 for any quantity indexed by the time step n. For instance, (2.3) can be written as The resulting set of nonlinear equations is solved by an iterative Newton-Raphson method with exact Jacobian. Hence the complete discrete system is given by For the subsequent analysis, we introduce for r ≥ 0 the broken H r (T h ) space as well as the following mesh-dependent broken norms where the stronger norm · 2,T h is used to show continuity. This norm is equivalent to · 1,T h on H 1 (T h ) (see [11]). Finally, adapting the argument used in [28,Proposition 4.5] we have the discrete Sobolev embedding: for r = 2, 4 there exists a constant C emb > 0 such that Using these norms, we can establish continuity of the trilinear and bilinear forms involved. The proof follows from [11,Section 4].
Lemma 2.1. The following properties hold: and for all u, v, w ∈ H 1 (T h ) and ϕ, ψ ∈ H 1 (Ω), there holds where the constantC Lip > 0 is independent of h (cf. [18]). Note that while the coercivity of the form a 2 (·, ·) in the discrete setting is readily implied by (1.5b), there also holds (cf. [29,Lemma 3 provided that the stabilisation parameter a 0 > 0 in (2.2) is sufficiently large and independent of the mesh size.
Let w ∈ H 0 (div 0 ; Ω), then due to the skew-symmetric form of the operators c h 1 and c 2 , and the positivity of the non-linear upwind term of c h 1 , we can write as well as the following relation (which is based on (2.5) and follows by the same method as in [28]): For any We also have Finally, we recall from [29] the following discrete inf-sup condition for b(·, ·), whereβ is independent of h: We will use the following algebraic relation: for any real numbers a n+1 , a n , a n−1 and defining Λa n := a n+1 − 2a n + a n−1 , we have 2(3a n+1 − 4a n + a n−1 , a n ) = |a n+1 | 2 + |2a n+1 − a n | 2 + |Λa n | 2 − |a n | 2 − |2a n − a n−1 | 2 .
(2.12)  We proceed in the same way taking ψ h = 4c n+1 h in the fourth equation of (2.4), to get the third result. We get the first result by substituting bounds for c h and s h into (2.14).
Theorem 2.2 (Existence of discrete solutions). The problem (2.4) admits at least one solution Proof of Theorem 2.2. To simplify the proof we introduce the following constants: ,Ω . We proceed by induction on n ≥ 2. We define the mapping

Note this map is well-defined and continuous on
and employ (2.9a), (2.9b) and (2.8), we obtain Next, using (2.13) we deduce that Then, setting Hence, the right-hand side is nonnegative on a sphere of radius r := C r /C R . Consequently, by Theorem 2. Let us denote by I h : C(Ω) → M h the nodal interpolator with respect to a unisolvent set of Lagrangian interpolation nodes associated with M h . Π h u will denote the BDM projection of u, and L h p is the L 2 -projection of p onto Q h . Under usual assumptions, the following approximation properties hold (see [29]): The following development follows the structure adopted in [5].
Lemma 2.2. Assume that u ∈ H 2 (Ω), p ∈ L 2 (Ω) and s, c ∈ H 1 (Ω). Then we have Proof. Since we assume u ∈ H 2 (Ω), integration by parts yields the required result. See also [11]. The third and fourth equations are a straightforward result from the continuous form.
Now we decompose the errors as follows: Assuming that u 0 h = Π h u(0), s 0 h = I h s(0) and c 0 h = I h c(0), we will use also the notation E n u = (u(t n )−Π h u(t n )) and ξ n u = (Π h u(t n ) − u n h ), and similar notation for other variables. Since for the first time iteration of system (2.1) we adopt a backward Euler scheme, an error estimate is required for this step.
and also that u L ∞ (0,t end ;W 1,∞ (Ω)) < M for a sufficiently small M > 0 (a precise condition can be found in Theorem 2.5). Then there exist positive constants C 1 u , C 1 s , C 1 c , independently of h and ∆t, such that Proof. As in the continuous case, we define the discrete kernel of the bilinear form b(·, ·) as and relying on the inf-sup condition (2.11), we can continue with an equivalent discrete problem without pressure.
Taking into account the assumed regularity for u, we have for all x, a γ(x) ∈ (0, 1) such that then u satisfies the following error equation , which results after choosing ξ 1 u as test function in the first equation of the reduced form of Lemma 2.2 and system (2.1), performing an Euler scheme step, subtracting both equations, and adding ±a h 1 (c 1 h ; Π h u(∆t), ξ 1 u ). Now, invoking the approximation estimates (2.15), Young's inequality, and the stability properties, we get (2.17) Next, we choose ξ 1 s as test function in (2.16c) and system (2.1); we follow the same steps as before, adding to the sum of both equations the term ±a 2 (I h s 1 , ξ 1 c ) to obtain 1 4 ξ 1 In the same way, choosing ξ 1 c as test function in (2.16d) and in (2.1) we obtain Now, from (2.17) we deduce that and u L ∞ (0,t end ;H 1 (Ω)) < M for a sufficiently small constant M > 0. Then there exist positive constants C, γ 1 ≥ 0 independent of h and ∆t such that for all m + 1 ≤ N , Proof. We appeal to the reduced form of the problem again, taking u h ∈ X h and u ∈ X, then we choose as test function v h = ξ n+1 u in the first equation of (2.4) and insert the terms We consider (2.16a) (see Lemma 2.2) at t = t n+1 and v = ξ n+1 u . Inserting the term ±(Du(t n+1 ), ξ n+1 u ) Ω /(2∆t), we readily deduce the identity We can then subtract (2.23) from (2.22) and multiply both sides by 4∆t, yielding I 1 + I 2 + · · · + I 7 = 0, with Du(t n+1 ), ξ n+1 u Ω , .
For the first term, using (2.12) we can assert that ,Ω . Using the ellipticity stated in (2.8), we readily get By using Taylor's formula with integral remainder we have then by combining Cauchy-Schwarz and Young's inequality, we obtain the bound Now we insert ±4∆tE u (t n+1 ) into the fourth term, which leads to Proceeding as before and using (2.15) on the first term of I 4 , we get . Then by (2.15), (2.7), Lemma 2.1 and Young's inequality we immediately have Adding and subtracting suitable terms within I 6 yields where we defineĨ The bound (2.10) and (2.15) imply that where C * is a positive constant coming from (2.15). We also have Hence, by choosing ε i = 3α a /5 for i = 1, . . . , 10, collecting the above estimates, and summing over 1 ≤ n ≤ m for all m + 1 ≤ N , we get where Finally, Theorem 2.4 yields the desired result.
Theorem 2.6. Let (u, p, s, c) be the solution of (1.2) under the assumptions of Section 1.4, and (u h , p h , s h , c h ) be the solution of (2.4). If Proof. Proceeding similarly as for Theorem 2.5, we choose as test function ϕ h = ξ n+1 s in the second equation of (2.4) and insert suitable additional terms to obtain the following identity (analogous to (2.22 From (2.16b), focusing on t = t n+1 , using ϕ = ξ n+1 s and proceeding as in the derivation of (2.23), we obtain Next we subtract (2.25) from (2.24) and multiply both sides by 4∆t. This yieldsÎ 1 + · · · +Î 6 = 0, wherê Inserting ±4∆tE s (t n+1 ) intoÎ 4 and using (2.15) leads to the bound Employing again (2.15) in combination with (1.3a) we have In order to derive a bound forÎ 6 we proceed as for the bound on I 7 in the proof of Theorem 2.5; namely adding and subtracting suitable terms in the definition ofÎ 6 , definingĨ 6 in this case bŷ and applying (2.9b), (2.6a), (2.15) and Young's Inequality to the result, we get In this manner, and after choosing ε i = 6α a /7 for i = 1, . . . , 7, we can collect the above estimates and sum over 1 ≤ n ≤ m, for all m + 1 ≤ N , to get This concludes the proof.  Proof. It follows along the same lines of the proof of Theorem 2.6, with constant γ 3 given by Proof. Owing to the inf-sup condition (2.11), there exists From (2.4) and Lemma 2.2, proceeding as in the proof of Theorem 2.5, we obtain Summing over 1 ≤ n ≤ m for all m + 1 ≤ N and substituting back into (2.27), we obtain , and the desired result readily follows from Theorem 2.8.
We next proceed to derive and analyse a posteriori error estimators. We split the presentation into three cases of increasing complexity, starting with an estimator focusing on the steady coupled problem.
Next, we decompose the H(div)-conforming velocity approximation uniquely into Proof. It follows straightforwardly from the decomposition u h = u c h + u r h and from the edge residual. Lemma 3.2. If u 1,∞ < M , s ∞ < M and c ∞ < M , then the following estimate holds: Proof. Using u h = u c h + u r h and the triangle inequality imply

Owing to the relation
, while using the properties and we readily see that after stating the discrete problem as and employing (3.4), the sought results follow. (3.5)

Proof. Using integration by parts
gives where we define the terms Applying the Cauchy-Schwarz inequality to T 1 implies Next, we rewrite T 2 in terms of a sum over interior edges and apply again the Cauchy-Schwarz inequality. Then Then, owing to the Cauchy-Schwarz inequality, it follows that Proceeding similarly, we may establish the following bounds for T 4 and T 5 : C ∇φ 0,Ω , where C > 0 is a constant independent of h.
Proof. It suffices to apply Lemmas 3.2 and 3.3.

3.2.
Efficiency. For each K ∈ T h , we can define the standard polynomial bubble function b K . Then, for any polynomial function v on K, the following results hold: where C is a positive constant independent of K and v.
Lemma 3.4. The following estimates hold, where C is a positive constant: Moreover, it also follows that Then, using (3.10), we have Using the Cauchy-Schwarz inequality and (3.10) we obtain , and combining these bounds leads to the first stated result. The other two bounds follow similarly.
Let e denote an interior edge that is shared by two elements K and K . Let ω e be the patch which is the union of K and K . Next, we define the edge bubble function ζ e on e with the property that it is positive in the interior of the patch ω e and zero on the boundary of the patch. From [36], the following results hold: Moreover, we also have Proof. Let e be an interior edge and let us define a rescaling of the edge bubble function in the form Using (3.11) gives (3.12) Using integration by parts on each element of patch ω e implies Note that (u, p, s, c) solves the underlying problem, so we then have Next, applying the Cauchy-Schwarz inequality together with Lemma 3.4 and (3.11) gives Combining the bounds of T 1 , T 2 and T 3 with (3.12) and (3.13) implies the first stated result. Similarly, we can prove the other two bounds. Let Ψ be defined as in (3.3). Then there exists a constant C > 0 that is independent of h such that Proof. Combining Lemmas 3.4 and 3.5 implies the stated result.

A posteriori error bound for the semidiscrete method
For each t ∈ (0, T ], let us consider the problem: find (ũ,p,c,s) Also, for each t ∈ (0, T ], we write the discrete weak formulation: where (4.1) remains in effect.
Next we introduce the error indicator Θ as where whereas Ψ is the global a posteriori error estimator for the steady problem with element and edge residual contributions defined in (3). In this case we now replace f and f 1 , f 2 by (4.1).  On the other hand, integrating by parts in Lemma 4.1 yields − a 1 (c,ũ, v) − c 1 (e u ;ũ, v), ∀v ∈ H 1 0 (Ω) (∂ t e s , φ) = −a 2 (ρ s , φ) − c 2 (u; ρ s , φ) − c 2 (e u ;s, φ), ∀φ ∈ H 1 0 (Ω), We apply Young's inequality and the definition of the dual norm. Then, we integrate in time the resulting expression. Finally, the second result is a consequence of Theorem 3.2 and (4.4).

A posteriori error analysis for the fully discrete method
In this section, we develop an a posteriori error estimator for the fully discrete problem and focus the presentation on the simpler case of a time discretisation by the backward Euler method. For each time step k (1 ≤ k ≤ N ), we define the (global in space) time indicator Ξ k as Next we define the accumulated time and spatial error indicators as where the terms Υ 2 k are constructed with the a posteriori error estimator contributions defined as in the steady case (3), but at a given time step k. That is, For each time step k, we can split again the H(div)-conforming discrete solution u k h into a conforming part u k hc and a non-conforming part u k hr such that u k h = u k h,c + u k h,r . For each t ∈ (t k−1 , t k ], we introduce a linear interpolant u h (t) in terms of t as where {l k , l k+1 } is the standard linear interpolation basis defined on [t k , t k+1 ]. Similarly, we may introduce u h,c (t) and u h,r (t). Then, setting e uc = u − u h,c , we have e u = u − u h = e uc − u h,r . For t ∈ (t k−1 , t k ), we define and for all t ∈ (t k−1 , t k ), we consider the problem of finding (ũ k ,p k ,s k , Lemma 5.1. The following estimates hold Proof. Combining (1.2) and (5.2) implies Choosing v = e uc τ , q = p − p k , φ = e s τ , ψ = e c τ and then combining the first two equations, we have (∂ t e uc τ , e uc τ ) + a 1 (c, e uc τ , e uc τ ) + c 1 (u, e uc τ , e uc τ ) − (∂ tũh,r , e uc τ ) + a 1 (c − c h ,ũ, e uc τ ) + c 1 (u − u h , u, e uc τ ) + a 1 (c,ũ h,c −ũ k , e uc τ ) + c 1 (u,ũ h,c −ũ k , e uc τ ) = 0, (∂ t e s τ , e s τ ) + a 2 (e s τ , e s τ ) + c 2 (u; e s τ , e s τ ) − c 2 (u − u h ;s k , e s τ ) + a 2 (s h −s k , e s τ ) − c 2 (u; s h −s k , e s τ ) = 0, where In light of the definition of uτ , sτ and cτ , we get Then we can apply triangle inequality, which gives Finally, applying integration by parts in (5.3a) yields Next we apply Young's inequality and the definition of the dual norm. Then, we integrate the whole expression in time between t k−1 and t k for each k = 1, 2, . . . , n and sum the expression for each k. Finally, we use (5.4), (5.5), (5.6) and (5.7) to establish the second stated result.
Proof. Using u k h = u k h,c + u k h,r together with the identity in [24, (5.59)-(5.60)] that in our context reads , we can invoke Lemma 5.1 and reuse the strategy applied in Theorem 4.1 to complete the proof.
As u is prescribed everywhere on ∂Ω, for sake of uniqueness we impose p ∈ L 2 0 (Ω) through a real Lagrange multiplier approach. To verify the a priori error estimates, we introduce the discrete norms where e,ẽ denote errors generated on two consecutive pairs of mesh size and time step (h, ∆t), and (h,∆t), respectively. Choosing ∆t = √ 2h and using scheme (2.4), the results in Table 6.1 confirm that the rates of convergence are optimal, coinciding with the theoretical bounds anticipated in Theorem 2.8. Approximate velocity magnitude (after 3 refinement steps), pressure (after 4 refinement steps), concentration s (after 5 refinement steps), and distribution of c after 6 steps of adaptive refinement.
6.2. Example 2: adaptive mesh refinement for the stationary problem. The classical Dörfler strategy [23] is employed for the adaptive algorithm based on the steps of solving, estimating, marking, and refining. Estimation is performed by computing the error indicators and using them to select/mark elements that contribute the most to the error [30]. The marking is done following the bulk criterion of selecting sufficiently many elements so that they represent a given fraction of the total estimated error. That is, one refines all elements K ∈ T h for which where 0 < γ ratio < 1 is a user-defined constant (that we tune in order to generate a similar number of degrees of freedom, or comparable errors, as those obtained under uniform refinement). And then the algorithm aims for equidistribution of the local error indicator in the updated mesh.
In the adaptive case, instead of (6.1) the convergence rates (for the spatial errors) are computed as where DoF and DoF are the number of degrees of freedom associated with each refinement level. The robustness of the global estimators is measured using the effectivity index (ratio between the total error and the indicator) .
These solutions exhibit a generic singularity towards the reentrant corner, and therefore one expects that the error decay is suboptimal when applying uniform mesh refinement. After solving the coupled stationary problem on sequences of uniformly and adaptively refined meshes, the aforementioned behaviour is indeed observed in Table 6.2, where the first part of the table shows deterioration of the convergence due to the high gradients of the exact solutions on the non-convex domain. The results shown in the bottom block of the table confirm that as more degrees of freedom are added, a restored error reduction rate is observed due to adaptive mesh refinement guided by the a posteriori error estimator Ψ. The second-last column of the table also indicates that the effectivity index oscillates under uniform refinement, while it is much more steady in the adaptive case. We tabulate as well the Newton-Raphson iteration count (needed to reach the relative residual tolerance of 1e-6), and this number is also systematically smaller for the adaptive case (about four steps in all instances) than for the uniform refinement case (up to seven nonlinear steps for certain refinement levels). As an example we plot  c(x, y, t) = 1 2 (1 + cos(π/4(xy))) exp(−t), s(x, y, t) = 1 2 (1 + sin(π/2(xy))) exp(−t).
Cumulative errors up to t final are computed as , and the resulting error history, after six steps of uniform mesh refinement, is collected in Table 6.3. To be consistent with the development in Section 5, the numerical verification in this set of tests was carried out using a backward Euler time discretisation. The a posteriori error estimator (5.1) is computed and the effectivity index is also tabulated, showing that the estimator is robust (and confirming the theoretical reliability bound as well as giving an heuristic indication of its efficiency). Note that in this case, since the mesh refinement is uniform, the auxiliary interpolation of the solutions at the last time step on the current mesh is not necessary. The average number of Newton-Raphson iterations required for convergence was 3.2.
6.4. Example 4: simulation of salinity-driven flow instabilities. To illustrate the behaviour of the model and the proposed method, we simulate a salinity-driven flow problem. Similar examples are found in [19] where direct numerical simulations (DNS) are applied to the version of (1.1) that has constant viscosity.
The configuration of layering in sedimentation is taken from [33]. We consider a rectangular domain with L x = 40, L y = 300 and an initial solid-particle concentration profile that is periodic in the horizontal direction, and periodic with Gaussian noise in the vertical direction: s(x, y, z, 0) = A 0 exp(−z 2 /σ 2 ) + A 1 sin(x),  with initial amplitudes A 0 , A 1 and width σ (see Figure 6.2). For the velocity field, we use non-slip boundary condition in all four walls and we choose ∆t = 0.1. As discussed in [33], simulations at low density ratios are extremely costly because of the large Reynolds numbers of fingering convection. In consequence we choose an initial density ratio R 0 = αs 0z /βc 0,z ≈ 4, and we carry out the simulations on tall, thin domains. Apart from the specifications above, the remaining constant parameters needed in the model take the following values A 0 = 2.86, A 1 = 0.5, σ = 0.35, ν = 1 × 10 −3 kg/m 3 , g = 9.8m/s 2 , Sc = 7.0, τ = 25, v p = 0.04m/s, α = −2.0, β = 0.5.
According to [33] a linear fingering instability occurs provided 1 < R 0 < τ , hence the instability shown in Figure 6.2 is expected. 6.5. Example 5: adaptive simulation of exothermic flows. To conclude this section, and to include an illustrative simulation exemplifying that the H(div)-conforming scheme along with the a posteriori error estimator perform well for an applicative problem, we address the computation of exothermic flows that develop fingering instabilities. The problem configuration is adapted as a simplification of the problem solved in [31], where the fields c, s represent solutal concentration and temperature, respectively. The model assumes an additional drag term due to porosity so that the momentum equation is of Navier-Stokes-Brinkman type. The domain is the rectangular region Ω = (0, L) × (0, H), and the initial solutal and temperature profiles are imposed as Boundary conditions are of mixed type for solutal and temperature. Both fields are prescribed to 0 and 1 on the bottom and top of the domain, respectively; while on the vertical walls we impose zero-flux boundary conditions. The velocity is of slip type on the whole boundary, and therefore a zero-mean condition for the pressure is considered using a real Lagrange multiplier. The solution algorithm, differently from the previous tests, is based on an inner fixed-point iteration between an Oseen and a transport system, rather than an exact Newton-Raphson method. An initial coarse mesh of 5300 elements is constructed, and an adaptive mesh refinement (only one iteration) guided by the estimator (5.1) is applied at the end of each time step. Figure 6.3 shows snapshots of adapted meshes at different times, and also samples of solute concentration, temperature distribution, velocity, and pressure at the final time.


Trajectory Similarity Analysis with the Weight of Direction and k -Neighborhood for AIS Data

: Automatic Identiﬁcation System (AIS) data have been widely used in many ﬁelds, such as collision detection, navigation, and maritime trafﬁc management. Similarity analysis is an important process for most AIS trajectory analysis topics. However, most traditional AIS trajectory similarity analysis methods calculate the distance between trajectory points, which requires complex and time-consuming calculations, often leading to substantial errors when processing AIS trajectory data characterized by substantial differences in length or uneven trajectory points. Therefore, we propose a cell-based similarity analysis method that combines the weight of the direction and k -neighborhood (WDN-SIM). This method quantiﬁes the similarity between trajectories based on the degree of proximity and differences in motion direction. In terms of its effectiveness and efﬁciency, WDN-SIM outperformed seven traditional methods for trajectory similarity analysis. Particularly, WDN-SIM has a high robustness to noise and can distinguish the similarities between trajectories under complex situations, such as when there are opposing directions of motion, large differences in length, and uneven point distributions. highest. The results indicate that WDN-SIM can sufﬁciently distinguish the similarities between trajectories with different spatial relationships.

Introduction
Rapid developments in wireless communication technology and continuous improvements to positioning accuracy have led to significant progress in terms of the collection, analysis, and application of trajectory data. Trajectory data contains substantial information, which provides strong support for determining motion patterns and feature extraction. In the maritime transport field, the Automatic Identification System (AIS) is a new type of digital navigation aid system and maritime safety equipment that can record and transmit the position, heading, speed, and other information of a ship in real time [1,2]. AIS has a high spatiotemporal resolution, large data volumes, and covers the ports and seas in most global regions [3,4]. It has been used by an increasing number of studies to analyze various maritime traffic problems [5][6][7].
With the aid of cluster analysis, neural networks, association analysis, feature analysis, and other data mining technologies, massive quantities of AIS data can be analyzed to potentially extract useful rules from chaotic ship trajectory points [8][9][10][11][12]. This has important implications for various applications, including maritime safety [13,14], vessel destination prediction [15,16], collision risk identification [17,18], and maritime traffic management [19]. Among the currently popular ship trajectory analysis topics (e.g., trajectory classification, clustering, and trajectory anomaly detection, etc.), calculating the similarity between trajectories is one of the most important and basic processes [20][21][22]. As different data types or application scenarios require distinct similarity analysis methods, the trajectory Figure 1. (a) Schematic illustration of the process to represent trajectory points in real space as a two-dimensional (2-D) grid cells (including neighborhoods). (b) Problem associated with the unreasonable weights set using the C-SIM method for neighborhoods. In (b), T1, T2, and T3 represent three trajectories; using the C-SIM method, the similarity result between T1 and T2 is the same as that between T1 and T3.
To solve these two problems, we propose a cell-based similarity analysis method that combines the weight of the direction and k-neighborhood methods (WDN-SIM). Similar to C-SIM, WDN-SIM calculates the trajectory similarity based on the number of overlapping cells, but WDN-SIM additionally considers the movement direction relationship between two trajectories at the grid cell scale. Moreover, based on distance attenuation, WDN-SIM creates multi-level neighborhoods for the cells and sets different weights for neighborhoods at different levels. WDN-SIM also considers the overall characteristics of the trajectory and converts the similarity measurement result to the range of [−1, 1] to facilitate the comparison of similarity.
The remainder of this paper is organized as follows. Section 2 introduces the theoretical basis and principles of the WDN-SIM similarity analysis method. Section 3 presents a series of experiments to show the effectiveness and advantages of WDN-SIM. Section 4 Figure 1. (a) Schematic illustration of the process to represent trajectory points in real space as a two-dimensional (2-D) grid cells (including neighborhoods). (b) Problem associated with the unreasonable weights set using the C-SIM method for neighborhoods. In (b), T 1 , T 2 , and T 3 represent three trajectories; using the C-SIM method, the similarity result between T 1 and T 2 is the same as that between T 1 and T 3 .
As an easily understandable and implementable method, C-SIM is based on the number of overlapping cells, thus avoiding the complicated distance calculation process. The computational load significantly decreases because the number of representative cells in the trajectory is usually less than the number of trajectory points. Most importantly, the use of grid cells instead of trajectory points reduces the impact that the distribution of sampling points has on the results, which addresses the issue of uneven sampling of ship trajectories. These are also the main reasons why we use cell sequences to represent trajectories in this study.
However, two challenges for similarity analyses of ship trajectories still remain. (1) The influence of diverse shipping directions on trajectory similarity is ignored. C-SIM only considers the spatial location of the track while ignoring the direction of movement. Two ships with similar routes may move in opposite directions which should not be treated as being similar trajectories. However, C-SIM fails to deal with this situation reasonably. (2) The equal and fixed importance of neighborhood cells does not adapt to complex analysis. C-SIM considers neighborhood cells and target cells to be equally important and fixed to eight neighborhoods. However, for tracks T 1 , T 2 , and T 3 , as shown in Figure 1b, the similarity (T 1 , T 2 ) = similarity (T 1 , T 3 ) can be calculated. It is obvious that the distance between T 1 and T 2 is less than the distance between T 1 and T 3 . Therefore, C-SIM is not suitable for analyzing trajectory similarity in the complex situation with high precision requirements.
To solve these two problems, we propose a cell-based similarity analysis method that combines the weight of the direction and k-neighborhood methods (WDN-SIM). Similar to C-SIM, WDN-SIM calculates the trajectory similarity based on the number of overlapping cells, but WDN-SIM additionally considers the movement direction relationship between two trajectories at the grid cell scale. Moreover, based on distance attenuation, WDN-SIM creates multi-level neighborhoods for the cells and sets different weights for neighborhoods at different levels. WDN-SIM also considers the overall characteristics of the trajectory and converts the similarity measurement result to the range of [−1, 1] to facilitate the comparison of similarity.
The remainder of this paper is organized as follows. Section 2 introduces the theoretical basis and principles of the WDN-SIM similarity analysis method. Section 3 presents a series of experiments to show the effectiveness and advantages of WDN-SIM. Section 4 discusses, in detail, the relevant parameters that affect the measurement results. Finally, Section 5 presents our conclusions and directions for future research.

Related Research
Researchers have recently developed many methods for trajectory similarity analyses. One of the most common methods is DTW [28], which obtains the cumulative distance between all optimally matched trajectory points through an iterative approach, thereby allowing the local expansion and contraction of the trajectory. DTW is widely used in AISrelated data analysis. For example, Li et al. [29] applied DTW to the robust ship trajectory clustering. Zhao et al. [30] improved DTW by considering the direction of trajectory point motion and the weight of endpoints. Liu et al. [31] combined DTW with the adaptive Douglas-Peucker (ADP) algorithm to improve the efficiency and accuracy of similarity measurements. However, the DTW method is sensitive to noise.
To overcome this noise sensitivity, Lachos et al. [32] combined the LCSS model with trajectory similarity measurements. The LCSS method quantifies the distance between 0 and 1. By setting a distance threshold, track point pairs, whose distance is less than the threshold, are added to the common trajectory sequence. The trajectory similarity is the distance for the longest common subsequence. LCSS has no strict requirements for the number of trajectory points or trajectory lengths and exhibits good robustness. However, the result is greatly affected by the distance threshold [33]. Furthermore, the EDR method was originally used to calculate the smallest number of operations (i.e., addition, deletion, and modification) required between two strings for complete consistency. This method has been widely used and was extended to similarity analyses of spatiotemporal trajectories [34][35][36][37][38]. However, similar to LCSS, this method cannot easily determine the distance threshold; LCSS and EDR only consider the similar or different components, respectively, such that the results of the similarity analysis are not ideal with uneven sampling point distributions or significant variations in the number of track points.
The methods mentioned above are all based on warping distance, but other studies have also proposed methods based on shape distance [24]. The Hausdorff distance method is a typical method based on shape distance. It calculates the maximum value of the shortest distance from a point on one trajectory to all points on the other trajectory. Wang et al. [39] found that the Hausdorff distance has an optimal effect when measuring the shape similarity between trajectories. However, as the Hausdorff distance must calculate the minimum value of the distance between a point and all points of the other trajectory, and does not consider the time sequence, the calculation efficiency is low and cannot identify the direction of motion. Therefore, Zhen et al. [22] proposed an improved Hausdorff distance method, which expresses the directional distance as the absolute value of the difference between the average course of two tracks. In some cases, this method can distinguish the trajectories of different motion directions but determining the weight of the directional distance remains difficult. Additionally, there are significant variations in the directions of movement when a ship moves at sea because of the influences of wind and seawater. Therefore, using the average direction of movement of all trajectory points as the direction of movement for the entire trajectory cannot yield the true direction of the ship at each point. Ma et al. [40] introduced OWD for the similarity calculation of ship trajectories when conducting research on ship motion pattern recognition. OWD calculates the average of all of shortest distances between two trajectory points [41]. OWD can reduce sensitivity to noise points, but still cannot identify the direction of movement.
Among shape-based methods, another commonly used technique is the Fréchet distance, which originates from the problem associated with the shortest leash length required when walking a dog, i.e., the shortest distance required for the intersection of two curves. Therefore, unlike the Hausdorff distance and OWD, the Fréchet distance considers the time relationship between the trajectory points, such that the resulting trajectories in the same and opposite directions are highly different [42][43][44]. However, as the distance between curves is difficult to obtain, the discrete Fréchet distance is now more commonly used, i.e., calculating the maximum value of the minimum distance between discrete point pairs [45,46].

Overall Idea
The method proposed in this study divides the 2-D space where the trajectory is located into regular grids while combining the spatial neighborhood and direction of the trajectory to analyze the similarity between trajectories. Figure 2 shows the process of this method, which was divided into three steps.   Quantify the direction and neighborhood of the trajectory. We assigned corresponding weights to various directional relationships for different trajectories on the same grid cell. The directional relationships included three types: same direction, inclined direction, and opposite direction. Meanwhile, different neighborhoods of the central cells were also given corresponding weights according to the degree of proximity to the central cell.

3.
Calculate similarity between trajectories. The similarity between the trajectories was measured by calculating the number and proportion of overlapping cells between representative trajectories, followed by assigning corresponding weights to the kneighborhood and motion direction characteristics.

Reconstructing the Representative Trajectory Based on Cell
A ship trajectory is a set of points arranged in chronological order by a series of trajectory points with time labels [23], which can be expressed as P = p 1 , p 2 , p 3 , p 4 , . . . , p n , where n is the number of track points in a certain track and p n is the nth track point present in an AIS record. The attributes of an AIS record include the ship name, call sign, International Maritime Organization (IMO) code, MMSI code, ship type, navigation status, length, width, draft, heading, course, speed, longitude, latitude, destination, estimated time of arrival, and time. The track points for different ships can usually be identified by their unique MMSI code [47]. Reconstructing a representative trajectory mainly includes two steps: trajectory identification and trajectory point mapping.
Trajectory segment identification. Owing to the large time intervals in AIS data, the AIS data for a ship includes all previous voyage information. After obtaining the AIS record of each ship, identifying the records for different voyages is thus necessary. We identified the stopping points in the trajectories to separate different trajectory segments. First, we sorted all trajectory points belonging to the same ship in chronological order, then divided the trajectory points of that ship into sailing and stopping points. The stopping points included situations where the navigation status was docking, anchoring, loading or unloading cargo, stranded, or maintenance, etc., which indicated that either the ship has ceased operation or that the time interval between one point and the following track point exceeds 12 h. The sailing points describe the condition when a ship operates normally at sea. Based on the above rules, we divided the trajectory points of ships into different trajectories according to the positions of stopping points. Figure 3 illustrates a diagram of trajectory identification based on stopping points. As shown in Figure 3b, when a stopping point appears, the original trajectory is divided into two new trajectories by the stopping point. Figure 3a shows multiple consecutive stopping points, in which the first stopping point is the end of the previous trajectory and the last stopping point is the starting point of the next trajectory.
When gridding the trajectory, we used a square structure. Storing and indexing was convenient because the division method for the square grid was more consistent with the latitudinal and longitudinal axes. In contrast, as our method required the calculation of all cells within a certain range surrounding the central cell, the square neighborhood divided the cell into more layers within a smaller distance. ferent trajectories according to the positions of stopping points. Figure 3 illustrates a gram of trajectory identification based on stopping points. As shown in Figure 3b, w a stopping point appears, the original trajectory is divided into two new trajectorie the stopping point. Figure 3a shows multiple consecutive stopping points, in which first stopping point is the end of the previous trajectory and the last stopping point i starting point of the next trajectory. When gridding the trajectory, we used a square structure. Storing and indexing convenient because the division method for the square grid was more consistent with latitudinal and longitudinal axes. In contrast, as our method required the calculatio all cells within a certain range surrounding the central cell, the square neighborhoo vided the cell into more layers within a smaller distance.
For trajectory point mapping, the speed of the ship usually varied; due to uns signals at sea, trajectory points were often missing. Simply assigning the attributes o trajectory points to the grid cells when mapping the trajectory points would produ trajectory with multiple moving directions in the same cell or interruptions in the structed representative trajectory. To solve these problems, we set the following rule As shown in Figure 4a, when a trajectory had consecutive trajectory points pas through the same cell, the average value of the directions of these trajectory points taken as the direction of the entire trajectory on the cell. (2) As shown in Figure 4b, w the representative trajectory was interrupted, the last trajectory point before the inter tion (d6) and the first trajectory point after the interruption (d7) form a line, the cell pas through the midpoint of the line was added to the representative trajectory. The ave value of the directions of d6 and d7 was used as the movement direction in cells with m ing information, until there was no interruption in the trajectory. For trajectory point mapping, the speed of the ship usually varied; due to unstable signals at sea, trajectory points were often missing. Simply assigning the attributes of all trajectory points to the grid cells when mapping the trajectory points would produce a trajectory with multiple moving directions in the same cell or interruptions in the constructed representative trajectory. To solve these problems, we set the following rules. (1) As shown in Figure 4a, when a trajectory had consecutive trajectory points passing through the same cell, the average value of the directions of these trajectory points was taken as the direction of the entire trajectory on the cell. (2) As shown in Figure 4b, when the representative trajectory was interrupted, the last trajectory point before the interruption (d 6 ) and the first trajectory point after the interruption (d 7 ) form a line, the cell passing through the midpoint of the line was added to the representative trajectory. The average value of the directions of d 6 and d 7 was used as the movement direction in cells with missing information, until there was no interruption in the trajectory.

Weight of Direction and Neighbor Cell
In this section, we introduce the method for measuring the direction relations a neighborhoods, which is not considered in C-SIM; however, these are important fact that affect the results of ship trajectory similarity measurements. As we calculated

Weight of Direction and Neighbor Cell
In this section, we introduce the method for measuring the direction relations and neighborhoods, which is not considered in C-SIM; however, these are important factors that affect the results of ship trajectory similarity measurements. As we calculated the similarity between the trajectories based on the number of overlapping cells, we quantified the relationship between the direction of motion between two trajectories in overlapping cells and the different levels of neighborhoods for each cell. This was then added to the final similarity calculation function in the form of a weight.

Weight of Direction
As the course of each point was recorded in the AIS data, we used it to indicate the direction of the track point. We used the average course for all points in the same cell as the trajectory direction in the cell. To obtain the movement direction relationship of two trajectories for each overlapping cell, we used the absolute value of subtracting their directions in the cell, as represented by T (T ∈ [0 • , 180 • ]). The larger the value of T, the greater the difference in direction.
To avoid excessive interference from the direction, as the direction of ship movement was unstable, we simplified the directional relationship between the trajectories into the following three relationships: same direction, inclined direction, or opposite direction. We divided T at equal intervals, using 60 • as the boundary, to allocate equal proportions to the three relations. (1) The same direction, with a T range of [0 • , 60 • ], represents the highest directional similarity, such that the highest weight of 1 was assigned. (2) The inclined direction, with a T range of (60 • , 120 • ], represents a relatively high direction similarity, such that the weight of 0.5 was assigned. (3) The opposite direction, with a T range of (120 • , 180 • ], represents the lowest direction similarity, such that the lowest weight of −1 was assigned. The weight was set to a negative value because an opposite directional relationship weakens the overall similarity. The mathematical expression is as follows: where w 1 D , w 2 D , and w 3 D are the weights corresponding to different directional relationships between the trajectories.

Weight of Neighbor Cells
As the distance increases, the degree of similarity between the trajectories slowly decreases, but due to the arbitrary division of grid cells, adjacent points may be assigned to different cells, which makes the similarity directly change from 1 (same) to 0 (completely different). Therefore, to fully consider the influence of the spatial relationship between trajectories and improve the rationality of the similarity measurement, we added the k-neighborhood of the trajectory grid cell, where k is the maximum number of neighborhoods that participate in the calculation. We set different weight values for neighborhoods at different distances for smoother change in similarity. Figure 5a shows the division method, where "a" is the central cell to which the neighborhood belongs and the remaining cells of the same color represent the same level of neighborhood for "a". The number in the cell corresponds to the level of the neighborhood: cell 1 is the nearest neighbor domain for central cell and cell 2 is the second nearest neighbor domain, extending outward in turn. The number of neighborhood levels involved in the calculation can be determined based on the size of the study area and the density of the trajectory points. at different distances for smoother change in similarity. Figure 5a shows the division method, where "a" is the central cell to which the neighborhood belongs and the remaining cells of the same color represent the same level of neighborhood for "a". The number in the cell corresponds to the level of the neighborhood: cell 1 is the nearest neighbor domain for central cell and cell 2 is the second nearest neighbor domain, extending outward in turn. The number of neighborhood levels involved in the calculation can be determined based on the size of the study area and the density of the trajectory points. To determine the weights of neighborhoods at different levels, as a greater distance between cells results in a lower similarity, we introduced Inverse Distance Weighted method (IDW). We used different neighborhood levels as "distances" and the largest neighborhood level k participating in the calculation as the "distance threshold" to assign corresponding weights to the neighborhoods of each level. The higher the neighborhood level, the smaller the weight. This method can be expressed as follows: where is the weight of the k-neighborhood and k is the level of the neighborhood. According to the above formula, the weights of 1-5 levels are 0.707, 0.577, 0.500, 0.447, 0.408, respectively.

Measuring Similarity between Trajectories
The key to trajectory similarity analysis is the determination of the distance. We used the number of overlapping cells between representative trajectories as the "trajectory distance", assigning corresponding weights to different movement direction relationships and neighborhood levels for each overlapping cell (see Section 3.3 for the weight calculation rules). The calculations for this method included two steps: (1) calculate the degree of overlap between two cell sequences and (2) calculate the proportion of the overlapping part of the trajectories in the total length. Finally, we multiplied the results of these two steps.
When calculating the degree of overlap, to make the similarity calculation results more comparable, we divided the number of overlapping cells, considering the weight of To determine the weights of neighborhoods at different levels, as a greater distance between cells results in a lower similarity, we introduced Inverse Distance Weighted method (IDW). We used different neighborhood levels as "distances" and the largest neighborhood level k participating in the calculation as the "distance threshold" to assign corresponding weights to the neighborhoods of each level. The higher the neighborhood level, the smaller the weight. This method can be expressed as follows: where w k N is the weight of the k-neighborhood and k is the level of the neighborhood. According to the above formula, the weights of 1-5 levels are 0.707, 0.577, 0.500, 0.447, 0.408, respectively.

Measuring Similarity between Trajectories
The key to trajectory similarity analysis is the determination of the distance. We used the number of overlapping cells between representative trajectories as the "trajectory distance", assigning corresponding weights to different movement direction relationships and neighborhood levels for each overlapping cell (see Section 3.3 for the weight calculation rules). The calculations for this method included two steps: (1) calculate the degree of overlap between two cell sequences and (2) calculate the proportion of the overlapping part of the trajectories in the total length. Finally, we multiplied the results of these two steps.
When calculating the degree of overlap, to make the similarity calculation results more comparable, we divided the number of overlapping cells, considering the weight of the direction and neighborhood, by the total number of overlapping cells. This result was converted to the range of [−1, 1], where 1 indicates that the two trajectories were identical and in the same direction, −1 indicates that the two trajectories were identical and in opposite directions, and 0 indicates that the two trajectories were irrelevant; the closer the result is to 0, the lower the degree of similarity.
When the total number of cells in the trajectories was identical, if the ratio of overlapping cells in the two trajectories to the total cells was high, the similarity between them was also high. Therefore, we considered the proportion of the overlapping part of the trajectories to the total length. This is achieved by calculating the ratio of the number of overlapping cells to the total number of cells in the two trajectories, where the similarity was expressed as follows: where S is the similarity for the two trajectories involved in the calculation; i is the i th overlapping cell between the trajectories; n is the total number of overlapping cells between the trajectories; w j N is the weight of the j-neighborhood and j is the level of the neighborhood; w 1 D , w 2 D , and w 3 D are the weights corresponding to different directional relationships between the trajectories; c ij 1 is the number of motion directions at the same interval in the j-neighborhood of the i th overlapping cell between the two trajectories; c ij 2 is the number of motion directions in adjacent intervals to the j-neighborhood of the i th overlapping cell between the two trajectories; c ij 3 is the number of motion directions in opposite intervals in the j-neighborhood of the i th overlapping cell between the two trajectories; C all is the total number of cells in the two trajectories; and C overlap is the number of overlapping cells.

Experimental Dataset
The dataset used in this study was AIS data purchased from the Shipping News Network of Elane, Inc. The collection period was from 1 to 7 January 2015, and the range was 94 • E-127 • E and 6 • S-26 • N. We first processed the data and deleted records with incorrect MMSI codes, incomplete attributes, and repeated records [3]. The ship type was then set to cargo. The data used in experiments to quantitatively analyze the effectiveness and efficiency of the WDN-SIM method were selected from this AIS dataset.

Trajectory Similarity of Different Positional Relationships
The experiment designed herein verified the rationality of the proposed method (WDN-SIM) by evaluating the calculated similarity results between the trajectories of different position relationships in real application scenarios. As WDN-SIM can distinguish whether the trajectory movement direction is the same or opposite, this experiment specifically compared the two cases. Furthermore, this experiment compared and analyzed the similarity of different parts of a trajectory to compare the calculation results of different position relationships using the proposed method. Figure 6 shows the trajectory data selected for this experiment, where Trajectories 1, 2, and 3 had 324, 265, and 310 trajectory points, respectively, and the average sampling interval was approximately 500 m; the direction of movement between trajectories 1 and 2 was the same, while the direction of movement between trajectories 1 and 3, and between trajectories 2 and 3 were opposite; A, B, and C were the further divided areas; the arrows represented the direction of movement of the trajectory.

Comparisons with Other Similarity Measurement Methods
In this section, we used four trajectory transformation experiments and compared the proposed method (WDN-SIM) with several typical ship trajectory similarity measurement methods to verify the effectiveness and efficiency of WDN-SIM [21,48]. The methods used for comparison included DTW, EDR, LCSS, the Fréchet distance, the Hausdorff distance, OWD, and HC-SIM.
Parameter Settings. The parameters involved in the WDN-SIM method included the cell length (L) and maximum neighborhood level (k). By achieving an optimal balance between the trajectory compression rate and feature point retention, after numerous comparison experiments, we set L and k as 2 km and 2, respectively (please see Section 4 for a detailed discussion on the methods of parameter setting and their influence on the results). To maintain consistency, the distance threshold of the LCSS and EDR was correspondingly set as 2 km. After testing, we found that grid cell sizes ≥24 km yielded similarity results with negligible changes, so the minimum cell length of HC-SIM was set to 0.75 km. Unless explicitly stated, the parameters in all experiments were set to the above values (Table 1). Additionally, except for WDN-SIM and HC-SIM, several other methods required the calculation of distances. Haversine is a formula that obtains the arc distance between two points on longitude and latitude; this method can approximate the shortest distance between points on the surface of Earth [49], allowing the measurement of distance as follows: where d is the Haversine distance between two points; r is the approximate radius of Earth (6371 km); ∆Lat denotes subtraction between the latitudes of two points; ∆Lon denotes subtraction between the longitudes of two points; and Lat 1 and Lat 2 are the latitudes of the two points, respectively. Effectiveness. Our experiments used 20 adjacent trajectory pairs in different regions from the AIS dataset, taking the average value as the final result. The average number of points for each trajectory was 350, and the sampling interval was about 500 m. The effectiveness of WDN-SIM was tested by comparing the measured results of the different methods in various transformations (including changing sampling, adding noise, changing direction, and deleting endpoints). Owing to the lack of benchmarks for the similarity of two trajectories, the measurement results for the various methods had different ranges of value and significance; for example, LCSS, EDR, HC-SIM, and WDN-SIM obtained comparable degrees of similarity, whereas DTW, Fréchet, Hausdorff, and OWD only calculated the distance, such that we applied the change rate in the results for comparison: where R is the change rate of similarity measurement result, r is the result before a change, and r is the result after a change.
The four types of transformations methods applied in this study are described in detail as follows.
Changing sampling rate. It includes two types: increasing sampling rate and decreasing sampling rate. Increasing sampling rate was achieved by randomly adding the midpoint, p m , of two continuous sampling points, p i and p i+1 , where the coordinates, time, speed, and course of p m were the mean of p i and p i+1 . Decreasing sampling rate was achieved by randomly deleting a portion of trajectory points from the original trajectory. In this experiment, the sampling rate range was set to 50-150%, where the robustness of the method to a change in sampling rate was evaluated by comparing the changing rate of the measurement result. Additionally, to avoid deleting important points and causing large errors, we used the Threshold-guided Sampling method, proposed by Zhang et al. [50], to extract feature points. The feature points of the trajectory were not changed when adding or deleting points. The specific method was as follows. First, the trajectory points in a trajectory were sorted in chronological order. Changes in the speed between two continuous trajectory points were then calculated, in addition to the course. Finally, based on experience, trajectory points with direction changes >5 • or speed changes >2 knots were set as feature points.
Adding noise. Noise points were added by randomly moving 1-10% of the original trajectory points in the range of 2-10 km. The change rate reflected the sensitivity of different methods to noise points.
Changing direction. In this experiment, we rotated one of the two trajectories moving in the same direction by 90 • and 180 • , where the direction of rotation was clockwise and the center of rotation was the midpoint of the trajectory. The course was also increased by 90 • and 180 • , such that the movement direction relationship between the trajectories after rotation became vertical and reverse, respectively. By comparing the similarity measurement results before and after the transformations, we tested the recognition ability of the different methods with respect to the three movement direction types.
Deleting endpoints. We deleted 5-20% of the points from one end of a trajectory; the resulting change rate was used to measure the ability of the method to identify the length differences between trajectories.
To sum up, during the experiment of changing sampling rate or adding noise, if the change rate after transformation was smaller, the method was less affected by such transformation, hence the stability was higher. However, during the experiment of deleting endpoints, if the change rate after transformation was smaller, the ability to recognize similarities between trajectories of different lengths was weaker.
Finally, we measured the similarity of trajectory pairs containing different numbers of trajectory points. The efficiency of the method was evaluated by comparing the time required for each method to obtain the similarity between two trajectories.  Figure 7 suggests that the similarity between trajectories 1 and 2 was the largest and positive, indicating that the similarity between trajectories 1 and 2 was the highest and the movement direction was identical overall. The similarity results of trajectories 3 and 1, and 3 and 2 were negative because the movement directions between the two trajectories were opposite. The results show that the WDN-SIM method can distinguish the similarity between trajectories with the same and opposite movement directions based on the positive and negative values of the measurement results. From the results of each part, the proximity of different parts between the trajectories was different and the similarity results were also different. From the similarity measurement result of trajectories 2 and 3 in part A (Figure 7b), the similarity between the mutually perpendicular trajectories tended to be 0. Comparing the similarity results in part B between each trajectory pair, the closer the trajectories, the greater the absolute value of the similarity results. In part C, the three trajectories were the closest to each other and the similarity result value was also the highest. The results indicate that WDN-SIM can sufficiently distinguish the similarities between trajectories with different spatial relationships. the positive and negative values of the measurement results. From the results of each part, the proximity of different parts between the trajectories was different and the similarity results were also different. From the similarity measurement result of trajectories 2 and 3 in part A (Figure 7b), the similarity between the mutually perpendicular trajectories tended to be 0. Comparing the similarity results in part B between each trajectory pair, the closer the trajectories, the greater the absolute value of the similarity results. In part C, the three trajectories were the closest to each other and the similarity result value was also the highest. The results indicate that WDN-SIM can sufficiently distinguish the similarities between trajectories with different spatial relationships.  Figure 8 shows the result of the experiment when changing the sampling rate. WDN-SIM was less affected by the sampling rate because the neighborhood cell had a buffering effect; the cells were filled even when the cells were interrupted. Among other methods, DTW was the most sensitive to changes in the sampling rate because it needed to match all trajectory points when calculating the distance, such that changes in the number of points and distance had a significant effect on the results of this method.  Figure 8 shows the result of the experiment when changing the sampling rate. WDN-SIM was less affected by the sampling rate because the neighborhood cell had a buffering effect; the cells were filled even when the cells were interrupted. Among other methods, DTW was the most sensitive to changes in the sampling rate because it needed to match all trajectory points when calculating the distance, such that changes in the number of points and distance had a significant effect on the results of this method.  Figure 8. Results of the experiment with a changing sampling rate. The horizontal axis is the rang of sampling rate changes, while the vertical axis is the change rate of similarity measurement results.

Results of Measurement Comparison Experiments
As shown in Figure 9, DTW remained the most sensitive to added noise points. OWD was the least sensitive to noise points because it used the mean of the minimum distance the change in the distance caused by a small number of noise points was weakened by th other tracking points. The other methods were not sensitive to the added noise points. As shown in Figure 9, DTW remained the most sensitive to added noise points. OWD was the least sensitive to noise points because it used the mean of the minimum distance; the change in the distance caused by a small number of noise points was weakened by the other tracking points. The other methods were not sensitive to the added noise points.

sults.
As shown in Figure 9, DTW remained the most sensitive to added noise points. OWD was the least sensitive to noise points because it used the mean of the minimum distance the change in the distance caused by a small number of noise points was weakened by th other tracking points. The other methods were not sensitive to the added noise points.   As shown in Figure 10, for two similar trajectories, when the length of one changed, the DTW, OWD, and WDN-SIM methods showed higher acuity than the other methods, indicating that these methods were better at identifying differences in the trajectory length.
As shown in Figure 11, when the trajectory had a small number of points, WDN-SIM was less efficient than other methods, mainly because WDN-SIM searched for cell neighborhoods, and was thus time-consuming. However, when there were many trajectory points, WDN-SIM was more efficient than other methods because they required traversing the trajectory points multiple times. For distance-based methods, more trajectory points require a longer calculation time. Cell-based methods did not require measuring the distance, such that HC-SIM and WDN-SIM were less affected by the increase in the number of trajectory points. Moreover, because HC-SIM needed to calculate the average of the similarity of six different cell lengths, its efficiency was lower than that of WDN-SIM. 1 Units for the DTW, Fréchet, Hausdorff, and OWD methods are km.
As shown in Figure 10, for two similar trajectories, when the length of one changed the DTW, OWD, and WDN-SIM methods showed higher acuity than the other methods indicating that these methods were better at identifying differences in the trajector length. Figure 10. Result of the experiment with deleted endpoints. The horizontal axis is the ratio of the deleted endpoints, while the vertical axis is the change rate of similarity measurement results.
As shown in Figure 11, when the trajectory had a small number of points, WDN-SIM was less efficient than other methods, mainly because WDN-SIM searched for cell neigh borhoods, and was thus time-consuming. However, when there were many trajector points, WDN-SIM was more efficient than other methods because they required travers ing the trajectory points multiple times. For distance-based methods, more trajector points require a longer calculation time. Cell-based methods did not require measurin the distance, such that HC-SIM and WDN-SIM were less affected by the increase in th number of trajectory points. Moreover, because HC-SIM needed to calculate the averag of the similarity of six different cell lengths, its efficiency was lower than that of WDN SIM.

Grid Cell Size Selection Problem
The proposed method extracted the representative trajectory by mapping the trajec tory points on a regular grid cell. Different cell sizes have a significant impact on the ac curacy of trajectory extraction results. If the cell size is too large, the trajectory is exces sively compressed and key features are blurred. In contrast, an excessively small grid cel size increases the number of cells in the representative trajectory and increases the com putational load. Therefore, this section designed an experiment to find the optimal cel size. Specifically, as the selection of the cell size is mainly related to the shape and distri bution density of the trajectory, we set different cell sizes to compare the relationship be tween the trajectory point compression rate and the feature point missing rate. We used the Threshold-guided Sampling method (discussed in Section 4.1.3) to extract the featur points of the trajectory [50]. Additionally, we used cell length to reflect the grid cell size Figure 11. Results of the efficiency experiment. The horizontal axis is the number of points included in the calculated two trajectories, while the vertical axis is the average time spent.

Grid Cell Size Selection Problem
The proposed method extracted the representative trajectory by mapping the trajectory points on a regular grid cell. Different cell sizes have a significant impact on the accuracy of trajectory extraction results. If the cell size is too large, the trajectory is excessively compressed and key features are blurred. In contrast, an excessively small grid cell size increases the number of cells in the representative trajectory and increases the computational load. Therefore, this section designed an experiment to find the optimal cell size. Specifically, as the selection of the cell size is mainly related to the shape and distribution density of the trajectory, we set different cell sizes to compare the relationship between the trajectory point compression rate and the feature point missing rate. We used the Threshold-guided Sampling method (discussed in Section 4.1.3) to extract the feature points of the trajectory [50]. Additionally, we used cell length to reflect the grid cell size. Based on Table 3, for the experimental data in this study, when the cell length was set to 2 km, the two indices reached an optimal balance.

k Setting Problem
To illustrate the influence that the k value has on the results of the similarity measurement, the following comparative experiment was conducted. We chose six adjacent trajectories moving in the same direction for the experiment (as shown in Figure 12). The change in the similarity measurement results was then compared for different k values. To compare the results, a certain trajectory was used as the target trajectory. Figure 12 shows the similarity results between the target trajectory and other trajectories, where the horizontal axis is arranged in descending order according to the degree of proximity. Overall, when considering the neighborhood, the value of the similarity measurement result was significantly increased, while the order of proximity remained unchanged. According to the results between trajectories 1 and 3, 1 and 4, and 3 and 5, some trajectory pairs that did not have a similar relationship when not considering the neighborhood, exhibited low similarity when considering the neighborhood. This shows that using neighbor cells can expand the range of the trajectory similarity comparison. However, the more neighborhood levels considered, the longer the execution time required by the program. Therefore, the running time must be considered when determining the k value. We also note that when k = 2, the neighbor grid involved in the calculation was exactly equivalent to the eight neighborhoods of the central cell, therefore, the diagonal cell was not ignored. Based on the above analysis, k was generally set as 2. Overall, when considering the neighborhood, the value of the similarity measurement result was significantly increased, while the order of proximity remained unchanged. According to the results between trajectories 1 and 3, 1 and 4, and 3 and 5, some trajectory pairs that did not have a similar relationship when not considering the neighborhood, exhibited low similarity when considering the neighborhood. This shows that using neighbor cells can expand the range of the trajectory similarity comparison. However, the more neighborhood levels considered, the longer the execution time required by the program. Therefore, the running time must be considered when determining the k value. We also note that when k = 2, the neighbor grid involved in the calculation was exactly equivalent to the eight neighborhoods of the central cell, therefore, the diagonal cell was not ignored. Based on the above analysis, k was generally set as 2.

Conclusions
We proposed a similarity analysis method for AIS ship trajectories in this study, referred to as WDN-SIM, which addresses the problem of C-SIM, which cannot identify the trajectory direction and yields an unreasonable neighborhood weight. WDN-SIM was compared with several traditional trajectory similarity analysis methods (i.e., DTW, EDR, LCSS, discrete Fréchet distance, Hausdorff distance, OWD, and HC-SIM), which showed that WDN-SIM was comparatively less affected by the trajectory sampling rate and noise points, could identify trajectories of different lengths, and performed well in terms of the efficiency, so WDN-SIM is suitable for the similarity analyses of AIS data.
WDN-SIM can not only recognize the similarity between trajectories with different motion direction relationships, but it can also obtain entirely different similarity results for trajectories with the same, perpendicular, and opposite directions. Existing methods can only identify two of these relationships. Therefore, WDN-SIM can obtain more finely resolved trajectory similarity results, which can improve the accuracy of subsequent trajectory classifications and clustering analyses. Additionally, the similarity results obtained with WDN-SIM ranged from −1 to 1. A fixed range such as this can improve the comparability of the results.
WDN-SIM provides new insights for improving the methods used to quantify the similarities between different ship trajectories. This is important in terms of logistics, shipping efficiency, and course plotting, among other issues, in addition to facilitating trajectory analyses, such as clustering, classification, or anomaly detection. To generalize the proposed method, our future research needs to explore the spatio-temporal similarities in AIS trajectories. In addition, since we judge the movement direction of the trajectory in each grid cell by calculating the average value of the direction, the direction will be inaccurate for complex situations such as zigzag movement or circular movement. Thus, further research should investigate the direction recognition method in complex scenes.
Author Contributions: Conceptualization, Pin Nie, Zhenjie Chen and Nan Xia; methodology, Pin Nie and Zhenjie Chen; validation, Qiuhao Huang and Feixue Li; formal analysis, Qiuhao Huang and Nan Xia; investigation, Zhenjie Chen and Feixue Li; resources, Zhenjie Chen; writing-original draft preparation, Pin Nie, Zhenjie Chen; writing-review and editing, Pin Nie, Zhenjie Chen, Nan Xia, Qiuhao Huang and Feixue Li; funding acquisition, Zhenjie Chen, Qiuhao Huang and Feixue Li. All authors have read and agreed to the published version of the manuscript.
Funding: This work was supported by the National Natural Science Foundation of China (No. 42171396, 42101415). Data Availability Statement: Some or all data or code generated or used during the study are available from the corresponding author by request.

Conflicts of Interest:
The authors declare no conflict of interest.