On the Horizontal Divergence Asymmetry in the Gulf of Mexico

1. Introduction

In large-scale oceanic flows, the Coriolis effect and the pressure gradient dominate the horizontal momentum balance [1,2]. When the Coriolis parameter is a constant, geostrophic flow is horizontally non-divergent. This horizontal non-divergent condition simplifies many theoretical analyses. For instance, the third-order structure function theory with the incompressible condition can be used to study the direction of energy cascades in the atmosphere [3] and ocean [4,5].

Geostrophic balance describes the oceanic mesoscale to the leading order of small Rossby number, while at submeso and smaller scales, horizontal divergence can be of leading order [6] and even mesoscale motion can have weak horizontal divergence. Second-order structure functions can be constructed to detect the amplitude of this compressibility [7], but the sign of horizontal velocity divergence is unknown.

We go beyond the velocity structure function and study the combined information from the velocity difference vector $\delta \mathit{u}$ and acceleration difference vector $\delta \mathit{a}$. In three-dimensional (3D) turbulence, the ensemble average of their dot product, denoted as $\langle \delta \mathit{u}·\delta \mathit{a}\rangle$, is found to be a constant under a fixed energy cascading rate with its sign corresponding to the direction of cascade in the inertial range, which is justified by experiments [8]. Thus, it is natural to study the angle $\theta$ between $\delta \mathit{u}$ and $\delta \mathit{a}$. In a simple model with constant strain rates and without rotation, the probability density function (p.d.f.) of $\cos \theta$ is related to the aspect ratio of the eigenvalues of the velocity strain rate tensor. The p.d.f. of $\cos \theta$ under this model fits well with the p.d.f. from experiments [9].Therefore, it is possible to obtain information about the eigenvalues of the strain rate tensor from the p.d.f. of $\cos \theta$.

This paper treats the oceanic surface flow as a projection of 3D flow. In this sense, the eigenvalues of the strain rate tensor reveal the compressibility of the surface flow. By analyzing drifter data in the Gulf of Mexico, we aim to detect the amplitude and sign of the horizontal divergence in oceanic surface flow.

The structure of the rest of this paper is as follows. Section 2 introduces the theoretical background and methods for analyzing the drifter data. Then we show the main results in Section 3. Finally, we discuss the phenomenon found from the results and propose the possible physics in Section 4. Details of data processing and error estimation will be shown in Appendix A and Appendix B. More fitting results are displayed in Appendix C.

2. Theories and Methods

We build a model that captures the relation between the angle of acceleration and velocity difference vectors and the divergence at a certain scale, which is inspired by Gibert et al. (2012) [9]. Since drifters may be distributed at strain-dominant regions around vortices in a diffusion process [10], we consider that at each scale, random strains dominantly impact pair separation. Thus we keep the model to be irrotational.

2.1. Random Strain Model

We can start from a simple case where the flow has constant strain rates and does not rotate. Consider a 3D incompressible flow, then the velocity gradient tensor $\mathit{M}$ is constant and symmetric with zero trace:

$$M_{ij}={\displaystyle \frac{\partial u_i}{\partial x_j}}={\displaystyle \frac{\partial u_j}{\partial x_i}}=M_{ji}=\mathrm{const}.$$

(1)

For simplicity, we consider a linear velocity profile and the corresponding acceleration can be expressed as

$$\begin{array}{cc}\hfill u_i& =M_{ij}{x}_j,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill a_i& =u_j{\displaystyle \frac{\partial u_i}{\partial x_j}}=M_{ij}{u}_j.\hfill \end{array}$$

(2b)

Here, we chose the Cartesian coordinates O-$xyz$, where the x-axis points to the East, the y-axis points to the North and the z-axis is vertically upwards. The indices $i,j$ range from 1 to 3, and each pair of repeating indices represents the corresponding summation from 1 to 3. The velocity is also denoted as $\mathit{u}=u\mathbf{i}+v\mathbf{j}+w\mathbf{k}$, where $\mathbf{i},\mathbf{j},\mathbf{k}$ represent the unit vectors of the axes $x,y,z$ respectively.

We can perform an orthogonal diagonalization since $\mathit{M}$ is assumed symmetric. Choosing the principal coordinates O-$\xi \eta \zeta$, we obtain

$$\mathit{M}=\left(\begin{array}{ccc}{\lambda }_1& & \\ & {\lambda }_3& \\ & & {\lambda }_2\end{array}\right)={\lambda }_1\left(\begin{array}{ccc}1& & \\ & -(1+\alpha )& \\ & & \alpha \end{array}\right),$$

(3)

where ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3$ are the eigenvalues of $\mathit{M}$. Here, we applied the incompressibility condition

$$\nabla ·\mathit{u}=\mathrm{tr}\mathit{M}=\sum _{i=1}^3{\lambda }_i=0,$$

(4)

and $\alpha ={\lambda }_2/{\lambda }_1$.

For the large-scale oceanic flow, due to the dominant geostrophic balance, the principal axes $\xi$ and $\eta$ lie close to the x-y plane. Therefore, it is reasonable to assume $\left|\alpha \right|\ll 1$ and the $\zeta$-axis is parallel to the z-axis. Under these assumptions, the horizontal velocity divergence is not zero:

$${\nabla }_z·\mathit{u}={\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}=-{\displaystyle \frac{\partial w}{\partial z}}=-{\lambda }_2,$$

(5)

where ${\nabla }_z:=\mathrm{i}{\partial }_x+\mathrm{j}{\partial }_y$ is the horizontal gradient.

Thus, a divergent flow would correspond to a negative ${\lambda }_2$ and a negative $\alpha$, and a convergent flow would correspond to a positive ${\lambda }_2$ and a positive $\alpha$.

Considering two-dimensional (2D) drifter pair dispersion in the ocean, the relative position vector between two drifters can be expressed as

$$\delta \mathit{r}=r_0({\mathbf{e}}_{\xi }\cos \delta \varphi +{\mathbf{e}}_{\eta }sin\delta \varphi ),$$

(6)

where ${\mathbf{e}}_{\xi }$ and ${\mathbf{e}}_{\eta }$ are the unit vectors of the axes $\xi$ and $\eta$, $r_0$ is the distance between the drifters, and $\delta \varphi$ represents the angle between $\delta \mathit{r}$ and ${\mathbf{e}}_{\xi }$. When $\delta \mathit{r}$ is fixed, the velocity and acceleration difference vectors can be expressed as

$$\delta \mathit{u}=\left(\begin{array}{ccc}{\lambda }_1& & \\ & {\lambda }_3& \\ & & {\lambda }_2\end{array}\right)·\delta \mathit{r}=r_0{\lambda }_1\left(\begin{array}{c}\cos \delta \varphi \\ -(1+\alpha )sin\delta \varphi \\ 0\end{array}\right),$$

(7)

and

$$\delta \mathit{a}=\left(\begin{array}{ccc}{\lambda }_1^2& & \\ & {\lambda }_3^2& \\ & & {\lambda }_2^2\end{array}\right)·\delta \mathit{r}=r_0{\lambda }_1^2\left(\begin{array}{c}\cos \delta \varphi \\ {(1+\alpha )}^{2}sin\delta \varphi \\ 0\end{array}\right).$$

(8)

Let $\theta$ be the angle between $\delta \mathit{u}$ and $\delta \mathit{a}$. We can obtain the exact expression of $\cos \theta$:

$$\cos \theta ={\displaystyle \frac{\delta \mathit{u}·\delta \mathit{a}}{\left|\delta \mathit{u}\right|\left|\delta \mathit{a}\right|}}={\displaystyle \frac{{\cos }^2\delta \varphi -{(1+\alpha )}^3{sin}^2\delta \varphi }{\sqrt{{\cos }^2\delta \varphi +{(1+\alpha )}^2{sin}^2\delta \varphi }\sqrt{{\cos }^2\delta \varphi +{(1+\alpha )}^4{sin}^2\delta \varphi }}},$$

(9)

which is independent of $r_0$ and ${\lambda }_1$.

Oceanic flows can be different when the scale changes. To apply the model to real oceanic flows, we assume ${\lambda }_i$ ($i=1,2,3$) and $\alpha$ to be functions of the scale, i.e. functions of $r_0$ in a 2D drifter pair dispersion. We further assume that the angle $\delta \varphi$ is a random variable (r.v.), to capture the changing direction of the principal axes and drifter pairs. When the flow is assumed isotropic, $\alpha$ is independent of $\delta \varphi$. Considering (9), for a fixed scale, the quantity $\cos \theta$ as a function of $\delta \varphi$ is also a r.v. Our aim is to obtain $\alpha$ from the distribution of $\cos \theta$ under different scales.

2.2. Probability Distribution of $\cos \theta$

The exact expression of $\cos \theta$ as a function of $\delta \varphi$ (9) has a complicated form which is not easily applicable to data analysis. As has been discussed above, the dominant geostrophic balance implies $\left|\alpha \right|\ll 1$, so we can expand $\cos \theta$ into a series of $\alpha$:

$$\cos \theta =\sum _{k=0}^{\infty }{\alpha }^k{f}_k\left(\delta \varphi \right).$$

(10)

At the first order, we obtain

$$\cos \theta =\cos 2\delta \varphi +{\displaystyle \frac{3}{2}}\alpha ({\cos }^{2}2\delta \varphi -1).$$

(11)

It is reasonable to reset $\cos 2\delta \varphi$ instead of $\delta \varphi$ to be the independent variable. We introduce the following notations for convenience:

$$x:=\cos 2\delta \varphi ,y:=\cos \theta .$$

(12)

The reverse function of (11) exists when $\left|\alpha \right|\le 1/3$, and the 1st-order relations become

$$\begin{array}{cc}\hfill y& =x-{\displaystyle \frac{3}{2}}\alpha (1-x^2),\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill x& =y+{\displaystyle \frac{3}{2}}\alpha (1-y^2).\hfill \end{array}$$

(13b)

To the second order, these relations become

$$\begin{array}{c}\hfill {\displaystyle y=x-\frac{3}{2}\alpha (1-x^2)+{\alpha }^2(1-x^2)\left(\frac{3}{4}-\frac{19}{8}x\right),}\\ \hfill {\displaystyle x=y+\frac{3}{2}\alpha (1-y^2)-{\alpha }^2(1-y^2)\left(\frac{3}{4}+\frac{17}{8}y\right).}\end{array}$$

(14)

The results of our study will show that $\alpha$ is on the order of ${10}^{-2}$ for oceanic surface flow, making the 1st-order expansion a good approximation.

The probability distribution of $\cos \theta$ can be expressed by the distribution of $\delta \varphi$ (or $\cos 2\delta \varphi$) in an asymptotic form. Consider a simple case when $\delta \varphi$ is uniformly distributed over $[0,2\pi ]$. The p.d.f. of $x=\cos 2\delta \varphi$ was derived in the 2D incompressible case with $\alpha =0$ [9]:

$$P\left(x\right)={\displaystyle \frac{1}{\pi \sqrt{1-x^2}}},x\in (-1,1).$$

(15)

To the 1st order of $\alpha$, We derive the expression of p.d.f. of $\cos \theta$ as

$$\tilde{P}\left(y\right)={\displaystyle \frac{1-{\displaystyle \frac{3}{2}}\alpha y}{\pi \sqrt{1-y^2}}},y\in (-1,1).$$

(16)

This expression introduces an asymmetric part with $\alpha$ as the leading order, which inspires us to describe the influence of non-zero $\alpha$ by even-odd decomposition.

2.3. Procedure for Analyzing Oceanic Surface Flow

Hinted by (16), we decompose the p.d.f.s into the even and odd parts

$$\begin{array}{ccc}\hfill P_{\mathrm{even}}\left(x\right)& ={\displaystyle \frac{P\left(x\right)+P(-x)}{2}},\hfill & \hfill {\tilde{P}}_{\mathrm{even}}\left(y\right)={\displaystyle \frac{\tilde{P}\left(y\right)+\tilde{P}(-y)}{2}},\\ \hfill P_{\mathrm{odd}}\left(x\right)& ={\displaystyle \frac{P\left(x\right)-P(-x)}{2}},\hfill & \hfill {\tilde{P}}_{\mathrm{odd}}\left(y\right)={\displaystyle \frac{\tilde{P}\left(y\right)-\tilde{P}(-y)}{2}}.\end{array}$$

(17)

The analysis will be based on the following assumptions:

Assumption 1. 

Assuming that in the even-odd decomposition $P\left(x\right)=P_{\mathrm{even}}\left(x\right)+P_{\mathrm{odd}}\left(x\right)$, $P_{\mathrm{even}}\left(x\right)\gg P_{\mathrm{odd}}\left(x\right)$. So we can introduce a bookkeeping parameter ϵ to mark the smallness of the odd part of the p.d.f., i.e.

$$P\left(x\right)=P_{\mathrm{even}}\left(x\right)+ϵP_{\mathrm{odd}}\left(x\right),0<ϵ\ll 1.$$

(18)

Assumption 2. 

$\alpha ={\lambda }_2/{\lambda }_1$ is small in oceanic flows, i.e. $\left|\alpha \right|\ll 1$. Hence the discussion in Section 2.2 holds valid.

Assumption 3. 

$\left|\alpha \right|\gg ϵ$, i.e. the odd p.d.f. of $\cos \theta$ is dominated by α instead of ϵ. Thus, α measures the skewness of p.d.f. of the angular distribution of drifter pairs.

The p.d.f. of $y=\cos \theta$ can be derived through a simple relation: $\tilde{P}\left(y\right)=P\left(x\left(y\right)\right)|\mathrm{d}x/\mathrm{d}y|$, where $\mathrm{d}x/\mathrm{d}y=1-3\alpha y$ always holds non-negative under assumption 2. Denoting $x=x_0\left(y\right)+\alpha x_1\left(y\right)$ with $x_0\left(y\right)=y$ and $x_1\left(y\right)=3(1-y^2)/2$, we have:

$$\begin{array}{cc}\hfill \tilde{P}\left(y\right)& =P\left(x\left(y\right)\right){\displaystyle \frac{\mathrm{d}x}{\mathrm{d}y}}=P(x_0+\alpha x_1)\left({\displaystyle \frac{\mathrm{d}{x}_0}{\mathrm{d}y}}+\alpha {\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}y}}\right)\hfill \\ \hfill & =\left[P\left(x_0\right)+\alpha x_1{P}^{\prime }\left(x_0\right)\right]\left({\displaystyle \frac{\mathrm{d}{x}_0}{\mathrm{d}y}}+\alpha {\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}y}}\right)\hfill \\ \hfill & =P\left(x_0\right){\displaystyle \frac{\mathrm{d}{x}_0}{\mathrm{d}y}}+\alpha {\displaystyle \frac{\mathrm{d}}{\mathrm{d}y}}\left[x_1\left(y\right)P\left(x_0\left(y\right)\right)\right]\hfill \\ \hfill & =P\left(y\right)+{\displaystyle \frac{3}{2}}\alpha {\displaystyle \frac{\mathrm{d}}{\mathrm{d}y}}\left[(1-y^2)P\left(y\right)\right],\hfill \end{array}$$

(19)

where we omitted the second-order term. Combined with assumption 1 and assumption 3, we get the even-odd decomposition for $\tilde{P}\left(y\right)$:

$$\tilde{P}\left(y\right)=\underset{\mathrm{even}}{\underbrace{P_{\mathrm{even}}\left(y\right)}}+\underset{\mathrm{odd}}{\underbrace{{\displaystyle \frac{3}{2}}\alpha {\displaystyle \frac{\mathrm{d}}{\mathrm{d}y}}\left[(1-y^2)P_{\mathrm{even}}\left(y\right)\right]}},$$

(20)

where the skewness of p.d.f. of the angular distribution is captured by $\alpha$.

Considering the decomposition in (20), we can use the even part to express the odd part, thus obtaining the $\alpha$ by fitting. Here we use the least square method for fitting.

After obtaining the odd and even parts of the p.d.f. from data, we introduce the loss function as

$$L\left(\alpha \right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left[{\displaystyle \frac{3}{2}}\alpha q^{\prime }\left(y_i\right)-p\left(y_i\right)\right]}^2,$$

(21)

where

$$\begin{array}{cc}\hfill p\left(y\right)& ={\tilde{P}}_{\mathrm{odd}}\left(y\right),\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill q\left(y\right)& =(1-y^2){\tilde{P}}_{\mathrm{even}}\left(y\right),\hfill \end{array}$$

(22b)

then we obtain $\alpha$ by minimizing $L\left(\alpha \right)$:

$$\alpha ={\displaystyle \frac{2{\sum }_{i=1}^{n}p\left(y_i\right)q^{\prime }\left(y_i\right)}{3{\sum }_{i=1}^n{\left[q^{\prime }\left(y_i\right)\right]}^2}}.$$

(23)

Here the sequence ${\left\{y_k\right\}}_{k=1}^n$ represents the data points, and it is chosen to be a finite arithmetic sequence with a common difference of $\Delta y$.

In our model, $\alpha$ is closely related to horizontal compressibility by (5). If $\alpha$ itself is a r.v. with its p.d.f. denonted as $f\left(\alpha \right)$, assuming that $f\left(\alpha \right)$ is negligible when $\left|\alpha \right|>1/3$, then the relation between x and y is 1-to-1, therefore the $\alpha$ in (20) can be understood as its expectation $\int \alpha f\left(\alpha \right)\mathrm{d}\alpha$. Thus, we may expect to obtain a statistical average of horizontal velocity divergence, i.e. the asymmetry towards divergence or convergence, by the sign of $\alpha$.

3. Results

In this study, we analyze two data sets of drifter motion from the Grand LAgrangian Deployment (GLAD) experiment and the LAgrangian Submesoscale ExpeRiment (LASER). Drifters were deployed in the Gulf of Mexico in Summer/July to August 2012 (GLAD) [11] and in Winter/January to February 2016 (LASER) [12]. The locations of drifters are tracked using the Global Positioning System (GPS) within a nominal position error less than 10 m [5]. We think these two sets of data can reflect the velocity of oceanic surface flows precisely at the scales studied here.

The p.d.f. of $\cos \theta$ to the zeroth order is (15) when $\delta \varphi$ is uniformly distributed over $[0,2\pi ]$. Therefore, we can do a rough comparison between the actual p.d.f. of $\cos \theta$ and (15). Choose the data in LASER and $\left|\delta \mathit{r}\right|$ around 33 km as an example, and the result is shown in Figure 1. Details of data processing are written in Appendix A. We find that (15) can already qualitatively describe the p.d.f. of $\cos \theta$, i.e. $\cos \theta$ distributing around $\pm 1$ has the maximum probability, while the p.d.f. changes little and lies below $0.5$ around $\cos \theta =0$. The dimensionless quantity $\alpha$, which may represent the compressibility by (5), appears in the first-order expansion. Then, a rough estimation of $\alpha$ can be done according to (16).

Section 2.3 provides a method to obtain $\alpha$ with a weaker assumption to the distribution of $\delta \varphi$, and we can take an example to show the whole process. Figure 2 shows the fitting process and its performance for $\left|\delta \mathit{r}\right|$ around 33 km in LASER. We first decompose the p.d.f. of $\cos \theta$ into the even and odd part, shown in Figure 2(a) and (c), respectively. Note that in Figure 2(a), q relates to the even part of p.d.f. of $\cos \theta$ through (22b). Then, using the discrete format in Appendix A, we take the derivative of q, shown in Figure 2(b). Thus, we fit p using $q^{\prime }$ to obtain $\alpha$ following the least-square procedure shown in (21)–(23). The red fitting curve in Figure 2(c) well captures p. Figure 2(d) shows the comparison between the p.d.f. of $\cos \theta$ obtained from data and the theoretical expression (20) obtained from our fitting procedure.

The above process can be carried out in other scales of $\left|\delta \mathit{r}\right|$, in GLAD and LASER. More fitting results are shown in Appendix C. Then we can get $\alpha$ as a function of scale, which is shown in Figure 3. Details of data processing and error estimation are written in Appendix A and Appendix B, which will tell us that Figure 3 is reliable when the scale is above 3 km. We can find that $\alpha$ is on the order of ${10}^{-2}$ in general, and stays negative when the scale ranges from about 10 km to about 300 km. This range of scale is larger in LASER than in GLAD, and so does the amplitude of negative $\alpha$.

4. Discussion

Figure 3 shows that $\alpha$ is negative (implying on average divergent flow) at most of the scales, with some positive (convergent flows) $\alpha$ at scales around 1-10km. The negative $\alpha$ is more prominent in LASER than in GLAD. This at least verifies that the surface flow is weakly compressible, as $\alpha$ is small and non-zero, and may involve some patterns in 3D flows. Besides, compared with theories of 3D homogeneous isotropic turbulence [13] with $\langle {\lambda }_1{\lambda }_2{\lambda }_3\rangle <0$ (which may imply a positive $\alpha$)1, the negative $\alpha$ implies the uniqueness of oceanic flow. We do not know the exact reason behind this negative $\alpha$, but we propose in the following sections two possible explanations.

4.1. Kinematics: Drifter Concentration Caused by Mesoscale Vortices and Surface Compressibility

Cressman et al. (2004) [14,15] analysed the motion of small particles on the surface of water in a square tank, and simulated the velocity distribution numerically. They claimed that surface flow is a compressible flow with low speed (lower than the speed of sound). Structures of source and sink will exist on the surface, as the motion is still 3D in essence. Particles will be trapped in narrow convergent regions (with negative 2D divergence of velocity), and be distribute around vortices with a larger scale.

For oceanic flows, the leading order vortices are mesoscale eddies (ranging from several km to several hundreds of km). It is possible that these mesoscale vortices are accompanied with narrow convergent fronts, so that the drifters will accumulate around these vortices [12]. Meanwhile, the vortices correspond to divergent regions due to mass conservation. This qualitative view agrees with the suggestion of divergence (negative $\alpha$) at mesoscales ($\sim 10-100$ km), with some suggestion of convergence (positive $\alpha$) at submesoscale fronts ($\sim 1-10$ km). We further speculate, that the negative alpha seen at even smaller scales may correspond to the divergent zones of Langmuir cells.

4.2. Dynamics: Weak Compressibility Caused by Ekman Pumping

Coriolis force is non-negligible for large scales in oceanic flows. The balance between Coriolis force and internal friction force is significant near the surface, resulting in the boundary layer between the ocean and the atmosphere, which is called the Ekman layer. Some contents in this section will use Vallis (2017) [1] as a reference.

For the Ekman layer on the surface of the ocean, friction at the bottom is zero (from definition), and wind stress ${\mathit{\tau }}_{\mathrm{T}}$ will happen at the top. The z-component w of the velocity stays zero at the top, but non-zero at the bottom:

$$w_{\mathrm{B}}={\displaystyle \frac{1}{{\rho }_0}}\left({\mathrm{curl}}_z{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{T}}}{f}}+\nabla ·{\mathit{M}}_g\right).$$

(24)

Here ${\mathrm{curl}}_z$ denotes the z-component of the curl, ${\mathit{M}}_g={\int }_{\mathrm{Ek}}{\rho }_0{\mathit{u}}_g\mathrm{d}z$, where ${\mathit{u}}_g$ is the horizontal component of the velocity under geostrophic balance. The term $\nabla ·{\mathit{M}}_g=-(\beta /f){\mathit{M}}_g·\mathrm{j}$ represents the divergence of geostrophic transport, which is relatively small in (24).

Therefore, wind stress with a curl consistent with the rotation of the earth (i.e. positive curl in the northern hemisphere) will cause the up-welling and surface divergence, which is likely to contribute to the negative $\alpha$. The Gulf of Mexico is a region associated with negative wind stress curl on average [16], which could not explain the observed positive $\alpha$ (divergence) at scales larger than $\sim 200$ km. However, this region is also subject to hurricanes, which are cyclones with a positive wind stress curl. In particular, the LASER experiment experienced a strong hurricane, which may explain the more negative $\alpha$ detected in the LASER dataset.

5. Conclusions

We constructed a pure-strain model with a parameter $\alpha$ to analyse the drifter data, where $\alpha$ can reflect the existence of compressibility in a 2D flow, and the horizontal divergence asymmetry in 2D turbulence. We find that $\alpha$ stays negative at the mesoscales, representing a divergent flow, and we conjecture that mesoscale eddies interspersed with convergent fronts or a wind stress with a cyclonic wind stress curl arising from hurricanes may explain this observed divergence.

This study analyzed two data sets, and our model better matches LASER than GLAD (shown in Appendix B). This difference may be caused by the difference in data volume, or our assumptions, such as homogeneity and isotropy, are better satisfied in LASER than in GLAD. The inhomogeneity and anisotropy in GLAD were already discussed in other studies [17].

A major caveat of this work is that we assumed the flow to be pure-strain; but oceanic flows are known to have comparably strong vorticity as well [18]. Thus, this strong assumption on our flow model may impact some of the results and conclusions. Also, our model is probably invalid when the velocity gradient tensor changes significantly with time, and the assumption of linearity is restrictive. Regardless, we provide an interesting and novel proof of concept for analyzing ocean observations, which is worthy of further investigation.

Appendix B.1. Noise in the Generated p.d.f.

We generated p.d.f. of $\cos \theta$ from drifter data directly. But these functions are not smooth, and contain irregular oscillation, which is called noise. The existence of noise will reduce the effect and reliability of our fitting. Observation tells us that noise will probably decrease (i.e. Signal to Noise Ratio (SNR) will increase) as the scale and data volume increase. For instance, there are more valid drifters in LASER (956) than in GLAD (297), and accordingly, the p.d.f.s appear smoother in LASER than in GLAD.

Moreover, SNR will also increase when $\Delta y$ appropriately increases. An example is taken for 33 km in LASER, shown in Figure A1. Obvious oscillations can be seen when $\Delta y=0.02$, and this phenomenon vanishes gradually as $\Delta y$ grows up. However, changes in $\Delta y$ will result in changes in the obtained $\alpha$, and the amplitude of this change also depends on the scale and data volume.

We can plot the $\alpha$-$\left|\delta \mathit{r}\right|$ curves at different values of $\Delta y$, which are shown in Figure A2. There should not be big changes in $\alpha$ if our method is proper enough. According to Figure A2, small scales correspond to a bigger amplitude of change in $\alpha$ and thus a lower reliability compared with the larger scales. The demarcation approximately locates at 3 km.

Figure A1. Odd part of the p.d.f. and the fitting curve when $\Delta y$ equals (a) $0.02$, (b) $0.05$ and (c) $0.1$, at the scale around 33 km in LASER.

Figure A1. Odd part of the p.d.f. and the fitting curve when $\Delta y$ equals (a) $0.02$, (b) $0.05$ and (c) $0.1$, at the scale around 33 km in LASER.

Figure A2. The $\alpha$-$\left|\delta \mathit{r}\right|$ curves at different values of $\Delta y$ in (a) GLAD and (b) LASER. The values of $\Delta y$ are set as $0.02$, $0.05$ and $0.1$, plotted in blue, red and black respectively.

Figure A2. The $\alpha$-$\left|\delta \mathit{r}\right|$ curves at different values of $\Delta y$ in (a) GLAD and (b) LASER. The values of $\Delta y$ are set as $0.02$, $0.05$ and $0.1$, plotted in blue, red and black respectively.

Appendix B.2. Root Mean Square Error

In the former section, We found that the fitting is effective and reliable at mesoscale and large scale. But we have not checked whether the fitting matches the p.d.f. well enough, and this is independent from the existence of noise. Therefore, we set $\Delta y=0.1$ to increase SNR, and define the bias of fitting as:

$${\sigma }_{\mathrm{Bias}}=\sqrt{{\displaystyle \frac{1}{n^{\prime }}}\sum _{i=1}^{n^{\prime }}{\left[p\left(y_i^{\prime }\right)-p^{\mathrm{fit}}\left(y_i^{\prime }\right)\right]}^2},$$

(A2)

where $n^{\prime }$ denotes the number of data points $y_i^{\prime }$ when $\Delta y=0.1$, $p\left(y_i^{\prime }\right)$ and $p^{\mathrm{fit}}\left(y_i^{\prime }\right)={\displaystyle \frac{3}{2}}\alpha q^{\prime }\left(y_i^{\prime }\right)$ denote the odd part of the p.d.f. of $\cos \theta$ and the fitting curve. Figure A3 shows the ${\sigma }_{\mathrm{Bias}}$-$\left|\delta \mathit{r}\right|$ curves. The value of ${\sigma }_{\mathrm{Bias}}$ reduces to a stable platform when the scale is over 500 m.

Figure A3. ${\sigma }_{\mathrm{Bias}}$-$\left|\delta \mathit{r}\right|$ curves in (a) GLAD and (b) LASER.

Figure A3. ${\sigma }_{\mathrm{Bias}}$-$\left|\delta \mathit{r}\right|$ curves in (a) GLAD and (b) LASER.

In addition, the amplitude of the odd part of the p.d.f. of $\cos \theta$ will vary as the scale changes. A quantity describing the relative error needs to be introduced, i.e.

$$e_{\mathrm{Bias}}={\displaystyle \frac{{\sigma }_{\mathrm{Bias}}}{\sqrt{\frac{1}{n^{\prime }}{\sum }_{i=1}^{n^{\prime }}{p}^2\left(y_i^{\prime }\right)}}}.$$

(A3)

The $e_{\mathrm{Bias}}$-$\left|\delta \mathit{r}\right|$ curves are shown in Figure A4. We can find that the fitting matches better in LASER than in GLAD, and the error is indeed smaller in larger scales (roughly larger than 3 km). The difference between the two data sets may be resulted by the difference of data volume, or GLAD is less likely to accord with our model.

Figure A4. $e_{\mathrm{Bias}}$-$\left|\delta \mathit{r}\right|$ curves in (a) GLAD and (b) LASER.

Figure A4. $e_{\mathrm{Bias}}$-$\left|\delta \mathit{r}\right|$ curves in (a) GLAD and (b) LASER.