Ice Spiral Patterns on the Ocean Surface

TL;DR

This paper introduces a novel hydrodynamic model to explain large-scale ice swirl patterns on the ocean surface, analyzing their formation, stability, and geometric properties through multiple mathematical approaches.

Contribution

The study develops a new two-dimensional compressible Navier-Stokes model for ice swirls, deriving analytical solutions and stability analysis that align with experimental observations.

Findings

Logarithmic spiral solutions are derived from the nonlinear equations.

Multiple spiral modes, including pure radial and azimuthal, are obtained.

The stability analysis reveals geometric phase transitions in the patterns.

Abstract

We investigate a new two-dimensional compressible Navier-Stokes hydrodynamic model design to explain and study large scale ice swirls formation at the surface of the ocean. The linearized model generates a basis of Bessel solutions from where various types of spiral patterns can be generated and their evolution and stability in time analyzed. By restricting the nonlinear system of equations to its quadratic terms we obtain swirl solutions emphasizing logarithmic spiral geometry. The resulting solutions are analyzed and validated using three mathematical approaches: one predicting the formation of patterns as Townes solitary modes, another approach mapping the nonlinear system into a sine-Gordon equation, and a third approach uses a series expansion. Pure radial, azimuthal and spiral modes are obtained from the fully nonlinear equations. Combinations of multiple-spiral solutions are also obtained, matching the experimental observations. The nonlinear stability of the spiral patterns is analyzed by Arnold's convexity method, and the Hamiltonian of the solutions is plotted versus some order parameters showing the existence of geometric phase transitions.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \vec{v})=0,$$

$$\frac{\partial (\rho \vec{v})}{\partial t}+\nabla \cdot (\rho \vec{v}\otimes \vec{v})=-\nabla P+\nabla \cdot \mathcal{T}+\vec{f},$$

$$\mathcal{T}^{ki}=\eta\biggl{(}\frac{\partial v^{i}}{\partial x_{k}}+\frac{\partial v^{k}}{\partial x_{i}}\biggr{)},$$

$$\Phi (x,y,t)=\frac{A_i}{A},0\le \Phi \le 1,$$

$$\rho (x,y,t)=\Phi {\rho }_i+(1-\Phi ){\rho }_w.$$

$$P(x,y,t)=\Phi P_i+(1-\Phi )P_w,$$

$$P_i={\rho }_{w}gT.$$

$$(1-a\Phi)\biggl{(}u_{r}+\frac{u}{r}+\frac{v_{\theta}}{r}\biggr{)}-a\biggl{(}\Phi_{t}+u\Phi_{r}+\frac{v\Phi_{\theta}}{r}\biggr{)}=0,$$

$$(1-a\Phi)\biggl{(}u_{t}+uu_{r}+\frac{vu_{\theta}}{r}-\frac{v^{2}}{r}\biggr{)}-a\biggl{(}u\Phi_{t}+u^{2}\Phi_{r}+\frac{uv\Phi_{\theta}}{r}\biggr{)}=\biggl{(}\frac{P_{w}}{\rho_{w}}-gT_{0}\biggr{)}\Phi_{r},$$

$$(1-a\Phi)\biggl{(}v_{t}+\frac{uv}{r}+uv_{r}+\frac{vv_{\theta}}{r}\biggr{)}-a\biggl{(}v\Phi_{t}+uv\Phi_{r}+\frac{v^{2}\Phi_{\theta}}{r}\biggr{)}=\biggl{(}\frac{P_{w}}{\rho_{w}}-gT_{0}\biggr{)}\frac{\Phi_{\theta}}{r},$$

$${\Psi }_t+(\Psi u{)}_r+\frac{(\Psi v{)}_{\theta }}{r}+\frac{\Psi u}{r}=0,$$

$$(\Psi u{)}_t+u(\Psi u{)}_r+\frac{v(\Psi u{)}_{\theta }}{r}-\frac{\Psi v^2}{r}+\frac{b}{a}{\Psi }_r=0,$$

$$(\Psi v{)}_t+u(\Psi v{)}_r+\frac{v(\Psi v{)}_{\theta }}{r}+\frac{\Psi uv}{r}+\frac{b}{ar}{\Psi }_{\theta }=0.$$

$$(\Psi u{)}_t+\frac{b}{a}{\Psi }_r=0,(\Psi v{)}_t+\frac{b}{ar}{\Psi }_{\theta }=0,$$

$$\frac{u}{v}=-\frac{i}{m}\frac{r{\Psi }_r}{\Psi }.$$

$$\frac{u}{v}=\frac{R_{\theta }}{R}.$$

$$\Psi =e^{i(m{g}_0\ln r+m\theta +\omega t)},$$

$$\Phi ={\Phi }^0+\sum _{n=1}^{\infty }{\epsilon }^n{\Phi }^n(r,\theta ,t),$$

$$u=\sum _{n=1}^{\infty }{\delta }^n{u}^n(r,\theta ,t),v=\sum _{n=1}^{\infty }{\delta }^n{v}^n(r,\theta ,t),$$

$$v_{gr}\sim \mathcal{O}(u)\sim \mathcal{O}(v)\sim \frac{\epsilon R}{\delta \Theta }\sim \frac{\epsilon \Theta v_{ph}^2}{\delta R}.$$

$$-a{\Phi }_t-au{\Phi }_r+(1-a\Phi )u_r+(1-a\Phi )\frac{u}{r}+(1-a\Phi )\frac{v_{\theta }}{r}-a\frac{v{\Phi }_{\theta }}{r}=0,$$

$$(1-a\Phi )u_t-au{\Phi }_t+(1-a\Phi )u{u}_r-a{u}^2{\Phi }_r+(1-a\Phi )\frac{v{u}_{\theta }}{r}$$

$$-a\frac{uv\Phi_{\theta}}{r}-(1-a\Phi)\frac{v^{2}}{r}=\biggl{(}\frac{P_{w}}{\rho_{w}}-gT_{0}\biggr{)}\Phi_{r},$$

$$(1-a\Phi )v_t-av{\Phi }_t+(1-a\Phi )u{v}_r-(1-a\Phi )\frac{uv}{r}-a\rho uv{\Phi }_r$$

$$+(1-a\Phi)\frac{vv_{\theta}}{r}-a\rho\frac{v^{2}\Phi_{\theta}}{r}=\biggl{(}\frac{P_{w}}{\rho_{w}}-gT_{0}\biggr{)}\frac{\Phi_{\theta}}{r}.$$

$$\frac{\varepsilon R}{\delta\Theta}\Phi_{\tau}^{1}-\biggl{(}\frac{1}{a}-\Phi_{0}\biggr{)}\biggl{(}\frac{u^{1}}{\zeta}+\frac{v_{\theta}^{1}}{\zeta}+u^{1}_{\zeta}\biggr{)}=0,$$

$$(1-a\Phi^{0})u^{1}_{\tau}=\frac{\varepsilon\Theta}{\delta R}\biggl{(}\frac{P_{w}}{\rho_{w}}-gT_{0}\biggr{)}\Phi^{1}_{\zeta},$$

$$(1-a\Phi^{0})v^{1}_{\tau}=\frac{\varepsilon\Theta}{\delta R}\biggl{(}\frac{P_{w}}{\rho_{w}}-gT_{0}\biggr{)}\frac{\Phi^{1}_{\theta}}{\zeta},$$

$$\frac{\partial }{\partial t}[u_{\theta }^1-(r{v}^1{)}_r]=-r{w}_t^1=0,$$

$$\Phi^{1}_{tt}-c^{2}\biggl{(}\Phi^{1}_{rr}+\frac{\Phi^{1}_{r}}{r}+\frac{\Phi^{1}_{\theta\theta}}{r^{2}}\biggr{)}=0,$$

Full text

††footnotetext: Emails: zongzhi@dlut.edu.cn, ludua@erau.edu

1 Introduction

The rapid decline of summer ice extent that has occurred in the Arctic Ocean over recent years has prompted a surge of research activity and, in particular, the role of sea-ice morphology, has been increasingly recognized. This is especially relevant in the context of climate change as the resulting ice melting and proliferation of open water promote further wave growth over increasing fetches, thus allowing long waves to propagate larger distances into the ice field. Of particular interest is the marginal ice zone which is the fragmented part of the ice cover closest to the open ocean and, as such, it is a very dynamic region strongly affected by ocean currents. The formation of sea-ice fragments increases the further penetration and damage of the ice cover, so it is highly important to understand the dynamics of these fragments on the ocean surface. In this respect, a lot of research was dedicated to linear and nonlinear ocean waves propagating through fragmented sea-ice [4, 1, 2, 3, 5, 6].

Sea-ice can directly collide with marine structures, [5], and induce wave scattering or attenuation, [7, 8], thus indirectly varying the hydrodynamic response of passing vessels and offshore platform. To safely manage the navigation and functioning of such structures, and to develop accurate models for their interaction with sea-ice, it is important to be able to predict the distribution and morphology of sea-ice fragments.

The mechanism for the formation of swirl patterns observed mostly in arctic ocean, see Figs. 1,5, [9, 10], has not yet been fully understood or modeled. There is little doubt, from the observational data available, that these huge scale phenomenon is associated with the ocean currents and possibly wind and waves [11]. Nevertheless, given the very slow time scale, it is unlikely for these spiral to be generated or enhanced by Coriolis force. A complete model should consider local viscosity generated by collisions of ice fragments, and subsequently the local phase transitions, clustering effect, the vertical motion of water and the elevation waves effect, the wind, and a certain probability distribution for the ice fragments size and shapes. Such a complete theory, taking into consideration all these components and forces into account, and put their relative importance into perspective is needed, but such a theory is not yet available. There is a considerable body of research on the interaction between waves, [3, 4, 5, 6], solitons and sea-ice, [1, 2], and quasi-granular aggregate models for the sea-ice [12]. The most natural explanation regarding such sea-ice swirl structure is the formation of a wave pattern, which either remains stationary, or at least quasi-stationary, in a frame of reference rotating around its center at a proper angular speed. Similar waves of patterns structures were used to explain the spiral galaxies formation [13]. In this manuscript we favor the point of view that the sea-ice can maintain a density wave through water currents interaction. This density wave provides a spiral density field which underlies the observable concentration of sea-ice. In this way, an observable swirl pattern can be maintained over the whole structure. In this paper we demonstrate the formation and stability of such ice swirls. We introduce a two-dimensional model describing large space-time scale spiral patterns of sea-ice fragments floating on water, like the ones observed in arctic ocean, [9, 10], see for example Fig. 1.

The paper is organized as follows: in section 2 we introduce a 2-dimensional two-phases compressible fluid model governed on mass and momentum conservation of water plus ice. In section 3 we expand the density and velocity fields in a two-scale series, controlled by two smallness parameters and the system of equations is linearized and solved exactly under initial and boundary conditions and the sea-ice swirls solutions are obtained and discussed. In section 4 we analyze the time evolution of various types of linearized solutions and spirals and their stability. In section 5 we introduce limiting situations for the ice patterns, namely the pure radial motion, the azimuthal (rotational) motion, and the spiral patterns. In section 6, following a qualitative discussion on the structure of the nonlinear system and its solutions, we find solutions by using three approaches: mapping the nonlinear system into a sine-Gordon equation, using an iterative procedure of partial differential operators, and using a quadratic truncation. In section 7 we study the nonlinear stability of the spiral solutions.

2 Inviscid two-fluid model

In our two-dimensional hydrodynamic model the water (density $\rho_{w}$) is assumed incompressible and inviscid, and elevation waves at the water surface are neglected. Multiple studies of both linear, [14], and nonlinear, [1, 2, 3], surface waves propagating through fragmented sea-ice show that the attenuation by scattering, damping, multiple wave reflections and ice viscosity is most effective for floe configuration representing a good compromise between ice concentration and ice fragmentation, combination which represents the typical distribution in the ice swirl patterns. Given the relative short attenuation length of such ocean waves through sea-ice fragments, [2, 3, 5], and the long life time of the observed ice swirls, it seems natural to neglect the influence of ocean waves in modeling the dynamics of ice swirl formation, while ocean currents and winds represent to main contributions.

We consider that the ice fragments (density $\rho_{i}$) have average size much smaller than the characteristic size of the pattern, so we can assume an almost uniform size, shape, and volume distribution among fragments, and almost constant ice draft (denoted $T_{0}$), except special local situation at the center of the swirl. We model the mixture between water and floating sea-ice fragments as a compressible two-dimensional fluid in the surfactant approximation. The flow takes places in the $(x,y)$-plane with $z$-axis oriented upwards, and $z=0$ being the water surface. The two-fluid system is considered in isothermal and isobaric equilibrium, without substantial amounts of phase transitions between ice and water. Other higher order of approximation effects, like vertical displacement of ice fragments induced by ocean waves, by wind or by adjacent ice fragments collision, ice elasticity, or strong collisions of ice fragments followed by local ice melting will be considered in a subsequent model.

For the extended ice swirl formations the average value of the water flow is relatively small $|\vec{v}|\sim 1\div 10$ mm/s, [15], while the sea-ice fragments horizontal size can range between $10\div 50$ m [12]. In this conditions the Reynolds number ranges Re$\sim 10^{4}\div 10^{7}$, the Froude number ranges Fr$\sim 10^{-6}\div 10^{-4}$ and Euler number is of the order of the unit. Consequently, the system is governed by the law of mass conservation for the mixture ice-water and by the compressible Cauchy linear momentum equation [16]

[TABLE]

[TABLE]

where $\vec{v}$ is the flow field, $\rho$ is the density, $P$ is the total pressure, the symbol $\otimes$ is the direct tensor product, $\vec{f}$ is the volume density of external forces, and $\mathcal{T}$ is the second order symmetric viscous stress tensor

[TABLE]

and $\eta$ is the dynamical viscosity. The high Reynolds number value and the low Froude number value for the ice swirl patterns suggest a sub-critical regime in which Newtonian viscosity can be neglected in the first approximation. Indeed, from the Navier-Stokes Eq. (2) written in dimensionless form is known that the Poiseuille number Poi$=$Eu$\cdot$Re$\sim PL/(v\eta)$ controls the importance of the viscosity term, where $L$ is the typical length scale of the system. The pressure of water upon the ice fragments can be evaluated by considering water uniformly pushing with a Stokes drag force upon the submerged regions of the ice fragment. For a regular prismatic ice fragment, for example, it results Poi $\sim 0.05$ which means that it is sufficient to consider viscosity in the second order of approximation, for the linearized equations. If for some reason the ice fragments tend to cluster (for example at large impact velocity) the friction between the mutual fragments can trigger local melting followed by re-solidification, which generate a local increase in the effective viscosity of the mixture. In such a case the floating ice mixture can be considered at large scale as a non-Newtonian fluid, specifically a share thickening (dilatant) mixture where viscosity increases with the rate of shear strain. In other words, the ice fragments and water can transition to coagulation, forming larger ice blocks. In this case the water-ice mixture can be modeled as a Bingham type of fluid with a viscous shear-thinning power flow [17]. Such situation occurs towards the central bulge of the spiral. In the following, we elaborate on a basic model for inviscid flow.

For the two-dimensional mixture of water sea-ice fragments we introduce the Ice Fraction Function (IFF) as the surface density of ice vs. surface, namely the ratio of ice in a unit surface

[TABLE]

where $A_{i}$ is the area occupied by brush ice in the unit area, and $A$ is the unit area under consideration. Correspondingly, the area occupied by free surface water is $A_{w}=A(1-\Phi)$.

Similarly to the introduction of fractional volume scalar field in air-water-ice-structure CFD modeling, [5], we introduce the density of the water-ice mixture by

[TABLE]

The movement of sea-ice and ocean currents (wind stress is neglected here) are combined in the pressure term in our model. Since we neglect the deformation of sea surface, the oceanic currents act on ice by ocean–ice interfacial stress which is related to the velocity difference between the ice movement and surface ocean currents. The momentum transfer is usually from the ice to the ocean, [18], case in which the pressure is given by the ice pressure on water through its buoyancy. Such pressure distribution acting on the sides of the ice fragments is considered positive because acts towards water. However, the pressure can become slightly negative if the momentum transfer is from water towards the ice [19]. Usually in sea-ice modeling, [20], the pressure is taken as the weighted sum of pressures of the two phases

[TABLE]

where $P_{i},P_{w}$ are the ice and water pressure, respectively. Neglecting the dynamic effects because the motion is very slow, the pressure on an ice fragment is approximately equal to its static pressure. We can express the ice pressure per unit of length of ice horizontal perimeter by

[TABLE]

As we mentioned above, we assume this draft to be in average independent of the ice fragment, so we treat it as a constant $T=T_{0}$. In the following we use the notation $\vec{v}=u\vec{e}_{r}+v\vec{e}_{\theta}$ for the velocity field, and $\hat{i},\hat{j},\vec{e}_{r}$ and $\vec{e}_{\theta}$ for the unit vectors, correspondingly. The derivatives with respect to $t,x,y,r,\theta$ are labeled with the corresponding letter subscript, while all other subscripts in this text (like for example $i,w$ for ice and water) do not represent derivatives.

In the following, by introducing Eq. (5) in Eq. (1) with $\eta=0$ and expressed in polar coordinates $(r,\theta)$, we obtain a nonlinear differential equation for the mass conservation, in terms of IFF and velocity components

[TABLE]

where $a=1-\rho_{i}/\rho_{w}$ is a constant. By implementing Eqs. (5-7) in Eq. (2) without the viscosity term, we obtain an equation for each polar component of the law of momentum conservation

[TABLE]

for radial component, and

[TABLE]

for the azimuthal one. This system can be written in a more compact form using the notation $\Psi=(1-a\Phi)$ and $b=P_{w}/\rho_{w}-gT_{0}$

[TABLE]

[TABLE]

[TABLE]

The system Eqs. (8-10) , or Eqs. (11-13), describes the dynamics of the three fields $\Phi$ (or $\Psi$), $u$ and $v$ depending on time $t\geq 0$ and polar coordinates $r\in(0,\infty),\theta\in(-\infty,\infty)$. In the following, we solve the system Eqs. (8-13) under various boundary and initial conditions and discuss exact and asymptotic solutions.

3 Linearized model

In this section we obtain exact solutions for the linear approximation of the conservative two-fluid model. A quick look at the linearization of Eqs. (12, 13)

[TABLE]

can give already a hint about the geometry of the solution’s patterns. For any of the coefficients of the double Fourier series expansion in time ($i\omega t$), and in the polar angle ($im\theta$) of the solution, we have

[TABLE]

The left hand side describes the equation for the stream lines expressed in polar coordinates by $r=R(\theta)$

[TABLE]

If the flow follows, for example, a logarithmic spiral ($R=R_{0}\exp{g_{0}\theta}$) we obtain

[TABLE]

which indeed represents a one-arm logarithmic spiral pattern in the ice density.

It is natural to ask whether distribution of sea-ice in the form of a rotating spiral is in a state of stable equilibrium. Instability can take the form of a warping of the spiral shape or dispersing the spiral pattern into a more uniform distribution of sea-ice. To obtain the states of equilibrium we consider the linearization of Eqs. (8-10) and we expect these linear solutions to generate wave type of equations and dispersion laws from where we can calculate the phase and group velocity, the latest being responsible for the dynamics of sea-ice patterns. Then, as in all stability problems, we consider the small perturbations from the equilibrium state. In the linear approximation spiral patterns are possible only through the effect of interacting waves of considerably different scales [41]. Consequently, we choose to use the method of multi-scale expansion in which we express the dependent variables in asymptotic formal series with respect to two small, dimensionless parameters $\varepsilon$ and $\delta$ as follows

[TABLE]

[TABLE]

where we are making the hypothesis that the zero orders of all quantities are time and position independent, and the velocity field has negligible zero order terms. These assumptions are validated by geophysical observations, [9, 10, 15, 12, 2, 3, 5], showing that such huge spiral sea-ice structures rotate slow and rather like a rigid pattern.

In order to separate different physical space-time scales in the mixed flow we are using the method of scaled parameters, [32], for the independent variables to $r=R\zeta$ and $t=\Theta\tau$, using $\zeta,\tau$ as dimensionless independent variables. Since the water flow is faster than the motion of ice fragments we expect the higher order corrections in the water flow field to dominate the corrections of the IFF density field, which implies $\varepsilon<\delta$, so we can choose $\mathcal{O}(\epsilon)=\mathcal{O}(\delta^{2})$ without any loss of generality. The ice fragments move together with the spiral pattern with the group velocity $v_{gr}$. It is simple to demonstrate that for this system the group velocity is always smaller than the phase velocity $\hbox{max}||\vec{v}||\simeq v_{gr}<v_{ph}$. Indeed, a simple estimation for the phase velocity for the pattern waves can be given by $v_{ph}\simeq\sqrt{\partial P/\partial\rho}$. From Eqs. (5-7) it results $v_{ph}\simeq\sqrt{gT_{0}}$. For regular size ice fragments floating freely we can assume $T_{0}\sim 10-100$ m which sets the phase velocity in the range $v_{ph}\sim 20$ m/s $\gg v_{gr}\simeq 0.01$ m/s. On the other hand, the phase velocity obtained from the linearized system of equations is of order $v_{ph}=\omega/k\sim R/\Theta$. From these observations it results $|u|\sim|v|\sim v_{gr}<v_{ph}$. From this scale hierarchy it results the balancing of same order types of terms

[TABLE]

By using the material conditions of our model Eqs. (5-7) in Eqs. (8-10) in absence of viscosity, we can re-write the system Eqs. (8-10) in the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the first order in $\varepsilon,\delta$ Eqs. (16-18) become

[TABLE]

[TABLE]

[TABLE]

where we kept the terms in the same order of magnitude according to Eq. (15).

In order to linearize the system Eqs. (19-20) we note that $\theta$ and $t$ are cyclic variables since they do not appear explicitly in the equations, so there are two constants of motion related to these two variables. Moreover, the equations are invariant to Galilean rotations around $\theta$ so the solution should depend on these variables only through phase expressions $m\theta+\omega t$, $m$ being integer in order to satisfy periodicity. From the rotational symmetry of the equations we expect that some of the solutions exhibit the same symmetry and provide rotational patterns in rotation. Any steady linearized solution in the first order has uniform density distribution with any divergence free velocity field. From the integrability conditions in the first order, $\Phi^{1}_{r\theta}=\Phi^{1}_{\theta r}$, and by using Eqs. (19-20) we re-obtain the vorticity $\vec{w}=\nabla\times\vec{v}=\vec{k}w$ conservation law in the first order

[TABLE]

conservation validates our choice for the smallness orders in the model. We differentiate Eqs (19) with respect to $\tau$, and differentiate Eqs. (20) with respect to $\zeta$ and $\theta$ respectively, and plug the results into the $\tau$ derivative of Eqs. (19). Eq. (19) decouples the density function $\Phi^{1}$ from velocity field and, back in dimensional variables, takes the form of a wave equation in cylindrical coordinates

[TABLE]

where the constant in front of the second term

[TABLE]

if positive, it can be interpreted as the phase velocity of the system’s linear waves. Eq. (21) can be solved by separation of variables, and the general solution has the form

[TABLE]

where $\omega=\omega_{r}+i\omega_{i}$ is a separation parameter with its real part $\omega_{r}$ representing the angular frequency of the waves, $m$ is an arbitrary integer, and $Z_{m}$ is a linear combination of Bessel functions of first or second kind of integer order $m$. In principle the solution must be regular at the origin $r=0$, which would limit the general solutions of Eq. (23) to be only Bessel functions of the first kind $J_{m}(\omega r/c)$. However, the real sea-ice swirl structure cannot be covered by the linearized model at its very central part where the ice always clusters into a more or less solid fragment, and hence the mixture cannot be treated as an idea compressible fluid. We therefore admit a singularity at $r=0$ in the solution Eq. (23) and request the solutions to describe a small radius $r_{0}$ ice core at the center of the swirl, see arrow 2 in Fig. 1. If $\omega$ is real and the water pressure exceeds the ice pressure ($P_{w}>\rho_{w}gT_{0}$) the argument of the Bessel functions is real, the oscillations of the radial function, and hence spiral shapes are allowed in the asymptotic region. If the water pressure is less than the ice pressure, we have $c^{2}<0$, the radial part of the solution becomes modified Bessel function, and oscillations are absent. The solution for IFF approaches asymptotically a uniform distributed value $\Phi_{\infty}\in[0,1]$ towards infinity, definitely outside the perimeter $L$ of the swirl, see see arrow 1 in Fig. 1. In other words the parameter $\omega$ will be determined from suitable chosen boundary conditions for $\lim\limits_{r\rightarrow 0}\Phi\sim\Phi(r_{0})=1$ and $\lim\limits_{r\rightarrow\infty}\Phi\sim\Phi(L)=\Phi_{\infty}$, where $L$ will be a suitable chosen parameter describing the radius of the swirl. These boundary conditions for our solutions request a discrete spectrum for $\omega$. Physically, this means that we are setting an exact solution as a representation of the swirl part of the mixture, and we leave the dynamics near the center, where ice concentration is larger and maybe even three-dimensional, to adjust itself to almost any requirement of the swirl part. We also expect that all the perturbations would decay to zero at infinity in a smooth manner.

With these comments being said, the the radial part of the solutions of Eq. (21) will be denoted generically by any Bessel function $Z_{m}$ and specified in more detail when is needed. These solutions involve damped wave oscillating evolution in the radial direction of the IFF, since in the far asymptotic range $r\rightarrow\infty$ the Bessel functions approach $1/\sqrt{r}$ times a periodic oscillating function in argument. This asymptotic behavior results in the coupling of the radial periodicity of the asymptotic expression of the Bessel solutions with the natural periodicity in the azimuthal angle and with time, hence generating solutions $\Phi^{1}$ with circular pattern symmetries, including spirals

[TABLE]

This asymptotic solution represents an Archimedean with $m$ arms, and these are trailing spiral arms if $\omega_{r}/c<0$, and leading spiral arms if $\omega_{r}/c>0$ [13]. The spirals represent stable modes if $\omega_{i}>0$. The Bessel solutions are bounded on $[r_{0},\infty)$ matching the necessary range for $\Phi\in[0,1]$, and they always decreases asymptotically to zero towards $r$ approaching infinity, meaning that far away from the center of the spiral the sea-ice degenerates in constant and uniform distribution, rather controlled at this point by random dynamics.

In order to build the general solution of the linearized system Eqs. (19-20) we use linear combinations of Eq. (23) summed over all ranges of parameters $m,\omega$

[TABLE]

These integral form generates solutions for given initial conditions, boundary conditions, or regularity conditions at $r=r_{0}$ and $r=L$ or even $r\rightarrow\infty$. Actually, even this general solution can describe spirals when the amplitude of the Bessel functions varies slowly with the radial distance, and if the phase varies quickly (large $\omega_{r}$ values). Indeed using the integral representation for Bessel functions, [16], we have

[TABLE]

which is a linear combination of Archimedean spirals of equations

[TABLE]

We built solutions within a disk of radius $L\gg r_{0}$. Physically, further away from the boundary of this disk $L$ we have only a uniform and sparse mixture of water with less and less ice, so we can impose the boundary condition $\Phi^{1}(L,\theta,t)=0$. Under this assumption, and inspired by Eq. (25), we can consider that the solution in Eqs. (23) form, for each $m$, an orthogonal system of functions labeled by $n$, complete over the space of continuous functions defined inside the disk $L$, and for $t>0$ in the form

[TABLE]

where $\xi_{m,n}$ is the root of order $n\geq 1$ of the Hankel function of the first kind and order $m\geq 0$, i.e. $H_{m}(\xi_{m,n})=0$. The order of the Hankel function is not relevant for the geometry of the spiral. We chose Hankel functions over other types of Bessel functions because these ones have the desired asymptotic behavior to generate spirals. The completeness of this orthogonal system of Bessel functions reduces the integral Eq. (25) over frequencies to a sum over the Bessel function roots, Eq. (27).

With this basis of functions we can solve the initial condition $\Phi^{1}(r,\theta,0)=\Phi^{1}_{0}(r,\theta)$ problem for the system Eqs. (19-20), by expanding the initial condition function in the a double Fourier and Fourier-Bessel series and thus obtain the continuous solution in the form

[TABLE]

[TABLE]

The basis of functions Eq. (27) built based on the boundary condition $\Phi^{1}(L,\theta,t)=0$ plus the initial condition provide the general solution for the IFF density field $\Phi^{1}$ in the first order linearized case.

From Eqs. (20) and (28) we obtain the linearized velocity field

[TABLE]

[TABLE]

All the partial waves in Eq. (27) have linear dispersion relations with $\omega_{radial}=ck$ and $\omega_{azimuthal}=c\xi_{m,n}k/(mL)$, so the solution Eq. (28) represents a linear combination of different rotational waves in the $\theta$ direction, and same type of waves in the radial direction. For large values of $r$, because of the asymptotic form of the Hankel functions, the azimuthal frequencies tend to be the same for any type of partial wave, $\xi_{m,n}\rightarrow m\pi/L$ so $\omega_{azimuthal}\rightarrow c\pi/(2L)$. That means that for large $r$ the solutions of Eqs. (19-20) the solution tends to become coherent in space-time, describing large-scale collective and coherent patterns of ice and water, including various number of arms spirals.

A solution can also describe a partial wave, with given $m,n$, and then the IFF field has the form

[TABLE]

The velocity field associated with this density field is given by

[TABLE]

and

[TABLE]

The properties of the general solution Eqs. (28-30) are in agreement with all the properties of observed patterns of ice. It is straightforward to check the IFF function fulfills the linear integrability condition, namely $-rw_{z}=u_{\theta}-(rv)_{r}=$constant. We also note that the only way to have steady motion solution is to cancel the coefficient in front of the time exponential. The only Bessel function fulfilling the case is for $m=0,n=1$ which, when implemented in the first order solution it reduces it to zero.

4 Time evolution of the Archimedean spiral

In this section we study the time evolution of some particular initial conditions for ice distribution in order to establish what types of patterns generate long time stable solutions and are more likely to develop and maintain in the sea. The most interesting pattern is the Archimedean spiral described in polar coordinates by the equation $ar+m\theta+\psi=0$, valid for any real $a$, any integer $m$, and any real parameter $\psi$ between [math] and $2\pi$. The Archimedean spiral has constant pitch $2\pi m/a$, that is the radial distance between successive turns. The integer $m$ describes the number of arms of the spiral, and the last term describes a rotation of angle $\psi$. If, for example, $\psi=cat$, the spiral rotates in time like a rigid pattern with constant angular speed $ac$. Moreover, any function defined in the plane which depends on space coordinates and time only through the variable $\chi=ar+m\theta+v_{0}t$ has contour levels with Archimedean spiral shape. Of course the function chosen in this study must be continuous and periodic in $\chi$, otherwise we will have multi-valued representation of the plane. Example of such function can be sine/cosine trigonometric functions, Bessel functions, hyperbolic functions, etc. We consider that a field defined on the plane has contour levels following Archimedean spiral if it can be written as $\Phi(r,\theta,t)=g(r)P(ar+m\theta+v_{0}t)$ with $g$ a decreasing function approaching [math] at infinity, and $P$ a continuous periodic function.

In the following, we discuss what are the initial conditions that can generate solutions with Archimedean pattern shape, and also the linear stability in time of such solutions. The radial extension of such observed spiral is very large (tens to hundred of Km) so in our calculation we can take at different stages the limit $L\gg R$ and even $L\rightarrow\infty$ which are realistic choices on one hand, and useful shortcuts in calculations since the integral over radius can be approximated with a Fourier transforms. The time scale of the phenomenon is also very large since the spiral pattern was observed moving very slow and persisting long time. Since in the absence of storms or strong sea currents, the distribution of ice fragments can be considered random, which means the IFF function is a constant. The sudden formation of such an Archimedean spiral pattern represents an interesting phenomenon of large scale collective behavior of ice fragments and water. Such coherent stable structures are possible only through a nonlinear dynamics which allows many scales interaction. Nevertheless, the linear solutions obtained in the previous section form a basis of complete functions, so we can expand any hypothetical nonlinear solution in series of the linear solutions and manage the series coefficients to assure stability of some nature.

Let us consider a standard Archimedean spiral $\Phi_{0}=\cos(a_{0}r+m_{0}\theta)$ as initial condition. The solution in Eq. (28) will retain only the term with $m=m_{0}$. In order to obtain the coefficients of the remaining series over $n$ we have to integrate terms of the general form

[TABLE]

Such integrals are exactly the coefficients of the Fourier-Bessel series in $H_{m_{0}}$ for the initial condition

[TABLE]

The solution for this initial condition, Eq. (28), becomes

[TABLE]

In the limit $L$ a very large number, and taking into account only the real part of the exponential involved, we can approximate the domain of integration with the positive real semi-axis, and for the even cosine function the numerator in Eq. (34) can be approximated with

[TABLE]

where $\hat{H}_{m_{0}}(\zeta)$ is the Fourier transform of the corresponding Bessel function. This Fourier transform is a polynomial of order $m_{0}$ in $\zeta$, with support only within the unit symmetric interval $\zeta\in[-1,1]$ and zero in the rest. Consequently, its derivative becomes the sum of two delta-Dirac distribution, out of which only $\delta(\zeta-a_{0}L/\xi_{m_{0},n})$ matters since $\zeta>0$. It means that in the sum over $n$ in Eq. (35) only the terms having $n$ given by the solutions of the transcendental equation $a_{0}L=\xi_{m_{0},n}$ have the dominant contribution, and let us denote this solution for $n$ by $n_{0}$. The rest of the terms with $n\neq n_{0}$ can be neglected, especially for very large values of $t$ since their contribution is densely chopped by the high frequency oscillations. The resulting solution for Eq. (35) has the form

[TABLE]

which in the asymptotic limit of very large distance, where the Bessel function can be approximated with cosine we have

[TABLE]

which means that any partial wave initial condition will be linearly stable in time, and it will keep its shape. Nevertheless, such solutions are nonphysical because the spiral extends to infinity and does not decay and smoothly connect to the surrounding water. We need to multiply such initial solutions with a decaying factor with respect to $r$.

4.1 Time evolution of partial waves

The linear stability of the partial waves can be demonstrated by using the asymptotic expansion of the Bessel function towards infinity, and by choose the initial condition the very same partial wave. In this case the coefficients $T_{m,n}$ are proportional to the expressions

[TABLE]

where $C_{1,2}$ are normalization constants and erf is the error function. The first factor in the above evaluation of $T_{m,n}$ has its Fourier transform the unit step function with support $|m_{0}-m|\leq 2\pi$ so it is only relevant in the IFF series if the summation index $m$ is in a $\pm 2\pi$ neighborhood of $m_{0}$. The first term in the parenthesis is subjected to the same behavior. Its contribution is relevant only if $|\xi_{m,n}-a_{0}L|\leq 1$ which means only those $(m,n)$ pairs (actually the $(m_{0},n)$ pairs) who fulfill this inequality are dominant. Finally, the last term in the parenthesis has a $\delta-$Dirac behavior. In a neighborhood of $x\sim 0$ it behaves like $x^{-1/2}$ and it has an exponential decay towards infinity, as one can see from its Taylor series around any point $x_{0}\neq 0$. These arguments demonstrate that out of the whole series from Eq. (28) only $m_{0}$ term and other very few $n$ terms fulfilling the condition $|\xi_{m_{0},n}-a_{0}L|$ are relevant. It means that towards large time values practically the solution reconstructs the initial condition partial wave.

4.2 Time evolution of the asymptotic solution

If we consider the limit $L\rightarrow\infty$, which is a natural limit given the large scale of observed circular ice formations, [9, 10], the solution Eq. (28) has the asymptotic expression

[TABLE]

where we used the asymptotic representation of the Bessel functions $J_{m}(x)\rightarrow\sqrt{2/(\pi x)}\cos(x-m\pi/2-\pi/4)$, and where $c_{m}$ represent the complex Fourier series coefficients of the function in the parenthesis (as a function of $\phi$), and $\mathcal{F}_{c}$ is the Fourier cosine of the argument as a function of $s$ evaluated at $\xi_{m,n}/L$. That is

[TABLE]

This result is a very good tool to classify the stability of different initial spiral configurations. According to the definition of a Archimedean spiral pattern in our physical system in the beginning of section 4, the initial distribution of ice fragments has the general form

[TABLE]

with $g(r)$ a decreasing function describing the diffusion of the spiral towards the boundaries into the sea, and $P$ o periodic function which governs the geometry of the spiral and its rotation. To be able to model we will consider for $g$ various powers of $1/r^{\alpha},\alpha>0$ and for $P$ a cosine. From the asymptotic expansion formula Eq. (37) we calculate the distribution of the coefficients of the series. There are only two situations of stable spirals. One is if $\alpha=0$, that is uniform distributed spiral. In this pattern the ice fragments do not decrease their density, or the width of spiral arms, larger $r$. The other case of stability is when $\alpha=1$ and the spiral radially decays as $1/r$. In the first case the Fourier cosine transform in Eq. (37) is applied to $r$ times trigonometric function. In the second case the transform is applied just to the trigonometric function. In both these situations the Fourier cosine transform result is the derivative of a $\delta-$Dirac distribution, or the $\delta-$Dirac distribution times constants. The result of the transform is a constant times $\delta(\xi_{m,n}-aL)$. This result restrict the sum over $n$ to only one term $n_{0}$, namely the solution of the transcendental equation $\xi_{m,n_{0}}=aL$, for each $m$. The next step is the calculation of complex Fourier coefficient $c_{m}$, but this turns to be trivial, since by orthogonality, the only contributing term filtered from the series over $m$ is $m_{0}$. This means that the time dependent solution $\Phi^{1}$ maintains the pattern of the initial condition, and stabilizes the Archimedean spirals of this types. It results that out of the double series in Eq. (37), only the term with $m=m_{0}$ and $n=n_{0}$ are contributing.

More importantly, this finding means that the only Archimedean spirals of the types

[TABLE]

which fulfill the special relation between their pitch and radial extension $L$ (size of the swirl)), falling in one of the zeroes of the Bessel functions, and given by

[TABLE]

are linearly stable.

Other types of Archimedean spirals, of the class we investigate here, may have a different type of decay law of the form $1/r^{\alpha}$ with $\alpha\in(0,1)$, which includes $1/\sqrt{r}$. These radial dependence generate Fourier cosine transform with rational functions dependence of the variable $\zeta$, with singularity point also at $\xi_{m_{0},n_{0}}=aL$. The Fourier cosine spectrum is concentrated in this singularity as in the previous cases, but there are also tails which mix the initial condition in time with other partial waves. The initial spirals are weakly stable, and tend to slowly disperse in time, without actually completely erasing the spiral pattern. Among these situations we have the decay type of the asymptotic expression for the Bessel functions as $1/\sqrt{r}$. This result explains why single partial waves, involving one Bessel function, cannot form sharp stable spirals and tend to mix modes and deform the pattern.

An example showing five types of decay law with $r$ are presented in Fig. 2 at initial moment of time (Archimedean spirals of the same pitch), and in Fig. 3 at a later moment in time, when the spiral with too slow or too steep decay law disintegrated. Finally, Archimedean spirals with steeper decay with a dependence of the form $1/r^{\alpha}$ with $\alpha>1$ generate a wide band Fourier cosine spectrum, and some higher power generate even divergent integrals. All these steeper decay laws are unstable and disperse quickly in time.

The results presented in this subsection demonstrate that, at least in the linear approximation of the dynamical equations and in the inviscid case, only some selected Archimedean spiral patterns tend to be stable for large scales of space and time.

4.3 Time Evolution of Archimedean spirals

In this subsection we choose for initial condition more general form of Archimedean spirals, namely

[TABLE]

where the function $h$ must decrease with $r$ and $g$ is periodic. We plug this initial condition in the general solution (not the asymptotic approximation) Eq. (28). The double series coefficients contain the form

[TABLE]

We make a change of variable $q=as+m_{0}\phi$, use the fact that $g$ is $2\pi$-periodic to translate the limits of integration with respect to $q$, and the above integral becomes

[TABLE]

We can expand the function $h$ in its Taylor series around a point between [math] and $L$, make the same supposition that $L$ is large enough to consider the upper limit of the integrals in $s$ as infinity, and consequently the inside integral is just a sum of inverse Fourier transforms of a power law time a Bessel function. These inverse Fourier transforms output as function of $\zeta$ the product between a finite support step function, and a polynomial. The integral over $\phi$ also reduces to a $\delta-$Kronecker symbol. We have

[TABLE]

where $H$ is the Heaviside distribution, $\delta$ is the Kronecker symbol, and $\mathcal{P}_{o(m)}$ is a polynomial of order $o(m)$ which is $m$ if $m$ is even, and $m+1$ is $m$ is odd. The hat means the Fourier transform. This term filters out of the double series the terms with $m=m_{0}$ and $n>n_{0}$ where $n_{0}$ is the solution of the transcendental equation $\xi_{m_{0},n}>aL$, exactly as in the result obtained in subsection 4.2. The solutions Eq. (28) for the initial condition $\Phi^{1}(r,\theta,0)=h(r)g(ar+m_{0}\theta)$ becomes

[TABLE]

We would like to evaluate the geometric parameters of the Archimedean spiral in a real situation. In all previous subsections we noticed that the solution series begins with a first term for $n=n_{0}$. Since the series must be convergent, we know the first term is the dominant one. Another noticed fact is that the real case spiral has one, maximum two arms, so $m_{0}=1$ or $2$. The core of the spiral, which usually is ice so $\Phi(0)\sim 1$ can be considered to extend radially from origin up to the first zero of the Bessel function involved in the solution. We can obtain the radius of the core from the equation $J_{m_{0}}(\xi_{m_{0},n_{0}}r_{c}/L)\simeq J_{m_{0}}(\xi_{m_{0},1})$ so we have $r_{c}\simeq\xi_{m_{0},1}L/\xi_{m_{0},n_{0}}$. We can take the radius of the core to be at half way between origin and first zero, so $r_{c}\simeq\xi_{m_{0},1}L/(2\xi_{m_{0},n_{0}})$. Also from the dominant term we can roughly estimate the radius of the swirl. Let us assume that the swirl ends when IFF becomes $\epsilon\Phi$ evaluated at the core. From here we obtain $L\simeq\xi_{m_{0},n_{0}}/(2\epsilon^{2})$. The pitch can also be evaluated as pitch$=2\pi m_{0}/\xi_{m_{0},n_{0}}\simeq(2\pi\xi_{m_{0},n_{0}})L$. The documentation shows spirals with $5-10$ turns, so $L/$pitch $\simeq 5-10$. If we choose $\epsilon=0.05$ we have $\xi_{m_{0},n_{0}}\simeq 2\pi m_{0}L/$pitch which ranges $\xi_{m_{0},n_{0}}$ between $30$ and $120$. It results for $m_{0}=1$ or $2$ the possible range $n_{0}=10-38$. From here we can estimate $L\simeq 510^{3}\xi_{m_{0},n_{0}}\simeq 50-230$Km with core $r_{c}\simeq(\xi_{m_{0},1}/\xi_{m_{0},n_{0}})L\simeq 2-40$Km, of course function of the radius of the swirl.

5 Limiting nonlinear modes

There are three limiting types of axially symmetric flows (that is neglecting plane density waves at the surface) for the ice and water mixture: the radial mode when we can neglect azimuthal flow and rotations, and the whole system performs radial oscillations of compression and dilation; the rotational mode in which the system is in rotational flow and we neglect radial components, and the spiral mode when the flow is concentrated mainly along a certain spiral direction field.

5.1 Radial limiting mode

In this situation we neglect the azimuthal velocity $v$ terms in Eqs. (8, 9) and, with the notations $\Psi=1-a\Phi$, $b=(P_{w}/\rho_{w}-gT_{0})/a$, we have

[TABLE]

[TABLE]

while the last Eq. (10) reduces to $\Psi_{\theta}=0$. It is the most natural in this case to assume that a the fields ($\Psi,u$) have radial symmetry. With the notation $\bar{u}=\Psi u$, and multiplying the second equation with $\Psi$ we have

[TABLE]

[TABLE]

equations representing conservation laws for surface mass density and mass flow

[TABLE]

respectively, under the assumption that the radial velocity decreases towards infinity more rapid that $1/r$. From a brief nonlinear dispersion relation analysis, we note that the system of equations is scale invariant to the magnitude of $\Psi$, the mean value of the radial velocity ranges $<u>\sim\sqrt{b}=\sqrt{(P_{w}-\rho_{w}gT_{0})/(\rho_{w}-\rho_{i})}$ and the radial oscillation frequency ranges $\omega/(2\pi)\sim<u>/R\simeq\sqrt{b}/R$, where $R$ is the radius of the ice-water mixture.

The radial mode system of equations in the form Eq. (40, 41) is separable, and the resulting PDE for $\bar{u}(r,t)$ has a very busy expression

[TABLE]

where $\Psi_{1}$ is an arbitrary function of $t$ only, and $Y$ is a nonlinear expression in $\bar{u}$ and its derivatives

[TABLE]

[TABLE]

[TABLE]

Nevertheless, since the first equation in the system Eq. (40) is originally linear, we apply perturbative methods to the above system depending only on the two unknown functions $\Psi,\bar{u}$, namely $\Psi\simeq\Psi_{0}+\Psi_{1},\bar{u}\simeq\bar{u}_{0}+\bar{u}_{1},\bar{u}_{r}\simeq\bar{u}_{0r}+\bar{u}_{1,r}$. We linearize of Eq. (41), differentiate it to $t$, and couple it with Eq. (40) differentiated to $r$, and we obtain in the zero order $\bar{u}_{0}=\sqrt{b}\Psi_{0}$, and in the first order

[TABLE]

By considering the Fourier decomposition and substituting $\partial_{t}=i\omega$Id, we obtain from Eq. (43) a generalized hypergeometric differential equation. The solution reads

[TABLE]

where $U$ is the confluent hypergeometric function, $L$ is the generalized Laguerre function, $C_{1,2}$ ar constants of integration, and we have the parameters

[TABLE]

[TABLE]

[TABLE]

To simplify further this solution, we can neglect the contribution of the term $\bar{u}_{1}\bar{u}_{1,r}$ from Eq. (41) and the corresponding simplified solution is expressible in term of the modified Bessel functions

[TABLE]

The calculations are complete by implementing the above solutions for $\bar{u}$ in Eq. (40), integrating with respect to time to obtain the IFF field $\Psi$

[TABLE]

and substituting this $\Psi$ solution $u=\bar{u}/\Psi$ to finally obtain the radial velocity.

From the linear stability analysis of the system Eqs. (40, 41) we obtain that always one of the linear equilibrium solutions is asymptotically stable, and the second always unstable. This fact however, does not guarantee the nonlinear stability of the system all together. For a large range of parameters the IFF solution is always around unit at origin and decreases parabolic towards large distances, while the radial velocity is zero at origin and increases almost linear towards the periphery of the ice disk, and both oscillate coherently in the radial direction. The resulting radial mode consists in very slowly periodic radial compression and dilation with the rate of change increasing with distance to center. Typical value are $\Phi(r=30\hbox{Km})\sim 0.85,u\simeq 100\div 2000$ m/day, and $2\pi/\omega\simeq 6\div 40$ hours.

5.2 Rotational limiting mode

In this limiting mode we neglect radial velocity $u\simeq 0$ and Eqs. (8) amd (10) acquire the simplified form

[TABLE]

[TABLE]

while Eq. (9) can be integrated into $v=\pm\sqrt{br\Psi_{r}/(a\Psi)}$. Such solution exists only if $\Psi_{r}\geq 0,\Phi_{r}\leq 0$ meaning that ice always accumulates at the center of the rotating pattern. By differentiating the equations above to $\theta,t$ respectively we obtain a nonlinear equation for $\Omega=\ln(\Psi)$

[TABLE]

This equation is similar to the Nonlinear Schrödinger equation, but it contains several other terms. A simple solution for this equation is the rigid rotation of the whole pattern, $v=\omega_{0}r$, which is also a solution of Eq. (45). In this case Eq. (46) can be easily integrated given the observation that $b\ll ar^{2}\omega_{0}^{2}$ and finally generates

[TABLE]

where $P$ is an arbitrary function. This rotational pattern and $r$ dependence confirms the accumulation of ice towards the center. While looking for other possible solution we notice that in the long range, the asymptotic form of Eq. (47) for $r\rightarrow\infty$ would have a blow-up solution $\Psi\sim t$. The equilibrium solutions $\partial_{t}=0$ of the system Eqs. (45, 46) represent un-physical situation. Indeed it results $(\ln\Omega_{\theta})_{r}+2\Omega_{r}+r=0$ and hence the corresponding solutions $\Psi\sim\sqrt{\theta}\exp(-r^{2})$ are not periodic in $\theta$, hence are strongly unstable.

5.3 Spiral limiting mode

If the flow follows the geometry of a family of time-independent logarithmic spirals of equations $R(\theta)=R_{0}\exp{g_{0}\theta}$ we have the simplifying relation $u=R_{\theta}v/R=g_{0}v$. With this substitution for $u$, by dividing Eq. (12) by $g_{0}$ and subtract it from Eq. (13) we obtain

[TABLE]

or

[TABLE]

If we assume the simplest periodic dependence of $\Psi$ by the azimuthal angle such that $\partial_{\theta}\Psi=im\Psi$ we have after one quadrature

[TABLE]

This is an integral representation of the solution, but it reveals a stable dependence on the logarithmic spiral argument, modulated with an exponential depending only on the tangent velocity. Most importantly, since the multiplicity $m$ factors out, only one-arm logarithmic spirals are possible

From the remaining equations Eq. (11) and Eq. (12), also under the substitution $u=g_{0}v$, and expressed as time Fourier transform ($\partial_{t}=i\omega$) we obtain

[TABLE]

[TABLE]

Once we proved that the ice pattern follows a stable logarithmic spiral, from here the procedure to obtain the velocity field is direct. By replacing $(\Psi v)_{r}$ in Eq. (52) from its expression in Eq. (51), we obtain an algebraic equation in the functions $\Psi$ and $v$

[TABLE]

We solve this algebraic equation for $v=v(\Psi,r)$ and substitute this solution in Eq. (50) transforming this equation into an integro-differential equation for $\Psi$ only, and with $\Psi$ obtained find $v$ and $u$.

As a final observation in this section, we mention that we can find a very stable and simple form for the logarithmic spiral solution in a particular equilibrium situation. If we assume $b\ll v^{2}$ in Eqs. (11, 13), both these equations become identical to

[TABLE]

which generates a solution in the form

[TABLE]

This is again the structure of a stable one-arm logarithmic spiral. The condition for small $b=P_{w}/\rho_{w}-gT_{0}$ represent a case when the phase velocity of the density waves through the sea-ice pattern is way smaller than the group velocities, and actually means that the water pressure is balanced by the ice buoyancy, and there is no strong interaction and exchange of momentum between ice and water. It also means that the gradient of pressure terms in the Navier-Stokes equations are negligible, which means that logarithmic spiral patterns tend to form and stabilize when the sea-ice is rather isolated from currents or winds.

6 The nonlinear equations

Eq. (11) has exact time-independent solutions generating logarithmic spiral patterns

[TABLE]

where

[TABLE]

and $g_{0},r_{0},q_{0},f_{0}$ are arbitrary real parameters and where the value of $a$ must be adjusted such that $\Phi\in[0,1]$. The IFF distribution continuously decreases from center towards its asymptotic value

[TABLE]

so the boundary conditions around the spiral determine the parameter $r_{0}\sqrt{f_{0}}$. We can choose a minimal radius $r_{min}$ above which the model validates, and since in all observations the center of the spiral appears to me solid ice we can determine another condition $\Phi(r\simeq 0)\sim 1\sim r_{0}q_{0}\sqrt{a(g_{0}^{2}+m^{2})}/(g_{0}r_{min}\sqrt{b})$. The ratio between the radial and azimuthal velocities is a constant determining uniform flow along spiral curve coordinates, $r=r_{0}\exp(m\theta/g_{0})$. The velocity is also independent of time and azimuthal angle.

Solution Eq. (54) is almost an exact solution for the momentum conservation equations, Eqs. (12, 13), within an approximation of order $\mathcal{O}(\sqrt{b}/r^{2})$ which is about $10^{-9}$ times smaller than the order of the $\Psi$ solution, for example, within a large range of parameters values. The parameter $r_{0}$ was introduced to allow dimension free argument in the logarithm, but also controls the smallness of the defect of this solutions in Eqs. (12, 13).

With the notations $b=(P_{w}/\rho_{w}-gT_{0})/a,\Psi=1-a\Phi$, and similarity substitutions $\xi=\ln r,s=t/r=t\ \hbox{exp}(-\xi)$, Eqs. (8-10) can be written in a coefficients-free symmetric form

[TABLE]

[TABLE]

[TABLE]

The system can be written in an almost symplectic form by using the covariant gradient operator $\mathcal{L}=\partial_{s}+u\partial_{\xi}+v\partial_{\theta}$ along the flow vector field, and define the gradient $\nabla=(-b\partial_{\xi},-b\partial_{\theta})$ [39]

[TABLE]

Eq. (55) for mass conservation contains linear and quadratic nonlinear terms, while the last two equations (Eqs. (56, 57) for momentum conservation) involve also cubic nonlinearity. We follow this idea to identify solutions for the full nonlinear system.

6.1 Spiral patterns from the full nonlinear system

By differentiating Eqs. (56, 57) with respect to $\xi$, and $\theta$, respectively, by substituting the first term from each resulting equation into the $s-$differentiated Eq. (55) we obtain the nonlinear extension of Eq. (21)

[TABLE]

where the velocities occur only in the higher orders terms, meaning that Eq. (55) describes mainly the geometry of the patterns through $\Psi$, while Eqs. (56, 57) describe the dynamics of these patterns. Without any loss of generality we look for solutions appropriate for description of the ice patterns with circular of spiral symmetry in the form $\Psi=f(\xi)\hbox{Exp}[i\chi_{1}(\xi,\theta,s)],u=g(\xi)\hbox{Exp}[i\chi_{2}(\xi,\theta,s)],v=h(\xi)\hbox{Exp}[i\chi_{3}(\xi,\theta,s)]$ with $f,g,h,\chi_{1,2,3}$ real functions.

By implementing these forms in Eq. (59), we obtain from its real part

[TABLE]

Under certain legitimate hypotheses we can use this equation to find the phase $\chi_{1}$ for the shape function $\Psi$. The remaining imaginary part of Eq. (55), and the real/imaginary parts of Eqs. (56, 57) provide a system of 5 equations for the remaining functions $f,g,h,\chi_{2},\chi_{3},$ to be determined.

If the amplitudes $f,g,h$ vary slowly with $\xi$, while the phases $\chi_{1,2,3}$ vary quickly with $\xi$ (see for example the models for spiral galaxies [13, 21]) Eq. (59) generates solutions with spiral symmetry in the IFF distribution at any instant of time, the shape of the spiral being given $\chi(\xi,\theta,s)=$const. In order to demonstrate the existence of such spiral solutions, and in accordance with to the above hypothesis on the rate of change of various functions with $\xi$, we can neglect the term $(fg)_{\xi}$ with respect to the $\chi$ derivatives. Moreover, we can approximate in this equation $\chi_{1}\simeq\chi_{2}\simeq\chi_{3}$ denoted in the following $\chi$, because differences between the phase terms occur in the $4^{th}$ order. The above equation becomes

[TABLE]

Because of the $\theta-$periodicity, the general spiral solution described by the equation $\chi(\xi,\theta,s)=$const. must be in the form $\chi=G(\xi,s)-m\theta=$const., hence $\chi_{\theta}=m$, in polar coordinates, with $m$ integer describing the number of arms in the spiral. The spiral pattern is coherent in time only if the velocity field $(u,v)$ is directed along the tangent to the spiral pattern at any point, which involves the condition

[TABLE]

otherwise the spiral would radially dilate, shrink or disperse. For example, an Archimedean spiral is generated if $u/v\sim 1/r$, and a logarithmic spiral is generated if this ratio is constant. In the logarithmic spiral case, $h/g=a_{0}=$ constant, Eq. (60) provides an implicit solution $\chi(\xi,\theta)=G(\xi)-m\theta$

[TABLE]

This solution generates a distribution of logarithmic spirals in the plane with centers along a straight line, Fig. 4, of equations $G(\ln r)=m\theta$. The constant value $a_{0}$ for the ratio $h/g$ is related to the asymptotic values of $G\rightarrow G_{\infty}$ when $r\rightarrow\infty$, $a_{0}=\cot^{-1}G_{\infty}/(2m)$ which describes the maximum angle of the spiral where the pattern stops $\theta\leq G_{\infty}/m$.

Solution Eq. (61) also predicts the limiting size of the logarithmic spirals as depending on the ratio of the flow velocities

[TABLE]

Also, this solution can explain the occurrence of multiple separated ice swirls, as they were observed in the ocean, Figs. 5, [11]

In the following we present an algorithm to construct exact series solutions of the full nonlinear system Eqs. (55, 56, 57) by using the scaled series in Eq. (14). The procedure is rooted in the observation that in the linearized version, Eq. (21), the system reduces to a Bessel equation in polar coordinates. When we implement the iterative algorithm to the full nonlinear system, we notice that the same Bessel differential operator occurs recurrently in all orders of smallness, property which induces a diagonal structure to the iteration algorithm, hence validating its efficiency.

Beginning from Eq. (59) we obtain the following expansion in the smallness parameter $\epsilon$

[TABLE]

[TABLE]

where the Bessel linear operator $B$ has the following action

[TABLE]

and $\mathcal{N}_{1,2}$ are specific nonlinear operators acting on the arguments placed in front of them. For example

[TABLE]

[TABLE]

Under the previously made scaling assumption $\mathcal{O}(\epsilon)=\mathcal{O}(\delta^{2})$, the structure of Eq. (62) becomes self-consistent, and it can be solved by successive iterations following the Lindstedt-Poincaré method, [32], since for any $k\geq 2$ and $j\geq 1$, $\mathcal{O}(\epsilon^{k})\neq\mathcal{O}(\delta\epsilon^{j})$ are independent orders. Consequently, each term from each of the three main terms in Eq. (62) generate independent equations. Namely, each pair of terms with coefficient $\delta\epsilon^{k-1}$ and $\epsilon^{k},k\geq 2$ generate system of two nonlinear complex differential equations in only 5 dependent variables $\{\Phi_{k},u_{2k-2},u_{2k-1},v_{2k-2},v_{2k-1}\}$, since the other dependent variables occurring in these terms, $\{\Phi_{0},\dots,\Phi_{k-1},u_{1},\dots,v_{2k-3}\}$, were solved in the previous pair of differential equations, corresponding to the coefficients $\delta\epsilon^{k-2},\epsilon^{k-1}$. The only exception is the order $\mathcal{O}(\epsilon)$ where the corresponding nonlinear equation has the form

[TABLE]

and it can be solved for $\Phi_{1},u_{1}$ and $v_{1}$ by using the substitutions

[TABLE]

where all the three physical fields assume the same phase dependence imposed by the cylindrical (and later on spiral) symmetry. From the real/imaginary, and even/odd separations (i.e. $\cos,\sin$) of the terms of this equation we obtain exactly four differential equations for four real functions $\alpha,\beta,\gamma,\chi$. In the next step, the resulting term of order $\epsilon\delta$ provides the differential equation from where one can integrate $u_{2},v_{2}$ from $\Phi_{1},u_{1},v_{1}$. The next order, the term $\mathcal{O}(\epsilon^{2})$ generates an equation for $\Phi_{2}$ function of the previously obtained $\Phi_{1},u_{1},v_{1}$, the next term in order $\mathcal{O}(\epsilon^{2}\delta)$ provides an equation for $u_{3},v_{3}$ function of $\Phi_{1},\Phi_{2},u_{1},\dots,v_{2}$, and so on.

An important observation is coming from the linearized Eq. (21) for mass conservation, which depends only on the $\Phi$ variable. When we consider higher orders, and other nonlinear terms, no matter of the scaling or approximation procedure used, this equation will have the general form

[TABLE]

where $\triangle_{cyl}$ is the Laplacian in cylindrical coordinates, $\mathit{N}$ is a nonlinear operator in the model dependent variables and their derivatives. Without any loss of generality we can request the IFF function to have the form

[TABLE]

with real functions $f,\chi$ such that $f$ varies slowly with the radial distance $r$, while $\chi$ varies quickly. Independently of the procedure to find the exact expression for the right hand term in Eq. (67) we can always expand this term in a Taylor series in all its arguments, around zero. Since the zero and first order terms were already considered in the left hand part, the right hand term begins with quadratic terms, like for example $C_{0}\Phi^{2}$. When we implement the form of solutions from Eq. (68), after some algebra, the imaginary part of the resulting system of equations reads

[TABLE]

which, in the approximation of slow variation of the amplitude $f$ in comparison to the phase $\chi$, is the sine-Gordon equation [24, 36, 26]. The remaining real part of the equation is used to determined the function $f$, and the integration of the momentum conservation system of equations determines the velocities. This result is independent of any procedure used to decompose the system of equations, hence it has a very high degree of generality and validity.

Multiple turns spirals are obtained by using multiple-fronts kinks (topological solitons) solutions of the sine-Gordon equation, built by gluing together shifted kinks on plateau of arbitrary width of other base kinks, and so on. Such solutions were obtained in the case of nonlinear dispersion equations, [31] (sometimes referred to as kovatons, [30]) whose stability is the same as for regular kinks, since multiple-fronts kinks have the same minimum value of energy as kinks, being calculated from their Hamiltonian (see e.g. [24], pp. 150).

6.2 Spiral solutions in the quadratic approximation

In order to study quantitatively the nonlinear solutions and their stability for our model we retained in Eqs. (8-10) only the quadratic terms. The final results in this section amply justify this truncation, because we are able to obtain spirals of various geometries that fit very well the observed patterns in this order of approximation. We search solutions for the IFF in the form $\Phi=f(r)\hbox{Exp}[i(G(r)+m\theta+\omega_{c}t)]$ where $f,G$ are real functions, and we use complex frequency $\omega_{c}=\omega_{r}+i\omega_{i}$ to take into account dissipation too, $m$ is integer, and we denoted like before $c^{2}=(P_{w}/\rho_{w}-gT_{0})$. By integrating Eqs. (9-10) with respect to time

[TABLE]

[TABLE]

and by plugging these velocities in Eq. (8) we can separate three main terms. The real part of the time dependent terms has the form

[TABLE]

The imaginary part of the time dependent terms has the form

[TABLE]

and the time-independent terms are

[TABLE]

where we mention that the $\theta$ derivative here does not reduce to $m$ because the stationary solutions do not generate spiral patterns. Eqs. (71-72) form a nonlinear differential system for the functions $f(r),G(r)$, $u_{s}(r,\theta)$ and $v_{s}(r,\theta)$ and from its solutions we can build the fields $\Phi,u,v$ and eventually predict occurrence of sea-ice patterns formed on the water surface from our nonlinear model in the quadratic approximation, Eqs. (8, 9, 10), or in its the compact version Eqs. (55, 56, 57).

Eqs. (71-72) are similar to very well studied NLS-type equations in literature. A very simplified version of our radial model Eq. (71) was derived by Benney and Roskes, [26, 27], for the case of a $(2+1)-$dimensional NLS-type equation for slowly varying envelope surface water waves in potential flow and finite depth. The BR equation can explain long-time similarity solutions of NLS with radiative tails. A more general form for NLS-type integrable system was obtained by including surface tension effects to the BR equation shallow water, [28, 27]. In addition to general nonlinear integrable models for shallow water waves using the BR equation, [36], similar equations to Eq. (71) are also described in nonlinear optics, where they are referred as the Davey-Stewartson (DS) equation, [34, 35], or for $\chi^{(2)}$-nonlinear optics [37, 38]. Moreover, it was shown, [26, 37], that the DS equation applied to shallow water waves with large coefficient of surface tension admits localized boundary induced pulse solutions, or weakly decaying lump-type solutions, including the so-called dromions, exhibiting interesting nonlinear behaviors including wave collapse, blow-up and dark-lump solitons. These results can be somehow condensely expressed through a multidimensional ($d\geq 1$) generalized NLS equation analysis, [29], of the form

[TABLE]

It was shown that Eq. (73) admits subcritical global solutions $\Psi(r,\theta,t)=T(r)\hbox{Exp}[i(G(r)+\omega t)]$, for which blow-up does not occur, if $d\sigma<2$. With this substitution, the resulting equation becomes

[TABLE]

which has the so-called Townes modes as pulses solutions, if $2\sigma\leq 2$, [33]. The Townes modes can become stable solitons, [40], in the case of self-focusing NLS solution collapse and if their energy is smaller than the critical energy $2\pi\int_{0}^{\infty}rT^{2}(r)dr$, [33]. In our case, Eq. (71), by neglecting dissipation, $\omega_{i}<<\omega_{r}$, the coefficient $2\sigma\sim 1$ which fulfills the subcriticality condition of existence of global stable solitons.

Eq. 72 represents a divergence free condition for the stationary part of the velocity flow, namely $\nabla\cdot(u_{s},v_{s})=0$, coming from the incompressible mass conservation of the stationary part. Simplified solutions $f(r$) where $G(r)=\pm m$ would not generate swirl patterns, but rather oscillating monomials in $r$, because such a choice means neglecting a strong nonlinear coupling between the IFF function and the velocities. The term $1-a\Phi(r,\theta,t)$ occurring systematically next to the components of the velocities is not too relevant for the dynamics because the small values ($a\simeq 0.3$, $\Phi\simeq 0.5$) in front of the time oscillating part under the complex logarithm. The time variation of this term is in second order of smallness, so occasionally the term can be considered constant for rough qualitative evaluations. The consideration of complex frequency helps in the evaluation of the dispersion relations, and in the study of the stability of the nonlinear solutions as it will be seen in continuation. Actually, for zero damping effects $\omega_{i}=0$ Eq. (72) becomes identical to the DS (or BR) equations. In this ideal case we expect that the spiral solutions, in appropriately chosen system of coordinates become a representation of Townes modes lump-solitons, or rather lump-breathers [33, 40].

Because the amplitude $f$ and the phase $G(r)$ are intertwined in the system Eqs. (70, 71) we obtain a more convenient version where the amplitude $f$ is decoupled from the phase, and fulfills the equation

[TABLE]

[TABLE]

where $C_{1}$ is a constant of integration. Once Eq. (75) is solved for $f$, we can implement that solution into Eq. (75) and obtain the phase solution, i.e.

[TABLE]

such that the coupling is now reduced to the shared constant $C_{1}$. Eq. (75) is a strongly nonlinear integro-differential equation with singularities at $r,f=0$ very difficult to solve and even to characterize qualitatively. Nevertheless, we note that in the initial version of the equation for the amplitude $f$, Eq. (71), the phase enters only as $G^{\prime}$, so from the phase solution in Eq. (76) we notice that the control term towards $r\rightarrow\infty$ is the denominator $rf^{2}$. If at large distance from its center the spiral pattern is surrounded by ice $(f=1)$ or by a uniform mixture of water and ice, so asymptotically $f\neq 0$, it results that the phase $G$ approaches a constant asymptotic value, which would generate logarithmic type of spirals. However, if the spiral pattern is surrounded by water and $f\rightarrow 0$ towards the boundaries, the asymptotic behavior of the phase $G$ can be unstable and depends on how fast the amplitude $f$ approaches zero value. In order to insure spiral stability we need $rf^{2}\rightarrow 0$ which requests $f(r\rightarrow\infty)\sim\mathcal{O}(r^{-\alpha})$ with $\alpha\geq 1/2$. Even in this case, the fourth term in Eq. (75) can approach infinity because of its higher power $f^{3}$ in denominator.

In order to avoid this ambiguity in our analysis, we consider here only solutions with $C_{1}=0$. In this case Eq. (75) reduces to its first three terms in the LHS, which are nothing but a Bessel type of Sturm-Liouville linear differential operator, and in the RHS the square of the integral. With the substitutions $r^{2}=2q$ and

[TABLE]

we can map the integro-differential equation into a quadratic nonlinear differential equation

[TABLE]

We thus obtained a nonlinear ODE model through Eqs. (70, 71), or equivalently Eqs. (75, 76), or finally more particular system of order 3 through Eqs. (75, 77), which generate an exact solution in the quadratic approximation for the IFF field as $\Phi=f\hbox{exp}[i(G(r)+m\theta+(\omega_{r}+i\omega_{i})t)]$.

In the following we present an example of stable ($\omega_{i}=0$, no dissipation) Archimedean spiral pattern $G_{0}r+m\theta+\omega_{r}t=$const., as a asymptotic solution of Eqs. (71, 70). Indeed, for $G(r)=G_{0}r$ we have $G^{{}^{\prime}2}=G_{0}^{2}$ constant and Eq. (71) provides the solution $\Phi=f\hbox{exp}[i(G_{0}r+m\theta+\omega_{r}t)]$ with

[TABLE]

where $a_{1,2}$ are constants. In the asymptotic region, for large enough $r$, regime which is very easy to meet given the space extension of the spiral solutions, the Bessel functions in the solution Eq. (78) behave like $f\rightarrow 1/\sqrt{r}$, which makes $G^{\prime}\sim 1/(rf^{2})$ a constant, exactly as it should be to close the system. In a similar way we find a solution generating a logarithmic spiral $G1\log r+m\theta+\omega_{r}t=$const. By implementing $G^{{}^{\prime}2}=G_{1}^{2}/r^{2}$ in Eq. (71) generates the solution $\Phi=f\hbox{exp}[i(G_{1}\log r+m\theta+\omega_{r}t)]$ with amplitude

[TABLE]

In this later case, in the asymptotic region for large $r$ the amplitude approaches the constant $a_{0}$ and the expression $1/(rf^{2})$ approaches $G^{\prime}\sim G_{1}/r$ as requested.

We close this section by presenting an exact IFF amplitude $f$ (or $\varphi$) solution for the third order nonlinear Eq. (77) in the non-dissipation case, $\omega_{i}=0$. The absence of dissipation involves the dropping of the integral term in this equation, or equivalently dropping of the non-derivative term in the $\varphi$ equation.

[TABLE]

[TABLE]

$a_{0,1}$ being constants and with the phase obtained by implementing this amplitude $f$ in Eq. (76). One example of such solution is presented in Fig. 8.

In addition to the specific spiral geometries, the present model can predict the average size of the sea-ice spiral patterns, as well as the number of turns, or their pitch. The general equation of a time independent spiral can be written in the general form $\cos[G(r)+m\theta]$ or $\cos[r+mG^{-1}(\theta)]$ where cosine or any other periodic function can be used, the integer $m$ is the number of arms, its sign identifies the helicity of the spiral. The spiral equation in polar coordinates becomes $r=\pm mG^{-1}(\theta)$. For one arm spirals, the pitch is given by $p=G^{-1}(\theta+2\pi)-G^{-1}(\theta)$. For example, the Archimedean spiral $\alpha r+\theta$ has $p=2\pi/\alpha$, and the logarithmic spiral $\alpha\ln r+\theta$ has variable pitch $p=r(e^{2\pi/\alpha}-1)$. Independent of the mathematical approach of the model equations, fully nonlinear analysis, quadratic approximation, linearized, etc., we notice that the solutions for the IFF distributions are always combinations of Bessel functions, see for example Eqs. (23,24,63,78,79,80), consequence of the basis generated by the linearized equations. All solutions contain the general term

[TABLE]

where $\alpha$ stands generically for the coefficient in front of $r$ for various solutions. From Eq. (76), and considering just the dominant term in the asymptotic region (where the large part of the spiral actually unfolds) for the phase $G(r)$ function, we can write

[TABLE]

which implies the following spiral parameterization

[TABLE]

Since $A$ is of the order of unity because of the range for $0\leq\Phi\leq 1$, and stable solutions have usually $C_{1}$ negligible and very small, the limiting radius of such spirals is given by $r_{max}=L=\pi/(2\alpha)$, see Fig. 7 right, and the pitch is almost constant within the whole active range of $\theta$, $p\simeq\pi/\alpha$. From the linearized model, Eq. (27), we know the scaling factor in front of $r$ inside the Bessel functions is $\xi_{m,n}/L$. If we consider the fundamental mode $m=1$, which is probably the most stable, knowing that $\xi_{1,n}\sim 3.83n$ we can consider the mean value of the pitch $p=L/\mathcal{N}$, where $\mathcal{N}$ is the mean number of full turns of such a spiral. Hence we have $\alpha\sim\xi_{1,n}/L\sim\pi/p\sim\pi\mathcal{N}/L$ where from $\mathcal{N}\sim\xi_{1,n}/\pi\sim 1.2n$. It results the spiral have at least one full turn, and the maximum number of turns is very close the the highest order term taken in the Bessel-Fourier series for the IFF solution. Moreover, for ice fragments with average draft $T_{0}\simeq 10$m, and water pressure upon ice blocks in the range of normal atmospheric pressure, we can evaluate from Eqs. (22, 28) and from $\alpha\sim\sqrt{a}\omega_{r}/c$ that $0.25\leq 2\pi\alpha/\omega_{r}=\alpha T\leq 0.45$. This relation together with $L\sim\xi_{1,n}/\alpha$ gives a fundamental general criterion to predict the size (and from here the velocity, energy and vorticity) of such spirals in relation to their period $T$ of revolution. For example a spiral with a period of revolution ten hours can be modeled by our solution with the scaling factor in from of $r$ variable $\alpha\sim 5\times 10^{-4}$m*-1* predicting a spiral/swirl size of $L\sim 30\div 40$Km. Moreover, the parameter $\alpha T$ shown above provides a sort of inverse of characteristic velocity of the spiral.

In order to connect the spiral evolution with possible solitary waves we start again from the quadratic approximation Eqs. (70, 71) and, in the case of weak dissipation such that we can neglect lower order terms containing the product $\omega_{r}\omega_{i}$, and by using the substitutions

[TABLE]

we can transform Eq. (71) in a nonlinear, variable coefficients integro-differential equation

[TABLE]

[TABLE]

where $b_{1,\dots,3}$ are constants of integration. By further substituting $T(\xi)=\ln S(\xi)$ we can further reduce Eq. (81) to a partially separated form

[TABLE]

where $\mathcal{F}_{1,2}$ are functions, and $\mathcal{I}$ is nonlinear integral operator. The first three non-integral terms in Eq. (82) form a known differential form which accepts exact solutions in parametric form (see for example [32] equation 13.6.3.8). Consequently, this modified form obtained for the model equation can be approached by integral equations technique to analyze existence and uniqueness of solutions.

The advantage of such a change of variables through the logarithmic substitution may favor the representation of solutions in terms of traveling localized perturbations along the $\xi$ axis.

7 Nonlinear Stability

The surface of the sea is a random nonlinear wave field, and coherent structures localized in time and space such as wave trains, rogue waves, and wave packages also get birth randomly. There is no correlation between waves at large distances except in cases where wind would blow constantly and uniformly on a large area, but even these events are transient and localized in time. For example, even in the case of the largest wave phase velocity (20 m/s) it would take about half an hour for information to travel from an initial point at the water surface, along 30 Km, time during which the conditions at the initial place have changed already, since the transient seas can change in matter of minutes. Studies on correlations of surface waves back up this conclusion [25]. The case becomes even more complicated when there are other systems in the water which interact with waves and currents, like ice fragments, oil spots, garbage, algae, etc. It is exceptional when a large space and time scale stable structure occurs spontaneously at the surface of the sea, like the spirals observed in the polar seas, [9, 10], Fig. 1. Such spiral structures extend over tens of kilometers, move very slow and survive for days. Such stable, large space-scale, long life time rotating patterns must be generated by a strong type of correlation and nonlinear coupling of modes. The ice and water spirals observed, [9, 10], are definitely the signature of coherent, nonlinear collective interaction, generated by spontaneous large-scale nonlinear collective sea modes.

The linear solutions $\Phi^{1}_{m,n}$ given by Eqs. (28-30) form a complete orthonormal basis of functions for the continuous functions defined on $[r_{0},L]\times[0,2\pi]$, which means that any solution of the fully nonlinear system Eqs. (16-18), or any solutions of some lower order of approximation containing quadratic, cubic, etc. nonlinear terms can be expanded in this basis.

In this section we consider nonlinear terms up to the second order of approximation, which of course includes the linearization Eqs. (19-20), already analyzed in section 3. The system Eqs. (8-10), written in dimensionless form in the second order of approximation with the terms slightly re-ordered, contains now cubic nonlinearities, and has the form

[TABLE]

[TABLE]

[TABLE]

where we denoted

[TABLE]

Before applying a rigorous approach to the nonlinear stability of the exact solution for the linearized system, Eqs. (28-30), we make some qualitative comments. From the $\theta$ and $t$ dependence of the solution Eq. (28), and knowing that the tangent velocity is usually about one order of magnitude less than the radial velocity, we make the hypothesis $v_{\theta}=\omega v_{t}$, with $\omega$ a parameter to be adjusted conveniently. The only term where $v$ occurs in Eqs. (83-85) is through $v_{\theta}/r\sim\omega v_{t}/r$. If we perform this substitution in Eq. (83) the resulting term has a smaller order compared to the rest of the terms in the equation, so it is enough to forward use the linear substitution given by the second equation in Eqs. (20). In this way the sub-system made by the first two equations, Eqs. (28-29), depends only on $\Phi$ and $u$, and once a solution is found for these two quantities, we can solve Eq. (85) for $v$, because this equation transforms now into a linear, first order PDE in $t$ and $r$ which is fully integrable. Under the above hypothesis Eq. (83) becomes

[TABLE]

with

[TABLE]

We note that at large distance from the origin $r\rightarrow\infty$ the terms containing $1/r$ can be neglected in Eq. (86). After some algebra we obtain for this asymptotic limit the condition $(1-a\Phi^{0})u_{t}=0$ which involves the that the radial speed becomes uniform asymptotically. Consequently, from Eq. (84), and in the same limit we obtain that all functions $\Phi,u$ and $v$ approach constant and uniform values in the far asymptotic region.

In the following we plug the linear analytic solutions Eqs. (28-30) in the nonlinear system Eqs. (84-86) containing terms up to order two in $\varepsilon$ and $\delta$. The result is that all linear terms will contain a double series over the order $m$ of the Bessel functions, and over the order $n$ of their zeroes. The quadratic nonlinear terms given be the second parenthesis left hand side of Eqs. (84-86), will involve products of such series. Consequently, the linear solution remain stable with the inclusion of the nonlinear terms if the products of series can be expressed in terms of the solutions series. Namely, if the following relation exists

[TABLE]

[TABLE]

[TABLE]

where $C_{m,n}$ are the coefficients of the solution Eq. (28). Obviously, because of completeness of the Fourier series we have the constraint $m=m_{1}+m_{2}$. This constraint transform the above relation into

[TABLE]

[TABLE]

[TABLE]

The completeness of the Fourier-Bessel series determines the range of the summations over $m_{1}$ and $m$, and then coefficients $C_{m,n}$ should fulfill the resulting nonlinear relation.

In the following we apply Arnold’s convexity method, [42], to prove stability estimates for smooth solutions in our nonlinear two-phases model. This estimation provides $L^{2}$ squared integrable bounds on perturbations of the IFF function and velocities. Our first observation is that the original system Eqs. (1, 2) , in the inviscid case ($\eta=0$), with density and pressure given by Eqs. (5,6) is a planar ideal barotropic fluid, and consequently Hamiltonian. Indeed, the pressure

[TABLE]

depends only on the density and thus warrants the barotropicity. From here we calculate the specific enthalpy $h(\rho)$, and specific energy $e(\rho)$ of our fluid by using [43]

[TABLE]

and obtain

[TABLE]

with $C_{1,2}$ arbitrary constants of integration. Because the flow is planar and barotropic, and we can express the flow vorticity as $\vec{\Omega}=\vec{k}(v_{2,x}-v_{1,y})$, it results the conservation of the quantity $\Omega/\rho$. Consequently, the volume ($\rho dxdy$) double integral on the plane of any smooth function of real variable $F$ depending only on this invariant, can act as a Lagrangian multiplier in sum to the Hamiltonian of the flow

[TABLE]

given here by the sum between the kinetic and internal energy densities. We consider the linear solutions of the linearized system obtained in section 3 for this two-phase, two-dimensional, barotropic flow as equilibrium states. For such equilibrium states the gradient vectors $\nabla[(u)^{2}/2+(v)^{2}/2+h(\rho)]$ and $\nabla(\Omega/\rho)$ can be shown to be orthogonal to the velocity [42]. Consequently, these gradient vectors must be collinear and hence there must be a functional relation between their potentials

[TABLE]

The functional relationship through the function $G$ represents Bernoulli’s law, and $G$ is called the Bernoulli function. This function qualifies for the role of the Lagrangian multiplier density function $F$ from Eq. (91). By applying to the functional in Eq. (91) the nonlinear Lyapunov stability criterion, [42, 44], we obtain the following conditions for the stability against small perturbations of our linear solutions

[TABLE]

and

[TABLE]

The first stability condition Eq. (93) is just a request that the flow must be subsonic, because it actually reduces to $(c^{2}-||\vec{v}^{1}||^{2})/\rho^{1}>0$, and this condition is obviously accomplished in the slow motion of water around the sea-ice. The stability constraint arises from the second stability condition, Eq. (94).

In this equation, it makes sense to write a ratio of gradients, because they are collinear. We calculated the IFF function $\Phi^{1}$, the density, the corresponding velocity components, their vorticity, and their enthalpy $h$ for $m=1$. In Figs. 9, 10 we plotted the functional from Eq. (94) for $m=1$ and for values of the parameters that match the geometrical parameters of the swirls observed in the ocean.

8 Conclusions

In the present work we investigate a two-dimensional compressible Navier-Stokes hydrodynamic model design to explain and study large scale ice swirls formation at the surface of the ocean. The regular motion of the ocean water and sea-ice is powered by winds, ocean currents, thermal convection, Coriolis effect. For the mixture ice and water which resides mainly at the surface, the thermal effects produced by internal convection is not a major cause. The observed swirls ot the surface of polar oceans have very large size yet they move very slow which eliminates the Coriolis force as a dominating cause. Also, the long life time of these structures demonstrates that their formation is rather governed by internal dynamics and ocean currents than waves and wind. The difference in density between ice and water makes this combination a compressible two-dimensional type of fluid mixture. Such two-dimensional structures with rotational symmetry as noticed by aerial photography may exercise three types of large amplitude motions: radial oscillations of compression and dilation, azimuthal motion including rotation and shear flow, and a third one consisting in the coupling of radial modes with azimuthal modes and may thus generate very large scale spiral patterns. The question is if indeed these modes tend to couple. Following day and night cycles of warming and cooling the surface ice may melt and freeze back even while the temperature is almost constant (by transfer of heat through radiation, wind and ocean currents). When melting, the ice trapped in the center induces a radial global two-dimensional compression mode. By the Kelvin theorem, this reduction in radius of the two-dimensional layers, and the corresponding conservation of circulation induces an increasing in azimuthal motion of the fluid. So, compression induces rotation, which is exactly the requested coupling. The circulation-conserving compression rate of $1/r$, while integrated generates exactly a sort of coupling inducing the logarithmic spiral pattern.

The three equations governing mass and momentum conservation form a partial differential system with cubic nonlinearity. A linearization procedure demonstrates the evidence the formation od linearly stable spiral patterns including Archimedean and logarithmic ones. In order to find fully nonlinear solutions we write the density field as a real magnitude and a phase function, and we demonstrate that with a sufficiently good physical approximation the dynamical equation for the phase can be mapped into a sine-Gordon equation whose stable multi-front solutions generate the ice patterns. By truncating the nonlinear system to its quadratic terms we obtain spiral solutions with logarithmic geometry. In addition to the existence of sine-Gordon solitons, the nonlinear solutions are analyzed using two additional mathematical approaches: one predicting the formation of patterns as Townes solitary modes, and another using a series expansion. Pure radial, azimuthal and spiral modes are obtained from the fully nonlinear equations. Combinations of multiple-spiral solutions are also obtained, matching the experimental observations. The nonlinear stability of the spiral patterns is analyzed by Arnold’s convexity method, and the Hamiltonian of the solutions is plotted versus some order parameters showing the existence of geometric phase transitions.

Acknowledgements One of the authors (AL) is grateful to Dalian University of Technology for sharing the preliminary results on this project, and for funding and hospitality during the accomplishment of this research project in 2018-2019.