Dynamics of inertial particles on the ocean surface with unrestricted reserve buoyancy

Abstract

The purpose of this note is to present an enhancement to a Maxey–Riley theory proposed in recent years for the dynamics of inertial particles on the ocean surface. This model upgrade removes constraints on the reserve buoyancy, defined as the fraction of the particle volume above the ocean surface. The refinement results in an equation that correctly describes both the neutrally buoyant and fully buoyant particle scenarios.

1 Introduction

In [7], a Maxey–Riley [18] framework was presented to describe the dynamics of floating finite-size or inertial particles on the ocean surface, hereafter called BOM model, expanding on a model for submerged particles in the ocean [5]. The fluid mechanics Maxey–Riley equation is a Newton’s second law that accounts for various forces acting on spherical inertial particles carried by fluid flow [8]. The BOM theory has been validated both under natural field conditions [20, 16] and under controlled lab settings [14], and a previous model [6] successfully explained garbage patches in the subtropical gyres of the ocean. Moreover, the BOM model has been extended to simulate networks of elastically interacting floating inertial particles [4], offering a simplified model for Sargassum seaweed transport. Several results pertaining to the BOM model and its extensions are reviewed in [3], additional laboratory experiments utilizing the BOM model are discussed in [19], and an application of the “elastic” BOM model to Sargassum dynamics in the Caribbean Sea is detailed in [1].

The purpose of this short paper is to introduce an enhanced version of the BOM model that removes limitations on reserve buoyancy, $\sigma$, which denotes the portion of the particle’s volume that stays above the ocean surface (Sec. 2). This improvement is motivated by the observation made in [20], who pointed out that the water-to-particle density ratio cannot be too large—a limitation that we seek to relax. In doing this, the dependence on $\sigma$ of the projected length of the submerged particle piece, a factor in the drag force formula, is retrospectively determined using existing data on the windage parameter (Sec. 3). We complement the enhanced model with a reduced-order model that is asymptotically valid over the long term (Sec. 4). Finally, we consider the incorporation of general particle interactions and demonstrate how to account for the effects of Earth’s curvature (Sec. 5). Remarks have been interspersed in the text to clarify ideas and correct propagated typographical errors in earlier works.

2 Setup

Let $\mathbf{x}=(x,y)$ denote the position in a domain $D$ of the $\beta$-plane, where $x$ (eastward) and $y$ (northward) represent its coordinates. The domain $D\subset\mathbb{R}^{2}$ rotates with an angular velocity of $\frac{1}{2}f$, where $f=f_{0}+\beta y$ denotes the Coriolis parameter. Let $z$ be the vertical coordinate, and let $t$ represent time. Imagine two homogeneous fluid layers stacked on top of each other, separated by an interface fixed at $z=0$. The bottom layer, which represents water, has a density $\rho$, while the upper lighter layer represents air and has a density $\rho_{\mathrm{a}}$. Denote the dynamic viscosities of water and air as $\mu$ and $\mu_{\mathrm{a}}$, respectively. The velocities of water and air near the interface are represented as $\mathbf{v}(\mathbf{x},t)$ and $\mathbf{v}_{\text{a}}(\mathbf{x},t)$ accordingly. Finally, consider a solid spherical particle with radius $a$, presumed sufficiently small, and density $\rho_{\mathrm{p}}$ floating at the air–sea interface.

Definition 1 (Reserve volume).

Let

$$\delta:=\frac{\rho}{\rho_{\mathrm{p}}},\quad\delta_{\mathrm{a}}:=\frac{\rho_{\mathrm{a}}}{\rho_{\mathrm{p}}}\equiv\delta\frac{\rho_{\mathrm{a}}}{\rho}.$$ (1)

Let $0\leq\sigma\leq 1$ be the reserve buoyancy, that is, the fraction of emergent (above the ocean surface) particle volume. Static vertical stability of the particle (Archimedes’ principle),

$$(1-\sigma)\delta+\sigma\delta_{\mathrm{a}}=1,$$ (2)

is satisfied for

$$\sigma=\frac{\delta-1}{\delta-\delta_{\mathrm{a}}}\equiv\frac{1-\delta^{-1}}{1-\rho_{\mathrm{a}}/\rho},$$ (3)

which requires

$$1\leq\delta\leq\frac{\rho}{\rho_{\mathrm{a}}},\quad\frac{\rho_{\mathrm{a}}}{\rho}\leq\delta_{\mathrm{a}}\leq 1.$$ (4)

The configurations for neutrally buoyant and maximally buoyant particles reside at opposite extremes. The neutrally buoyant condition is characterized by $\sigma=0$, which translates to $\delta=1$ and $\delta_{\mathrm{a}}=0$. Conversely, the maximally buoyant case is characterized by $\sigma=1$, formally achieved by making $\delta=\rho/\rho_{\mathrm{a}}$ and $\delta_{\mathrm{a}}=1$.

Remark 1.

Two points should be highlighted:

1. Given that $\rho_{\mathrm{a}}/\rho$ is significantly smaller than 1, the reserve buoyancy $\sigma$ can be approximated by $1-\delta^{-1}$ (Fig. 1). Observe the distinct behaviors of $\sigma$ and $\delta$, the preferred parameter in the BOM theory. When $\delta$ equals 1, $\sigma$ vanishes, representing the neutrally buoyant configuration. In contrast, as $\delta$ increases, $\sigma$ approaches 1, corresponding to the maximally buoyant configuration. This cannot be described by the BOM theory as it imposes a constraint on $\delta$, preventing it from becoming too large [20].

2. The value of the parameter $\delta_{\mathrm{a}}$ does not have to be small in all cases. It is close to 0 when $\delta$ is finite (meaning $\sigma$ remains below 1), but it nears 1 as $\delta$ increases significantly (or in other words, as $\sigma$ approaches 1).

Figure 1: Reserve buoyancy as a function of the particle-to-water density ratio, which, unlike in the BOM theory, can be taken as large as desired.

The height ($h_{\mathrm{a}}$) of the emerged spherical cap can be described using $\sigma$. This results from matching its volume formula, expressed with $a$ and $h_{\mathrm{a}}$, to the volume of the emerged spherical cap. Specifically,

$$\frac{\pi h_{\mathrm{a}}^{2}}{3}(3a-h_{\mathrm{a}})=\sigma\frac{4}{3}\pi a^{3},$$ (5)

whose only physically meaningful root is

$$h_{\mathrm{a}}/a=\Phi:=\frac{\mathrm{i}\sqrt{3}}{2}\left(\frac{1}{\varphi}-\varphi\right)-\frac{1}{2\varphi}-\frac{\varphi}{2}+1$$ (6)

where

$$\varphi:=\sqrt[3]{\mathrm{i}\sqrt{1-(1-2\sigma)^{2}}+1-2\sigma}$$ (7)

(Fig. 2, left panel). The depth ($h$) of the submerged spherical cap,

$$h=(2-\Phi)a.$$ (8)

To be referenced in the future, the area of the emerged spherical cap projected in the direction of the flow, $A_{\mathrm{a}}$, is clearly given by

$$A_{\mathrm{a}}=\pi\Psi a^{2},\quad\Psi:=\pi^{-1}\cos^{-1}(1-\Phi)-\pi^{-1}(1-\Phi)\sqrt{1-(1-\Phi)^{2}},$$ (9)

a function of $\sigma$ exclusively (Fig. 2, right panel). In turn, the immersed projected area, denoted $A$, is equal to

$$A=\pi a^{2}-A_{\mathrm{a}}=\pi(1-\Psi)a^{2}.$$ (10)

Figure 2: As a function of reserve buoyancy, height of the emerged spherical cap normalized by radius (left) and area of the of the emerged spherical cap normalized by maximal sectional area (right).

3 Maxey–Riley equation for unrestricted reserve buoyancy

Fluid variables and parameters, denoted with a subscript f, differ for water and air, for example,

$$\mathbf{v}_{\text{f}}(\mathbf{x},z,t)=\begin{cases}\mathbf{v}_{\text{a}}(\mathbf{x},t)&\text{if }z\in(0,h_{\mathrm{a}}],\\ \mathbf{v}(\mathbf{x},t)&\text{if }z\in[-h,0).\end{cases}$$ (11)

Following [7], we write

$$\dot{\mathbf{v}}_{\text{p}}+f\mathbf{v}_{\text{p}}^{\perp}=\langle\mathbf{F}_{\mathrm{flow}}\rangle+\langle\mathbf{F}_{\mathrm{mass}}\rangle+\langle\mathbf{F}_{\mathrm{lift}}\rangle+\langle\mathbf{F}_{\mathrm{drag}}\rangle,$$ (12)

where $\perp$ denotes $+\frac{\pi}{2}$-rotation and $\langle\,\rangle$ is an average over $z\in[-h,h_{\mathrm{a}}]$. Here, $\mathbf{v}_{\text{p}}(t)$ represents the particle velocity, $\mathbf{F}_{\text{flow}}$ is the force per unit density of the undisturbed flow, $\mathbf{F}_{\text{mass}}$ is the added mass force, $\mathbf{F}_{\text{lift}}$ is the lift force from a horizontally sheared flow, and $\mathbf{F}_{\text{drag}}$ is the drag force from fluid viscosity. All forcing terms in Eq. (12) are part of the Maxey–Riley equation [18]. Exceptions are the Coriolis force term [5] and the lift force term [17]. For details, see [7, 3].

In contrast to [7], we calculate the vertical averages in Eq. (12) retaining all terms irrespective of their insignificance relative to the air-to-particle density ratio ($\delta_{\mathrm{a}}$) when considering a finite water-to-particle density ratio ($\delta$). We begin with the flow force:

$$\langle\mathbf{F}_{\mathrm{flow}}\rangle={} \frac{1}{2a}\int_{-h}^{h_{\mathrm{a}}}\frac{m_{\mathrm{f}}}{m_{\mathrm{p}}}\left(\frac{\operatorname{D}\!{\mathbf{v}_{\text{f}}}}{\operatorname{D}\!{t}}+f\mathbf{v}_{\text{f}}^{\perp}\right)\operatorname{d}\!{z}$$

$$={} \frac{1}{2a}\int_{(\Phi-2)a}^{0}\frac{\sigma\frac{4}{3}\pi a^{3}\rho}{\frac{4}{3}\pi a^{3}\rho_{\mathrm{p}}}\left(\frac{\operatorname{D}\!{\mathbf{v}}}{\operatorname{D}\!{t}}+f\mathbf{v}^{\perp}\right)\operatorname{d}\!{z}$$

$$+\frac{1}{2a}\int_{0}^{\Phi a}\frac{\sigma\frac{4}{3}\pi a^{3}\rho_{\mathrm{a}}}{\frac{4}{3}\pi a^{3}\rho_{\mathrm{p}}}\left(\frac{\operatorname{D}\!{\mathbf{v}_{\text{a}}}}{\operatorname{D}\!{t}}+f\mathbf{v}_{\text{a}}^{\perp}\right)\operatorname{d}\!{z}$$

$$={} \left(1-\frac{\Phi}{2}\right)(1-\sigma)\delta\left(\frac{\operatorname{D}\!{\mathbf{v}}}{\operatorname{D}\!{t}}+f\mathbf{v}^{\perp}\right)+\frac{\Phi}{2}\sigma\delta_{\mathrm{a}}\left(\frac{\operatorname{D}\!{\mathbf{v}_{\text{a}}}}{\operatorname{D}\!{t}}+f\mathbf{v}_{\text{a}}^{\perp}\right),$$ (13)

where

$$\frac{\operatorname{D}\!{\mathbf{v}_{\text{f}}}}{\operatorname{D}\!{t}}=\partial_{t}\mathbf{v}_{\text{f}}+v^{x}_{\mathrm{f}}\partial_{x}\mathbf{v}_{\text{f}}+v^{y}_{\mathrm{f}}\partial_{y}\mathbf{v}_{\text{f}}$$ (14)

is the material derivative of $\mathbf{v}_{\text{f}}$ along a fluid trajectory, i.e., obtained by solving $\dot{\mathbf{x}}=\mathbf{v}_{\text{f}}$.

Similarly, we compute

$$\langle\mathbf{F}_{\mathrm{mass}}\rangle =\frac{1}{2a}\int_{-h}^{h_{\mathrm{a}}}\frac{\frac{1}{2}m_{\mathrm{f}}}{m_{\mathrm{p}}}\left(\frac{\operatorname{D}\!{\mathbf{v}_{\text{f}}}}{\operatorname{D}\!{t}}+f\mathbf{v}_{\text{f}}^{\perp}-\dot{\mathbf{v}}_{\mathrm{p}}-f\mathbf{v}_{\text{p}}^{\perp}\right)\operatorname{d}\!{z}$$

$$=\frac{1}{2}\left(1-\frac{\Phi}{2}\right)(1-\sigma)\delta\left(\frac{\operatorname{D}\!{\mathbf{v}}}{\operatorname{D}\!{t}}+f\mathbf{v}^{\perp}-\dot{\mathbf{v}}_{\mathrm{p}}-f\mathbf{v}_{\mathrm{p}}^{\perp}\right)+\frac{1}{2}\frac{\Phi}{2}\sigma\delta_{\mathrm{a}}\left(\frac{\operatorname{D}\!{\mathbf{v}_{\text{a}}}}{\operatorname{D}\!{t}}+f\mathbf{v}_{\text{a}}^{\perp}-\dot{\mathbf{v}}_{\mathrm{p}}-f\mathbf{v}_{\text{p}}^{\perp}\right)$$ (15)

and, letting

$$\omega_{\mathrm{f}}=\partial_{x}v^{y}_{\mathrm{f}}-\partial_{y}v^{x}_{\mathrm{f}},$$ (16)

that is, the vorticity of the fluid, we likewise compute

$$\langle\mathbf{F}_{\mathrm{lift}}\rangle =\frac{1}{2a}\int_{-h}^{h_{\mathrm{a}}}\frac{\frac{1}{2}m_{\mathrm{f}}}{m_{\mathrm{p}}}\omega_{\mathrm{f}}\left(\mathbf{v}_{\text{f}}-\mathbf{v}_{\text{p}}\right)^{\perp}\operatorname{d}\!{z}$$

$$=\frac{1}{2}\left(1-\frac{\Phi}{2}\right)(1-\sigma)\delta\omega\left(\mathbf{v}-\mathbf{v}_{\text{p}}\right)^{\perp}+\frac{1}{2}\frac{\Phi}{2}\sigma\delta_{\mathrm{a}}\omega_{\mathrm{a}}\left(\mathbf{v}_{\text{a}}-\mathbf{v}_{\text{p}}\right)^{\perp}.$$ (17)

To determine the drag force, it is necessary to select suitable characteristic projected length scales for the submerged and emergent sections of the particle, denoted by $\ell$ and $\ell_{\mathrm{a}}$, respectively. We reasonably assume that $\ell$ and $\ell_{\mathrm{a}}$ both are smaller or equal to $2a$, with $\ell_{\mathrm{a}}=2a-\ell$ as a function of $\sigma$, constrained by $\ell(0)=2a$ and $\ell(1)=0$. Recognizing the heuristic nature of this approach, given the absence of an explicit formula for the drag force on spherical caps, the natural way forward is to leave $\ell(\sigma)$ for empirical determination from appropriate measurements. With this in mind, we compute:

$$\langle\mathbf{F}_{\mathrm{drag}}\rangle= {}\frac{1}{2a}\int_{-h}^{h_{\mathrm{a}}}\frac{12\mu_{\mathrm{f}}\frac{A_{\mathrm{f}}}{\ell_{\mathrm{f}}}}{m_{\mathrm{p}}}(\mathbf{v}_{\text{f}}-\mathbf{v}_{\text{p}})\operatorname{d}\!{z}$$

$$= {}\frac{1}{2a}\int_{(\Phi-2)a}^{0}\frac{12\mu\frac{\pi(1-\Psi)a^{2}}{\ell(\sigma)}}{\frac{4}{3}\pi a^{3}\rho_{\mathrm{p}}}(\mathbf{v}-\mathbf{v}_{\text{p}})\operatorname{d}\!{z}$$

$$+\frac{1}{2a}\int_{0}^{\Phi a}\frac{12\mu_{\mathrm{a}}\frac{\pi\Psi a^{2}}{2a-\ell(\sigma)}}{\frac{4}{3}\pi a^{3}\rho_{\mathrm{p}}}(\mathbf{v}_{\text{a}}-\mathbf{v}_{\text{p}})\operatorname{d}\!{z}$$

$$= {}\frac{3}{2}\left(\frac{2}{3}+\frac{(1-\sigma)\delta}{3}+\Big{(}\sigma\delta_{\mathrm{a}}-\sigma\delta\Big{)}\frac{\Phi}{6}\right)\frac{\mathbf{u}-\mathbf{v}_{\text{p}}}{\tau},$$ (18)

where

$$\mathbf{u}:=(1-\alpha)\mathbf{v}+\alpha\mathbf{v}_{\text{a}},$$ (19)

and the parameters

$$\tau :=\frac{\frac{2}{3}+\frac{(1-\sigma)\delta}{3}+\Big{(}\sigma\delta_{\mathrm{a}}-(1-\sigma)\delta\Big{)}\frac{\Phi}{6}}{3\left((2-\Phi)(1-\Psi)+\frac{\ell(\sigma)}{2a-\ell(\sigma)}\gamma\Phi\Psi\right)\delta}\cdot\frac{a\ell(\sigma)}{\mu/\rho},$$ (20)

$$\alpha :=\frac{\gamma\Phi\Psi}{(2-\Phi)(1-\Psi)\frac{2a-\ell(\sigma)}{\ell(\sigma)}+\gamma\Phi\Psi},$$ (21)

$$\gamma :=\frac{\mu_{\mathrm{a}}}{\mu}.$$ (22)

Finally, plugging Eqs. (13)–(18) into Eq. (12) and after some algebraic manipulation, we obtain the following Maxey–Riley equation:

$$\dot{\mathbf{v}}_{\mathrm{p}}+\left(f+\tfrac{1}{3}R\omega+\tfrac{1}{3}R_{\mathrm{a}}\omega_{\mathrm{a}}\right)\mathbf{v}_{\text{p}}^{\perp}+\tau^{-1}\mathbf{v}_{\text{p}}=R\frac{\operatorname{D}\!{\mathbf{v}}}{\operatorname{D}\!{t}}+R\left(f+\tfrac{1}{3}\omega\right)\mathbf{v}^{\perp}+R_{\mathrm{a}}\frac{\operatorname{D}\!{\mathbf{v}_{\text{a}}}}{\operatorname{D}\!{t}}+R_{\mathrm{a}}\left(f+\tfrac{1}{3}\omega_{\mathrm{a}}\right)\mathbf{v}_{\text{a}}^{\perp}+\tau^{-1}\mathbf{u},$$ (23)

where

$$R :=\frac{\left(1-\frac{\Phi}{2}\right)(1-\sigma)\delta}{\frac{2}{3}+\frac{(1-\sigma)\delta}{3}+\Big{(}\sigma\delta_{\mathrm{a}}-(1-\sigma)\delta\Big{)}\frac{\Phi}{6}},$$ (24)

$$R_{\mathrm{a}} :=\frac{\frac{\Phi}{2}\sigma\delta_{\mathrm{a}}}{\frac{2}{3}+\frac{(1-\sigma)\delta}{3}+\Big{(}\sigma\delta_{\mathrm{a}}-(1-\sigma)\delta\Big{)}\frac{\Phi}{6}}.$$ (25)

Typically, the parameter $\gamma$ is approximately 0.0167. The parameter $\tau$ is positive and small, because $a$ is assumed to be small. The parameter $\alpha$ increases from 0 to 1 as $\sigma$ increases from 0 to 1, which can be expected regardless of the specific form of the projected length scale of the submerged spherical cup $\ell(\sigma)$. We only need to assume that $\ell$ decreases from $2a$ to 0 as $\sigma$ increases from 0 to 1. Finally, the parameter $R$ decreases from 1 to 0 as $\sigma$ increases from 0 to 1, whereas $R_{\mathrm{a}}$ increases from 0 to 1.

Remark 2.

The parameter $\tau$ can be identified with Stokes’ time, which represents the time required for a particle to respond to inertial effects. The parameter $\alpha$ can be called windage parameter because it quantifies the contribution of air velocity to the carrying flow velocity $\mathbf{u}$, given by the convex combination of water and air velocities in Eq. (19).

Interestingly, $\alpha$ measures windage similarly to the “leeway parameter” commonly used in the search-and-rescue literature [2]. However, unlike $\alpha$, the leeway parameter is based primarily on ad hoc considerations rather than derived from first principles.

In order to complete the formulation of the Maxey–Riley model, as specified in Eq. (23), it is required to evaluate $\ell(a)$. An evaluation of $\ell(a)$ can be achieved by using results from controlled laboratory measurements of the windage parameter $\alpha$ within the range $0\leq\delta\lessapprox 4$, as documented by [15]. This $\delta$-range corresponds to $0\leq\sigma\lessapprox 0.75$, as it follows from Eq. (3) as a consequence of the smallness of $\rho_{\mathrm{a}}/\rho$ in the noted finite $\delta$-range. These empirical observations suggest that, within the inferred $\sigma$-range,

$$\alpha\approx\frac{\gamma\Psi}{1+(\gamma-1)\Psi}.$$ (26)

Consequently,

$$\ell\approx(2-\Phi)a=h,$$ (27)

indicating that, over the range $0\leq\sigma\lessapprox 0.75$, the appropriate characteristic projected length scale of the submerged spherical cap is simply its depth. While this might appear self-evident, the drag force formula applies to complete spheres rather than spherical caps. Consistent with Eq. (27),

$$\tau\approx\frac{\frac{2}{3}+\frac{(\sigma-1)\delta}{3}+\Big{(}\sigma\delta_{\mathrm{a}}-(1-\sigma)\delta\Big{)}\frac{\Phi}{6}}{3\left(1+(\gamma-1)\Psi\right)\delta}\cdot\frac{a^{2}}{\mu/\rho}\approx\frac{1-\frac{1}{6}\Phi}{3(1+(\gamma-1)\Psi\delta}\cdot\frac{a^{2}}{\mu/\rho},$$ (28)

where the last approximation holds because $\sigma\approx 1-\delta^{-1}$ independent of the value of $\delta$. The derived formula for $\tau$ differs by a factor of $\delta^{-3}$ from the one obtained through empirical fitting to observational data reported in [20]. However, these observations dealt with objects of various geometries, not just spheres, which required applying the heuristic shape correction introduced in [10]. Extending the closure in Eq. (27) beyond the interval $0\leq\sigma\lessapprox 0.75$ will require further measurements. It is important to note that the interval $0\leq\sigma\lessapprox 0.75$ encompasses a substantial portion of the possible range of reserve buoyancy values. Therefore, it appears that Eq. (27) already offers an empirically determined projected length well-supported by data.

Remark 3.

In the Appendix of [15], the expression in Eq. (26) contains a typo where $1-\gamma$ is used instead of $\gamma-1$. This typographical error has been propagated to [16, 4], and we wish to bring attention to this issue.

Lastly, all terms related to air in the Maxey–Riley model, Eq. (23), can be safely ignored when $\delta$ is finite (i.e., $\sigma$ is sufficiently smaller than 1), in which case $R_{\mathrm{a}}$ in Eq. (25) is small due to the smallness of $\delta_{\mathrm{a}}$. These air-related terms are excluded from the BOM equation. However, these terms will gain significance as $\delta$ becomes sufficiently large (or equivalently as $\sigma$ approaches 1), thereby expanding the applicability of the BOM model.

4 Relevant limits and slow manifold reduction

For neutrally buoyant particles, i.e., when $\sigma=0$, which corresponds to $\delta=1$ and $\delta_{\mathrm{a}}=\rho_{\mathrm{a}}/\rho$, Eq. (23) reduces to

$$\dot{\mathbf{v}}_{\mathrm{p}}+\left(f+\tfrac{1}{3}\omega\right)\mathbf{v}_{\text{p}}^{\perp}+\tau^{-1}\mathbf{v}_{\text{p}}=\frac{\operatorname{D}\!{\mathbf{v}}}{\operatorname{D}\!{t}}+\left(f+\tfrac{1}{3}\omega\right)\mathbf{v}^{\perp}+\tau^{-1}\mathbf{v},$$ (29)

with

$$\tau=\frac{a\ell(0)}{6\mu/\rho}=\frac{a^{2}}{3\mu/\rho}.$$ (30)

Conversely, in the case of maximally buoyant particles, i.e., for $\sigma=1$, which translates to $\delta=\rho/\rho_{\mathrm{a}}$ and $\delta_{\mathrm{a}}=1$, Eq. (23) simplifies to

$$\dot{\mathbf{v}}_{\mathrm{p}}+\left(f+\tfrac{1}{3}\omega_{\mathrm{a}}\right)\mathbf{v}_{\text{p}}^{\perp}+\tau^{-1}\mathbf{v}_{\text{p}}=\frac{\operatorname{D}\!{\mathbf{v}_{\text{a}}}}{\operatorname{D}\!{t}}+\left(f+\tfrac{1}{3}\omega_{\mathrm{a}}\right)\mathbf{v}_{\text{a}}^{\perp}+\tau^{-1}\mathbf{v}_{\text{a}},$$ (31)

with

$$\tau=\frac{a\ell_{\mathrm{a}}(1)}{6\mu_{\mathrm{a}}/\rho_{\mathrm{a}}}=\frac{a^{2}}{3\mu_{\mathrm{a}}/\rho_{\mathrm{a}}}.$$ (32)

The above equations exactly match the Maxey–Riley equation in fluid mechanics [18] when considering particles that have the same density as the carrying fluid (with the assumptions of motion on a $\beta$-plane, and the inclusion of lift force).

The movement of neutrally buoyant particles generally differs from the movement of water particles. Under specific conditions, it can be demonstrated to be unstable [12], irrespective of the Coriolis and lift effects [7, 3]. This applies mutatis mutandis to the motion of maximally buoyant particles with respect to that of air particles. These results are relevant in the short run as, over time, the movement of both neutrally buoyant and maximally buoyant particles aligns with the Lagrangian movement, meaning they align with the movement of water and air particles, respectively.

Remark 4.

The motion of fully exposed particles is significantly more restricted by the assumption of a flat air–sea interface compared to the movement of partially or completely submerged particles. For these submerged particles, the influence of a wavy interface can be considered by incorporating a representation of the Stokes drift in the near-surface water velocity. This does not have an equivalent for the near-surface air velocity. Thus, the maximally buoyant limit should be considered primarily of theoretical interest rather than of practical value.

To understand how aligment with Lagrangian motion is realized in the long run, we note that for sufficiently small inertial response time $\tau$, Eq. (23) reduces time-asymptotically to

$$\dot{\mathbf{x}}\sim\mathbf{u}+\tau\mathbf{u}_{\tau},$$ (33)

where

$$\mathbf{u}_{\tau}:=R\frac{\operatorname{D}\!{\mathbf{v}}}{\operatorname{D}\!{t}}+R_{\mathrm{a}}\frac{\operatorname{D}\!{\mathbf{v}_{\text{a}}}}{\operatorname{D}\!{t}}+R\left(f+\tfrac{1}{3}\omega\right)\mathbf{v}^{\perp}+R_{\mathrm{a}}\left(f+\tfrac{1}{3}\omega_{\mathrm{a}}\right)\mathbf{v}_{\text{a}}^{\perp}-\frac{\operatorname{D}\!{\mathbf{u}}}{\operatorname{D}\!{t}}-\left(f+\tfrac{1}{3}R\omega+\tfrac{1}{3}R_{\mathrm{a}}\omega_{\mathrm{a}}\right)\mathbf{u}^{\perp}.$$ (34)

This equation governs the motion on the slow manifold, namely, the (2+1)-dimensional manifold $M_{\tau}:=\{(\mathbf{x},\mathbf{v}_{\text{p}},t):\mathbf{v}_{\text{p}}=\mathbf{u}(\mathbf{x},t)+\tau\mathbf{u}_{\tau}(\mathbf{x},t)\}$ embedded in the (4+1)-dimensional extended phase space $(\mathbf{x},\mathbf{v}_{\text{p}},t)$. Rigorously speaking, $M_{\tau}$ represents a normally hyperbolic manifold that has the ability to exponentially attract all solutions of the complete Maxey–Riley equation, Eq. (23). A technical observation is that such an attraction cannot be expected to be monotonic when the carrying flow velocity $\mathbf{u}$ varies rapidly. To derive the reduce equation (33), one must recognize that the complete equation, Eq. (23), comprises slow ($\mathbf{x}$) and fast ($\mathbf{v}_{\text{p}}$) variables, enabling the application of geometric singular perturbation theory [9, 13], extended to non-autonomous systems [11]. This approach was utilized in [7, 3] for the BOM equation and can be straightforwardly applied on Eq. (23) to obtain Eq. (33).

Note that when particles have neutral buoyancy, $\mathbf{u}_{\tau}\equiv 0$ and therefore Eq. (33) reduces to $\dot{\mathbf{x}}\sim\mathbf{v}$. In a similar fashion, for particles exhibiting maximum buoyancy, $\mathbf{u}_{\tau}\equiv 0$, reducing Eq. (33) to $\dot{\mathbf{x}}\sim\mathbf{v}_{\text{a}}$. This elucidates the earlier mention of synchronization with Lagrangian motion over time.

5 Final considerations

We note that a Maxey–Riley theory can be developed using Eq. (23) for systems of inertial particles interacting through a force, such as elastic coupling, which has been proposed for modeling Sargassum seaweed transport [4]. Various other types of interactions are also possible. For example, they may represent cohesion required to simulate the movement of oil or sediment particles. Regardless of the interaction type, the $i$th particle in a system of $N$ particles will evolve as described by Eq. (23), and in the long term by Eq. (33), with the coupling force

$$\mathbf{F}_{i}=\sum_{j\in\operatorname{neighbor}(i)}\mathbf{F}(|\mathbf{x}_{ij}|),$$ (35)

where $|\mathbf{x}_{ij}|$ is the distance between particles $i$ and $j$. The term $\mathbf{F}_{i}$ enters as an additional term in the right-hand-side of Eq. (23) and of Eq. (33), but multiplied by $\tau$, under the implicit assumption that $\mathbf{F}_{i}=O(1)$ as measured by $\tau$.

Another remark involves the derivation of both the complete, Eq. (23), and simplified, Eq. (33), Maxey–Riley equations on a spherical surface. Consider the rescaled longitude–latitude coordinates $x=a_{\odot}\cos\vartheta_{0}(\lambda-\lambda_{0})$ and $y=a_{\odot}(\vartheta-\vartheta_{0})$, where $a_{\odot}$ is the Earth’s average radius, and $(\lambda_{0},\vartheta_{0})$ is a reference point. We define $\gamma_{\odot}(y):=\sec\vartheta_{0}\cos\vartheta$ and $\tau_{\odot}(y):=a_{\odot}^{-1}\tan\vartheta$, which account for the geometric effects of Earth’s spherical shape [21]. The fluid vorticity

$$\omega_{\mathrm{f}}=\gamma_{\odot}^{-1}\partial_{x}v^{y}_{\mathrm{f}}-\partial_{y}v^{x}_{\mathrm{f}}+\tau_{\odot}v^{x}_{\mathrm{f}}.$$ (36)

The covariant derivative of the fluid velocity $\mathbf{v}_{\text{f}}$ along a fluid trajectory,

$$\frac{\operatorname{D}\!{\mathbf{v}_{\text{f}}}}{\operatorname{d}\!{t}}=\partial_{t}\mathbf{v}_{\text{f}}+\gamma_{\odot}^{-1}v^{x}_{\mathrm{f}}\partial_{x}\mathbf{v}_{\text{f}}+v^{y}_{\mathrm{f}}\partial_{y}\mathbf{v}_{\text{f}}+\tau_{\odot}v^{x}_{\mathrm{f}}\mathbf{v}_{\text{f}}^{\perp}.$$ (37)

To obtain the full Maxey–Riley equation on the sphere, one needs to:

1. replace the vorticities in Cartesian coordinates in Eq. (23), given by Eq. (16), with those in spherical coordinates specified in Eq. (36), and substitute the material derivatives in Eq. (14) with covariant derivatives from Eq. (37); and

2. add the term $\tau_{\odot}v^{x}_{\mathrm{p}}\mathbf{v}_{\text{p}}^{\perp}$ to the left-hand side of Eq. (23).

For deriving the reduced Maxey–Riley equation on the sphere, one should:

1. similarly replace the material derivatives in Eq. (33) with covariant derivatives from Eq. (37), extending $\mathbf{v}_{\text{f}}$ to also represent $\mathbf{u}$;

2. subtract from the right-hand side of Eq. (33) the term $\tau_{\odot}u^{x}\mathbf{u}^{\perp}$; and

3. multiply the left-hand side of Eq. (33) by $\smash{\sqrt{\mathsf{m}}}$, where the matrix $\mathsf{m}$ provides a representation of the spherical metric in the rescaled coordinates $(x,y)$ introduced above, with elements $\mathsf{m}_{11}=\gamma_{\odot}^{2}$, $\mathsf{m}_{11}=1$, and $\mathsf{m}_{12}=0=\mathsf{m}_{21}$.

Remark 5.

We have identified three inaccuracies in [7] that we take the chance to rectify: 1) the initial equality in Eq. (A5) should be omitted; 2) the term $\dot{x}$ in Eq. (A8) must be scaled by $\sqrt{\mathsf{m}}$; and 3) the term $\tau_{\odot}v^{1}$ preceding $\frac{1}{3}R\omega$ in Eq. (A8) should be replaced with $2\tau_{\odot}u^{1}$ (the discrepancy in the total derivative definition with respect that used in this paper results in the factor of 2). The first error was propagated to [3], as can be seen in Eq. (39) in that paper. The third error was also propagated to [16]; however, in that study, the complete BOM equation was employed, as opposed to its reduced form applicable asymptotically on the slow manifold, thereby not affecting the results.

Author declarations

Conflict of interest

The author has no conflicts to disclose.

Author contributions

This paper is authored by a single individual who entirely carried out the work.

Data availability

Data sharing is not applicable to this article as no new data were created or analyzed in this study.