# Matter-wave interference with particles selected from
a molecular library with masses exceeding 10 000 amu

Sandra Eibenberger  University of Vienna, Faculty of Physics, VCQ, QuNaBioS,
Boltzmanngasse 5, 1090 Vienna (Austria)  Stefan Gerlich  University of Vienna,
Faculty of Physics, VCQ, QuNaBioS, Boltzmanngasse 5, 1090 Vienna (Austria)
Markus Arndt<EMAIL_ADDRESS>University of Vienna, Faculty of
Physics, VCQ, QuNaBioS, Boltzmanngasse 5, 1090 Vienna (Austria)  Marcel Mayor
<EMAIL_ADDRESS>Department of Chemistry, University of Basel, St.
Johannsring 19, 4056 Basel (Switzerland),
Karlsruhe Institute of Technology (KIT), Institute of Nanotechnology, P.O. Box
3640, 76021 Karlsruhe (Germany)  Jens Tüxen  Department of Chemistry,
University of Basel, St. Johannsring 19, 4056 Basel (Switzerland)

###### Abstract

The quantum superposition principle, a key distinction between quantum physics
and classical mechanics, is often perceived as a philosophical challenge to
our concepts of reality, locality or space-time since it contrasts our
intuitive expectations with experimental observations on isolated quantum
systems. While we are used to associating the notion of localization with
massive bodies, quantum physics teaches us that every individual object is
associated with a wave function that may eventually delocalize by far more
than the body’s own extension. Numerous experiments have verified this concept
at the microscopic scale but intuition wavers when it comes to delocalization
experiments with complex objects. While quantum science is the uncontested
ideal of a physics theory, one may ask if the superposition principle can
persist on all complexity scales. This motivates matter-wave diffraction and
interference studies with large compounds in a three-grating interferometer
configuration which also necessitates the preparation of high-mass
nanoparticle beams at low velocities. Here we demonstrate how synthetic
chemistry allows us to prepare libraries of fluorous porphyrins which can be
tailored to exhibit high mass, good thermal stability and relatively low
polarizability, which allows us to form slow thermal beams of these high-mass
compounds, which can be detected in electron ionization mass spectrometry. We
present successful superposition experiments with selected species from these
molecular libraries in a quantum interferometer, which utilizes the
diffraction of matter waves at an optical phase grating. We observe high-
contrast quantum fringe patterns with molecules exceeding a mass of 10 000 amu
and 810 atoms in a single particle.

††preprint: PCCP, DOI: 10.1039/C3CP51500A

##  I Introduction

Quantum physics has long been regarded as the science of small things , but
experimental progress throughout the last two decades has led to the insight
that it can also be observable for mesoscopic or even macroscopic objects.
This applies for instance to the superposition of macroscopic numbers of
electrons in superconducting quantum devices Mooij et al. (1999), the
realization of large quantum degenerate atomic clouds in Bose-Einstein
condensates Anderson et al. (1995), the cooling of micromechanical oscillators
to their mechanical ground state.

Quantum superposition studies with complex molecules Arndt et al. (1999)
became possible with the advent of new matter wave interferometers Brezger et
al. (2003); Gerlich et al. (2007); Haslinger et al. (2013) and techniques for
slow macromolecular beams Deachapunya et al. (2008). These interferometers
were for instance practically used to enable quantum enhanced measurements of
internal molecular properties. The quantum fringe shift of a molecular
interference pattern in the presence of external electric fields resulted in
information for example on electric polarizabilities Berninger et al. (2007),
dipole moments Eibenberger et al. (2011), or configuration changes Gring et
al. (2010). Molecule interferometry can complement mass spectrometry Gerlich
et al. (2008) and help to distinguish constitutional isomers Tüxen et al.
(2010). In addition to their applications in chemistry, quantum interference
experiments with massive molecules currently set the strongest bound on
certain models that challenge the linearity of quantum mechanics Bassi et al.
(2013).

Further exploration of the frontiers of de Broglie coherence now profits from
new capabilities in tailoring molecular properties to the needs of quantum
optics. Our quantum experiment dictates the design of the molecules and the
challenges increase with the number of atoms involved. In order to realize a
molecular beam of sufficient intensity, the model compounds must be volatile,
thermally stable and accessible in quantities of several hundred milligrams.
Moreover, in order to minimize absorption at the wavelength of the optical
diffraction (see below) they need to feature low absorption but sufficient
polarizability at 532 nm. In reply to these needs, a dendritic library concept
is particularly appealing since it can be scaled up to complex particles, once
a suitable candidate has been identified.

To meet these requirements we have functionalized organic chromophores with
extended perfluoroalkyl chains. Such compounds show low inter-molecular
binding and relatively high vapor pressures Stock et al. (2004); Krusic et al.
(2005). They possess strong intra-molecular bonds and therefore sufficient
thermal stability. In addition we start with a porphyrin core which is
compatible with the required optical and electronic properties Meot-Ner et al.
(1973). In the past, monodisperse fluorous porphyrins were generated by
substituting the four para-fluorine substituents of tetrakis
pentafluorophenylporphyrin (TPPF20) by dendritically branched fluorous
moieties Tüxen et al. (2011). With this approach, molecules composed of 430
atoms were successfully synthesized and applied in quantum interference
experiments Gerlich et al. (2011).

With increasing complexity it becomes more challenging to purify monodisperse
particles in sufficient amounts. Here we profit from the fact that our
interferometer arrangement allows us to work with compound mixtures since each
molecule interferes only with itself. By substituting some of the twenty
fluorine atoms of TPPF20 with a branched, terminally perfluorinated alkylthiol
(1), we obtain a mixture of compounds with molecular masses that differ
exactly by an integer multiple of a particular value as molecular library. The
molecular beam density is sufficiently low for the molecules not to interact
with each other and the individual library compounds can be mass-specifically
detected in a quadrupole mass spectrometer (QMS Extrel, 16 000 amu).

Our synthetic approach is based on the fact that pentafluorophenyl moieties
can be used to attach up to five polyfluoroalkyl substituents in nucleophilic
aromatic substitution reactions. Substitutions at TPPF20 with its 20
potentially reactive fluorine substituents lead to a molecular library of
derivatives with a varying number of fluorous side chains.

##  II Results and Discussion

We used sodium hydride as a base, microwave radiation as a heating source and
diethylene glycol dimethyl ether (diglyme) as fluorophilic solvent. TPPF20,
and a large excess of the thiol 1 (60 equivalents) and sodium hydride in
diglyme were heated in a sealed microwave vial to 220 °C for 5 minutes.111
Synthetic protocol and analytical data of the porphyrin libraries L: General
Remarks: All commercially available starting materials were of reagent grade
and used as received. Microwave reactions were carried out in an Initiator 8
(400 W) from Biotage. Glass coated magnetic stirring bars were used during the
reactions. The solvents for the extractions were of technical grade and
distilled prior to use. Matrix Assisted Laser Desorption Ionization Time of
Flight (MALDI-ToF) mass spectra were performed on an Applied Bio Systems
Voyager-De™ Pro mass spectrometer or a Bruker microflex mass spectrometer.
Significant signals are given in mass units per charge (m/z) and the relative
intensities are given in brackets.
Porphyrin library L: The thiol 1 was synthesized in seven reaction steps in an
overall yield of 70% as reported elsewhere Tüxen et al. (2011).
5,10,15,20-Tetrakis(pentafluoro-phenyl)-porphyrin (TPPF20, 4.0 mg, 4.10
$\rm{mu}$mol, 1.0 eq.), thiol 1 (193 mg, 246 $\rm{mu}$mol, 60 eq.) and sodium
hydride (60% dispersion in mineral oil, 14.8 mg, 369 $\rm{mu}$mol, 90 eq.)
were added to diglyme (4 mL) in a microwave vial. The sealed tube was heated
under microwave irradiation to 220 °C for 5 minutes. After cooling to room
temperature the reaction mixture was quenched with water and subsequently
extracted with diethylether. The organic layer was washed with brine and
water, dried over sodium sulfate and evaporated to dryness. The resulting
product mixture (183 mg) was analyzed by MALDI-ToF mass spectrometry. MS
(MALDI-ToF, m/z): 12 403 (29%), 11 645 (52%), 10 884 (100%), 10 121 (78%), 9
339 (15%), 8 597 (7%). After aqueous workup the resulting mixture was analyzed
by MALDI-ToF mass spectrometry (Figure 1) and subsequently used in our quantum
interference experiments without further purification. We find up to 15
substituted fluorous thiol chains reaching to a molecular weight well beyond
10 000 amu.

Figure 1: Synthetic scheme and MALDI-ToF mass spectrum of the fluorous
porphyrin library textitL. High-mass matter-wave experiments were performed
with component L12 of the library L. This structure is composed of 810 atoms
and has a nominal molecular weight of 10 123 amu.

In order to study the delocalized quantum wave nature of compounds in the
fluorous library we utilize a Kapitza-Dirac-Talbot-Lau interferometer (KDTLI),
which has already proven to be a viable tool with good mass scalability in
earlier studies Gerlich et al. (2007, 2011). The interferometer is sketched in
Figure 2: a molecular beam is created by thermal evaporation of the entire
library in a Knudsen cell. The mixture traverses three gratings G1, G2, and
G3, which all have the same period of d $\cong$ 266 nm. The molecules first
pass the transmission grating G1, a SiNx mask with a slit opening of
$s\approx$ 110 nm, where each molecule is spatially confined to impose the
required spatial coherence by virtue of Heisenberg’s uncertainty principle
Nairz et al. (2002). This is sufficient for the emerging quantum wavelets to
cover several nodes of the optical phase grating G2, 10.5 cm further
downstream. The standing light wave G2 is produced by retro-reflection of a
green laser ($\lambda_{\rm{L}}$= 532 nm) at a plane mirror.

Figure 2: KDTL interferometer setup: The molecules are evaporated in a
furnace. Three height delimiters, D1 \- D1, define the particle velocity by
selecting a flight parabola in the gravitational field. The interferometer
consists of three gratings with identical periods of $d$ = 266 nm. G1 and G3
are SiNx gratings, whereas G2 is a standing light wave which is produced by
retro-reflection of a green laser at a plane mirror. A phase modulation
$\Phi\propto(\alpha\cdot P)/(v\cdot\omega_{y})$ is imprinted onto the
molecular matter wave via the optical dipole force which is exerted by the
light grating onto the molecular optical polarizability $\alpha_{opt}$. Here
$P$ is the laser power, $v$ the molecular velocity,
$\omega_{x}\cong\rm{18\,\mu m}$ and $\omega_{y}\cong\rm{945\,\mu m}$ the
Gaussian laser beam waists. The transmitted molecules are detected using
electron ionization quadrupole mass spectrometry after their passage through
G3, which can be shifted along the z-axis to sample the interference pattern.

When the molecular matter wave traverses the standing light wave, the dipole
interaction between the electric field of power $P$ and the molecular optical
polarizability $\alpha_{\rm{opt}}$ entails a periodic phase modulation
$\Phi=\Phi_{0}\cdot\sin^{2}{(2\pi z/\lambda_{L})}$ with the maximum phase
shift $\Phi_{0}=8\sqrt{2\pi}\alpha_{opt}P/(\hbar c\omega_{y}v_{x})$. Here
$v_{x}$ is the forward directed molecular velocity, $\omega_{y}\cong$ 945
$\rm{\mu}$m is the Gaussian laser beam waist along the grating slits and $z$
is the coordinate along the laser beam.

Interference of the molecular wavelets behind G2 leads to a molecular density
pattern of the same period $d$ in front of the third grating. As long as we
can neglect photon absorption by the molecules in the standing light wave G2
acts effectively as a pure phase grating. The expected fringe visibility $V$
is then given by Hornberger et al. (2009) $V=2(\sin(\pi f)/\pi
f)^{2}J_{2}(-sgn(\Phi_{0}\sin(\pi L/L_{T}))\Phi_{0}\sin(\pi L/L_{T}))$. Here
$f$ designates the grating open fraction, i.e. the ratio between the open slit
width $s$ and the grating period $d$. $J_{2}$ is the Bessel function of second
order, $L$ the separation of the gratings, $L_{T}=d^{2}/\lambda_{dB}$ the
Talbot length and $\lambda_{dB}=h/(mv)$ the de Broglie wavelength with $h$ as
Planck s quantum of action, $m$ the molecular mass and $v$ the modulus of its
velocity.

G3 is again a SiNx structure and lends spatial resolution to the detector. The
interferogram is sampled by tracing the transmitted particle beam intensity as
a function of the lateral ($z$) position of G3 and the mass selection is
performed in the mass spectrometer behind this stage.

Talbot-Lau interferometers Clauser and Li (1997) offer the important advantage
over simple grating diffraction that the required grating period $d$ only
weakly depends on the molecular de Broglie wavelength:
$d\propto\sqrt{\lambda_{dB}}$. The setup accepts a wide range of velocities
and low initial spatial coherence. This facilitates the use of dilute thermal
molecular beams. Diffraction at the standing light wave G2 avoids the
dephasing caused by the van der Waals interaction between the molecules and a
dielectric wall. This is indeed present in G1 and G3 but can be neglected
there since the molecular momentum distribution at G1 is wider than that
caused by the grating and any phase shift in G3 is irrelevant if we are only
interested in counting the particles that reach the mass spectrometer.

Indistinguishability in all degrees of freedom is the basis for quantum
interference Dirac (1958) and naturally given if a single molecule interferes
only with itself. We only have to make sure that every molecule contributes to
the final pattern in a similar way, which is true for all members of the
library at about the same mass, independent of their internal state.

Differences between various molecules, such as their isotopic distribution or
the addition of a single atom are still acceptable. The KDTLI can tolerate a
wavelength distribution of $\Delta\lambda_{dB}/\lambda_{dB}\leq\rm{20\,\%}$
and still produce a quantum fringe visibility in excess of the classical
threshold.

We here present quantum interference collected at the mass of one specific
library compound, particularly for L12$=\rm{C_{284}H_{190}F_{320}N_{4}S_{12}}$
which has 12 fluorous side chains, a mass of 10 123 amu and 810 atoms bound in
a single hot nanoparticle.

All molecules of the library were evaporated at a temperature of about 600 K.
We selected the velocity class around $v$ = 85 m/s ($\Delta v_{\rm{FWHM}}$ =
30 m/s) corresponding to a most probable de Broglie wavelength of
approximately 500 fm. This is about four orders of magnitude smaller than the
diameter of each individual molecule. We detected the signal by electron
ionization quadrupole mass spectrometry. During the interference measurements
the mass filter was set to the target mass of L12 and only this compound
contributed to the collected interference pattern.

The molecular beam was dilute enough to prevent classical interactions between
any two molecules within the interferometer. Given that 80 mg of library L
molecules were evaporated in 45 minutes, we estimate a flux at the source exit
of $2\times 10^{15}$ particles per second. Including the acceptance angle of
the instrument, the velocity selection as well as the grating transmission we
estimate a molecular density inside the interferometer of 30 mm-3. This
corresponds to a mean particle distance of about 300 $\rm{\mu}$m which is
sufficient to exclude interactions with other neutral molecules in the beam.

The average flight time of a molecule through the tightly focused standing
light wave amounts to about 400 ns, i.e. much longer than the time scale of
molecular vibrations (10-14 \- 10-12 s) and rotations ( 10-10 s). Therefore,
the mean scalar polarizability governs the interaction with the standing light
wave although the optical polarizability is generally described by a tensor.
Thermal averaging also occurs for the orientation of any possibly existing
molecular electric dipole moment Eibenberger et al. (2011). The internal
molecular states are decoupled from de Broglie interference as long as we
exclude effects of collisional or thermal decoherence Hackermüller et al.
(2003, 2004) or external force fields Berninger et al. (2007).

The thermal mixture of internal states is another reason why two-particle
interference, i.e. mutual coherence of two macromolecules, is excluded in our
experiments. The chances of finding two of them in the same indistinguishable
set of all internal states electronic, vibrational and rotational levels,
configuration, orientation and spin is vanishingly small.

Figure 3: (a) Quantum interference pattern of L12 recorded at a laser power of
$P\cong\rm{1\,W}$. The circles represent the experimental signal $s$ as a
function of the position $z$ of the third grating. The solid line is a
sinusoidal fit to the data, with a quantum fringe visibility of $V$ = 33(2) %.
The shaded area represents the background signal of the detector. A classical
picture predicts a visibility of only 8% for the same experimental parameters.
b) Measured fringe visibility $V$ as a function of the diffracting laser power
$P$. The expected contrast according to the quantum and the classical model
are plotted as the blue and red lines, respectively Hornberger et al. (2009).
The dashed blue lines correspond to the expected quantum contrast when the
mean velocity is increased (reduced) by 5 ms-1.

In Figure 3(a) we show a high contrast quantum interference pattern of L12. In
contrast to far-field diffraction where the fringe separation is governed by
the de Broglie wave-length of the transmitted molecules,Nairz et al. (2003);
Juffmann et al. (2012) near-field interferometry of the Talbot-Lau type
generates fringes of a fixed period, which are determined by the experimental
geometry. Specifically, the expected interference figure in our KDTLI
configuration is a sine curve whose contrast varies with the phase-shifting
laser power as well as with the molecular beam velocity and polarizability. We
distinguish the genuine quantum interferogram Brezger et al. (2003);
Hornberger et al. (2009) from a classical shadow image by comparing the
expected and experimental interference fringe visibility (contrast) with a
classical model.

The far off-resonance optical polarizability is assumed to be well
approximated by the static value
$\alpha_{opt}\cong\alpha_{stat}\cong\rm{410}\rm{\AA^{3}}\times\rm{4}\pi\epsilon_{0}$
as estimated using Gaussian G09 Frisch et al. (2009) with the 6-31 G basis
set. The absorption cross section of L12 at 532 nm was estimated using the
value of pure tetraphenylporphyrine dissolved in toluene Dixon et al. (2005)
assuming that the perfluoroalkyl chains contribute at least an order of
magnitude less to that value Gotsche et al. (2007). We thus find
$\sigma_{532}\cong 1.7\times\rm{10^{-21}\,m^{2}}$.

In Figure 3(b) we show the expected classical and quantum contrast as a
function of the diffracting laser power. Our experimental contrast is derived
from the recorded signal curves, such as shown in Figure 3(a) by
$V=(S_{max}-S_{min})/(S_{max}+S_{min})$, where $S_{max}$ and $S_{min}$ are the
maxima and minima of the sine curve fitted to the data. The contrast it is
plotted as the black dots in Figure 3(b). These points are in good agreement
with the quantum prediction (blue line). The classical effect, describing the
shadow image by two material gratings and one array of dipole force lenses
Hornberger et al. (2009) is shown as the red line.

The experiment clearly excludes this classical picture. Since the amount of
precious molecular material and thermal degradation processes in the source
limited the measurement time we performed our experiments only in the
parameter regime that maximizes the fringe visibility, to show the clear
difference between the quantum and classical contrast. The experimental fringe
visibility reproduces well the maximally expected quantum contrast.

##  III Conclusion

We have shown that a library approach towards stable and volatile high-mass
molecules can substantially extend the complexity range of molecular de
Broglie coherence experiments to masses in excess of 10 000 amu. Our data
confirm the fully coherent quantum delocalization of single compounds composed
of about 5 000 protons, 5 000 neutrons and 5 000 electrons. The internal
complexity, number of vibrational modes and also the internal energy of each
of these particles is higher than in any other matter-wave experiment so far.

###### Acknowledgements.

This work was supported by the FWF programs Wittgenstein (Z149-N16) and CoQuS
(W1210-2), the European Commission in the project NANOQUESTFIT (304 886), the
Swiss National Science Foundation, the NCCR Nanoscale Science , and the Swiss
Nanoscience Institute (SNI). The authors thank Michel Rickhaus for the artwork
of figure 2.

## References

  * Mooij et al. (1999) J. E. Mooij, T. P. Orlando, L. Levitov, L. Tian, C. H. v. d. Wal, and S. Lloyd, 285, 1036 (1999), ISSN 0036-8075, 1095-9203, URL http://www.sciencemag.org/content/285/5430/1036.
  * Anderson et al. (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, 269, 198 (1995), ISSN 0036-8075, 1095-9203, URL http://www.sciencemag.org/content/269/5221/198.
  * Arndt et al. (1999) M. Arndt, O. Nairz, J. Voss-Andreae, C. Keller, G. van der Zouw, and A. Zeilinger, 401, 680 (1999), ISSN 0028-0836, URL http://dx.doi.org/10.1038/44348.
  * Brezger et al. (2003) B. Brezger, M. Arndt, and A. Zeilinger, 5, 82 (2003), ISSN 1464-4266, URL http://iopscience.iop.org/1464-4266/5/2/362.
  * Gerlich et al. (2007) S. Gerlich, L. Hackermüller, K. Hornberger, A. Stibor, H. Ulbricht, M. Gring, F. Goldfarb, T. Savas, M. Müri, M. Mayor, et al., 3, 711 (2007), ISSN 1745-2473, URL http://dx.doi.org/10.1038/nphys701.
  * Haslinger et al. (2013) P. Haslinger, N. Dörre, P. Geyer, J. Rodewald, S. Nimmrichter, and M. Arndt, 9, 144 (2013), ISSN 1745-2473, URL http://dx.doi.org/10.1038/nphys2542.
  * Deachapunya et al. (2008) S. Deachapunya, P. J. Fagan, A. G. Major, E. Reiger, H. Ritsch, A. Stefanov, H. Ulbricht, and M. Arndt, 46, 307 (2008), ISSN 1434-6060, 1434-6079, URL http://www.springerlink.com/content/v543345622317652/.
  * Berninger et al. (2007) M. Berninger, A. Stefanov, S. Deachapunya, and M. Arndt, 76, 013607 (2007), URL http://link.aps.org/doi/10.1103/PhysRevA.76.013607.
  * Eibenberger et al. (2011) S. Eibenberger, S. Gerlich, M. Arndt, J. Tüxen, and M. Mayor, 13, 043033 (2011), URL http://stacks.iop.org/1367-2630/13/i=4/a=043033.
  * Gring et al. (2010) M. Gring, S. Gerlich, S. Eibenberger, S. Nimmrichter, T. Berrada, M. Arndt, H. Ulbricht, K. Hornberger, M. M., M. Mayor, et al., 81, 031604 (2010).
  * Gerlich et al. (2008) S. Gerlich, M. Gring, H. Ulbricht, K. Hornberger, J. Tüxen, M. Mayor, and M. Arndt, 47, 6195 (2008), ISSN 1521-3773, URL http://onlinelibrary.wiley.com/doi/10.1002/anie.200801942/abstract.
  * Tüxen et al. (2010) J. Tüxen, S. Gerlich, S. Eibenberger, M. Arndt, and M. Mayor, 46, 4145 (2010).
  * Bassi et al. (2013) A. Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht, 85, 471 (2013), URL http://link.aps.org/doi/10.1103/RevModPhys.85.471.
  * Stock et al. (2004) N. L. Stock, D. A. Ellis, L. Deleebeeck, D. C. G. Muir, and S. A. Mabury, 38, 1693 (2004), URL http://dx.doi.org/10.1021/es034773+.
  * Krusic et al. (2005) P. J. Krusic, A. A. Marchione, F. Davidson, M. A. Kaiser, C.-P. C. Kao, R. E. Richardson, M. Botelho, R. L. Waterland, and R. C. Buck, 109, 6232 (2005), URL http://dx.doi.org/10.1021/jp0502961.
  * Meot-Ner et al. (1973) M. Meot-Ner, J. H. Green, and A. D. Adler, 206, 641 (1973), ISSN 1749-6632, URL http://onlinelibrary.wiley.com/doi/10.1111/j.1749-6632.1973.tb43242.x/abstract.
  * Tüxen et al. (2011) J. Tüxen, S. Eibenberger, S. Gerlich, M. Arndt, and M. Mayor, pp. 4823–4833– (2011), ISSN 1099-0690, URL http://onlinelibrary.wiley.com/doi/10.1002/ejoc.201100638/abstract.
  * Gerlich et al. (2011) S. Gerlich, S. Eibenberger, M. Tomandl, S. Nimmrichter, K. Hornberger, P. J. Fagan, J. Tüxen, M. Mayor, and M. Arndt, 2, 263 (2011), URL http://dx.doi.org/10.1038/ncomms1263.
  * Nairz et al. (2002) O. Nairz, M. Arndt, and A. Zeilinger, 65, 65 0321091 (2002).
  * Hornberger et al. (2009) K. Hornberger, S. Gerlich, H. Ulbricht, L. Hackermüller, S. Nimmrichter, I. V Goldt, O. Boltalina, and M. Arndt, 11, 043032 (2009), ISSN 1367-2630, URL http://iopscience.iop.org/1367-2630/11/4/043032.
  * Clauser and Li (1997) J. Clauser and S. Li, in _Atom Interferometry_ , edited by P. R. Berman (Academic Press, 1997), vol. 11, pp. 121–151–.
  * Dirac (1958) P. A. M. Dirac, _The principles of quantum mechanics_ (Clarendon Pr., Oxford, 1958).
  * Hackermüller et al. (2003) L. Hackermüller, K. Hornberger, B. Brezger, A. Zeilinger, and M. Arndt, 77, 781 (2003).
  * Hackermüller et al. (2004) L. Hackermüller, K. Hornberger, B. Brezger, A. Zeilinger, and M. Arndt, 427, 711 (2004).
  * Nairz et al. (2003) O. Nairz, M. Arndt, and A. Zeilinger, 71, 319 (2003), ISSN 0002-9505, URL <GotoISI>://000181806900006.
  * Juffmann et al. (2012) T. Juffmann, A. Milic, M. Müllneritsch, P. Asenbaum, A. Tsukernik, J. Tüxen, M. Mayor, O. Cheshnovsky, and M. Arndt, 7, 297 (2012), ISSN 1748-3387, URL http://www.nature.com/nnano/journal/v7/n5/full/nnano.2012.34.html.
  * Frisch et al. (2009) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. J. Montgomery, T. Vreven, K. N. Kudin, J. C. Burant, et al., _Gaussian 09_ , Wallingford, CT (2009).
  * Dixon et al. (2005) J. M. Dixon, M. Taniguchi, and J. S. Lindsey, 81, 212 (2005), ISSN 1751-1097, URL http://onlinelibrary.wiley.com/doi/10.1111/j.1751-1097.2005.tb01544.x/abstract.
  * Gotsche et al. (2007) N. Gotsche, H. Ulbricht, and M. Arndt, 17, 583 (2007), ISSN 1054-660X, URL <GotoISI>://000245807500043.


# High-order Compact Gas-kinetic Scheme for Two-layer Shallow Water Equations
on Unstructured Mesh

Fengxiang Zhao<EMAIL_ADDRESS>Jianping Gan<EMAIL_ADDRESS>Kun Xu
<EMAIL_ADDRESS>Department of Mathematics, Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, HongKong  Center for Ocean Research in
Hong Kong and Macau (CORE), Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong  Shenzhen Research Institute, Hong Kong
University of Science and Technology, Shenzhen, China

###### Abstract

For the two-layer shallow water equations, a high-order compact gas-kinetic
scheme (GKS) on triangular mesh is proposed. The two-layer shallow water
equations have complex source terms in comparison with the single layer
equations. The main focus of this study is to construct a time-accurate
evolution solution at a cell interface and to design a well-balanced scheme.
The evolution model at a cell interface provides not only the numerical
fluxes, but also the flow variables. The time-dependent flow variables at the
closed cell interfaces can be used to update the cell-averaged gradients for
the discretization of the the source terms inside each control volume in the
development of the well-balanced scheme. Based on the cell-averaged flow
variable and their gradients, high-order initial data reconstruction can be
achieved with compact stencils. The compact high-order GKS has advantages to
simulate the flow evolution in complex domain covered by unstructured mesh.
Many test cases are used to validate the accuracy and robustness of the scheme
for the two-layer shallow water equations.

###### keywords:

Two-layer shallow water equations; Gas-kinetic scheme; High-order compact
reconstruction; Unstructured mesh

##  1 Introduction

The shallow water equations (SWE) are useful in studying both large-scale
ocean circulations and small-scale coastal and channel flows, such as
tsunamis, pollutant transport, tidal waves, and dam break problems. However,
real flows often exhibit stratification, which cannot be captured accurately
by a single-layer SWE. For instance, the injection of freshwater into seawater
creates plumes that are important for the coastal marine environment, with
salinity stratification being a possible feature. In addition, the flow
velocity at coastal area may vary significantly or exhibit stratification
along the depth. To model the stratified water flow, the multi-layer SWE will
be used with the superposition of coupled layers with the force interaction
between them. This paper will focus on the development of high-order compact
scheme for the two-layer SWE (TLSWE), which is the basis for the multi-layer
SWE. In particular, the numerical scheme developed in this study for TLSWE can
be naturally extended to multi-layer SWE with the inclusion of the interaction
between layers as the source term and their dynamic effect in the calculation
of numerical fluxes.

Numerous numerical schemes have been developed for solving SWE with second-
order accuracy [1, 2, 3]. High-order numerical methods have gained popularity
in recent years due to their advantages in accuracy and computational
efficiency [4, 5, 6, 7]. As a result, several high-order numerical schemes
have been proposed for solving SWE [8, 9, 4]. However, there are few works on
numerical methods for the two-layers SWE. Most of them are still based on the
1-D model [10, 11, 12, 13] or 2-D model on structured mesh [14, 15]. Second-
order schemes for 2D TLSWE on unstructured mesh have been developed [16], with
great difficulties due to the loss of hyperbolicity under certain conditions
and the stiff coupling between layers with the product of flow variables and
their derivatives [11, 17, 18, 16].

Unstructured mesh is highly adaptable to complex geometries, making it a
popular choice for numerical studies on real flow simulations [19]. This is
particularly relevant for coastal hydrodynamics simulation, given the
irregular and multiscale nature of coastal boundary geometries. However, to
construct a high-order finite volume scheme on unstructured mesh presents a
challenge due to the use of large stencils in the reconstruction [20]. Most
high-order schemes for single-layer SWE on unstructured mesh are based on the
discontinuous Galerkin (DG) formulation [4, 19], and there are few high-order
finite volume schemes for solving TLSWE. The DG method updates the inner
degrees of freedom (DOFs) from its weak formulation, and is widely used to
solve compressible gas dynamics equations [21, 22] due to its compact spatial
discretization. However, for flows with discontinuities, additional numerical
treatments, such as identifying troubling cells and limiting procedures, must
be designed within the DG framework [23, 24, 25]. In this study, a high-order
compact gas-kinetic scheme (GKS) will be constructed. The finite volume GKS
updates both cell-averaged flow variables and their gradients from the moments
of the time-accurate gas distribution function at a cell interface and compact
initial reconstruction can be obtained. At the same time, the multistage and
multiderivative method will be adopted for achieving high-order temporal
accuracy with less stages [26].

The structure of this paper is as follows. Section 2 introduces the GKS for
TLSWE. Section 3 discusses the high-order compact reconstruction on
unstructured mesh and temporal discretization. In Section 4, the compact GKS
is validated by studying shallow water flow in various cases. Finally, Section
5 is the conclusion.

##  2 Two-layer shallow water equations and gas-kinetic evolution model

This section will present the gas-kinetic evolution model for solving TLSWE.
The corresponding GKS for TLSWE will be constructed based on the extension of
the scheme for the single-layer SWE [27], where the interaction between layers
will be explicitly included in the scheme.

###  2.1 Two-layer shallow water equations

In [18], three equivalent forms of TLSWE are presented. In this study, the
conservative form of TLSWE will be adopted and the interaction between layers
is included in the source term,

$\frac{\partial\textbf{W}}{\partial
t}+\frac{\partial\textbf{F}^{x}(\textbf{W})}{\partial
x}+\frac{\partial\textbf{F}^{y}(\textbf{W})}{\partial
y}=\textbf{S}(\textbf{W}),$ (1)

where

${\textbf{W}}=\left(\begin{array}[]{c}h_{2}\\\ h_{2}U_{2}\\\ h_{2}V_{2}\\\
h_{1}\\\ h_{1}U_{1}\\\ h_{1}V_{1}\\\ \end{array}\right),\\\
{\textbf{F}^{x}}=\left(\begin{array}[]{c}h_{2}U_{2}\\\
h_{2}U_{2}^{2}+\frac{1}{2}Gh^{2}_{2}\\\ h_{2}U_{2}V_{2}\\\ h_{1}{U_{1}}\\\
h_{1}U_{1}^{2}+\frac{1}{2}Gh^{2}_{1}\\\ h_{1}U_{1}V_{1}\\\
\end{array}\right),\\\ {\textbf{F}^{y}}=\left(\begin{array}[]{c}h_{2}V_{2}\\\
h_{2}U_{2}V_{2}\\\ h_{2}V_{2}^{2}+\frac{1}{2}Gh^{2}_{2}\\\ h_{1}V_{1}\\\
h_{1}U_{1}V_{1}\\\ h_{1}V_{1}^{2}+\frac{1}{2}Gh^{2}_{1}\\\
\end{array}\right),$

and

${\textbf{S}}=\left(\begin{array}[]{c}0\\\ -Gh_{2}B_{x}-Gh_{2}h_{1,x}\\\
-Gh_{2}B_{y}-Gh_{2}h_{1,y}\\\ 0\\\ -Gh_{1}B_{x}-\chi Gh_{1}h_{2,x}\\\
-Gh_{1}B_{y}-\chi Gh_{1}h_{2,y}\\\ \end{array}\right).$

Here W are the flow variables, and $\textbf{F}^{x}$ and $\textbf{F}^{y}$ are
the corresponding fluxes in the $x$ and $y$ directions. $B$ is the bottom
topography, $G$ is the gravitational acceleration, and $\chi$ is the density
ratio defined as $\chi=\rho_{2}/\rho_{1}$, where $\rho_{1}$ and $\rho_{2}$ are
the densities of the first and second fluid layer. The flow variables of the
lower and upper layers are denoted as $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$,
respectively. The fluxes of the two layers are
$(\mathbf{F}_{1}^{x},\mathbf{F}_{1}^{y})$ and
$(\mathbf{F}_{2}^{x},\mathbf{F}_{2}^{y})$ with the corresponding source terms
$\mathbf{S}_{1}$ and $\mathbf{S}_{2}$. Fig.1 presents a schematic of the two-
layer shallow water flow.

Figure 1: Schematic of the two-layer shallow water flow.

By adopting the form of TLSWE in Eq. (1), the equations for each layer are
similar as the single-layer SWE except the additional source term related to
the interaction between layers. The source term makes the TLSWE conditionally
hyperbolic [18], which may cause difficulty in the construction of numerical
scheme based on Riemann solver and flux splitting method. In addition, the
source terms related to the interaction between layers are nonlinear, which
have challenges in the discretization for the high-order schemes. In the gas-
kinetic scheme, the dynamics in the TLSWE will be recovered by the time
evolution of gas distribution function, and the effect of source term will be
incorporated into the particle transport process. The numerical fluxes will be
directly evaluated from the time-dependent gas distribution function. Since
the governing equations of each layer in the TLSWE have similar forms, a
general formulation for one of the layers will be presented in the following.

###  2.2 Gas-kinetic evolution model

The GKS is based on the time evolution solution of the gas distribution
function for the flux evaluation [28]. The gas-kinetic BGK model can be
written as [3]

$f_{t}+\textbf{u}\cdot\nabla_{\textbf{x}}f+\nabla\Phi\cdot\nabla_{\textbf{u}}f=\frac{g-f}{\tau},$
(2)

where $f$ is the distribution function $f(\textbf{x},t,\textbf{u})$,
$\textbf{u}=(u,v)$ is the particle velocity, and $g$ is the equilibrium state
approached by $f$. $\tau$ is the relaxation time. $\nabla\Phi$ is the
acceleration of particle due to external force and is related to the source
term in TLSWE, such as the force from bottom topography and the friction. The
equilibrium state $g$ is a Maxwellian distribution function [3],

$\begin{split}g=h\big{(}\frac{\lambda}{\pi}\big{)}e^{-\lambda(\mathbf{u}-\mathbf{U})^{2}},\end{split}$
(3)

where $\lambda$ is defined by $\lambda=1/Gh$. Due to the conservation in
relaxation process from $f$ to $g$, $f$ and $g$ satisfy the compatibility
condition,

$\int\frac{g-f}{\tau}\boldsymbol{\psi}\mathrm{d}\Xi=\textbf{0},$ (4)

where $\boldsymbol{\psi}=(\psi_{1},\psi_{2},\psi_{3})^{T}=(1,u,v)^{T}$ and
$\text{d}\Xi=\text{d}u\text{d}v$.

Based on the moments of the gas distribution function, the flow variables and
their fluxes can be obtained. Due to the similar equations for different
layers, the schemes for layer 1 and layer 2 can be formulated similarly. In
the general scheme, the macroscopic flow variables and the fluxes can be
obtained from the distribution function $f$ as

${\textbf{W}}=\int f\boldsymbol{\psi}\mathrm{d}\Xi,$ (5)

and

${\big{(}\textbf{F}^{x},\textbf{F}^{y}\big{)}^{T}}=\int
f\boldsymbol{\psi}\mathbf{u}\mathrm{d}\Xi.$ (6)

The source term S becomes

${\textbf{S}}=-\int\nabla\Phi\cdot\nabla_{\textbf{u}}f\boldsymbol{\psi}\mathrm{d}\Xi,$
(7)

and $\nabla\Phi$ is determined by

$\nabla\Phi=\mathbf{S}/h,$ (8)

where $\mathbf{S}$ takes $\mathbf{S}_{1}=h_{1}(0,-GB_{x}-G\chi
h_{2,x},-GB_{y}-G\chi h_{2,y})^{\mathrm{T}}$ and
$\mathbf{S}_{2}=h_{2}(0,-GB_{x}-Gh_{1,x},-GB_{y}-Gh_{1,y})^{\mathrm{T}}$ for
layer 1 and layer 2, respectively.

The formal solution of the BGK model in Eq. (2) with external forcing term is

$f(\textbf{x},t,\textbf{u})=\frac{1}{\tau}\int_{0}^{t}g(\textbf{x}^{{}^{\prime}},t^{\prime},\textbf{u}^{{}^{\prime}})e^{-(t-t^{\prime})/\tau}\mathrm{d}t^{\prime}+e^{-t/\tau}f_{0}(\textbf{x}_{0},\textbf{u}_{0}),$
(9)

where x is the numerical quadrature point on the cell interface for flux
evaluation, and x can be set as $(0,0)$ for simplicity in a local coordinate
system with both normal and tangential directions as the x- and y-directions.
The formal solution describes an evolution process for the distribution
function. The trajectory of fluid particle is given by
$\textbf{x}=\textbf{x}^{{}^{\prime}}+\textbf{u}^{{}^{\prime}}(t-t^{{}^{\prime}})+\frac{1}{2}\nabla\Phi(t-t^{{}^{\prime}})^{2}$,
and the velocity of the particle is
$\textbf{u}=\textbf{u}^{{}^{\prime}}+\nabla\Phi(t-t^{{}^{\prime}})$. The
acceleration has a second-order effect ($\sim t^{2}$) on the particle
trajectory, but has the first-order contribution ($\sim t$) to the particle
velocity.

The second-order in time and the ell-balanced explicit evolution solution $f$
is obtained for SWE [27]. In this paper, the same evolution solution of $f$ is
used for the individual layer. The solution of $f$ is

$\begin{split}f(\textbf{x},t,\textbf{u})&=\overline{g}(\mathbf{x},0,\mathbf{u})\big{[}C_{1}+C_{2}\big{(}\overline{\mathbf{a}}^{l}\cdot\mathbf{u}H(u)+\overline{\mathbf{a}}^{r}\cdot\mathbf{u}(1-H(u))\big{)}+C_{3}\overline{A}\big{]}\\\
&+C_{2}\overline{g}(\mathbf{x},0,\mathbf{u})\big{[}-2\alpha_{k,m}\overline{\lambda}\big{(}\nabla\Phi^{l}H(u)+\nabla\Phi^{r}(1-H(u))\big{)}\cdot(\mathbf{u}-\overline{\mathbf{U}})\big{]}\\\
&+C_{4}\big{[}g^{l}(\mathbf{x},0,\mathbf{u})H(u)+g^{r}(\mathbf{x},0,\mathbf{u})(1-H(u))\big{]}\\\
&+C_{5}g^{l}(\mathbf{x},0,\mathbf{u})\big{[}\mathbf{a}^{l}\cdot\mathbf{u}-2\alpha_{k,m}\lambda^{l}\nabla\Phi^{l}\cdot(\mathbf{u}-\mathbf{U}^{l})\big{]}H(u)\\\
&+C_{5}g^{r}(\mathbf{x},0,\mathbf{u})\big{[}\mathbf{a}^{r}\cdot\mathbf{u}-2\alpha_{k,m}\lambda^{r}\nabla\Phi^{r}\cdot(\mathbf{u}-\mathbf{U}^{r})\big{]}(1-H(u)),\end{split}$
(10)

where $\alpha_{k,m}~{}(k=1,2,~{}m=1,2,3)$ are constants for a well-balanced
scheme, $(\alpha_{1,1},\alpha_{1,2},\alpha_{1,3})=(1,3/4,1/4)$ and
$\alpha_{2,m}=1$, $m$ and $k$ is related to taking moment, and the details are
given in the Appendix of [27]. The coefficients $C_{i}~{}(i=1,2,\cdots,5)$ are

$\begin{split}C_{1}&=1-e^{-t/\tau},~{}C_{2}=-\tau(1-e^{-t/\tau})+te^{-t/\tau},~{}C_{3}=-\tau(1-e^{-t/\tau})+t,\\\
C_{4}&=e^{-t/\tau},~{}C_{5}=-te^{-t/\tau}.\end{split}$

The fluxes at the cell interface are evaluated by taking moments of the above
gas distribution function and the total transport of mass and momentum within
a time step can be further integrated in time. More details in the formulation
can be found in [28].

###  2.3 Acceleration force modeling at the interface between two water layers

The interaction between layers is modeled as the acceleration term in the
kinetic equation. The spatial derivatives of the water column height determine
the acceleration, where the values of the height derivatives can be obtained
by the compact reconstruction at the cell interface. However, the possible
discontinuity of the interface can trigger a sudden “pull” or “push” between
water layers. For cases with discontinuities, the spatial derivatives of the
water height from the reconstruction will not be used to calculate the force,
and the “step effect” due to the discontinuity needs to be considered. The
acceleration from a discontinuous interface will be modeled.

Figure 2:  A 1-D schematic for the modeling discontinuous interface between
water layers to evaluate the acceleration term in the BGK model. In the 1-D
discretized space, for the reconstructed distributions of the water interface
inside each cell with the possible discontinuities at the cell interface, the
topography for evaluating the source term between layers is represented by the
connected solid line.

Without loss of generality, for the momentum equation of layer 2 as an
example, the corresponding acceleration is given by

$\displaystyle\nabla\Phi_{2}=-GB_{x}-\chi Gh_{1,x}.$

For the continuous bottom topography $B$ and water height $h_{1}$, the
acceleration can be directly evaluated based on the functions of $B_{x}$ and
the reconstructed $h_{1,x}$. However, when the water height is discontinuous
at a cell interface, such as the reconstructed dash lines in Fig. 2, the
corresponding forcing term between layers will be modeled from a re-
constructed continuous profile at the cell interface. The construction of this
continuous profile will take into account the forcing interaction between
neighboring cells.

Firstly, let’s construct the continuous line at the cell interface. In each
cell, the continuous line connects the respective unique values of the water
height on the cell interfaces $x_{j\pm 1/2}$, which are denoted by
$\widehat{h}_{1}(x_{j\pm 1/2})$ as the black dots in Fig. 2 with the values
given later. The continuous line in the cell is obtained as

$P_{j}^{1}(x)=\frac{1}{2}(\widehat{h}_{1}(x_{j+1/2})+\widehat{h}_{1}(x_{j-1/2}))+\frac{\widehat{h}_{1}(x_{j+1/2})-\widehat{h}_{1}(x_{j-1/2})}{x_{j+1/2}-x_{j-1/2}}(x-x_{j}),$

where $P_{j}^{1}(x)$ is a linear interpolation based on the values at the cell
interfaces of the cell. $\widehat{h}_{1}(x_{j+1/2})$ are modeled based on the
discontinuous left and right states

$\widehat{h}_{1}(x_{j+1/2})=\xi
h^{l}_{1}(x_{j+1/2})+(1-\xi)h^{r}_{1}(x_{j+1/2}),$

where $h^{l}_{1}(x_{j+1/2})$ and $h^{r}_{1}(x_{j+1/2})$ are the reconstructed
values at the cell interface. $\xi$ is a coefficient for the convex
combination, and it is defined as

$\xi=\frac{1}{2}\mathrm{erfc}\big{(}(U^{l}_{1}(x_{j+1/2})+U^{r}_{1}(x_{j+1/2}))/2\big{)},$

where the function $\mathrm{erfc}(\cdots)$ is the complementary error
function, and $U^{l,r}_{1}(x_{j+1/2})$ are the left and right values of the
velocity at the interface. $\mathrm{erfc}(\cdots)$ makes a smooth transition
from $2$ to $0$ when the independent variable covers $(-\infty,+\infty)$ with
a value $\mathrm{erfc}(0)=1$. The above linear distribution in the cell has
dynamically upwind-biased slope.

In the smooth case, the updated derivative of the water height can be used in
the evaluation of acceleration inside each cell. In order cope with both
discontinuous and smooth cases, the final derivative of the water height is
determined by the following nonlinear convex combination method

$\displaystyle\begin{split}\widetilde{h}^{l}_{1,x}(x_{j+1/2})&=w_{j+1/2}h^{l}_{1,x}(x_{j+1/2})+(1-w_{j+1/2})P^{1}_{j,x}(x_{j+1/2}),\\\
\widetilde{h}^{r}_{1,x}(x_{j+1/2})&=w_{j+1/2}h^{r}_{1,x}(x_{j+1/2})+(1-w_{j+1/2})P^{1}_{j+1,x}(x_{j+1/2}),\end{split}$
(11)

where $w_{j+1/2}$ is a nonlinear weighting function to identify the smoothness
of the solution. $w_{j+1/2}$ tends to $1$ in the smooth region and to $0$ in
the discontinuous region. The value of $w_{j+1/2}$ is the same nonlinear
weight as that in the high-order time stepping reconstruction scheme of [29].

In two dimensions, similar modeling of the derivative of the water height can
be done. Different from the one-dimensional one, the modeled continuous line
in Fig. 2 is extended to a 2-D continuous plane. A smooth linear interpolation
in the cell is determined by the following constraints.

$P^{1}(x^{c}_{k},y^{c}_{k})=\widehat{h}_{1}(x^{c}_{k},y^{c}_{k}),~{}k=1,2,3,$

where $(x^{c}_{k},y^{c}_{k})$ is the center of the cell interface,
$\widehat{h}_{1}(x^{c}_{k},y^{c}_{k})$ can be obtained by taking the
arithmetic average of the values $\widehat{h}_{1}(x_{m},y_{m})$, where
$m=1,2$, on the Gaussian quadrature points of the corresponding cell side.

##  3 Compact GKS based on high-order compact reconstruction

In this section, the compact GKS for the TLSWE will be constructed, where the
high-order compact reconstruction to obtain the initial values of flow
distributions is implemented and the two-stage fourth-order (S2O4) temporal
discretization is used. Since the two layers in the shallow water equations
can be numerically treated in the same way, the evolutions for
$\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ will be presented by the discretization
of $\mathbf{W}$ below.

###  3.1 Finite volume discretization

Taking moments $\boldsymbol{\psi}$ on Eq. (2), the flow variables in a cell
$\Omega_{j}$ are updated by

$\frac{\partial\textbf{W}_{j}}{\partial
t}=-\frac{1}{\big{|}\Omega_{j}\big{|}}\int_{\partial\Omega_{j}}\textbf{F}\cdot\textbf{n}\mathrm{d}l+\frac{1}{\big{|}\Omega_{j}\big{|}}\int_{\Omega_{j}}\textbf{S}\mathrm{d}\Omega_{j},$
(12)

where $\textbf{W}_{j}$ is the cell-averaged flow variable,
$\textbf{F}=(\textbf{F}^{x},\textbf{F}^{y})$ is the time-dependent flux at
cell interface, which can be obtained from the moments of the gas distribution
function in Eq. (10). The $\textbf{W}_{j}$ is defined as

$\displaystyle\textbf{W}_{j}$
$\displaystyle\equiv\frac{1}{\big{|}\Omega_{j}\big{|}}\int_{\Omega_{j}}\textbf{W}(\textbf{x})\text{d}\Omega.$
(13)

The line integral of the flux in Eq. (12) can be discretized by a q-point
Gaussian quadrature formula,

$\displaystyle\begin{split}-\frac{1}{\big{|}\Omega_{j}\big{|}}\int_{\partial\Omega_{j}}\textbf{F}\cdot\textbf{n}\mathrm{d}l&=-\frac{1}{\big{|}\Omega_{j}\big{|}}\sum_{l=1}^{l_{0}}\big{(}\big{|}\Gamma_{l}\big{|}\sum_{k=1}^{q}\omega_{k}\textbf{F}(\textbf{x}_{k})\cdot\textbf{n}_{l}\big{)}\\\
&\equiv\mathcal{L}^{F}_{j}(\textbf{W}),\end{split}$ (14)

where $\big{|}\Gamma_{l}\big{|}$ is the side length of the cell, $l_{0}$ is
the total number of cell sides, such as $l_{0}=3$ for a triangular mesh,
$\textbf{n}_{l}$ is the unit outer normal vector, and $q$ and $\omega_{k}$ are
the total number of integration points and weights of the Gaussian integration
formula. In order to evaluate the above numerical flux, the initial data
$\textbf{W}(\textbf{x}_{k})$ is reconstructed using the compact spatial
stencil, which are presented in Section 3.3. The cell-averaged S becomes

$\begin{split}\frac{1}{\big{|}\Omega_{j}\big{|}}\int_{\Omega_{j}}\textbf{S}\mathrm{d}\Omega_{j}\equiv\mathcal{L}^{S}_{j}(\textbf{W}).\end{split}$
(15)

###  3.2 Discretization for source term

The source term in the momentum equations includes two parts, the first one
depends on the bottom topography, and the second one is related to the
variation of the interface between layers and the water height of the up
layer.

The first part of the source term depending on the bottom topography is
defined as

$\begin{split}\frac{1}{\big{|}\Omega_{j}\big{|}}\iint_{\Omega_{j}}\textbf{S}\mathrm{d}\Omega_{j}&=h_{j}(0,-GB_{j,x},-GB_{j,y})^{T}\\\
&\equiv\mathcal{L}^{S_{1}}_{j}(\textbf{W}),\end{split}$ (16)

where $h_{j}$ is the cell average of $h$ in $\Omega_{j}$. High-order spatial
and temporal discretizations of the first part can be implemented directly, as
in the single-layer SWE in [27].

The second part of $\mathcal{L}^{S}_{j}(\textbf{W})$ is related to the
variation of the water height. Taking the source term in the equation of
$h_{1}U_{1}$ as an example, the spatial discretization becomes

$\displaystyle\begin{split}\mathcal{L}^{S_{2}}_{j}(\textbf{W})&\equiv\frac{1}{|\Omega_{j}|}\int_{\Omega_{j}}-\chi
Gh_{1}h_{2,x}\mathrm{d}x\mathrm{d}y\\\ &=-\chi
G\frac{1}{|\Omega_{j}|}\int_{\Omega_{j}}h_{1}\mathrm{d}x\mathrm{d}y\cdot\frac{1}{|\Omega_{j}|}\int_{\Omega_{j}}h_{2,x}\mathrm{d}x\mathrm{d}y+O(\Delta
X^{2})\\\ &=-\chi
G\frac{1}{|\Omega_{j}|}\int_{\Omega_{j}}h_{1}\mathrm{d}x\mathrm{d}y\cdot\frac{1}{|\Omega_{j}|}\int_{\partial\Omega_{j}}h_{2}n_{x}\mathrm{d}\Gamma+O(\Delta
X^{2})\\\ &=-\chi
G\sum_{k=1}^{3}\sum_{l=1}^{2}\widetilde{w}_{k,l}h_{1}(\mathbf{x}_{k,l})\cdot\frac{1}{|\Omega_{j}|}\sum_{k=1}^{3}\big{(}\sum_{l=1}^{2}w_{k,l}h_{2}(\mathbf{x}_{k,l})\big{)}n_{k,x}|\Gamma_{k}|+O(\Delta
X^{2}).\end{split}$ (17)

where $|\Omega_{j}|$, $\big{|}\Gamma_{k}\big{|}$, $\omega_{l}$, $n_{x}$ and
$\mathbf{x}_{k,l}$ have the same definition as those in Eq. (14), $\Delta X$
is the mesh cell size, $\widetilde{w}_{k,l}$ is the weight to obtain the
numerical integration over $\Omega_{j}$ based on $h_{1}(\mathbf{x}_{k,l})$,
and $\widetilde{w}_{k,l}=1/6$. The second-order spatial discretizations is
implemented in Eq. (17). High-order discretization of
$\mathcal{L}^{S_{2}}_{j}(\textbf{W})$ can be achieved by introducing more
numerical integration points. However, considering the balance between
accuracy and efficiency, the simple method given in Eq. (17) is adopted for
the spatial discretization of the second part of the source term in this
paper.

The compact GKS of the TLSWE is a well-balanced scheme. The well-balanced
property is achieved through the balance of the time-accurate flux function at
the cell interface and the spatial discretization of the source terms inside
the control volume. In the previous study [30], the well-balanced GKS for the
single-layer SWE on triangular mesh has been developed, where a corresponding
well-balanced evolution solution of the gas distribution function shown in Eq.
(10) is obtained. For the TLSWE, with the well-balanced initial conditions

$\displaystyle\begin{split}&h_{1}+B=\mathrm{Const},\\\
&h_{2}=\mathrm{Const},\end{split}$

and $(U_{1},V_{1})=(U_{2},V_{2})=(0,0)$, at the quadrature points on the cell
interface, the initial conditions should be $\nabla h_{2}=\mathbf{0}$ and
$\nabla(h_{1}+B)=\mathbf{0}$. With the adoption of water level reconstruction
technique [2], this initial condition can be preserved numerically. As a
result, the compact GKS for the TLSWE can keep such a solution and the scheme
is a well-balanced one. In the following, the solution update in the compact
GKS on the triangular mesh will be presented.

###  3.3 The time evolutions of flow variables and their derivatives

By adopting the S2O4 time stepping method [26, 31], the fully discretized form
of the TLSWE in Eq. (16) over the cell $\Omega_{j}$ in a time step
$[t^{n},t^{n+1}]$ is given by

$\begin{split}\textbf{W}^{n+1/2}_{j}=&\textbf{W}^{n}_{j}+\frac{1}{2}\Delta
t\mathcal{L}_{j}(\textbf{W}^{n})+\frac{1}{8}\Delta
t^{2}\frac{\partial}{\partial t}\mathcal{L}_{j}(\textbf{W}^{n}),\\\
\textbf{W}^{n+1}_{j}=&\textbf{W}^{n}_{j}+\Delta
t\mathcal{L}_{j}(\textbf{W}^{n})+\frac{1}{6}\Delta
t^{2}\frac{\partial}{\partial
t}\mathcal{L}_{j}(\textbf{W}^{n})+\frac{1}{3}\Delta
t^{2}\frac{\partial}{\partial
t}\mathcal{L}_{j}(\textbf{W}^{n+1/2}),\end{split}$ (18)

where $\mathcal{L}_{j}=\mathcal{L}^{F}_{j}+\mathcal{L}^{S}_{j}$ includes the
flux and source term contribution.

Figure 3:  A schematic of flow variables and fluxes on the interface given by
the time-dependent evolution solution in the compact GKS. The dots represent
time accurate flow variables at the interface with possible discontinuities,
but corresponding a single valued flux function.

In the current compact GKS, besides the update of cell-averaged flow variables
in Eq. (18), the cell-averaged derivatives can be updated as well by the
Gauss’s theorem as

$\nabla\mathbf{W}_{j}(t^{n+1})=\frac{1}{|\Omega_{j}|}\int_{\partial\Omega_{j}}{\bf
W}(\mathbf{x},t^{n+1}){\bf n}\mathrm{d}S,$ (19)

with the discretized form

$\displaystyle\nabla\textbf{W}_{j}^{n+1}$
$\displaystyle=\frac{1}{\big{|}\Omega_{j}\big{|}}\sum_{l=1}^{l_{0}}\big{(}|\Gamma_{l}|\mathbf{n}_{l}\sum_{k=1}^{q}\omega_{k}\textbf{W}^{n+1}(\textbf{x}_{k})\big{)},$
(20)

where $|\Omega_{j}|$, $\big{|}\Gamma_{l}\big{|}$, $l_{0}$, $\omega_{k}$ and
$\textbf{n}_{l}$ have the same definition as those in Eq. (14). The flow
variables $\mathbf{W}(\mathbf{x},t^{n+1})$ should be provided at the inner
sides of the cell boundary of the control volume at the time step $t^{n+1}$.
Fig.3 shows the time-accurate flow variables and fluxes on the cell interface
from the evolution solution of the gas distribution function in the compact
GKS. In the discrete scheme, the discontinuous evolution solution
$\mathbf{W}^{l,r}(\mathbf{x},t^{n+1})$ at the cell interface have been
obtained in the GKS for the highly compressible Navier-Stokes solutions [29].
However, in the current study for the shallow water equations, a continuous
evolution solution, namely $\mathbf{W}^{l}=\mathbf{W}^{r}$, for the update the
cell-averaged derivatives within the cell by Eq. (19) seems work very well. In
order to obtain a high-order time-accurate flow variable at the quadrature
point in Eq. (20), the macroscopic flow variable is evolved by two stages

$\begin{split}\textbf{W}^{n+1/2}(\textbf{x})=\textbf{W}^{n}(\textbf{x})+\frac{1}{2}\Delta
t\textbf{W}_{t}^{n}(\textbf{x}),\\\
\textbf{W}^{n+1}(\textbf{x})=\textbf{W}^{n}(\textbf{x})+\Delta
t\textbf{W}_{t}^{n+1/2}(\textbf{x}).\end{split}$ (21)

###  3.4 High-order compact reconstruction

In this section, the high-order compact spatial reconstruction for flow
variables will be presented. Based on the cell averages and their derivatives,
the high-order reconstruction can be obtained compactly with the stencils
involving the closest neighboring cells only, as shown in Fig. 4. The compact
stencil provides consistent domains of dependence between the numerical and
physical ones. The reconstruction with the accuracy from fourth-order to
sixth-order can be obtained on the compact stencils [30]. The fourth-order
reconstruction will be used in this study.

Figure 4:  A schematic of reconstruction stencil of compact GKS. The dotted
circle is the physical domain of dependence of the cell 0, which is covered
compactly by the cells surrounding 0. In each cell, the flow variables and
their derivatives in the x- and y-directions are known.

For the fourth-order reconstruction, $P^{3}(\bm{x})$ polynomial is constructed
as

$\displaystyle\begin{split}P^{3}(\bm{x})=\sum_{k=0}^{9}a_{k}\varphi_{k}(\bm{x}),\end{split}$
(22)

where $a_{k}$ is the degrees of freedom (DOFs) of $P^{3}(\bm{x})$, the total
number of $a_{k}$ is $10$ and the complete polynomial basis with the highest
order of $3$ are included, and $\bm{x}=(x,y)$ is the coordinate. The basis
function $\varphi_{k}(\bm{x})$ can take the zero-averaged basis as

$\displaystyle 1,~{}\delta x-\overline{\delta x}^{(0)},~{}\delta
y-\overline{\delta y}^{(0)},~{}\frac{1}{2}\delta
x^{2}-\overline{\frac{1}{2}\delta x^{2}}^{(0)},~{}\delta x\delta
y-\overline{\delta x\delta y}^{(0)},~{}\frac{1}{2}\delta
y^{2}-\overline{\frac{1}{2}\delta y^{2}}^{(0)},~{}\cdots.$ (23)

To fully determine $P^{3}(\bm{x})$, the DOFs on the cells of the compact
stencil is selected to give the constraints on $P^{3}(\bm{x})$.

$\displaystyle\begin{split}&\big{(}\frac{1}{\big{|}\Omega_{l}\big{|}}\int_{\Omega_{l}}\varphi_{k}(\bm{x})\mathrm{d}x\mathrm{d}y\big{)}a_{k}=Q_{l},\\\
&\big{(}\frac{1}{\big{|}\Omega_{l}\big{|}}\int_{\Omega_{l}}\varphi_{k,x}(\bm{x})\mathrm{d}x\mathrm{d}y\big{)}a_{k}=Q_{l,x},\\\
&\big{(}\frac{1}{\big{|}\Omega_{l}\big{|}}\int_{\Omega_{l}}\varphi_{k,y}(\bm{x})\mathrm{d}x\mathrm{d}y\big{)}a_{k}=Q_{l,y},\end{split}$
(24)

where the same subscript $k$ of $\varphi_{k}$ and $a_{k}$ on the left-hand
side of the equations follow the Einstein summation. $Q_{l}$, $Q_{l,x}$ and
$Q_{l,y}$ are the DOFs in the cells for any component of $\mathbf{W}$.

Due to arbitrary geometrical triangular mesh, the number of equations $M$ in
Eq. (24) should be greater than the number of DOFs $a_{k}$ to avoid an ill-
conditioned system. For the fourth-order reconstruction, the set of DOFs
$S_{0}$ is given by

$\displaystyle
S_{0}=\\{Q_{l_{1}},Q_{l_{2},x},Q_{l_{2},y}\\},~{}l_{1}=0,i,j,k,i_{1},i_{2},\cdots,k_{2},~{}l_{2}=0,i,j,k.$
(25)

Eq. (24) can determine a linear system of $a_{k}$, and it is written as

$a_{0}=Q_{0},$ (26)

and

$\left(\begin{array}[]{cccc}A_{1,1}&A_{1,2}&\cdots&A_{1,9}\\\
A_{2,1}&A_{2,2}&\cdots&A_{2,9}\\\ \vdots&\vdots&\vdots&\vdots\\\
A_{9,1}&A_{9,2}&\cdots&A_{9,9}\\\
\widetilde{A}^{x}_{0,1}&\widetilde{A}^{x}_{0,2}&\cdots&\widetilde{A}^{x}_{0,9}\\\
\widetilde{A}^{y}_{0,1}&\widetilde{A}^{y}_{0,2}&\cdots&\widetilde{A}^{y}_{0,9}\\\
\vdots&\vdots&\vdots&\vdots\\\
\widetilde{A}^{y}_{3,1}&\widetilde{A}^{y}_{3,2}&\cdots&\widetilde{A}^{y}_{3,9}\\\
\end{array}\right)\left(\begin{array}[]{c}a_{1}\\\ a_{2}\\\ \vdots\\\ a_{9}\\\
\end{array}\right)=\left(\begin{array}[]{c}Q_{1}-Q_{0}\\\ Q_{2}-Q_{0}\\\
\vdots\\\ Q_{9}-Q_{0}\\\ Q_{0,x}h\\\ Q_{0,y}h\\\ \vdots\\\ Q_{3,y}h\\\
\end{array}\right),$ (27)

where $A_{l,k}$ and $\widetilde{A}_{l,k}$ are defined as

$\displaystyle\begin{split}A_{l,k}&=\frac{1}{\big{|}\Omega_{l}\big{|}}\int_{\Omega_{l}}\varphi_{k}(\bm{x})~{}\mathrm{d}x\mathrm{d}y,\\\
\widetilde{A}^{i}_{l,k}&=\frac{h}{\big{|}\Omega_{l}\big{|}}\int_{\Omega_{l}}\partial\varphi_{k}(\bm{x})/\partial
r~{}\mathrm{d}x\mathrm{d}y,~{}r=x,y,~{}k=1,2,\cdots,9.\\\ \end{split}$ (28)

The system can be solved by the least square (LS) method. The solution of
$a_{k},~{}(k=1,2,\cdots,9)$ is given by

$\displaystyle\mathbf{a}=\big{[}(\mathbf{A}^{\mathrm{T}}\mathbf{A})^{-1}\mathbf{A}^{\mathrm{T}}\big{]}\mathbf{Q},$
(29)

where $\mathbf{a}$ is the vector of DOFs without $a_{0}$, $\mathbf{A}$ is the
coefficient matrix in Eq. (27), and $\mathbf{Q}$ is the vector of the RHS in
Eq. (27).

To deal with discontinuities in the solution, the nonlinear reconstruction is
needed. The nonlinear compact reconstruction is obtained based on the WENO
method by nonlinearly combining the high-order polynomial $P^{3}$ and several
lower-order polynomials, where the lower-order polynomials are determined
based on the sub-stencils by using the LS method. The nonlinear reconstruction
in the compact GKS has been developed in [30], and the same techniques will be
used here.

##  4 Numerical validations

The compact GKS for the two-layer SWE will be validated by the cases of two-
layer shallow flow in this section. All the computations in this section are
performed on 2-D triangular mesh. The time step used in the computation is
determined by the CFL condition as $\Delta t=CFL\frac{\Delta X}{U_{max}}$,
where $\Delta X$ is the size of the mesh cell,
$U_{max}=\mathrm{max}\\{\sqrt{U_{1}^{2}+V_{1}^{2}}+\sqrt{Gh_{1}},\sqrt{U_{2}^{2}+V_{2}^{2}}+\sqrt{Gh_{2}}\\}$,
and $CFL$ number takes $0.5$. The gravitational acceleration is taken as
$G=9.81$ if not specified.

The collision time $\tau$ in the BGK model for inviscid flow at a cell
interface is defined by

$\displaystyle\tau=\varepsilon\Delta
t+\varepsilon_{num}\displaystyle|\frac{h_{l}^{2}-h_{r}^{2}}{h_{l}^{2}+h_{r}^{2}}|\Delta
t,$

where $\varepsilon=0.05$, $\varepsilon_{num}=5$, and $h^{2}_{l}$ and
$h^{2}_{r}$ are the pressures at the left and right sides of a cell interface.
The reason for including the pressure jump term in the relaxation time is to
enhance the artificial dissipation in case of bore wave.

###  4.1 Accuracy test

The accuracy of the compact GKS with high-order compact reconstruction is
tested. In order to calculate the error in the numerical solution, an initial
condition with analytical evolution solution is used directly

$\displaystyle\begin{split}&h_{1}=0.9+0.02e^{-50((x-1)^{2}+(y-1)^{2})},\\\
&h_{2}=1-h_{1},\end{split}$

with a uniform velocity $(U_{1},V_{1})=(U_{2},V_{2})=(1,1)$. The density ratio
is taken as $\chi=1.0$. The gravitational acceleration is $G=9.81$. The free
boundary condition is taken. The analytical solution of this problem is given
by

$\displaystyle h_{1}(t)=0.9+0.02e^{-50((x-1-t)^{2}+(y-1-t)^{2})},$
$\displaystyle h_{2}(t)=1-h_{1}(t),$
$\displaystyle(U_{1}(t),V_{1}(t))=(U_{2}(t),V_{2}(t))=(1,1).$

The computational domain is taken as $[0,2]\times[0,2]$. The triangular mesh
is used.

The $L^{1}$ errors of $h_{1}$ and $h_{2}$ at $t=0.1$ and the convergence
orders are presented in Table 1. The convergence order of current compact GKS
does not keep the 4th order, which is due to the second-order approximation is
used when discretizing the source term in Eq. (17). Although the optimal 4th-
order convergence is not realized, the advantages of high resolution from the
compact spatial reconstruction will be demonstrated in other complex flow
problems.

$h_{re}$ $Error_{L^{1}}(h_{1})$ $\mathcal{O}_{L^{1}}(h_{1})$
$Error_{L^{1}}(h_{2})$ $\mathcal{O}_{L^{1}}(h_{2})$ 1/8 1.2664e-04 1.2015e-04
1/16 1.9987e-05 2.66 1.8203e-05 2.72 1/32 1.0744e-06 4.22 9.4509e-07 4.27 1/64
1.4263e-07 2.91 1.2522e-07 2.92

Table 1:  Accuracy test: errors and convergence orders of compact GKS with
high-order reconstruction.

###  4.2 Well-balanced property

The well-balanced property of the compact GKS on unstructured mesh is
validated in the following. The initial condition is a two-dimensional steady
state solution with non-flat bottom topography. The bottom topography is

$\displaystyle B(x,y)=0.5e^{-50[(x-1)^{2}+(y-1)^{2}]}.$

The steady state is

$\displaystyle\begin{split}&h_{1}=0.8-B(x,y),\\\ &h_{2}=0.2,\end{split}$

and all the velocities are $0$. The density ratio and the gravitational
acceleration are taken as $\chi=1.0$ and $G=9.81$, respectively. The
computational domain is $[0,2]\times[0,2]$. The triangular mesh with cell size
$\Delta X=0.05$ is used. The wall boundary condition is imposed on all the
boundaries.

The discretized bottom topography is shown in Fig. 5. The errors history of
flow variables is plotted in Fig. 6. The error remains at the same level at
different computational time. At very long computation times, the errors of
water surface level and momentum are less than $1.0\times 10^{-8}$. The
current compact GKS is able to maintain an initial balanced steady state
solution.

Figure 5:  Well-balanced property study: the left figure shows the bottom
profile and the unstructured mesh with cell size $\Delta X=0.05$, and the
right figure presents the water level contours at $t=100$.

Figure 6:  Well-balanced property study: the left figure shows the water level
distributions along the horizontal centerline in the computational domain at
$t=100$ and the right figures presents the time evolution of the numerical
error.

Figure 7:  Riemann problem I: the left figure shows the unstructured mesh with
cell size $\Delta X=1/100$, and the right figure is the 3-D water surface
distributions at $t=0.1$.

Figure 8:  Riemann problem I: the left figure is the 1-D distributions of
water level and the right figure is about the discharge along the horizontal
centerline at $t=0.1$.

Figure 9:  Riemann problem I: comparison of the water levels obtained by the
compact GKS and from the reference solution.

###  4.3 Riemann problems of TLSWE

In this section, the Riemann problems with a discontinuity at the interface
between two fluid layers are studied to validate the compact GKS for TLSWE.
Due to unequal densities of the two layers of fluid, the discontinuity at the
interface will evolve and propagate.

The first test was introduced to verify the stability of the numerical schemes
for unsteady two-layer exchange flows [10]. It can also be used to evaluate
the accuracy of different numerical schemes in computing unsteady solutions
over a flat bottom. The initial water level is set as

$(h_{1},h_{2})=\begin{cases}(0.5,0.5),~{}~{}~{}~{}0\leq x<0.3,\\\
(0.55,0.45),0.3\leq x\leq 1,\end{cases}$

and the uniform velocity $(U_{1},V_{1})=(U_{2},V_{2})=(2.5,0)$ is given in the
whole domain. In the computation, the 2-D computational domain is taken as
$[0,1]\times[0,0.5]$, and the triangular mesh is used. The computational time
is $t=0.1$. The density ratio is $\chi=0.98$. The gravitational acceleration
is taken as $G=10$ in this case.

The coarse mesh with $\Delta X=1/100$ used in the computation and the 3-D
water surface obtained by the compact GKS are shown in Fig. 7. The solution of
the evolved free surface has a square-wave structure with small variation. The
current compact GKS captures this solution with no obvious numerical
oscillations. In Fig. 8 and Fig. 9, the water levels along the horizontal
centerline of the computational domain is plotted, where the results on a
finer mesh with $\Delta X=1/400$ are also given to verify the mesh convergence
solution from the current compact GKS. To quantitatively verify the
correctness of the results obtained by the current scheme, the reference
solution obtained by the 1-D model with a cell size of $\Delta X=1/10000$ in
[14] is also plotted. The compact GKS gives consistent solutions on both
coarse and fine meshes. The resolution of the local solution structure on the
fine mesh by the compact GKS is comparable to the reference solution.

Figure 10:  Riemann problem II: the left figure shows the 3-D water level
contours of $h_{1}$ and the right figure gives the 1-D distribution of $h_{1}$
along horizontal centerline at $t=1$.

Figure 11:  Riemann problem II: the left figure is the 1-D distributions of
water level and the right figure is about the discharge along the horizontal
centerline at $t=1$.

The second case is the Riemann problem with a large discontinuity at the
interface between the two layers [14]. The initial value of water levels is
given by

$(h_{1},h_{2})=\begin{cases}(0.2,1.8),~{}0\leq x<5,\\\ (1.8,0.2),~{}5\leq
x<10.\end{cases}$

The initial velocity is $0$. The water density ratio is $\chi=0.98$. The
gravitational acceleration is taken as $G=9.81$. The computational domain is
set as $[0,10]\times[0,1]$. The triangular mesh with a cell size of $\Delta
X=1/40$ is used in the computation.

The evolved results at $t=1.0$ obtained by the compact GKS is presented in
Fig. 10 and Fig. 11. In Fig. 10 the result of $h_{1}$ on the 2-D triangular
mesh is compared with the reference solution presented in [14]. Good agreement
has been obtained. The water levels of the first layer together with the water
surface and discharge are plotted in Fig. 11.

Figure 12:  Dam-break flow at $\chi=1$: the 3-D contours of $h_{1}+h_{2}$
(left) and $h_{1}$ (right) at $t=0.08$. The cell sizes of the coarse and fine
meshes are $\Delta X_{CM}=1/100$ and $\Delta X_{FM}=1/400$.

Figure 13:  Dam-break flow at $\chi=1$: the water levels (left)and discharge
distributions (tight) along the horizontal centerline at $t=0.08$. The cell
sizes of the coarse and fine meshes are $\Delta X_{CM}=1/100$ and $\Delta
X_{FM}=1/400$.

###  4.4 Dam-break problems at different density ratios

The two-layer dam-break flows are used to validate the compact GKS. The
initial state is given as

$(h_{1},h_{2})=\begin{cases}(0.357,1),0\leq x<0.5,\\\ (0.357,0),0.5\leq x\leq
1.\end{cases}$

The velocity is set as $(U_{1},V_{1})=(U_{2},V_{2})=(0,0)$ in the whole
domain, and the computational domain is $[0,1]\times[0,0.5]$. The
gravitational acceleration is $G=9.81$. Dam-break flows at two density ratios
are studied. In the computation, a coarse triangular mesh with $\Delta
X=1/100$ and a fine triangular mesh with $\Delta X=1/400$ are used.

The first case is the dam-break flow at same density of the two layers, i.e.,
the density ratio with $\chi=1$. The 3-D water level distributions of
$h_{1}+h_{2}$ and $h_{1}$ at $t=0.08$ obtained by the compact GKS on the
coarse mesh are shown in Fig. 12. The water levels and discharge distributions
along the horizontal centerline are given in Fig. 13. The results on the
coarse mesh are consistent with those on the fine mesh, and the water levels
obtained by the current compact scheme are consistent with those in [18].

Figure 14:  Dam-break flow of a light fluid over a dense fluid with
$\chi=0.2$: the 3-D height distributions of $h_{1}+h_{2}$ (left) and $h_{1}$
(right) at $t=0.08$. The cell sizes of the coarse and fine meshes are $\Delta
X_{CM}=1/100$ and $\Delta X_{FM}=1/400$.

Figure 15:  Dam-break flow of a light fluid over a dense fluid with
$\chi=0.2$: the height (left) and discharge (right) along the horizontal
centerline at $t=0.08$. The cell sizes of the coarse and fine meshes are
$\Delta X_{CM}=1/100$ and $\Delta X_{FM}=1/400$.

The second case is the dam-break flow of a light fluid over a dense one. The
density ratio is $\chi=0.2$. The 3-D water level distributions of
$h_{1}+h_{2}$ and $h_{1}$ at $t=0.08$ obtained by the compact GKS are shown in
Fig.14. The 1-D water levels and discharge distributions along the horizontal
centerline are given in Fig.15. Due to the complexity of the solution, the
fine mesh result has a better spatial resolution and gives the solution close
to the reference ones in [18].

Figure 16:  Channel flow with non-flat bottom: the 2-D contours of the water
levels (up) and the water level distributions along the horizontal centerline
(down) at $t=0.1$ (left) and $t=1.0$ (right).

###  4.5 Channel flow with non-flat bottom

This case is about the two-layer flow through a channel with non-flat bottom
topography. The bottom topography is defined by

$B(x,y)=0.5e^{-100(x-0.5)^{2}}.$

The initial condition is given as

$\begin{split}h_{1}&=0.8-B(x,y),~{}~{}h_{2}=0.4,\\\
U_{1}&=-0.2,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}U_{2}=0.15.\end{split}$

The channel covers a domain $[0,1]\times[0,0.25]$. The reflecting boundary
condition is applied at the channel walls. The free boundary condition is used
on the left and right boundaries. The triangular mesh with a cell size of
$\Delta X=1/200$ is used in the computation.

Fig. 16 shows the results of water levels at $t=0.1$ and $t=1.0$,
respectively. Due to the non-flat bottom topography, the interface between two
layer fluids evolves from an initial smooth interface to a discontinuous one.
The reference solution comes from solving the 1-D TLSWE on a uniform mesh with
$1000$ cells in [13]. At the early time, a smooth interface evolves, such as
the left figures in Fig. 16, and the solution has good agreement with the
reference solution. At a later time, a discontinuous interface emerges, such
as the right figures in Fig. 16, and the position of the discontinuity
obtained by the compact GKS has a good match with the reference solution.

###  4.6 2-D interface propagation

The 2-D circular interface propagation Riemann problem is studied. The initial
condition of the test case is given by

$(h_{1},h_{2})=\begin{cases}(1.8,0.2),~{}~{}(x-5)^{2}+(y-5)^{2}<4.0,\\\
(0.2,1.8),~{}~{}\mathrm{otherwise}.\end{cases}$

The initial velocity is $(U_{1},V_{1})=(U_{2},V_{2})=(0,0)$ in the
computational domain $[0,10]\times[0,10]$. The gravitational acceleration is
$G=9.81$. The density ratio between layers is $\chi=0.98$. The free boundary
condition is adopted on all boundaries. The triangular mesh with a cell size
of $\Delta X=1/10$ is used in the computation.

The 3-D water level distributions of $h_{1}$ and its distributions along the
horizontal centerline at $t=0$, $t=2.0$ and $t=4.0$ are presented in Fig. 17
and Fig. 17, respectively. The results show the circular propagation of the
water column.

Figure 17:  2D interface propagation: the 3-D water level distributions of
$h_{1}$ at $t=0$, $t=2.0$ and $t=4.0$. The cell size of the triangular mesh is
$\Delta X=1/10$.

Figure 18:  2D interface propagation: the water level distribution of $h_{1}$
along the horizontal centerline at $t=0$, $t=2.0$ and $t=4.0$. The cell size
of the triangular mesh is $\Delta X=1/10$.

Figure 19:  2-D dam-break in an irregular domain: the left figure is the
computational domain and mesh, and the right figure is the enlarged view of
the mesh around the dam breach. The mesh size far away from the dam is $\Delta
X=2.5$, and the mesh size is refined by $3.3$ times in the region close to dam
breach.

Figure 20:  2-D dam-break in an irregular domain with a wet bed: the 3-D
contours of water levels and the distributions along the horizontal centerline
of the breach. The computational time is $t=7.2$.

Figure 21:  2-D dam-break in an irregular domain with a dry bed: the 3-D
contours of water levels and the distributions along the horizontal centerline
of the breach. The computational time is $t=7.2$.

###  4.7 2-D dam-break in an irregular domain

The 2-D dam-break problem in [32, 9] is used in the current study to validate
the compact GKS. Fig. 19 shows the computational domain and the mesh. The
length of the dam breach is $75$ and it starts at $y=95$. The dam itself has a
width of $10$ and its left side is located at $x=95$. At $t=0$ the stationary
water surface has a discontinuity with $h_{l}=10$ and $h_{r}=\epsilon$ across
the breach, and two values of $\epsilon=5$ and $\epsilon=1\times 10^{-3}$ are
used to simulate the wet and dry bed cases, respectively. For the wet case,
the individual water levels of layer 1 and layer 2 are set as

$(h_{1},h_{2})=\begin{cases}(9,1),~{}0\leq x<95,\\\ (5,0),~{}95\leq
x.\end{cases}$

For the dry case, the individual water levels of layer1 and layer 2 are set as

$(h_{1},h_{2})=\begin{cases}(9,1),~{}0\leq x<95,\\\
(\epsilon,\epsilon),~{}95\leq x.\end{cases}$

The boundary condition on the far right is the free boundary, and the other
boundary conditions are the non-penetration slip wall boundaries. The mesh
size far from the breach is $h_{mesh}=2.5$, and is locally refined by $3.3$
times around the dam breach.

The 3-D water surface heights at $t=7.2$ are shown in Fig.20 and Fig.21. The
discontinuous bore waves are captured without spurious oscillation. It clearly
shows that the wave propagating speed is higher in the dry bed case.

##  5 Conclusion

In this study, we have developed a compact high-order Gas-Kinetic Scheme (GKS)
on a triangular mesh to solve the Two-Layer Shallow Water Equations (TLSWE).
The compact scheme is highly accurate and robust in capturing discontinuous
solutions.

The gas evolution model at the cell interface in the kinetic scheme explicitly
captures the dynamics from the particle free transport, collisions, and
acceleration from the external forcing term on the particle trajectory. The
time-accurate evolution solution provides not only the flow variable update
inside each cell, but also the gradients of the flow variables. As a result,
based on the updated flow variables and their gradients, compact stencil can
be used in the reconstruction and the design of the compact scheme.

The compact GKS has several key features in solving TLSWE. The high-order
compact reconstruction on a triangular mesh is naturally obtained. The
existence of the time derivative of the flux function uses less stages to get
a scheme with high-order accuracy in time, such as the two-stages for the
fourth-order time accuracy. This compact GKS provides accurate numerical
solutions for the TLSWE and is ready for its engineering application in the
coastal area ocean flow.

## Acknowledgments

The current research is supported by CORE as a joint research centre for ocean
research between QNLM and HKUST through the project QNLM20SC01-A and
QNLM20SC01-E, the National Natural Science Foundation of China (No. 12172316),
and Hong Kong research grant council 16208021 and 16301222.

## References

## References

  * [1] R. J. LeVeque, “Balancing source terms and flux gradients in high-resolution godunov methods: the quasi-steady wave-propagation algorithm,” Journal of computational physics, vol. 146, no. 1, pp. 346–365, 1998.
  * [2] J. G. Zhou, D. M. Causon, C. G. Mingham, and D. M. Ingram, “The surface gradient method for the treatment of source terms in the shallow-water equations,” Journal of Computational physics, vol. 168, no. 1, pp. 1–25, 2001.
  * [3] K. Xu, “A well-balanced gas-kinetic scheme for the shallow-water equations with source terms,” Journal of Computational Physics, vol. 178, no. 2, pp. 533–562, 2002.
  * [4] D. Wirasaet, E. Kubatko, C. Michoski, S. Tanaka, J. Westerink, and C. Dawson, “Discontinuous galerkin methods with nodal and hybrid modal/nodal triangular, quadrilateral, and polygonal elements for nonlinear shallow water flow,” Computer methods in applied mechanics and engineering, vol. 270, pp. 113–149, 2014.
  * [5] Z. J. Wang, “High-order methods for the euler and navier–stokes equations on unstructured grids,” Progress in Aerospace Sciences, vol. 43, no. 1-3, pp. 1–41, 2007.
  * [6] M. Alhawwary and Z. Wang, “On the mesh resolution of industrial les based on the dns of flow over the t106c turbine,” Advances in Aerodynamics, vol. 1, no. 1, pp. 1–18, 2019.
  * [7] Z. Wang and E. Jourdan, “Benchmark for scale-resolving simulation with curved walls: the taylor couette flow,” Advances in Aerodynamics, vol. 3, no. 1, pp. 1–18, 2021.
  * [8] Y. Xing and C.-W. Shu, “High order well-balanced finite volume weno schemes and discontinuous galerkin methods for a class of hyperbolic systems with source terms,” Journal of Computational Physics, vol. 214, no. 2, pp. 567–598, 2006.
  * [9] M. Ricchiuto, R. Abgrall, and H. Deconinck, “Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes,” Journal of Computational Physics, vol. 222, no. 1, pp. 287–331, 2007.
  * [10] M. Castro, J. Macías, and C. Parés, “A q-scheme for a class of systems of coupled conservation laws with source term. application to a two-layer 1-d shallow water system,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 35, no. 1, pp. 107–127, 2001.
  * [11] R. Abgrall and S. Karni, “Two-layer shallow water system: a relaxation approach,” SIAM Journal on Scientific Computing, vol. 31, no. 3, pp. 1603–1627, 2009.
  * [12] F. Bouchut and V. Zeitlin, “A robust well-balanced scheme for multi-layer shallow water equations,” Discrete and Continuous Dynamical Systems-Series B, vol. 13, no. 4, pp. 739–758, 2010.
  * [13] N. Krvavica, M. Tuhtan, and G. Jelenić, “Analytical implementation of roe solver for two-layer shallow water equations with accurate treatment for loss of hyperbolicity,” Advances in water resources, vol. 122, pp. 187–205, 2018.
  * [14] A. Kurganov and G. Petrova, “Central-upwind schemes for two-layer shallow water equations,” SIAM Journal on Scientific Computing, vol. 31, no. 3, pp. 1742–1773, 2009.
  * [15] M. Dudzinski and M. Lukáčová-Medvid’ová, “Well-balanced bicharacteristic-based scheme for multilayer shallow water flows including wet/dry fronts,” Journal of Computational Physics, vol. 235, pp. 82–113, 2013.
  * [16] X. Liu, “A new well-balanced finite-volume scheme on unstructured triangular grids for two-dimensional two-layer shallow water flows with wet-dry fronts,” Journal of computational physics, vol. 438, p. 110380, 2021.
  * [17] M. J. Castro-Díaz, E. D. Fernández-Nieto, J. M. González-Vida, and C. Parés-Madroñal, “Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system,” Journal of Scientific Computing, vol. 48, pp. 16–40, 2011.
  * [18] B. Spinewine, V. Guinot, S. Soares-Frazão, and Y. Zech, “Solution properties and approximate riemann solvers for two-layer shallow flow models,” Computers & Fluids, vol. 44, no. 1, pp. 202–220, 2011.
  * [19] S. Brus, D. Wirasaet, E. Kubatko, J. Westerink, and C. Dawson, “High-order discontinuous galerkin methods for coastal hydrodynamics applications,” Computer Methods in Applied Mechanics and Engineering, vol. 355, pp. 860–899, 2019.
  * [20] M. Dumbser, D. S. Balsara, E. F. Toro, and C.-D. Munz, “A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes,” Journal of Computational Physics, vol. 227, no. 18, pp. 8209–8253, 2008.
  * [21] B. Cockburn and C.-W. Shu, “The Runge–Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems,” Journal of Computational Physics, vol. 141, no. 2, pp. 199–224, 1998.
  * [22] Z. J. Wang and H. Gao, “A unifying lifting collocation penalty formulation including the discontinuous galerkin, spectral volume/difference methods for conservation laws on mixed grids,” Journal of Computational Physics, vol. 228, no. 21, pp. 8161–8186, 2009.
  * [23] J. Qiu and C.-W. Shu, “Hermite WENO schemes and their application as limiters for Runge–Kutta discontinuous Galerkin method II: Two dimensional case,” Computers & Fluids, vol. 34, no. 6, pp. 642–663, 2005.
  * [24] J. Zhu and J. Qiu, “Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method, III: Unstructured meshes,” Journal of Scientific Computing, vol. 39, no. 2, pp. 293–321, 2009.
  * [25] C.-W. Shu, “High order WENO and DG methods for time-dependent convection-dominated PDEs: A brief survey of several recent developments,” Journal of Computational Physics, vol. 316, pp. 598–613, 2016.
  * [26] J. Li, “Two-stage fourth order: temporal-spatial coupling in computational fluid dynamics (cfd),” Advances in Aerodynamics, vol. 1, pp. 1–36, 2019\.
  * [27] F. Zhao, J. Gan, and K. Xu, “The study of shallow water flow with bottom topography by high-order compact gas-kinetic scheme on unstructured mesh,” Physics of Fluids, vol. 33, no. 8, p. 083613, 2021.
  * [28] K. Xu, “A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method,” Journal of Computational Physics, vol. 171, no. 1, pp. 289–335, 2001.
  * [29] F. Zhao, X. Ji, W. Shyy, and K. Xu, “Direct modeling for computational fluid dynamics and the construction of high-order compact scheme for compressible flow simulations,” Journal of Computational Physics, p. 111921, 2023.
  * [30] F. Zhao, X. Ji, W. Shyy, and K. Xu, “A compact high-order gas-kinetic scheme on unstructured mesh for acoustic and shock wave computations,” Journal of Computational Physics, vol. 449, p. 110812, 2022.
  * [31] L. Pan, K. Xu, Q. Li, and J. Li, “An efficient and accurate two-stage fourth-order gas-kinetic scheme for the Euler and Navier–Stokes equations,” Journal of Computational Physics, vol. 326, pp. 197–221, 2016\.
  * [32] C. Mingham and D. Causon, “High-resolution finite-volume method for shallow water flows,” Journal of Hydraulic Engineering, vol. 124, no. 6, pp. 605–614, 1998.