11institutetext:  Departament d’Astronomia i Astrofísica, Universitat de
València, C/ Dr. Moliner, 50, 46100, Burjassot, València, Spain
Observatori Astronòmic, Universitat de València, C/ Catedràtic José Beltrán 2,
46980, Paterna, València, Spain

# Shocks, clouds and atomic outflows in active galactic nuclei hosting
relativistic jets

Manel Perucho 11221122

(Received xx / Accepted yy)

###### Abstract

Context. A number of observations have revealed atomic and/or molecular lines
in active galaxies hosting jets and outflows. Line widths indicate outward
motions of hundreds to few thousands of kilometers per second. They appear
associated to the presence of radio emission in Gigahert-peaked spectrum (GPS)
and/or compact steep spectrum (CSS) sources, with linear sizes $\leq 10\,{\rm
kpc}$. Numerical simulations have shown that the bow shocks triggered by
relativistic jets in their host galaxies drive ionisation and turbulence in
the interstellar medium (ISM). However, the presence of atomic lines requires
rapid recombination of ionised gas, which seems to be hard to explain from the
physical conditions revealed so far by numerical simulations of powerful jets.

Aims. The aim of this paper is to provide a global frame to explain the
presence of lines in terms of jet and shock evolution, and fix the parameter
space in which the atomic and molecular outflows might occur.

Methods. This parameter space is inspired by numerical simulations and basic
analytical models of jet evolution as a background.

Results. Our results show that a plausible, general explanation involves
momentum transfer and heating to the interstellar medium gas by jet triggered
shocks within the inner kiloparsecs. The presence of post-shock atomic gas is
possible in the case of shocks interacting with dense clouds that remain
relatively stable after the shock passage.

Conclusions. According to our results, current numerical simulations cannot
reproduce the physical conditions to explain the presence of atomic and
molecular outflows in young radio-sources. However, I show that these outflows
might occur in low-power jets at all scales, and predict a trend towards
powerful jets showing lines at CSS scales, when clouds have cooled to
recombination temperatures.

###### Key Words.:

galaxies – individual – jets

##  1 Introduction

Relativistic jets in AGN are produced in the environment of the central black
holes. Our current understanding explains their formation in terms of magneto-
hydrodynamical processes, via extraction of rotational energy either from a
Kerr black-hole (Blandford & Znajek 1977), or from the inner accretion disk
(Blandford & Payne 1982). The jets are launched as sub-Alfvénic Poynting-flux
dominated and accelerated by magnetic fields (Vlahakis & Königl 2004;
Komissarov et al. 2007; Komissarov 2012) and internal energy (Perucho & Martí
2007, Ricci et al., submitted) to relativistic, superfast-magnetosonic speeds.
It has been shown that the jets not only contribute to feedback at large
scales and may participate in regulating star formation within their host
galaxies (see,e.g. Shabala & Alexander 2009; Raouf et al. 2019), but also may
play a significant role in driving atomic and molecular outflows within the
inner kiloparsecs from the AGN (e.g., Morganti & Oosterloo 2018; Santoro et
al. 2020; Morganti et al. 2021; Tadhunter et al. 2021; Kukreti et al. 2023).

VLBI observations of compact, young radio sources revealed advance velocities
of $\sim 0.2\,c$ (e.g., Owsianik & Conway 1998; Polatidis & Conway 2003). Such
velocities probably correspond to powerful jets. Numerical simulations of jet
propagation within the inner kiloparsec of active galaxies have shown that
they do expand at supersonic velocities and trigger strong shocks in the
interstellar medium (see, e.g., Wagner & Bicknell 2011; Wagner et al. 2012;
Mukherjee et al. 2016; Bicknell et al. 2018; Mukherjee et al. 2018; Perucho et
al. 2021; Meenakshi et al. 2022). Simulations also show that lateral expansion
is slower by a factor of 3-5, as derived from aspect ratios, so the radial
expansion velocity can be $\simeq 0.04-0.06\,c$. Such strong shocks increase
the temperature of the gas to values high enough to completely ionise hydrogen
atoms. Post-shock temperatures of $10^{8}-10^{10}$ K can be easily achieved in
this region, as we also justify in this work. In contrast, simulations of low
power jets ($L_{j}\simeq 10^{42}\,{\rm erg/s}$) present advance velocities
$\leq 10^{-3}\,c$, and temperature jumps $\leq 2$ (Perucho et al. 2014a) to
values $\sim 10^{7}$ K (in the hot ISM component).

However, radio observations of atomic lines show the presence of HI outflows
in compact radio sources, i.e., Gigahertz-peaked spectrum (GPS) sources (Holt
et al. 2011), with linear sizes $\leq 1\,{\rm kpc}$, and compact steep
spectrum (CSS) sources, $\leq 10\,{\rm kpc}$ in size (see, e.g. O’Dea et al.
2002; Labiano, A. et al. 2005; Shih et al. 2013; Reynaldi & Feinstein 2016;
Morganti et al. 2016; Zovaro et al. 2019; Schulz et al. 2021, and Morganti &
Oosterloo (2018) for a review). The presence of HI requires gas temperatures
$\leq 1.5\times 10^{5}\,{\rm K}$ in the case of collisional equilibrium, or
even smaller in the case of thermodynamical equilibrium, i.e., at least four
orders of magnitude below the shocked gas temperature in powerful jets (e.g.,
Perucho et al. 2021). The inferred velocities range from several hundreds to
few thousands of km/s. In Holt et al. (2011), the authors point out that the
line widths reveal larger velocities in two observed GPS galaxies than in
CSSs, and that they are compatible with an in situ shock-cloud collision,
which has also been reported by Couto et al. (2020) for detected OIII lines.
However, this is not obvious in all sources.

Although the aforementioned numerical simulations have tackled the problem of
jet evolution inside the host galaxy, the authors have been forced to focus
either on the inner kiloparsec (GPS region) or beyond (CSS region) by
computational limitations, so the transition has barely been discussed in
terms of heating and cooling of the ISM gas. Furthermore, none of the
simulations had the required resolution to study the evolution of shocked
clouds (see, however, Mandal et al., in preparation).

Bicknell et al. (1997, 2003) presented a model based on Begelman (1996) that
successfully explained a number of properties of GPS sources. The authors show
that radiative shocks can create an ionised shocked ISM that free-free absorbs
low frequency radiation and explain the evolution of the peak frequency with
size in these young radio-galaxies (Fanti et al. 1990; O’Dea & Baum 1997;
Stanghellini et al. 1998). Those theoretical papers focus on high power jets
($10^{45-46}\,{\rm erg/s}$).

This paper takes a simple analytical approach to tackle the discrepancy
between physical conditions derived from numerical simulations and those
required for the presence of HI outflows in host galaxies (see, e.g., Snow &
Hillier 2021a, b, and references therein for a detailed study of the physics
of two-fluid shocked plasmas). The model presented here is motivated by
results from numerical simulations of relativistic outflows (see, e.g.
Hardcastle 2018; Turner et al. 2023, for other analytical models also driven
by simulations and Turner & Shabala (2023) for a review) and focuses on the
role of shocks on the environment, extending the study to low power radio-
galaxies. The shocks are shown to be non radiative in the case of powerful
jets, but to become radiative at a few kpc for low power jets. The model shows
that atomic outflows are compatible with conditions produced by low power jets
at all scales, but require further propagation time and cloud stability to
allow for gas cooling in the case of powerful sources.

These results would point towards more powerful radio-sources showing atomic
lines when they have propagated longer distances, but low power radio-sources
presenting them at all scales, including smaller linear sizes. Altogether,
this conclusion sets a plausible general scenario to explain observations of
atomic outflows at all scales, depending on the jet power, either within or
beyond the inner kiloparsec.

The paper is structured as follows. In section 2, I present the physical
arguments that could explain the presence of atomic lines in active galaxies
based on shock-heating and estimated cooling times. In section 3, I discuss
the implications for both low and high power radio-galaxies, the caveats of
the simple analytic approach and the implications for our understanding of jet
physics and feedback at larger scales. Finally, in section 4 I summarise the
main conclusions of this paper.

##  2 Shocks: heating and cooling

###  2.1 The heating

We focus here on the lateral expansion of jet-triggered shocks (i.e.,
misaligned with the jet direction), which affects a large volume of the ISM.
The shock propagates into the two-phase ISM, with a hotter, dilute medium with
temperatures $\sim 10^{6}-10^{7}\,{\rm K}$ and number densities $\sim
0.1-1\,{\rm cm^{-3}}$, and a colder, denser component in clouds, with
temperatures $\sim 10-1000\,{\rm K}$ and densities $\sim 10^{2}-10^{3}\,{\rm
cm^{-3}}$. The shocked ISM medium surrounds the shocked jet gas (or cocoon),
with very high temperatures and lower densities (see, e.g., Wagner & Bicknell
2011; Mukherjee et al. 2016, 2018; Perucho et al. 2021). Nevertheless, our
interest is focused on the shocked ISM-ISM interaction region (see Krause &
Alexander 2007; Krause 2008, for the case of a cloud entering and interacting
with the jet cocoon).

The following expressions (jump conditions in the shock reference frame) give
the post-shock pressures and temperatures:

$\rho_{1}\,v_{1}\,=\,\rho_{2}\,v_{2},$ (1)
$v_{2}\,=\,\frac{P_{1}\,-\,P_{2}}{\rho_{1}\,v_{1}}\,+\,v_{1},$ (2)
$v_{1}^{2}\,=\,\frac{1}{2\,\rho_{1}}\left[(\Gamma+1)\,P_{2}+(\Gamma-1)\,P_{1}\right]$
(3)

with subscripts $1$ and $2$ standing for pre and post-shock states of the
interstellar medium (ISM), respectively, and where $P$ represents pressure,
$\rho$ the gas rest-mass density, $v$ the velocity in the shock propagation
direction, and $\Gamma$ the adiabatic exponent of an ideal gas.

In the reference frame of the shocked ISM, $v_{1}=v_{s}$, the shock velocity.
For an assumed ideal gas, we also have $P/\rho\,=\,k\,T/\mu m_{H}$, where $k$
is the Boltzmann constant, $\mu$ is the mean molecular weight, and $m_{H}$ is
the atomic hydrogen mass. The previous expressions allow us to derive the
shock advance velocity by imposing the ambient density and both pre and post-
shock pressures.

####  2.1.1 The hot phase

In the case of AGN bow-shocks, the post-shock pressure is rapidly homogenised
in the shocked hot phase because of the large sound speed achieved by the gas,
as it has been shown by numerical simulations. Perucho et al. (2017) also
showed that the following expression gives an order-of-magnitude estimate of
post-shock pressure in relativistic jets:

$P_{2}\,\leq\,\frac{L_{j}\,t}{V},$ (4)

where $L_{j}$ is the jet power, and $V$ is the shocked volume.111The
expression implicitly assumes that all injected power is converted into
internal energy driving the lobe expansion, which is valid for relativistic
(hot and fast) outflows. Thus, this represents an upper limit for pressure,
which propagates as an upper limit for the derived post-shock temperatures.
The ISM pressure can be taken to be the usually estimated value
$P_{1}\,\simeq\,10^{-10}\,{\rm dyn\,cm^{-2}}$. Therefore, if we are able to
estimate the shocked volume with time, we can get an approximate value for the
pressure jump. We can do this in a very raw way, taking the shocked volume to
be a cylinder, with radius $r_{s}$ and length $l_{s}$ (a usual approximation
used in analytical models, e.g., Begelman & Cioffi 1989):
$V\,\sim\,\pi\,r_{s}^{2}\,l_{s}$. Considering the mean head advance velocity
at small evolution times (within the host galaxy), we can write
$r_{s}\,\sim\,v_{s}\,t$ and $l_{s}\,\sim\,v_{h}\,t$, where $v_{s}$ is the
lateral shock velocity, which is the one we are interested in, and $v_{h}$ is
the head advance velocity. As stated above, $v_{s}$ is smaller than $v_{h}$.
Self-similarity appears as a result of keeping the ratio $v_{h}/v_{s}$
constant. Kaiser & Alexander (1997) showed that this is a consequence of
adiabatic expansion of jet-inflated cocoons. This is supported by observed
lobe/hot-spot$-$distance evolution throughout the inner kpc (Perucho & Martí
2002; Polatidis & Conway 2003; Kawakatu et al. 2008) in the case of Compact
Symmetric Objects. With this information, we can write the expression for the
volume as a function of time:

$V(t)\,\sim\,\pi\,v_{0,s}^{2}\,v_{0,h}\,t^{3-3\alpha}.$ (5)

For the first time step we can safely assume strong shock conditions
$P_{2}\,=\,\rho(r)\,v_{s}^{2}$ (where $\rho(r)$ is the ambient density). With
this, we obtain

$v_{s,0}\,=\,\left(\frac{L_{j}}{3\,\pi\rho_{0}t_{0}^{2}}\,\right)^{1/5}.$ (6)

Figure 1: Left panel: radial shock velocity ($v_{s}$, violet) and axial
position (black) for jet powers $10^{46}$ (solid lines), $10^{44}$ (dashed),
and $10^{42}\,{\rm erg/s}$ (dotted). The thin, dark red lines indicate the ISM
sound speed at the axial position. Right panel: post-shock temperatures for
same jet powers as the left panel in the case of a medium density following
Eq. 7 (black lines), i.e., a hot, diffuse ISM, and a medium with constant
number density $10^{3}\,{\rm cm^{-3}}$, i.e., cold clouds. Horizontal lines
indicate hydrogen ionization temperatures in local thermodynamical
equilibrium, $3500\,$K (solid), and collisional equilibrium, $1.58\times
10^{5}\,$K (dashed). In these clouds, the sound speed is $3\,{\rm km/s}$. The
temperatures have been derived for ionised hydrogen in Eq. 10.

Finally, we assume a typical ambient medium density profile

$\rho(r)\,=\,\rho_{0}\,\frac{1}{1\,+\,\left(\frac{r}{r_{c}}\right)^{2}},$ (7)

and $\rho_{0}\,=\,1.67\times 10^{-24}\,\rm{g/cm^{3}}$.

We can provide a general expression to estimate post-shock gas temperature by
starting from Eq. 4. Taking an initial volume of $\sim 100\,{\rm pc}^{3}$ as a
reference at $t\sim 10^{3}\,$yr, we get:

$P_{2}(t_{0})\,\simeq\,3\times 10^{-7}\,\left(\frac{L_{j}}{10^{45}\,{\rm
erg/s}}\right)\,\left(\frac{t}{10^{3}\,{\rm
yr}}\right)\,\left(\frac{V}{(100\,{\rm pc})^{3}}\right)^{-1}\,{\rm
dyn\,cm^{-2}},$ (8)

i.e., a factor $10^{3}$ over the ISM pressure (validating the strong shock
hypothesis at $t_{0}$). The evolution of pressure with time can then be
written from the previous expression and the time-dependency obtained for the
volume (see Eq. 5):

$P_{2}(t)\,\simeq\,3\times 10^{-7}\,\left(\frac{L_{j}}{10^{45}\,{\rm
erg/s}}\right)\,\left(\frac{t}{10^{3}\,{\rm yr}}\right)^{-2+3\alpha}\,{\rm
dyn\,cm^{-2}},$ (9)

where the constants in Eq. 5 have been absorbed in the initial volume taken to
estimate $P_{2}(t_{0})$.

Post-shock temperature can be estimated from (for a strong, non-relativistic
shock):

$T_{2}(t)\,=\,\frac{P_{2}(t)}{\rho_{2}}\,\frac{\mu\,m_{H}}{k}\,\simeq\frac{P_{2}(t)}{4\,\rho_{1}}\,\frac{\mu\,m_{H}}{k}.$
(10)

Taking $\rho_{1}\,=\,\rho_{0}\simeq 1\,m_{p}\,{\rm cm^{-3}}$ as a reference
value for the unperturbed ISM, we obtain:

$T_{2}(t)\,\simeq\,2\times 10^{9}\left(\frac{L_{j}}{10^{45}\,{\rm
erg/s}}\right)\left(\frac{\rho_{1}}{1\,m_{p}\,{\rm
cm^{-3}}}\right)^{-1}\,\left(\frac{t}{10^{3}\,{\rm
yr}}\right)^{-2+3\alpha}\,\rm{K}.$ (11)

The equation shows that the post-shock temperature is linearly dependent on
the jet power and inversely proportional to ISM gas density. The dependence on
time is controlled by exponent $\alpha$ for velocity evolution. Plausible
values range from $\alpha=0$, which is a consequence of a very steep density
profile ($\rho(r)\propto r^{-2}$) or a buoyantly propagating bubble with
constant velocity, and $\alpha\leq 1$ for low-power and or mass-loaded jets
(Perucho et al. 2011, 2014a). Self-similarity implies $\alpha\sim 0.4$
(Scheuer 1974; Kaiser & Alexander 1997), so post-shock temperature would fall
with time as $t^{-0.8}$. The sensitivity of $T_{2}$ to source age would be
damped for $\alpha\simeq 0.5-0.6$. On the contrary, it would be stronger if
velocity is constant ($\alpha=0$ would make $T_{2}\propto t^{-2}$, owing to a
fast pressure decrease) or if it fell rapidly ($\alpha\simeq 1$ would make
post-shock pressure become proportional to $t$, due to the increase in
pressure if the advance is slowed down). For moderate values of $\alpha$ as
that given by self-similarity, post-shock temperature might fall approximately
as $t^{-1}$ when propagating through an ambient medium with constant density,
i.e., the galactic core.

Therefore, beyond this initial set of values, we could impose self-similar
expansion and strong-shock conditions all along the studied evolution. In this
case $\alpha=2/5$ and the evolution is determined by the previous equations.
However, as we will see, the latter assumption could be not fulfilled for low
power jets. In order to avoid this problem, we can proceed in a different way:
the shock advance velocity in the radial direction at time $t_{i}$,
$v_{s}(t_{i})$, can be estimated from Eq. 3, i.e., skipping the strong-shock
assumption. Once $v_{s}(t_{i})$ is derived (for time $t_{i}$), we can obtain
$V(t_{i+1})$ assuming constant expansion velocity for a sufficiently small
time-step.

However, Eq. 5 requires not only $v_{s}$ but also $v_{h}$ to be known. We have
seen that self-similarity can be a good assumption through the inner region
and, in addition, the lobe length-to-radius ratio is $\leq 5$, so we can
assume a constant relation between both expansion velocities, e.g.,
$v_{h}\simeq 3v_{s}$. Although this assumption is probably weaker later on,
when the ratio can decrease due to jet head deceleration, $v_{h}/v_{s}$ must
be always larger than 1, so we can impose the same velocity ratio along the
whole studied distance range. Detailed analytical modelling by Turner et al.
(2018) showed that, in realistic environments, the ratio of jet head to
transverse expansion velocity only increases moderately (factor 1.5 on scales
of 100s of kpc; see Turner et al.’s Figures 8 and 9). As shown in that paper,
these analytical results are consistent with numerical simulations, and also
observations of lobed 3C radio sources. Furthermore, this estimate cannot
introduce significantly larger errors than, for instance, considering a
cylindrical volume for the shocked region.

Finally, once we obtain the value of $V(t_{i+1})$, we find $P_{2}(t_{i+1})$
from Eq. 4, and proceed to compute $v_{s}(t_{i+1})$. In this way, we can avoid
the dependency on $\alpha$ for the expansion velocity.

Figure 1 shows the evolution of shock position and velocity along the jet
propagation direction (left panel) and post-shock temperatures (right panel)
for different jet powers. The right panel also displays post-shock
temperatures for media with constant number density $10^{3}\,{\rm cm^{-3}}$
(and pressure $P_{1}$) in order to show the values attained inside clouds
shocked at different locations within the host galaxy (see Sect. 2.1.2). We
have set $t_{0}=10^{3}$ yr and have used constant ambient pressure
$P_{1}=10^{-10}\,{\rm dyn/cm^{2}}$.

Figure 2: Lobe overpressure for jet powers $10^{46}$ (solid lines), $10^{44}$
(dashed), and $10^{42}\,{\rm erg/s}$. The black lines for the hot, diffuse
ISM, and red lines correspond to cold, dense clouds with number density
$10^{3}\,{\rm cm^{-3}}$ are overlapped.

The right panel in Fig. 1 shows that post-shock temperatures reached by the
hot, dilute ISM agree with those obtained for powerful radio-galaxies from
numerical simulations. The ratios $P_{2}/P_{1}$ for pressure jumps across the
shocks (see Fig. 2) indicate that overpressure goes below a factor 10 only in
the case of the low-power radio galaxies after $500\,{\rm kyr}$, when the jet
head is at $\simeq 1\,{\rm kpc}$, and reach values $\simeq 7$ at $1{\rm Myr}$
($D\simeq 1.2\,{\rm kpc}$). Therefore, the strong shock approximation would
remain valid throughout most of the shown profiles for powerful jets, but low
power radio galaxies may transit to weaker (and radiative, see Sect. 2.2)
shocks within the inner 2 kpc of propagation.

####  2.1.2 The cold phase

The red lines in the right panel of Fig. 1 and Fig. 2 show the estimated
values of temperature and overpressure for shock-cloud interactions, assuming
$n\,=\,10^{3}\,{\rm cm^{-3}}$, for the three considered jet powers.

We have assumed $v_{1,c}\sim(\rho_{1}/\rho_{1,c})^{1/2}v_{1}$ (subscript $c$
stands for clouds, and no subscript for the values derived for the hot phase
in the previous section). Taking $P_{1,c}\,=\,P_{1}$ and $\Gamma=5/3$ for a
cold, monoatomic gas (hydrogen), we can derive $v_{2,c}$ by rearranging the
jump Eqs. 1-3, and then $P_{2,c}$.

Figure 3: HI fraction in shocked clouds in LTE for jets with powers $10^{44}$
(dashed), and $10^{42}\,{\rm erg/s}$ (dotted). The intermediate case of
$10^{43}\,{\rm erg/s}$ (dash-dot line) has been also included in this plot.

Assuming collisional ionization, which is expected to be a valid approximation
for diffuse, astrophysical plasmas (Dopita & Sutherland 2003), hydrogen
ionising temperatures are $\sim 10^{5}\,$K. Figure 1 shows that, for
$L_{j}\,=\,10^{46}\,{\rm erg/s}$ temperatures $\simeq 10^{5}\,$K are reached
in shocked clouds up to $\simeq 100\,{\rm kyr}$ ($D\simeq 2\,{\rm kpc}$). In
the case of $L_{j}\,=\,10^{44}\,{\rm erg/s}$, this happens before $50\,{\rm
kyr}$ ($D\simeq 300\,{\rm pc}$). Finally, low power jets may barely reach
ionising post-shock temperatures in dense clouds within the inner tens of
parsecs.

Although local thermodynamical equilibrium (LTE), in which the ionised states
follow Boltzmann’s law (Dopita & Sutherland 2003), is probably not fulfilled
in diffuse astrophysical plasmas, we can explore this possibility, too, for
the case of clouds. Using the post-shock values derived for clouds and Saha
equation, we can derive the degree of ionization in shocked clouds. In this
case, ionising temperatures are $\simeq 3000\,$K. The results are shown in
Fig. 3 for $L_{j}\,=\,10^{42}$, $10^{43}$ and $10^{44}\,{\rm erg/s}$. Low
power jets would be unable to ionise hydrogen in cold, dense clouds, beyond a
few hundred parsecs, whereas the first traces of HI in shocked clouds would
appear after $D\simeq 3\,{\rm kpc}$ in $10^{44}\,{\rm erg/s}$ jets. Post-shock
temperatures are incompatible with the presence of non-ionised HI in shocked
clouds in LTE for higher powers. Therefore, higher power jets would require
efficient cooling to explain the presence of atomic outflows within the inner
kiloparsecs of propagation, if LTE is assumed.

These values are only a basic approximation because shocks in clouds are
probably radiative (see Sect. 2.2), allowing for fast cooling. As an opposing,
heating mechanism, the development of instabilities and mixing with the hotter
environment could change this picture, as we later discuss.

###  2.2 The cooling

Figure 4: Cooling times for the same cases given in Fig. 1 (left panel) and
relative values with respect to the evolution times (right panel). Solid lines
represent the powerful, $10^{46}\,{\rm erg/s}$ case, dashed lines the
$10^{44}\,{\rm erg/s}$ case, and dotted lines the $10^{42}\,{\rm erg/s}$.
Calculations are shown for the shocked hot, dilute ISM described by Eq. 7.

As we have stated in the Introduction, outflows of atomic gas have been
observed both in GPS ($l_{s}\leq 1\,{\rm kpc}$) and CSS sources ($l_{s}\leq
10\,{\rm kpc}$). The values obtained for temperature in the previous section
(see Eq. 11) and in numerical simulations of powerful jets certainly prohibit
the presence of atomic hydrogen when the radio-galaxy extends over the inner
hundreds of parsecs. Efficient cooling must take place to explain the presence
of atomic outflows correlated with the presence of more powerful ($L_{j}\geq
10^{43}\,{\rm erg/s}$) relativistic jets.

Dense clouds are expected to be located within the inner kiloparsec of the
host galaxy. In the case of powerful jets that ionise most of the shocked ISM,
we then expect the observed atomic outflows to be formed by the recombination
of those ionised atoms, once the radio-source has evolved to larger scales.

According to cooling theory (see, e.g., Myasnikov et al. 1998; Dopita &
Sutherland 2003; Draine 2011) the cooling function, $\Lambda/n_{H}n_{e}$ can
be approximated to $1.1\times 10^{-22}\,(T/10^{6})^{-0.7}\,{\rm
erg\,cm^{3}\,s^{-1}}$ for gas temperatures between $10^{5}\,$K and $2\times
10^{7}\,$K, beyond which it is dominated by Bremsstrahlung ($2.3\times
10^{-24}\,(T/10^{6})^{1/2}{\rm erg\,cm^{3}\,s^{-1}}$), for number density
$n_{H}=1\,{\rm cm^{-3}}$, and, at lower temperatures (between $10^{4}\,$ and
$10^{5}\,$K), it becomes $7\times 10^{-27}\,T{\rm erg\,cm^{3}\,s^{-1}}$. Below
$10^{4}\,$K, we use the cooling function for recombination of hydrogen atoms
($2.39\times 10^{-27}\,T^{1/2}{\rm erg\,cm^{3}\,s^{-1}}$, see Vaidya et al.
2015). Using the values for post-shock densities and temperatures obtained in
Sect. 2.1, we can compute the cooling times of the shocked gas.

The cooling times are computed as follows:

$\tau\,=\,\frac{(\Gamma-1)\,P_{2}}{n_{2}^{2}\,2.3\times
10^{-24}(T_{2}/10^{6})^{0.5}}\,{\rm s},\quad T_{2}>2\times 10^{7}{\rm K},$
$\tau\,=\,\frac{(\Gamma-1)\,P_{2}}{n_{2}^{2}\,1.1\times
10^{-22}(T_{2}/10^{6})^{-0.7}}\,{\rm s},\quad 10^{5}\,{\rm K}<T_{2}<2\times
10^{7}{\rm K},$ $\tau\,=\,\frac{(\Gamma-1)\,P_{2}}{n_{2}^{2}\,7.0\times
10^{-27}T_{2}}\,{\rm s}\,\quad 10^{4}\,{\rm K}<T_{2}<10^{5}{\rm K},$
$\tau\,=\,\frac{(\Gamma-1)\,P_{2}}{n_{2}^{2}\,2.39\times
10^{-27}T_{2}^{0.5}}\,{\rm s}\quad T_{2}<10^{4}{\rm K}.$

The results are shown in Fig. 4, where we observe that low-power jets present
cooling times for the shocked hot, dilute ISM of the order or less than
evolution times after $10^{5}\,$yr, when they become radiative. The cooling
times derived for cold, dense clouds are shorter than 1 kyr and are thus not
shown. This means that cooling in the clouds is basically instantaneous in
relative terms. Therefore, although fast jet expansion and high post-shock
temperatures would make it difficult to find atomic hydrogen within the inner
kiloparsec in powerful jets, such short cooling times could make recombination
plausible in shocked clouds even before a linear size of 2 kpc is reached.

##  3 Discussion

###  3.1 Low power radio galaxies

We have seen that the bow shock can weaken and cease to be ’strong’ beyond a
few kpc for jet powers $L_{j}\,<\,10^{43}\,{\rm erg/s}$ (see Figs. 1 and 2).
Even though post-shock temperatures are high for the hot, dilute ISM phase,
cooling times become similar to source age when the tip of the bow shock is
still located within the inner kpc, and keep decreasing farther on. In the
case of cold, dense clouds, however, the post-shock temperatures are always
$\leq 10^{5}$ K. This means that 1) it is possible that part of the atomic gas
is not ionised by the shock, and 2) cooling times are basically instantaneous
compared to the evolution time-scales.

As the shock weakens, post-shock temperature and velocity become smaller for
the hot, dilute gas surrounding the clouds. Because the difference in
conditions inside and outside shocked clouds would be then alleviated, shocked
clouds could be more stable within this environment than in that generated by
stronger shocks. In addition, even if mixing occurred, the lower surrounding
temperatures would still allow for rapid cooling. This is, however, to be
studied by devoted numerical simulations.

Nevertheless, the simple analytical approach that we have used already
indicates that the presence of atomic gas with outward motions in shocked
regions within radio galaxies could be relatively easily explained.

###  3.2 Powerful radio galaxies

In Sect. 2.2 we have seen that, despite the high temperatures achieved by the
ISM after the passage of shocks, large density regions could allow fast enough
cooling and, therefore, the presence of atomic (and molecular) gas within the
inner kiloparsec(s) from the active nucleus. In addition, this gas has an
outward momentum that can be aligned with the jet propagation direction or be
perpendicular to it because the bow shocks propagate in all directions, in
accordance with observations. The observed lines reveal radial velocities of
hundreds up to few thousand km/s, and these velocity fields are indeed
reproduced in numerical simulations (Perucho et al. 2021; Meenakshi et al.
2022; Murthy et al. 2022).

The higher temperatures achieved by the shocked gas in powerful radio-sources
and thus longer cooling times could lead us to expect a bias in the detection
of atomic outflows towards larger distances to the nucleus in this case.
Furthermore, (shocked) cloud stability has to be taken into account:
disruption and mixing with the hotter environment would significantly delay
cooling times/distance. An _a fortiori_ argument favouring cloud stability
conditions is given by observations. Clouds must retain a certain degree of
coherence after shock passage if atomic and molecular lines are observed
within shock affected regions in radio-galaxies. Otherwise, this would be
difficult to explain for relatively powerful radio-galaxies.

In a recent paper, Snow & Hillier (2021a) explore the physics of shock
propagation through partially ionised plasmas and show that post-shock
equilibrium in such systems can lead to rapid cooling and lower post-shock
than pre-shock temperatures. If this scenario is confirmed and generalised in
the context of AGN shocks and clouds, it could provide an easy explanation to
the presence of atomic lines associated with young radio-galaxies.

Unstable processes evolve on timescales of $\sim R/c_{s}$, where $R$ is the
cloud radius and $c_{s}$ its sound speed. The relatively low temperatures kept
by cloud material (both unshocked and shocked) ensures relatively low sound
speeds. For a cloud size of 10 pc and a sound speed (after shock passage) of
$10^{2}-10^{3}\,{\rm km/s}$ dynamical times become $\sim 10^{4}-10^{5}$ yr,
i.e. of the order of typical evolution times shown in Fig. 1 and certainly
longer than cooling times. There could thus be time for cooling and
recombination prior to eventual disruptive processes. So, even if mixing would
induce further ionization, there could be a long enough time interval between
shock passage and disruption when atomic lines could be observed. Another
physical mechanism to be taken into account is enhanced clumping induced by
self-gravity could also contribute to further cooling and stability (Mandal et
al., in preparation).

As stated above, the development of instabilities would favour mixing with the
hotter, shocked surrounding gas, acting against cooling and recombination in
the region. Numerical simulations of radio source evolution (see, e.g.
Mukherjee et al. 2016; Perucho et al. 2021) do not show any hints of post-
shock gas colder than $10^{8}\,$K, although the densest clouds in set-ups do
have densities $\geq 100\,{\rm cm^{-3}}$. The reason is numerical: the clouds
in such simulations do not involve many cells ($\leq$ 20 cells, in general),
favouring rapid mixing by numerical diffusion. We therefore need devoted
numerical simulations of shock-cloud interactions (evolved versions of the one
presented here) to study this process in detail. It would be of particular
interest to study the dependence of cloud stability on the properties of the
surrounding shocked gas, which shows very different thermodynamic properties
for different jet powers (see, e.g., Figs. 2 and 4). These runs are, however,
costly owing to the relatively small shock propagation velocity in clouds and
out of the scope of this theoretical approach.

###  3.3 Implications for large-scale feedback

The presence of atomic outflows with large mass fluxes has been interpreted
sometimes in the literature as jets being cold, strongly mass loaded and slow
at large scales. This conclusion has been applied to study jet feedback at
large scales (beyond the host galaxy Vernaleo & Reynolds 2006; Cattaneo &
Teyssier 2007; Sternberg & Soker 2008; Gaspari et al. 2011), even for powerful
outflows ($L_{j}\sim 10^{46}\,{\rm erg/s}$). In contrast, numerical
simulations of powerful, relativistic jets show that it is very difficult to
affect the global dynamics of powerful jets by stellar wind mass-load (Perucho
et al. 2014a; Anglés-Castillo et al. 2021). In addition, observational studies
have shown that such powerful jets must be mildly relativistic at large
scales, on the basis of jet-to-counter-jet brightness ratios (e.g., Arshakian
& Longair 2004).

Furthermore, the scenario we propose in this paper, which is based on results
obtained from numerical simulations and is compatible with the physical scales
of the problem, provides a reasonable frame that indicates that cold, massive
outflows observed in active galaxies are produced as a result of the
development of relativistic jets, but do not represent the nature of jets as
it was assumed by a number of works (see, e.g., the references above) for
studying jet feedback at hundreds of kiloparsecs. The same conclusion has also
been derived from observational results (see, e.g. Dasyra et al. 2014, who
show that massive outflows can be explained from a small fraction of the total
jet power or Laing & Bridle (2002) for the case of the FRI radio-galaxy 3C
31). Altogether, the possibility of massive, cold outflows with powers $\sim
10^{46}\,{\rm erg/s}$ not only not playing a role in large-scale feedback, but
also not existing at all in the Universe, would have to be taken into serious
consideration in cosmological or large-scale feedback simulations. In
contrast, low power sources can be efficiently mass-loaded (Bowman et al.
1996; Perucho et al. 2014a; Anglés-Castillo et al. 2021), losing collimation
and decelerating. However, the resulting outflows are not expected to reach
hundreds of kiloparsecs as classical radio-galaxies do.

##  4 Summary

In this work, I present a general physical scenario to explain the presence of
atomic gas outflows in GPS and CSS radio-sources. I have based the model in
basic knowledge obtained from both numerical simulations and analytical models
that have been systematically applied to analyse them (radio-galaxy evolution,
e.g., Perucho & Martí 2007; Perucho et al. 2014b, 2019).

The model predicts that atomic outflows could easily be detected in the GPS
regime ($l_{s}\leq 1\,{\rm kpc}$) for low-power jets ($L_{j}\leq 10^{43}\,{\rm
erg/s}$), as post-shock temperatures in clouds could easily fall below
$10^{4}\,{\rm K}$. In addition, the cooling times of the shocked hot phase can
be shorter than the evolution time-scales.

On the contrary, the post-shock temperatures reached in the case of powerful
jets could ionise a large percentage of the existing atomic gas at these
scales. Thus, the presence of atomic gas in powerful radio-galaxies at those
scales requires rapid cooling in shocked clouds. This relies on cloud
stability for sufficiently long periods of time, which seems to be granted by
the small sound speeds expected for post-shock temperatures of
$10^{4}-10^{5}\,$K (see Fig. 1). In this case, the cooling times allowing
recombination to take place are compatible with the dynamical times required
by jets to reach the CSS scales ($\geq 2\,{\rm kpc}$). In conclusion, a
positive correlation between linear size and jet power should be expected in
radio-galaxies showing atomic outflows.

Devoted numerical simulations of shock-cloud interaction must be run in order
to study the stability and evolution of shocked cloud gas, also including the
role of magnetic fields on the stability of the system (see, e.g., Snow &
Hillier 2021a, for a recent work on the physics of shocked plasmas that could
be relevant in this context). It has actually been shown that magnetic fields
can be relevant to damp instability growth even if they are energetically
irrelevant (López-Miralles et al. 2022). If shocked clouds retain at least
partly coherence, these regions can cool faster because of their larger
densities and lower post-shock temperatures; in contrast, if instabilities and
turbulent mixing with the shocked hotter ISM component, the resulting
temperature of the mixture could delay cooling and put constraints to the
validity of this simple model. Nevertheless, observational evidence seems to
require shocked clouds to be fairly stable.

## References

  * Anglés-Castillo et al. (2021) Anglés-Castillo, A., Perucho, M., Martí, J. M., & Laing, R. A. 2021, MNRAS, 500, 1512
  * Arshakian & Longair (2004) Arshakian, T. G. & Longair, M. S. 2004, MNRAS, 351, 727
  * Begelman (1996) Begelman, M. C. 1996, in Cygnus A – Studay of a Radio Galaxy, ed. C. L. Carilli & D. E. Harris, 209
  * Begelman & Cioffi (1989) Begelman, M. C. & Cioffi, D. F. 1989, ApJ, 345, L21
  * Bicknell et al. (1997) Bicknell, G. V., Dopita, M. A., & O’Dea, C. P. O. 1997, ApJ, 485, 112
  * Bicknell et al. (2018) Bicknell, G. V., Mukherjee, D., Wagner, A. Y., Sutherland, R. S., & Nesvadba, N. P. H. 2018, MNRAS, 475, 3493
  * Bicknell et al. (2003) Bicknell, G. V., Saxton, C. J., & Sutherland, R. S. 2003, PASA, 20, 102
  * Blandford & Payne (1982) Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433
  * Bowman et al. (1996) Bowman, M., Leahy, J. P., & Komissarov, S. S. 1996, MNRAS, 279, 899
  * Cattaneo & Teyssier (2007) Cattaneo, A. & Teyssier, R. 2007, MNRAS, 376, 1547
  * Couto et al. (2020) Couto, G. S., Storchi-Bergmann, T., Siemiginowska, A., Riffel, R. A., & Morganti, R. 2020, MNRAS, 497, 5103
  * Dasyra et al. (2014) Dasyra, K. M., Combes, F., Novak, G. S., et al. 2014, A&A, 565, A46
  * Dopita & Sutherland (2003) Dopita, M. A. & Sutherland, R. S. 2003, Astrophysics of the diffuse universe
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium, Princeton University Press
  * Fanti et al. (1990) Fanti, R., Fanti, C., Schilizzi, R. T., et al. 1990, A&A, 231, 333
  * Gaspari et al. (2011) Gaspari, M., Brighenti, F., D’Ercole, A., & Melioli, C. 2011, MNRAS, 415, 1549
  * Hardcastle (2018) Hardcastle, M. J. 2018, MNRAS, 475, 2768
  * Holt et al. (2011) Holt, J., Tadhunter, C. N., Morganti, R., & Emonts, B. H. C. 2011, MNRAS, 410, 1527
  * Kaiser & Alexander (1997) Kaiser, C. R. & Alexander, P. 1997, MNRAS, 286, 215
  * Kawakatu et al. (2008) Kawakatu, N., Nagai, H., & Kino, M. 2008, ApJ, 687, 141
  * Komissarov (2012) Komissarov, S. 2012, in Relativistic Jets from Active Galactic Nuclei, 81–114
  * Komissarov et al. (2007) Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51
  * Krause & Alexander (2007) Krause, M. & Alexander, P. 2007, MNRAS, 376, 465
  * Krause (2008) Krause, M. G. H. 2008, Mem. Soc. Astron. Italiana, 79, 1162
  * Kukreti et al. (2023) Kukreti, P., Morganti, R., Tadhunter, C., & Santoro, F. 2023, A&A, 674, A198
  * Labiano, A. et al. (2005) Labiano, A. et al. 2005, A&A, 436, 493
  * Laing & Bridle (2002) Laing, R. A. & Bridle, A. H. 2002, MNRAS, 336, 1161
  * López-Miralles et al. (2022) López-Miralles, J., Perucho, M., Martí, J. M., Migliari, S., & Bosch-Ramon, V. 2022, A&A, 661, A117
  * Meenakshi et al. (2022) Meenakshi, M., Mukherjee, D., Wagner, A. Y., et al. 2022, MNRAS, 511, 1622
  * Morganti & Oosterloo (2018) Morganti, R. & Oosterloo, T. 2018, A&ARv, 26, 4
  * Morganti et al. (2021) Morganti, R., Oosterloo, T., Murthy, S., & Tadhunter, C. 2021, Astronomische Nachrichten, 342, 1135
  * Morganti et al. (2016) Morganti, R., Veilleux, S., Oosterloo, T., Teng, S. H., & Rupke, D. 2016, A&A, 593, A30
  * Mukherjee et al. (2016) Mukherjee, D., Bicknell, G. V., Sutherland, R., & Wagner, A. 2016, MNRAS, 461, 967
  * Mukherjee et al. (2018) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544
  * Murthy et al. (2022) Murthy, S., Morganti, R., Wagner, A. Y., et al. 2022, Nature Astronomy, 6, 488
  * Myasnikov et al. (1998) Myasnikov, A. V., Zhekov, S. A., & Belov, N. A. 1998, MNRAS, 298, 1021
  * O’Dea & Baum (1997) O’Dea, C. P. & Baum, S. A. 1997, AJ, 113, 148
  * O’Dea et al. (2002) O’Dea, C. P., de Vries, W. H., Koekemoer, A. M., et al. 2002, AJ, 123, 2333
  * Owsianik & Conway (1998) Owsianik, I. & Conway, J. E. 1998, A&A, 337, 69
  * Perucho et al. (2021) Perucho, M., López-Miralles, J., Reynaldi, V., & Labiano, Á. 2021, Astronomische Nachrichten, 342, 1171
  * Perucho & Martí (2002) Perucho, M. & Martí, J. M. 2002, ApJ, 568, 639
  * Perucho & Martí (2007) Perucho, M. & Martí, J. M. 2007, MNRAS, 382, 526
  * Perucho et al. (2014a) Perucho, M., Martí, J. M., Laing, R. A., & Hardee, P. E. 2014a, MNRAS, 441, 1488
  * Perucho et al. (2019) Perucho, M., Martí, J.-M., & Quilis, V. 2019, MNRAS, 482, 3718
  * Perucho et al. (2017) Perucho, M., Martí, J.-M., Quilis, V., & Borja-Lloret, M. 2017, MNRAS, 471, L120
  * Perucho et al. (2014b) Perucho, M., Martí, J.-M., Quilis, V., & Ricciardelli, E. 2014b, MNRAS, 445, 1462
  * Perucho et al. (2011) Perucho, M., Quilis, V., & Martí, J.-M. 2011, ApJ, 743, 42
  * Polatidis & Conway (2003) Polatidis, A. G. & Conway, J. E. 2003, PASA, 20, 69
  * Raouf et al. (2019) Raouf, M., Silk, J., Shabala, S. S., et al. 2019, MNRAS, 486, 1509
  * Reynaldi & Feinstein (2016) Reynaldi, V. & Feinstein, C. 2016, MNRAS, 455, 2242
  * Santoro et al. (2020) Santoro, F., Tadhunter, C., Baron, D., Morganti, R., & Holt, J. 2020, A&A, 644, A54
  * Scheuer (1974) Scheuer, P. A. G. 1974, MNRAS, 166, 513
  * Schulz et al. (2021) Schulz, R., Morganti, R., Nyland, K., et al. 2021, A&A, 647, A63
  * Shabala & Alexander (2009) Shabala, S. & Alexander, P. 2009, ApJ, 699, 525
  * Shih et al. (2013) Shih, H.-Y., Stockton, A., & Kewley, L. 2013, ApJ, 772, 138
  * Snow & Hillier (2021a) Snow, B. & Hillier, A. 2021a, A&A, 645, A81
  * Snow & Hillier (2021b) Snow, B. & Hillier, A. 2021b, MNRAS, 506, 1334
  * Stanghellini et al. (1998) Stanghellini, C., O’Dea, C. P., Dallacasa, D., et al. 1998, A&AS, 131, 303
  * Sternberg & Soker (2008) Sternberg, A. & Soker, N. 2008, MNRAS, 384, 1327
  * Tadhunter et al. (2021) Tadhunter, C., Morganti, R., Santoro, F., & Bernhard, E. 2021, Astronomische Nachrichten, 342, 1200
  * Turner & Shabala (2023) Turner, R. J. & Shabala, S. S. 2023, Galaxies, 11, 87
  * Turner et al. (2018) Turner, R. J., Shabala, S. S., & Krause, M. G. H. 2018, MNRAS, 474, 3361
  * Turner et al. (2023) Turner, R. J., Yates-Jones, P. M., Shabala, S. S., Quici, B., & Stewart, G. S. C. 2023, MNRAS, 518, 945
  * Vaidya et al. (2015) Vaidya, B., Mignone, A., Bodo, G., & Massaglia, S. 2015, A&A, 580, A110
  * Vernaleo & Reynolds (2006) Vernaleo, J. C. & Reynolds, C. S. 2006, ApJ, 645, 83
  * Vlahakis & Königl (2004) Vlahakis, N. & Königl, A. 2004, ApJ, 605, 656
  * Wagner & Bicknell (2011) Wagner, A. Y. & Bicknell, G. V. 2011, ApJ, 728, 29
  * Wagner et al. (2012) Wagner, A. Y., Bicknell, G. V., & Umemura, M. 2012, ApJ, 757, 136
  * Zovaro et al. (2019) Zovaro, H. R. M., Nesvadba, N. P. H., Sharp, R., et al. 2019, MNRAS, 489, 4944

###### Acknowledgements.

The author acknowledges support by the Spanish Ministry of Science through
Grants PID2019- 105510GB-C31/AEI/10.13039/501100011033 and
PID2022-136828NB-C43, from the Generalitat Valenciana through grant
CIPROM/2022/49, and from the Astrophysics and High Energy Physics programme
supported by MCIN and Generalitat Valenciana with funding from European Union
NextGenerationEU (PRTR-C17.I1) through grant ASFAE/2022/005. The author also
acknowledges the referee of this paper, Stas Shabala for his constructive
criticism, which has improved the resulting text, and Pau Beltran-Palau for
critical reading of the paper.


D. Mickaël / Ocean Modelling

Mickaël Delcey et al.

mode = title]Enhancing Gravity Currents Analysis through Physics-Informed Neural Networks: Insights from Experimental Observations

1]Mickaël Delcey[orcid = 0000-0003-2225-8531]

[1]organization=LEMTA, Université de Lorraine, CNRS,
addressline=2, Avenue de la Forêt de Haye, B.P. 160,

1]Yoann Cheny
[1]

1]Jean Schneider[

1]Simon Becker[]

1]Yvan Dossmann

[cor1]<EMAIL_ADDRESS>
[2]Institut Universitaire de France (IUF)

Sébastien Kiesgen De Richter

Gravity currents in oceanic flows require simultaneous measurements of pressure and velocity to assess energy flux, which is crucial for predicting fluid circulation, mixing, and overall energy budget. In this paper, we apply Physics Informed Neural Networks (PINNs) to infer velocity and pressure field from Light Attenuation Technique (LAT) measurements for gravity current induced by lock-exchange. In a PINN model, physical laws are embedded in the loss function of a neural network, such that the model fits the training data but is also constrained to reduce the residuals of the governing equations. PINNs are able to solve ill-posed inverse problems training on sparse and noisy data, and therefore can be applied to real engineering applications.
The noise robustness of PINNs and the model parameters are investigated in a 2 dimensions toy case on a lock-exchange configuration , employing synthetic data. Then  we train a PINN with experimental LAT measurements and quantitatively compare the velocity fields inferred to PIV measurements performed simultaneously on the same experiment. Finally, we study the energy flux field $J=p \boldsymbol{u}$ derived from the model. The results state that accurate and useful quantities can be derived from a PINN model trained on real experimental data which is encouraging for a better description of gravity currents and improve models of ocean circulation.

Physics Informed Neural Networks Gravity Currents Lock-Exchange Deep Learning

§ INTRODUCTION

Gravity currents are common in oceanic flows like ocean overflows, or hydrothermal vents for example. It is crucial to determine energy transport in such density varying flows to predict fluid circulation, mixing and overall budget. For this, it is necessary to evaluate the energy flux fields $J=p\boldsymbol{u}$ which requires simultaneous measurements of the pressure $p$ and velocity $\boldsymbol{u}$ simultaneously in $3$ dimensions ($3D$) in order to obtain an upper bound on the energy available for dissipation and turbulent mixing.

In recent years, Particule Image Velocimetry (PIV) ([Perez-Diaz et al., 2018]), Planar Laser Induced Fluorescence (PLIF) ([Crimaldi, 2008]) and Light Attenuation Technique (LAT) [Dossmann et al., 2016] have been widely used and combined to get simultaneous measurements of velocity and density fields in 2D ([Balasubramanian et al., 2018]). 3D measurements have also been performed by using 3D-LIF ([Tian et al., 2003]), Tomo-BOS ([Raffel, 2015]) and PIV-3D ([Elsinga et al., 2006], [Partridge et al., 2019]) techniques.
Although well established, these techniques are very costly and difficult to implement in practice in a simultaneous manner in the field. Furthermore, it is necessary to use complex inverse methods to extract the pressure field from the density measurements. Thus, achieving a complete reconstruction of the hydrodynamic field remains a significant challenge.
In their work, [Nash et al., 2005] proposed a method to extract the pressure field in the case of a linear perturbation of the latter by the density field. It succeeds to predict the time average energy flux but assumes that there is no contribution of the dynamic pressure. [Clark et al., 2010] managed to overcome this limitation by using the Boussinesq polarization relations for calculating the pressure but requires nearly periodic perturbations of the density field. These methods provide leading order approximations for the averaged energy flux from measurements but cannot capture transient features contrary to Green function methods recently proposed by [Allshouse et al., 2016] but which remains valid only in the case of linear or very weakly non-linear perturbations of the density field. To account for non-linear perturbations of the density field, and achieve a more realistic reconstruction, data assimilation techniques can be employed  ([Fourés et al., 2014], [Mons et al., 2022]) but this requires to run multiple CFD simulations which is prohibitively time consuming.

Recent advances in deep learning ([Goodfellow et al., 2016]) and the universal approximation nature of Neural Networks offer new possibilities for reconstruction problems in fluid mechanics ([Brunton et al., 2020]) even for non linear perturbations of the hydrodynamic fields. However, the chaotic nature of the flows investigated makes the reconstruction problem most of the time out of the reach of purely data driven algorithms. In that context, Raissi et al. ([Raissi et al., 2020]) recently proposed to use Physics-Informed Neural Network (PINNs) for solving forward and inverse partial differential equation problems which overcomes the low data availability in physical systems. It consists of adding prior knowledge about physical laws by constraining a neural network to respect a set of partial differential equations (PDEs). This method has quickly proven itself on various non-linear forward ([Mao, 2020]) , and inverse ([Cai et al., 2021]) problems for synthetic data (see [Cuomo et al., 2022] for a complete review). Since PINNs can handle noisy and gappy data, as well as ill-posed problems, it seems to be the ideal candidate to tackle inverse problems with experimental data, however few applications on this subject are available. In [Cai et al., 2021] Kaniadakis demonstrated the ability of PINNs to solve inverse problems using experimental Background-riented schlieren (BOS) data, however, velocity fields inferred have only been validated by a qualitative comparison on PIV  measurements in a companion experiment. Recently, PINNs have also been applied BOS data to reconstruct high speed flows fields ([Molnar et al., 2023]), but no validation of the velocity field is available.

This paper has two main objectives, the first is to provide for the first time a quantitative validation of the velocity fields obtained by PINNs. To do this, we combine LAT and PIV measurements on the same lock-exchange experiments, then we build a Physics Informed Neural Network model to infer velocity and pressure from density measurements only, finally, we compare the velocity field inferred with PIV measurements kept for validation.
The second objective is to demonstrate that this method can be used to directly calculate the energy flux $J=pu$ and to draw up energy budgets for strongly non-linear perturbations of the density field.

The principles of Physics informed Neural Networks and the governing equations of Gravity Currents flow will be discussed in section <ref> and <ref> respectively.
In section <ref> we study the noise robustness of PINNs on synthetical data for the lock exchange configuration.
We apply PINNs model on LAT experimental measurements from lock-exchange experiment in section <ref>, velocity fields inferred are quantitatively compared to PIV measurements performed simultaneously on the same experience.
Finally, we discuss the energy flux and vorticity derived from the inferred pressure and velocity fields.

§ PHYSICS INFORMED NEURAL NETWORKS.

Artificial Neural Networks (ANNs) are known as universal approximators ([Cybenko, 1989]) so they can in theory learn any measurable function from observed data.
A neural network consists of an input layer, one or more hidden layers, and an output layer. The layers are composed of nodes, also known as neurons, and each neuron is associated with a set of learnable parameters called weights and biases which are used to transform the input signals and compute the output of the model.
A network of $N$  layers performs a series of computations, known as forward propagation, that transform the input vector $x$ through each layer of the network until a final prediction $a^{[N]}$ is obtained. :

\begin{align}
& a^{[1]} = x \in \mathbb{R}^{n_{1}}, \\
& a^{[k]} = \sigma \left( W^{[k]}a^{[k-1]}+b^{[k]} \right) ,~ 2 \leq k \leq N ,\\
& a^{[N]} =  W^{[N]}a^{[N-1]}+b^{[N]},~ ~k = N ,
\end{align}

Here $W^{[k]}\in \mathcal{M}_{n_{k}n_{k-1}}(\mathbb{R})$, $b^{[k]} \in \mathbb{R}^{n_{k}}$ are the weights and biases of the $k^{th}$-layers and $\sigma$ is a non linear transformation called activation function. $n_k$ is the number of neurons of the layer $k$.

The learnable parameters of the networks : $\boldsymbol{\theta} = (W_{k},b_{k})_{k = 1\dots N}$ are adjusted through the training process until they reach an optimal value $\boldsymbol{\theta^{*}} = (W_{k}^*,b_{k}^*)_{k = 1\dots N}$ that minimizes a loss function $\mathcal{L}$  which penalizes the mismatch between the network predictions and observed data :

\begin{align}
\label{argmin}
~~~~~~ ~~~~~~  ~~~~~~ \huge \theta^{*}= \underset{\theta}{\argmin}  ~\mathcal{L}(\theta)
\end{align}

most of the time, for regression tasks, the loss function $\mathcal{L}$ is the Mean Squared Error between the observations $\{a_{i}\}_{i = 1 \dots N_{\textrm{data}}}$ and the corresponding predictions of the network $\{a^{[N]}(x_{i})\}_{i = 1 \dots N_{\textrm{data}}}$ where $N_\textrm{data}$ is the size of the observed dataset.

Weights are updated at each training step, using the gradient descent :

\begin{align}
\label{gradient}
~~~~~~ ~~~~~~  ~~~~~~  \theta_{n+1} = \theta_{n} - \eta \nabla_{\theta} \mathcal{L}(\theta),
\end{align}

where $\eta$ is called the learning rate.
The key point in the success of neural networks is backpropagation. Backpropagation ([Higham et al., 2019]) is an algorithm for computing the gradients of the loss function with respect to the learnable parameters in a neural network.

Nowadays deep learning libraries such as Tensorflow/Keras ([Géron, 2019]) and Pytorch ([Paszke et al., 2019]) include automatic differentiation (AD) tools ([Griewank et al., 1989]) that allow an efficient computation of these gradients.
In the context of physics-informed neural networks (PINNs) ([Raissi et al., 2020]), a neural network $\boldsymbol{v}_{\mathcal{N}}$ is used to approximate the solution $\boldsymbol{v}$ of partial differential equations like :

\begin{align}
&\partial_{t}\boldsymbol{v} + \mathcal{D}_{\textrm{\textbf{x}}}[\boldsymbol{v}],~~\textrm{\textbf{x}} \in \Omega,~~t\in [0,T], \label{eq:PINN:1}\\
&\boldsymbol{v}(\textrm{\textbf{x}},0) = h(x), x \in \Omega , \\
&\boldsymbol{v}(\textrm{\textbf{x}},t) = g(x,t),~~ \textrm{\textbf{x}}\in \partial \Omega, t\in [0,T],
\end{align}

where $\mathcal{D}_{\textrm{\textbf{x}}}$ is a linear or non linear differential operator , $ \Omega,~ \partial \Omega, [0,T]$ are respectively the spatial domain and his frontier, and the temporal domain. $h|_{\Omega}, g|_{\partial\Omega,[0,T]}$ are the initial and boundary conditions.

We define $f$ that represents the left-hand-side of <ref>:

\begin{align}
~~~~~~~~~~~~~~~f:=\partial_{t}\boldsymbol{v}+ \mathcal{D}_{\textrm{\textbf{x}}}[\boldsymbol{v}]
\end{align}

$v_{\mathcal{N}}$ takes as input $(\textrm{\textbf{x}},t)$ and AD enables the computation of the derivatives of $\boldsymbol{v}_{\mathcal{N}}$ with respect to $\textrm{\textbf{x}}$ and $t$, thus we can get an explicit expression of $f(\textrm{\textbf{x}},t),~\forall (\textrm{\textbf{x}},t)$.
This stage is the key instance of PINNs because automatic differentiation does not introduce truncation errors due to numerical approximations as with conventional derivation methods.
The objective function to minimize is expressed as :

\begin{align}
& \mathcal{L} = \lambda_{1}\mathcal{L}_{f}+\lambda_{2}\mathcal{L}_{b}+\lambda_{3}\mathcal{L}_{0}+\lambda_{4}\mathcal{L}_{\textrm{data}}, \\
& \mathcal{L}_{f} = \dfrac{1}{N_f}\sum_{i=1}^{N_{f}}||f(\textrm{\textbf{x}}^{f}_{i},t^{f}_{i})||^{2}, \\
& \mathcal{L}_{b} = \dfrac{1}{N_b}\sum_{i=1}^{N_{b}} ||\boldsymbol{v}_{\mathcal{N}}(\textrm{\textbf{x}}^{b}_{i},t^{b}_{i})-g(\textrm{\textbf{x}}^{b}_{i},t^{b}_{i})||^{2}, \\
& \mathcal{L}_{0} = \dfrac{1}{N_0}\sum_{i=1}^{N_{0}} ||\boldsymbol{v}_{\mathcal{N}}(\textrm{\textbf{x}}^{0}_{i},0)-h(\textrm{\textbf{x}}^{0}_{i})||^{2}, \\
& \mathcal{L}_{\textrm{data}} = \dfrac{1}{N_\textrm{data}}\sum_{i=1}^{N_{\textrm{data}}} ||\boldsymbol{v}_{\mathcal{N}}(\textrm{\textbf{x}}^{\textrm{data}}_{i},t^{\textrm{data}}_{i})-\boldsymbol{v}(\textrm{\textbf{x}}^{\textrm{data}}_{i},t^{\textrm{data}}_{i})||^{2}.
\end{align}

Here $\mathcal{L}_{\textrm{data}}$ is due to additional observations inside the domain, and $\mathcal{L}_{f},~\mathcal{L}_{b},~\mathcal{L}_{0}$ penalize the residual equations, the boundary, and the initial conditions.

The parameters $\lambda_{1\dots 4} \in \mathbb{R}_{+}$ are non trainable weights set manually to manage the ratio between all loss terms.

Collocation points for data, residuals and boundary conditions : $\{\textrm{\textbf{x}}^{*}_{i},t^{*}_{i} \}_{i=1 \dots N_{*}}$ are often selected on a discrete Cartesian grid.
In most of the cases the algorithm selected to optimize the parameter $\theta$ is ADAM ([Kingma and Ba, 2014]), a first-order gradient-based optimization algorithm for stochastic objective functions .

In the next section, we will apply PINNs to gravity currents generated by lock exchange, to infer the density, velocity and pressure field with synthetical data from 2D simulation. We can see a schematic representation of PINN for this problem Fig. <ref>.

Physic informed neural network structure for gravity currents governing equations: a fully connected neural network take as input $(\textrm{\textbf{x}},t) \in \mathbb{R}^{3}$ and predicts $\boldsymbol{v}_{\mathcal{N}}(\textrm{\textbf{x}},t) = (\rho_{\mathcal{N}},u_{\mathcal{N}},v_{\mathcal{N}},p_{\mathcal{N}})$. The residuals of the governing equations $e_{1,\dots,4}$  are computed by automatic differentiation and $e_0$ denotes the mismatch between the observational data $\boldsymbol{v}$ and the predictions $\boldsymbol{v}_{\mathcal{N}}$, which are combined in the loss function $\mathcal{L}_0, \mathcal{L}_b$ and $\mathcal{L}_\textrm{data}$. Finally, the weights and biases of the network are updated iteratively with gradient descent.

§ GOVERNING EQUATIONS OF LOCK-EXCHANGE GRAVITY CURRENTS

In the lock exchange flow ([Shin et al., 2004]), fluids of different densities initially at
rest are separated by a vertical membrane, in a tank. The fluids in the left and right hand compartments have the respective densities $\tilde{\rho_1}$ and $\tilde{\rho_2}$ with  $\tilde{\rho_1} > \tilde{\rho_2}$, the density gradient
is caused for instance by salinity difference. Once the membrane is withdrawn the density difference $\Delta\rho = \tilde{\rho_1} - \tilde{\rho_2}$ between the two compartments induces the generation and propagation of a gravity current at the tank bottom. The intensity of turbulent structures and induced mixing in the current is controlled by the Reynolds number $Re$. The flow is governed by Navier-Stokes, mass transport and incompressibility equations. Moreover, in this paper, we study Gravity Currents induced by small difference in density, so that Boussinesq Aproximation ([Gray and Giorgini, 1976]) is  effective.

Lock exchange tank : The lock-exchange configuration. At the initial moment, the wall is removed and the system evolves under the action of gravity only. The difference of density is achieved by a difference of salinity.

In order to render the equations dimensionless (see  [Hartel et al., 2000]), we use $\tL_y$ as the characteristic
length scale and the characteristic velocity is the buoyancy velocity defined as :

\begin{equation}\label{eq:buyvel}
\tilde{u}_b=\sqrt{g\frac{\tilde{\rho_1}-\tilde{\rho_2}}{\tilde{\rho_a}}\tilde{L}_y}, \quad \text{with}\quad\tilde{\rho}_a=\frac{\tilde{\rho_1}+\tilde{\rho_2}}{2}.
\end{equation}

We consider small density difference $\tilde{\rho}_1 - \tilde{\rho}_2\ll \tilde{\rho}_1$ for which the Boussinesq approximation can be adopted, thus the governing equations read in dimensionless form :

\begin{align}
& \dfrac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \cdot \nabla )~ \boldsymbol{u}  +\dfrac{\partial p}{\partial x} - \dfrac{1}{Re}\nabla ^{2} \boldsymbol{u}  \label{eq:momx} + \rho ~ \boldsymbol{e_{y}}=0,\\
& \dfrac{\partial \rho}{\partial t}+(\boldsymbol{ u} \cdot \nabla )~ \rho - \dfrac{1}{Re\,Sc}\nabla^{2}\rho \label{eq:mass}=0,\\
& \nabla \cdot \boldsymbol{ u} = 0\label{eq:cont},
\end{align}

where  $\boldsymbol{u} = (u,v)^\intercal$ denotes the dimensionless velocity vector. The non-dimensional pressure $p$ and density $\rho$ are given by :

\begin{equation}\label{eq:prhoadim}
p=\frac{\tilde{p}}{\tilde{\rho}_a\tilde{u}_b^2},\quad \rho=\frac{\tilde{\rho}-\tilde{\rho}_2}{\tilde{\rho_1}-\tilde{\rho_2}}.
\end{equation}

The Reynolds number $Re$ and the Schmidt number $Sc$ arising in the dimensionless equations (<ref>) are defined by :

\begin{equation}\label{eq:prhoadim}
Re=\frac{\tilde{u}_b\tL_y}{\tilde{\nu}},\quad Sc=\frac{\tilde{\nu}}{\tilde{\kappa}},
\end{equation}

where $\tilde{\nu}$ is the kinematic viscosity and $\tilde{\kappa}$ denotes the molecular diffusivity of the chemical specie producing the density difference.

§ APPLICATION OF PINNS ON SYNTHETICAL DATA

§.§ synthetic data for the lock-exchange problem

An accurate representation of the interface between the miscible fluids is crucial in gravity currents simulations, which requires high-order numerical methods to compute steep gradients in the vicinity of the interface. To this end, the governing equations (<ref>) are solved with the code Nek5000 ([Fischer et al., 2008]) that has been for example successfully employed in numerical investigations by  [Özgökmen et al., 2004]. The discretization scheme in Nek5000 is based on the spectral element method proposed by  [Patera, 1984] with exponential convergence in space and $3^\text{rd}$-order timestepping scheme.

We consider the dimensionless computational domain $\Omega = [0,6.0] \times [0,1.0]$ and $T = [0,10]$ with $Re=5500$. While for liquids such as salt water $Sc\simeq 700$, we consider here $Sc=1$. This assumption is commonly made to reduce computational costs and has low influence on the gravity current front as shown by  [Marshall et al., 2021].
A free-slip condition (no normal flow) is applied at $y=1.0$ and no slip-conditions are applied on the remaining domain boundaries. Grid independent results were obtain on a computational mesh of $90\times45$ elements where unknowns are represented by $7^\text{th}$-order Lagrange interpolating polynomials, and with a fixed timestep $\Delta t = 2.10^{-2}$ giving a Courant–Friedrichs–Lewy number less than $0.5$. We show in Fig. <ref> the time evolution of the density field obtain with Nek5000.

Heatmap of the density field $\rho$  for different time steps during the flow.

We build from this simulation a training dataset $\mathcal{D}$ where the collocation points are selected in the sub domain : $\Omega_{obs} \times T_{obs}= [0.0, 4.5] \times [0,1.0] \times [4.0,7.5] \subset \Omega \times T $,  which is discretized into a uniform Cartesian grid of $N_x \times N_y \times N_t$ points , where $N_{x}=750$, $N_{y}= 150$ and $N_{t} = 400$.

Finally, to mimic experimental measurements of LAT, we extract from this simulation the observed quantity of the density $\rho$ on each collocation points, and build our training dataset $\mathcal{D} = (\textbf{\textrm{x}}_i,t_i,\rho_i)_{i=1 \dots N_xN_yN_t}$ while the other fields $u,~v,~p$ are used for validation purpose.

§.§ Noise robustness

In this section, noises of varying magnitudes were added to $\mathcal{D}$.
We investigate 4 different cases :

* Case 1 :  Original dataset $\mathcal{D}$.

* Case 2 :  We add a gaussian noise to $\rho$ field with standard deviation $\sigma = 5\%$ (experimental precision of LAT).

* Case 3 :  same as Case 2 but with $\sigma =25\%$.

* Case 4 :  Same as case 3 but in addition, we remove a part of the information in all the snapshots. (we draw 200 squares of $10 \times 10$ pixels randomly from Case 3 dataset

A visualization of training data for $\rho$ for cases 1-4 is available Fig. <ref>

Density field at t = 10 for the different cases studied. In case number 4, 200 squares of size 10 by 10 are randomly selected in the mesh $\Omega_{obs}$

For all cases, the PINN loss function is expressed as:

\begin{align}
& \mathcal{L} = \lambda\mathcal{L}_{\textrm{data}} + \mathcal{L}_{f}\\
& \mathcal{L}_{f} = \dfrac{1}{N_f}\sum_{i=1}^{N_{f}}||f(\textrm{\textbf{x}}^{f}_{i},t^{f}_{i})||^{2} \\
& \mathcal{L}_{\textrm{data}} = \dfrac{1}{N_{\textrm{data}}} \sum_{i=1}^{N_{\textrm{data}}} ||\rho_{\mathcal{N}}(\textrm{\textbf{x}}^{\textrm{data}}_{i},t^{\textrm{data}}_{i})-\rho(\textrm{\textbf{x}}^{\textrm{data}}_{i},t^{\textrm{data}}_{i})||^{2}  \label{eq:Loss2D1:data}
\end{align}

The neural networks $\mathcal{N}$ take as input $(\textbf{\textrm{x}},t) = (x,y,t)$ and outputs the prediction $\boldsymbol{v}_{\mathcal{N}}(\textbf{\textrm{x}},t) = (\rho_{\mathcal{N}},u_{\mathcal{N}},v_{\mathcal{N}},p_{\mathcal{N}})$. The term $\mathcal{L}_{f}$ penalizes the governing equations  (<ref>) and $\mathcal{L}_{\textrm{data}}$ penalizes mismatch between network prediction for $\rho_{\mathcal{N}}$ and observed $\rho$.  The parameter $\lambda$ is a weighting  coefficient used to reinforce the influence of <ref>, in practice we set $\lambda = 500$.

The neural network is composed of 8 layers of 190 neurons each. The stochastic Adam optimizer [Kingma and Ba, 2014] is employed for solving <ref> : Glorot normal initializer [Glorot and Bengio, 2010] is employed to initialize the biases and the weights that are computed iteratively with a gradient update <ref> on a subset $\mathcal{B} \subset \mathcal{D}$ called a batch. One training round over $\mathcal{D}$ is called an epoch.  Training is done for 100 epochs with a batch size of 4096. We use a decreasing learning rate schedule. The learning rate starts at 1e-3 and ends at 1e-5. swish  activation function [Ramachandran et al., 2017] is selected for all layers except the last one.  Finally, training is done with two GPU A6000 in parallel.
The choice of these hyperparameters provides a goodsatisfactory compromise between quality of the results and computation time.

§.§ Results on synthetical data

For quantitative results, we define the following relative L2-norm :

\begin{equation}\label{eq:L2gene}
\epsilon_V =\frac{100}{\sup_{\bfx \in  \Omega_{obs} \times T_{obs}} \left| V_\bfx \right|}\sqrt{\frac{\sum_{\bfx\in\bfX}	\left|V_{\mathcal{N}}(\bfx) - V_\bfx\right|^{2}}{\left|\bfX\right|}}
\end{equation}

where $V \in (\rho,u,v,w,p)$ and $\bfX \subseteq \Omega_{obs} \times T_{obs}  $ defines the domain where the error is computed. We choose $\sup_{\bfx \in \Omega_{obs} \times T_{obs}} \left| V_\bfx \right|$ as reference value to avoid the division near to zero issue reported in [Delcey et al., 2023].

In Fig. <ref>, we can see that the models for cases 1,2 and 3 are very efficient and accurate for the inference of all fields. Case 4 also performs great mostly everywhere, but we can see a big drop in performance near $x = 1.5$ which may be due to a lack of training data in this area. For all the cases, the results for $u$ and $v$ close to $x = 4.49$ are less accurate, but the steep variations near the front are well inferred.

Fig. <ref>  shows the evolution in time of the error in $\Omega_{obs}$ for the fields $\rho,~u,~v$. We can see the $L_2$-errors for the density and pressure field are almost constant in time except for density in Case 4. For all cases, the accuracy of the predictions $u$ and $v$ improves with time. Since gravity currents propagate from left to right, at the initial stages of the flow, a large portion of the area covered by the flow contains zero values that might not be useful for training the neural network, making it more challenging to effectively learn and approach the velocity. However, as the gravity current progresses and reaches the end of the flow, the entire current is present in $\Omega_{obs}$, and more useful values such as non zero variations of density can be utilized by the neural network, possibly leading to improved performance.

Table 1 illustrates the global error for Cases 1 to 4 between the density, velocity and pressure fields predicted by PINN, and the corresponding ground thruth. The first three cases have highly accurate predictions , indicating that the model is robust to noise. Despite the higher errors in the fourth case, the predictions still fall within an acceptable range, suggesting that the model can handle incomplete data with high noise level.

These results are very promising and the robustness of PINNs to noise in the 2D lock exchange configuration suggests the potential of this framework for real experimental data.

In the following section, we apply the PINN model to experimental data obtained by LAT.

Horizontal profiles comparison. Profiles of $(\rho,u,v)$ reconstructed by PINNs for Cases 1-4 along the x-direction at y = 0.24 and t = 7.5 against the reference solution (first column).
Corresponding location of the profile in $\Omega_{obs}$ (second column).

Evolution of the L2-relative error through time for density, velocity and pressure fields.

Case     $\epsilon_\rho$     $\epsilon_u$      $\epsilon_v$       $\epsilon_p$

Case $n^{o}$1     0.035%     1.17%      1.81%      1.03%

Case $n^{o}$2     0.085%     1.29%       2.41%      1.07%

Case $n^{o}$3    0.53%     2.65%     5.30%    2.48%

Case $n^{o}$4     7.98%     4.93%     5.12%     4.48%

Relative $L_2$-norm errors of the fields reconstructed by PINN computed on the spatio-temporal domain $\Omega_{obs} \times T_{obs}$.

§ APPLICATION OF PINNS TO LAT EXPERIMENTAL MEASUREMENTS

§.§ Experimental setup

The experimental apparatus is summarized in Fig.  <ref>. Experiments are conducted in a $150$ cm long, $10$ cm wide and $35$ cm high plexiglas tank. After the tank is filled with freshwater to an height $25$ cm, a gate with insulating sides is inserted at 20 cm from the left end of the tank.

Salt is dissolved in the lock to reach a homogeneous density $\tilde{\rho}_1= 1001.2~\textrm{kg.m}^{-3}$ slightly larger than the freshwater density $\tilde{\rho}_2=999.8~\textrm{kg.m}^{-3}$ which gives us $\Delta \rho = 1.4~\textrm{kg.m}^{-3}$, $Re= 14000$ and $Sc=1000$ according to the definitions in section (<ref>).

Measurements are performed using two cameras and two light sources.
The PIV and LAT cameras are respectively mounted with a $180$ mm f2-8 and $25$ mm lens.
The cameras are mounted on tripods and placed at a distance of $250$ cm from the tank. Images are recorded at a frame rate of $10$ fps and the spatial resolution of the raw images is 8Mpx for both cameras.

The PIV (particle image velocimetry [Raffel et al., 2007]) camera is synchronized with a vertical LASER sheet that illuminates the central section of the tank, providing two-dimensional vertical velocity field measurements as commonly performed in stratified flows.

The LAT (light attenuation technique) camera is synchronized with a white LED panel place behind the tank. The latter technique allows to measure high resolution cross-averaged density fields by adding a known concentration of tracer in the dense compartment. In the present experiment, a red food dye is used as a tracer for salinity and hence density. The reader may refer to [Dossmann et al., 2016] for a detailed description of the LAT protocol.

The two cameras and light sources are computer
-controlled, which permits to synchronize velocity and density measurements. As the LED panel interferes with the PIV measurements, it is required to trigger the PIV and LAT sources in a slightly phase-shifted manner. Measurements are latter recombined on an identical timestep during the PIV post-process. A common spatial origin is used to match the recorded velocity and density field.

(a) Sketch of the experimental setup, including the tank,
barrier, light panel, fluids and cameras 1 (LAT) and 2 (PIV). Combination of LAT and PIV measurement devices with $\tilde{L_y}$ = 25 cm,  $\tilde{L_x}$ = 150 cm and l = 10 cm.  The cameras are located at a distance of $D$=250 mm from the tank. (b) Time scale used to represent, in chronological order, the duration of PIV time with 8000 $\mu s$ correlation time between two images and LAT time with 2500 $\mu s$ exposure time.

As the studied flow evolves in 3D, using the canonical PINN model would imply forcing a neural network to impose the Navier-Stokes equations in 3D and having available measurements of the density field in the whole tank which is prohibitive in our case.

So in the following, we assume that the flow is homogeneous in the z-direction, this assumption is commonly made for the lock release configuration [La Rocca and Bateman, 2010]. As a consequence, we reduce the problem size and consider only the 2D Navier-Stokes equations <ref> and, under this assumption, the spanwise-averaged density $\overline{\rho}$ measured by LAT is assimilated to 2D density measurements $\rho$.
In Fig. <ref>, we can see a sample of the fields measured by the combination of LAT and PIV.

§.§ Training data

The data are acquired on the spatial and temporal dimensionless domain $\Omega_{exp} \times T_{exp} = [3.0,4.65] \times [0.05,0.5] \times [19.5,21.5]$.

this domain is discretized into a uniform grid of $N_x \times N_y \times N_t$ points , where $N_{x}=1050$, $N_{y}= 300$ and $N_{t} = 90$.

measures of density at $t = 20.45$ and two components of velocity at $t = 20.80. ~ \rho $ is measured by LAT and $(u, v)$ are measured by PIV in the middle plane of the tank

Our objective is to infer the hydrodynamic field in the middle plane of the domain given LAT measurements of density. We keep the available PIV measurements for validation purpose.

We apply canonical PINN method with LAT experimental measures as observation data in the domain $\Omega_{exp} \times T_{exp}$ with Navier-Stokes equation <ref> and loss fonction <ref>.

The neural network is composed of 8 layers of 250 neurons each which, in our case is the structures which gives the best results. Training is done for 150 epochs with a batch size of 4096. We use Adam optimizer with a decreasing learning rate schedule. The learning rate starts at 5e-4 and ends at 1e-5. Finally we use $swish$ activation function on each layers transition, except for the last layers where no activation function is selected.

Comparison between measured and predicted density and velocity for model $\mathcal{N}$ at $t=20.8$. First line refers to $\rho$ obtained by LAT and $(u,v)$ measured by PIV . Second line refers to corresponding fields inferred by PINN , and third line corresponds to horizontal profiles along x-axis for $y= 0.11$ (red), $y = 0.28$ (blue) and $y = 0.48 $ (green) for both inferred and measured fields.

§.§ Results on experimental data

In Fig <ref>, we can see a comparison between predicted $\rho$ $u$ and $v$ fields and the corresponding measured fields by LAT/PIV. Since the model is trained on density, it is not surprising that the predicted density values are highly accurate, but we can notice that PINN has a smoothing effect on the density field. What is the most interesting here is the comparison on the velocity fields. Indeed, the NN has succeeded without prior information on $u$ and $v$ to reproduce the right amplitudes of values as well as the location of the current front. On the other hand, there is a bias in the values for $u$ and the small variations above the current (blue profile) are poorly captured. Moreover, for $v$, the values on the right of the front seem to be outliers.

After a thorough study of the problem parameters, we came to the conclusion that the residual related to the continuity equation (<ref>) in 2D degraded the results. So, to account for the spanwise variations of the flow, we propose a modification of the continuity equation (<ref>) :

\begin{equation}
\label{eq:cont1}
~~~~~~~~~~~~~~~~~~~~\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \xi = 0
\end{equation}

where $\xi$ which depends on $(x,y,t) \in \Omega_{exp} \times T_{exp}$ is a new output of the neural network.
Thus, the new model $\mathcal{M}$ take as input $(\textbf{\textrm{x}},t) = (x,y,t)$ and outputs the prediction $\boldsymbol{v}_{\mathcal{M}}(\textbf{\textrm{x}},t) = (\rho_{\mathcal{M}},u_{\mathcal{M}},v_{\mathcal{M}},p_{\mathcal{M}},\xi_{\mathcal{M}})$.
Since the 3D Navier-Stokes equations involve $\dfrac{\partial w}{\partial z}$ in the continuity equation, we expect $\xi_{\mathcal{M}}$ to converge during training to the restriction of $\dfrac{\partial w}{\partial z}$ to the PIV plane.
In the following, we denote by $\mathcal{N}$ or $\mathcal{M}$ either if the model was trained with the equation (<ref>) or with (<ref>).

Comparison between measured and predicted density and velocity for model $\mathcal{M}$ at $t=20.8$.  at $t = 20.8$. First line refers to $\rho$ obtained by LAT and $(u,v)$ measured by PIV . Second line refers to corresponding fields inferred by PINN , and third line corresponds to horizontal profiles along x-axis for $y= 0.11$ (red), $y = 0.28$ (blue) and $y = 0.48 $ (green) for both inferred and measured fields.

The modification made to the equation has led to significant improvements in the accuracy of the predictions.
The results Fig. <ref> show the efficiency of this new formulation of the problem.
For the $u$ component of velocity, although a small bias can still be observed on the values above the stream, the predicted quantities inside the stream and near the blue zone are suitable. The predicted  profiles for the $v$ component of velocity in the blue and green zones are now accurate, and the model is able to capture steep variations, such as the large peak at $x=4.3$, with good accuracy.

Model     $\epsilon_\rho$     $\epsilon_u$      $\epsilon_v$

$\mathcal{N}$     1.72%     20.90%       31.5%

$\mathcal{M}$       1.35%     14.6%     12.65%

Relative $L_2$-norm errors on the spatio-temporal domain $\Omega_{exp} \times T_{exp}$ between reconstructed and measured fields. model name is $\mathcal{N}$ or $\mathcal{M}$ either if training was done with standard continuity equation (<ref>) or variation (<ref>)

The table <ref> shows the overall improvement in predictions brought by the simple change in the continuity equation <ref>. The model achieve great results on the global space $\Omega_{exp} \times T_{exp}$.

Model $\mathcal{M}$ doesn't take into account all the 3D effects, but correction (<ref>) makes it more precise.
The surrogate model is fully differentiable and enables the computation of various quantities of interest, including vorticity, pressure, and energy fluxes and budgets.
In the following, we examine the pressure field inferred by analyzing energy fluxes and the residual of the energy conservation equation given by :

\begin{equation}
\label{eq:energy}
\begin{aligned}
& \left| e_{5} \right|= \bigg| \frac{\partial}{\partial t}\left( \frac{|\boldsymbol{u}|^{2}}{2}\right)  + \nabla \cdot \left[ \left( \dfrac{|\boldsymbol{u}|^{2}}{2} + p\right)\boldsymbol{u}\right]  \\
& -\nabla \cdot (\dfrac{1}{\textrm{Re}} \nabla \boldsymbol{u} \cdot \boldsymbol{u} ) + \dfrac{1}{\textrm{Re}} \nabla\boldsymbol{u} :  \nabla\boldsymbol{u} + \rho\boldsymbol{u} \cdot \boldsymbol{e_y} \bigg|.
\end{aligned}
\end{equation}

In Fig. <ref>, we represent the quantity $J = p \boldsymbol{u}$ where $\boldsymbol{u} = (u,v)^{T}$, and respectively the isovalues of the separated components $pv$, $pu$ given by $\mathcal{M}$ at three different times during the flow. These quantities provide information about the energy flux within the gravity current.

We clearly evidence a zone where the current transfers energy to the surrounding fluid, which are the zone where $pv$ reach its maximum values in the right of the current front. In the same way, one identifies a zone where the current captures energy from the ambient fluid which are the zone of circulation above the current where $ pu $ reach its minimum values. We can notice a very sharp transition between these two zones (red dotted line) where $pv$ changes sign. This criterion is consistent through time and represents a quantitative criterium to divide the density field in two parts, the front and the head of the current. In the last raw of Fig. <ref> vorticity isovalues inferred $\omega$ are plotted and superposed to the energy flux vector field. We observe that the transition between the two previously identified zones is also associated with the appearance of re-circulation rolls where the mixing becomes efficient.

Although we have no direct experimental measurement of the pressure field for validation, it is interesting to discuss its relevance in regards to the energy conservation equation (see eq (<ref>)).
For that, we plot in Fig. <ref> the magnitude of energy residual $\lvert e_{5} \rvert$ defined as:

We find that $\lvert e_{5} \rvert$ is around $3 \times 10^{-2}$. This error is small from a physical point of view, especially considering the experimental errors of measurements for such unsteady flows but is significant from a mathematical point of view. This can be explained by the nature of the physical model used to model the flow here.
Indeed, we use the 2D equations to train the neural network. This approximation is justified in view of the mainly two-dimensional character of the flow but the small fluctuations of the transverse component of the velocity can carry a non-negligible fraction of the energy which would explain a significantly high value of the residuals $e_{5}$. An improvement in the future would be to develop the same type of approach but using 3D flow reconstruction methods from 2D measurements (see [Delcey et al., 2023]). The implementation of the 3D equations would then allow to significantly reduce the error on the energy transport equation.

Energy flux and vorticity. Isovalues of $\rho$ (first line), $pu, pv$ (seconde and third line) and vorticity $\omega$ (fourth line) at times : $t = 20.7, 21.14$ and $21.6$ in addition in lines 2,3,4, the vector field $p\boldsymbol{u}$ is added.

Time evolution of the residual of the energy equation (<ref>).

§ CONCLUSION

In this work, we quantitatively demonstrate that PINNs are able to solve inverse problems from experimental data. We first investigated the robustness of this model to noisy and gappy data. This study was made from synthetic data obtained by the CFD solver Nek5000 for the lock-exchange configuration, and allowed us to study the parametric aspect of the model. Then we trained a physics informed model from LAT measurements acquired on the lock-exchange configuration, the inferred velocity fields have been compared to PIV measurements made simultaneously on the same experiment. Despite the assumption of homogeneity in z-direction, the results are convincing and the model succeeds in predicting the pressure field. Based on the energy flux computed from the results, we discuss the structure of the gravity current and show the potential of Physics Informed Neural Networks for ocean modeling.

§ BIBLIOGRAPHY

[Perez-Diaz et al., 2018]
authorB. Perez-Diaz, authorPalomar, et al.,
titlePiv-plif characterization of nonconfined saline
density currents under different flow conditions,
journalJournal of Hydraulic Engineering
volume144 (year2018) pages04018063.

[Crimaldi, 2008]
authorJ. Crimaldi,
titlePlanar laser induced fluorescence in aqueous flows,
journalExperiments in fluids volume44
(year2008) pages851–863.

[Dossmann et al., 2016]
authorDossmann, et al.,
titleExperiments with mixing in stratified flow over a
topographic ridge,
journalJournal of Geophysical Research: Oceans
volume121 (year2016) pages6961–6977.

[Balasubramanian et al., 2018]
authorBalasubramanian, et al.,
titleEntrainment and mixing in lock-exchange gravity
currents using simultaneous velocity-density measurements,
journalPhysics of Fluids volume30
(year2018) pages056601.

[Tian et al., 2003]
authorTian, et al.,
titleA 3d lif system for turbulent buoyant jet flows,
journalExperiments in fluids volume35
(year2003) pages636–647.

[Raffel, 2015]
authorM. Raffel,
titleBackground-oriented schlieren (bos) techniques,
journalExperiments in Fluids volume56
(year2015) pages1–17.

[Elsinga et al., 2006]
authorG. E. Elsinga, et al.,
titleTomographic particle image velocimetry,
journalExperiments in fluids volume41
(year2006) pages933–947.

[Partridge et al., 2019]
authorPartridge, et al.,
titleA versatile scanning method for volumetric
measurements of velocity and density fields,
journalMeasurement Science and Technology
volume30 (year2019) pages055203.

[Nash et al., 2005]
authorNash, et al.,
titleEstimating internal wave energy fluxes in the ocean,
journalJournal of Atmospheric and Oceanic Technology
volume22 (year2005) pages1551–1570.

[Clark et al., 2010]
authorClark, et al.,
titleGeneration, propagation, and breaking of an internal
wave beam,
journalPhysics of Fluids volume22
(year2010). note076601.

[Allshouse et al., 2016]
authorAllshouse, et al.,
titleInternal wave pressure, velocity, and energy flux
from density perturbations,
journalPhysical Review Fluids volume1
(year2016) pages014301.

[Fourés et al., 2014]
authorFourés, et al.,
titleA data-assimilation method for reynolds-averaged
navier–stokes-driven mean flow reconstruction,
journalJournal of Fluid Mechanics volume759
(year2014) pages404–431.

[Mons et al., 2022]
authorV. Mons, et al.,
titleDense velocity, pressure and eulerian acceleration
fields from single-instant scattered velocities through navier–stokes-based
data assimilation,
journalMeasurement Science and Technology
volume33 (year2022) pages124004.

[Goodfellow et al., 2016]
authorI. Goodfellow, et al., titleDeep Learning,
publisherMIT Press, year2016.

[Brunton et al., 2020]
authorS. L. Brunton, et al.,
titleMachine learning for fluid mechanics,
journalAnnual review of fluid mechanics
volume52 (year2020) pages477–508.

[Raissi et al., 2020]
authorM. Raissi, authorA. Yazdani,
authorG. E. Karniadakis,
titleHidden fluid mechanics: Learning velocity and
pressure fields from flow visualizations,
journalScience volume367
(year2020) pages1026–1030.

[Mao, 2020]
titlePhysics-informed neural networks for high-speed
journalComputer Methods in Applied Mechanics and
Engineering volume360 (year2020)

[Cai et al., 2021]
authorCai, et al.,
titlePhysics-informed neural networks for heat transfer
journalJournal of Heat Transfer volume143

[Cuomo et al., 2022]
authorCuomo, et al.,
titleScientific machine learning through physics-informed
neural networks: Where we are and what's next,
journalJournal of Scientific Computing
volume92 (year2022).

[Cai et al., 2021]
authorCai, et al.,
titleFlow over an espresso cup: inferring 3-d velocity and
pressure fields from tomographic background oriented schlieren via
physics-informed neural networks,
journalJournal of Fluid Mechanics volume915

[Molnar et al., 2023]
authorJ. Molnar, et al.,
titleEstimating density, velocity, and pressure fields in
supersonic flows using physics-informed bos,
journalExperiments in Fluids volume64
(year2023) pages14.

[Cybenko, 1989]
authorG. Cybenko,
titleApproximation by superpositions of a sigmoidal
function. math cont sig syst (mcss) 2:303-314,
journalMathematics of Control, Signals, and Systems
volume2 (year1989) pages303–314.

[Higham et al., 2019]
authorHigham, et al.,
titleDeep learning: An introduction for applied
journalSiam review volume61
(year2019) pages860–891.

[Géron, 2019]
authorA. Géron, titleHands-on machine learning with
Scikit-Learn, Keras, and TensorFlow: Concepts, tools, and techniques to build
intelligent systems, publisherO'Reilly Media,

[Paszke et al., 2019]
authorPaszke, et al., titlePytorch: An imperative style,
high-performance deep learning library, year2019. <http://arxiv.org/abs/1912.01703>, notecite
arxiv:1912.01703Comment: 12 pages, 3 figures, NeurIPS 2019.

[Griewank et al., 1989]
authorA. Griewank, et al.,
titleOn automatic differentiation,
journalMathematical Programming: recent developments and
applications volume6 (year1989)

[Kingma and Ba, 2014]
authorD. P. Kingma, authorJ. Ba,
titleAdam: A method for stochastic optimization,
journalarXiv preprint arXiv:1412.6980

[Shin et al., 2004]
authorShin, et al.,
titleGravity currents produced by lock exchange,
journalJournal of Fluid Mechanics volume521
(year2004) pages1–34.

[Gray and Giorgini, 1976]
authorD. D. Gray, authorA. Giorgini,
titleThe validity of the boussinesq approximation for
liquids and gases,
journalInternational Journal of Heat and Mass Transfer
volume19 (year1976) pages545–551.

[Hartel et al., 2000]
authorHartel, et al.,
titleAnalysis and direct numerical simulation of the flow
at a gravity-current head. part 1. flow topology and front speed for slip and
no-slip boundaries,
journalJournal of Fluid Mechanics volume418
(year2000) pages189–212.

[Fischer et al., 2008]
authorP. F. Fischer, et al., titlenek5000 Web page,
year2008. noteHttp://nek5000.mcs.anl.gov.

[Özgökmen et al., 2004]
authorÖzgökmen, et al.,
titleThree-dimensional turbulent bottom density currents
from a high-order nonhydrostatic spectral element model,
journalJournal of Physical Oceanography
volume34 (year2004) pages2006–2026.

[Patera, 1984]
authorA. T. Patera,
titleA spectral element method for fluid dynamics: laminar
flow in a channel expansion,
journalJournal of computational Physics
volume54 (year1984) pages468–488.

[Marshall et al., 2021]
authorMarshall, et al.,
titleThe effect of schmidt number on gravity current
flows: The formation of large-scale three-dimensional structures,
journalPhysics of Fluids volume33
(year2021) pages106601.

[Glorot and Bengio, 2010]
authorX. Glorot, authorY. Bengio,
titleUnderstanding the difficulty of training deep
feedforward neural networks,
in: booktitleProceedings of the thirteenth international
conference on artificial intelligence and statistics,
organizationJMLR Workshop and Conference Proceedings,
year2010, pp. pages249–256.

[Ramachandran et al., 2017]
authorRamachandran, et al.,
titleSearching for activation functions,
journalarXiv preprint arXiv:1710.05941

[Delcey et al., 2023]
authorDelcey, et al.,
titlePhysics-informed neural networks for gravity
currents reconstruction from limited data,
journalPhysics of Fluids volume35
(year2023). note027124.

[Raffel et al., 2007]
authorM. Raffel, authorC. Willert,
authors. wereley, authorJ. Kompenhans,
titleParticle Image Velocimetry: A Practical Guide, Experimental
Fluid Mechanics, publisherSpringer Berlin Heidelberg,
year2007. <https://books.google.fr/books?id=WH5w2TvK2bQC>.

[La Rocca and Bateman, 2010]
authorM. La Rocca, authorA. Bateman,
titleExperimental and Theoretical Modelling of 3D Gravity
Currents, year2010. 10.5772/13251.