# Event topology and constituent-quark scaling of elliptic flow in heavy-ion
collisions at the Large Hadron Collider using a multiphase transport model

Neelkamal Mallick1 Sushanta Tripathy2 Raghunath Sahoo1,3111Corresponding
author<EMAIL_ADDRESS>1Department of Physics, Indian Institute of
Technology Indore, Simrol, Indore 453552, India  2INFN - sezione di Bologna,
via Irnerio 46, 40126 Bologna BO, Italy  3CERN, CH 1211, Geneva 23,
Switzerland

###### Abstract

Transverse spherocity is an event shape observable, which has got a unique
capability to separate the events based on their geometrical shapes. In this
work, we use transverse spherocity to study the identified light flavor
production in heavy-ion collisions using A Multi-Phase Transport (AMPT) model.
We obtain the elliptic flow coefficients for pions, kaons and protons in Pb-Pb
collisions at $\sqrt{s_{\rm{NN}}}=5.02$ TeV as a function of transverse
spherocity and collision centrality. Also, we study the number of constituent-
quark (NCQ) scaling of elliptic flow which interprets the dominance of the
quark degrees of freedom at early stages of the collision. We observe a clear
dependence of the elliptic flow for identified particles on transverse
spherocity. It is found that the NCQ-scaling is strongly violated in events
with low transverse spherocity compared to transverse spherocity-integrated
events.

###### pacs:

##  I Introduction

The ultra-relativistic heavy-ion collisions at particle accelerators like the
Large Hadron Collider (LHC) at CERN, Switzerland and Relativistic Heavy-Ion
Collider (RHIC) at BNL, USA, provide signatures of possible formation of a
deconfined state of quarks and gluons, known as Quark Gluon Plasma (QGP).
Although, one does not observe any direct evidence of possible QGP formation
due to its very short lifetime, instead several indirect signatures such as
direct photon measurements, jet-quenching, elliptic flow ($v_{2}$),
strangeness enhancement, quarkonia suppression etc. suggest that formation of
QGP is highly probable in such collisions. Due to the unprecedented energy
scale of a few TeV, some of the QGP signatures ALICE:2017jyt ;
Khachatryan:2016txc  are also observed in smallest collision systems (i.e.
proton-proton collisions) at the LHC. To understand the underlying dynamics of
small collision systems, an event shape observable, transverse spherocity, has
been introduced recently Cuautle:2014yda ; Cuautle:2015kra ; Salam:2009jx ;
Bencedi:2018ctm ; Banfi:2010xy  at the LHC energies. From these studies, it
was observed that transverse spherocity has a very unique capability to
separate the events based on their geometrical shapes, i.e. isotropic and
jetty events Tripathy:2019blo ; Khuntia:2018qox . After its successful
implementation in small collision systems, transverse spherocity can be used
as a tool for heavy-ion collisions to differentiate events based on their
event topology, even the environment is much more complex. It might reveal new
and unique results from heavy-ion collisions where the production of a QGP
medium is already established. In one of our recent publications
Mallick:2021rsd , we have explicitly used transverse spherocity in heavy-ion
collisions for the first time to study the azimuthal anisotropy of
unidentified charged particles as a function of transverse spherocity in A
Multi-Phase Transport (AMPT) model. We observe a strong anti-correlation of
transverse spherocity with the ellipticity of the events in heavy-ion
collisions. It was found that high transverse spherocity events have nearly
zero $v_{2}$ while low transverse spherocity events contribute significantly
to the $v_{2}$.

$v_{2}$ is generally caused by the initial spatial anisotropy in the system
produced in a non-central collision. It plays a significant role to understand
the collective motion and bulk property of the QGP. $v_{2}$ is defined as the
second-order Fourier component of the particle azimuthal distribution, which
provides information of the transport properties of created medium in heavy-
ion collisions v2 . The comparison of $v_{2}$ measurements to hydrodynamic
calculations reveal that the elliptic flow is built up mainly during the early
partonic stages of QGP and it is governed by the evolution of QGP
Abelev:2014pua ; Arsene:2004fa ; Adcox:2004mh ; Back:2004je ; Adams:2005dq .
However, the hadronic rescattering in the hadronic phase could also contribute
to the evolution of elliptic flow Hirano:2005xf . Thus, to understand the
interplay of partonic versus hadronic phase in the evolution of elliptic flow,
it is important to study the elliptic flow of different identified particles.

Figure 1: (Color Online) Transverse momentum space correlation ($p_{\rm y}$
vs. $p_{\rm x}$) for $\pi^{+}+\pi^{-}$ in different spherocity classes for
(0-10)% (top) and (40-50)% (bottom) central Pb-Pb collisions at
$\sqrt{s_{\rm{NN}}}=5.02$ TeV in AMPT model.

In addition, RHIC and LHC experimental results found that in the intermediate
transverse momentum range, the $v_{2}$ of baryons is larger than that of
mesons Abelev:2014pua . Ref. Voloshin:2002wa  suggests that this phenomenon
could be explained by the development of hydrodynamic flow in the partonic
phase and then quarks coalesce into hadrons during the hadronization. This
mechanism leads to number of constituent quarks (NCQ) scaling, a hierarchy
observed in the values of $v_{2}$ as a function of transverse momentum, where
flow might no longer be dominant and this may compete with the contribution
from jet fragmentation. This scaling property has been studied extensively in
experiments Adams:2003am ; Adler:2003kt ; Abelev:2014pua ; v2_ALICE2  and in
theoretical studies Molnar:2003ff ; Greco:2003mm ; Fries:2003kq . Some of
these results suggest that the NCQ scaling of $v_{2}$ at the LHC energies is
violated Abelev:2014pua ; v2_ALICE2  for intermediate or high transverse
momenta. It would be very interesting to study these properties of $v_{2}$ in
different event shapes. In this work, we obtain the transverse momentum
spectra, two-particle correlations and elliptic flow of identified particles
as a function of transverse spherocity in Pb-Pb collisions at
$\sqrt{s_{\rm{NN}}}=5.02$ TeV using AMPT model. We believe that the results
would motivate experimentalists to pursue such a novel method in experiments
at the LHC. In fact, in a recent study it has been shown that how machine
learning tools can be implemented in experiments to obtain transverse
spherocity Mallick:2021wop .

The paper is organised as follows. We begin with a brief introduction and
motivation for the study in Section I. In Section II, the detailed analysis
methodology along with brief description about AMPT are given. Section III
discusses about the results and finally they are summarized in Section IV.

##  II Event Generation and Analysis Methodology

In this section, we begin with a brief introduction on AMPT model and then we
proceed with defining transverse spherocity and elliptic flow.

###  II.1 A Multi-Phase Transport (AMPT) Model

AMPT model contains four components namely, initialization, parton transport,
hadronization mechanism and hadron transport AMPT2 . The initialization of the
collisions is performed using HIJING ampthijing , where the differential
cross-section of the produced mini-jets in $pp$ collisions is calculated.
Then, the produced partons calculated in $pp$ collisions is converted into
$p$-A and heavy-ion collisions by parametrized shadowing function and nuclear
overlap function with inbuilt Glauber model. The initial low-momentum partons
are produced from parametrized coloured string fragmentation mechanisms and
they are separated from high momentum partons by a momentum cut-off. The
produced partons are propagated into parton transport part via Zhang’s Parton
Cascade (ZPC) model amptzpc . In the String Melting version of AMPT, melting
of colored strings into low momentum partons takes place at the start of the
ZPC. It is calculated using Lund FRITIOF model of HIJING. The resulting
partons undergo multiple scatterings which take place when any two partons are
within a distance of minimum approach. In AMPT-SM, the transported partons are
finally hadronized using spatial coalescence mechanism Lin:2001zk ; He:2017tla
. The produced hadrons further undergo final evolution in a relativistic
transport mechanism amptart1 ; amptart2  via meson-meson, meson-baryon and
baryon-baryon interactions. There is also a default version of AMPT, where
instead of coalescing the partons, fragmentation mechanism using Lund
fragmentation parameters $a$ and $b$ are used for hadronizing the transported
partons. However, the particle flow and spectra at the mid-$p_{\rm T}$ regions
are well explained by quark coalescence mechanism for hadronization
ampthadron1 ; ampthadron2 ; ampthadron3 . Thus, we have used string melting
mode of AMPT (AMPT version 2.26t9b) for all of our calculations in the current
work. The AMPT settings in the current work, are the same as reported in Refs.
Tripathy:2018bib ; Mallick:2021rsd . One should note that, to exactly match
with the experimental data, one can vary the tunes of the AMPT model, which is
currently out of the scope of this manuscript. For the input of impact
parameter values for different centralities in Pb-Pb collisions, we have used
Ref. Loizides:2017ack .

###  II.2 Transverse Spherocity

Transverse spherocity is an event property which is defined for a unit vector
$\hat{n}(n_{T},0)$ that minimizes the ratio Cuautle:2014yda ; Cuautle:2015kra
:

$\displaystyle
S_{0}=\frac{\pi^{2}}{4}\bigg{(}\frac{\Sigma_{i}~{}|\vec{p}_{T_{i}}\times\hat{n}|}{\Sigma_{i}~{}p_{T_{i}}}\bigg{)}^{2}.$
(1)

By restricting it to the transverse plane, transverse spherocity becomes
infrared and collinear safe Salam:2009jx . By construction, the extreme limits
of transverse spherocity are related to specific configurations of events in
transverse plane i.e. jetty and isotropic. The value of transverse spherocity
ranges from 0 to 1, which is ensured by multiplying the normalization constant
$\pi{{}^{2}}/4$ in Eq. 1. Here onwards, for the sake of simplicity the
transverse spherocity is referred as spherocity. To disentangle the
low-$S_{0}$ and high-$S_{0}$ events from the spherocity-integrated events, we
have applied spherocity cuts on the generated events. The spherocity
distributions are selected in the pseudo-rapidity range of $|\eta|<0.8$ with a
minimum constraint of 5 charged particles with $p_{\rm{T}}>$ 0.15 GeV/$c$ to
recreate conditions similar to ALICE experiment at the LHC. The low-$S_{0}$
events are those events having spherocity values in the lowest 20% and the
high-$S_{0}$ events are those occupying the highest 20% in the spherocity
distribution. The spherocity cuts are shown in Table 1.

Table 1: Lowest 20 % (low-$S_{0}$) and highest 20% (high-$S_{0}$) cuts on spherocity distribution. Centrality (%) |  low-$S_{0}$ |  high-$S_{0}$
---|---|---
0-10 |  0-0.88025 |  0.95285-1
40-50 |  0-0.71475 |  0.86365-1
60-70 |  0-0.70775 |  0.87355-1

Figure 1 shows the transverse momentum space correlation ($p_{\rm y}$ vs.
$p_{\rm x}$) for $\pi^{+}+\pi^{-}$ in different spherocity classes for (0-10)%
(top) and (40-50)% (bottom) central Pb-Pb collisions at
$\sqrt{s_{\rm{NN}}}=5.02$ TeV in AMPT model. The spherocity integrated events
for (40-50)% centrality class clearly show anisotropy in the momentum space
which is linked to the initial state geometrical anisotropy. However such an
anisotropy is not observed in the (0-10)% centrality class. This behavior is
understood as semi-central collisions develop an almond shape of nuclear
overlap region in the initial stage while for most central collisions the
nuclear overlap region is nearly circular. When one studies the transverse
momentum space correlation in different spherocity classes, the low-$S_{0}$
events show higher anisotropy compared to spherocity-integrated events while
high-$S_{0}$ events show nearly zero anisotropy. This is one of the
testimonies that spherocity is a unique tool for the event classification to
study the particle production in heavy-ion collisions.

Figure 2: (Color Online) Top plot: $p_{\rm T}$-spectra for pions in (0-10)%
central Pb-Pb collisions for high-$S_{0}$, $S_{0}$-integrated and low-$S_{0}$
classes. Bottom plot: Ratio of $p_{\rm T}$-spectra for high-$S_{0}$ and
low-$S_{0}$ events to the $S_{0}$-integrated events.

###  II.3 Elliptic Flow

In non-central heavy-ion collisions, the initial spatial anisotropies of the
nuclear overlap region gets converted to the final state momentum space
azimuthal anisotropy of the produced hadrons due to the strong differential
pressure gradients of the QGP medium. The observed non-zero finite azimuthal
anisotropy (evident from Fig. 1) can be expressed as a Fourier series in the
azimuthal angle, $\phi$:

$\displaystyle E\frac{d^{3}N}{dp^{3}}=\frac{d^{2}N}{2\pi p_{\rm T}dp_{\rm
T}dy}\bigg{(}1+2\sum_{n=1}^{\infty}v_{n}\cos[n(\phi-\psi_{n})]\bigg{)}\,.$ (2)

Here, $v_{n}$ is the magnitude of anisotropic flow of different order $n$, and
$\psi_{n}$ is the n${}^{\text{th}}$ harmonic event plane angle v2eventplane .
$v_{n}$ is usually studied as a function of transverse momentum ($p_{\rm T}$)
and pseudorapidity ($\eta$).

Figure 3: (Color Online) One dimensional azimuthal correlation of pions for
low-$S_{0}$ (red triangles), high-$S_{0}$ (blue squares) and spherocity-
integrated (black circles) events in 40-50% central Pb-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV using AMPT model.  Figure 4: (Color Online) Two
particle elliptic flow coefficients of pions for low-$S_{0}$ (red triangles),
high-$S_{0}$ (blue squares) and spherocity-integrated (black circles) events
in 40-50% central Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV using AMPT
model.  Figure 5: (Color Online) Single particle elliptic flow coefficients of
pions for low-$S_{0}$ (red triangles), high-$S_{0}$ (blue squares) and
spherocity-integrated (black circles) events in 40-50% central Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV using AMPT model.

Figure 6: (Color Online) Centrality dependence of elliptic flow co-efficient
($v_{2}(p_{\rm T}^{a}$)) as a function of $p_{\rm T}^{a}$ of pions (blue
circles), kaons (black squares) and protons (red triangles) for spherocity-
integrated (top panel) and low-$S_{0}$ (bottom panel) events in Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV using AMPT model.

For $n$ = 2, Eq. 2 gives the second order harmonics in the expansion and its
coefficient, $v_{2}$ provides the numerical measure of the elliptic flow which
is defined as:

$\displaystyle v_{2}=\langle\cos(2(\phi-\psi_{2}))\rangle$ (3)

Due to the almond shaped nuclear overlap region in a non-central heavy-ion
collision, the azimuthal anisotropy is expected to be dominated by the $v_{2}$
component. In addition, the signatures of $v_{1}$ (directed flow) and $v_{3}$
(triangular flow) are observed but their contribution is quite less as
compared to $v_{2}$. Another interesting observable owing to anisotropic flow
and properties of the medium is the correlation function between two particles
in relative pseudorapidity ($\Delta\eta=\eta_{a}-\eta_{b}$) and azimuthal
angle ($\Delta\phi=\phi_{a}-\phi_{b}$). The labels $a$ and $b$ denote the two
particles in the pair, which are chosen from different (or same) transverse
momentum intervals. In this study, we have taken particles in $|\eta|<$ 2.5 to
include wider range of particles mostly contributing to $v_{2}$ and $p_{\rm
T}>$ 0.5 GeV/c to be consistent with Refs. Aad:2015lwa ; Mallick:2021rsd . The
two particle correlation function, $C(\Delta\eta,\Delta\phi)$ can be
constructed as the ratio of distributions for same-event pairs,
$S(\Delta\eta,\Delta\phi)$ and mixed-event pairs, $B(\Delta\eta,\Delta\phi)$,
which is given by:

$\displaystyle
C(\Delta\eta,\Delta\phi)=\frac{S(\Delta\eta,\Delta\phi)}{B(\Delta\eta,\Delta\phi)}$
(4)

We have chosen five randomly selected events for the mixed-event pairing and
hence the pairs contain no physical correlations whatsoever. To obtain the
anisotropy in azimuthal direction, our analysis focuses mainly on the shape of
the 1D correlation function $\Delta\phi$, given as;

$\displaystyle C(\Delta\phi)=\frac{dN_{\rm
pairs}}{d\Delta\phi}=A\times\frac{\int
S(\Delta\eta,\Delta\phi)d\Delta\eta}{\int
B(\Delta\eta,\Delta\phi)d\Delta\eta}.$ (5)

Here, the normalization constant $A$ ensures that the number of pairs to be
the same between the signal ($S$) and background ($B$) events for a given
$\Delta\eta$ interval. The $\Delta\eta$ interval is chosen carefully by
excluding the jet peak region observed in the $C(\Delta\eta,\Delta\phi)$
distribution. For our case, the interval is taken to be
$2.0<|\Delta\eta|<4.8$. By choosing particle pairs with a proper
pseudorapidity cut, the two-particle correlation method removes the residual
non-flow effects in the estimation of elliptic flow. Non-flow azimuthal
correlations usually arise from jets and short range resonance decays, which
are not associated with the symmetry planes at the early stage of the
collisions. Although the non-flow effects are very small in heavy-ion
collisions at very high energies, but our study based on event shape
dependence being sensitive to jets may get affected by these effects and hence
we proceed with this method. To ensure the consistency of the method, one can
have a look at Fig. 3 of Ref. Mallick:2021rsd , where the comparison of
calculated $v_{2,2}$ between PYTHIA8 (Angantyr) Bierlich:2018xfw  and AMPT
model is shown. It is shown that the method removes residual non-flow
contributions without affecting the event shape studies based on transverse
spherocity. Now, the pair distribution ($N_{\rm pairs}$) can be expanded in
$\Delta\phi$ into a Fourier series:

$\displaystyle\frac{dN_{\rm
pairs}}{d\Delta\phi}\propto\bigg{[}1+2\sum_{n=1}^{\infty}v_{n,n}(p_{\rm
T}^{a},p_{\rm T}^{b})\cos n\Delta\phi\bigg{]}.$ (6)

Here, $v_{n,n}$ is called as the two particle flow coefficient. By rewriting
Eq. 6 with the help of Eq. 5 one obtains:

$\displaystyle
C(\Delta\phi)\propto\bigg{[}1+2\sum_{n=1}^{\infty}v_{n,n}(p_{\rm T}^{a},p_{\rm
T}^{b})\cos n\Delta\phi\bigg{]}.$ (7)

Figure 7: (Color Online) Centrality dependence of elliptic flow co-efficient
with number of quark participant scaling ($v_{2}(p_{\rm T}^{a}$)/$n_{q}$) as a
function of $p_{\rm T}^{a}$ of pions (blue circles), kaons (black squares) and
protons (red triangles) for spherocity-integrated (top panel) and low-$S_{0}$
(bottom panel) events in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV
using AMPT model.

Now, $v_{n,n}$ can be calculated as:

$\displaystyle v_{n,n}(p_{\rm T}^{a},p_{\rm T}^{b})=\langle
cos(n\Delta\phi)\rangle$ (8)

We have taken very narrow bins in $-\pi/2<\Delta\phi<3\pi/2$. $v_{n,n}$ are
symmetric functions with respect to $p_{\rm T}^{a}$ and $p_{\rm T}^{b}$. The
definition of harmonics defined in Eq. 2 also enters in Eq. 6, which can be
rewritten as;

$\displaystyle\frac{dN_{\rm
pairs}}{d\Delta\phi}\propto\bigg{[}1+2\sum_{n=1}^{\infty}v_{n}(p_{\rm
T}^{a})v_{n}(p_{\rm T}^{b})\cos n\Delta\phi\bigg{]}.$ (9)

By this definition, the event plane angle drops out in convolution. So, if the
azimuthal anisotropy is driven by collective expansion then $v_{n,n}$ should
factorize into the product of two single-particle harmonic coefficients.

$\displaystyle v_{n,n}(p_{\rm T}^{a},p_{\rm T}^{b})=v_{n}(p_{\rm
T}^{a})v_{n}(p_{\rm T}^{b}).$ (10)

From Eq. 10 we derive the single particle flow coefficient as:

$\displaystyle v_{n}(p_{\rm T}^{a})=v_{n,n}(p_{\rm T}^{a},p_{\rm
T}^{b})/\sqrt{v_{n,n}(p_{\rm T}^{b},p_{\rm T}^{b})}$ (11)

Figure 8: (Color Online) Centrality dependence of elliptic flow co-efficient
with transverse mass ($m_{\rm T}$) scaling as a function of $p_{\rm T}^{a}$ of
pions (blue circles), kaons (black squares) and protons (red triangles) for
spherocity-integrated (top panel) and low-$S_{0}$ (bottom panel) events in Pb-
Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV using AMPT model.

We now proceed for the estimation of azimuthal anisotropy of
$\pi^{+}+\pi^{-}$, $K^{+}+K^{-}$ and $p+\bar{p}$ in different spherocity
classes for Pb-Pb collisions at the LHC energies from AMPT. For better
readability, here onwards we refer $\pi^{+}+\pi^{-}$, $K^{+}+K^{-}$ and
$p+\bar{p}$ as pions, kaons and protons, respectively.

##  III Results and Discussions

Top panel of Fig. 2 shows the $p_{\rm T}$-spectra for pions in (0-10)% Pb-Pb
collisions for high-$S_{0}$, $S_{0}$-integrated and low-$S_{0}$ events and
bottom panel shows the ratio of $p_{\rm T}$-spectra for high-$S_{0}$ and
low-$S_{0}$ events to the $S_{0}$-integrated events. In the low-$p_{\rm T}$
region, we observe that pion production is dominated by high-$S_{0}$ events.
But as we move on to higher $p_{\rm T}$, the relative yield from high-$S_{0}$
and low-$S_{0}$ events are similar. These results show that spherocity
successfully differentiates the heavy-ion collisions event topology based on
their geometrical shapes i.e. high-$S_{0}$ and low-$S_{0}$.

Figure 3 shows one dimensional azimuthal correlation of pions for low-$S_{0}$,
high-$S_{0}$ and spherocity-integrated events in 40-50% central Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. As described in the previous
section, the correlation function is calculated in the $p_{\rm T}$ range of
0.5 to 5 GeV/c for $2.0<|\Delta\eta|<4.8$. We see a strong correlation of
modulation magnitude with spherocity in Fig. 3. Also, double peak is observed
for away side ($\Delta\phi\sim\pi$) for high-$S_{0}$ events, which may arise
from the residual contribution of the triangular flow. The observed behavior
is analogous to the behavior observed in Refs. Aad:2015lwa , where flow vector
is used as an event classifier for heavy-ion collisions. Figures 4 and 5 show
the two particle and single particle elliptic flow coefficients of pions for
low-$S_{0}$, high-$S_{0}$ and spherocity-integrated events in 40-50% central
Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV obtained using Eqs. 8 and 11,
respectively. We observe that the contribution towards elliptic flow
coefficients is dominated by low-$S_{0}$ events while high-$S_{0}$ events show
nearly zero elliptic flow coefficients. In a similar fashion, we obtain the
elliptic flow coefficients for pions, kaons and protons and then study the
number of constituent quark scaling (NCQ) of the elliptic flow coefficients.

Figure 6 shows the single particle elliptic flow co-efficient ($v_{2}(p_{\rm
T}^{a}$)) as a function of $p_{\rm T}^{a}$ for pions, kaons and protons in
(0-10)%, (40-50)% and (60-70)% central Pb-Pb collisions at $\sqrt{s_{\rm
NN}}=5.02$ TeV for spherocity-integrated and low-$S_{0}$ classes. The values
of $v_{2}(p_{\rm T}^{a}$) increases from (0-10)% to (40-50)% central
collisions while slightly decreases for (60-70)% central collisions for all
particle species. The increasing behavior is consistent with the picture of
geometry of the collision driving the final state anisotropy. The slight
decrease for (60-70)% central collisions might have originated from an
interplay of different effects like the smaller lifetime of the fireball in
peripheral collisions with respect to most central collisions which does not
allow the elliptic flow to further develop and less significant contribution
of eccentricity fluctuations and to final state hadronic effects Song:2013qma
. For (40-50)% and (60-70)% central collisions, we see a clear mass hierarchy
in the low-$p_{\rm T}$ region ($p_{\rm T}<$ 2 GeV/c) and in the
intermediate-$p_{\rm T}$ the elliptic flow for protons are higher than that of
pions and kaons. As per hydrodynamic calculations, the mass hierarchy in the
low-$p_{\rm T}$ could be an interplay between radial and elliptic flow. The
baryon and meson separation in the intermediate-$p_{\rm T}$ could arise from
the quark coalescence mechanism of AMPT string melting mode. As described
earlier, in the string melting mode of AMPT, when a quark and an anti-quark
are close in phase space with their momenta very close to each other, they
coalesce to form a meson. When three quarks come closer in phase space, they
combine to form a baryon. Thus the probability of recombination of three
quarks becomes higher in the intermediate-$p_{\rm T}$ range.

Owing to the quark coalescence mechanism, one may expect elliptic flow to
follow the constituent quark number (NCQ) scaling behavior. If elliptic flow
follows NCQ-scaling, then it can be expressed as,

$\displaystyle v_{2}^{h}(p_{\rm T}^{a})=n_{q}\times v_{2}^{q}(p_{\rm
T}^{a}/n_{q}).$ (12)

Here, $n_{q}$ is the number of constituent quarks for h-hadron. In RHIC
experimental observations Adams:2003am ; Abelev:2007qg ; Adler:2003kt , it was
found that the identified particles follow the NCQ scaling behavior. The
PHENIX collaboration extends the scaling property to the lower transverse
momentum by obtaining elliptic flow as a function of the transverse kinetic
energy, where the kinetic energy is defined by $m_{\rm T}-m_{0}$ Adler:2003kt
; Afanasiev:2007tv ; Adare:2006ti . Here, $m_{\rm T}$ is the transverse mass
($\sqrt{p_{\rm T}^{2}+m_{0}^{2}}$) and $m_{0}$ is the rest mass of the
particle. The kinetic energy representation of scaling was found to be valid
in the RHIC results. Also, this scaling behavior was successfully reproduced
by models which invokes quark coalescence as the dominant hadronization
mechanism Molnar:2003ff ; Greco:2003mm ; Fries:2003kq ; Hwa:2003ic . Thus, the
NCQ scaling of the elliptic flow has been interpreted as evidence of the
dominance of quark degrees of freedom in the initial stages of heavy-ion
collisions Molnar:2003ff ; Greco:2003mm ; Fries:2003kq ; Hwa:2003ic . In this
work, we have studied both the NCQ-scaling behavior for pions, kaons and
protons. Figure 7 shows centrality dependence of elliptic flow co-efficient
with number of quark participant scaling ($v_{2}(p_{\rm T}^{a}$)/$n_{q}$) as a
function of $p_{\rm T}^{a}/n_{q}$ of pions, kaons and protons for spherocity-
integrated and low-$S_{0}$ events. Figure 8 shows the transverse kinetic
energy, $v_{2}(p_{\rm T}^{a}$)/$n_{q}$ scaling as a function of $(m_{\rm
T}-m_{0})/n_{q}$. In the intermediate transverse momentum range, where the
quark coalescence mechanism is dominant, both the scaling of $v_{2}(p_{\rm
T}^{a}$)/$n_{q}$ show a significant deviation. Results from LHC suggests that
both NCQ and transverse kinetic energy scaling of $v_{2}(p_{\rm
T}^{a}$)/$n_{q}$ at the LHC energies show a deviation of 20% Abelev:2014pua ;
v2_ALICE2  for intermediate/high transverse momenta. The observed deviations
seems to be similar for all centrality classes.

Let us now discuss how the results vary when such a study is performed in
different spherocity classes. In Fig. 6, for (0-10)% centrality we see that
the values of $v_{2}$ is similar in both spherocity-integrated and low-$S_{0}$
events. For (40-50)% and (60-70)% centrality classes we see significantly
higher $v_{2}$ in low-$S_{0}$ events compared to the spherocity-integrated
events for all particles. For (40-50)% centrality class, the baryon-meson
separation in intermediate-$p_{\rm T}$ is more prominent in low-$S_{0}$ events
compared to spherocity-integrated events. While studying the NCQ scaling and
transverse kinetic energy scaling of the elliptic flow for different
spherocity classes in Figs. 7 and 8, it was found that the low-$S_{0}$ events
show a larger deviation compared to spherocity-integrated events in all
centrality classes. The deviation is highly significant, where one can draw a
conclusion of violation in the NCQ scaling of elliptic flow at the LHC
energies in AMPT model. This essentially indicates that using spherocity not
only one can study the events with enhanced elliptic flow but also one can
directly probe into the scaling properties of the elliptic flow. Thus,
spherocity can separate out the events where the initial stage of the systems
is dominated by the quark degrees of freedom from total events. Such an
observation strengthens our previous publications using simulation from AMPT
model Mallick:2021rsd ; Mallick:2021wop , where it was argued that the
transverse spherocity can be used as an event classifier in experimental
analyses for heavy-ion collisions at RHIC and LHC.

##  IV Summary

In summary,

  * •

The transverse momentum space correlation as a function of transverse
spherocity indicates a strong correlation of transverse spherocity with the
azimuthal anisotropy in the produced system in heavy-ion collisions. This
correlation is further confirmed with the study of transverse spherocity
dependent one dimensional azimuthal correlation.

  * •

The contribution towards elliptic flow coefficients is dominated by
low-$S_{0}$ events while high-$S_{0}$ events show nearly zero elliptic flow.

  * •

The baryon-meson separation at the intermediate transverse momentum, the NCQ
scaling and the transverse kinetic energy scaling of $v_{2}(p_{\rm
T}^{a}$)/$n_{q}$ from AMPT show a similar behavior for pions, kaons and
protons as observed at the LHC.

  * •

The baryon-meson separation in intermediate-$p_{\rm T}$ was found to be more
prominent in low-$S_{0}$ events compared to spherocity-integrated events.

  * •

The study on the NCQ scaling and the transverse kinetic energy scaling of the
elliptic flow in different spherocity classes shows that the low-$S_{0}$
events violate the scaling properties of the elliptic flow. This essentially
indicates that using spherocity not only one can study the events with
enhanced elliptic flow but also one can directly probe into the scaling
properties of the elliptic flow.

One should note that these results are obtained from an event generator (AMPT)
which successfully describes the elliptic flow for spherocity-integrated
events in heavy-ion collisions. Thus it is expected that the observed
behaviours to mimic the experimental data. However, any deviation of
experimental results from AMPT results on elliptic flow could be very
interesting to understand the system dynamics. We believe that in coming days
of Run-3 at the LHC, transverse spherocity can serve as an important event
classifier which could give a better insight to the possible QGP formation in
heavy-ion collisions.

## Acknowledgements

R.S. acknowledges the financial supports under the CERN Scientific
Associateship and the financial grants under DAE-BRNS Project No.
58/14/29/2019-BRNS of the Government of India. S.T. acknowledges the supports
under the INFN postdoctoral fellowship. The authors would like to acknowledge
the usage of resources of the LHC grid Tier-3 computing facility at IIT
Indore.

## References

  * (1) J. Adam et al. [ALICE Collaboration], Nature Phys. 13, 535 (2017).
  * (2) V. Khachatryan et al. [CMS Collaboration], Phys. Lett. B 765, 193 (2017).
  * (3) E. Cuautle, R. Jimenez, I. Maldonado, A. Ortiz, G. Paic and E. Perez, arXiv:1404.2372 [hep-ph].
  * (4) A. Ortiz, G. Paic and E. Cuautle, Nucl. Phys. A 941, 78 (2015).
  * (5) G. P. Salam, Eur. Phys. J. C 67, 637 (2010).
  * (6) G. Bencédi [ALICE Collaboration], Nucl. Phys. A 982, 507 (2019)
  * (7) A. Banfi, G. P. Salam and G. Zanderighi, JHEP 1006, 038 (2010).
  * (8) S. Tripathy, A. Bisht, R. Sahoo, A. Khuntia and M. P. Salvan, Adv. High Energy Phys. 2021, 8822524 (2021).
  * (9) A. Khuntia, S. Tripathy, A. Bisht and R. Sahoo, J. Phys. G 48, 035102 (2021).
  * (10) N. Mallick, R. Sahoo, S. Tripathy and A. Ortiz, J. Phys. G 48, 045104 (2021).
  * (11) S. Voloshin and Y. Zhang, Z. Phys. C 70, 665 (1996).
  * (12) B. B. Abelev et al. [ALICE Collaboration], JHEP 06, 190 (2015).
  * (13) I. Arsene et al. [BRAHMS Collaboration], Nucl. Phys. A 757, 1 (2005).
  * (14) K. Adcox et al. [PHENIX Collaboration], Nucl. Phys. A 757, 184 (2005).
  * (15) B. B. Back et al. [PHOBOS Collaboration], Nucl. Phys. A 757, 28 (2005).
  * (16) J. Adams et al. [STAR Collaboration], Nucl. Phys. A 757, 102 (2005).
  * (17) T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Lett. B 636, 299 (2006).
  * (18) S. A. Voloshin, Nucl. Phys. A 715, 379 (2003).
  * (19) J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 92, 052302 (2004).
  * (20) S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 91, 182301 (2003).
  * (21) K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 105, 252302 (2010).
  * (22) D. Molnar and S. A. Voloshin, Phys. Rev. Lett. 91, 092301 (2003).
  * (23) V. Greco, C. M. Ko and P. Levai, Phys. Rev. C 68, 034904 (2003).
  * (24) R. J. Fries, B. Muller, C. Nonaka and S. A. Bass, Phys. Rev. C 68, 044902 (2003).
  * (25) N. Mallick, S. Tripathy, A. N. Mishra, S. Deb and R. Sahoo, Phys. Rev. D (2021: In Press), arXiv:2103.01736 [hep-ph].
  * (26) Z. W. Lin, C. M. Ko, B. A. Li, B. Zhang and S. Pal, Phys. Rev. C 72, 064901 (2005).
  * (27) X. N. Wang and M. Gyulassy, Phys. Rev. D 44, 3501 (1991).
  * (28) B. Zhang, Comput. Phys. Commun. 109, 93 (1998).
  * (29) Z. w. Lin and C. M. Ko, Phys. Rev. C 65, 034904 (2002).
  * (30) Y. He and Z. W. Lin, Phys. Rev. C 96, 014910 (2017).
  * (31) B. Li, A. T. Sustich, B. Zhang and C. M. Ko, Int. J. Mod. Phys. E 10, 267 (2001).
  * (32) B. A. Li and C. M. Ko, Phys. Rev. C 52, 2037 (1995).
  * (33) V. Greco, C. M. Ko and P. Levai, Phys. Rev. C 68, 034904 (2003).
  * (34) R. J. Fries, B. Muller, C. Nonaka and S. A. Bass, Phys. Rev. Lett. 90, 202303 (2003).
  * (35) R. J. Fries, B. Muller, C. Nonaka and S. A. Bass, Phys. Rev. C 68, 044902 (2003).
  * (36) S. Tripathy, S. De, M. Younus and R. Sahoo, Phys. Rev. C 98, 064904 (2018).
  * (37) C. Loizides, J. Kamin and D. d’Enterria, Phys. Rev. C 97, 054910 (2018) Erratum: [Phys. Rev. C 99, 019901 (2019)]
  * (38) B. B. Abelev et al. [ALICE Collaboration], JHEP 1506, 190 (2015).
  * (39) G. Aad et al. [ATLAS Collaboration], Phys. Rev. C 92, 034903 (2015).
  * (40) C. Bierlich, G. Gustafson, L. Lönnblad and H. Shah, JHEP 10, 134 (2018).
  * (41) H. Song, S. Bass and U. W. Heinz, Phys. Rev. C 89, 034919 (2014).
  * (42) B. I. Abelev et al. [STAR Collaboration], Phys. Rev. C 75, 054906 (2007).
  * (43) S. Afanasiev et al. [PHENIX Collaboration], Phys. Rev. Lett. 99, 052301 (2007).
  * (44) A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 98, 162301 (2007).
  * (45) R. C. Hwa and C. B. Yang, Phys. Rev. C 70, 024904 (2004).


# GD-CAF: Graph Dual-stream Convolutional Attention Fusion for Precipitation
Nowcasting

Lóránd Vatamány  Siamak Mehrkanoon<EMAIL_ADDRESS>Department of
Information and Computing Sciences, Utrecht University, Utrecht, The
Netherlands

###### Abstract

Accurate precipitation nowcasting is essential for various purposes, including
flood prediction, disaster management, optimizing agricultural activities,
managing transportation routes and renewable energy. While several studies
have addressed this challenging task from a sequence-to-sequence perspective,
most of them have focused on a single area without considering the existing
correlation between multiple disjoint regions. In this paper, we formulate
precipitation nowcasting as a spatiotemporal graph sequence nowcasting
problem. In particular, we introduce Graph Dual-stream Convolutional Attention
Fusion (GD-CAF), a novel approach designed to learn from historical
spatiotemporal graph of precipitation maps and nowcast future time step ahead
precipitation at different spatial locations. GD-CAF consists of spatio-
temporal convolutional attention as well as gated fusion modules which are
equipped with depthwise-separable convolutional operations. This enhancement
enables the model to directly process the high-dimensional spatiotemporal
graph of precipitation maps and exploits higher-order correlations between the
data dimensions. We evaluate our model on seven years of precipitation maps
across Europe and its neighboring areas collected from the ERA5 dataset,
provided by Copernicus. The model receives a fully connected graph in which
each node represents historical observations from a specific region on the
map. Consequently, each node contains a 3D tensor with time, height, and width
dimensions. Experimental results demonstrate that the proposed GD-CAF model
outperforms the other examined models. Furthermore, the averaged seasonal
spatial and temporal attention scores over the test set are visualized to
provide additional insights about the strongest connections between different
regions or time steps. These visualizations shed light on the decision-making
process of our model.

###### keywords:

Precipitation nowcasting , High dimensional graph precipitation data , Graph
Attention Networks , Deep Learning

††journal: arXiv

##  1 Introduction

Precipitation nowcasting involves forecasting the forthcoming intensity of
rainfall typically on a timescale ranging from minutes to a few hours.
Nowcasting can help the operations of several weather dependent sectors
including energy management, retail, flood, traffic control and emergency
services [1]. To serve these sectors effectively, the accuracy of nowcasting
must extend across a range of spatial and temporal scales. Two primary
approaches are commonly employed for precipitation nowcasting. The first one
involves ensemble numerical weather prediction (NWP) systems, which rely on
the physical properties of the atmosphere to generate multiple realistic
precipitation forecasts. However, these methods are not suitable for short-
term predictions due to their high computational expense, sensitivity to
noise, and dependence on the initial conditions of the event [2]. The second
approach is optical flow, a technique for deriving a velocity field from
consecutive images. This method typically involves two steps: feature tracking
and extrapolation [3]. Unlike supervised deep learning approaches, optical
flow methods are unsupervised algorithms and often serve as baseline
predictions [4].

In contrast to NWP models, data-driven approaches do not rely on the physical
properties of the atmosphere. Instead, they utilize historical weather
observations to train models capable of mapping input data to target outputs
[5]. Among these data-driven models, deep neural network architectures stand
out, as they are trained in an end-to-end fashion and possess the ability to
extract complex underlying patterns from data by incorporating multiple
nonlinear layers. Recent advances in deep learning have showcased remarkable
progress in the fields of weather element forecasting and related nowcasting
tasks [6, 7, 8, 9, 10, 11, 12].

In particular, Convolutional Neural Networks (CNNs) based models have already
demonstrated success in addressing the weather forecasting challenge [13, 14].
However, it’s worth noting that CNN-based methods typically do not account for
the spatial relationships between weather stations. In a prior study [5], the
approach involved transforming historical data into a tensor format
(comprising weather stations, weather variables, and time steps), which was
subsequently fed into the model, with convolution operations applied across
the data volume. Consequently, the neighborhood relationships between weather
stations were primarily determined by their order in the dataset rather than
explicitly considering their spatial proximity.

Weather patterns are inherently spatial, with various meteorological factors
interacting across geographical regions. Therefor, graph neural networks
(GNNs) based models that can generalize CNNs to work on graphs rather than on
regular grids are among promising architecture for the weather elements
nowcasting. In particular, GNNs can capture the intricate spatial dependencies
by modeling data as a graph, where nodes represent locations or weather
stations, and edges represent the connections between them. This enables GNN
based models to account for the influence of neighboring regions on each
other’s weather conditions, making them effective at modeling spatial
correlations. In addition, GNNs can also be extended to incorporate temporal
information, allowing them to model how weather conditions change over time.
This is crucial for short-term weather predictions and nowcasting.

In this paper we propose a novel Graph Dual-stream Convolutional Attention
Fusion (GD-CAF)111https://github.com/wendig/GD-CAF which encompasses
spatiotemporal convolutional attention and gated fusion modules, enriched with
depthwise-separable convolutional operations. This augmentation empowers the
model to effectively capture and leverage the inherent correlations and
dependencies present in the spatiotemporal graph sequences, leading to
improved nowcasting performance. In contrast to many other graph based models
such as [15, 16] that require the node features to be one dimensional, the
proposed GD-CAF model can directly analyze high-dimensional, tensorial node
features within the spatiotemporal graph of precipitation maps. This
versatility enables GD-CAF to work with richer and more complex data
representations, enhancing its ability to make accurate predictions in the
context of weather nowcasting. Comprehensive experiments and comparative
evaluations demonstrate that the GD-CAF model surpasses the widely-used SmaAt-
UNet [9] and persistence model [17] in precipitation nowcasting task,
highlighting the significance of this contribution in advancing the field and
addressing the challenges associated with high-dimensional data in weather
forecasting/nowcasting.

This paper is organized as follows. A brief overview of the related research
works is given in Section 2. Section 3 introduces the proposed GD-CAF model.
The experimental settings and description of the used datasets are given in
Section 4. The obtained results are discussed in Section 5 and the conclusion
is drawn in Section 6.

##  2 Related Work

Figure 1: The input is a spatiotemporal graph sequence, and the output is a
spatiotemporal graph.

Weather element forecasting based on deep-learning architectures has recently
gained a lot of attention due to the availability of large amount of weather
data and the rapid advances in neural network techniques. The literature has
already witnessed successful application of different architectures including
Recurrent Neural Network (RNN) [18], Long short-term memory (LSTM) [19],
Convolutional LSTM (ConvLSTM) [20], Convolutional Neural Network (CNN),
encoder-decoder [21], UNet [22] and graph neural networks [23] in weather
forecasting and nowcasting related tasks.

For instance, the authors in [20], introduced a convolutional LSTM model to
predict future rainfall intensity in Hong Kong over a relatively short period.
In [4], the authors proposed Convcast, an embedded convolutional LSTM-based
architecture. Additionally, [24] introduced a dynamic convolutional layer for
short-range weather prediction. Furthermore, in [25], a deep convolutional
neural network is employed for predicting thunderstorms and heavy rains. A
CNN-based wind speed prediction model that effectively captures spatiotemporal
patterns in wind data using real weather datasets from Denmark and the
Netherlands is presented in [14]. The author in [5] proposed different CNN
architectures, including 1-D, 2-D, and 3-D convolutions, to accurately predict
wind speed and temperature for few hours ahead.

The UNet architecture, initially a successful model primarily used in the
field of medical image analysis, has found application in precipitation
nowcasting as well, as demonstrated in the study by Lebedev et al. [26]. In a
subsequent work by Trebing et al. [9], a model known as SmaAt-UNet was
introduced as an extension of the core UNet model. This extended model
significantly reduces the number of parameters in the UNet without
compromising its performance, as outlined in [9]. Furthermore, in another
study by Fernandez et al. [10], a modification called Broad-UNet was
introduced, which enhances the UNet architecture by incorporating asymmetric
parallel convolutions and the Atrous Spatial Pyramid Pooling (ASPP) module.

Although the previous works exploit spatio-temporal correlations, they do not
fully leverage the spatial information from multiple weather stations
(regions). Graph neural networks (GNNs) [27] have recently attracted a lot of
attention due to their expressive power and ability to infer information from
complex data, such as brain signals, social network interactions, and weather
prediction [28, 29, 30, 31]. They propagate information through the graph
nodes and edges, enabling the model to capture the underlying structure and
dependencies within the graph. For instance, the authors in [32] utilized a
graph neural network for global weather forecasting, where the system learns
to project the current 3D atmospheric state six hours ahead, yielding improved
results. The authors in [33] used graph neural networks to predict power
outages based on current weather conditions.

One category of Graph Neural Networks (GNNs) includes Graph Convolution
Networks (GCNs) as described in [34]. GCNs extend the capabilities of
Convolutional Neural Networks (CNNs) to operate on graphs rather than regular
grids. They are particularly adept at integrating neighbor relationships,
often through the adjacency matrix of a graph. In the study by Stanczyk et al.
[30], the authors applied graph convolutional networks (GCNs) to tackle the
challenge of wind speed prediction using data from multiple weather stations.
Their model outperformed existing baseline methods when tested with real
datasets from weather stations in Denmark and the Netherlands. Similarly, in
[23], the authors introduced GCLSTM and GCTrafo, graph convolutional models
used to address solar power generation forecasting from multi-site
photovoltaic production data represented as signals on a graph. These models,
solely reliant on production data, surpassed existing multi-site forecasting
methods, especially for a six-hour prediction horizon.

Another category of GNNs comprises Graph Attention Networks (GAT) [35], which
are designed to work with graph data and leverage attention mechanisms. An
extension of GAT, known as the Heterogeneous Graph Attention Network, was
introduced by Wang et al. in [36]. This approach handles the complexities
associated with heterogeneous graphs containing different types of nodes and
links by incorporating node-level and semantic-level attentions. Aykas et al.
extended GAT further in [31], introducing Multistream Graph Attention
Networks. This model incorporates a learnable adjacency matrix and a novel
attention mechanism, which they applied to predict wind speeds for multiple
cities. In [37], a Hierarchical Graph Attention Network was proposed by the
authors to capture dependencies at both the object-level and triplet-level,
allowing the model to represent interactions between objects and dependencies
among relation triplets. For spatial-temporal analysis, Spatial-Temporal Graph
Attention Networks were introduced in [38], where graph attention mechanisms
were used to capture spatial dependencies among road segments, and LSTM
networks were employed to extract temporal features.

##  3 Methods

###  3.1 Preliminaries

$\bullet$ Graph precipitation maps:

We denote a network of precipitation maps as a fully connected graph
$G=(V,E)$. Here, V is a set of $N=|V|$ vertices (nodes), and $E$ represents
the edges connecting them. Each node holds historical observations from a
particular region, see the input graph in Fig. 1. Given the precipitation maps
with $H$ (height) and $W$ (width) dimensions at $N$ vertices over $T$ time
steps, we denote one sample as $\omega$ in the dataset $\mathcal{D}$ with
input $X$ and output $Y$ as follows:

$\omega=\biggl{\\{}X,Y\biggr{\\}}=\biggl{\\{}x_{v_{i},t_{j}},y_{v_{i},\,t_{\Delta}}\biggr{\\}}_{i=1,\,j=1}^{N,T},$
(1)

where $x_{v_{i},t_{j}}\in\mathbb{R}^{H\times W}$ is a precipitation map at
time step $t_{j}$ for node $v_{i}$. The target
$y_{v_{i},\,t_{\delta}}\in\mathbb{R}^{H\times W}$ is the precipitation map in
the future time step ($t_{\Delta=T+\delta}$) for node $v_{i}$. Given the input
graph, the goal is to predict the precipitation maps of all nodes for a single
time step into the future as illustrated in Fig. 1.

###  3.2 Proposed Model

Here, we introduce the Graph Dual-stream Convolutional Attention Fusion (GD-
CAF) model that leverages higher-order correlations among the node dimensions
of the historical spatiotemporal graph of precipitation maps and nowcast
precipitation for a time step ahead at various spatial locations. GD-CAF is
composed of spatio-temporal convolutional attention and gated fusion modules,
both of which incorporate depthwise-separable convolutional operations. Unlike
other competing graph-based models such as [15, 16, 36] that necessitate low-
dimensional node input representation, our proposed model can directly process
high-dimensional nodes. As a result, the fundamental structure of the graph
information remains unchanged.

Figure 2: (a) GD-CAF architecture overview (b) Double convolution block with
pooling (c) ST-Attention block in the $\ell$-th block. (d) Double convolution
block with upsampling. The numbers above the double convolutional blocks
indicate the number of input and output channels respectively.

####  3.2.1 Architecture overview

An overview of the proposed GD-CAF model is illustrated in Fig. 2 (a). The
input to the model, i.e. $X=\left\\{x_{v_{i},t_{j}}\in\mathbb{R}^{H\times
W}\right\\}_{i=1,j=1}^{N,T}$, consists of historical observations at $N$ nodes
for $T$ time steps. The input $X$ is initially transformed into
$\tilde{X}^{(0)}$ using a double convolutional operation, during which the
temporal depth increases $K$ times, corresponding to the number of heads used
in the attention mechanisms. Therefore, $\tilde{X}^{(0)}$ is defined as
$\left\\{\tilde{x}_{v_{i},t_{j}}\in\mathbb{R}^{H\times
W}\right\\}_{i=1,j=1}^{N,T\times K}$. Next, $\tilde{X}^{(0)}$ is passed
through a sequence of $L$ ST-Attention blocks, producing the output
$\tilde{X}^{(L)}$ at the $L$-th block, with the same dimension as
$\tilde{X}^{(0)}$. The depth dimension ($T\times K$) of the output of the ST-
Attention blocks are then reduced using a double convolutional operation to
obtain a single time step for all nodes, represented as
$\hat{Y}=\left\\{y_{v_{i},t_{\Delta}}\in\mathbb{R}^{H\times
W}\right\\}_{i=1}^{N}$. The proposed GD-CAF is trained in an end-to-end
fashion by minimizing the mean squared error (MSE) between the predicted and
the ground-truth precipitation maps. It’s worth highlighting that when
striving to establish a universal representation across all nodes, shared
filters are utilized for every graph node throughout all convolutional
operations.

$\bullet$ ST-Attention Block:

As illustrated in Fig. 2 (c), the ST-Attention block consists of pooling,
upsampling as well as spatial and temporal attention modules which are
combined through a gated fusion. Furthermore, multi head attention is used to
stabilize the learning process for spatial and temporal attention. Therefore,
the attentions are computed $K$ times with different learnable nonlinear
projections. As opposed to graph attention based models [15, 16, 36] that use
fully-connected layers for computing the queries and keys, here we introduce
convolutional operations as nonlinear projections. This enables us to learn
directly from high-dimensional node representation without the need of
flattening them, therefore the structure of data remains unchanged. In
addition the total number of trainable parameters are also reduced. In
particular, we use depthwise-separable convolution [39] where convolutions are
applied separately on individual input channels and then combined. This
results in a substantial decrease in the number of trainable parameters when
compared to standard convolution, leading to a more efficient model with lower
computational complexity and memory demands.

In situations where the quantity of nodes greatly surpasses the number of
features within each node, numerous approaches have been suggested in the
literature to mitigate computational complexity. These techniques include node
grouping and contextual attention approaches, as discussed in [40, 41, 42]. In
spatiotemporal graphs, where the number of features exceeds the number of
nodes, reducing computational complexity can for instance be achieved through
pooling operations. Pooling involves downsampling or aggregating information
thereby reducing spatial dimension and accelerating overall computation.

$\bullet$ Spatial attention:

(a)

(b)

Figure 3: Spatial and temporal attention on 3D tensors with only one attention
head. (a) Spatial attention is calculated between different nodes, at the same
time step. (b) Temporal attention is calculated within one node, but between
different time steps.

Let us denote the input representation corresponding to the $k$-th head in the
$\ell$-th ST-Attention block at node $v_{i}$ as follows:

$\tilde{x}^{(\ell-1),(k)}_{v_{i}}=[\tilde{x}^{(\ell-1)}_{v_{i},\,t_{{}_{kT+1}}},\tilde{x}^{(\ell-1)}_{v_{i},\,t_{{}_{kT+2}}},\ldots,\tilde{x}^{(\ell-1)}_{v_{i},\,t_{{}_{(k+1)T}}}]\in\mathbb{R}^{T\times
H\times W},$ (2)

and

$\left\\{\begin{array}[]{l}i:1\leq i\leq N\textrm{ is the node index}\\\
k:0\leq k\leq K-1\textrm{ is the head index}\\\ \ell:1\leq\ell\leq L\textrm{
is index of the ST-Attention block.}\end{array}\right.$

Here, $\tilde{x}^{(\ell-1)}_{v_{i},\,t_{{}_{kT+j}}}\in\mathbb{R}^{H\times W}$
is the input representation of the perception map of node $v_{i}$ for the
$k$-th head at time step $t_{kT+j}$ in the $\ell$-th ST-Attention block.

The precipitation in a particular area is influenced by precipitation in other
areas to varying degrees. To model this highly dynamic relationship, we extend
the Scaled Dot-Product Attention [43] to dynamically assign weights to
tensorial node representation, at each time steps, also illustrated in Fig. 3
(a). This relationship between tensorial node $v_{i}$ and $v$ is expressed as
follows:

$s_{v_{i},v}^{(k)}=\frac{\langle
f_{query}^{(k)}*\tilde{x}_{v_{i}}^{(\ell-1),(k)},f_{key}^{(k)}*\tilde{x}_{v}^{(\ell-1),(k)}\rangle}{\sqrt{d}}\in\mathbb{R}^{T}.$
(3)

Here ‘$*$’ represents a double depthwise-separable convolution, consisting of
the following sequential operations: depthwise convolution, pointwise
convolution, normalization and a ReLU activation function.
$\langle\cdot,\cdot\rangle$ is the inner product operator, and $d$ is the size
of a single precipitation map, which is used as a scaling factor. The
relevance $s_{v_{i},v}^{(k)}$ is a $T$ dimensional vector that represents the
relationship between tensorial node $v$ and node $v_{i}$ in the $k$-th head.
The attention scores are computed as follows:

$\alpha_{v_{i},v}^{(k)}=\frac{exp\Big{(}\textrm{LeakyReLu}(s_{v_{i},v}^{(k)})\Big{)}}{\sum_{v_{m}\in
V}exp\Big{(}\textrm{LeakyReLu}(s_{v_{i},v_{m}}^{(k)})\Big{)}}\in\mathbb{R}^{T},$
(4)

where $\alpha_{v_{i},v}^{(k)}$ is the attention score between tensorial node
$v_{i}$ and $v$ in the $k$-th head. By normalizing we ensure that for each
node $\sum_{v\in V}\alpha_{v_{i},v}^{(k)}=1$. The output of the spatial
attention module for node $v_{i}$ in $\ell$-th ST-Attention block is then
computed as follows:

$p_{v_{i}}^{(\ell)}=\bigparallel_{k=1}^{K}\bigg{(}\sum_{v\in
V}\alpha_{v_{i},v}^{(k)}\Big{(}f_{value}^{(k)}*\tilde{x}_{v_{i}}^{(\ell-1),(k)}\Big{)}\bigg{)}\in\mathbb{R}^{T\times
H\times W},$ (5)

where $\parallel$ is the concatenation operation and ‘$*$’, defined as
previously, is applied to the transformed node representation
$\tilde{x}_{v_{i}}^{(\ell-1),(k)}$ with learnable parameters
$f_{value}^{(k)}$. Next, the obtained output of Eq. (5) is fed into a double
convolutional operation to further enhance the feature representation by
capturing more local patterns and structures.

$\bullet$ Temporal attention:

Precipitation in one area is influenced by the past values in the same area.
To model this highly dynamic relationship, we extend the Scaled Dot-Product
Attention [43], also illustrated in Fig. 3 (b). Let us denote the input
representation of all the nodes corresponding to the $k$-th head in the
$\ell$-th ST-Attention block at time step $t_{i}$ as follows:

$z^{(\ell-1),(k)}_{t_{i}}=[\tilde{x}^{(\ell-1),(k)}_{v_{1},\,t_{{}_{i}}},\tilde{x}^{(\ell-1),(k)}_{v_{2},\,t_{{}_{i}}},\ldots,\tilde{x}^{(\ell-1),(k)}_{v_{N},\,t_{{}_{i}}}]\in\mathbb{R}^{N\times
H\times W},$ (6)

where

$\left\\{\begin{array}[]{l}i:1\leq i\leq T\textrm{ is the time step},\\\
k:0\leq k\leq K-1\textrm{ is the head index},\\\ \ell:1\leq\ell\leq L\textrm{
index of the ST-Attention block.}\end{array}\right.$

Here, $\tilde{x}^{(\ell-1),(k)}_{v,\,t_{i}}\in\mathbb{R}^{H\times W}$ is the
input representation of the perception map, corresponding to the $k$-th head,
of node $v$ at time step $t_{i}$ in the $\ell$-th ST-Attention block. The
relevance between all the nodes at time step $t_{i}$, and $t$ in the $k$-th
head can now be expressed as follows:

$u_{t_{i},\,t}^{(k)}=\frac{\langle
g_{query}^{(k)}*z_{t_{i}}^{(\ell-1),(k)}\,,\,g_{key}^{(k)}*z_{t}^{(\ell-1),(k)}\rangle}{\sqrt{d}}\in\mathbb{R}^{N}.$
(7)

Here, the relevance $u_{t_{i},t}^{(k)}$ is an $N$ dimensional vector. The
symbol ‘$*$’ is defined as previously, and it is now applied using learnable
parameters associated with $g_{query}^{(k)}$, and $g_{key}^{(k)}$.
$\langle\cdot,\cdot\rangle$ is defined as previously, and $d$ is the size of a
single precipitation map, serving as a scaling factor. Subsequently, the
attention scores are obtained as follows:

$\beta_{t_{i},\,t}^{(k)}=\frac{exp\Big{(}\textrm{LeakyReLu}(u_{t_{i},\,t}^{(k)})\Big{)}}{\sum_{t_{m}\in
T}exp\Big{(}\textrm{LeakyReLu}(u_{t_{i},\,t_{m}}^{(k)})\Big{)}}\in\mathbb{R}^{N},$
(8)

where the attention score $\beta_{t_{i},t}^{(k)}$ is an $N$ dimensional vector
in the $k$-th head, indicating the significance between time step $t$ and
$t_{i}$. The output of the temporal attention module in the $\ell$-th block is
then computed as follows:

$o_{t_{i}}^{(\ell)}=\bigparallel_{k=1}^{K}\bigg{(}\sum_{t\in
T}\beta_{t_{i},\,t}^{(k)}\Big{(}g_{value}^{(k)}*z_{t_{i}}^{(\ell-1),(k)}\Big{)}\bigg{)}\in\mathbb{R}^{N\times
H\times W}.$ (9)

Here, $o_{t_{i}}^{(\ell)}$ is the output representation of the temporal
attention in the $\ell$-th ST-Attention block, and ‘$*$’ is defined as
previously and applied to $\tilde{z}_{t_{i}}^{(\ell-1),(k)}$ using learnable
parameters associated with $g_{value}^{(k)}$. Next, in order to enrich the
learned representations, the output of Eq. (9) is fed into a double
convolutional operation.

The above equations show how temporal attention in GD-CAF captures the
correlations between different time steps for a given node. The convolution
operation extracts local information from the node tensors, and the attention
mechanism assigns weights to each time step based on their relevance to the
target time step.

$\bullet$ Gated Fusion:

We use gated fusion to combine the obtained representations of spatial and
temporal attention modules. To this end, in the $\ell$-th ST-attention block,
the outputs $P^{(\ell)}$ and $O^{(\ell)}$ are first concatenated followed by a
double depthwise-separable convolution operation to create a new
representation $\tilde{G}^{(\ell)}$ as follows:

$\tilde{G}^{(\ell)}=h_{gated}*\big{(}P^{(\ell)}\parallel O^{(\ell)}\big{)}.$
(10)

Here, ‘$*$’ is defined as previously with learnable parameters $h_{gated}$,
and concatenation is performed along the temporal depth. Gated fusion controls
the flow between spatial and temporal attention at each node and time step.
The output of the gated fusion is upsampled, in case pooling is used. Then the
output of the $(\ell)$-th ST-Attention block is obtained as follows (see Fig.
2(c)):

$\tilde{X}^{(\ell)}=\tilde{X}^{(\ell-1)}+\textrm{upsample}(G^{(\ell)}).$ (11)

###  3.3 Training

Figure 4: Precipitation data captured from the $4$-th region over a 6-hour
time interval.

All models are trained using PyTorch Lightning [44] framework for a maximum of
150 epochs. However, an early stopping criterion that keeps track of the
validation loss is used. The training stops when the validation loss no longer
decreases in the last 15 epochs. This criterion was met in all training
iterations and the maximum of 100 epochs was never reached. Additionally, we
used a learning rate scheduler that reduced the learning rate to a tenth of
the previous learning rate when the validation loss did not improve for four
consecutive epochs. The initial learning rate was set to 0.001 and we used the
Adam optimizer [45] with default values. The training was done on Google Colab
Notebook [46] on a single Nvidia Tesla T4 graphics card that had 16 GB of
memory.

We configured the hyperparameters of our model based on our empirical
observations and also took into account that the compared models should have
similar computational times. As a result, we opted for two ST-Attention blocks
($L=2$), each with four attention heads ($K=4$). For both GD-CAF and SmaAt-
UNet, we used two kernels per channel within the depthwise-separable
convolutional operations. The persistence model uses the most recent available
observation for each region and treats it as the model prediction. An example
of the nowcastings obtained by our proposed GD-CAF model is shown in Fig. 11.

###  3.4 Model evaluation

In order to evaluate the performance of our proposed model as well as the
other examined models, we use the same metrics that are used in [9]. Our main
metric, the loss function used in this study is the mean squared error (MSE)
between the predicted and the ground truth precipitation maps.

$MSE=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\hat{Y}_{i})^{2},$ (12)

where $n$ is the number of samples, $Y_{i}$ is the ground truth for the $i$-th
sample, which contains $N$ nodes, and $\hat{Y}_{i}$ is the predicted value for
the $i$-th sample, also containing $N$ nodes.

Following the lines of [9], in addition to the MSE, we also compute other
metrices such as Precision, Recall (probability of detection), Accuracy and
F1-score, critical success index (CSI), false alarm rate (FAR) and Heidke
Skill Score (HSS). Similar to [9], these scores are calculated for rainfall
bigger than a threshold of $0.5mm/h$. To do this, we convert each pixel of the
predicted output and target images to a boolean mask using this threshold.
From this, one can calculate the true positives (TP) (prediction = 1, target =
1), false positives (FP) (prediction = 1, target = 0), true negatives (TN)
(prediction = 0, target = 0) and false negatives (FN) (prediction = 0, target
= 1). Subsequently, the CSI, FAR and HSS metrics can be computed as follows:

$CSI=\frac{TP}{TP+FP+FN},$ (13)

$FAR=\frac{FP}{TP+FP},$ (14)

$\displaystyle HSS=\frac{2(TP\times TN)-2(FP\times
FN)}{(TP+FN)(FN+TN)+(TP+FP)(FP+TN)}.$ (15)

##  4 Experiments

Figure 5: Randomly placed areas on the map containing local precipitation maps. |  Table 1: 16 selected regions, and their coordinates. LU means the left upper corner of the bounding box. Region id LU* latitude LU* longitude 1 66.58 18.41 2 57.31 2.44 3 51.37 -20.27 4 64.03 -29.25 5 73.19 -23.02 6 68.34 -8.54 7 62.18 54.59 8 44.19 53.59 9 44.19 80.04 10 80.70 17.91 11 77.54 2.19 12 72.44 39.36 13 79.93 -27.26 14 78.39 48.85 15 52.63 26.64 16 77.76 80.79
---|---

###  4.1 Precipitation map dataset

Copernicus Emergency Management Service (CEMS) offers a variety of
meteorological data. These datasets are accessible for various periods on a
global and regional scale. The meteorological information obtained from CEMS,
namely the ERA5 hourly data on single levels [47], is used in this work to
collect precipitation data from the Euro-Asian Region. The longitudinal and
latitudinal boundaries of the selected area are West: $-31^{\circ}$, East:
$100^{\circ}$, South: $15^{\circ}$, and North: $82^{\circ}$.

We have collected hourly precipitation maps over seven years, from January
2016 to December 2022. The first 6 years are used for training, and the last
year is used for testing. When training the models, first we shuffle the
indices of our dataset, then create a training and validation split with a
$0.9/0.1$ ratio. The training set ($47348$ samples) is used for training the
models, and the validation set ($5260$ samples) is used for model selection.
After training the models, we chose the one with the lowest validation loss
for each model type. These top-performing models are then tested on the test
set ($8760$ samples). Our goal is to study precipitation maps of disjoint
regions. Therefore, we have selected sixteen disjoint regions on the map, and
used their precipitation data. Fig. 5 shows the selected sixteen regions
placed on the map. These regions form the nodes of the graph and are referred
to as R1, R2, $\ldots$, R16, and their positions on the map are tabulated in
Table 1.

###  4.2 Studied scenarios

We performed a series of experiments to quantify the nowcasting performance of
the proposed GD-CAF model, comparing it with the SmaAt-UNet and persistence
models. Specifically, we examined three scenarios, which are detailed as
follows:

  * 1.

Ablation study:

To explore the enhanced impact of pooling within GD-CAF, we assessed the
model’s performance by testing the inclusion of pooling layers after the input
data and/or within the ST-Attention blocks. The detailed cases studied are
described in Table 2.

  * 2.

Changing graph size:

To explore the additional benefits of incorporating multiple disjoint regions,
we assessed the models’ performance across varying graph sizes. While the
models produced predictions for all nodes during training, the testing phase
was conducted in two distinct cases: one where Mean Squared Error (MSE) and
other metrics were computed exclusively for the target region R1, and another
where we evaluated all regions collectively. To attain this objective, we
generated five graphs of varying sizes: 1, 2, 4, 8, and 16. In the case of a
graph size of 1, it exclusively comprises the region R1. For a graph size of
2, it encompasses both R1 and R2. Similarly, a graph size of 4 includes
regions R1 through R4. As we progressively increase the graph size to 8, it
spans from R1 to R8. Finally, with a graph size of 16, it encompasses an
extensive range, including regions R1 through R16. Our predictions were
centered on the upcoming 6 hours, leveraging data from the preceding 6 hours.

Figure 6: Changing the number of nodes in the graph, and calculating MSE on all regions. |  Figure 7: Changing the number of nodes in the graph, and calculating MSE only on one target region (R1).
---|---
  * 3.

Changing number of past observations and prediction time:

We investigated the performance of the models under variations in both input
sample size and prediction time. Specifically, we examined four different
input durations, i.e., 6 hours, 9 hours, 12 hours, and 15 hours, and conducted
a series of nowcasting tasks for each input size, predicting outcomes at
1-hour, 3-hour, and 6-hour ahead.

##  5 Results and discussion

###  5.1 Ablation study

Unlike the UNet architecture, which commonly reduces input size through
downsampling, our model has the tendency to increase input size by
incorporating multiple attention heads. This behavior is akin to graph
attention models. To alleviate this computational overhead, we implement
pooling techniques to counteract the size increase and uphold computational
efficiency.

The performance evaluation of different GD-CAF variants is depicted in Table
2. Specifically, GD-CAF (case id=1) exhibits the lowest MSE, while GD-CAF
(case id=4) emerges as the most efficient option. In the upcoming experiments,
we will utilize the latter variant (case id=4) due to its comparable
computational time with SmaAt-UNet, while still maintaining an MSE similar to
other GD-CAF variants. Integrating pooling into either the input layer or the
computation of query ($Q$), key ($K$), and value ($V$) matrices leads to a
slight reduction in performance, albeit still surpassing benchmark models.
However, this integration enhances computational speed by 2.13 ($=1920/900$)
when pooling is applied to $Q$, $K$, $V$, or by 3.33 $(=1920/576)$ when
directly utilizing pooling on the input.

Table 2: An ablation study on the effects of pooling in GD-CAF with regards to
MSE. Computational time is the number of seconds one training epoch took, and
gain represents the percentage increase in MSE compared to the persistence
model.

Model Case id Pooling Q, K and V Pooling Input Training Computational Time
(s/epoch) MSE ↓ Gain ↑ (%) Persistence 0 0.00259167 - SmaAt-UNet 390
0.00151368 1.712 GD-CAF 1 1920 0.00122645 2.113 GD-CAF 2 ✓ 900 0.00122928
2.108 GD-CAF 3 ✓ 576 0.00123289 2.102 GD-CAF 4 ✓ ✓ 312 0.00123959 2.091

Table 3: Test all three models and change the number of nodes in the graph.
However, MSE and other metrics were calculated on target region: R1. A ↑
indicates that higher values for that score are better whereas a ↓ indicates
that lower scores are better.

Nr. nodes Model MSE ↓ Accuracy ↑ Precision ↑ Recall ↑ F1 ↑ CSI ↑ FAR ↓ HSS ↑ 1
Persistence 0.00022440 0.92906 0.23448 0.23437 0.23442 0.13277 0.76552 0.09861
SmaAt-UNet 0.00013329 0.95643 0.15328 0.00021 0.00041 0.00021 0.84672 0.00015
GD-CAF 0.00014201 0.95366 0.88889 0.00004 0.00007 0.00004 0.11111 0.00004 2
Persistence 0.00022440 0.92906 0.23448 0.23437 0.23442 0.13277 0.76552 0.09861
SmaAt-UNet 0.00012864 0.94740 0.38652 0.22994 0.28835 0.16846 0.61348 0.13143
GD-CAF 0.00011324 0.95381 0.50362 0.24009 0.32517 0.19415 0.49638 0.15217 4
Persistence 0.00022440 0.92906 0.23448 0.23437 0.23442 0.13277 0.76552 0.09861
SmaAt-UNet 0.00012223 0.95227 0.45983 0.17133 0.24964 0.14262 0.54017 0.11514
GD-CAF 0.00011282 0.95309 0.48886 0.26972 0.34764 0.21039 0.51114 0.16270 8
Persistence 0.00022440 0.92906 0.23448 0.23437 0.23442 0.13277 0.76552 0.09861
SmaAt-UNet 0.00013153 0.94973 0.40385 0.17796 0.24705 0.14094 0.59615 0.11254
GD-CAF 0.00010970 0.95531 0.53641 0.26345 0.35336 0.21459 0.46359 0.16650 16
Persistence 0.00022440 0.92906 0.23448 0.23437 0.23442 0.13277 0.76552 0.09861
SmaAt-UNet 0.00013031 0.95098 0.39373 0.10682 0.16805 0.09173 0.60627 0.07563
GD-CAF 0.00010874 0.95438 0.51513 0.26706 0.35176 0.21341 0.48487 0.16529

Figure 8: MSE per region for GD-CAF and SmaAt-UNet. Figure 9: MSE of all three models, while changing the past and future window sizes. The performance is averaged on different input window sizes. |  Figure 10: MSE of all three models, while changing the past and future window sizes. The performance is averaged on different prediction steps.
---|---

###  5.2 Changing graph size

The obtained overall test MSE values encompassing all regions are shown in
Fig. 6. It can be seen that as the number of regions increases, the overall
MSE increases, however, the proposed GD-CAF model surpasses the baseline
persistence as well as the SmaAt-UNet model.

Fig. 7 visualizes the obtained test MSE of only target region $R1$, while the
models are trained using different number of nodes. In addition, the obtained
results of other tested metrics are shown in Table 3. When examining a single
target region R1, GD-CAF surpasses SmaAt-UNet when the input graph exceeds a
size of 1. Notably, the Mean Squared Error (MSE) for SmaAt-UNet gradually
decreases as the number of nodes or regions grows, hitting its minimum with 4
input nodes. Yet, post this threshold, the MSE begins to rise. Conversely, for
GD-CAF, the MSE continues to decline with an increasing number of nodes, even
surpassing four regions. This suggests that GD-CAF effectively learns and
leverages the correlations among these distinct regions for enhanced
nowcasting. Meanwhile, the performance of the persistence model remains
consistent, unaffected by changes in graph size.

Additionally, Fig. 8 illustrates the individual improvement achieved for each
region when employing GD-CAF compared to SmaAt-UNet. The most significant
improvements are observed in $R2$, $R3$, $R4$, and $R6$.

###  5.3 Changing input amount and prediction time

The obtained results of the test set, corresponding to various amounts of past
observations utilized by the models to nowcast multiple steps ahead are
tabulated in Table 4. It can be observed that all models outperform the
persistence baseline model. However, it’s noteworthy that within a 1-hour
timeframe, the weather doesn’t undergo significant changes at this resolution.
Hence, persistence effectively yields accurate predictions for such scenarios.

Fig. 9 shows the test MSE averaged over different input size for all examined
models. It can be observed that MSE increases as the future prediction steps
increases. Moreover, GD-CAF consistently achieves lower MSE than SmaAt-UNet.
Fig. 10 visualizes the test MSE averaged over prediction steps for all
examined models. Furthermore Fig. 11 shows examples of predictions from
various regions.

Figure 11: True label, baseline models (Persistence, and SmaAt-UNet), GD-CAF
prediction from six selected regions. Figure 12: The top 20 strongest spatial
attentions are displayed for two time step in a circular graph format. Columns
represent seasons, displaying averaged attention matrices across seasons and
heads, with prominent correlations observed in $R1$, $R4$, $R6$ and $R11$ with
other regions. Figure 13: Averaged temporal attentions are displayed for
specific regions ($R2,R3$), revealing notable correlations between $t=1$ and
$t=2$, and $t=2$ and $t=3$, gradually decreasing afterward. Additionally,
correlations between non-consecutive time steps are also observed. Table 4:
Test all three models and change the input amount and prediction time. We used
all 16 regions during this experiment. A ↑ indicates that higher values for
that score are better whereas a ↓ indicates that lower scores are better.

Input Amount Prediction Model MSE ↓ Accuracy ↑ Precision ↑ Recall ↑ F1 ↑ CSI ↑
FAR ↓ HSS ↑ 6 hour 1 hour Persistence 0.00064499 0.97717 0.71488 0.71486
0.71487 0.55626 0.28512 0.35149 SmaAt-UNet 0.00041348 0.98155 0.79804 0.72163
0.75792 0.61020 0.20196 0.37417 GD-CAF 0.00034831 0.98341 0.80340 0.77519
0.78904 0.65158 0.19660 0.39020 3 hour Persistence 0.00177373 0.95620 0.45288
0.45285 0.45286 0.29271 0.54712 0.21503 SmaAt-UNet 0.00119339 0.96742 0.61723
0.48978 0.54617 0.37568 0.38277 0.26475 GD-CAF 0.00079560 0.97134 0.66401
0.57462 0.61609 0.44518 0.33599 0.30064 6 hour Persistence 0.00259167 0.94233
0.27954 0.27952 0.27953 0.16247 0.72046 0.12474 SmaAt-UNet 0.00151368 0.95901
0.47165 0.20024 0.28112 0.16355 0.52835 0.13178 GD-CAF 0.00123959 0.96280
0.55945 0.33177 0.41653 0.26305 0.44055 0.19930 9 hour 1 hour Persistence
0.00064499 0.97717 0.71488 0.71486 0.71487 0.55626 0.28512 0.35149 SmaAt-UNet
0.00052450 0.98096 0.77444 0.73972 0.75668 0.60860 0.22556 0.37339 GD-CAF
0.00033724 0.98385 0.81423 0.77267 0.79291 0.65688 0.18577 0.39226 3 hour
Persistence 0.00177373 0.95620 0.45288 0.45285 0.45286 0.29271 0.54712 0.21503
SmaAt-UNet 0.00096567 0.96810 0.63735 0.47065 0.54146 0.37123 0.36265 0.26266
GD-CAF 0.00079051 0.97178 0.67401 0.57117 0.61834 0.44754 0.32599 0.30190 6
hour Persistence 0.00259167 0.94233 0.27954 0.27952 0.27953 0.16247 0.72046
0.12474 SmaAt-UNet 0.00150384 0.95922 0.47699 0.19653 0.27836 0.16169 0.52301
0.13055 GD-CAF 0.00123802 0.96257 0.55102 0.34927 0.42754 0.27189 0.44898
0.20460 12 hour 1 hour Persistence 0.00064499 0.97717 0.71488 0.71486 0.71487
0.55626 0.28512 0.35149 SmaAt-UNet 0.00039728 0.98180 0.79094 0.74110 0.76521
0.61971 0.20906 0.37788 GD-CAF 0.00033237 0.98385 0.80460 0.78770 0.79606
0.66121 0.19540 0.39383 3 hour Persistence 0.00177373 0.95620 0.45288 0.45285
0.45286 0.29271 0.54712 0.21503 SmaAt-UNet 0.00110666 0.96587 0.60281 0.43151
0.50298 0.33599 0.39719 0.24290 GD-CAF 0.00078344 0.97176 0.67074 0.57811
0.62099 0.45032 0.32926 0.30320 6 hour Persistence 0.00259167 0.94233 0.27954
0.27952 0.27953 0.16247 0.72046 0.12474 SmaAt-UNet 0.00152703 0.95789 0.44970
0.23396 0.30779 0.18189 0.55030 0.14415 GD-CAF 0.00125883 0.96258 0.55336
0.33652 0.41852 0.26464 0.44664 0.20019 15 hour 1 hour Persistence 0.00064499
0.97717 0.71488 0.71486 0.71487 0.55626 0.28512 0.35149 SmaAt-UNet 0.00043228
0.98127 0.77626 0.74734 0.76153 0.61489 0.22374 0.37589 GD-CAF 0.00032502
0.98410 0.81325 0.78223 0.79744 0.66312 0.18675 0.39458 3 hour Persistence
0.00177373 0.95620 0.45288 0.45285 0.45286 0.29271 0.54712 0.21503 SmaAt-UNet
0.00100494 0.96761 0.63306 0.45344 0.52840 0.35907 0.36694 0.25605 GD-CAF
0.00079395 0.97193 0.67516 0.57532 0.62126 0.45060 0.32484 0.30339 6 hour
Persistence 0.00259167 0.94233 0.27954 0.27952 0.27953 0.16247 0.72046 0.12474
SmaAt-UNet 0.00153682 0.95880 0.46486 0.19871 0.27841 0.16172 0.53514 0.13035
GD-CAF 0.00125927 0.96270 0.55536 0.33808 0.42030 0.26606 0.44464 0.20110

Fig. 12 displays spatial attentions in the last ST-Attention block in a
circular graph format. Each graph shows the top 20 strongest attention values
between different regions. The first row contains spatial attention matrices
from the first time step, while the second row contains those from the last
time step. Each column represents a different season, with the spatial
attention matrices averaged for each season and across attention heads.
Notably, when $t=1$, $R1$ exhibited a strong correlation with $R4$, $R10$, and
$R14$. $R11$ showed a lot of strong connections in spring, autumn, and winter.
When $t=6$, we noticed one region in each season that has high cardinality;
these are $R14$, $R4$, $R13$, and $R6$. Overall, a strong correlation does not
necessarily imply that the cells are physically close to each other. However,
in the cases of $R1-R15$, $R4-R5$, and $R4-R6$, this is indeed the case.

Fig. 13 presents the temporal attentions for selected regions ($R2,R3$). These
temporal attention values were averaged across each season in the test set.
Notably, a strong correlation is observed between $t=1$ and $t=2$, as well as
between $t=2$ and $t=3$. However, after these time steps, the correlation
gradually decreases. It’s also worth mentioning that there is some correlation
between non-consecutive time steps.

##  6 Conclusion

In this study, we introduce an innovative Graph Dual-stream Convolutional
Attention Fusion (GD-CAF) model tailored for precipitation nowcasting tasks.
The incorporation of spatiotemporal convolutional attention and gated fusion
modules, enhanced by depthwise-separable convolutional layers, enables our
model to extract valuable insights from high-dimensional spatiotemporal graph
nodes that encapsulate historical precipitation maps. Leveraging direct
exploration of higher-order correlations among input data dimensions, our
approach proves particularly advantageous in scenarios constrained by limited
access to local information across multiple regions. Rigorous evaluations
against benchmark models such as SmaAt-UNet and persistence models attest to
the superior performance of our GD-CAF model under varying conditions.
Notably, our model surpasses these benchmarks, showcasing its efficacy in
handling the complexities of nowcasting tasks. To offer additional insights
into the predictive capabilities of our model, we visualize the strongest
connections between regions or time periods. This visualization is derived
from the averaged spatial and temporal attention scores computed across each
season within test set. The implementation of our proposed model, including
the trained models, is available on GitHub222https://github.com/wendig/GD-CAF.

## References

  * [1] J. W. Wilson, Y. Feng, M. Chen, R. D. Roberts, Nowcasting challenges during the beijing olympics: Successes, failures, and implications for future nowcasting systems, Weather and Forecasting 25 (6) (2010) 1691–1714.
  * [2] A. Heye, K. Venkatesan, J. Cain, Precipitation nowcasting: Leveraging deep recurrent convolutional neural networks, in: Proceedings of the Cray User Group (CUG), Vol. 2017, Cray User Group Washington, 2017.
  * [3] G. Ayzel, M. Heistermann, T. Winterrath, Optical flow models as an open benchmark for radar-based precipitation nowcasting (rainymotion v0. 1), Geoscientific Model Development 12 (4) (2019) 1387–1402.
  * [4] A. Kumar, T. Islam, Y. Sekimoto, C. Mattmann, B. Wilson, Convcast: An embedded convolutional lstm based architecture for precipitation nowcasting using satellite data, Plos one 15 (3) (2020) e0230114.
  * [5] S. Mehrkanoon, Deep shared representation learning for weather elements forecasting, Knowledge-Based Systems 179 (2019) 120–128.
  * [6] P. Hess, N. Boers, Deep learning for improving numerical weather prediction of heavy rainfall, Journal of Advances in Modeling Earth Systems 14 (3) (2022) e2021MS002765.
  * [7] Y. Han, L. Mi, L. Shen, C. Cai, Y. Liu, K. Li, G. Xu, A short-term wind speed prediction method utilizing novel hybrid deep learning algorithms to correct numerical weather forecasting, Applied Energy 312 (2022) 118777.
  * [8] D. Cho, C. Yoo, J. Im, D.-H. Cha, Comparative assessment of various machine learning-based bias correction methods for numerical weather prediction model forecasts of extreme air temperatures in urban areas, Earth and Space Science 7 (4) (2020) e2019EA000740.
  * [9] K. Trebing, T. Stanczyk, S. Mehrkanoon, SmaAT-UNet: Precipitation nowcasting using a small attention-unet architecture, Pattern Recognition Letters 145 (2021) 178–186.
  * [10] J. G. Fernández, S. Mehrkanoon, Broad-UNet: Multi-scale feature learning for nowcasting tasks, Neural Networks 144 (2021) 419–427.
  * [11] Y. Yang, S. Mehrkanoon, AA-TransUNet: Attention augmented transunet for nowcasting tasks, in: IEEE International Joint Conference on Neural Networks (IEEE-IJCNN), 2022, pp. 01–08.
  * [12] J. G. Fernández, I. A. Abdellaoui, S. Mehrkanoon, Deep coastal sea elements forecasting using unet-based models, Knowledge-Based Systems 252 (2022) 109445.
  * [13] R. Fukuoka, H. Suzuki, T. Kitajima, A. Kuwahara, T. Yasuno, Wind speed prediction model using lstm and 1d-cnn, Journal of Signal Processing 22 (4) (2018) 207–210.
  * [14] K. Trebing, S. Mehrkanoon, Wind speed prediction using multidimensional convolutional neural networks, in: IEEE Symposium Series on Computational Intelligence (IEEE-SSCI), 2020, pp. 713–720.
  * [15] C. Zheng, X. Fan, C. Wang, J. Qi, Gman: A graph multi-attention network for traffic prediction, in: Proceedings of the AAAI conference on artificial intelligence, Vol. 34, 2020, pp. 1234–1241.
  * [16] P. Velickovic, G. Cucurull, A. Casanova, A. Romero, P. Lio, Y. Bengio, et al., Graph attention networks, stat 1050 (20) (2017) 10–48550.
  * [17] M. Aiolfi, A. Timmermann, Persistence in forecasting performance and conditional combination strategies, Journal of Econometrics 135 (1-2) (2006) 31–53.
  * [18] Q. Cao, B. T. Ewing, M. A. Thompson, Forecasting wind speed with recurrent neural networks, European Journal of Operational Research 221 (1) (2012) 148–154.
  * [19] M. A. Zaytar, C. El Amrani, Sequence to sequence weather forecasting with long short-term memory recurrent neural networks, International Journal of Computer Applications 143 (11) (2016) 7–11.
  * [20] X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W.-K. Wong, W.-c. Woo, Convolutional lstm network: A machine learning approach for precipitation nowcasting, Advances in neural information processing systems 28 (2015).
  * [21] P. R. Larraondo, L. J. Renzullo, I. Inza, J. A. Lozano, A data-driven approach to precipitation parameterizations using convolutional encoder-decoder neural networks, arXiv preprint arXiv:1903.10274 (2019).
  * [22] S. Agrawal, L. Barrington, C. Bromberg, J. Burge, C. Gazen, J. Hickey, Machine learning for precipitation nowcasting from radar images, arXiv preprint arXiv:1912.12132 (2019).
  * [23] J. Simeunović, B. Schubnel, P.-J. Alet, R. E. Carrillo, Spatio-temporal graph neural networks for multi-site pv power forecasting, IEEE Transactions on Sustainable Energy 13 (2) (2021) 1210–1220.
  * [24] B. Klein, L. Wolf, Y. Afek, A dynamic convolutional layer for short range weather prediction, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 4840–4848.
  * [25] S. Guastavino, M. Piana, M. Tizzi, F. Cassola, A. Iengo, D. Sacchetti, E. Solazzo, F. Benvenuto, Prediction of severe thunderstorm events with ensemble deep learning and radar data, Scientific Reports 12 (1) (2022) 20049.
  * [26] V. Lebedev, V. Ivashkin, I. Rudenko, A. Ganshin, A. Molchanov, S. Ovcharenko, R. Grokhovetskiy, I. Bushmarinov, D. Solomentsev, Precipitation nowcasting with satellite imagery, in: Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, 2019, pp. 2680–2688.
  * [27] M. Gori, G. Monfardini, F. Scarselli, A new model for learning in graph domains, in: Proceedings. 2005 IEEE international joint conference on neural networks, Vol. 2, 2005, pp. 729–734.
  * [28] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, M. Sun, Graph neural networks: A review of methods and applications, AI open 1 (2020) 57–81.
  * [29] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, S. Y. Philip, A comprehensive survey on graph neural networks, IEEE transactions on neural networks and learning systems 32 (1) (2020) 4–24.
  * [30] T. Stanczyk, S. Mehrkanoon, Deep graph convolutional networks for wind speed prediction, in: ESANN 2021 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning, 2021, pp. 147–152.
  * [31] D. Aykas, S. Mehrkanoon, Multistream graph attention networks for wind speed forecasting, in: IEEE Symposium Series on Computational Intelligence (SSCI), IEEE, 2021, pp. 1–8.
  * [32] R. Keisler, Forecasting global weather with graph neural networks, arXiv preprint arXiv:2202.07575 (2022).
  * [33] D. Owerko, F. Gama, A. Ribeiro, Predicting power outages using graph neural networks, in: 2018 ieee global conference on signal and information processing (globalsip), IEEE, 2018, pp. 743–747.
  * [34] T. N. Kipf, M. Welling, Semi-supervised classification with graph convolutional networks, arXiv preprint arXiv:1609.02907 (2016).
  * [35] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, Y. Bengio, Graph attention networks, arXiv preprint arXiv:1710.10903 (2017).
  * [36] X. Wang, H. Ji, C. Shi, B. Wang, Y. Ye, P. Cui, P. S. Yu, Heterogeneous graph attention network, in: The world wide web conference, 2019, pp. 2022–2032.
  * [37] L. Mi, Z. Chen, Hierarchical graph attention network for visual relationship detection, in: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2020, pp. 13886–13895.
  * [38] C. Zhang, J. James, Y. Liu, Spatial-temporal graph attention networks: A deep learning approach for traffic forecasting, IEEE Access 7 (2019) 166246–166256.
  * [39] F. Chollet, Xception: Deep learning with depthwise separable convolutions, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 1251–1258.
  * [40] N. Kitaev, Ł. Kaiser, A. Levskaya, Reformer: The efficient transformer, arXiv preprint arXiv:2001.04451 (2020).
  * [41] S. Wang, B. Z. Li, M. Khabsa, H. Fang, H. Ma, Linformer: Self-attention with linear complexity, arXiv preprint arXiv:2006.04768 (2020).
  * [42] Z. Shen, M. Zhang, H. Zhao, S. Yi, H. Li, Efficient attention: Attention with linear complexities, in: Proceedings of the IEEE/CVF winter conference on applications of computer vision, 2021, pp. 3531–3539.
  * [43] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin, Attention is all you need, Advances in neural information processing systems 30 (2017).
  * [44] W. Falcon, The PyTorch Lightning team, PyTorch Lightning (Mar. 2019).  doi:10.5281/zenodo.3828935.
URL https://github.com/Lightning-AI/lightning

  * [45] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [46] E. Bisong, E. Bisong, Google colaboratory, Building machine learning and deep learning models on google cloud platform: a comprehensive guide for beginners (2019) 59–64.
  * [47] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, et al., The era5 global reanalysis, Quarterly Journal of the Royal Meteorological Society 146 (730) (2020) 1999–2049.