# Stochastic Transport in Upper Ocean Dynamics II

Bertrand Chapron, Dan Crisan, Darryl Holm, Etienne Mémin, Anna Radomska

► To cite this version:

Bertrand Chapron, Dan Crisan, Darryl Holm, Etienne Mémin, Anna Radomska. Stochastic Transport in Upper Ocean Dynamics II. STUOD 2022 - Workshop on Stochastic Transport in Upper Ocean Dynamics II, 11, Springer Nature Switzerland, pp.XIV, 338, 2024, Mathematics of Planet Earth, 978-3-031-40094-0. 10.1007/978-3-031-40094-0 . hal-04354387

HAL Id: hal-04354387

https://inria.hal.science/hal-04354387v1

Submitted on 19 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Bertrand Chapron · Dan Crisan ·
Darryl Holm · Etienne Mémin ·
Anna Radomska *Editors*

Stochastic
Transport in
Upper Ocean
Dynamics II

STUOD 2022 Workshop, London, UK,
September 26–29

<figure><img src="image_2.png" /></figure>

OPEN ACCESS

<figure><img src="image_3.png" /></figure>

Springer
Mathematics of Planet Earth

Volume 11

**Series Editors**

Dan Crisan, Imperial College London, London, UK

Ken Golden, University of Utah, Salt Lake City, UT, USA

Darryl D. Holm, Imperial College London, London, UK

Mark Lewis, University of Victoria, Victoria, BC, Canada

Yasumasa Nishiura, Tohoku University, Sendai, Miyagi, Japan

Joseph Tribbia, National Center for Atmospheric Research, Boulder, CO, USA

Jorge Passamani Zubelli, Instituto de Matemática Pura e Aplicada, Rio de Janeiro, Brazil
This series provides a variety of well-written books of a variety of levels and styles, highlighting the fundamental role played by mathematics in a huge range of planetary contexts on a global scale. Climate, ecology, sustainability, public health, diseases and epidemics, management of resources and risk analysis are important elements. The mathematical sciences play a key role in these and many other processes relevant to Planet Earth, both as a fundamental discipline and as a key component of cross-disciplinary research. This creates the need, both in education and research, for books that are introductory to and abreast of these developments.

Springer’s MoPE series will provide a variety of such books, including monographs, textbooks, contributed volumes and briefs suitable for users of mathematics, mathematicians doing research in related applications, and students interested in how mathematics interacts with the world around us. The series welcomes submissions on any topic of current relevance to the international Mathematics of Planet Earth effort, and particularly encourages surveys, tutorials and shorter communications in a lively tutorial style, offering a clear exposition of broad appeal.

**Responsible Editors:**

Martin Peters, Heidelberg (martin.peters@springer.com)

Robinson dos Santos, São Paulo (robinson.dossantos@springer.com)

**Additional Editorial Contacts:**

Donna Chernyk, New York (donna.chernyk@springer.com)

Masayuki Nakamura, Tokyo (masayuki.nakamura@springer.com)
Stochastic Transport in
Upper Ocean Dynamics II

STUOD 2022 Workshop, London, UK,
September 26–29

<figure><img src="image_4.png" /></figure>
Editors

Bertrand Chapron
Ifremer - Institut Français de Recherche
pour l'Exploitation de la Mer
Plouzané, France

Darryl Holm
Imperial College London
London, UK

Anna Radomska
Imperial College London
London, UK

Dan Crisan
Imperial College London
London, UK

Etienne Mémin
Campus Universitaire de Beaulieu
Inria - Institut National de Recherche en
Sciences et Technologies du Numérique
Rennes, France

<figure><img src="image_5.png" /></figure>

ISSN 2524-4264
Mathematics of Planet Earth
ISBN 978-3-031-40093-3
https://doi.org/10.1007/978-3-031-40094-0

ISSN 2524-4272 (electronic)
ISBN 978-3-031-40094-0 (eBook)

Mathematics Subject Classification: 60Hxx, 60H17, 70L10, 35R60, 37M05, 37-11, 35Qxx, 65Pxx, 00B25

This work was supported by Horizon 2020 Framework Programme (856408)

© The Editor(s) (if applicable) and The Author(s) 2024. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG
The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

Paper in this product is recyclable.
# Preface

This volume contains the Proceedings of the 3rd Stochastic Transport in Upper Ocean Dynamics Annual Workshop held on 26–29 September 2022. The workshop is part of the project “Stochastic Transport in Upper Ocean Dynamics” (STUOD) supported by an European Research Council Synergy Grant. The four principal investigators of the projects, Prof. Dan Crisan (ICL), Prof. Bertrand Chapron (IFREMER) Prof. Darryl Holm (ICL) and Prof. Etienne Memin (INRIA), were delighted to organize the annual workshop: the event was largest to date with an impressive range of inspiring talks, engaging theoretical and applied sessions, as well as networking opportunities.

The STUOD project is supported by an ERC Synergy Grant, led by three world-class institutions: Imperial College London (ICL), National Institute for Research in Digital Science and Technology (INRIA) and the French Research Institute for Exploitation of the Sea (IFREMER). The project aims to deliver new capabilities for assessing variability and uncertainty in upper ocean dynamics and provide decision makers a means of quantifying the effects of local patterns of sea level rise, heat uptake, carbon storage and change of oxygen content and pH in the ocean. The project will make use of multimodal data and will enhance the scientific understanding of marine debris transport, tracking of oil spills and accumulation of plastic in the sea.

As in the previous years, the 3rd STUOD Annual Workshop 2022 focused on a range of fundamental topical areas, including:

1. Observations at high resolution of upper ocean properties such as temperature, salinity, topography, wind, waves and velocity

2. Large-scale numerical simulations

3. Data-based stochastic equations for upper ocean dynamics that quantify simulation errors

4. Stochastic data assimilation to reduce uncertainty

Each chapter in the present volume illustrates one or several of these topical areas. Many chapters offer new mathematical frameworks that are intended to enhance future research in the STUOD project. The workshop was held in the
hybrid mode and brought together participants from countries such as: UK, France, USA, Italy, Germany, Saudi Arabia, the Netherlands and Switzerland. It was well attended by early-career academics, post-graduate students, senior members of the community and other invited guests.

The scientific program of the four-day conference was divided into five sessions covering Data Assimilation, Physics Models, Data, Numerics for Ocean Models, and Theoretical Analysis of SPDEs. Several members of the STUOD External Advisory Board, Prof Sebastian Reich (University of Potsdam), Prof Baylor Fox-Kemper (Brown University), Prof Rosemary Morrow (Laboratoire d'Études en Géophysique et Océanographie Spatiales), Prof Laurent Debreu (INRIA) and Prof Arnaud Debussche Ecole Normale Supérieure de Rennes, gave invited talks. The programme also included individual presentations by the STUOD principal investigators and post-doctoral researchers and snapshot presentations by invited speakers such as Dr Ali Mashayekhi (Imperial College London) and Dr Hamza Ruzayqat (King Abdullah University of Science and Technology). The workshop also provided opportunities for investigators at an early stage of their career to have discussions with established scientists, fostering potential future research collaborations, networking as well as inclusion and training of the next generation of researchers.

<figure><img src="image_6.png" /><figcaption>Some of the in-person participants attending the event</figcaption></figure>
The following is a brief description of the 15 contributions included in the proceedings:

The submitted manuscripts include the chapter by Adrien Bella, Noé Lahaye, and Gilles Tissot entitled "Internal Tides Energy Transfers and Interactions with the Mesoscale Circulation in Two Contrasted Areas of the North Atlantic". The chapter investigates the energy budget of the internal tide and its life cycle with a high-resolution numerical simulation and a vertical normal mode decomposition. Two areas of interest are considered: the Azores Islands over the mid-Atlantic ridge and the Gulf Stream offshore the North of the US East coast shelf break. Low mode (1 and 2) internal tides are found to propagate from 100 km (mode 2) to 1000 km (mode 1) away from their generation sites. Waves lose a significant portion of their energy as they propagate through the Gulf Stream, in contrast to the Azores domain. In the Gulf Stream domain, the mesoscale circulation is responsible for energy transfers from low to high mode internal tides, while the topographic scattering is dominant in the Azores area. This transfer of energy toward high modes favours energy dissipation. The mesoscale is significant in the energy budget of modes higher than mode 1 for both domains, and for all baroclinic modes in the Gulf Stream area. The internal tide is found to extract or lose energy toward the mesoscale circulation, but this accounts for less than 14% of the energy scattered from low internal tide modes to higher ones once summed over all contributions of the modal energy budget.

The contribution of Paolo Cifani, Sagy Ephrati and Milo Viviani entitled "Sparse-Stochastic Model Reduction for 2D Euler Equations" introduces a reduction technique for ideal fluids.

The 2D Euler equations are a simple but rich set of non-linear PDEs that describe the evolution of an ideal inviscid fluid, for which one dimension is negligible. Solving numerically these equations can be extremely demanding. Several techniques to obtain fast and accurate simulations have been developed during the last decades. In this chapter, the authors present a novel approach which combines recent developments in the stochastic model reduction and conservative semi-discretization of the Euler equations. In particular, starting from the Zeitlin model on the 2-sphere, they derive reduced dynamics for large scales and we close the equations either deterministically or with a suitable stochastic term. Numerical experiments show that, after an initial turbulent regime, the influence of small scales to large scales is negligible, even though a non-zero transfer of energy among different modes is present.

Franco Flandoli, Silvia Morlacchi, and Andrea Papini present in their work, "Effect of Transport Noise on Kelvin-Helmholtz Instability", a numerical investigation of the dissipation properties of very small-scale transport noise. As a test problem, the authors consider the Kelvin-Helmholtz instability and compare the inviscid case, the viscous one, both without noise, and the inviscid case perturbed by transport noise. They observe a partial similarity with the viscous case, namely a delay of instability.

The chapter by Daniel Goodair and Dan Crisan, entitled "On the 3D Navier-Stokes Equations with Stochastic Lie Transport", investigates the existence and uniqueness of maximal solutions to the 3D SALT Navier-Stokes Equation in
velocity and vorticity form, on the torus and the bounded domain, respectively. The work partners an earlier paper of the authors as an application of the abstract framework presented there. They demonstrate the efficacy of the earlier work in showing well-posedness for both the velocity and vorticity form of the equation, as well as obtaining the first analytically strong existence result for a fluid equation perturbed by Lie transport noise on a bounded domain.

Darryl Holm, Ruiao Hu, and Oliver Street in their paper entitled "On the Interactions Between Mean Flows and Inertial Gravity Waves" derive a generalized Lagrangian mean (GLM) theory as a phase-averaged Eulerian Hamilton variational principle expressed as a composition of two smooth invertible maps. Following an earlier work of Holm, they consider 3D inertial gravity waves (IGWs) in the Euler-Boussinesq fluid approximation. The authors provide both deterministic and stochastic closure models for GLM IGWs at leading order in 3D complex vector WKB wave asymptotics. The chapter brings the earlier results of Holm at leading order into an easily assimilable short form and proposes a stochastic generalization of the wave mean flow interaction (WMFI) equations for IGWs.

The contribution of Quentin Jamet, Etienne Mémin, Franck Dumas, Long Li, and Pierre Garreau entitled "Toward a Stochastic Parameterization for Oceanic Deep Convection" investigates parametrizations for ocean convections. Current climate models are known to systematically overestimate the rate of deep-water formation at high latitudes in response to too deep and too frequent deep convection events. The authors propose in this study to investigate a misrepresentation of deep convection in Hydrostatic Primitive Equation (HPE) ocean and climate models due to the lack of constraints on vertical dynamics. They discuss the potential of the Location Uncertainty (LU) stochastic representation of geophysical flow dynamics to help in the process of re-introducing some degree of non-hydrostatic physics in HPE models through a pressure correction method. The authors then test these ideas with idealized Large Eddy Simulations (LES) of buoyancy-driven free convection with the CROCO modelling platform. Preliminary results at LES resolution exhibit a solution obtained with our Quasi-nonhydrostatic (Q-NH) model that tends toward the reference non-hydrostatic (NH) model. As compared to a pure hydrostatic setting, our Q-NH solution exhibits vertical convective plumes with a larger horizontal structure, a better spatial organization, and a reduced intensity of their associated vertical velocities. The simulated Mixed Layer Depth (MLD) deepening rate is however too slow in our Q-NH experiment as compared to the reference NH, a behaviour that opposes to that of hydrostatic experiments of producing too fast MLD deepening rate. These preliminary results are encouraging, and support future efforts in the direction of enriching coarse resolution, hydrostatic ocean, and climate models with a stochastic representation of non-hydrostatic physics.

The contribution of Alexander Lobbe, Dan Crisan, Darryl Holm, Etienne Mémin, Oana Lang, and Bertrand Chapron, entitled "Comparison of Stochastic Parametrization Schemes Using Data Assimilation on Triad Models", investigates stochastic parametrization schemes. In recent years, stochastic parametrizations have been ubiquitous in modelling uncertainty in fluid dynamics models. One source
of model uncertainty comes from the coarse-graining of the fine-scale data and is
in common usage in computational simulations at coarser scales. In this chapter,
the authors look at two such stochastic parametrizations: the Stochastic Advection
by Lie Transport (SALT) parametrization and the Location Uncertainty (LU)
parametrization. Whilst both parametrizations are available for full-scale models,
the authors study their reduced order versions obtained by projecting them on a
complex vector Fourier mode triad of eigenfunctions of the curl. Remarkably, these
two parametrizations lead to the same reduced order model, termed the helicity-
preserving stochastic triad (HST). This reduced order model is then compared with
an alternative model which preserves the energy of the system, and which is termed
the energy preserving stochastic triad (EST). These low-dimensional models are
ideal benchmark models for testing new Data Assimilation algorithms: they are easy
to implement, exhibit diverse behaviours depending on the choice of the coefficients,
and come with natural physical properties such as the conservation of energy and
helicity.

Erwin Luesink and Bernard Geurts consider in the chapter “An Explicit Method to Determine Casimirs in 2D Geophysical Flows” a new method to construct Casimirs for geophysical flows.

Conserved quantities in geophysical flows play an important role in the charac-
terization of geophysical dynamics and aid the development of structure-preserving
numerical methods. A significant family of conserved quantities is formed by
the Casimirs: These are integral conservation laws that are in the kernel of the
underlying Poisson bracket. The Casimirs hence determine the geometric structure
of the geophysical fluid equations among which the enstrophy is well known. Often
Casimirs are proposed on heuristic grounds and later verified to be part of the kernel
of the Poisson bracket. In this work, the authors explicitly construct Casimirs by
rewriting the Poisson bracket in vorticity-divergence coordinates thereby providing
explicit construction of Casimirs for 2D geophysical flow dynamics.

The work of Igor Maingonnat, Gilles Tissot, and Noé Lahaye is entitled
“Correlated Structures in a Balanced Motion Interacting with an Internal Wave”
and investigates the correlations between a balanced motion and the incoherent part
of a wave in an idealized configuration.

Characterizing the loss of coherence of an internal tide propagating through
mesoscale turbulence has been a major challenge in oceanography, particularly due
to its implications for the interpretation of satellite data. The authors introduce a new
modal decomposition technique, named broad-band proper orthogonal decomposi-
tion (BBPOD), which consists in performing a proper orthogonal decomposition
(POD) on complex demodulated variables. After connecting BBPOD to the standard
SPOD, they show that BBPOD, coupled with the extended POD technique, enables
them to associate the principal components of the incoherent field to the slow flow
structures responsible for this loss of coherence through triadic interactions with the
incident wave.

The chapter by Étienne Mémin, Long Li, Noé Lahaye, Gilles Tissot, and Bertrand
Chapron, entitled “Linear Wave Solutions of a Stochastic Shallow Water Model”
investigates the wave solutions of a stochastic rotating shallow water model. The
approximate model covered in the chapter provides an interesting simple description of the interplay between waves and random forcing ensuing either from the wind or coming as the feedback of the ocean on the atmosphere and leading in a very fast way to the selection of some wavelength. This interwoven, yet simple, mechanism explains the emergence of typical wavelength associated with near-inertial waves. Ensemble-mean waves that are not in phase with the random forcing are damped at an exponential rate, whose magnitude depends on the random forcing variance. Geostrophic adjustment is also interpreted as a statistical homogenization process in which, in order to conserve potential vorticity, the small-scale component tends to align to the velocity fields to form a statistically homogeneous random field.

The contribution of Said Ouala, Bertrand Chapron, Fabrice Collard, Lucile Gaultier, and Ronan Fablet entitled "Analysis of Sea Surface Temperature variability Using Machine Learning" analyses sea surface temperature (SST). SST is a critical factor in the global climate system and plays a key role in many marine processes. Understanding the variability of SST is therefore important for a range of applications, including weather and climate prediction, ocean circulation modelling, and marine resource management. In this study, the authors use machine learning techniques to analyse SST anomaly (SSTA) data from the Mediterranean Sea over a period of 33 years. The objective is to best explain the temporal variability of the SSTA extremes. These extremes are revealed to be well explained through a non-linear interaction between multi-scale processes. The results contribute to better unveil factors influencing SSTA extremes, and the development of more accurate prediction models.

The contribution of Sebastian Reich entitled "Data Assimilation: A Dynamic Homotopy-Based Coupling Approach" covers a new approach for data assimilation based on homotopy. Homotopy approaches to Bayesian inference have found widespread use especially if the Kullback-Leibler divergence between the prior and the posterior distribution is large. The author extends one of these homotopy approach to include an underlying stochastic diffusion process. The underlying mathematical problem is closely related to the Schrödinger bridge problem for given marginal distributions. He shows that the proposed homotopy approach provides a computationally tractable approximation to the underlying bridge problem. In particular, the implementation builds upon the widely used ensemble Kalman filter methodology and extends it to Schrödinger bridge problems within the context of sequential data assimilation.

Valentin Resseguier, Yicun Zhen, and Bertrand Chapron show in their work entitled "Constrained Random Diffeomorphisms for Data Assimilation" that both the Stochastic Advection by Lie Transport (SALT) and the Location Uncertainty (LU) equations can be recovered using a prescribed definition of a random diffeomorphism used to perturb the physical space. However, unlike the SALT and LU settings, they propose a perturbation scheme does not directly rely on a particular physics. Hence, the random mapping is more flexible and can be applied to any PDE.

The work of Gilles Tissot, Étienne Mémin, and Quentin Jamet entitled "Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty" provides
a stochastic version under location uncertainty of the compressible Navier-Stokes equations. To that end, some clarifications of the stochastic Reynolds transport theorem are given when stochastic source terms are present in the right-hand side. The authors apply this conservation theorem to density, momentum, and total energy in order to obtain a transport equation of the primitive variables, i.e. density, velocity, and temperature. They show that performing low Mach and Boussinesq approximations to this more general set of equations allows to recover the known incompressible stochastic Navier-Stokes equations and the stochastic Boussinesq equations, respectively. Finally, they provide some research directions for using this general set of equations in the perspective of relaxing the Boussinesq and hydrostatic assumptions for ocean modelling.

Francesco Tucciarone, Étienne Mémin, and Long Li present in their work entitled "Data-Driven Stochastic Primitive Equations with Dynamic Modes Decomposition" an implementation of a stochastic version of the primitive equations into the NEMO community ocean model to assess the capability of the so-called Location Uncertainty framework in representing the small scales of the ocean flows. The work is important as planetary flows are characterized by interaction of phenomena in a huge range of scales, and it is unaffordable today to resolve numerically the complete ocean dynamics.

Finally, the STUOD Organizing Committee would like to acknowledge the financial and in-kind support received from several sources: the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (ERC, Grant Agreement No 856408)—for providing funds to cover the travel expenses of the invited speakers, catering costs, and administrative support; Imperial College London—for offering the conference venue.

Plouzané, France
London, UK
London, UK
Rennes, France
London, UK

May 2023

Bertrand Chapron

Dan Crisan

Darryl Holm

Etienne Mémin

Anna Radomska
# Contents

<table><tbody><tr><td><b>Internal Tides Energy Transfers and Interactions with the Mesoscale Circulation in Two Contrasted Areas of the North Atlantic</b><br>Adrien Bella, Noé Lahaye, and Gilles Tissot</td><td>1</td></tr><tr><td><b>Sparse-Stochastic Model Reduction for 2D Euler Equations</b><br>Paolo Cifani, Sagy Ephrati, and Milo Viviani</td><td>17</td></tr><tr><td><b>Effect of Transport Noise on Kelvin-Helmholtz Instability</b><br>Franco Flandoli, Silvia Morlacchi, and Andrea Papini</td><td>29</td></tr><tr><td><b>On the 3D Navier-Stokes Equations with Stochastic Lie Transport</b><br>Daniel Goodair and Dan Crisan</td><td>53</td></tr><tr><td><b>On the Interactions Between Mean Flows and Inertial Gravity Waves in the WKB Approximation</b><br>Darryl D. Holm, Ruiao Hu, and Oliver D. Street</td><td>111</td></tr><tr><td><b>Toward a Stochastic Parameterization for Oceanic Deep Convection</b><br>Quentin Jamet, Etienne Mémin, Franck Dumas, Long Li, and Pierre Garreau</td><td>143</td></tr><tr><td><b>Comparison of Stochastic Parametrization Schemes Using Data Assimilation on Triad Models</b><br>Alexander Lobbe, Dan Crisan, Darryl Holm, Etienne Mémin, Oana Lang, and Bertrand Chapron</td><td>159</td></tr><tr><td><b>An Explicit Method to Determine Casimirs in 2D Geophysical Flows</b><br>Erwin Luesink and Bernard Geurts</td><td>193</td></tr><tr><td><b>Correlated Structures in a Balanced Motion Interacting with an Internal Wave</b><br>Igor Maingonnat, Gilles Tissot, and Noé Lahaye</td><td>207</td></tr></tbody></table>
<table><tr><td>Linear Wave Solutions of a Stochastic Shallow Water Model</td><td>223</td></tr><tr><td>Etienne Mémin, Long Li, Noé Lahaye, Gilles Tissot, and Bertrand Chapron</td><td></td></tr><tr><td>Analysis of Sea Surface Temperature Variability Using Machine Learning</td><td>247</td></tr><tr><td>Said Ouala, Bertrand Chapron, Fabrice Collard, Lucile Gaultier, and Ronan Fablet</td><td></td></tr><tr><td>Data Assimilation: A Dynamic Homotopy-Based Coupling Approach</td><td>261</td></tr><tr><td>Sebastian Reich</td><td></td></tr><tr><td>Constrained Random Diffeomorphisms for Data Assimilation</td><td>281</td></tr><tr><td>Valentin Resseguier, Yicun Zhen, and Bertrand Chapron</td><td></td></tr><tr><td>Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty</td><td>293</td></tr><tr><td>Gilles Tissot, Étienne Mémin, and Quentin Jamet</td><td></td></tr><tr><td>Data Driven Stochastic Primitive Equations with Dynamic Modes Decomposition</td><td>321</td></tr><tr><td>Francesco L. Tucciarone, Etienne Mémin, and Long Li</td><td></td></tr><tr><td>Index</td><td>337</td></tr></table>
# Internal Tides Energy Transfers and Interactions with the Mesoscale Circulation in Two Contrasted Areas of the North Atlantic

<figure><img src="image_7.png" /></figure>

Check for
updates

Adrien Bella, Noé Lahaye, and Gilles Tissot

**Abstract** The energy budget of the internal tide and its life cycle is investigated with a high resolution numerical simulation and a vertical normal mode decomposition. Two areas of interest are considered: the Azores Islands over the mid Atlantic ridge and the Gulf Stream offshore the North of the US East coast shelf break. Low mode (1 and 2) internal tides are found to propagate from 100 km (mode 2) to 1000 km (mode 1) away from their generation sites. Waves lose a significant portion of their energy as they propagate through the Gulf Stream, in contrast to the Azores domain. In the Gulf Stream domain, the mesoscale circulation is responsible for energy transfers from low to high modes internal tides, while the topographic scattering is dominant in the Azores area. This transfer of energy toward high modes favours energy dissipation. The mesoscale is significant in the energy budget of modes higher than mode 1 for both domains, and for all baroclinic modes in the Gulf Stream area. The internal tide is found to extract or lose energy toward the mesoscale circulation, but this accounts for less than 14%, of the energy scattered from low internal tide modes to higher ones once summed over all contributions of the modal energy budget.

**Keywords** Internal tide · Mesoscale flow · Energy budget · Vertical modal decomposition · High resolution numerical simulation

## 1 Introduction

Internal tides are a category of oceanic internal inertia-gravity waves that are generated when the astronomical tide interacts with topographic features such as shelf breaks, ridges or seamounts. They are encountered in many areas in the ocean [Zaron et al. 2022, Zhao et al. 2016, among others] and are of crucial importance for

<footnote>A. Bella (✉) · N. Lahaye · G. Tissot
INRIA Rennes Bretagne Atlantique and IRMAR – UMR CNRS 6625, Rennes, France
e-mail: adrien.bella@inria.fr</footnote>
the mixing of deep water and the closure of the Meridional Overturning circulation [Munk and Wunsch 1998]. The role of mesoscale currents for the energy pathway starting from the barotropic tide input and leading to this mixing is not yet fully understood. As a different thread of motivation, internal tides exhibit the same spatial scale as (sub-)mesoscale features, which renders challenging the estimation of the geostrophic velocities and the disentanglement of the two types of motions from satellite altimetry [Ponte and Klein 2015]. This problem is exacerbated by the ability of the new SWOT mission to obtain a resolution of 20 km to 30 km, compared to the approximately 100 km of current altimeter products [Ballarotta et al. 2019]. It is therefore needed to better understand the dynamics of the internal tide if we hope to disentangle it from the geostrophic field via data assimilation for instance.

Since a decade or so [e.g. Arbic et al. 2010], global and regional realistic simulations have reached a sufficient resolution to be able to represent both internal tides and mesoscale features, opening the way to studying their interactions. Most studies focusing on these interactions use the vertical modal decomposition framework [Gill 1982], which allows to get a clear separation of the barotropic and baroclinic tides as well as alleviating the computational cost associated with processing 3D simulation outputs.

The importance of the mesoscale circulation and buoyancy field have been shown by Kelly and Lermusiaux [2016] in a realistic configuration over the Gulf Stream area for the energy budget of the first baroclinic mode (denoted mode 1). They quantified the energy transfers between the mode 1 and the mesoscale circulation as well as its relative importance compared to topographic effects, showing that the mesoscale explains first order interactions pattern visible on the mode 1 energy flux divergence. Similarly, Pan et al. [2021] have analysed the energy budget of the mode 1 in a realistic setup including the mesoscale circulation and heterogeneous buoyancy field in the Gulf Stream area and the Mallau Island, leading to the same conclusion concerning the importance of mesoscale—internal tides interactions. Both studies have also considered realistic simulations with baroclinic currents, and have used a modal decomposition approach to define their energy budget. The deviation of internal tide energy flux have been analysed by Duda et al. [2018], indicating a refracting behaviour of a Gulf Stream like current on beams of internal tides energy fluxes. Dunphy and Lamb [2014] have shown in an idealised setup that the interactions of the first baroclinic internal tide mode with the barotropic eddy and the first baroclinic eddy mode lead to a modification of the energy flux propagation and a scattering of energy toward higher modes, respectively.

In the present paper, our focus will be on the energy budget yielded by a modal analysis and we will in particular extend the work of other studies to the 10 first baroclinic modes, as well as studying the couplings between modes.

The contribution of the mesoscale flow and buoyancy field will be taken into account and the relative importance of the coupling between modes as well as their spatial pattern will be quantified in an area where topographic features are prominent, and another area with a strong mesoscale circulation.

The document is organised as follows: Sect. 2 will explain the theoretical framework with the modal decomposition and the governing dynamical equations.
Section 3 will detail the dataset used and computational techniques. Section 4 will present results. Section 5 closes this article with a conclusion.

## 2 Governing Equations and Energy Budget

To investigate the modal internal tide energy budget, we rely on linearized equations of motions projected on vertical normal modes. The derivation follows previous studies (e.g. Kelly and Lermusiaux [2016]) and is only briefly depicted here. It starts from the hydrostatic primitive equations under the Boussinesq approximation with a free surface [e.g. Vallis 2017]:

$$ \partial_t \mathbf{u}_h + \mathbf{u} \cdot \nabla \mathbf{u}_h + f \vec{e}_z \wedge \mathbf{u}_h = -\nabla_h p - \nabla \Pi_{\text{tide}}, \quad (1a) $$

$$ \partial_z p - b = 0, \quad (1b) $$

$$ \partial_t b + \mathbf{u} \cdot \nabla b = 0, \quad (1c) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (1d) $$

where $\mathbf{u}$ is the velocity, $b$ the buoyancy, $p$ the reduced pressure (divided by $\rho_0$, the reference density), and the index $h$ denotes horizontal component of a 3D vector. Other variable names follow standard notations. Here, the forcing and dissipation terms have been omitted on purpose, except for the tidal potential term $\Pi_{\text{tide}}$.

Next, the flow is decomposed into a low-frequency mesoscale flow and a high-frequency component that includes the internal tide: i.e., for the velocity, we split $\mathbf{u} = \mathbf{U} + \tilde{\mathbf{u}}$. We will further decompose the low frequency part of the buoyancy into a long-time mean and a slowly variable part: $\mathbf{B} = \bar{\mathbf{B}} + \mathbf{B}'$ ($\bar{\mathbf{B}}$ denotes long time average), and introduce the Brunt-Väisälä frequency: $\mathbf{N}^2 = \partial_z \mathbf{B}$ and $\bar{\mathbf{N}}^2 = \partial_z \bar{\mathbf{B}}$. Subtracting the low frequency equations (1) to the initial system and further assuming that nonlinearities amongst the high-frequency flow are negligible, one obtains the following linear system of equations:

$$ \partial_t \tilde{\mathbf{u}} + \mathbf{U} \cdot \nabla \tilde{\mathbf{u}} + \tilde{\mathbf{u}} \cdot \nabla \mathbf{U} + f \vec{e}_z \wedge \tilde{\mathbf{u}} = -\nabla \tilde{\mathbf{p}} - \nabla \Pi_{\text{tide}}, \quad (2a) $$

$$ \partial_z \tilde{\mathbf{p}} - \tilde{\mathbf{b}} = 0, \quad (2b) $$

$$ \partial_t \tilde{\mathbf{b}} + \mathbf{U} \cdot \nabla \tilde{\mathbf{b}} + \tilde{\mathbf{u}}_h \cdot \nabla_h (\bar{\mathbf{B}} + \mathbf{B}') + \tilde{\mathbf{w}} (\bar{\mathbf{N}}^2 + \mathbf{N}'^2) = 0, \quad (2c) $$

$$ \nabla \cdot \tilde{\mathbf{u}} = 0. \quad (2d) $$

Note that, since the above system of equations is linear in the high-frequency variables, and assuming that most of the mesoscale variability has a timescale longer or equal to the cutoff period used in the internal tides complex demodulation (3 days), band-pass filtering of these equations leaves the system unchanged. Hence, in our study, we will focus on the semi diurnal (including M2 and S2 frequencies)
internal tide with a period centred on 12.2 h, relying on the above set of equations where the high-frequency variables describe this frequency band. These equations are completed by the usual linearized boundary conditions at the mean surface $\bar{\eta}$ and the bottom at $z = -H$:

$$ \tilde{p}(\bar{\eta}) = g\tilde{\eta}, \qquad (3) $$

$$ \tilde{w}(\bar{\eta}) = \partial_t(\bar{\eta}) + \mathbf{U}_h(\bar{\eta}) \cdot \nabla_h(\bar{\eta}), \qquad (4) $$

$$ \tilde{w} = -\tilde{u}_h(-H) \cdot \nabla H. \qquad (5) $$

We thus have the internal tide part of the primitives equations that are now ready to be projected on a set of vertical modes. The latter are given locally by the standard Sturm-Liouville problem for internal waves in a rotating stratified ocean [e.g. Gill 1982, Vallis 2017], assuming a flat ocean, a horizontally-homogeneous stratification profile and no background flow:

$$ \left( \frac{\Phi'_n}{\bar{N}^2} \right)' + \frac{\Phi_n}{c_n^2} = 0, \quad \text{with } \Phi'_n = 0 \text{ at } z = -H, \quad \text{and } g\Phi'_n + \bar{N}^2\Phi_n = 0 \text{ at } z = \bar{\eta}, \qquad (6) $$

where prime denotes vertical derivative. The obtained set of vertical modes $\Phi_n$ is complemented by another set of vertical functions which obey:

$$ \varphi'_n = \Phi_n, \quad \Phi_n' = -\frac{\bar{N}^2}{c_n^2} \varphi_n, \qquad (7) $$

and these modes follow the orthogonality conditions:

$$ \int_{-H}^{\bar{\eta}} \Phi_m \Phi_n \, dz = \int_{-H}^{\bar{\eta}} \frac{\bar{N}^2}{c_n^2} \varphi_m \varphi_n \, dz + \frac{g}{c_n^2} \varphi_m(\bar{\eta}) \varphi_n(\bar{\eta}) = H \delta_{mn}. \qquad (8) $$

Physical fields can then be expanded over these bases through:

$$ [\mathbf{u}_h, p] = \sum_n [\mathbf{u}_n, p_n] \Phi_n, \qquad (9) $$

$$ [w, b] = \sum_n [w_n, \bar{N}^2 b_n] \varphi_n, \qquad (10) $$

and the modal projection coefficients are obtained as follows:

$$ [\mathbf{u}_n, p_n] = \langle [\mathbf{u}_h, p], \Phi_n \rangle, \qquad (11) $$

$$ w_n = \frac{1}{c_n^2} \left( \langle \varphi_n, w \bar{N}^2 \rangle + \frac{g}{H} w(\bar{\eta}) \varphi_n(\bar{\eta}) \right), \qquad (12) $$
$$b_n = \frac{1}{c_n^2} \left( \langle \varphi_n, b \rangle + \frac{g}{H} \frac{b(\bar{\eta})}{\bar{N}^2(\bar{\eta})} \varphi_n(\bar{\eta}) \right), \quad (13)$$

with the inner product defined as such: $\langle f, g \rangle = \frac{1}{H} \int_{-H}^{\bar{\eta}} f(z)g(z) dz$. Note that the two bases have a parametric dependency on the horizontal coordinates because of the varying topography $H(x, y)$ and stratification profile $\bar{N}^2(z; x, y)$. In this paper, we are going to consider the modes 0 to 10, with the mode 0 being the barotropic tide, and all others being the baroclinic internal tide.

The dynamical equations are then obtained by projecting the tidal momentum equations and continuity on a mode $\Phi_m$, and buoyancy equations on the mode $\varphi_m$, using the relations exposed above plus integration by part and Leibniz formula, and substituting the buoyancy modal amplitude by the pressure anomaly (after projection of the hydrostatic balance on $\varphi_m$) [Kelly and Lermusiaux 2016, Kelly 2016, Pan et al. 2021]. Omitting the ~ notation from now on, one obtains—after a tedious but straightforward derivation:

$$\begin{aligned} \partial_t \boldsymbol{u}_m + \nabla_h p_m + f \ddot{\boldsymbol{e}}_z \wedge \boldsymbol{u}_m = & - \sum_n U_{mn} \cdot \nabla \boldsymbol{u}_n - \sum_n \boldsymbol{u}_n U_{mn}^{\Phi} \\ & - \sum_n \boldsymbol{u}_n \cdot U_{mn}^{\nabla} - \sum_n w_n U_{mn}^{z} \\ & - \sum_n p_n T_{mn} - \frac{1}{H} \nabla \Pi_{\text{tide}}(\varphi_m(\bar{\eta}) - \varphi_m(-H)), \end{aligned} \quad (14a)$$

$$\begin{aligned} \partial_t p_m - c_m^2 w_m = & -\frac{g}{H} \varphi_m(\bar{\eta}) U_h(\bar{\eta}) \cdot \nabla \eta - \sum_n U_{mn}^p \cdot \nabla p_n \\ & - \sum_n p_n \left\langle \varphi_m, U_h \cdot \nabla \left( \frac{\bar{N}^2}{c_n^2} \varphi_m \right) \right\rangle \\ & + \sum_n \boldsymbol{u}_n \cdot (B_{mn} - \bar{B}_{mn}) + \sum_n w_n \left\langle \varphi_m, \varphi_n N^2 \right\rangle, \end{aligned} \quad (14b)$$

$$\nabla_h \cdot (H \boldsymbol{u}_m) + H w_m = H \sum_n \boldsymbol{u}_n T_{mn}, \quad (14c)$$

with:

$$T_{mn} = \langle \Phi_m, \nabla_h \Phi_n \rangle, \quad U_{mn} = \langle \Phi_m, U_h \Phi_n \rangle, \quad U_{mn}^{\Phi} = \langle \Phi_m, U_h \cdot \nabla_h \Phi_n \rangle,$$

$$U_{mn}^{\nabla} = \langle \Phi_m, \Phi_n \nabla U_h \rangle, \quad U_{mn}^{z} = \langle \Phi_m, \varphi_n \partial_z U_h \rangle, \quad U_{mn}^{p} = \left\langle \varphi_m, U_{tot} \frac{\bar{N}^2}{c_n^2} \varphi_n \right\rangle,$$

$$B_{mn} = \langle \varphi_m, \Phi_n \nabla B \rangle, \quad \bar{B}_{mn} = \langle \varphi_m, \Phi_n \nabla_h \bar{B} \rangle.$$
We then take the scalar product of the projected momentum equation (14a) with $\mathbf{u}_n$, multiply the pressure equation (14b) with $p_m$ and substitute $w_n$ using the continuity equation (14c) to obtain the energy budget of the internal tide mode $m$ in presence of a slowly varying mesoscale circulation and buoyancy field:

$$
\begin{aligned}
& \overbrace{\frac{\partial}{\partial t} \left( \frac{\mathbf{u}_m^2}{2} + \frac{p_m^2}{2c_m^2} \right)}^{\text{temporal variation of energy}} + \overbrace{\frac{1}{H} \nabla_h \cdot (H \mathbf{u}_m p_m)}^{\text{divergence of energy flux}} = \\
& \quad - \sum_n \left[ \overbrace{\underbrace{\mathbf{u}_m \cdot (U_{mn} \cdot \nabla \mathbf{u}_n) - \mathbf{u}_m \cdot \mathbf{u}_n U_{mn}^{\Phi}}_{\text{background flow shear}} + \frac{p_m}{c_m^2} U_{mn}^p \cdot \nabla p_n} \right. \\
& \qquad \left. - \underbrace{\mathbf{u}_n \cdot U_{mn}^{\nabla} \cdot \mathbf{u}_m + \mathbf{u}_m \cdot U_{mn}^{z}}_{\text{horizontal gradient of background buoyancy}} w_n \right. \\
& \qquad + \overbrace{\frac{p_m}{H c_m^2} \mathbf{u}_n \cdot (\tilde{B}_{mn} + B_{mn})}^{\text{Topography and mean stratification: } C_{mn}} + \overbrace{\frac{p_m}{c_m^2} w_n \langle \varphi_m, N^2 \varphi_n \rangle}^{\text{slowly variable stratification}} \\
& \qquad + \left. \underbrace{p_m \mathbf{u}_n \cdot T_{nm} - p_n \mathbf{u}_m \cdot T_{mn}}_{\text{Tidal potential work}} \right] + R_m
\end{aligned}
\quad (15)
$$

with

$$
\begin{aligned}
R_m = & - \overbrace{\frac{p_m}{H c_m^2} g \varphi_m(\bar{\eta}) \mathbf{U}_h(\bar{\eta}) \cdot \nabla \tilde{\eta}}^{\text{advection of ssh by the background flow}} \\
& - \sum_n \overbrace{\frac{p_m}{c_m^2} \sum_n p_n \left\langle \varphi_m, \mathbf{U}_h \cdot \nabla_h \left( \frac{\tilde{N}^2}{c_n^2} \varphi_n \right) \right\rangle}^{\text{advection of mean stratification by the background flow}} \\
& - \overbrace{\mathbf{u}_m \cdot \nabla \Pi_{\text{tide}}(\varphi_m(\bar{\eta}) - \varphi_m(-H))}^{\text{Tidal potential work}}.
\end{aligned}
\quad (16)
$$

By anticipation of the results in Fig. 3 in Sect. 4.2, the term $R_m$ gathers all contributions that are not of primary importance (and will not be discussed in this paper).

This equation describes the different energy exchanges (coupling) between modes as well as energy sources and sinks. Most coupling terms can be represented with coupling matrices –noted $K_{mn}$ here—that will be shown below in Sect. 4.2 (Fig. 2). Generally, this coupling matrix contains an anti-symmetric part, describing net exchanges of energy between modes, and a symmetric part, with the diagonal included, that describes gain/loss of energy for the internal tide.
The first three terms in the right-hand side (RHS) arise from the presence of a mesoscale flow, while the next two ones include contributions from associated buoyancy perturbations as well as seasonal variations. The term labelled as slowly variable stratification arises during the derivation of the CSW buoyancy equation (14) because of our choice to construct the modal bases by using a stationary profile of stratification in the problem (6). In case of a definition of the vertical modal bases with the instantaneous stratification profile, this term would vanish. Indeed, for a linear propagation of an internal wave without topography, mesoscale and dissipation, the vertical modes diagonalise the operator, thus decouples the components from each other. By time-averaging the stratification profile, the associated vertical mode basis becomes inadequate to perform this decoupling. The slowly variable stratification term then arises, taking into account in the energy equation coupling between modes associated with this linear propagation. This term can be interpreted as energy exchanges between modes during a linear propagation of an internal wave without topography, mesoscale and dissipation. Since vertical modes are eigenfunctions of this associated idealised operator, we associate the variable stratification term to the inadequacy of the basis to represent internal wave dynamics (through the time-averaged stratification profile). We do not expect, however, that using a basis defined with a slowly time-varying stratification would change significantly the magnitude of the dominant interaction terms discussed in this paper (see e.g. Sect. 4.2), and in particular the mesoscale and topographic contributions.

Finally, $C_{mn} = -C_{nm}$ is the topographic scattering, which includes the generation of baroclinic tide from the barotropic tide ($C_{m0}$). Note that this term has no symmetric part—it can only redistribute energy amongst the vertical modes.

Once averaged in time over a sufficiently long period (typically, one month), the time-variation of the energy vanishes and the above modal energy equation reduces to the following, simplified form:

$$ \nabla_h \mathbf{F}_m = - \sum_n K_{mn} + R_m, \qquad (17) $$

where $\mathbf{F}_m$ denotes the modal energy flux, $K_{mn}$ gathers the energy transfers, sources and sinks that are not negligible and $R_m$ gathers all the neglected terms. These include a contribution from the free surface and the tidal potential, which would both vanish under a rigid-lid assumption for the baroclinic modes—in agreement with the fact that we find they have negligible magnitude in our diagnostics (cf. Sect. 4). The remaining term in $R_m$ involves horizontal variations of the modes and associated background stratification, which is small since a year average is used.
# 3 Data and Method

## 3.1 eNATL60 Simulation

The high-resolution realistic simulation of the North Atlantic Ocean eNATL60 [Brodeau et al. 2020] is used to diagnose the energy budget of the internal tide following Eq. (15). It uses the Nucleus for European Modelling of the Ocean (NEMO) model, which solves the primitive equations under the Boussinesq and hydrostatic approximations with an Arakawa C grid and z-coordinate with partial step [Madec and the NEMO team 2008]. The eNATL60 run includes astronomical tidal forcing with the M2, S2, N2, O1 and K1 frequencies. In addition, surface forcing from the 3-hourly ERA-interim (ECMWF) reanalysis is used, which enables the simulation to develop a realistic mesoscale field. The horizontal resolution of the simulation is 1/60° (about 1.5 km in the mid latitudes) and features 300 vertical levels with a thickness starting from less than 1 m at the top of the ocean to 100 m at the bottom. We processed the hourly outputs of the horizontal velocity **u**, sea level, temperature and salinity fields from which we computed the pressure and stratification.

We will use the 2009 October month over two subdomains representative of different dynamical regimes of the ocean. The first one is centred over the Azores Islands and the North Mid-Atlantic Ridge. The second one is located offshore the North of the US East coast and includes a portion of the Gulf Stream. The Azores domain is characterised by the predominance of topographic features such as the mid Atlantic ridge exhibiting low amplitudes topographic variation over 100 km, and a group of seamount with strong slope and a typical scale of 10 km. Comparatively to the Gulf Stream area, the mesoscale flow is weak. The Gulf Stream domain is marked by a strong mesoscale activity and a large continental slope leading to a flat abyssal plain with a few surrounding seamounts (see Fig. 1 discussed below).

## 3.2 Filtering and Computing Methods

In order to obtain the vertical normal modes bases, the mean stratification $\bar{N}^2$ is computed from the time average of temperature and salinity using the nonlinear equation of state from TEOS-10. The Sturm-Liouville problem (6) is then discretised and solved on the staggered vertical grid in each horizontal cell, giving the two modal bases on the T grid, at the centre (in the horizontal) of each cell. From there, the modal amplitude $u_n$, $v_n$ and $p_n$ are obtained by projecting the corresponding fields $\mathbf{u}$ and $\mathbf{p}$, following (11)–(13). Since the horizontal grid is staggered, the horizontal velocity is located at the edge of each cell. We thus interpolate the base $\Phi_n$ before projecting the velocities. Unfortunately, this step induces a loss of orthogonality for the newly interpolated basis. The coefficients determined
<figure><img src="image_8.png" /><figcaption>**Fig. 1** One-month average of the horizontal energy flux divergence for mode 1 (left column) and mode 3 (right column), for the Azores domain (top) and the Gulf Stream domain (bottom). A positive value indicates a gain of energy for the mode. Black arrows show the square root of the corresponding modal horizontal energy flux. A spatial Gaussian filter with a half modal wavelength kernel has been applied on these fields in order to remove the local alternate of positive and negative values and make the plot cleaner. Green streamlines indicate, qualitatively, the mesoscale horizontal currents at 52 m depth, and isobaths at 1000, 2000, and 3000 m are superimposed in grey contours</figcaption></figure>

by innerproduct between the horizontal velocity fields and the interpolated basis are corrected *a posteriori* by inverting the cross-correlation matrix between the interpolated basis and the local basis issued from the Sturm-Liouville problem (6) on the *u*, *v*-grid. This allows us to obtain the projection coefficients *u*<sub>n</sub> and *v*<sub>n</sub> defined on the *u*, *v*-points. Vertical velocity is then computed using the continuity equation formulated with the discrete scheme employed in the NEMO code [Madec and the NEMO team 2008].

The next step consists in separating the mesoscale contribution from the semi diurnal internal tide. The latter is extracted by means of a complex demodulation at a central semi-diurnal frequency ω = 1.415×10⁻⁴ rad s⁻¹ and a low-pass cutoff cutting period of 3 days. The mesoscale flow is directly low-pass filtered with the same cutting period. The mesoscale slowly-variable buoyancy fields is computed from the daily-averaged temperature and salinity, and the corresponding stratification is then estimated from this buoyancy field and a mean sound speed profile [using standard formula, see e.g. Vallis 2017].

Since our modal decomposition is inaccurate in areas of shallow water, a mask is applied to remove all locations shallower than 300 m in the Azores domain, and
250 m in the Gulf Stream domain. The mask is slightly shallower in the Gulf Stream domain in order to retain the topographic generation occurring at the shelf break.

Since high modes are expected to propagate over shorter distances from their generation site compared to low modes, we usually group the modes 4–10 together, and interpret them as part of dissipative processes in the energy budget. Finally, all terms of Eq. (15) display small scale patterns, particularly pronounced when the mesoscale and internal tides interact. These patterns are not of particular physical significance for the diagnostics discussed in this paper. Therefore, we apply a Gaussian spatial filter before displaying any term in (15), with a kernel width equal to 25 grid points (a little less than 40 km). This procedure is only applied for map plotting.

# 4 Results

The energy budget of internal tides in both domains are examined in this section, first by considering the geographical distribution of their generation, propagation and sinks, and secondly by quantifying the spatially integrated and temporally averaged energy budget.

## 4.1 Life Cycle of the Internal Tide

Both the Azores Islands and Gulf Stream domains are sites of powerful internal tides generation, exposed to contrasted background conditions. The divergence of the energy flux averaged in time over one month (for the period of October 2009) is shown in Fig. 1, along with the mesoscale currents at 52 m depth and the topography. In the Azores area, mode 1 energy flux originates from the seamounts and propagates far away (at least 1000 km), indicating a moderate energy loss during the propagation (Fig. 1, top left panel). Modes 2 (not shown) and 3 propagate over shorter distances. All dominant sources and sinks are located around the topography (visible in Fig. 1). For the third mode, we find a local loss of energy just next to the generation at the seamount, of the same order of magnitude but smaller value. This property holds for higher modes as well (not shown). In the Gulf Stream domain, in contrast, the input of energy for mode 1 is located at the shelf break, producing a strong offshore beam. This beam manifests strong energy loss as it propagates as can be seen in the negative patch displayed in the mode 1 energy flux divergence, Fig. 1, bottom left panel. A first hot spot of sinks for the mode 1 is directly next to the shelf-break, another being at the encounter with the strong north-eastward current. This beam is also refracted by the Gulf Stream, a behaviour previously reported in Duda et al. [2018]. As we will show in the next section (cf. Figs. 2 and 3), this loss of energy is mainly caused by the advection terms in Eq. (15) and is largely converted in mode 2 energy. Mode 3 (Fig. 1, bottom right panel) displays a more
complex repartition of sources and sinks, with an overall gain of energy at and near the shelfbreak and an overall loss of energy further away from the land, where it encounters the current. The sources and sinks of energy from modes 1 and 3 are therefore much less localised around the prominent topographic feature than in the Azores domain, and clearly exhibit strong interactions with the mesoscale currents.

## 4.2 Importance of the Different Contributions in the Energy Transfers

We now focus on quantifying the causes of the multiple interactions of the internal tide with its surrounding which results in the pattern in the energy flux divergence previously described (Fig. 1).

### 4.2.1 Detailed View of Coupling Terms

To this aim, the modal energy budget (15) is integrated over the two domains considered, averaged over the October month and the different contributions are displayed under their matrix form in Fig. 3. The matrix are decomposed into their anti-symmetric part (upper triangle) denoting a transfer of energy from mode to mode, referred afterwards as scattering, and symmetric part (diagonal plus lower triangle) showing sources or sinks for the internal tide. We have only plotted terms of (15) that are of first order for at least one mode in one domain.

Amongst these terms, the topographic contribution is at least one order of magnitude higher than all others for the barotropic mode in both domains and for the first baroclinic mode in the Azores domain. However, for the Gulf Stream area, the advection coupling becomes of similar magnitude for the first mode and exceeds the topographic scattering for modes 2 and 3. In the Azores domain, this is non-negligible for modes 2 and 3 but never exceeds the topographic contribution. In this region, the second largest contribution is given by the variable stratification coupling for the baroclinic modes. In The Gulf Stream area, the mesoscale vertical shear production and the horizontal buoyancy gradient loss are also of first order. However, as already found in Kelly and Lermusiaux [2016], these two terms seem to compensate each other leaving a near zero energy gain for mode 1. The same phenomenon approximately occurs for higher modes, cf. Fig. 2: the lower triangles of buoyancy and shear matrix cancel out.

These dominant processes are dominated by their anti-symmetric part. The topographic contribution is purely anti-symmetric by construction, while the sum of the symmetric part for a given mode for the advection coupling is always smaller than 10% of the sum of the anti-symmetric part when this term is significant in regard to the topographic scattering. This quantity never exceeds 20% for the variable stratification coupling in the Azores domain. For the Gulf Stream domain
<figure><img src="image_9.png" /><figcaption>**Fig. 2** Matrices $K_{nm}$ summarising the mean energy transfer integrated over the whole subdomains between a mode *n* (row) and a mode *m* (column). The upper-right triangle is the anti-symmetric part, and the lower-left triangle (diagonal included) the symmetric part. The sign convention follows Eq. (15): a positive value indicates a gain of energy for the reference—"receiving"— mode *m* from *n* mode and conversely. Modes 4 to 10 are grouped together. From top to bottom are the coupling caused by the topography, the internal tide advection by the background flow, the terms implying the horizontal and vertical shear of the background flow and the term implying the gradient of the buoyancy field (both variable and annual mean). Notice the different range of values in the colorbar for the topographic coupling terms compared to the other ones</figcaption></figure>

however, it is close to one for the first baroclinic mode, indicating that symmetric and antisymmetric parts have comparable magnitudes. The background flow shear and buoyancy gradient couplings have a dominant symmetric part, but these two terms are cancelling each others. Last, if we sum all the contributions of the modal energy budget (15), the symmetric transfer of energy is found to be worth less or equal to 14% for modes 1, 3 in both domains and 2 in the Azores one. Therefore, the interaction between the mesoscale circulation and the internal tides mainly
contribute to a scattering of modes equivalent to the one produced by the topography without exchanges between the internal tides and the mesoscale circulation and buoyancy field.

As a last comment, Fig. 2 shows that topographic scattering is the only process that significantly transfers energy between non-neighbouring modes. For the other dominant processes, the energy transfer toward non-neighbouring modes is one order of magnitude lower than the one toward a neighbouring mode, except for the variable stratification coupling in the Gulf Stream area. For the mesoscale contributions, this implies that only the third mode is able to transfer energy toward what we consider as part of dissipative processes.

### 4.2.2 Modal Energy Budget

We are now looking at a more aggregated view of (15) thanks to Fig. 3 which shows the modal energy budget of mode 1, 2 and 3 for both domains. Here we have separated the energy transfers with the following decomposition:

**Topographic Scattering :** $C_{nm} = p_m u_n \cdot T_{nm} - p_n u_m \cdot T_{mn}$

**Advection of the internal tide :** $A_{nm} = u_m \cdot (U_{mn} \cdot \nabla u_n) - u_m \cdot u_n U_{mn}^{\Phi} + \frac{p_m}{c_m^2} U_{mn}^{p} \cdot \nabla p_n$

**Shear of background flow :** $\nabla U_{nm} = u_n \cdot U_{mn}^{\nabla} \cdot u_m + u_m \cdot U_{mn}^{z} w_n$

**Horizontal gradient of buoyancy field :** $\nabla B_{nm} = \frac{p_m}{H c_m^2} u_n \cdot (B_{mn} + B_{mn})$

**Surface contribution and advection of the mean stratification profile :** $S_{nm} = -\frac{p_m}{H c_m^2} g \varphi_m(\bar{\eta}) U_h(\bar{\eta}) \cdot \nabla \tilde{\eta} - \frac{p_m}{c_m^2} \sum_n p_n \langle \varphi_m, U_h \cdot \nabla_h \left( \frac{\tilde{\eta}^2}{c_m^2} \varphi_n \right) \rangle$

**Variable stratification :** $N_{nm} = \frac{p_m}{c_m^2} w_n \langle \varphi_m, N^2 \varphi_n \rangle$

We also have separated the symmetric part (*) of the antisymmetric one, and have furthermore divided the latter into an energy exchange coming from lower modes, and an energy exchange going from higher modes (**). The flux coming from lower modes is again divided into an energy transfer coming from the barotropic mode (****), and transfers coming from the baroclinic ones (***). A positive value means a transfer of energy toward the mode of interest.

If we consider the index of the mode of interest and $A_{nm}$ the energy transfer matrix already separated into symmetric and antisymmetric part as in Fig. 2, then we have:

$$ **** = A_{0i} $$

$$ *** = \sum_{n=1}^{i-1} A_{ni} $$

$$ ** = - \sum_{n=i+1}^{4} A_{in} $$

$$ * = \sum_{n=0}^{i} A_{in} - \sum_{n=i+1}^{4} A_{ni} $$

Comparing the integrated value of the energy divergence flux, the Gulf-Stream domain is overall less energetic than the Azores domain: mode 1 has a divergence around 4 time larger in the last domain than in the first and modes 2 and 3 are also
<figure><img src="image_10.png" /><figcaption>**Fig. 3** Mean modal budget for modes 1 (top), 2 (middle row) and 3 (bottom) for the Azores domain (left) and the Gulf Stream (right) domains, as described by Eq. (15). The divergence of the modal flux is plotted in red, while others contributions are separated into the symmetric part (green) and the anti-symmetric part where interaction with lower modes and with higher modes (orange) are separated. The anti-symmetric part coming from low modes is furthermore separated into fluxes coming from the barotropic mode (blue) and fluxes coming from the baroclinic modes (fuchsia). A positive value means a gain of energy for the mode. Processes taken into account are Topo : $C_{nm}$, Adv : $A_{nm}$, shear : $\nabla U_{nm}$, ssh + Ns : $S_{nm}$, Buo : $\nabla B_{nm}$, Strat : $N_{nm}$</figcaption></figure>

more energetic in the Azores area. In the results displayed in Fig. 3, the mode 1 in the Gulf-Stream domain is the only one to have an overall negative energy flux divergence.

The barotropic mode has an importance in the energy budget of mode 1 to 3 only through the topographic scattering, and is dominant in the energy budget of mode 1 for both areas. However, it becomes less important for modes 2 and 3, with a
contribution to the topographic scattering toward mode 3 inferior to the one caused by modes 1 and 2, but still significant.

Among the first order energy transfers exposed in all panels of Fig. 3, none transfers energy from small scale to large scale. Mesoscale and topographic couplings therefore create a forward energy cascade.

Last, Fig. 3 enables identifying the terms of the modal energy budget (15) that are of negligible importance. The advection of ssh and stationary stratification are negligible, which justifies gathering these terms in Eq. (15) in a small residual term. Similarly, the astronomical tidal potential (not included in Fig. 3) only projects weakly on the baroclinic modes. The mean of temporal variations of potential and kinetic energy are also small: one month is enough to flatten the two weeks periodic variation introduced by the spring neap cycle.

# 5 Conclusion

The semidiurnal internal tide energy budget was characterised based on a vertical mode decomposition, allowing a detailed investigation of the energy transfers amongst different vertical scales—and associated horizontal scales—of the internal tide field. We showed that the Azores area, a region of weak mesoscale activity and strong topographic features, is essentially dominated by the topographic induced energy transfer from low to higher modes. While this effect is dominant for low modes, variable stratification and—to a lesser extent—advection of the internal tide by the mesoscale flow become of similar importance for high modes (equal or larger than 2). Some questions remain concerning the precise role of the variable stratification, in relation with the definition of the modal basis expansion and how this definition would change the energy exchanges between modes. In the Gulf Stream area, the importance of the advection of the internal tide by the background flow is comparable to the topographic scattering, and even greater for modes 2 and higher. The variable stratification term is less important. However, the net impact of the mesoscale circulation on the energetic of internal tide is mainly limited to a scattering of energy from the low to high-order modes. In comparison, exchange of energy between the internal tide and the mesoscale fields are negligible: the energy exchanged by modes 1 to 3 with the background is always less than 14% of the total energy transfer for the respective mode.

The type of modal analysis conducted in the present paper on a high resolution primitive equation simulation is able to give some insights on the energetic impact of mesoscale circulation and its associated buoyancy field on the internal tide and the transfer of energy from the astronomical tide toward dissipation. This work could be extended by investigating the deep water mixing caused by these interactions. We have in particular overlooked for now the internal tide dissipation, but it is worth investigating in the near future.

Finally, the temporal variability of mesoscale—internal tides interactions will also be studied in the near future.
## References

B. K. Arbic, A. Wallcraft, and E. Metzger. Concurrent simulation of the eddying general circulation and tides in a global ocean model. *Ocean Modelling*, 32:175–187, 2010.

M. Ballarotta, C. Ubelmann, M-I. Pujol, G. Taburet, F. Fournier, J-F. Legeais, Y. Faugère, A. Delepoulle, D. Chelton, G. Dibarbouré, and N. Picot. On the resolutions of ocean altimetry maps. *Ocean science*, page 1091–1109, 2019.

L. Brodeau, J. Le Sommer, and A. Albert. Ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. Zenodo, September 2020.

T. F. Duda, Y.-T. Lin, M. Buijsman, and A. E. Newhall. Internal tidal modal ray refraction and energy ducting in baroclinic gulf stream currents. *Journal of Physical Oceanography*, 48, 2018.

M. Dunphy and K. G. Lamb. Focusing and vertical mode scattering of the first mode internal tide by mesoscale eddy interaction. *Journal of Geophysical Research: Oceans*, 33:523–536, 2014.

A. E. Gill. *Atmosphere-Ocean dynamics*. International Geophysics Series, Academic Press, Vol 30, 1982.

S. M. Kelly. The vertical mode decomposition of surface and internal tides in the presence of a free surface and arbitrary topography. *American Meteorological Society*, Volume 46, Issue 12: 3777–3788, 2016.

S. M. Kelly and P. F. J. Lermusiaux. Internal-tide interactions with the gulf stream and middle atlantic bight shelfbreak front. *Journal of Geophysical Research: Oceans*, 121:6271–6294, 2016.

G. Madec and NEMO ocean engine the NEMO team. *Note du Pôle de modélisation*. Institut Pierre-Simon Laplace (IPSL), France, No 27, ISSN No 1288-1619, 2008.

W. Munk and C. Wunsch. Abyssal recipes ii: energetics of tidal and wind mixing. *Deep-Sea Research*, Volume 45, Issue 12:1977–2010, 1998.

Y. Pan, P. J. Haley, and P. J. F. Lermusiaux. Interactions of internal tides with a heterogeneous and rotational ocean. *Journal of Fluid Mechanics*, 920, 2021.

A. L. Ponte and P Klein. Incoherent signature of internal tides on sea level in idealized numerical simulations. *Geophysical Research Letters*, Volume 42:1520–1526, 2015.

G. K. Vallis. *Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. Cambridge University Press, 2017.

E. D. Zaron, R. C. Musgrave, and G. D. Egbert. Baroclinic tidal energetics inferred from satellite altimetry. *Journal of Physical Oceanography*, page 1015–1032, 2022.

Z. Zhao, M. H. Alford, J. B. Girton, L. Rainville, and H. L. Simmons. Global observations of open-ocean mode-1 M2 internal tides. *Journal of Physical Oceanography*, 46:1657–1684, 2016.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_11.png" /></figure>
# Sparse-Stochastic Model Reduction for 2D Euler Equations

<figure><img src="image_12.png" /></figure>

Check for
updates

Paolo Cifani, Sagy Ephrati, and Milo Viviani

**Abstract** The 2D Euler equations are a simple but rich set of non-linear PDEs that describe the evolution of an ideal inviscid fluid, for which one dimension is negligible. Solving these equations numerically can be extremely demanding. Several techniques to obtain fast and accurate simulations have been developed during the last decades. In this paper, we present a novel approach that combines recent developments in stochastic model reduction and conservative semi-discretization of the Euler equations. In particular, starting from the Zeitlin model on the 2-sphere, we derive reduced dynamics for large scales and we close the equations either deterministically or with a suitable stochastic term. Numerical experiments show that, after an initial turbulent regime, the influence of small scales to large scales is negligible, even though a non-zero transfer of energy among different modes is present.

## 1 Introduction

The 2D Euler equations are a fundamental model for ideal fluids [10]. During the last two centuries, these equations have stimulated an intense activity both in terms of mathematics and physics (see for example the seminal works of Helmholtz and Arnol'd [14, 2]). In computational science and numerical analysis, retaining at a discrete level the rich non trivial structure of these equations is still

**Supplementary Information** The online version contains supplementary material available at
https://github.com/cifanip/GLIFS.

P. Cifani · S. Ephrati
University of Twente, Enschede, The Netherlands
e-mail: p.cifani@utwente.nl; s.r.ephrati@utwente.nl

M. Viviani (✉)
Scuola Normale Superiore of Pisa, Pisa, Italy
e-mail: milo.viviani@sns.it
a challenging problem [1, 19]. One main computational issue is the “curse of dimensionality”. Indeed, turbulent phenomena vary in different spatial and time scales and the distribution of energy over a vast range of scales of motion makes it computationally infeasible to fully resolve the flow in a numerical simulation. A well-established technique to mitigate large computational costs is large-eddy simulation (LES), where a spatial filter is applied to the governing equations after which only the large scales of motion are resolved [13, 21]. The filter may be defined explicitly, as a smoothing function where the filter width determines the level of detail left in the filtered solution and is chosen by the user, or implicitly, by coarsening the discretization operator. Either approach requires a closure model to represent the effect of the filtered scales on the unfiltered resolved scales. In this paper, we apply an explicit spectral cut-off filter to the 2D Euler equations, where the cut-off frequency is based on observed intrinsic scale separation. This filter leads to a discrete problem formulation in which matrix sparsity can be exploited to reduce computational costs. At the same time, this formulation still allows for an explicit representation of the effects of small scales on large scales. The equations are closed by a deterministic or stochastic model term, based on high-resolution measurements. Here, we choose to model small-scale flow features using a stochastic term mimicking high-resolution numerical data. Subsequently, we analyze the model performance by means of energy fluxes between the scales of motion. The proposed stochastic closures and assessment of energy fluxes in the high-resolution data can serve as a point of departure for further development of stochastic closure models.

A peculiar aspect of 2D ideal fluids is the presence of infinitely many conservation laws. In particular, as firstly described by Kraichnain [18], the conservation of energy and enstrophy (the $L^2$ norm of the curl of the velocity field) implies a double cascade phenomenon: the energy tends to move from small scales to large scales, whereas the enstrophy tends to follow the opposite direction. Hence, in terms of the curl of velocity, or vorticity, it is possible to clearly separate two regimes: one slowly evolving at large scales and one fast at small scales. This was shown to hold numerically in [20] for the Euler–Zeitlin model on the sphere. Theoretically, the study of non-deterministic fluid models for different regimes have gained interest in the SPDE community [12]. The equations studied in [12] and the results proved therein, show a precise connection between different space-time regimes with a reduced model for large scales. Indeed, it is shown that a suitable model for large scales is given by the so-called *Stochastic Advection by Lie Transport* (SALT) equations [15], in which a transport noise term models the infinitesimal action of the small scales on the large ones. Several numerical tests have shown the usefulness of the SALT equations as a powerful tool for model reduction [6, 9].

However, defining precisely what large and small scales are is still an open problem. In this paper, we present a criterion for defining large scales in terms of truncation of the Fourier expansion. We point out that other choices and interpretations of large and small scales are possible (see for example [20]). Let us first introduce the governing equations for the vorticity field $\omega$, defined on the
2-sphere $\mathbb{S}^2$ embedded in $\mathbb{R}^3$:

$$ 
\begin{aligned}
\dot{\omega} &= \{\psi, \omega\} \\
\Delta\psi &= \omega.
\end{aligned}
\qquad (1) $$

The Poisson bracket is defined as

$$ \{\psi, \omega\} := \nabla\psi \cdot \nabla^{\perp}\omega $$

and the Laplacian is the Laplace-Beltrami operator on $\mathbb{S}^2$. As mentioned above, equations (1) have infinitely many first integrals: energy $H(\omega) = \frac{1}{2} \int_{\mathbb{S}^2} \psi \omega$, Casimirs $C_n(\omega) = \int_{\mathbb{S}^2} \omega^n$, for $n \ge 1$, and angular momentum. Understanding the role played by these invariants is still an open problem, especially for the long-time evolution of the fluid [7].

In order to gain numerical insight on this question, V. Zeitlin proposed a spatial discretization of (1), which retains many of the first integrals above [22, 23]. The Euler-Zeitlin equations are defined as follows:

$$ 
\begin{aligned}
\dot{W} &= [P, W] \\
\Delta_N P &= W.
\end{aligned}
\qquad (2) $$

Here $W$ is a $N \times N$ skew-Hermitian matrix with zero trace, that is, an element of the Lie algebra $\mathfrak{su}(N)$. The bracket $[P, W]$ is the usual matrix commutator and the discrete Laplacian $\Delta_N$ is defined such that its spectrum is a truncation of the spectrum of $\Delta$ [17]. As mentioned above, the Euler-Zeitlin equations possess the following integral of motions: energy $H(W) = \frac{1}{2} \text{Tr}(PW)$, Casimirs $C_n(W) = \text{Tr}(W^n)$, for $n = 2, \dots, N$, and angular momentum. The core of the Zeitlin model is how the original vorticity $\omega$ and the discrete one $W$ are linked. Indeed, the representation theory of $SU(2)$ provides a deep connection between the discrete Laplacian $\Delta_N$ and a particular basis $\{T_{lm}\}$ of $\mathfrak{su}(N)$, for $l = 1, \dots, N-1$ and $m = -l, \dots, m$ [17, 4]:

* each $T_{lm}$ is an eigenvector of $\Delta_N$, with eigenvalue $-l(l+1)$,

* for each $N \ge 1$, there exists a linear map $p_N : C^\infty(\mathbb{S}^2) \to \mathfrak{su}(N)$, defined via the (real) spherical harmonics basis $\{Y_{lm}\}$ as $p_N(Y_{lm}) = T_{lm}$, if and only if $l \le N-1$,

* $\|p_N\{\psi, \omega\} - N^{3/2}[p_N\psi, p_N\omega]\| \to 0$, for $N \to \infty$, where the norm is the operator norm.

The classical way to determine large and small scales is to choose a wave number $\bar{l}$ as a threshold for the large scales (see for example [3, 6]). In this work, we propose the following criterion to set the threshold $\bar{l}$. Consider a time scale in which the fluid's energy spectrum profile has reached a stationary state. Then, typically (that is, out of equilibrium) the spectrum exhibits a double slope, which determines a kink at a certain wave number $\bar{l}$. Then, we define the large-scales $\bar{W}$ as the filtered
vorticity with modes up to $\bar{l}$, obtaining a banded matrix. We propose three possible ways, both deterministic and stochastic, of closing the equations for $\bar{W}$, by choosing different interactions with the small scales. Finally, we provide numerical tests to assess the different models introduced.

## 2 Sparse-Stochastic Model Reduction

The Euler–Zeitlin equations (2) allow to study some typical features of the 2D fluids in the matrix language. In this section, we propose a way to reduce the complexity of the Eq. (2), by defining from $W$ a sparse matrix $\bar{W}$ which retains the relevant large-scale information. Then, we show different ways of closing the equations for $\bar{W}$, by adding a suitable stochastic term.

In the Zeitlin model, the basis element $T_{lm}$ of $\mathfrak{S}(N)$ has non-zero entries only in the lower and upper $\pm m$ diagonal. If we look at the anti-diagonals, instead, we are looking at the components determining the value of the vorticity field at certain latitude bandwidth on $\mathbb{S}^2$, as shown in Fig. 1.

The large scales are typically chosen to be the modes such that $l$ is smaller than a threshold level $\bar{l}$. In the Euler–Zeitlin model, this corresponds to considering the banded matrices limited in the diagonals $\pm l \le \bar{l}$ and then removing the components corresponding to $l > \bar{l}$. The Poisson equation which defines the stream matrix $P$ preserves this sparsity structure, since the basis elements $T_{lm}$, the eigenvectors of the Laplacian, are themselves sparse. However, the Lie bracket does not restrict to this space. Indeed at each time-step we have to project the vector field into the right space.

Usually, we do not have any chance to guess the contribution of the small scales to the evolution of the large ones. However, we expect that after an initial turbulent transition, the fluid exhibits two clearly separated spatial scales. The hint for such a scenario is due to several numerical simulations of the Euler–Zeitlin equations [3, 20]. Eventually, the energy profile reaches a fixed configuration with two slopes. The first part of the spectrum represents the distribution of energy at large scales, whereas the second part the distribution of energy at small scales. Typically, the separation between large and small scales occurs at a wave number $\bar{l} \approx \sqrt{N}$.

<figure><img src="image_13.png" /><figcaption>**Fig. 1** Structure of the discrete vorticity $W$ in the Zeitlin model</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 2** Filtering of large-scale and definition of random small scale vorticity $\tilde{W}$, via the independent Brownian motions $\beta^{lm}$</figcaption></figure>

For wave numbers larger than $\bar{l}$ the energy spectrum has the characteristic slope of $l^{-1}$, which is the one of white noise, see Fig. 3. The universal nature of the small scales suggests a model reduction in terms of large-scale evolution combined with a stochastic term contribution. In Fig. 2, we show the procedure to get the two new fields $\overline{W}$ and $\tilde{W}$. To define $\overline{W}$, we introduce the orthogonal projection $\pi$ onto the modes $l \le \bar{l}$. The small-scale field $\tilde{W}$ is defined as the linear combination of the basis elements $T_{lm}$, for $l > \bar{l}$ with coefficients $\beta^{lm}$ as independent Brownian motions, with mean and variance obtained from the high-resolution DNS. This choice of coefficients ensures that the random field yields the same mean energy spectrum at the small scales as measured from the DNS. Additionally, applying the Kolmogorov–Smirnov and Anderson–Darling tests for normality to the high-resolution data suggests that the distribution of the basis coefficients for $T_{lm}$, for $l > \bar{l}$, is Gaussian.

Hence, we define $\overline{W} := \pi W$ and $\tilde{W} := \sum_{l=\bar{l}+1}^{N-1} \sum_{m=-l}^{l} \beta^{lm} T_{lm}$. With these new fields, we essentially have three possible choices. The first one consists of a deterministic closure simply via the projection of the vector field onto the large scales:

$$ \begin{aligned} \frac{\dot{\overline{W}}}{\overline{W}} &= \pi[\overline{P}, \overline{W}] \\ \Delta_N \overline{P} &= \overline{W}, \end{aligned} \qquad (3) $$
which we refer to as the no-model closure. The second model is the large scale enstrophy-preserving stochastic closure, which is up to the projection $\pi$ a type of SALT equation (see [15]):

$$ \frac{d\overline{W}}{dt} = \pi[\overline{P}, \overline{W}]dt + \sum_{l=\overline{l}+1}^{N-1} \sum_{m=-l}^{l} \frac{1}{-l(l+1)} \pi[T_{lm}, \overline{W}] \circ d\beta^{lm} \quad (4) \\ \Delta_N \overline{P} = \overline{W}. $$

We recall that the symbol $\circ$ denotes the Stratonovich integral and indeed we find that the large scale enstrophy is conserved, via

$$ d\frac{1}{2}\mathrm{Tr}(\overline{W}^2) = \mathrm{Tr}(\overline{W}d\overline{W}) = \mathrm{Tr}(\overline{W}\pi[\overline{P}, \overline{W}])dt + \sum_{l=\overline{l}+1}^{N-1} \sum_{m=-l}^{l} \frac{1}{-l(l+1)}\mathrm{Tr}(\overline{W}\pi[T_{lm}, \overline{W}]) \circ d\beta^{lm} = 0, $$

being $\pi \overline{W} = \overline{W}$ and $[\overline{W}, \overline{W}] = 0$. Notice that the other large-scale Casimirs are not preserved, since in general $\pi \overline{W}^k \neq (\pi \overline{W})^k$, for $k > 1$.

Finally, the third model is a large-scale energy-preserving stochastic closure (see [11] for its analysis and [16, 8] for more recent applications). We note that [16] introduces this stochastic closure as *Stochastic Forcing by Lie Transport* (SFLT) and provides its general definition. Here, we refer to this closure as energy-preserving noise (EPN) and adopt the following definition:

$$ \frac{d\overline{W}}{dt} = \pi[\overline{P}, \overline{W}]dt + \sum_{l=\overline{l}+1}^{N-1} \sum_{m=-l}^{l} \pi[\overline{P}, T_{lm}] \circ d\beta^{lm} \quad (5) \\ \Delta_N \overline{P} = \overline{W}. $$

The proof of conservation of the large scale energy for (5) is identical to the one of large scale enstrophy conservation by noticing that $\pi$ commutes with $\Delta_N$. The benefit of the stochastic models (4) and (5) compared to the no-model closure is that the stochastic models provide an explicit representation of the small-scale flow features and, doing so, aim to truthfully affect the evolution of the large scales of motion. No such representation is included in the no-model closure. In the next section, we perform a numerical test for the three different models (3), (4), (5), comparing them with the high-resolution DNS.

## 3 Numerical Simulations

In this section, we carry out a numerical experiment to study the performance of the models proposed in the previous section. The numerical experiment is conducted as follows. We set the high-resolution level at $N = 128$. Then we generate a random initial condition and we run a high-resolution DNS. We stop the simulation once
<figure><img src="image_15.png" /><figcaption>**Fig. 3** Initial vorticity obtained via high-resolution DNS. Top left, the field $W$, top right, the filtered field $\bar{W}$, bottom left, $W - \bar{W}$, bottom right, the energy spectrum of $W$. Note the change of slope in the energy profile at $l \approx \sqrt{N}$</figcaption></figure>

a stationary energy profile is reached (see Fig. 3). The solution at this point in time defines the initial condition of the reference solution and the ensuing model simulations.

From the DNS we select the large-scale threshold as wave number $\bar{l} \approx \sqrt{N}$, at which the kink in the energy spectrum appears. In our numerical simulation the kink is found to be at $\bar{l} = 14$.

**Remark 1** The kink in the energy spectrum must depend on the truncation level $N$. Indeed, it was shown that the tail of the energy distribution at small scales for conservative schemes, like the Euler–Zeitlin one, has a characteristic slope of $k^{-1}$. This would imply unbounded energy for $N \to \infty$, which contradicts the fact that we only consider vorticities with bounded energy. Therefore, the kink wave number between the two slopes in the energy profile must increase with $N$. Numerically, we have observed that it grows like $\sqrt{N}$.

Then, we define our large-scale field as $\bar{W} := \pi W$, where $\pi$ denotes the orthogonal projection onto the modes for $l \le \bar{l}$. The projection consists of two steps: first we extract the components up to $\bar{l}$ and then we generate the field $\bar{W}$. The cost of calculating each component is $O(N)$ and since we need to repeat this operation $\bar{l}^2 - 1 \approx N$ times, the total cost of extracting the components is $O(N^2)$. Clearly, to construct the field $\bar{W}$ we have to perform $O(N^2)$ operations. Hence, the total computational cost of the projection $\pi$ is $O(N^2)$. Hence, given two matrices $A, B \in \text{su}(N)$, the cost of evaluating $\pi[\pi A, \pi B]$ is given by the evaluation of $\pi$ plus the cost of multiplying $\pi A$ and $\pi B$. Since we are interested only in the $\pm \bar{l}$
diagonals we need to perform $O(N\bar{l})$ vector-vector multiplications of the cost $O(\bar{l})$, which implies a total cost of $O(N\bar{l}^2) \approx O(N^2)$. The same cost $O(N^2)$ holds for the stochastic term, since we can actually consider only $m = -\bar{l}, \dots, \bar{l}$. We also define $\tilde{W}$, as explained in the previous section. Finally, we define the reference solution as the high-resolution numerical simulation using the initial condition previously defined. The large-scale field obtained by projecting this initial condition serves as a starting point for the simulations where the small scales are modeled as described in Eqs. (3), (4), (5). All simulations are run for 250 time units from the initial condition. In our numerical simulations, the time integration is done via the Heun-type scheme adapted for the SDEs, with time-step $h = 0.25$. In the following, we perform an ensemble of realizations for the stochastic closures. Each ensemble consists of 25 realizations, which is sufficient to show the qualitative difference between the proposed models.

We notice from Figs. 4 and 5 that the no-model solution and the mean SALT solution perform well compared to the reference solution, in terms of qualitative vorticity evolution as well as the kinetic energy spectrum. On the contrary, the energy-preserving scheme loses accuracy and a cascade of energy to lower wave numbers occurs. We attribute the observed difference between SALT and EPN to two causes. Firstly, the orders of magnitude of the stochastic terms entering the equations differ between the methods. The stochastic forcing types are parametrizations of components of the effect of small scales on large scales, respectively given by $\pi[\tilde{P}, \tilde{W}]$ and $\pi[\bar{P}, \bar{W}]$. Figure 6 shows these quantities for a high-resolution snapshot, which illustrates that the term parametrized by the EPN is significantly larger than the term that SALT parametrizes. Therefore, it is reasonable to expect that the proposed use of EPN leads to more substantial deviations from the no-model

<figure><img src="image_16.png" /><figcaption>**Fig. 4** Evolution of the large scales over 250 time units, via the different models proposed and the high-resolution simulation. No model corresponds to (3), SALT to (4) and EPN to (5). The shown results using SALT and EPN are the mean of an ensemble each consisting of 25 independent realizations</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 5** Energy spectra at $t = 250$ using SALT (left) and EPN (right), both computed from an ensemble of 25 independent realizations and compared to the reference and no-model energy spectra. Note that the standard deviation of the ensemble obtained using SALT is too small to discern in this figure</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 6** Components of the evolution of the reference vorticity projected onto the large scales, generated from the reference vorticity field at $t = 250$. Note that the ranges of the color bars vary per field, which highlights the difference in magnitudes between the fields</figcaption></figure>

simulation than the use of SALT. In fact, the proposed use of SALT only leads to very small changes compared to the no-model setting. Secondly, in the energy-preserving scheme, no energy can leave the large scales. Hence, if the transfer of energy between different modes is non-zero, the conservation of the large-scale energy prevents the energy to flow from large scales to small scales, causing an extra accumulation of energy $l \approx \bar{l}$.
In order to check this thesis, we compute the energy transfer among different modes in the high-resolution DNS. Let us consider the energy at a level $l$:

$$E(l) = \frac{1}{2} \sum_{m=-l}^{l} \frac{\omega_{lm}^2}{l(l+1)}.$$

Then, the energy variation in time is given by

$$\frac{dE(l)}{dt} = \sum_{m=-l}^{l} \frac{\omega_{lm}[P, W]_{lm}}{l(l+1)}.$$

Let $F(l) := |\frac{dE(l)}{dt}|$ be the absolute value of the energy transfer due to the nonlinearity of the vector field $[P, W]$. In Fig. 7, we plot the energy transfer contributions of the four possible couplings of large and small scales. We notice that the transfer of energy between large and small scales is non-zero. In particular, the main driver of the energy for the components of $\tilde{W}$ is the vector field $[\tilde{P}, \tilde{W}]$, whereas for small scales it is $[\tilde{P}, \tilde{W}]$.

The plots in Fig. 7 show that the term $[\tilde{P}, \tilde{W}]$ becomes more and more relevant in the energy flux at large scales, while approaching the threshold level $\tilde{l}$. Therefore, one might expect the stochastic model (5), which takes into account this term too, to be more accurate than (3) or (4). However, form Figs. 4 and 5 this does not seem the case. Indeed, the way we model the term $[\tilde{P}, \tilde{W}]$ in (5) prevents the energy from flowing from large scales to small scales. However, also in (3) there is no energy flow from large scales to small scales. Therefore, we suggest that at large scales the term $[\tilde{P}, \tilde{W}]$ is responsible for redistributing the energy from lower to higher frequencies and in absence of a dissipation mechanism the rearranging of energy diverges from the correct spectrum.

The fact that (3) and (4) perform equally well is quite surprising and indicates that the small scales do not affect the large scales much when a stationary energy profile

<figure><img src="image_19.png" /><figcaption>Fig. 7 Energy transfer among different modes computed from a snapshot of the high-resolution simulation. Left, energy transfer among modes at large scales, right, at small scales</figcaption></figure>
is reached. However, for much longer times than those we have run, it is possible
that the effect of small scales on large scales become more relevant and (4) becomes
more accurate than (3) in terms of large-scale dynamics. Further investigations on
this aspect are ongoing research and will be presented in future work.

4 Conclusions and Outlook

In this paper, we have presented a possible strategy to reduce the complexity of the Euler–Zeitlin model, while performing long-time simulations. Numerical evidence shows that the Euler–Zeitlin equations exhibit a clear separation of scales such that the large-scale dynamics are quite robust to different couplings with small scales, either deterministic or stochastic. Interestingly, the energy-preserving scheme we have defined shows that the energy at large scales cannot be exactly conserved. This means that large and small scales are never completely decoupled, even when the energy spectrum profile reaches a stationary regime. This indicates that for very long times a non-zero transfer of energy among different scales is present.

The Zeitlin model has been criticized for unrealistic conservation of enstrophy and other Casimirs at a finite level of truncation *N*. Our result shows that this issue can be understood such that the Euler–Zeitlin equations are quite robust and precise in describing large scales, which means for wave numbers *l* ≈ √*N*. On the other hand, the remaining modes are themselves a model for the small scales, which correctly mimic the energy flux among different modes.

In conclusion, we have shown that the Zeitlin model can be a useful tool for
simulating long-time large-scale dynamics. In future work, we aim to perform more
systematic simulations using the parallelized code developed in [5] and available
on https://github.com/cifanip/GLIFS. Additionally, further analysis of energy and
enstrophy transfers between the scales of motion may serve to derive tailored data-
driven stochastic closure models for the Euler-Zeitlin equations.

References

[1] A. Arakawa. Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. *J. Comp. Phys.*, 1(1):119–143, 1966.

[2] V. I. Arnold. Sur la géometrie différentielle des groupes de lie de dimension infinie et ses applications à l'hydrodynamique des fluids parfaits. *Ann. Fourier*, 16(1):319–361, 1966.

[3] G. Boffetta and R. E. Ecke. Two-dimensional turbulence. *Annu. Rev. Fluid Mech.*, 44(1): 427–451, 2012. doi: 10.1146/annurev-fluid-120710-101240.

[4] M. Bordemann, J. Hoppe, P. Schaller, and M. Schlichenmaier. gl(∞) and geometric quantization. *Comm. Math. Phys.*, 138(2):209–244, 1991.

[5] P. Cifani, M. Viviani, and K. Modin. An efficient geometric method for incompressible hydrodynamics on the sphere. *Journal of Computational Physics*, 473:111772, 2023. ISSN 0021-9991.

[6] Colin Cotter, Dan Crisan, Darryl D. Holm, Wei Pan, and Igor Shevchenko. A particle filter for stochastic advection by lie transport: A case study for the damped and forced incompressible
two-dimensional euler equation. *SIAM/ASA Journal on Uncertainty Quantification*, 8(4): 1446-1492, 2020.

[7] M. Dolce and T. Drivas. On maximally mixed equilibria of two-dimensional perfect fluids. *Arxiv*, 2022.

[8] Theodore D. Drivas and Darryl D. Holm. Circulation and energy theorem preserving stochastic fluids. *Proceedings of the Royal Society of Edinburgh Section A: Mathematics*, 150(6):2776-2814, 2020. doi: 10.1017/prm.2019.43.

[9] S.R. Ephrati, P. Cifani, E. Luesink, and B.J. Geurts. Data-driven stochastic lie transport modeling of the 2d euler equations. *Journal of Advances in Modeling Earth Systems*, 51, 2023.

[10] L. Euler. Principes généraux de l'état d'équilibre d'un fluide. *Académie Royale des Sciences et des Belles-Lettres de Berlin, Mémoires*, 11:217-273, 1757.

[11] Franco Flandoli and Dariusz Gatarek. Martingale and stationary solutions for stochastic navier-stokes equations. *Probability Theory and Related Fields*, 102(3):367-391, 1995.

[12] Franco Flandoli and Umberto Pappalettera. 2d euler equations with stratonovich transport noise as a large-scale stochastic model reduction. *Journal of Nonlinear Science*, 31, 2021.

[13] Bernardus J Geurts. *Elements of direct and large eddy simulation*. RT Edwards, Inc, 2003.

[14] H. Helmholtz. Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen. *J. Reine Angew. Math.*, 25-55:25-55, 1858.

[15] D. D. Holm. Variational principles for stochastic fluid dynamics. *Proc. R. Soc. A.*, page 471(20140963), 2015.

[16] Darryl D Holm and Ruiao Hu. Stochastic effects of waves on currents in the ocean mixed layer. *Journal of Mathematical Physics*, 62(7):073102, 2021.

[17] J. Hoppe and S.-T. Yau. Some properties of matrix harmonics on $S^2$. *Comm. Math. Phys.*, 195:66-77, 1998.

[18] R. H. Kraichnan. Inertial ranges in two-dimensional turbulence. *Phys. Fluid.*, 10(7):1417-1423, 1967.

[19] K. Modin and M. Viviani. A Casimir preserving scheme for long-time simulation of spherical ideal hydrodynamics. *J. Fluid Mech.*, 884:A22, 2020.

[20] Klas Modin and Milo Viviani. Canonical scale separation in two-dimensional incompressible hydrodynamics. *Journal of Fluid Mechanics*, 943:A36, 2022.

[21] Pierre Sagaut. *Large eddy simulation for incompressible flows: an introduction*. Springer Science & Business Media, 2006.

[22] V. Zeitlin. Finite-mode analogues of 2D ideal hydrodynamics: Coadjoint orbits and local canonical structure. *Physica D*, 49(3):353-362, 1991.

[23] V. Zeitlin. Self-consistent-mode approximation for the hydrodynamics of an incompressible fluid on non rotating and rotating spheres. *Phys. Rev. Lett.*, 93(26):353-362, 2004.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
# Effect of Transport Noise on Kelvin-Helmholtz Instability

<figure><img src="image_20.png" /></figure>

Check for
updates

Franco Flandoli, Silvia Morlacchi, and Andrea Papini

**Abstract** The effect of transport noise on a 2D fluid may depend on the space-scale of the noise. We investigate numerically the dissipation properties of very small-scale transport noise. As a test problem we consider the Kelvin-Helmholtz instability and we compare the inviscid case, the viscous one, both without noise, and the inviscid case perturbed by transport noise. We observe a partial similarity with the viscous case, namely a delay of the instability.

## 1 Introduction

Stochastic transport is a new fundamental perspective on fluid dynamics, see e.g. [29, 39, 13, 8] and [22]. A transport type noise in a fluid dynamic model may be seen, loosely speaking, as a simplified description of small-medium space-scales of motion. In the numerical simulations (see for instance [8]) we may observe the way it perturbs large scale motion; in general, this perturbation destabilizes large scales producing smaller eddies.

This has been observed for homogeneous noise where small scale perturbations with possibly local backscattering of energy are observed. While there are cases in which transport noise has been shown also to trigger secondary circulation at large scales for inhomogeneous noises [3].

In this note we want to explore how transport noise may affect large scales, somewhat opposite to the one mentioned above. The key difference is the assumption that it is very-small-space-scale. The noise used below is made of very small, low intensity, vortex structures. In such a case it may happen that the transport noise acts as a dissipation, an additional viscosity.

It corresponds to Joseph Boussinesq intuition [6] that “turbulent small scales may be dissipative on the mean flow”. The physical intuition, beyond the specific

<footnote>F. Flandoli (✉) · S. Morlacchi · A. Papini
Scuola Normale Superiore, Pisa, Italy
e-mail: franco.flandoli@sns.it; silvia.morlacchi@sns.it; andrea.papini@sns.it</footnote>
mathematical derivation, is that fluid particles move so erratically to produce effects similar to the molecular motion. In a sense, it is like a macroscopic version of the molecular dynamics: as macroscopic kinetic energy transfers to molecular kinetic energy (heat), reducing the macroscopic motion, similarly large-scale kinetic energy moves to small scale turbulence and reduce the intensity of the mean flow. There have been many attempts to prove the validity of this picture, see for instance [5, 16, 21, 32, 30, 44], but a final conclusion is not clear, also because the empirical validity is sometimes moderate [32, 41]; in particular, the inverse cascade in 2D seems to break the viscosity effect due to small-scale turbulence [14].

However, this idea is certainly very useful for numerical simulations of fluid dynamics model, being the basis of the LES method [5], and, in case of vorticity equations, the vortex blob method [11].

Theoretically, this kind of turbulent-small-scale transport noise has been investigated in rigorous works. This line of research was initiated in [23] and developed in several works, see e.g [16, 19] and many others (see [22]).

However, its stabilizing power has not been tested numerically yet. Here we observe its action on one of the strongest and most common instabilities: the Kelvin-Helmholtz one. This is an instability that occurs when a velocity difference across the interface between two fluids is present, producing macroscopic vortex-like structures. For a description of Kelvin-Helmholtz instability see [31] or [36], Chapter 6, and our description below; as a typical instability of shear flows, see also [14].

We compare the inviscid case, the viscous one, both without transport noise, and the inviscid case perturbed by transport noise. The results are described in Sect. 4. Since we approximate the fluid dynamic equations by vortex methods, which fits quite well and in a unified way with the three cases analyzed here (inviscid, viscous and stochastic transport), we describe some preliminaries on this topic in Sects. 2 and 3.

## 2 Model Formulation

We consider the two-dimensional flow of an incompressible fluid on a 2D domain, either the full plane or $T^2$, the 2D torus with periodic boundary conditions.

As usual, the equations for the motion of the fluid are the conservation of mass and linear momentum, see e.g. [2], expressed through the Navier-Stokes equations in the null-divergence formulation. The main interest for our numerical simulations is the equation for the evolution of vorticity $\omega(t, x)$, which can be derived from the Navier-Stokes equations:

$$ \partial_t \omega + u \cdot \nabla \omega = \nu \Delta \omega. \qquad (1) $$

Here $\nu$ is the kinematic viscosity, and $u$ is the velocity of the fluid solving the NS-equation. In the 2D case, we can express $\omega := \nabla \times u$ as the curl of the velocity
field, a vector normal to the flow plane. We denote $x = (x_1, x_2)$ to be any point in the domain.

## 2.1 Point Vortex Method for Inviscid Flows

In order to state the problem we investigate numerically, we first recall several well-known facts; for a general introduction, see e.g. [36, 37, 43]. For the sake of simplicity, we focus on an inviscid fluid; the vorticity equation (1) reduces to,

$$ \partial_t \omega + u \cdot \nabla \omega = 0, \qquad (2) $$

From the null-divergence hypothesis, if $u$ has zero average we can define a stream function associated to the fluid $\psi(t, x)$, see [2], such that the velocity $u$ is given by $u = -\nabla^{\perp} \psi$. We obtain the stream function by solving the Poisson equation

$$ \Delta \psi = -\omega. \qquad (3) $$

To compute the solution, we need to express the Green's function and the convolution with the variations of constants methods, which for flows in the full plane $\mathbb{R}^2$ reads as

$$ \psi(t, x) = \int G(x - y) \omega(t, y) dy := G * \omega, \qquad (4) $$

where $G$ is the Green's function, the fundamental solution of the Laplace equation. The explicit form of $G$ in the full plane $\mathbb{R}^2$ is

$$ G(x, y) = \frac{1}{2\pi} \log(|x - y|). \qquad (5) $$

Using (4) we obtain the velocity field

$$ u(t, x) = \int K(x - y) \omega(t, y) dy := K * \omega, \qquad (6) $$

where $K$ is given by

$$ K(x) = \nabla^{\perp} G(x). \qquad (7) $$

Equations (6) and (7) for the velocity are known as the Biot-Savart law, and $K$ is the Biot-Savart Kernel. Note that we must correct this velocity for flows over bounded or periodic boundaries domains to satisfy the boundary condition (See Sect. 4.1, for details). Therefore, we have to change the Green's function according to the Poisson equation.
Consider now a fluid particle $X_t$, moving in the velocity field; from (2), the path of the particle is

$$ \frac{d}{dt} X_t = u(t, X_t) \\ \omega(t, X_t) = \omega_0(X_0). $$

Therefore, considering the fluid as distinct “fluid particles” of constant vorticity, these abstract objects’ motion determines the scalar field’s evolution. This is the premise of the point vortex method to compute (2), (see [36, 37] and [43] for details).

To this end, consider $N$ point vortices, idealizing a 2D inviscid fluid and occupying positions $X_t^1, ..., X_t^N$, with intensities (circulations) $\Gamma_1, ..., \Gamma_N$ respectively. They move accordingly to the following set of ordinary differential equations derived from the previous computations (8):

$$ \frac{dX_t^i}{dt} = \sum_{j \neq i} \Gamma_j K(X_t^i, X_t^j) \qquad (8) $$

where the vector-valued kernel $K(x, y)$ is the Biot-Savart kernel, equal to $\frac{1}{2\pi} \frac{(x-y)^\perp}{|x-y|^2}$ in full space, suitably modified on a torus or in a bounded domain (see Sect. 4.1). We refer to [15] for well-posedness of the ODE system for almost all initial configurations (including the case when the system is perturbed by additive noise). One can prove that the empirical measure

$$ \omega(t, \cdot) := \sum_{i=1}^{N} \Gamma_i \delta_{X_t^i} $$

is a weak solution of 2D Euler equations in vorticity form (suitably interpreted for distributional fields as in [42]), see [24] and [26, 27] for a discussion of well-posedness in this context):

$$ \begin{gathered} \partial_t \omega + u \cdot \nabla \omega = 0 \\ u(t, x) = \int K(x, y) \omega(t, y) dy \qquad (9) \\ \omega|_{t=0} = \omega_0 \end{gathered} $$

with $\omega_0 := \sum_{i=1}^N \Gamma_i \delta_{X_0^i}$. Here $\omega$ is the (scalar) vorticity, $u$ is the (vector) velocity.

Moreover, consider the empirical measure parametrized by $N$ and rescaled by $1/N$:

$$ \omega_N(t, \cdot) = \sum_{i=1}^{N} \frac{\Gamma_i}{N} \delta_{X_t^i}. $$
In [36] a convergence result to Euler equations was proved under some restrictions (specifically, in that case $X_t^i$ solves a vortex dynamics with mollified Biot-Savart Kernel, with vanishing mollification as $N \to \infty$, with rates of convergence of the parameters a bit far from optimal). The best result in the literature of convergence to Euler equations was then [42], see also other references therein.

Notice that the velocity field associated to the distributional vorticity $\omega_N(t, \cdot)$ is

$$u_N(t, x) := \frac{1}{N} \sum_{i=1}^{N} K(x, X_t^i).$$

This is a well defined vector field, of class $L_{loc}^p$ for every $p < 2$ but not for $p = 2$.

## 2.2 Point Vortex Method for Viscous Flows

To investigate viscous flows, following [10, 33, 38], we modifying the previous scheme by adding independent 2D Brownian motions $W_t^1, ..., W_t^N$ to the equations of point vortices

$$dX_t^i = \sum_{j \neq i} \Gamma_j K(X_t^i, X_t^j) dt + \sqrt{2\nu} dW_t^i. \quad (10)$$

Then, the empirical measure $\omega_N(t, \cdot)$ converges weakly, in probability, to the unique solution of the 2D Navier-Stokes equations in vorticity form

$$\partial_t \omega + u \cdot \nabla \omega = \nu \Delta \omega$$

$$u(t, x) = \int K(x, y) \omega(t, y) dy \quad (11)$$

$$\omega|_{t=0} = \omega_0$$

We refer to [7] for convergence to 2dNS on the full plane, and to [33, 38, 25] for results in the case of domains with boundary.

## 3 Point Vortex Method with Environmental Noise

As mentioned above in the Introduction, several works, e.g. [29, 39, 13, 28, 9, 8] and [22], indicated an interest in the following stochastic modification of the Euler equations
$$d\omega + u \cdot \nabla \omega dt = \sum_{k \in K} \sigma_k \cdot \nabla \omega \circ dB_t^k \quad (12)$$

where $\sigma_k = \sigma_k(x)$ are given vector fields, that we assume divergence-free, $(B_t^k)_{k \in K}$ are independent 1D Brownian motions and the stochastic operation $\circ$ stands for the Stratonovich integral. Due to this, formally, vorticity is conserved (it is transported randomly by the field $u \, dt + \sum_{k \in K} \sigma_k \circ dB_t^k$). This model bears similarities with the viscous flows, and the objective of this paper is to show differences and similarities of the elliptic operator obtained from such a transport-advection noise.

## 3.1 Transport Noise and Deterministic Scaling Limit

The point vortex dynamics associated to the model (12) is then given by the following expression:

$$dX_t^i = \frac{1}{N} \sum_{j \neq i} K(X_t^i, X_t^j) dt + \sum_{k \in K} \sigma_k (X_t^i) \circ dB_t^k. \quad (13)$$

Notice that this is a model of common noise (also called environmental noise): the BM's $B_t^k$ are the same for all particles, in contrast to the model (11) where each particle $X_t^i$ was affected by an independent BM $W_t^i$. See [18] for an example of theoretical results on this model. For models similar to this one, it has been proved (see e.g. [12]) that the empirical measure converges to the solution of the SPDE (12). At the same time, following [23] and subsequent works, if the noise is parameterized in such a way to become more and more small scale, the SPDE (12) converges to the deterministic equation with additional viscosity (see Sect. 3.2 for details on possible noise selection)

$$\partial_t \omega + u \cdot \nabla \omega = \nu \Delta \omega. \quad (14)$$

Inspired by [20], we consider a sort of mixed scaling limit: we take the point vortex dynamics with common noise (13), which for given fields $\sigma_k$ would converge to the SPDE (12), and we choose more and more small scale coefficients $\sigma_k$ in order to be close to the deterministic equation (14).

The scheme just described has two parallel interpretations. If we look at the point vortex dynamics just as a numerical method of discretization of Euler equation, the transport noise in the point vortex dynamics (13) is just like the transport noise in the Euler vorticity dynamics (12). Therefore the scaling limit when the coefficients $\sigma_k$ more concentrated is a numerical realization of the theoretical scaling limit investigated in [16], from the stochastic Euler equation to the deterministic Navier-Stokes equation and, as such, it is a form of proof of the validity of
Boussinesq hypothesis that small scale turbulence enhances the viscosity. That is the first interpretation. The second one is more related to particle systems and noise. Consider the particle system (13), having common noise and compare it with the more classical particle system with independent noise acting on each particle (Sect. 2.2 above). In the scaling limit of $\sigma_k$ described above, the common noise becomes more and more spatially de-correlated, close to independent noise at each space position, and thus it acts almost independently on each particle. Since point vortices with independent noise converge to the Navier-Stokes equations, it is natural to expect that also with common noise but very concentrated $\sigma_k$ the empirical measure of the particle system is close to the solution of the Navier-Stokes equations. This viewpoint has been investigated theoretically in [20], while the subject of linking the turbulence stresses to the mean flow, which is in synthesis the Boussinesq idea, was explored in appendix A and B of [40] and theoretically in [17], providing a link to eddy viscosity and an environmental small scale noise.

## 3.2 A Digression on the Theoretical Selection of the Noise

In this section, in the same spirit as [20, 22, 40], we explore some property of the environmental noise that we use in a simplified way in our numerical simulations. Following the works on modeling of passive scalars [34], when considering the scaling limit of (13) to $\omega(\mathbf{x})$ solution of the viscous Euler equation (14), we consider a model of noise in the fluid which is delta-correlated in time, namely a white noise with a precise space dependence.

$$ \mathbf{W}(t, \mathbf{x}) dt = \sum_{k \in K} \sigma_k(\mathbf{x}) dB_t^k \quad (15) $$

where $(\sigma_k(\mathbf{x}))_k$ is a family of smooth divergence free vector fields on the 2D domain of the equation, and $B_t^k$ are independent one-dimensional Brownian motions; $K$ is, usually, a finite index set, but with suitable assumption we could consider also the case of countable family of smooth fields.

In this case, the term $\mathbf{W}(t, \mathbf{x}) \cdot \nabla \omega(\mathbf{x})$ obtained in the convergence result of the point vortex empirical measure, must be interpreted as a Stratonovich integral

$$ \sum_{k \in K} \sigma_k(\mathbf{x}) \cdot \nabla \omega(\mathbf{x}) \circ d B_t^k. \quad (16) $$

Assume that the solution is sufficiently smooth so that the Stratonovich integral makes sense, then this is given by an Itô-Stratonovich corrector plus an Itô integral; precisely, is given by:
$$-\frac{1}{2} \sum_{k \in K} \sigma_k(\mathbf{x}) \cdot \nabla (\sigma_k(\mathbf{x}) \cdot \nabla \omega(\mathbf{x}, \mathbf{v})) dt + dM(t, \mathbf{x})$$

where $M(t, \mathbf{x})$ is a (local) martingale. Follows that the Itô-Stratonovich corrector takes the form of an elliptic operator:

$$-\frac{1}{2} \operatorname{div} (C(\mathbf{x}, \mathbf{x}) \nabla \omega(\mathbf{x})) dt$$

where $C(\mathbf{x}, \mathbf{y})$ is the space-covariance function of the noise

$$C(\mathbf{x}, \mathbf{y}) = \sum_{k \in K} \sigma_k(\mathbf{x}) \otimes \sigma_k(\mathbf{y}).$$

As an example, we take the noise [34], which is relevant to our numerical investigation in the choice of the divergence-free field in the point vortex model. For simplicity, assume the domain to be $\mathbb{R}^2$, but modifications on $\mathbb{T}^2$ are possible, see for example [16, 20].

Its covariance function is space-homogeneous, i.e. $C(\mathbf{x}, \mathbf{y}) = C(\mathbf{x} - \mathbf{y})$, with the form

$$C(\mathbf{z}) = \nu k_0^{\zeta} \int_{k_0 \le |\mathbf{k}| < k_1} \frac{1}{|\mathbf{k}|^{d+\zeta}} e^{i \mathbf{k} \cdot \mathbf{z}} \left( I - \frac{\mathbf{k} \otimes \mathbf{k}}{|\mathbf{k}|^2} \right) d\mathbf{k}.$$

The famous Kolmogorov 41 case follows if we take $\zeta = 4/3$. Taking $k_1 = +\infty$, then $C(\mathbf{0}) = K\sigma^2$ where the constant $K$ is given by

$$K = \int_{1 \le |\mathbf{k}| < \infty} \frac{1}{|\mathbf{k}|^{d+\zeta}} \left( I - \frac{\mathbf{k} \otimes \mathbf{k}}{|\mathbf{k}|^2} \right) d\mathbf{k}.$$

We consider small-scale turbulent velocity fields depending on a scaling parameter and taking the scaling limit in (12), as in [20, 23]. In the case of [34] we have

$$k_0 = k_0^N \to \infty$$

The result $C(\mathbf{0}) = K\nu$ is independent of $N$, so that the Itô-Stratonovich corrector becomes equal to

$$\nu \Delta \omega(\mathbf{x}),$$

and simultaneously, we may have that the Itô term goes to zero, hence recovering (14).

Let us sketch the argument (see [16] for details) which explains why the Itô term may go to zero, in spite of the convergence to a finite non-zero limit of the Itô-Stratonovich corrector. Let $\phi$ be a smooth test function. One has
$$ \mathbb{E} \left[ \left( \sum_{k \in K} \int_0^T \langle \sigma_k \cdot \nabla \omega_t, \phi \rangle_{L^2} dB_t^k \right)^2 \right] = \mathbb{E} \left[ \sum_{k \in K} \int_0^T \langle \sigma_k \cdot \nabla \omega_t, \phi \rangle_{L^2}^2 dt \right] $$

by the isometry formula of Itô integrals,

$$ = \mathbb{E} \left[ \sum_{k \in K} \int_0^T \langle \omega_t, \sigma_k \cdot \nabla \phi \rangle_{L^2}^2 dt \right] $$

since $\operatorname{div} \sigma_k = 0$,

$$ = \mathbb{E} \left[ \int_0^T \int \int \sum_{k \in K} \sigma_k(x) \cdot \nabla \phi(x) \sigma_k(y) \cdot \nabla \phi(y) \omega(t, x) \omega(t, y) dx dy dt \right] $$

$$ = \mathbb{E} \left[ \int_0^T \int \int \nabla \phi(y)^T \cdot C(x, y) \cdot \nabla \phi(x) \omega(t, x) \omega(t, y) dx dy dt \right] $$

$$ = \mathbb{E} \int_0^T \langle \mathbf{C} \theta_t, \theta_t \rangle_{L^2} dt $$

where is the linear operator on vector fields with kernel $C(x, y)$ and $\theta_t(x) = \nabla \phi(x) \omega(t, x)$,

$$ \leq \|\mathbf{C}\|_{L^2 \to L^2} \mathbb{E} \int_0^T \|\theta_t\|_{L^2}^2 dt. $$

Now, one can prove uniform bounds on $\mathbb{E} \int_0^T \|\theta_t\|_{L^2}^2 dt$ with respect to the scaling of the noise and one can choose a noise such that $\|\mathbf{C}\|_{L^2 \to L^2}$ goes to zero. Notice that in the Itô-Stratonovich corrector only the diagonal $C(x, x)$ counts, while the smallness of $\|\mathbf{C}\|_{L^2 \to L^2}$ is related to the smallness of $C(x, y)$ when $x \neq y$.

# 4 Numerical Results

## 4.1 Setting: Kelvin–Helmholtz Instability

In this section, we investigate classical results on the shear flow model in the setting of point vortices, analyzing the Kelvin–Helmholtz instability and the possibility of delaying the structure formation. In this way, we can both test the goodness of our point vortex models and, at the same time, set a benchmark for which we will show delayed instability. In order to test the point vortex model in the classical cases (8) and (10), we choose the particular fluid configuration of a shear flow because of its
fundamental property: developing instability without viscosity and delaying it when viscosity is present [36].

We work on the torus $\mathbb{T}^2$ equal to the set $[-1, 1]^2 / \sim$ with coordinates $x = (x_1, x_2)$ and identified boundaries at $x_1, x_2 = \pm 1$; all fields are periodic in the $x_1$ and $x_2$-direction. We take an initial velocity $u_0$ of the form

$$u^0(x_1, x_2) = (u_1^0(x_2), 0)$$

and corresponding vorticity $\omega_0 = \partial_{x_2} u_1^0(x_2)$. We choose, in particular, the function

$$u_1^0(x_2) = \begin{cases} -1 & \text{if } x_2 \le -\delta \\ \frac{x_2}{\delta} & \text{if } -\delta \le x_2 \le \delta \\ 1 & \text{if } \delta \le x_2 \end{cases} \quad (17)$$

where we fix $\delta = 0.02$ in our numerical simulations.

To compute the vorticity measure of our point vortices, we use vortex blobs, obtained by spreading the circulation of a point vortex over a chosen small area, the vortex core (see e.g. [43]). In this formulation, the vorticity field is approximated by

$$\omega_{\varepsilon}^{N}(\mathbf{x}, t) = \sum_{i} \Gamma_{i} \phi_{\varepsilon}(\mathbf{x} - X_{t}^{i}), \quad (18)$$

where the mollifier $\phi_{\varepsilon}$ (in our numerical simulations a Gaussian kernel with width dependent on the subscript $\varepsilon$, the characteristic size of the vortex core) describes the vorticity distribution in the vortex core. Following standard numerical techniques (see e.g. [4, 1]), the core size $\varepsilon$ of the vortices has to be much larger than the average spacing $d$ between the vortices; the core size is usually taken to be $\varepsilon = d^q$, with $q \ll 1$.

In (18), the vorticity distribution at any time depends on the point vortices $X_t^i$ through the vortex blobs. In our numerical simulations, we take $N \sim 10^4$ point vortices; following a mean-field approach, the initial circulation for every given point vortex is derived from $u_1^0$ and is equal to $\Gamma_0^i = \frac{1}{2\delta N}$. We solved the point vortex model (8) using a 2nd order Runge Kutta scheme coupled with a Heun techniques for the noise, for a second-order time discrete approximation; the time step for our simulations was selected to be $\Delta t \sim 10^{-3}$ to ensure a trade-off between the stability of our method and the generation of vortex-like structures in the shear flow model.

In the usual way, we also recall that in our numerical framework, the kernel $K$ in (13) corresponds to the Biot-Savart kernel. We have that $K = \nabla^{\perp} G = (\partial_2 G, -\partial_1 G)$, where $G$ is the Green function on $\mathbb{T}^2$. In the whole plane we have the simple expression $G_{\mathbb{R}^2} = \frac{1}{2\pi} \log |x|$; while for our domain we know that

$$G(x) = \frac{1}{2\pi} \log |x| + s(x), \quad \forall x \in \mathbb{T}^2 \setminus \{0\}, \quad (19)$$
and $s(x)$ is a smooth function on $\mathbb{T}^2$. Thus, $K$ is divergence-free, smooth away from the origin, and symmetrical; moreover, it holds the following behaviour:

$$|K(x)| \sim 1/|x|, \text{ as } |x| \to 0,$$

which we extensively use to approximate our kernel with $K_{\mathbb{R}^2}(x-y) = \frac{1}{2\pi} \frac{(x-y)^\perp}{|x-y|^2}$, throughout the numerical simulations. Without ambiguity, from here on, we consider the horizontal and vertical axes as our reference frame, naming them the x-axis and y-axis, as usual.

### 4.1.1 The Role of Intrinsic Instability

We know that, at least formally, the vector field $u = u^0$ is a solution of Euler equation (14) with $v = 0$. This system is unstable: small perturbations rapidly develop vortex blobs. As such, we consider the system of point vortices $(X_t^i)_i$ with initial vorticity, derived from (17), and expressed as

$$\omega_0(x_1, x_2) := \frac{1}{N} \sum_i \frac{1}{2\delta} \delta_{x_0^i}(x_1, x_2), \qquad (20)$$

where the circulation for each point is equal to $\frac{1}{2\delta N}$ and the initial positions of the vortices $X_0^i \forall i = 1, ..., N$ are uniformly distributed on the strip $[-1, 1] \times [-\delta, \delta]$. The randomly generated initial condition represents small perturbations in the system and is responsible for the different pattern formation.

The measure $\omega_t^N := \frac{1}{N} \sum_i \delta_{X_t^i}$ converges, in distribution, to the scalar vorticity field solution of the Euler equation (14); analogously with the continuous case, we see in Fig. 1a, c the development of instability in the form of macroscopic vortex-like structure on the boundary of the two fluid layers.

Note that the number of such macroscopic vortex-like structures and their position is entirely dependent on the initial condition: small perturbations on the randomly generated point vortices can produce entirely different macroscopic vortex-like structures, hence the instability of the two laminar fluids' profile.

### 4.1.2 The Role of Viscosity and Stability Restoration

The exact solution of the Navier-Stokes equation (14), with $v > 0$ and initial condition $u^0$, is given by

$$u(t, x_1, x_2) = (u_1(t, x_2), 0),$$

where $u_1(t, x_2)$ solves the heat equation,
$$ \partial_t u_1 = \nu \partial_{x_2}^2 u_1 \qquad (21) \\ u_1(0, x_2) = u_1^0(x_2). $$

Due to the spreading of the profile $u_1^0$, the solution becomes more stable; namely, the development of macroscopic vortex-like structures is delayed. In our numerical simulations, we reproduce this phenomenon by perturbing the system (8) through independent Brownian motions $B_t^i$, $i = 1, ..., N$, with variance linked to the viscosity parameter: $Var(B_t^i) \sim \sqrt{\nu}$. This system converges to the exact solution when $N \to \infty$; however, in our numerical study we are dealing with a finite system. For this reason, the profile of the strip remains quite stable for short times, with just a spread along the y-axes. We report in Fig. 2a, b a single configuration at two different

<figure><img src="image_21.png" /><figcaption>**Fig. 1** $\nu = 0$. (a): initial configuration, approximating a shear flow fluid dynamics; (b): iteration $t = 50$, formation of macroscopic vortex structures; (c): iteration $t = 100$, perfectly developed macroscopic vortex structures</figcaption></figure>
timesteps, $t = 50$ and $t = 100$; we take $\sqrt{\nu} = 0.095$ to better focus on the stability restoration. Our results are in agreement with the theory: from a comparison with Fig. 1b, c, we see that when $\nu > 0$, the profile is much more stable and diffused than in the deterministic case, and macroscopic vortex-like structures appear only at large times.

## 4.2 Numerical Results on Environmental Noise

### 4.2.1 Selection of Divergence Free Field

Starting with the same initial condition (17), we consider $N + M$ point vortices, each of which we associate with a position in $\mathbb{T}^2$ in the following way:

$$X_t^1, \dots, X_t^N, Y^1, \dots, Y^M.$$

Here, the vortices $Y^i, i = 1, \dots, M$ do not move, and when their interactions with point vortices became not negligible, they represent the feedback of small-scale turbulence acting on the fluid itself on large scales.

The new simulated vortex dynamics for $X_t^i$ in $\mathbb{T}^2$ as in (13), reads

<figure><img src="image_22.png" /><figcaption>**Fig. 2** $\nu > 0$. (a): iteration $t = 50$, preservation of strip profile; (b): iteration $t = 100$, first development of the strip profile's instability, formation of large rotating structures</figcaption></figure>
$$dX_t^i = \frac{1}{N} \sum_{i' \neq i} \Gamma_{i'} K (X_t^i - X_{t'}^{i'}) dt + \sum_j \sigma_j (X_t^i) \circ dW_t^j$$

with periodic boundary conditions, where $W_t^j$ are Brownian motions, all independent and uni-dimensional, and they are acting simultaneously on all the particles $i = 1, ..., N$. The environmental noise follows the Stratonovich integral prescription, automatically implemented in Heun's method [35].

We choose the divergence-free vector fields $\sigma_j$ as

$$\sigma_j(X_t^i) := a_j^{N,M} K (X_t^i - Y^j), \quad j = 1, ..., M$$

following the theoretical analysis performed in [22, 20]. Here, the intensities $a_j^{N,M}$ are linked to the scaling limit process, which produces a viscosity term on the large scales, and $K$, the Biot-Savart kernel, simulates the action of such small vortices. The idea behind such a selection is that we want to exploit the same features of the vortex model, with the formation of small-scale vortex structures generating feedback on the entire configuration. In the limit, the dynamics of such small structures, modulated through a Brownian motion, rebound on large scales, perturbing their motion with their dissipative properties and delaying the formation of the instability.

### 4.2.2 Positions and Intensities of Fixed Vortices

In order to connect with previous studies [20], we choose the positions of the fixed vortices $Y^j$ and their intensity $a_j^{N,M}$ according to the convergence of the scaling limit (4.2.1). More precisely, at each timestep, we generate $Y^j$, $j = 1, ..., M$, uniformly distributed point vortices; their position on the y-axis is apriori selected in the interval $[-\delta_{FX}, \delta_{FX}]$. In this setup, the vortices $Y^j$ are generated in a strip of variable height $2\delta_{FX}$; this strip contains the moving vortices $X_t^i$, and it is taken to be of the same height of our boundary fluid layers, or one order of magnitude greater. This choice emphasizes that our proposed “small-scale” structures should act on all points of the fluid in all directions: the average contribution of the $Y^j$ on the $X_t^i$ along every direction should mimic a Brownian motion. We explored different setups of positions and intensity; we selected meaningful realizations, as reported in Table 1.

<table><caption>**Table 1** Parameters of the discussed realizations</caption><thead><tr><th><i>M</i></th><th><i>δ</i><sub>FX</sub></th><th><i>m</i></th><th><i>a</i></th></tr></thead><tbody><tr><td>200000</td><td>0.1</td><td>0.0014</td><td>0.0005</td></tr><tr><td>132000</td><td>0.07</td><td>0.0014</td><td>0.0005</td></tr><tr><td>1000</td><td>0.07</td><td>0.0017</td><td>0.005</td></tr></tbody></table>
We choose the intensity of the “small-scale” perturbations following heuristic considerations. We consider the mean inter-particle distance between two fixed point vortices, $<r> := \sqrt{1/m}$, where $m := M/A$ is the particle density and $A$ is the total area occupied by the $M$ vortices.

Then, we focus on a single moving vortex, $X_t^i$; we compute the magnitude of its velocity when $X_t^i$ is at a distance $d = <r>/2$ from the nearest fixed vortex, i.e., its position is halfway between two fixed vortices. As a result, we obtain the following estimate for the velocity of $X_t^i$:

$$ \sum_j a_j^{N,M} \frac{1}{4\pi} \frac{|X_t^i - Y^j|^{\perp}}{\|X_t^i - Y^j\|^2} \sim \frac{1}{4\pi} \sum_j \frac{a_j^{N,M}}{d} \sim \frac{a^N}{4\pi} \sum_j \frac{1}{d}. $$

where we suppose that the coefficients depend on the configuration $(X_t^i)_i$, and are equal for each $j = 1, ..., M$. What is left is to estimate the number of fixed vortices such that the interaction with $X_t^i$ is not negligible: let us call this number $K$, giving us $\sim \frac{Ka^N}{4\pi d}$.

Thus, using the theory from the scaling limit of environmental transport noise (see e.g. [22, 20, 23]) and the construction of Sect. 3.2 for point vortices with transport noise, we assume that

$$ v \sim \frac{1}{2} \left( \frac{Ka^N}{4\pi d} \right)^2 . $$

This leads us to the estimate for the intensity of the fixed vortices

$$ a^N \sim 2\sqrt{2\pi} \frac{\sqrt{A}}{\sqrt{M}} \frac{\sqrt{v}}{K} $$

It remains to estimate $K$, the number of the nearest fixed vortices: consider a ball centered in $X_t^i$ with radius $d$, so that the area is $A_{near} = \pi d^2$. We recall that $m$ is the density of the fixed vortices, then the nearest vortices are:

$$ m \times A_{near} = \frac{M}{A} A_{near} = \frac{\pi}{4}. $$

Taking into account only the nearest vortices, we are underestimating the actual contribution of all the vortices. In particular, we should compute such contribution by considering a radius dependent on the effective range of the image of the Biot-Savart kernel. In fact, being $K \sim \frac{\Gamma_i}{|x|^2}$, contribution for distant particle is negligible. For this reason, we empirically selected a wider radius $\alpha d$, with $\alpha \sim 3$. Concluding we get our estimates for the intensity
$$a^N \sim \frac{8\sqrt{2}}{3} \sqrt{\nu} \frac{\sqrt{A}}{\sqrt{N}}$$

### 4.2.3 Effect of Small Scale Common Noise

As recalled in the previous paragraph's heuristics, the procedure of the scaling limit is preserved when both $N$ and $M$ are large, and the intensity $a_j^{N,M}$ is small. For this reason, we do not search for the same exact solution of the Navier-Stokes equation (14), $\nu > 0$, with initial condition $u^0$, as per the case of the independent noise. However, since the regime tends, in the limit, to the same solution, we expect a diffusive effect on the strip of the point vortex. More precisely, we expect to see a delay in the formation of macroscopic structures and a more dispersed displacement of such small vortex blobs that, on average, should maintain the strip configuration for a longer time.

In the first of our simulations, we generate at each time step $M \sim 2 \cdot 10^5$ fixed vortices, with intensity $a_j^{N,M} \sim 5 \cdot 10^{-4}$ which follows from the heuristics. The fixed vortices are uniformly distributed in a strip $[-1, 1] \times [-0.1, 0.1]$, containing the initial point vortices configuration. This particular setup captures the feedback effect of small scales on large vorticity structures, as the contribution of all the low-intensity perturbations on the dynamics of the point vortices averages in every direction. In the proposed numerical simulation, we show that the transport noise model reproduces the desired instability delay, even if it is slightly less effective than in the independent noise case; we illustrate a snapshot of a configuration for time $t = 50$ and $t = 100$ in Fig. 3a, b.

By comparison with Figs. 1b and 2a, we see that the initial strip configuration is preserved for a longer time than in the deterministic case and rotation of the fluid is milder, but the profile is less stable than in the viscous regime. If we focus on the deterministic case, we see blob-like structures formation already at $t = 50$; in contrast, in the transport noise regime, those structures are less visible and appear more prominently only at the end of our simulation ($t = 100$). This delay of the instability is evident in the realizations in Fig. 3b, compared with Figs. 1c and 2b: we notice a more diffused and homogeneous profile and a delayed formation of rotational structures due to the noise spreading the particles along the y-axis. A difference with the viscous case is that the compression in the x-axis is stronger than in the case of the independent noise, resulting in a more prominent stretch, which could resemble more the deterministic formations, placing the transport noise as a midpoint between the two regimes.

In the second of our simulations, the strip of fixed points $Y^j$ is generated in the same region as the point vortices at each timestep; we selected $M \sim 1.32 \cdot 10^5$, the fixed vortices are uniformly distributed in $[-1, 1] \times [-\delta - \epsilon, \delta + \epsilon]$, with $\epsilon = 0.03$, and their intensity is derived from the heuristics $a_j^{N,M} \sim 5 \cdot 10^{-4}$. The results are shown in Fig. 4: diffusion on the y-axis is present for short times, and preservation of the strip profile is guaranteed. However, the drawback of such a configuration
<figure><img src="image_23.png" /><figcaption>**Fig. 3** Environmental noise. (a): iteration *t* = 50, diffusive behaviour of strip profile; (b): iteration *t* = 100, degradation of profile, formation of macroscopic structures due to the stretch</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 4** Environmental noise case, iteration *t* = 50, diffusion of the strip is present for a short time with preserved configuration</figcaption></figure>

<figure><img src="image_25.png" /></figure>

is that analysis can be performed only for a short time: boundary effects of the fixed vortex strip can deteriorate the configuration, making the results unrealistic. In future works, we expect to overcome this obstacle by proposing a new method, now in the study, to generate small vortices only in regions activated by the shear flow's movement.

We performed a final simulation, in which we take the density of fixed vortices to be smaller than the density of the point vortices. In particular, the strip of fixed points $Y^j$ is generated in the same region as the point vortices at each timestep; we selected $M \sim 10^3$, the fixed vortices are uniformly distributed in $[-1, 1] \times [-\delta - \epsilon, \delta + \epsilon]$, with $\epsilon = 0.05$, and their intensity derived from the heuristics $a_j^{N,M} \sim 5 \cdot 10^{-3}$. While the diffusive behaviour is lost, as shown in Fig. 5, the strip is already broken at time $t = 50$, showing rotating structures. The fixed vortices' low density and higher intensity seem to produce new formations and medium-scale structures, showing
<figure><img src="image_26.png" /><figcaption>**Fig. 5** Environmental noise case, iteration $t = 50$, low density ratio between fixed vortices and point vortices showing emergence of medium-scale structures</figcaption></figure>

a completely different behaviour than the deterministic and viscous counterparts. For this reason, we need to investigate further the link between the ratio of vortices densities and the formation of new independent medium structures.

## 4.3 Diagnostics

In this section, we perform statistical analysis on the three configurations proposed in this study to highlight the differences and the reconstructed stability, or the emergence of new structures, in the Kelvin-Helmholtz instability problem.

As a first step, we compute the vorticity $\omega_e$ obtained through the same mollification as in Majda [4], through the vortex blob method applied to each of the point vortices $X_t^i$. We report our results for the vorticity computed in the deterministic case at $t = 100$ in Fig. 6a. We see that the vorticity measure concentration near the fluid's boundary layer is located in the newly developed macroscopic structures. Moreover, a displacement from the initial configuration where the laminar fluid started its evolution is present.

In contrast, the vortex blob solution with viscosity $\nu > 0$ retains its structure for longer times than in the inviscid case. The instability delay is graphically evident both from the configuration reported in Fig. 2a and the vorticity intensity reported in Fig. 6b: the vorticity measure concentration at time $t = 100$ is similar to the one of the initial strip but with a more diffused profile on the horizontal line.

From configuration Fig. 3b, finally, in Fig. 6c, we show the vorticity in the environmental noise regime at $t = 100$. In contrast with the inviscid case, we see no development of macroscopic structures in the profile. However, even though a more diffused profile, with lower density overall, is present, the stability of the strip profile is lost at larger times compared to Figs. 6b and 2b. This instability at larger times suggests that different behavior, dependent on the density of fixed vortex $Y^j$ and selection of transport noise fields $\sigma_j$, could arise in applying this kind of small-scale approximation. For this reason, we focus our analysis on small-time behaviour.

Concerning the formation of large rotating structures, we see that the particles spread in the horizontal direction when a forcing term, either an independent or transport noise, acts on the fluid, in contrast to the solution of the Euler equation
<figure><img src="image_27.png" /><figcaption>**Fig. 6** Vorticity $\omega_{\epsilon}$ at $t = 100$ in the case (a) inviscid, (b) viscous and (c) transport noise, showing macroscopic structures formation or delay of instability</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 7** Histograms of x-positions and empirical density at $t = 100$ in the case (a) inviscid, (b) viscous and (c) transport noise, showing formation of macroscopic structures</figcaption></figure>

with $v = 0$. In particular, we focus on the empirical density obtained from the x-axis in the three configurations at $t = 100$, Fig. 7. The deterministic case of Fig. 7a shows a complete formation of separate blobs with peaks in the exact locations of
the macroscopic structures, as in Fig. 1c. On the contrary, in the viscous Fig. 7b and transport noise case of Fig. 7c, the distribution of the vortices is more uniform, and it delays the instability of the fluid layers.

Following theoretical results, we know that with our initial condition $ω_0$, the velocity $u_t$ solves (21) when $v > 0$. This is crucial to understanding our system's short-time behaviour and the preservation of the initial configuration. This result states that the empirical density obtained from the y-position, when viscosity is present, maintains a Gaussian profile through time. In particular, in the case of $v = 0$, the viscosity follows a classic Euler equation. As such, the y-position profile is far from a Gaussian: it behaves like a multi-modal distribution, concentrated in the proximity of the center of the large structures. This result is supported both by the profile of the particle system and by Figs. 8a, and 9a; the qq-plot shows a distinct behaviour for small quantiles. The Kolmogorv-Smirnov test estimates a D-statistic of 0.031, with a p-value less than $10^{-9}$ confirming the rejection of the Gaussianity hypothesis.

Vice versa, when viscosity is present, i.e. $v > 0$, as in the case of independent Brownian motions, the noise's diffusive behaviour allows the profile's restoration in the y-direction: the strip configuration is preserved for a longer time. From the profile of the point vortex system and Figs. 8b, and 9b, we see that the empirical density approximates well the one of a Gaussian kernel. Moreover, the qq-plot suggests a perfect match with a normal distribution, suggesting the preservation of the strip through time, trading it with more spread on the y-axis. Finally, performing a Kolmogorv-Smirnov test, we see, in fact, a D-statistic of 0.005 and a p-value greater than 0.9 suggesting to accept the normality hypothesis. This behaviour is

<figure><img src="image_29.png" /><figcaption>**Fig. 8** Histograms of y-positions and empirical density at $t = 50$ in the case (a) inviscid, (b) viscous and (c) transport noise, showing different density profiles</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Fig. 9** qq-plot of the empirical density in the y-positions, at $t = 50$ in the case (a) inviscid, (b) viscous and (c) transport noise</figcaption></figure>

preserved throughout the simulation, degrading only at longer times when a few large formations start to rise, as in Fig. 2b.

In the transport noise case, we know that we recover the same viscous Euler equation solved in the independent Brownian motion case by applying a particular scaling limit procedure. To this end, we expected that the more stable profile shown in Fig. 3a presents the same diffusion on the density of the y-position as the case of the independent noise. This is the case as reported in Figs. 8c, and 9c: the behaviour at a short times, $t = 50$, in which the profile still approximates a Gaussian kernel. We perform a qq-plot and a Kolmogorov-Smirnov test on the y-position; we obtain a D-statistic of 0.011, and a p-value of 0.312; those results suggest the profile stability and the validity of our hypothesis.

However, later on ($t = 100$), even though the quantity obtained from the KS test is still preserved as in the viscous case, with a D-statistics of 0.008 and a p-value of 0.48, the profile degrades as shown in Fig. 3b. Those results show that the strip configuration is not preserved for longer times due to transport noise stretching
acting on the point vortices. This can be seen in the tail of the distribution in Fig. 9c, compared to the viscous case in Fig. 9b, which shows at $t = 50$ already a different behaviour. This suggests that the environmental noise's effect is not only related to the diffusivity of the strip, but is also responsible for the stretching and formation of different structures.

## 5 Concluding Remarks

In the present work, we have produced numerical simulations of 2D incompressible fluids, also perturbed by transport noise, using the point vortex method. We focused on the special case of shear flow formation that produced a Kelvin-Helmholtz instability, in order to test the dissipativity properties of small-space-scale transport noise. We confronted the intrinsic instability generated in the deterministic case with the possible recovery of the stability through injected noise in the system in the form of transport noise. We showed that, for short times, with a degree less intense than the viscous case, we can maintain the stability of the strip at the expense of a more small-scale irregularity and diffusion of the profile.

**Acknowledgments** We greatly thank anonymous referees for several improvements suggested in their reports.

The research of the first author is funded by the European Union (ERC, NoisyFluid, No. 101053472). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

## References

[1] C. R. Anderson and C. Greengard, On vortex methods, SIAM J. Numer. Anal. 22, 413–440 (1986).

[2] G. K. Batchelor, An Introduction to Fluid Dynamics Cambridge University Press, Cambridge, (1987).

[3] Bauer, W., Chandramouli, P., Chapron, B., Li, L. and Mémin, E., Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach, Journal of Physical Oceanography, 50 (4), 983–1003, 2020

[4] J. T. Beale and A. Majda, Higher order accurate vortex methods with explicit velocity kernels, J. Comput. Phys. 58, 188–208 (1985).

[5] L. C. Berselli, T. Iliescu, W. J. Layton. Mathematics of Large Eddy Simulation of Turbulent Flows, Springer, Berlin (2005).

[6] J. Boussinesq, Essai sur la théorie des eaux courantes, Mémoires présentés par divers savants à l'Académie des Sciences XXIII (1877), 1–680.

[7] Bresch, D., and Jabin, P.E., and Wang, Z. “On mean-field limits and quantitative estimates with a large class of singular kernels: Application to the Patlak–Keller–Segel model.” Comptes Rendus Mathematique, Volume 357, Issue 9, (2019).

[8] B. Chapron, D. Crisan, D. Holm, E. Mémin, A. Radomska Editors, Stochastic Transport in Upper Ocean Dynamics, Mathematics of Planet Earth 10, Springer, 2023.
[9] B. Chapron, E. Mémin, V. Resseguier, Geophysical flows under location uncertainty Geo-physical & Astrophysical Fluid Dynamics, 11(3), 149-176 (2017).

[10] Chorin, A.J., Numerical study of slightly viscous flow, J.Fluid Mech.,57,p.783, 1973

[11] R. R. Clements and D. J. Maull, The representation of sheets of vorticity by discrete vortices, Prog. Aerospace Sci. 16, 129-146 (1975).

[12] M. Coghi, F. Flandoli, Propagation of chaos for interacting particles subject to environmental noise, Ann. Appl. Probab. 26 (2016).

[13] C.J. Cotter, G.A. Gottwald, D.D. Holm, Stochastic partial differential fluid equations as a diffusive limit of deterministic lagrangian multi-time dynamics, Proc. R. Soc. A 473 (2017), 20170388.

[14] P. F. Cummins, G. Holloway. Reynolds stress and eddy viscosity in direct numerical simulations of sheared two-dimensional turbulence. J. Fluid Mech. 657 (2010), 394-412.

[15] Durr, D., and M. Pulvirenti. "On the vortex flow in bounded domains." Communications in Mathematical Physics 85.2 (1982): 265-273.

[16] F. Flandoli, L. Galeati, D. Luo, Scaling limit of stochastic 2D Euler equations with transport noises to the deterministic Navier-Stokes equations, J. Evol. Equ. 21 (2021), no. 1, 567-600.

[17] Flandoli, Franco, Lucio Galeati, and Dejun Luo. "Eddy heat exchange at the boundary under white noise turbulence." Philosophical Transactions of the Royal Society A 380.2219 (2022): 20210096.

[18] F. Flandoli, M. Gubinelli, E. Priola, Full well-posedness of point vortex dynamics corresponding to stochastic 2D Euler equations. Stochastic Process. Appl. 121 (2011), no. 7, 1445-1463.

[19] F. Flandoli, R. Huang, A. Papini, Turbulence enhancement of coagulation: the role of eddy diffusion in velocity, arXiv:2209.14387.

[20] F. Flandoli, D. Luo, Mean field limit of point vortices with environmental noises to deterministic 2D Navier-Stokes equations, arXiv:2101.06934.

[21] F. Flandoli, D. Luo, E. Luongo, 2D Smagorinsky type large eddy models as limits of stochastic PDEs, arXiv:2302.13614.

[22] F. Flandoli, E. Luongo, Stochastic Partial Differential Equations in Fluid Mechanics, Springer, to appear.

[23] L. Galeati, On the convergence of stochastic transport equations to a deterministic parabolic one, Stoch. Partial Differ. Equ. Anal. Comput. 8 (2020), no. 4, 833-868.

[24] Grotto, F. "Essential self-adjointness of Liouville operator for 2D Euler point vortices." Journal of Functional Analysis 279, Issue 6, (2020)

[25] Grotto F., Luongo E. and Maurelli M. Uniform Approximation of 2D Navier-Stokes Equations with Vorticity Creation by Stochastic Interacting Particle Systems, 2023, arXiv, math.AP, https://doi.org/10.48550/arXiv.2212.12753

[26] Grotto, F., Pappalettera, U. Burst of Point Vortices and Non-uniqueness of 2D Euler Equations. Arch Rational Mech Anal 245, 89-125 (2022).

[27] Grotto, F., Pappalettera, U. Correction to: Burst of Point Vortices and Non-uniqueness of 2D Euler Equations. Arch Rational Mech Anal 246, 139-140 (2022).

[28] S. K. Harouna, E. M'emin, Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling. Computers & Fluids, 156 (2017), 456-469.

[29] D. D. Holm, Variational principles for stochastic fluid dynamics, Proc. R. Soc. A. 471 (2015), 20140963.

[30] A. Labovsky, W. Layton. Magnetohydrodynamic flows: Boussinesq conjecture. J. Math. Anal. Appl. 434 (2016), 1665-1675.

[31] L. D. Landau, E. M. Lifshitz. Fluid Mechanics: Volume 6. Vol. 6. Elsevier, (1987).

[32] N. Jiang, W. Layton, M. McLaughlin, Y. Rong, H. Zhao. On the foundations of eddy viscosity models of turbulence. Fluids 2020, 5, 167.

[33] Jourdain, B., Meleard, S. "Probabilistic Interpretation and Particle Method for vortex equations with Neumann's boundary condition." Proceedings of the Edinburgh Mathematical Society, 47(3), 597-624. (2004).

[34] Kraichnan, R. (1974). Convection of a passive scalar by a quasi-uniform random straining field. Journal of Fluid Mechanics, 64(4), 737-762.
[35] R. Mannella, Integration Of Stochastic Differential Equations On A Computer, International Journal of Modern Physics, C 2002 13:09, 1177-1194.

[36] C. Marchioro, M. Pulvirenti, Vortex Methods in Two-Dimensional Fluid Dynamics, Lecture Notes in Physics (LNP, volume 203), 1984

[37] C. Marchioro, M. Pulvirenti, Mathematical Theory of Incompressible Nonviscous Fluids, Applied Mathematical Sciences, 96. Springer-Verlag, New York, 1994.

[38] Méléard, S., Monte-Carlo Approximations for 2D Navier-Stokes Equations with Measure Initial Data, Probability Theory and Related Fields, vol. 121, no. 3, pp. 367-388, 2001.

[39] E. Mémin, Fluid flow dynamics under location uncertainty, Geophysical & Astrophysical Fluid Dynamics, 108 (2014), no. 2, 119-146.

[40] Papini, A., Franco F., and Ruojun H. "Turbulence enhancement of coagulation: The role of eddy diffusion in velocity." Physica D: Nonlinear Phenomena 448 (2023): 133726.

[41] F. G. Schmitt, About Boussinesq's turbulent viscosity hypothesis: historical remarks and a direct evaluation of its validity, C. R. Mecanique 2007.

[42] Schochet, S. (1996), The point-vortex method for periodic weak solutions of the 2-D Euler equations. Comm. Pure Appl. Math., 49: 911-965.

[43] Subramaniam, S., A New Mesh-free Vortex Method, A New Mesh-free Vortex Method, Florida State University, (1996).

[44] A. Wirth, S. Gama, U. Frisch. Eddy viscosity of three-dimensional flows. J. Fluid Mech. 288 (1995), 249-264.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_31.png" /></figure>
# On the 3D Navier-Stokes Equations with Stochastic Lie Transport

<figure><img src="image_32.png" /></figure>

Check for
updates

Daniel Goodair and Dan Crisan

**Abstract** We prove the existence and uniqueness of maximal solutions to the 3D SALT (Stochastic Advection by Lie Transport) Navier-Stokes Equation in velocity and vorticity form, on the torus and the bounded domain respectively. In particular we demonstrate the efficacy of Goodair et al. (Existence and Uniqueness of Maximal Solutions to SPDEs with Applications to Viscous Fluid Equations, 2023. Stochastics and Partial Differential Equations: Analysis and Computations, pp.1-64) in showing the well-posedness for both the velocity and vorticity form of the equation, as well as obtaining the first analytically strong existence result for a fluid equation perturbed by Lie transport noise on a bounded domain.

## 1 Introduction

The theoretical analysis of Stochastic Navier-Stokes Equations dates back to the work of Bensoussan and Temam [4] in 1973, where the problem of existence of solutions is addressed in the presence of a random forcing term. The well-posedness question for additive and multiplicative noise has since seen significant developments, for example through the works [1, 7, 24, 29, 39, 42, 43] and references therein. The interest in this problem has expanded into analytical properties of these solutions, particularly along the lines of ergodicity, which can be seen in [17, 18, 22, 26, 27, 34]. In the present work our concern is the Navier-Stokes Equations with Stochastic Lie Transport, derived through the principle of Stochastic Advection by Lie Transport (SALT) introduced in [35]. We consider the equation

$$u_t - u_0 + \int_0^t \mathcal{L}_{u_s} u_s ds - \nu \int_0^t \Delta u_s ds + \int_0^t B u_s \circ dW_s + \nabla \rho_t = 0 \quad (1)$$

<footnote>D. Goodair (✉) · D. Crisan
Imperial College London, London, UK
e-mail: djg116@ic.ac.uk</footnote>
where $u$ represents the fluid velocity, $\rho$ the pressure,<fnref n="1" /> $W$ is a Cylindrical Brownian Motion, $\mathcal{L}$ represents the nonlinear term and $B$ is a first order differential operator (the SALT Operator) formally addressed in Sect. 2.3. Intrinsic to this stochastic methodology is that $B$ is defined relative to a collection of functions $(\xi_i)$ which physically represent spatial correlations. These $(\xi_i)$ can be determined at coarse-grain resolutions from finely resolved numerical simulations, and mathematically are derived as eigenvectors of a velocity-velocity correlation matrix (see [10, 11, 12]). We pose the equation (1) in 3 dimensions and impose the divergence free constraint on $u$. We shall consider the problem both over the torus $\mathbb{T}^3$ and a smooth bounded domain $\mathcal{O} \subset \mathbb{R}^3$. In the case of the torus we supplement the equation with the zero-average condition (as is classical), whilst for the bounded domain we impose the boundary condition

$$u \cdot n = 0, \quad w = 0 \qquad (2)$$

where $n$ represents the outwards unit normal at the boundary, and $w$ the fluid vorticity. These are the so called Lions boundary conditions, considered in [38] and shown to be a particular case of the Navier boundary conditions in [36] (note that this is done in 2D, whilst a treatment of the Navier boundary conditions in 3D can be found in [28]). The significance of such a boundary condition is well documented in that work by Kelliher, and can be seen in other works such as [28] where the boundary layer is explicitly addressed. The precise mathematical interpretation of these conditions, and the operators of (1), are explicated in Sect. 2.2. A complete derivation of this equation can be found in [45].

This work continues the theoretical development of fluid models perturbed by a transport type noise, the significance of which was posed as early as 1992 in the paper [6]. The area has garnered substantial attention in more recent years with the seminal works [35, 41], in which the authors establish a new class of stochastic equations driven by transport type noise which serve as fluid dynamics models by adding uncertainty in the transport of fluid parcels to reflect the unresolved scales. This paper partners that of [32] where we showed the existence and uniqueness of maximal solutions to Stochastic Partial Differential Equations satisfying an abstract framework, built to cope with a general transport type noise as we see in (1). The importance of such equations in modelling, numerical schemes and data assimilation is reviewed there, along with the theoretical developments of these equations: let us briefly mention some interesting results [3, 5, 21, 23]. We only draw particular attention here to the Navier-Stokes Equations, and results on a bounded domain. The Navier-Stokes Equations have been studied with transport type noise, for example in the works [14, 19, 25, 43], though typically solutions are analytically weak and where strong solutions are considered major concessions in the noise are made. In these cases a cancellation property is evident in the noise

<footnote n="1">¹ The pressure term is a semimartingale, and an explicit form for the SALT Euler Equation is given in [45] Subsection 3.3.</footnote>
term, so the resulting energy balance is formally the same as for the deterministic
equation without noise. These difficulties have been addressed on the torus, in the
likes of the papers [13, 37] and those further addressed in [32], but extending a
control of this noise term to a bounded domain remains open. Furthermore, whilst
the presence of viscosity does improve the solution theory, it invites additional
challenges in controlling the noise. Energy methods require non-standard Sobolev
inner products to conduct the required integration by parts for the viscous term in the
bounded domain, so we must provide novel estimates on the transport type noise in
these inner products. The problem of analytically strong solutions to fluid equations
perturbed by a transport type noise in the bounded domain has been considered in
[8], though the authors assume that the gradient dependency is of a small enough
size to be directly controlled and that the noise terms are traceless under Leray
Projection; such assumptions are designed to circumvent the technical difficulties
of a first order noise operator on a bounded domain, and is a luxury that the SALT
equations do not have.

The goal of this paper is to apply the abstract framework established in [32],
providing a rigorous justification of the results first announced in [31] and extending
them to the vorticity form of Eq. (1) on a bounded domain. The purpose of this is
twofold:

1. To demonstrate the efficacy of the criteria from [32]. This is most pertinent in our treatment of the velocity form on the torus; whilst the well-posedness results that we present are new, one would expect them to hold. This belief is compounded by the well-posedness results for both viscous and inviscid fluid equations perturbed by Lie transport noise on the torus [2, 13, 15, 16, 37] for which the Navier-Stokes equations do not present any intrinsic additional difficulty, though still contributing a necessary piece to the somewhat sparse literature. We hope to convey that verifying the assumptions of [32] is now the simplest and most efficient way to prove the well-posedness of this equation and, by extension, associated stochastic viscous fluid equations.

2. To obtain the well-posedness of the equation in the presence of a boundary, which does present intrinsic additional difficulty beyond the existing results.

In the interests of brevity we provide a shorter manuscript here, though greater
detail can be found in the former arXiv version [33]. In Sect. 2 we establish
the stochastic and functional framework necessary to understand (1), along with
fundamental properties of the operators involved. In Sect. 3 we make precise how
equation (1) fits into the framework of [32], as a problem posed on the torus $T^3$.
In Sect. 4 we consider the vorticity form of equation (1) as a problem posed on a
bounded domain of $\mathbb{R}^3$. We again justify the assumptions in [32] to prove existence
and uniqueness of maximal solutions to this equation. Additional details for the
proofs are given in Sect. 5, along with the results of the partnering paper [32].
# 2 Preliminaries

## 2.1 Elementary Notation

In the following $\mathcal{O}$ can represent both the 3-dimensional torus $\mathbb{T}^3$ and a smooth bounded domain $\mathcal{O} \subset \mathbb{R}^3$. We consider Banach Spaces as measure spaces equipped with the Borel $\sigma$-algebra, and use $\lambda$ to represent the Lebesgue Measure. Let $(\mathcal{X}, \mu)$ denote a general measure space, $(\mathcal{Y}, \|\cdot\|_{\mathcal{Y}})$ and $(\mathcal{Z}, \|\cdot\|_{\mathcal{Z}})$ be Banach Spaces, and $(\mathcal{U}, \langle \cdot, \cdot \rangle_{\mathcal{U}})$, $(\mathcal{H}, \langle \cdot, \cdot \rangle_{\mathcal{H}})$ be general Hilbert spaces. $\mathcal{O}$ is equipped with Euclidean norm.

* $L^p(\mathcal{X}; \mathcal{Y})$ is the class of measurable $p$-integrable functions from $\mathcal{X}$ into $\mathcal{Y}$, $1 \le p < \infty$, which is a Banach space with norm

$$ \|\phi\|_{L^p(\mathcal{X}; \mathcal{Y})}^p := \int_{\mathcal{X}} \|\phi(x)\|_{\mathcal{Y}}^p \mu(dx). $$

In particular $L^2(\mathcal{X}; \mathcal{Y})$ is a Hilbert Space when $\mathcal{Y}$ itself is Hilbert, with the standard inner product

$$ \langle \phi, \psi \rangle_{L^2(\mathcal{X}; \mathcal{Y})} = \int_{\mathcal{X}} \langle \phi(x), \psi(x) \rangle_{\mathcal{Y}} \mu(dx). $$

In the case $\mathcal{X} = \mathcal{O}$ and $\mathcal{Y} = \mathbb{R}^3$ note that

$$ \|\phi\|_{L^2(\mathcal{O}; \mathbb{R}^3)}^2 = \sum_{l=1}^{3} \|\phi^l\|_{L^2(\mathcal{O}; \mathbb{R})}^2, \quad \phi = (\phi^1, \dots, \phi^3), \quad \phi^l : \mathcal{O} \to \mathbb{R}. $$

We denote $\|\cdot\|_{L^p(\mathcal{O}; \mathbb{R}^3)}$ by $\|\cdot\|_{L^p}$ and $\|\cdot\|_{L^2(\mathcal{O}; \mathbb{R}^3)}$ by $\|\cdot\|_{\mathcal{O}}$.

* $L^\infty(\mathcal{X}; \mathcal{Y})$ is the class of measurable functions from $\mathcal{X}$ into $\mathcal{Y}$ which are essentially bounded, which is a Banach Space when equipped with the norm

$$ \|\phi\|_{L^\infty(\mathcal{X}; \mathcal{Y})} := \inf\{C \ge 0 : \|\phi(x)\|_{\mathcal{Y}} \le C \text{ for } \mu\text{-a.e. } x \in \mathcal{X}\}. $$

* $L^\infty(\mathcal{O}; \mathbb{R}^3)$ is the class of measurable functions from $\mathcal{O}$ into $\mathbb{R}^3$ such that $\phi^l \in L^\infty(\mathcal{O}; \mathbb{R})$ for $l = 1, \dots, N$, which is a Banach Space when equipped with the norm

$$ \|\phi\|_{L^\infty} := \sup_{l \le N} \|\phi^l\|_{L^\infty(\mathcal{O}; \mathbb{R})}. $$

* $C(\mathcal{X}; \mathcal{Y})$ is the space of continuous functions from $\mathcal{X}$ into $\mathcal{Y}$.

* $C^m(\mathcal{O}; \mathbb{R})$ is the space of $m \in \mathbb{N}$ times continuously differentiable functions from $\mathcal{O}$ to $\mathbb{R}$, that is $\phi \in C^m(\mathcal{O}; \mathbb{R})$ if and only if for every $N$ dimensional multi index
$α = α_1, ..., α_N$ with $|α| ≤ m$, $D^αφ ∈ C(O; R)$ where $D^α$ is the corresponding classical derivative operator $\partial_{X_1}^{α_1} ... \partial_{X_N}^{α_N}$.

* $C^∞(O; R)$ is the intersection over all $m ∈ N$ of the spaces $C^m(O; R)$.

* $C_0^m(O; R)$ for $m ∈ N$ or $m = ∞$ is the subspace of $C^m(O; R)$ of functions which have compact support.

* $C^m(O; \mathbb{R}^3)$, $C_0^m(O; \mathbb{R}^3)$ for $m ∈ N$ or $m = ∞$ is the space of functions from $O, O$ to $\mathbb{R}^3$ whose $N$ component mappings each belong to $C^m(O; \mathbb{R})$, $C_0^m(O; \mathbb{R})$.

* $W^{m,p}(O; R)$ for $1 ≤ p < ∞$ is the sub-class of $L^p(O, R)$ which has all weak derivatives up to order $m ∈ N$ also of class $L^p(O, R)$. This is a Banach space with norm

$$ \|\phi\|_{W^{m,p}(O, \mathbb{R})}^p := \sum_{|\alpha| \le m} \|D^\alpha \phi\|_{L^p(O, \mathbb{R})}^p $$

where $D^α$ is the corresponding weak derivative operator. In the case $p = 2$ the space $W^{m,2}(O, \mathbb{R})$ is Hilbert with inner product

$$ \langle \phi, \psi \rangle_{W^{m,2}(O; \mathbb{R})} := \sum_{|\alpha| \le m} \langle D^\alpha \phi, D^\alpha \psi \rangle_{L^2(O; \mathbb{R})}. $$

* $W^{m,∞}(O; R)$ for $m ∈ N$ is the sub-class of $L^∞(O, R)$ which has all weak derivatives up to order $m ∈ N$ also of class $L^∞(O, R)$. This is a Banach space with norm

$$ \|\phi\|_{W^{m,∞}(O, \mathbb{R})} := \sup_{|\alpha| \le m} \|D^\alpha \phi\|_{L^∞(O; \mathbb{R}^3)}. $$

* $W^{m,p}(O; \mathbb{R}^3)$ for $1 ≤ p < ∞$ is the sub-class of $L^p(O, \mathbb{R}^3)$ which has all weak derivatives up to order $m ∈ N$ also of class $L^p(O, \mathbb{R}^3)$. This is a Banach space with norm

$$ \|\phi\|_{W^{m,p}}^p := \sum_{l=1}^3 \|\phi^l\|_{W^{m,p}(O; \mathbb{R})}^p. $$

In the case $p = 2$ the space $W^{m,2}(O, \mathbb{R}^3)$ is Hilbertian with inner product

$$ \langle \phi, \psi \rangle_{W^{m,2}} := \sum_{l=1}^{3} \langle \phi^l, \psi^l \rangle_{W^{m,2}(O; \mathbb{R})}. $$

* $W^{m,∞}(O; \mathbb{R}^3)$ is the sub-class of $L^∞(O, \mathbb{R}^3)$ which has all weak derivatives up to order $m ∈ N$ also of class $L^∞(O, \mathbb{R}^3)$. This is a Banach space with norm
$$ \|\phi\|_{W^{m,\infty}(O, \mathbb{R}^3)} := \sup_{l \le N} \|\phi^l\|_{W^{m,\infty}(O; \mathbb{R})}. $$

• $\dot{L}^2(\mathbb{T}^3; \mathbb{R}^3)$ is the subset of $L^2(\mathbb{T}^3; \mathbb{R}^3)$ of functions $\phi$ such that

$$ \int_{\mathbb{T}^3} \phi \, d\lambda = 0. $$

• $\dot{W}^{m,2}(\mathbb{T}^3; \mathbb{R}^3)$ is simply the intersection $W^{m,2}(\mathbb{T}^3; \mathbb{R}^3) \cap \dot{L}^2(\mathbb{T}^3; \mathbb{R}^3)$.

• $W_0^{m,p}(\mathcal{O}; \mathbb{R})$, $W_0^{m,p}(\mathcal{O}; \mathbb{R}^3)$ for $m \in N$ and $1 \le p \le \infty$ is the closure of $C_0^\infty(\mathcal{O}; \mathbb{R})$, $C_0^\infty(\mathcal{O}; \mathbb{R}^3)$ in $W^{m,p}(\mathcal{O}; \mathbb{R})$, $W^{m,p}(\mathcal{O}; \mathbb{R}^3)$.

• $\mathcal{L}(\mathcal{Y}; \mathcal{Z})$ is the space of bounded linear operators from $\mathcal{Y}$ to $\mathcal{Z}$. This is a Banach Space when equipped with the norm

$$ \|F\|_{\mathcal{L}(\mathcal{Y}; \mathcal{Z})} = \sup_{\|y\|_{\mathcal{Y}}=1} \|Fy\|_{\mathcal{Z}} $$

and is simply the dual space $\mathcal{Y}^*$ when $\mathcal{Z} = \mathbb{R}$, with operator norm $\|\cdot\|_{\mathcal{Y}^*}$.

• $\mathcal{L}^2(\mathcal{U}; \mathcal{H})$ is the space of Hilbert-Schmidt operators from $\mathcal{U}$ to $\mathcal{H}$, defined as the elements $F \in \mathcal{L}(\mathcal{U}; \mathcal{H})$ such that for some basis $(e_i)$ of $\mathcal{U}$,

$$ \sum_{i=1}^{\infty} \|Fe_i\|_{\mathcal{H}}^2 < \infty. $$

This is a Hilbert space with inner product

$$ \langle F, G \rangle_{\mathcal{L}^2(\mathcal{U}; \mathcal{H})} = \sum_{i=1}^{\infty} \langle Fe_i, Ge_i \rangle_{\mathcal{H}} $$

which is independent of the choice of basis.

We will consider a partial ordering on the 3-dimensional multi-indices by $\alpha \le \beta$ if and only if for all $l = 1, 2, 3$ we have that $\alpha_l \le \beta_l$. We extend this to notation $<$ by $\alpha < \beta$ if and only if $\alpha \le \beta$ and for some $l = 1, 2, 3$, $\alpha_l < \beta_l$.

We also now introduce some less familiar spaces in slightly greater detail. We recall notation that $\mathcal{O}$ represents a smooth bounded domain in $\mathbb{R}^3$ which we now fix, $\mathbb{T}^3$ is the 3-dimensional torus, and $\mathcal{O}$ freely denotes both $\mathbb{T}^3$ and $\mathcal{O}$.

**Definition 2.1** We define $C_{0,\sigma}^\infty(\mathcal{O}; \mathbb{R}^3)$ as the subset of $C_0^\infty(\mathcal{O}; \mathbb{R}^3)$ of functions which are divergence-free. $L_\sigma^2(\mathcal{O}; \mathbb{R}^3)$ is defined as the completion of $C_{0,\sigma}^\infty(\mathcal{O}; \mathbb{R}^3)$ in $L^2(\mathcal{O}; \mathbb{R}^3)$, whilst we introduce $W_\sigma^{1,2}(\mathcal{O}; \mathbb{R}^3)$ as the intersection of $W_0^{1,2}(\mathcal{O}; \mathbb{R}^3)$ with $L_\sigma^2(\mathcal{O}; \mathbb{R}^3)$ and $W_\sigma^{2,2}(\mathcal{O}; \mathbb{R}^3)$ as the intersection of $W^{2,2}(\mathcal{O}; \mathbb{R}^3)$ with $W_\sigma^{1,2}(\mathcal{O}; \mathbb{R}^3)$.

Recall that any function $f \in L^2(\mathbb{T}^3; \mathbb{R}^3)$ admits the representation
$$f(x) = \sum_{k \in \mathbb{Z}^3} f_k e^{ik \cdot x} \qquad (3)$$

whereby each $f_k \in \mathbb{C}^3$ is such that $f_k = \overline{f}_{-k}$ and the infinite sum is defined as a limit in $L^2(\mathbb{T}^3; \mathbb{R}^3)$, see e.g. [44] Subsection 1.5 for details.

**Definition 2.2** We define $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$ as the subset of $\dot{L}^2(\mathbb{T}^3; \mathbb{R}^3)$ of functions $f$ whereby for all $k \in \mathbb{Z}^3$, $k \cdot f_k = 0$ with $f_k$ as in (3). For general $m \in \mathbb{N}$ we introduce $W_{\sigma}^{m,2}(\mathbb{T}^3; \mathbb{R}^3)$ as the intersection of $W^{m,2}(\mathbb{T}^3; \mathbb{R}^3)$ respectively with $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$.

Note that $W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$ is precisely the subspace of $W^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$ consisting of zero-average divergence free functions. Similarly $W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^3)$ is precisely the subspace of $W_0^{1,2}(\mathcal{O}; \mathbb{R}^3)$ consisting of divergence free functions. Moreover, $W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^3)$ is the completion of $C_{0,\sigma}^{\infty}(\mathcal{O}; \mathbb{R}^3)$ in $W^{1,2}(\mathcal{O}; \mathbb{R}^3)$. The general space $W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^3)$ thus incorporates the divergence-free and zero-average/zero-trace condition.

As for the stochastic set up, let $(\Omega, \mathcal{F}, (\mathcal{F}_t), \mathbb{P})$ be a fixed filtered probability space satisfying the usual conditions of completeness and right continuity. We take $\mathcal{W}$ to be a cylindrical Brownian Motion over some Hilbert Space $\mathcal{U}$ with orthonormal basis $(e_i)$. Recall ([30], Subsection 1.4) that $\mathcal{W}$ admits the representation $\mathcal{W}_t = \sum_{i=1}^{\infty} e_i W_t^i$ as a limit in $L^2(\Omega; \mathcal{U}')$ whereby the $(W^i)$ are a collection of i.i.d. standard real valued Brownian Motions and $\mathcal{U}'$ is an enlargement of the Hilbert Space $\mathcal{U}$ such that the embedding $J: \mathcal{U} \to \mathcal{U}'$ is Hilbert-Schmidt and $\mathcal{W}$ is a $JJ^*$-cylindrical Brownian Motion over $\mathcal{U}'$. Given a process $F: [0, T] \times \Omega \to \mathcal{L}^2(\mathcal{U}; \mathcal{H})$ progressively measurable and such that $F \in L^2(\Omega \times [0, T]; \mathcal{L}^2(\mathcal{U}; \mathcal{H}))$, for any $0 \le t \le T$ we define the stochastic integral

$$\int_0^t F_s d\mathcal{W}_s := \sum_{i=1}^{\infty} \int_0^t F_s(e_i) dW_s^i$$

where the infinite sum is taken in $L^2(\Omega; \mathcal{H})$. We can extend this notion to processes $F$ which are such that $F(\omega) \in L^2([0, T]; \mathcal{L}^2(\mathcal{U}; \mathcal{H}))$ for $\mathbb{P}-\text{a.e. } \omega$ via the traditional localisation procedure. In this case the stochastic integral is a local martingale in $\mathcal{H}^2$.

<footnote n="2">² A complete, direct construction of this integral, a treatment of its properties and the fundamentals of stochastic calculus in infinite dimensions can be found in [30] Section 1.</footnote>
## 2.2 Functional Framework

We now recap the classical functional framework for the study of the deterministic Navier-Stokes Equation. Firstly we briefly comment on the pressure term $\nabla \rho$, which will not play any role in our analysis. $\rho$ does not come with an evolution equation and is simply chosen to ensure the incompressibility (divergence-free) condition; moreover we will ignore this term via a suitable projection (in Sect. 3 we even consider a different form of the equation) and treat the projected equation, with the understanding to append a pressure to it later. This procedure is well discussed in [44] Sect. 5 and [40], and an explicit form for the pressure for the SALT Euler Equation is given in [45] Subsection 3.3.

The mapping $\mathcal{L}$ is defined for sufficiently regular functions $f, g : \mathcal{O} \to \mathbb{R}^3$ by $\mathcal{L}_f g := \sum_{j=1}^3 f^j \partial_j g$. Here and throughout the text we make no notational distinction between differential operators acting on a vector valued function or a scalar valued one; that is, we understand $\partial_j g$ by its component mappings $(\partial_l g)^l := \partial_j g^l$. We now give some clarification as to ‘sufficiently regular’, by stating basic properties of this mapping. For any $m \ge 1$, the mapping $\mathcal{L} : W^{m+1,2}(\mathcal{O}; \mathbb{R}^3) \to W^{m,2}(\mathcal{O}; \mathbb{R}^3)$ defined by $f \mapsto \mathcal{L}_f f$ is continuous. Additionally there exists a constant $c$ such that for any $f, g \in W^{k,2}(\mathcal{O}; \mathbb{R}^3)$ for $k \in \mathbb{N}$ as appropriate, we have the bounds:

$$ \|\mathcal{L}_f g\| + \|\mathcal{L}_g f\| \le c \|g\|_{W^{1,2}} \|f\|_{W^{2,2}}; \quad (4) $$

$$ \|\mathcal{L}_g f\|_{W^{1,2}} \le c \|g\|_{W^{1,2}} \|f\|_{W^{3,2}}; \quad (5) $$

$$ \|\mathcal{L}_g f\|_{W^{1,2}} \le c \|g\|_{W^{2,2}} \|f\|_{W^{2,2}}. \quad (6) $$

We introduce the Leray Projector $\mathcal{P}$ as the orthogonal projection in $L^2(\mathcal{O}; \mathbb{R}^3)$ onto $L_\sigma^2(\mathcal{O}; \mathbb{R}^3)$. It is well known (see e.g. [47] Remark 1.6.) that for any $m \in \mathbb{N}$, $\mathcal{P}$ is continuous as a mapping $\mathcal{P} : W^{m,2}(\mathcal{O}; \mathbb{R}^3) \to W^{m,2}(\mathcal{O}; \mathbb{R}^3)$. In fact, the complement space of $L_\sigma^2(\mathcal{O}; \mathbb{R}^3)$ can be characterised (this is the so called Helmholtz-Weyl decomposition), a result that we state explicitly as we will need to exploit the precise structure in future arguments.

**Lemma 2.3** *Define the space*

$$ L_{\sigma}^{2,\perp}(\mathcal{O}; \mathbb{R}^3) := \{\psi \in L^2(\mathcal{O}; \mathbb{R}^3) : \psi = \nabla g \text{ for some } g \in W^{1,2}(\mathcal{O}; \mathbb{R})\}. $$

Then indeed $L_\sigma^{2,\perp}(\mathcal{O}; \mathbb{R}^3)$ is orthogonal to $L_\sigma^2(\mathcal{O}; \mathbb{R}^3)$ in $L^2(\mathcal{O}; \mathbb{R}^3)$, i.e. for any $\phi \in L_\sigma^2(\mathcal{O}; \mathbb{R}^3)$ and $\psi \in L_\sigma^{2,\perp}(\mathcal{O}; \mathbb{R}^3)$ we have that

$$ \langle \phi, \psi \rangle = 0. $$

Moreover, every $f \in L^2(\mathcal{O}; \mathbb{R}^3)$ has the unique decomposition

$$ f = \phi + \psi \quad (7) $$
for some $\phi \in L_{\sigma}^{2}(\mathcal{O}; \mathbb{R}^{3})$, $\psi \in L_{\sigma}^{2, \perp}(\mathcal{O}; \mathbb{R}^{3})$ and every $f \in L^{2}(\mathbb{T}^{3}; \mathbb{R}^{3})$ has the unique decomposition

$$f = \phi + \psi + c \qquad (8)$$

where $\phi \in L_{\sigma}^{2}(\mathbb{T}^{3}; \mathbb{R}^{3})$, $\psi \in L_{\sigma}^{2, \perp}(\mathbb{T}^{3}; \mathbb{R}^{3})$ and $c$ is a constant function: that is, there exists $k \in \mathbb{R}^{3}$ such that each component mapping $c^{j}$ is identically equal to $k^{j}$, $j = 1, 2, 3$.

**Proof** See [47] Theorems 1.4, 1.5 and [44] Theorem 2.6.

**Corollary 2.3.1** Every $f \in L^{2}(\mathcal{O}; \mathbb{R}^{3})$ admits the representation

$$f = \mathcal{P}f + \nabla g \qquad (9)$$

for some $g \in W^{1,2}(\mathcal{O}; \mathbb{R})$. Similarly every $f \in L^{2}(\mathbb{T}^{3}; \mathbb{R}^{3})$ admits the representation

$$f = \mathcal{P}f + \nabla g + c \qquad (10)$$

for some $g \in W^{1,2}(\mathbb{T}^{3}; \mathbb{R})$ and constant function $c$.

**Proof** It is an immediate property of the orthogonal projection that $\mathcal{P}f$ is the unique element $\phi \in L_{\sigma}^{2}(\mathcal{O}; \mathbb{R}^{3})$ of (7) and (8).

Through $\mathcal{P}$ we define the Stokes Operator $A : W^{2,2}(\mathcal{O}; \mathbb{R}^{3}) \to L_{\sigma}^{2}(\mathcal{O}; \mathbb{R}^{3})$ by $A := -\mathcal{P}\Delta$. Once more we understand the Laplacian as an operator on vector valued functions through the component mappings, $(\Delta f)^{l} := \Delta f^{l}$. From the continuity of $\mathcal{P}$ we have immediately that for $m \in \{0\} \cup \mathbb{N}$, $A : W^{m+2,2}(\mathcal{O}; \mathbb{R}^{3}) \to W^{m,2}(\mathcal{O}; \mathbb{R}^{3})$ is continuous. We remark that $\Delta$ leaves the complement space $L_{\sigma}^{2, \perp}(\mathcal{O}; \mathbb{R}^{3})$ invariant, so $\mathcal{AP}$ is equal to $A$ on $W^{2,2}(\mathcal{O}; \mathbb{R}^{3})$. Moreover (see [44] Theorem 2.24) there exists a collection of functions $(a_{k})$, $a_{k} \in W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^{3}) \cap C^{\infty}(\bar{\mathcal{O}}; \mathbb{R}^{3})$ such that the $(a_{k})$ are eigenfunctions of $A$, are an orthonormal basis in $L_{\sigma}^{2}(\mathcal{O}; \mathbb{R}^{3})$ and an orthogonal basis in $W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^{3})$. The eigenvalues $(\lambda_{k})$ are strictly positive and approach infinity as $k \to \infty$. Therefore every $f \in L_{\sigma}^{2}(\mathcal{O}; \mathbb{R}^{3})$ admits the representation

$$f = \sum_{k=1}^{\infty} f_{k}a_{k} \qquad (11)$$

where $f_{k} = \langle f, a_{k} \rangle$, as a limit in $L^{2}(\mathcal{O}; \mathbb{R}^{3})$.

**Definition 2.4** For $m \in \mathbb{N}$ we introduce the spaces $D(A^{m/2})$ as the subspaces of $L_{\sigma}^{2}(\mathcal{O}; \mathbb{R}^{3})$ consisting of functions $f$ such that
$$ \sum_{k=1}^{\infty} \lambda_k^m f_k^2 < \infty $$

for $f_k$ as in (11). Then $A^{m/2} : D(A^{m/2}) \to L_{\sigma}^2(O; \mathbb{R}^3)$ is defined by

$$ A^{m/2} : f \mapsto \sum_{k=1}^{\infty} \lambda_k^{m/2} f_k a_k. $$

We present some fundamental properties regarding these spaces, which are justified in [9] Proposition 4.12, [44] Exercises 2.12, 2.13 and the discussion in Subsection 2.3.

1. $D(A^{m/2}) \subset W^{m,2}(O; \mathbb{R}^3) \cap W_{\sigma}^{1,2}(O; \mathbb{R}^3)$ and the bilinear form

$$ \langle f, g \rangle_m := \langle A^{m/2} f, A^{m/2} g \rangle $$

is an inner product on $D(A^{m/2})$;

2. For $m$ even the induced norm $\|\cdot\|_m^2 = \langle \cdot, \cdot \rangle_m$ is equivalent to the $W^{m,2}(O; \mathbb{R}^3)$ norm, and for $m$ odd there is a constant $c$ such that

$$ \|\cdot\|_{W^{m,2}} \le c \|\cdot\|_m; $$

3. $D(A) = W_{\sigma}^{2,2}(O; \mathbb{R}^3)$ and $D(A^{1/2}) = W_{\sigma}^{1,2}(O; \mathbb{R}^3)$ with the additional property that $\|\cdot\|_1$ is equivalent to $\|\cdot\|_{W^{1,2}}$ on this space.

It can be directly shown that for any $p, q \in \mathbb{N}$ with $p \le q$, $p + q = 2m$ and $f \in D(A^{m/2})$, $g \in D(A^{q/2})$ we have that

$$ \langle f, g \rangle_m = \langle A^{p/2} f, A^{q/2} g \rangle. \quad (12) $$

From here we can also see that the collection of functions $(a_k)$ form an orthogonal basis of $W_{\sigma}^{1,2}(O; \mathbb{R}^3)$ equipped with the $\langle \cdot, \cdot \rangle_1$ inner product. In addition to using these spaces defined by powers of the Stokes Operator, we also use the basis to consider finite dimensional approximations of these spaces.

**Definition 2.5** We define $\mathcal{P}_n$ as the orthogonal projection onto $\text{span}\{a_1, \dots, a_n\}$ in $L^2(O; \mathbb{R}^3)$. That is $\mathcal{P}_n$ is given by

$$ \mathcal{P}_n : f \mapsto \sum_{k=1}^{n} \langle f, a_k \rangle a_k $$

for $f \in L^2(O; \mathbb{R}^3)$.

The restriction of $\mathcal{P}_n$ to $D(A^{m/2})$ is self-adjoint for the $\langle \cdot, \cdot \rangle_m$ inner product, and there exists a constant $c$ independent of $n$ such that for all $f \in D(A^{m/2})$,
$$ \|P_n f\|_{W^{m,2}} \le c \|f\|_{W^{m,2}}, \qquad (13) $$

see [44] Lemma 4.1 for details. Similar ideas justify that for all $f \in W_{\sigma}^{1,2}(O; \mathbb{R}^3)$, $g \in W_{\sigma}^{2,2}(O; \mathbb{R}^3)$,

$$ \|(I - P_n) f\|^2 \le \frac{1}{\lambda_n} \|f\|_1^2, \qquad \|(I - P_n)g\|_1^2 \le \frac{1}{\lambda_n} \|g\|_2^2 $$

where $I$ represents the identity operator in the relevant spaces. To conclude this subsection we discuss briefly bounds related to the nonlinear term, which will be used in our analysis. For every $\phi \in W_{\sigma}^{1,2}(O; \mathbb{R}^3)$ and $f, g \in W^{1,2}(O; \mathbb{R}^3)$, we have that

$$ \langle \mathcal{L}_{\phi} f, g \rangle = -\langle f, \mathcal{L}_{\phi} g \rangle \qquad (14) $$

$$ \langle \mathcal{L}_{\phi} f, f \rangle = 0. \qquad (15) $$

## 2.3 The SALT Operator

Having established the relevant function spaces and some fundamental properties of the operators involved in the deterministic Navier-Stokes Equation, we now address the operator $B$ appearing in the Stratonovich integral of (1). As in [30] Subsection 2.2, the operator $B$ is defined by its action on the basis vectors $(e_i)$ of $\mathfrak{U}$. We shall show in Sect. 3.3 that $B$ does indeed satisfy Assumption 2.2.2 of [30] for the spaces to $V$, $H$, $U$, $X$ to be defined. With the notation of [30], each $B_i$ is defined relative to the correlations $\xi_i$ for sufficiently regular $f$ by the mapping

$$ B_i : f \mapsto \mathcal{L}_{\xi_i} f + \mathcal{T}_{\xi_i} f $$

where $\mathcal{L}$ is as before, and $\mathcal{T}$ is a new operator that we introduce defined by

$$ \mathcal{T}_g f := \sum_{j=1}^{3} f^j \nabla g^j. $$

We shall assume throughout that each $\xi_i$ belongs to the space $W_{\sigma}^{1,2}(O; \mathbb{R}^3)$. If for some fixed $m \in \mathbb{N}$ we have $\xi_i \in W^{m+2,\infty}(O; \mathbb{R}^3)$ then for all $k = 0, \dots, m+1$,

$$ \| \mathcal{T}_{\xi_i} f \|_{W^{k,2}}^2 \le c \| \xi_i \|_{W^{k+1,\infty}}^2 \| f \|_{W^{k,2}}^2 \qquad (16) $$

$$ \| \mathcal{L}_{\xi_i} f \|_{W^{k,2}}^2 \le c \| \xi_i \|_{W^{k,\infty}}^2 \| f \|_{W^{k+1,2}}^2 \qquad (17) $$

$$ \| B_i f \|_{W^{k,2}}^2 \le c \| \xi_i \|_{W^{k+1,\infty}}^2 \| f \|_{W^{k+1,2}}^2. \qquad (18) $$
Moreover $\mathcal{T}_{\xi_i}$ is a bounded linear operator on $L^2(\mathcal{O}; \mathbb{R}^3)$ so has adjoint $\mathcal{T}_{\xi_i}^*$ satisfying the same boundedness. In conjunction with property (14), $\mathcal{L}_{\xi_i}$ is a densely defined operator in $L^2(\mathcal{O}; \mathbb{R}^3)$ with domain of definition $W^{1,2}(\mathcal{O}; \mathbb{R}^3)$, and has adjoint $\mathcal{L}_{\xi_i}^*$ in this space given by $-\mathcal{L}_{\xi_i}$ with same dense domain of definition. Likewise then $B_i^*$ is the densely defined adjoint $-\mathcal{L}_{\xi_i} + \mathcal{T}_{\xi_i}^*$. Our techniques centre around energy estimates, where the key idea as to how we preserve these estimates in the case of a transport type noise owes to the following proposition.

**Proposition 2.6** There exists a constant $c$ such that for each $i$ and for all $f \in W^{k+2,2}(\mathcal{O}; \mathbb{R}^3)$ with $k = 0, \dots, m$, we have the bounds

$$ \langle B_i^2 f, f \rangle_{W^{k,2}} + \|B_i f\|_{W^{k,2}}^2 \le c \|\xi_i\|_{W^{k+2,\infty}}^2 \|f\|_{W^{k,2}}^2, \quad (19) $$

$$ \langle B_i f, f \rangle_{W^{k,2}}^2 \le c \|\xi_i\|_{W^{k+1,\infty}}^2 \|f\|_{W^{k,2}}^4. \quad (20) $$

**Proof** See Sect. 5.1. $\square$

<figure><img src="image_33.png" /></figure>

Another valuable result is given now, which will be necessary in showing comparable estimates to Proposition 2.6 in the $\langle \cdot, \cdot \rangle_k$ inner product for appropriate $k$.

**Lemma 2.7** We have that

$$ B_i : L_{\sigma}^{2,\perp}(\mathcal{O}; \mathbb{R}^3) \cap W^{1,2}(\mathcal{O}; \mathbb{R}^3) \to L_{\sigma}^{2,\perp}(\mathcal{O}; \mathbb{R}^3) $$

and moreover that $\mathcal{P}B_i = \mathcal{P}B_i\mathcal{P}$ on $W^{1,2}(\mathcal{O}; \mathbb{R}^3)$.

**Proof** See Sect. 5.1. $\square$

<figure><img src="image_34.png" /></figure>

We note that this result holds true only in the presence of the additional $\mathcal{T}_{\xi_i}$ term in the operator, highlighting the significance of considering a noise which is not purely transport. The Leray Projector does pose difficulties in the presence of a boundary though, which we state here.

**Remark 2** The Leray Projector does not preserve the space $W_0^{1,2}(\mathcal{O}; \mathbb{R}^3)$, and so we cannot say that $\mathcal{P}B_i : W_{\sigma}^{2,2}(\mathcal{O}; \mathbb{R}^3) \to W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^3)$. The issues arising from this operator not satisfying the zero-trace property are fundamentally why we only treat the Torus for the velocity form in Sect. 3.

# 3 The Velocity Equation on the Torus

In this section we restrict ourselves to the Torus $\mathbb{T}^3$, leaving a treatment of the bounded domain to Sect. 4. We also now fix our assumptions on the $(\xi_i)$, assuming that each $\xi_i \in W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3) \cap W^{3,\infty}(\mathbb{T}^3; \mathbb{R}^3)$ and they collectively satisfy
$$ \sum_{i=1}^{\infty} \|\xi_i\|_{W^{3,\infty}}^2 < \infty. \qquad (21) $$

## 3.1 Definitions and Results

Here we state the key definitions and results of this section. To facilitate our analysis we work with an equation projected by the Leray Projector as discussed at the start of Sect. 2.2. Thus we consider the new equation

$$ u_t - u_0 + \int_0^t \mathcal{P}L_{u_s} u_s ds + v \int_0^t Au_s ds + \int_0^t \mathcal{P}Bu_s \circ dW_s = 0 \qquad (22) $$

obtained at a heuristic level by projecting all terms of (1). Having not defined solutions of (1) we cannot be too formal here, but the idea is that we require solutions in $L_{\sigma}^2(O; \mathbb{R}^3)$ with initial condition also in this space so they are invariant under $\mathcal{P}$, and $\mathcal{P}$ is a bounded linear operator so can be taken through the integrals (see [30] Corollary 1.6.12.1 for this result in Itô integration, understanding the Stratonovich integral as the sum of an Itô and a Bochner integral).

**Theorem 3.1** For any given $\mathcal{F}_0$-measurable $u_0 : \Omega \to W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3)$ there exists a pair $(u, \tau)$ such that: $\tau$ is a $\mathbb{P}$-a.s. positive stopping time and $u$ is a process whereby for $\mathbb{P}$-a.e. $\omega$, $u.(\omega) \in C([0, T]; W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3))$ and $u.(\omega)1_{\le \tau(\omega)} \in L^2([0, T]; W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3))$ for all $T > 0$ with $u.1_{\le \tau}$ progressively measurable in $W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$, and moreover satisfying the identity

$$ u_t = u_0 - \int_0^{t \wedge \tau} \mathcal{P}L_{u_s} u_s ds - v \int_0^{t \wedge \tau} Au_s ds - \int_0^{t \wedge \tau} \mathcal{P}Bu_s \circ dW_s $$

$\mathbb{P}$-a.s. in $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$ for all $t \ge 0$.

Theorem 3.1 will be proved as a consequence of the following proposition.

**Proposition 3.2** Suppose that $(u, \tau)$ are such that: $\tau$ is a $\mathbb{P}$-a.s. positive stopping time and $u$ is a process whereby for $\mathbb{P}$-a.e. $\omega$, $u.(\omega) \in C([0, T]; W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3))$ and $u.(\omega)1_{\le \tau(\omega)} \in L^2([0, T]; W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3))$ for all $T > 0$ with $u.1_{\le \tau}$ progressively measurable in $W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$, and moreover satisfying the identity

$$ u_t = u_0 - \int_0^{t \wedge \tau} \mathcal{P}L_{u_s} u_s ds - v \int_0^{t \wedge \tau} Au_s ds \\ + \frac{1}{2} \int_0^{t \wedge \tau} \sum_{i=1}^{\infty} \mathcal{P}B_i^2 u_s ds - \int_0^{t \wedge \tau} \mathcal{P}Bu_s dW_s $$
P - a.s. in $W_{\sigma}^{1,2}(\mathcal{O}; \mathbb{R}^3)$ for all $t \ge 0$. Then the pair $(u, \tau)$ satisfies the identity

$$u_t = u_0 - \int_0^{t \wedge \tau} \mathcal{P}L_{u_s} u_s \, ds - \nu \int_0^{t \wedge \tau} A u_s \, ds - \int_0^{t \wedge \tau} \mathcal{P}B u_s \circ dW_s$$

P - a.s. in $L_{\sigma}^2(\mathcal{O}; \mathbb{R}^3)$ for all $t \ge 0$.

Proposition 3.2 motivates studying the converted equation

$$u_t = u_0 - \int_0^t \mathcal{P}L_{u_s} u_s \, ds - \nu \int_0^t A u_s \, ds + \frac{1}{2} \int_0^t \sum_{i=1}^{\infty} \mathcal{P}B_i^2 u_s \, ds - \int_0^t \mathcal{P}B u_s \, dW_s. \quad (23)$$

Once we convert to the Itô Form though, starting with an initial condition in $W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3)$ is not optimal in the sense that, at least according to the deterministic theory, we should be able to construct a solution (satisfying the identity in $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$ as is natural) for only a $W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$ initial condition. To this end we give the following definitions and the main result of this section. Definition 3.3 is stated for an arbitrary $\mathcal{F}_0$-measurable $u_0: \Omega \to W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$.

**Definition 3.3** A pair $(u, \tau)$ where $\tau$ is a $P$-a.s. positive stopping time and $u$ is a process such that for $P$-a.e. $\omega$, $u.(\omega) \in C([0, T]; W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3))$ and $u.(\omega)1_{\le \tau(\omega)} \in L^2([0, T]; W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3))$ for all $T > 0$ with $u.1_{\le \tau}$ progressively measurable in $W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3)$, is said to be a local strong solution of the equation (23) if the identity

$$u_t = u_0 - \int_0^{t \wedge \tau} \mathcal{P}L_{u_s} u_s \, ds - \nu \int_0^{t \wedge \tau} A u_s \, ds \\ + \frac{1}{2} \int_0^{t \wedge \tau} \sum_{i=1}^{\infty} \mathcal{P}B_i^2 u_s \, ds - \int_0^{t \wedge \tau} \mathcal{P}B u_s \, dW_s \quad (24)$$

holds $P$-a.s. in $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$ for all $t \ge 0$.

**Remark 3** If $(u, \tau)$ is a local strong solution of the equation (23), then $u. = u_{\wedge \tau}$. The time integrals in (24) are well defined as Bochner integrals in $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$, and the Itô integral in $W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$. These integrals are shown to be well defined in the abstract case in [30] Definition 2.2.3.

**Definition 3.4** A pair $(u, \Theta)$ such that there exists a sequence of stopping times $(\theta_j)$ which are $P$-a.s. monotone increasing and convergent to $\Theta$, whereby $(u_{\wedge \theta_j}, \theta_j)$ is a local strong solution of the equation (23) for each $j$, is said to be a maximal strong solution of the equation (23) if for any other pair $(v, \Gamma)$ with this property then $\Theta \le \Gamma P$-a.s. implies $\Theta = \Gamma P$-a.s..

**Definition 3.5** A maximal strong solution $(u, \Theta)$ of the equation (23) is said to be pathwise unique if for any other such solution $(v, \Gamma)$, then $\Theta = \Gamma P$-a.s. and for
all $t \in [0, \Theta)$,

$$ \mathbb{P}(\{\omega \in \Omega : u_t(\omega) = v_t(\omega)\}) = 1. $$

**Theorem 3.6** For any given $\mathcal{F}_0$- measurable $u_0 : \Omega \to W_\sigma^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$, there exists a pathwise unique maximal strong solution $(u, \Theta)$ of the equation (23). Moreover at $\mathbb{P}$ - a.e. $\omega$ for which $\Theta(\omega) < \infty$, we have that

$$ \sup_{r \in [0, \Theta(\omega))} \|u_r(\omega)\|_1^2 + \int_0^{\Theta(\omega)} \|u_r(\omega)\|_2^2 dr = \infty. \quad (25) $$

## 3.2 Operator Bounds

In this subsection we state some intermediary results regarding control on the operators involved. In the following and throughout this section $c$ will represent a generic constant changing from line to line, $c(\varepsilon)$ will be a generic constant dependent on a fixed $\varepsilon$, $f$ and $g$ will be arbitrary elements of $W_\sigma^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$ and $f_n \in \text{span}\{a_1, \dots, a_n\}$. The proofs of these lemmas can be found in Sect. 5.1.

**Lemma 3.7** For any $\varepsilon > 0$, we have that

$$ |\langle \mathcal{P}_n \mathcal{P} \mathcal{L}_{f_n}, f_n \rangle_2| \le c(\varepsilon) \|f_n\|_2^4 + \varepsilon \|f_n\|_3^2; \quad (26) $$

$$ \langle \mathcal{P}_n \mathcal{P} B_i^2 f_n, f_n \rangle_1 + \langle \mathcal{P}_n \mathcal{P} B_i f_n, \mathcal{P}_n \mathcal{P} B_i f_n \rangle_1 \le c(\varepsilon) \|\xi_i\|_{W^{3,\infty}}^2 \|f_n\|_1^2 \\ + \varepsilon \|\xi_i\|_{W^{3,\infty}}^2 \|f_n\|_2^2; \quad (27) $$

$$ \langle \mathcal{P}_n \mathcal{P} B_i^2 f_n, f_n \rangle_2 + \langle \mathcal{P}_n \mathcal{P} B_i f_n, \mathcal{P}_n \mathcal{P} B_i f_n \rangle_2 \le c(\varepsilon) \|\xi_i\|_{W^{3,\infty}}^2 \|f_n\|_2^2 \\ + \varepsilon \|\xi_i\|_{W^{3,\infty}}^2 \|f_n\|_3^2. \quad (28) $$

**Remark 4** The algebra property of $W^{k,2}(\mathbb{T}^3; \mathbb{R}^3)$ when $k = 2$ is fundamental in proving (26); a result of the form

$$ |\langle \mathcal{P}_n \mathcal{P} \mathcal{L}_{f_n}, f_n \rangle_1| \le c(\varepsilon) \|f_n\|_1^4 + \varepsilon \|f_n\|_2^2 $$

is unavailable to us given that the algebra property does not hold for $k = 1$. This is revisited in Remark 5.

**Lemma 3.8** For any $\varepsilon > 0$, we have that

$$ |\langle \mathcal{P} \mathcal{L}_f f - \mathcal{P} \mathcal{L}_g g, f - g \rangle_1| \le c(\varepsilon) (\|g\|_1^4 + \|f\|_2^2) \|f - g\|_1^2 + \varepsilon \|f - g\|_2^2. $$
**Lemma 3.9** For any $\epsilon > 0$, we have that

$$ \left| \langle \mathcal{P} \mathcal{L}_f f - \mathcal{P} \mathcal{L}_g g, f - g \rangle \right| \leq c(\epsilon) \|f\|_2^2 \|f - g\|^2 + \epsilon \|f - g\|_1^2. $$

## 3.3 Proof of Proposition 3.2

We prove this result through the abstract procedure used in [30] Subsections 2.2 and 2.3, the result of which is stated in Sect. 5.2. Towards this goal we define the quartet of spaces

$$ V := W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3), \quad H := W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3), $$

$$ U := W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3), \quad X := L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3). $$

We equip $L_{\sigma}^2(\mathbb{T}^3; \mathbb{R}^3)$ with the usual $\langle \cdot, \cdot \rangle$ inner product, but then equip $W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$ and $W_{\sigma}^{2,2}(\mathbb{T}^3; \mathbb{R}^3)$ with the $\langle \cdot, \cdot \rangle_1$ and $\langle \cdot, \cdot \rangle_2$ inner products respectively. In fact we also have that $D(A^{3/2}) = W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$ and that the $\langle \cdot, \cdot \rangle_3$ inner product is equivalent to the usual $\langle \cdot, \cdot \rangle_{W^{3,2}}$ one on this space (see [44] Theorem 2.27), so we endow $W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$ with $\langle \cdot, \cdot \rangle_3$. Our SPDE (22) takes the form of (68) for the operators

$$ Q := -(\mathcal{P}\mathcal{L} + \nu A), \quad G := -\mathcal{P}B $$

where we rewrite $(-\mathcal{P}B_i)^2$ as $\mathcal{P}B_i^2$ firstly from the linearity of $\mathcal{P}B_i$ to deal with the minus and secondly using the property that $\mathcal{P}B_i = \mathcal{P}B_i\mathcal{P}$. It is worth appreciating here that we chose to project the equation and then convert it into Itô Form, but we may equally have chosen to convert the unprojected Stratonovich Form and then project the resulting Itô Equation. Without addressing the conversion of the unprojected equation in complete detail, we would directly arrive at (23) taking this approach. Indeed before projection our correction term would be of the form $\sum_{i=1}^{\infty} B_i^2$ plus a function in $L_{\sigma}^{2\perp}(\mathbb{T}^3, \mathbb{R}^3)$ as defined in Lemma 2.7. Under projection this is $\mathcal{P}\sum_{i=1}^{\infty} B_i^2$ which is just $\sum_{i=1}^{\infty} \mathcal{P}B_i^2$ from the linearity and continuity. Therefore the property $\mathcal{P}B_i = \mathcal{P}B_i\mathcal{P}$ from Lemma 2.7 establishes the consistency between these approaches.

To prove the result we check the Assumptions 5.3 and 5.4. Starting with 5.3, we first of all have that $\nu A$ is continuous from $W_{\sigma}^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$ into $W_{\sigma}^{1,2}(\mathbb{T}^3; \mathbb{R}^3)$, therefore it is measurable and as a linear operator too satisfies the boundedness. As for $\mathcal{P}\mathcal{L}$, measurability is satisfied in the same way and for the boundedness we have that

$$ \|\mathcal{P}\mathcal{L}_f f\|_1 \leq c \|\mathcal{P}\mathcal{L}_f f\|_{W^{1,2}} \leq c \|\mathcal{L}_f f\|_{W^{1,2}} \leq c \|f\|_{W^{1,2}} \|f\|_{W^{3,2}} \leq c \|f\|_1 \|f\|_3 $$
for any $f \in W_g^{3,2}(\mathbb{T}^3; \mathbb{R}^3)$ where $c$ is a generic constant, critically applying (5). This verifies Assumption 5.3 so we move on to Assumption 5.4, which is immediate from (18) and the linearity of $\mathcal{P}B$ to show continuity in all relevant spaces.

## 3.4 Proofs of Theorems 3.1 and 3.6.

We use the abstract results of [32] stated in Sects. 5.3 and 5.4. Definitions 3.3, 3.4, 3.5 and Theorem 3.6 are precisely Definitions 5.19, 5.20, 5.21 and Theorem 5.22 for the equation (23) with respect to the spaces $V, H, U, X$ as defined in Sect. 3.3. Indeed we would also prove Theorem 3.1 through Proposition 3.2 by showing the existence of a local solution with the regularity of Definition 5.12. Therefore we prove both Theorem 3.1 and 3.6 by showing that the assumptions of Sects. 5.3 and 5.4, are satisfied. We work with the operators

$$ \mathcal{A} := -(\mathcal{P}\mathcal{L} + \nu A) + \frac{1}{2} \sum_{i=1}^{\infty} \mathcal{P}B_i^2, \quad \mathcal{G} := -\mathcal{P}B $$

which were addressed to be measurable mappings into the required spaces in Sect. 3.3. We now prove Theorem 3.1 by justifying the assumptions of Sect. 5.3.

**Proof of Theorem 3.1** First note that the density of the spaces is immediately inherited from the density of the usual Sobolev Spaces and the equivalence of the norms. The bilinear form satisfying (70) is chosen to be

$$ \langle f, g \rangle_{U \times V} := \langle A^{1/2} f, A^{3/2} g \rangle $$

which reduces to the $\langle \cdot, \cdot \rangle_2$ inner product from (12). The remainder of the proof is in treating the numbered assumptions.

**Assumption 5.6:** We use the system $(a_k)$ of eigenfunctions of the Stokes Operator.

**Assumption 5.7:** Once more (74) follows from the discussion in Sect. 3.3. For (75) we treat the different operators in $\mathcal{A}$ individually, starting from the nonlinear term:

$$ \begin{aligned} \|\mathcal{P}\mathcal{L}_f f - \mathcal{P}\mathcal{L}_g g\|_1 &= \|\mathcal{P}\mathcal{L}_f (f-g) + \mathcal{P}\mathcal{L}_{f-g} g\|_1 \\ &\le \|\mathcal{P}\mathcal{L}_f (f-g)\|_1 + \|\mathcal{P}\mathcal{L}_{f-g} g\|_1 \\ &\le c \|\mathcal{L}_f (f-g)\|_{W^{1,2}} + c \|\mathcal{L}_{f-g} g\|_{W^{1,2}} \\ &\le c \|f\|_{W^{1,2}} \|f-g\|_{W^{3,2}} + c \|f-g\|_{W^{1,2}} \|g\|_{W^{3,2}} \\ &\le c (\|f\|_{W^{1,2}} + \|g\|_{W^{3,2}}) \|f-g\|_{W^{3,2}} \end{aligned} $$