[a]Johannes Blümlein

#  DESY 19–061, DO-TH 20/09, TTP22–048, RISC Report Series 22–09, SAGEX 20–20,
PoS (LL2022) 048
The 3-loop anomalous dimensions from off-shell operator matrix elements111
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska–Curie grant
agreement No. 764850, SAGEX and from the Austrian Science Fund (FWF) grants
SFB F50 (F5009-N15) and P33530.

Peter Marquard  Carsten Schneider  Kay Schönwald

###### Abstract

We report on the calculation of the three–loop polarized and unpolarized
flavor non–singlet and the polarized singlet anomalous dimensions using
massless off–shell operator matrix elements in a gauge–variant framework. We
also reconsider the unpolarized two–loop singlet anomalous dimensions and
correct errors in the foregoing literature.

##  1 Introduction

In these proceedings we report on recent results of the calculation of
unpolarized and polarized three–loop anomalous dimensions based on massless
off–shell operator matrix elements (OMEs) [1, 2] and the unpolarized two–loop
case [3] correcting previous results in Refs. [4, 5].

The unpolarized and polarized non–singlet and singlet anomalous dimensions
have been calculated at one– [6], two– [7, 3, 8], and three–loop order [9, 10,
12, 11, 13, 14, 1, 2, 15]. Here different techniques as off–shell massless
OMEs, the forward Compton amplitude, in part also with scalar and
gravitational currents, and massive on–shell OMEs have been used. In the
latter case one obtains the contributions $\propto T_{F}$, which are the
complete anomalous dimensions in the cases $(\Delta)\gamma_{qq}^{(2),\rm PS}$
and $(\Delta)\gamma_{qg}^{(2)}$.

The method of massless off–shell OMEs is the traditional way to calculate the
anomalous dimensions, cf. [1a-c]. However, it is a gauge–dependent environment
implying in the unpolarized singlet case new operator mixings [8, 16, 3]. It
is our goal to set up a program chain allowing for a fully automated
calculation of the three–loop anomalous dimensions without making structural
assumptions motivated by QCD or the expected representation of the final
results in terms of harmonic sums [17, 18]. The paper is organized as follows.
We discuss first the basic formalism and describe then the calculation method,
before we present some examples for the three–loop anomalous dimensions.

##  2 Basic Formalism

We form the expectation values of the local twist–2 operators in the
unpolarized and polarized case, cf. e.g. [19], between off–shell quark and
gluon states and consider the physical projections, which are calculated to
three–loop order. There are also other contributions due to the violation of
the equation of motion (EOM), which, however, are not related to the anomalous
dimensions. From the physical projection we extract the three–loop anomalous
dimensions from the $O(1/\varepsilon)$ terms, where $\varepsilon=D-4$ denotes
the dimensional parameter. The corresponding amplitudes are gauge–variant,
i.e. they depend on the gauge parameter $\xi$ in the $R_{\xi}$ gauges. The
structure of the pole terms are fully predicted by the renormalization group
and do depend on lower order expansion coefficients to higher orders in
$\varepsilon$. In the unpolarized case, the mixing of more local operators has
to be considered [8, 16, 4, 5, 3], which we studied up to two–loop order in
[3].

In the polarized singlet case the anomalous dimensions are first calculated in
the Larin scheme [20], which is a consistent scheme. Finally, we transform the
anomalous dimensions to the M–scheme [5, 11].

##  3 Details of the calculation

The calculation of the unrenormalized massless off–shell OMEs is performed in
the following way. Diagram generation, the performance of Lorentz/Dirac and
color algebra are performed by using the packages QGRAF, FORM and color [21,
22, 23, 24]. The local operators are resummed into propagators by observing
the current crossing relations, cf. [26, 25], as has been described in Ref.
[27],

$\displaystyle\sum_{N=0}^{\infty}(\Delta.k)^{N}\left(t^{N}\pm(-t)^{N}\right)\rightarrow\left[\frac{1}{1-\Delta.k~{}t}\pm\frac{1}{1+\Delta.k~{}t}\right].$
(1)

The three–loop anomalous dimensions are obtained from the contributions of
$O(1/\varepsilon)$, determining the $N$th moment analytically.

In the individual channels there are up to $O(1600)$ Feynman diagrams. They
are reduced to up to $O(250)$ master integrals using the code Crusher [28] by
applying the integration–by–parts relations [29, 30]. Coupling constant and
wave function renormalizations are performed [31, 32] and results from lower
order factorizing diagrams [3] are accounted for. The method of arbitrary high
Mellin moments [33], implemented in the package SolveCoupledSystem [34], is
used to provide a necessary input set, here of 3000 moments. The necessary
initial values for the difference equations can be obtained from [30, 35]. We
then used the method of guessing [36, 37] and its implementation in Sage [38,
39] to obtain the recurrences for the different color and multiple zeta value
factors [40]. It turned out that $O(1600)$ moments suffice and the largest
difference equations was of order o = 16 and degree d = 304. These difference
equations are solved using methods from difference ring theory [41]
implemented in the package Sigma [42, 43]. Functions of the package
HarmonicSums [44, 47, 17, 18, 45, 46] are used to compactify the final
results. The automated calculation at Intel(R) Xeon(R) CPU E5-2643 v4
processors amounted to about 40 days of CPU time for the projects [1, 2].

All anomalous dimensions can be expressed in terms of harmonic sums [17, 18]

$\displaystyle S_{b,\vec{a}}(N)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{N}\frac{({\rm
sign}(b))^{k}}{k^{|b|}}S_{\vec{a}}(k),~{}~{}~{}S_{\emptyset}=1,~{}~{}~{}b_{,}a_{i}\in\mathbb{Z}\backslash\\{0\\},N\in\mathbb{N}\backslash\\{0\\}.$
(2)

Accordingly, the corresponding $z$-space expressions are given by harmonic
polylogarithms ${\rm H}_{\vec{a}}(z)$, [45]. Here $z\in[0,1]$ denotes the
momentum fraction w.r.t. the incoming nucleon momentum in the deep–inelastic
process.

##  4 Anomalous Dimensions

We have calculated all the non–singlet anomalous dimensions
$\gamma_{qq}^{\pm,\rm NS},\gamma_{qq}^{s,\rm NS},\Delta\gamma_{qq}^{s,\rm NS}$
for unpolarized and polarized deep-inelastic scattering and also for
transversity $\gamma_{qq}^{\pm,\rm tr,NS}$, as well as the singlet polarized
anomalous dimensions, Ref. [1, 2]. They can all be expressed in terms of
harmonic sums, for which we apply the algebraic relations [46]. If also the
structural relations are applied [47] only 10 harmonic sums contribute.

$\displaystyle\\{S_{1},S_{2,1},S_{-2,1},S_{-3,1},S_{-4,1},S_{2,1,1},S_{-2,1,1},S_{2,1,-2},S_{-3,1,1},S_{-2,1,1,1}\\}.$
(3)

As examples we show one of the transversity anomalous dimensions

$\displaystyle\gamma_{\rm NS}^{(2),\rm
tr,+}=\frac{1}{2}\left[1+(-1)^{N}\right]$
$\displaystyle\times\Biggl{\\{}{\color[rgb]{0,0,1}C_{F}}\Biggl{\\{}{\color[rgb]{0,0,1}T_{F}^{2}N_{F}^{2}}\Biggl{[}\frac{8\big{(}-8+17N+17N^{2}\big{)}}{9N(1+N)}-\frac{128}{27}S_{1}-\frac{640}{27}S_{2}+\frac{128}{9}S_{3}\Biggr{]}$
$\displaystyle+{\color[rgb]{0,0,1}C_{A}T_{F}N_{F}}\Biggl{[}-\frac{16\big{(}-22+45N+45N^{2}\big{)}}{9N(1+N)}+\Biggl{(}-\frac{16\big{(}9+209N+209N^{2}\big{)}}{27N(1+N)}+64S_{3}$
$\displaystyle+\frac{256}{3}S_{-2,1}-128\zeta_{3}\Biggr{)}S_{1}+\frac{5344}{27}S_{2}-\frac{448}{3}S_{3}+\frac{320}{3}S_{4}+\Biggl{(}-\frac{1280}{9}S_{1}+\frac{128}{3}S_{2}\Biggr{)}S_{-2}$
$\displaystyle+\Biggl{(}-\frac{640}{9}+\frac{128}{3}S_{1}\Biggr{)}S_{-3}+\frac{128}{3}S_{-4}-\frac{256}{3}S_{3,1}+\frac{1280}{9}S_{-2,1}+\frac{128}{3}S_{-2,2}-\frac{512}{3}S_{-2,1,1}$
$\displaystyle+96\zeta_{3}\Biggr{]}+{\color[rgb]{0,0,1}C_{A}^{2}}\Biggl{[}\frac{-968+1657N+1657N^{2}}{18N(1+N)}+\Biggl{(}\frac{4P_{14}}{3(-1+N)N(1+N)(2+N)}-176S_{3}$
$\displaystyle-256S_{4}+512S_{3,1}-\frac{704}{3}S_{-2,1}-1024S_{-2,2}-1024S_{-3,1}+2048S_{-2,1,1}\Biggr{)}S_{1}+\Biggl{(}-128S_{3}$
$\displaystyle-512S_{-2,1}\Biggr{)}S_{1}^{2}+\Biggl{(}-\frac{8344}{27}+384S_{3}+1536S_{-2,1}\Biggr{)}S_{2}+\frac{3112}{9}S_{3}-\frac{880}{3}S_{4}+64S_{5}$
$\displaystyle+\Biggl{(}\frac{16P_{2}}{(-1+N)N(1+N)(2+N)}+\frac{32\big{(}-241+134N+134N^{2}\big{)}S_{1}}{9(-1+N)(2+N)}-\frac{352}{3}S_{2}-64S_{3}$
$\displaystyle-1536S_{2,1}+128S_{-2,1}-192\zeta_{3}\Biggr{)}S_{-2}+\Biggl{(}48-192S_{1}\Biggr{)}S_{-2}^{2}+\Biggl{(}256S_{1}^{2}-768S_{2}-320S_{-2}$
$\displaystyle+\frac{32\big{(}-107+67N+67N^{2}\big{)}}{9(-1+N)(2+N)}-\frac{352}{3}S_{1}\Biggr{)}S_{-3}+\Biggl{(}-\frac{208}{3}+320S_{1}\Biggr{)}S_{-4}-704S_{-5}-384S_{2,3}$
$\displaystyle-768S_{2,-3}+\frac{704}{3}S_{3,1}+384S_{4,1}-\frac{64\big{(}-107+67N+67N^{2}\big{)}S_{-2,1}}{9(-1+N)(2+N)}-\frac{352}{3}S_{-2,2}$
$\displaystyle+1088S_{-2,3}-448S_{-4,1}+1536[S_{2,1,-2}+S_{-2,2,1}+S_{-3,1,1}]-768S_{3,1,1}$
$\displaystyle+\frac{1408}{3}S_{-2,1,1}+512S_{-2,1,-2}-3072S_{-2,1,1,1}-\frac{24\big{(}-6+5N+5N^{2}\big{)}\zeta_{3}}{(-1+N)(2+N)}\Biggr{]}\Biggr{\\}}$
$\displaystyle+{\color[rgb]{0,0,1}C_{F}^{2}}\Biggl{\\{}{\color[rgb]{0,0,1}T_{F}N_{F}}\Biggl{[}92+\Biggl{(}-\frac{8\big{(}-8+55N+55N^{2}\big{)}}{3N(1+N)}+\frac{1280}{9}S_{2}-\frac{512}{3}S_{3}-\frac{512}{3}S_{-2,1}+128\zeta_{3}\Biggr{)}$
$\displaystyle\times
S_{1}-\frac{80}{3}S_{2}-\frac{128}{3}S_{2}^{2}+\frac{1856}{9}S_{3}-\frac{512}{3}S_{4}+\Biggl{(}\frac{2560}{9}S_{1}-\frac{256}{3}S_{2}\Biggr{)}S_{-2}+\Biggl{(}\frac{1280}{9}$
$\displaystyle-\frac{256}{3}S_{1}\Biggr{)}S_{-3}-\frac{256}{3}S_{-4}+\frac{256}{3}S_{3,1}-\frac{2560}{9}S_{-2,1}-\frac{256}{3}S_{-2,2}+\frac{1024}{3}S_{-2,1,1}-96\zeta_{3}\Biggr{]}$
$\displaystyle+{\color[rgb]{0,0,1}C_{A}}\Biggl{[}-\frac{151}{2}+\Biggl{(}-\frac{8\big{(}-206+211N+211N^{2}\big{)}}{3(-1+N)N(1+N)(2+N)}-\frac{4288}{9}S_{2}+\frac{1984}{3}S_{3}+320S_{4}-1024S_{3,1}$
$\displaystyle+\frac{1984}{3}S_{-2,1}+3712S_{-2,2}+3840S_{-3,1}-7168S_{-2,1,1}\Biggr{)}S_{1}+\Biggl{(}256S_{3}+1792S_{-2,1}\Biggr{)}S_{1}^{2}$
$\displaystyle+\Biggl{(}\frac{604}{3}-832S_{3}-5248S_{-2,1}\Biggr{)}S_{2}+\frac{352}{3}S_{2}^{2}-\frac{6160}{9}S_{3}+\frac{2416}{3}S_{4}$
$\displaystyle+\Biggl{(}-\frac{48P_{2}}{(-1+N)N(1+N)(2+N)}+\Biggl{(}-\frac{64P_{7}}{9(-1+N)N(1+N)(2+N)}-256S_{2}\Biggr{)}S_{1}$
$\displaystyle+\frac{992}{3}S_{2}+64S_{3}+5376S_{2,1}-384S_{-2,1}+576\zeta_{3}\Biggr{)}S_{-2}+\Biggl{(}-96+512S_{1}\Biggr{)}S_{-2}^{2}$
$\displaystyle+\Biggl{(}-\frac{32\big{(}-187+134N+134N^{2}\big{)}}{9(-1+N)(2+N)}+\frac{992}{3}S_{1}-1152S_{1}^{2}+2624S_{2}+960S_{-2}\Biggr{)}S_{-3}$
$\displaystyle+\Biggl{(}\frac{560}{3}-1472S_{1}\Biggr{)}S_{-4}+2304S_{-5}+768S_{2,3}+2688S_{2,-3}-\frac{1856}{3}S_{3,1}-768S_{4,1}$
$\displaystyle+\frac{64\big{(}-187+134N+134N^{2}\big{)}S_{-2,1}}{9(-1+N)(2+N)}+\frac{992}{3}S_{-2,2}-3648S_{-2,3}+1728S_{-4,1}$
$\displaystyle-5376[S_{2,1,-2}+S_{-2,2,1}+S_{-3,1,1}]+1536[S_{3,1,1}-S_{-2,1,-2}]-\frac{3968}{3}S_{-2,1,1}$
$\displaystyle+10752S_{-2,1,1,1}+\frac{72\big{(}-6+5N+5N^{2}\big{)}\zeta_{3}}{(-1+N)(2+N)}\Biggr{]}\Biggr{\\}}$
$\displaystyle+{\color[rgb]{0,0,1}C_{F}^{3}}\Biggl{\\{}-29+\Biggl{(}\frac{384\big{(}-1+N+N^{2}\big{)}}{(-1+N)N(1+N)(2+N)}+128S_{2}^{2}-384S_{3}+128S_{4}+512S_{3,1}$
$\displaystyle-384S_{-2,1}-3328S_{-2,2}-3584S_{-3,1}+6144S_{-2,1,1}\Biggr{)}S_{1}-256S_{-2}^{2}S_{1}+\Biggl{(}12+512S_{3}$
$\displaystyle+4352S_{-2,1}\Biggr{)}S_{2}-96S_{2}^{2}+104S_{3}-480S_{4}+\Biggl{(}\frac{32P_{2}}{(-1+N)N(1+N)(2+N)}$
$\displaystyle+\Biggl{(}\frac{384}{(-1+N)N(1+N)(2+N)}+512S_{2}\Biggr{)}S_{1}-192S_{2}+128S_{3}-4608S_{2,1}+256S_{-2,1}$
$\displaystyle-384\zeta_{3}\Biggr{)}S_{-2}+\Biggl{(}\frac{192}{(-1+N)(2+N)}-192S_{1}+1280S_{1}^{2}-2176S_{2}-640S_{-2}\Biggr{)}S_{-3}$
$\displaystyle+\Biggl{(}-96+1664S_{1}\Biggr{)}S_{-4}-1792S_{-5}+384[-S_{2,3}+S_{3,1}+S_{4,1}]-2304S_{2,-3}$
$\displaystyle-\frac{384S_{-2,1}}{(-1+N)(2+N)}-1536S_{1}^{2}S_{-2,1}-192S_{-2,2}+2944S_{-2,3}-1664S_{-4,1}+4608S_{2,1,-2}$
$\displaystyle-768S_{3,1,1}+768S_{-2,1,1}+1024S_{-2,1,-2}+4608[S_{-2,2,1}+S_{-3,1,1}]-9216S_{-2,1,1,1}$
$\displaystyle-\frac{48\big{(}-6+5N+5N^{2}\big{)}\zeta_{3}}{(-1+N)(2+N)}\Biggr{\\}}\Biggr{\\}}$
(4)

and one of the polarized singlet anomalous dimensions,

$\displaystyle\Delta\gamma_{gg}^{(2)}$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,1}C_{A}T_{F}^{2}N_{F}^{2}}\Biggl{[}-\frac{16P_{8}}{27N^{2}(1+N)^{2}}S_{1}-\frac{4P_{48}}{27N^{3}(1+N)^{3}}\Biggr{]}+{\color[rgb]{0,0,1}C_{F}}\Biggl{[}{\color[rgb]{0,0,1}T_{F}^{2}N_{F}^{2}}\Biggl{[}-\frac{8P_{59}}{27N^{4}(1+N)^{4}}$
(5)
$\displaystyle+\frac{64(N-1)(2+N)\big{(}-6-8N+N^{2}\big{)}}{9N^{3}(1+N)^{3}}S_{1}+\frac{32(N-1)(2+N)}{3N^{2}(1+N)^{2}}S_{1}^{2}$
$\displaystyle-\frac{32(N-1)(2+N)}{N^{2}(1+N)^{2}}S_{2}\Biggr{]}+{\color[rgb]{0,0,1}C_{A}T_{F}N_{F}}\Biggl{[}\frac{8P_{6}}{N^{3}(1+N)^{3}}S_{2}-\frac{8P_{9}}{3N^{3}(1+N)^{3}}S_{1}^{2}$
$\displaystyle+\frac{2P_{77}}{27(N-1)N^{5}(1+N)^{5}(2+N)}+\Biggl{(}-\frac{8P_{67}}{9(-1+N)N^{4}(1+N)^{4}(2+N)}$
$\displaystyle-\frac{32(N-1)(2+N)}{N^{2}(1+N)^{2}}S_{2}+128\zeta_{3}\Biggr{)}S_{1}+\frac{32(N-1)(2+N)}{3N^{2}(1+N)^{2}}S_{1}^{3}-\frac{32\big{(}34+N+N^{2}\big{)}}{3N^{2}(1+N)^{2}}$
$\displaystyle\times
S_{3}+\Biggl{(}\frac{128P_{2}}{(N-1)N^{2}(1+N)^{2}(2+N)}S_{1}-\frac{32P_{23}}{(N-1)N^{2}(1+N)^{3}(2+N)}\Biggr{)}S_{-2}$
$\displaystyle-\frac{192\big{(}4-N-N^{2}\big{)}}{N^{2}(1+N)^{2}}S_{-3}+\frac{64(N-1)(2+N)}{N^{2}(1+N)^{2}}S_{2,1}-\frac{128\big{(}-8+N+N^{2}\big{)}}{N^{2}(1+N)^{2}}S_{-2,1}$
$\displaystyle-\frac{64(-3+N)(4+N)}{N^{2}(1+N)^{2}}\zeta_{3}\Biggr{]}\Biggr{]}+{\color[rgb]{0,0,1}C_{A}^{3}}\Biggl{[}\frac{64P_{16}}{9N^{2}(1+N)^{2}}S_{-2,1}-\frac{32P_{18}}{9N^{2}(1+N)^{2}}S_{3}$
$\displaystyle+\frac{P_{74}}{27(N-1)N^{5}(1+N)^{5}(2+N)}+\Biggl{(}\frac{4P_{69}}{9(N-1)N^{4}(1+N)^{4}(2+N)}$
$\displaystyle-\frac{64P_{17}}{9N^{2}(1+N)^{2}}S_{2}+128S_{2}^{2}+\frac{16\big{(}-96+11N+11N^{2}\big{)}}{3N(1+N)}S_{3}+192S_{4}$
$\displaystyle+\frac{1024}{N(1+N)}S_{-2,1}-640S_{-2,2}-768S_{-3,1}+1024S_{-2,1,1}\Biggr{)}S_{1}$
$\displaystyle+\Biggl{(}-\frac{256\big{(}1+3N+3N^{2}\big{)}}{N^{3}(1+N)^{3}}+128S_{3}-256S_{-2,1}\Biggr{)}S_{1}^{2}+\Biggl{(}-\frac{16P_{41}}{9N^{3}(1+N)^{3}}$
$\displaystyle+64S_{3}+640S_{-2,1}\Biggr{)}S_{2}-\frac{256}{N(1+N)}S_{2}^{2}-\frac{384}{N(1+N)}S_{4}+64S_{5}$
$\displaystyle+\Biggl{(}\frac{32P_{52}}{9(N-1)N^{3}(1+N)^{3}(2+N)}+\big{(}-\frac{64P_{32}}{9(-1+N)N(1+N)^{2}(2+N)}+256S_{2}\Biggr{)}$
$\displaystyle\times
S_{1}-\frac{512}{N(1+N)}S_{2}+128S_{3}-768S_{2,1}\Biggr{)}S_{-2}+\Biggl{(}-\frac{16\big{(}24+11N+11N^{2}\big{)}}{3N(1+N)}$
$\displaystyle+64S_{1}\Biggr{)}S_{-2}^{2}+\Biggl{(}-\frac{32P_{15}}{9N^{2}(1+N)^{2}}-\frac{1536}{N(1+N)}S_{1}+384S_{1}^{2}-320S_{2}\Biggr{)}S_{-3}$
$\displaystyle+\Biggl{(}-\frac{1024}{N(1+N)}+512S_{1}\Biggr{)}S_{-4}-192S_{-5}-384S_{2,-3}+\frac{1280}{N(1+N)}S_{-2,2}$
$\displaystyle+384S_{-2,3}+\frac{1536}{N(1+N)}S_{-3,1}-384S_{-4,1}+768S_{2,1,-2}-\frac{2048}{N(1+N)}S_{-2,1,1}$
$\displaystyle+768[S_{-2,2,1}+S_{-3,1,1}]-1536S_{-2,1,1,1}\Biggr{]}$
$\displaystyle+{\color[rgb]{0,0,1}C_{F}^{2}T_{F}N_{F}}\Biggl{[}-\frac{4P_{75}}{(N-1)N^{5}(1+N)^{5}(2+N)}+\Biggl{(}\frac{32(N-1)(2+N)S_{2}}{N^{2}(1+N)^{2}}$
$\displaystyle-\frac{16P_{42}}{N^{4}(1+N)^{4}}\Biggr{)}S_{1}+\frac{8(N-1)(2+N)\big{(}2+3N+3N^{2}\big{)}}{N^{3}(1+N)^{3}}S_{1}^{2}-\frac{32(N-1)(2+N)}{3N^{2}(1+N)^{2}}$
$\displaystyle\times
S_{1}^{3}-\frac{8(2+N)\big{(}2-11N-16N^{2}+9N^{3}\big{)}}{N^{3}(1+N)^{3}}S_{2}+\frac{32\big{(}10+7N+7N^{2}\big{)}}{3N^{2}(1+N)^{2}}S_{3}$
$\displaystyle+\Biggl{(}-\frac{64\big{(}10+N+N^{2}\big{)}}{(N-1)N(1+N)(2+N)}+\frac{512}{N^{2}(1+N)^{2}}S_{1}\Biggr{)}S_{-2}+\frac{256}{N^{2}(1+N)^{2}}S_{-3}$
$\displaystyle-\frac{64(N-1)(2+N)}{N^{2}(1+N)^{2}}S_{2,1}-\frac{512}{N^{2}(1+N)^{2}}S_{-2,1}+\frac{192\big{(}-2-N-N^{2}\big{)}}{N^{2}(1+N)^{2}}\zeta_{3}\Biggr{]}$
$\displaystyle+{\color[rgb]{0,0,1}C_{A}^{2}T_{F}N_{F}}\Biggl{[}\frac{32P_{4}}{9N^{2}(1+N)^{2}}S_{2}+\frac{32P_{11}}{9N^{2}(1+N)^{2}}S_{-3}-\frac{64P_{11}}{9N^{2}(1+N)^{2}}S_{-2,1}$
$\displaystyle+\frac{16P_{13}}{9N^{2}(1+N)^{2}}S_{3}+\frac{2P_{76}}{27(N-1)N^{5}(1+N)^{5}(2+N)}+\Biggl{(}\frac{1280}{9}S_{2}-\frac{64}{3}S_{3}$
$\displaystyle-\frac{8P_{68}}{27(-1+N)N^{4}(1+N)^{4}(2+N)}-128\zeta_{3}\Biggr{)}S_{1}+\frac{64}{3}S_{-2}^{2}$
$\displaystyle+\Biggl{(}\frac{64P_{45}}{9(N-1)N^{2}(1+N)^{2}(2+N)}S_{1}-\frac{32P_{50}}{9(N-1)N^{3}(1+N)^{3}(2+N)}\Biggr{)}S_{-2}$
$\displaystyle+\frac{128\big{(}-3+2N+2N^{2}\big{)}}{N^{2}(1+N)^{2}}\zeta_{3}\Biggr{]}.$

The polynomials are given in Refs. [1, 2]. The small $z$ and large $z$ limits
of the anomalous dimensions have also been considered in explicit form. The
latter are related to the cusp anomalous dimensions. There are various partial
predictions on the small $z$ behaviour, however, not specifying the
factorization scheme used. These are theoretically interesting, but are not of
quantitative importance since it is known for long that subleading terms do
strongly modify the leading order effects over about three additional
subleading orders [48].

We have also recalculated the unpolarized two–loop anomalous dimensions using
the method of massless off–shell OMEs, [3]. Here new local operators
contribute [8, 16, 3] to cancel the gauge–variant contributions to obtain the
anomalous dimensions. We corrected a series of errors contained in [4, 5] and
provided all expansion coefficients emerging at two–loop order [3], which are
relevant for upcoming four–loop calculations.

##  5 Conclusions

We have calculated the three–loop polarized and unpolarized non–singlet
anomalous dimensions, including transversity, and as well the polarized
three–loop singlet anomalous dimensions applying the traditional method of
massless off–shell operator matrix elements, a gauge–variant method, which
requires special projections. The calculations have been performed without
assuming special conditions concerning the structure of the final result. The
method of arbitrary high Mellin moments plays a central role in these
computations in establishing the corresponding difference equations through
the method of guessing, after having exploited the equations for the master
integrals. These equations are solved subsequently using algorithms of
difference ring theory. We agree with all the results previously obtained in
the literature. The polarized non–singlet anomalous dimension
$\Delta\gamma_{qq}^{s,\rm NS}$, related to the polarized structure function
$g_{5}^{-}$, [25], has been calculated using the associated forward Compton
amplitude. In the unpolarized singlet case also new OMEs contribute. Here we
have calculated all contributions emerging at the two–loop level and corrected
results in the literature. The present method is suited to be expanded to the
four–loop level. The knowledge of the higher–loop anomalous dimensions form
one important asset for the precise description of the scaling violations of
the deep–inelastic structure functions, which provide an important way to
measure the QCD coupling constant $\alpha_{s}(M_{Z}^{2})$ [49] at highest
precision possible.

## References

  * [1] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald, Nucl. Phys. B 971 (2021) 115542 [arXiv: 2107.06267 [hep-ph]].
  * [2] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald, JHEP 01 (2022) 193 [arXiv:2111.12401 [hep-ph]].
  * [3] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald, Nucl. Phys. B 980 (2022) 115794 [arXiv:2202.03216 [hep-ph]].
  * [4] Y. Matiounine, J. Smith and W.L. van Neerven, Phys. Rev. D 57 (1998) 6701–6722 [arXiv:hep-ph/9801224 [hep-ph]];
  * [5] Y. Matiounine, J. Smith and W.L. van Neerven, Phys. Rev. D 58 (1998) 076002 [hep-ph/9803439].
  * [6] D.J. Gross and F. Wilczek, Phys. Rev. D 8 (1973) 3633–3652; Phys. Rev. D 9 (1974) 980–993;
H. Georgi and H.D. Politzer, Phys. Rev. D 9 (1974) 416–420;
K. Sasaki, Prog. Theor. Phys. 54 (1975) 1816–1827;
M.A. Ahmed and G.G. Ross, Phys. Lett. B 56 (1975) 385–390.

  * [7] E.G. Floratos, D. A.Ross and C.T. Sachrajda, Nucl. Phys. B 129 (1977) 66–88 Erratum: [Nucl. Phys. B 139 (1978) 545–546]; 152 (1979) 493–520;
A. Gonzalez-Arroyo and C. Lopez, Nucl. Phys. B 166 (1980) 429–459;
A. Gonzalez-Arroyo, C. Lopez and F. J. Yndurain, Nucl. Phys. B 159 (1979)
512–527; Nucl. Phys. B 153 (1979) 161–186;
G. Curci, W. Furmanski and R. Petronzio, Nucl. Phys. B 175 (1980) 27–92;
W. Furmanski and R. Petronzio, Phys. Lett. B 97 (1980) 437–442;
E.G. Floratos, C. Kounnas and R. Lacaze, Nucl. Phys. B 192 (1981) 417–462;
R. Mertig and W.L. van Neerven, Z. Phys. C 70 (1996) 637–654 [hep-
ph/9506451v2];
W. Vogelsang, Phys. Rev. D 54 (1996) 2023–2029 [hep-ph/9512218]; Nucl. Phys. B
475 (1996) 47–72 [hep-ph/9603366];
R.K. Ellis and W. Vogelsang, The Evolution of parton distributions beyond
leading order: The Singlet case, arXiv:hep-ph/9602356 [hep-ph];
S. Moch and J.A.M. Vermaseren, Nucl. Phys. B 573 (2000) 853–907 [hep-
ph/9912355].

  * [8] R. Hamberg and W.L. van Neerven, Nucl. Phys. B 379 (1992) 143–171.
  * [9] S. Moch, J.A.M. Vermaseren and A. Vogt, Nucl. Phys. B 688 (2004) 101–134 [hep-ph/0403192].
  * [10] A. Vogt, S. Moch and J.A.M. Vermaseren, Nucl. Phys. B 691 (2004) 129–181 [hep-ph/0404111].
  * [11] S. Moch, J.A.M. Vermaseren and A. Vogt, Nucl. Phys. B 889 (2014) 351–400 [arXiv: 1409.5131 [hep-ph]].
  * [12] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, Nucl. Phys. B 890 (2014) 48–151 [arXiv:1409.1135 [hep-ph]].
  * [13] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, Nucl. Phys. B 922 (2017) 1–40 [arXiv:1705.01508 [hep-ph]].
  * [14] A. Behring, J. Blümlein, A. De Freitas, A. Goedicke, S. Klein, A. von Manteuffel, C. Schneider and K. Schönwald, Nucl. Phys. B 948 (2019) 114753 [arXiv:1908.03779 [hep-ph]].
  * [15] T. Gehrmann et al., these Proceedings and in preparation.
  * [16] J.A. Dixon and J.C. Taylor, Nucl. Phys. B 78 (1974) 552–560;
H. Kluberg-Stern and J. B. Zuber, Phys. Rev. D 12 (1975) 482–488; 3159–3180;
S. Sarkar, Nucl. Phys. B 82 (1974) 447–460;
S. Sarkar and H. Strubbe, Nucl. Phys. B 90 (1975) 45–51;
S.D. Joglekar and B.W. Lee, Annals Phys. 97 (1976) 160–215;
S.D. Joglekar, Annals Phys. 108 (1977) 233–241; 109 (1977) 210–287;
R. Hamberg, Second order gluonic contributions to physical quantities, Thesis,
Univ. Leiden (1991);
J.C. Collins and R.J. Scalise, Phys. Rev. D 50 (1994) 4117–4136 [arXiv:hep-
ph/9403231 [hep-ph]];
B.W. Harris and J. Smith, Phys. Rev. D 51 (1995) 4550–4560 [arXiv:hep-
ph/9409405 [hep-ph]];
G. Falcioni and F. Herzog, JHEP 05 (2022) 177 [arXiv:2203.11181 [hep-ph]] and
these Proceedings.

  * [17] J.A.M. Vermaseren, Int. J. Mod. Phys. A 14 (1999) 2037–2076 [hep-ph/9806280].
  * [18] J. Blümlein and S. Kurth, Phys. Rev. D 60 (1999) 014018 [hep-ph/9810241].
  * [19] J. Blümlein, Prog. Part. Nucl. Phys. 69 (2013) 28–84 [arXiv:1208.6087 [hep-ph]].
  * [20] S.A. Larin, Phys. Lett. B 303 (1993) 113–118 [hep-ph/9302240].
  * [21] P. Nogueira, J. Comput. Phys. 105 (1993) 279–289.
  * [22] I. Bierenbaum, J. Blümlein and S. Klein, Nucl. Phys. B 820 (2009) 417–482 [arXiv:0904.3563 [hep-ph]].
  * [23] J.A.M. Vermaseren, New features of FORM, math-ph/0010025;
M. Tentyukov and J. A. M. Vermaseren, Comput. Phys. Commun. 181 (2010)
1419–1427 [hep-ph/0702279].

  * [24] T. van Ritbergen, A.N. Schellekens and J.A.M. Vermaseren, Int. J. Mod. Phys. A 14 (1999) 41–96 [arXiv:hep-ph/9802376 [hep-ph]].
  * [25] J. Blümlein and N. Kochelev, Nucl. Phys. B 498 (1997) 285–309 [arXiv:hep-ph/9612318 [hep-ph]].
  * [26] H.D. Politzer, Phys. Rept. 14 (1974) 129–180
  * [27] J. Ablinger, J. Blümlein, C. Raab, C. Schneider and F. Wißbrock, Nucl. Phys. B 885 (2014) 409–447 [arXiv:1403.1137 [hep-ph]].
  * [28] P. Marquard and D. Seidel, The Crusher algorithm, unpublished.
  * [29] J. Lagrange, Nouvelles recherches sur la nature et la propagation du son, Miscellanea Taurinensis, t. II, 1760-61; Oeuvres t. I, p. 263;
C.F. Gauß, Theoria attractionis corporum sphaeroidicorum ellipticorum
homogeneorum methodo novo tractate, Commentationes societas scientiarum
Gottingensis recentiores, Vol III, 1813, Werke Bd. V pp. 5-7;
G. Green, Essay on the Mathematical Theory of Electricity and Magnetism,
Nottingham, 1828 [Green Papers, pp. 1-115];
M. Ostrogradski, Mem. Ac. Sci. St. Peters., 6 (1831) 39–53;
S. Laporta, Int. J. Mod. Phys. A 15 (2000) 5087–5159 [hep-ph/0102033].

  * [30] K.G. Chetyrkin and F.V. Tkachov, Nucl. Phys. B 192 (1981) 159–204.
  * [31] D.J. Gross and F. Wilczek, Phys. Rev. Lett. 30 (1973) 1343–1346;
H.D. Politzer, Phys. Rev. Lett. 30 (1973) 1346–1349;
W.E. Caswell, Phys. Rev. Lett. 33 (1974) 244–246;
D.R.T. Jones, Nucl. Phys. B 75 (1974) 531–538;
O.V. Tarasov, A.A. Vladimirov and A.Y. Zharkov, Phys. Lett. B 93 (1980)
429–432;
S.A. Larin and J.A.M. Vermaseren, Phys. Lett. B 303 (1993) 334–336 [hep-
ph/9302208];
T. van Ritbergen, J.A.M. Vermaseren and S.A. Larin, Phys. Lett. B 400 (1997)
379–384 [hep-ph/9701390];
M. Czakon, Nucl. Phys. B 710 (2005) 485–498 [hep-ph/0411261];
P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, Phys. Rev. Lett. 118 (2017) no.8,
082002 [arXiv:1606.08659 [hep-ph]];
P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, Phys. Rev. Lett. 118 (2017) no.8,
082002 [arXiv:1606.08659 [hep-ph]];
F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, JHEP 02 (2017)
090 [arXiv:1701.01404 [hep-ph]];
T. Luthe, A. Maier, P. Marquard and Y. Schröder, JHEP 10 (2017) 166
[arXiv:1709.07718 [hep-ph]];
K.G. Chetyrkin, G. Falcioni, F. Herzog and J.A.M. Vermaseren, JHEP 10 (2017)
179 [arXiv:1709.08541 [hep-ph]].

  * [32] K.G. Chetyrkin, G. Falcioni, F. Herzog and J.A.M. Vermaseren, JHEP 1710 (2017) 179 Addendum: [JHEP 1712 (2017) 006] [arXiv:1709.08541 [hep-ph]];
K.G. Chetyrkin, Nucl. Phys. B 710 (2005) 499–510 [arXiv:hep-ph/0405193 [hep-
ph]]. E. Egorian and O.V. Tarasov, Teor. Mat. Fiz. 41 (1979) 26–32 [Theor.
Math. Phys. 41 (1979) 863–867];
T. Luthe, A. Maier, P. Marquard and Y. Schröder, JHEP 1701 (2017) 081
[arXiv:1612.05512 [hep-ph]].

  * [33] J. Blümlein and C. Schneider, Phys. Lett. B 771 (2017) 31–36 [arXiv:1701.04614 [hep-ph]].
  * [34] J. Blümlein, P. Marquard and C. Schneider, PoS (RADCOR2019) 078 [arXiv:1912.04390 [cs.SC]].
  * [35] F.V. Tkachov, Phys. Lett. B 100 (1981) 65–68; Theor. Math. Phys. 56 (1983) 866–870; Teor. Mat. Fiz. 56 (1983) 350–356;
See also: T. Gehrmann, E.W.N. Glover, T. Huber, N. Ikizlerli and C. Studerus,
JHEP 06 (2010) 094 [arXiv:1004.3653 [hep-ph]];
R.N. Lee, A.V. Smirnov and V.A. Smirnov, Nucl. Phys. B Proc. Suppl. 205-206
(2010) 308–313 [arXiv:1005.0362 [hep-ph]];
R.N. Lee and V.A. Smirnov, JHEP 02 (2011) 102 [arXiv:1010.1334 [hep-ph]].

  * [36] M. Kauers, Guessing Handbook, JKU Linz, Technical Report RISC 09-07.
  * [37] J. Blümlein, M. Kauers, S. Klein and C. Schneider, Comput. Phys. Commun. 180 (2009) 2143–2165 [arXiv:0902.4091 [hep-ph]].
  * [38] Sage, http://www.sagemath.org/
  * [39] M. Kauers, M. Jaroschek, and F. Johansson, in: Computer Algebra and Polynomials, Editors: J. Gutierrez, J. Schicho, Josef, M. Weimann, Eds.. Lecture Notes in Computer Science 8942 (Springer, Berlin, 2015) 105–125 [arXiv:1306.4263 [cs.SC]].
  * [40] J. Blümlein, D.J. Broadhurst and J.A.M. Vermaseren, Comput. Phys. Commun. 181 (2010), 582-625 [arXiv:0907.2557 [math-ph]].
  * [41] M. Karr, J. ACM 28 (1981) 305–350;
M. Bronstein, J. Symbolic Comput. 29 (2000) no. 6 841–877;
C. Schneider, Symbolic Summation in Difference Fields, Ph.D. Thesis RISC,
Johannes Kepler University, Linz technical report 01–17 (2001);
C. Schneider, An. Univ. Timisoara Ser. Mat.-Inform. 42 (2004) 163–179;
C. Schneider, J. Differ. Equations Appl. 11 (2005) 799–821;
C. Schneider, Appl. Algebra Engrg. Comm. Comput. 16 (2005) 1–32;
C. Schneider, J. Algebra Appl. 6 (2007) 415–441;
C. Schneider, Clay Math. Proc. 12 (2010) 285–308 [arXiv:0904.2323 [cs.SC]];
C. Schneider, Ann. Comb. 14 (2010) 533–552 [arXiv:0808.2596];
C. Schneider, in: Computer Algebra and Polynomials, Applications of Algebra
and Number Theory, J. Gutierrez, J. Schicho, M. Weimann (ed.), Lecture Notes
in Computer Science (LNCS) 8942 (2015) 157–191 [arXiv:1307.7887 [cs.SC]];
C. Schneider, J. Symb. Comput. 72 (2016) 82–127 [arXiv:1408.2776 [cs.SC]]. C.
Schneider, J. Symb. Comput. 80 (2017) 616–664 [arXiv:1603.04285 [cs.SC]]. S.A.
Abramov, M. Bronstein, M. Petkovšek, Carsten Schneider, J. Symb. Comput. 107
(2021) 23–66 [arXiv:2005.04944 [cs.SC]].

  * [42] C. Schneider, Sém. Lothar. Combin. 56 (2007) 1–36 article B56b.
  * [43] C. Schneider, Simplifying Multiple Sums in Difference Fields, in: Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions Texts and Monographs in Symbolic Computation eds. C. Schneider and J. Blümlein (Springer, Wien, 2013) 325–360 [arXiv:1304.4134 [cs.SC]].
  * [44] J. Ablinger, J. Blümlein and C. Schneider, J. Phys. Conf. Ser. 523 (2014) 012060 [arXiv:1310.5645 [math-ph]];
J. Ablinger, PoS (LL2014) 019 [arXiv:1407.6180[cs.SC]]; A Computer Algebra
Toolbox for Harmonic Sums Related to Particle Physics, Diploma Thesis, JKU
Linz, 2009, arXiv:1011.1176[math-ph]; Computer Algebra Algorithms for Special
Functions in Particle Physics, Ph.D. Thesis, Linz U. (2012)
arXiv:1305.0687[math-ph]; PoS (LL2016) 067; Experimental Mathematics 26 (2017)
[arXiv:1507.01703 [math.CO]]; PoS (RADCOR2017) 001 [arXiv:1801.01039 [cs.SC]];
arXiv:1902.11001 [math.CO]; PoS (LL2018) 063;
J. Ablinger, J. Blümlein and C. Schneider, J. Math. Phys. 54 (2013) 082301
[arXiv:1302.0378 [math-ph]]; J. Math. Phys. 52 (2011) 102301 [arXiv:1105.6063
[math-ph]];
J. Ablinger, J. Blümlein, C.G. Raab and C. Schneider, J. Math. Phys. 55 (2014)
112301 [arXiv:1407.1822 [hep-th]].

  * [45] E. Remiddi and J.A.M. Vermaseren, Int. J. Mod. Phys. A 15 (2000) 725–754 [hep-ph/9905237].
  * [46] J. Blümlein, Comput. Phys. Commun. 159 (2004) 19–54 [arXiv:hep-ph/0311046 [hep-ph]].
  * [47] J. Blümlein, Comput. Phys. Commun. 180 (2009) 2218–2249 [arXiv:0901.3106 [hep-ph]].
  * [48] J. Blümlein and A. Vogt, Phys. Lett. B 370 (1996) 149–155 [hep-ph/9510410]; Phys. Lett. B 386 (1996), 350-358 [arXiv:hep-ph/9606254 [hep-ph]]; Phys. Rev. D 58 (1998) 014020 [arXiv:hep-ph/9712546 [hep-ph]];
J. Blümlein, Lect. Notes Phys. 546 (2000) 42–57, [arXiv:hep-ph/9909449 [hep-
ph]].

  * [49] S. Bethke et al., Workshop on Precision Measurements of $\alpha_{s}$, arXiv:1110.0016 [hep-ph];
S. Moch et al., High precision fundamental constants at the TeV scale,
arXiv:1405.4781 [hep-ph];
S. Alekhin, J. Blümlein and S.O. Moch, Mod. Phys. Lett. A 31 (2016) no.25,
1630023.


# Non-Markovian reversible diffusion-influenced
reactions in two dimensions

Thorsten Prüstel  Laboratory of Systems Biology
National Institute of Allergy and Infectious Diseases
National Institutes of Health  Martin Meier-Schellersheim  Laboratory of
Systems Biology
National Institute of Allergy and Infectious Diseases
National Institutes of Health

###### Abstract

We investigate the reversible diffusion-influenced reaction of an isolated
pair in the presence of a non-Markovian generalization of the backreaction
boundary condition in two space dimensions. Following earlier work by Agmon
and Weiss, we consider residence time probability densities that decay slower
than an exponential and that are characterized by a parameter $0<\sigma\leq
1$. We calculate an exact expression for the probability $S(t|\ast)$ that the
initially bound particle is unbound, which is valid for arbitrary $\sigma$ and
for all times. Furthermore, we derive an approximate solution for long times.
We show that the ultimate fate of the bound state is complete dissociation, as
in the 2D Markovian case. However, the limiting value is approached quite
differently: Instead of a $\sim t^{-1}$ decay, we obtain $1-S(t|\ast)\sim
t^{-\sigma}\ln t$.

11footnotetext: Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

It has been shown that the ultimate fate of an isolated pair for the
reversible diffusion-influenced reaction is independent from the
dimensionality of space. In fact, the ultimate escape probability is unity in
one, two and three space dimensions [2, 3, 9, 12]. However, the dimensionality
does have an influence on how the limiting value for $S(t|\ast)$, which
denotes the probability that the initially bound particle is unbound, is
approached: Denoting the number of dimensions by $d$, the decay goes as
$1-S(t|\ast)\sim t^{-d/2}$, which implies for instance considerable
ramifications for the off-rate. This dependence from the dimension can be
traced back to the return to the origin probability of the underlying random
walk. In one and two dimensions (2D), Brownian motion is recurrent, i.e. the
return to the origin happens with probability one, while in 3D it is
transient, meaning that the non-return probability does not vanish.

It is reasonable that, apart from the diffusive motion of the particles, the
type of their interaction potentially influences the time-dependence of the
probability $S(t|\ast)$ and, in particular, its long-time behavior. Indeed, in
most cases one studied the problem by tacitly assuming a Markovian residence
time probability density for the bound state. However, Agmon and Weiss [4]
demonstrated that abandoning the Markovian assumption changes the time-
dependence of $S(t|\ast)$ in a substantial way. They found a phase transition
as a function of the parameter $\sigma$, which is a measure for the deviation
from the Markovian density (where $\sigma=1$ corresponds to the Markovian
limit). In one dimension, they found a first-order phase transition in which
the long-time behavior of $S(t|\ast)$ undergoes a discontinuous change at
$\sigma=1/2$. In contrast, in 3D they observed a second-order phase transition
at $\sigma=1$ and the ultimate fate is again complete dissociation,
independent from $\sigma$.

In what follows, we will extend the analysis presented in Ref. [4] to the 2D
case. We will investigate non-Markovian dissociation for an isolated pair of
reversibly binding particles which move in an infinitely extended two
dimensional plane. The case of two dimensions is special, because it is the
critical dimension with regard to recurrence and transience [14]. We will show
whether even weak recurrence is sufficient to obtain a first-order phase
transition.

We would like to point out that the consideration of non-Markovian
dissociation is not only theoretically interesting. In particular in view of
the intricate and heterogeneous composition of biological cell membranes [6],
generalizations of the Markovian behavior could potentially be significant for
more realistic descriptions of the corresponding local biochemistries.

The organization of the paper is as follows. First, we will introduce the
general theoretical context, especially the Smoluchowski equation and
backreaction boundary condition (BC). Next, we will discuss memory effects and
a non-Markovian generalization of the backreaction BC. In theses parts of our
presentation, we shall closely follow Agmon and Weiss [4]. Then, we will apply
the general formalism to the 2D case. We will obtain an integral
representation for $S(t|\ast)$ in the time domain, which is valid for all
times. Furthermore, we will present an approximate solution for the long-time
behavior of the probability $S(t|\ast)$. Finally, we discuss similarities and
differences to the non-Markovian 1D and 3D cases.

##  2 Smoluchowski equation and backreaction BC

Solutions of the diffusion (Smoluchowski) equation [15, 7, 13] can be employed
to study diffusion-influenced reactions of an isolated pair. The different
types of chemical reactions are taken into account by imposing certain
boundary conditions at the encounter distance. Typically, one considers the
following scenario [2, 3, 9, 12]. Two spherical (or rather disklike in 2D)
particles $A$ and $B$ with diffusion constants $D_{A}$ and $D_{B}$,
respectively, may associate when their separation equals the encounter
distance $a$ to form a bound molecule $AB$. Such a system may equivalently be
described as the diffusion of a point-like particle with diffusion constant
$D=D_{A}+D_{B}$ around a static disk with radius $a$.

The irreversible association reaction is described by the radiation BC that is
characterized by an intrinsic association constant $\kappa_{a}$ [8]. In the
reversible case, the bound state may dissociate again to form an unbound pair
$A+B$ and the radiation boundary condition has to be generalized
correspondingly. The so-called backreaction BC [2, 10, 11, 5, 3, 9, 12] has
been used extensively to take into account reversible association-dissociation
reactions, but in its conventional form it assumes a Markovian probability
density for the residence time [4]

$\psi(t)=\kappa_{d}e^{-\kappa_{d}t},$ (2.1)

or, equivalently, in the Laplace domain

$\tilde{\psi}(s)=\frac{\kappa_{d}}{\kappa_{d}+s}.$ (2.2)

Eq. (2.1) corresponds to the situation where the bound molecule can be
sufficiently well described by a single, discrete state. However, the bound
state is typically comprised of a large number of different internal states.
It follows that the associated probability density $\psi(t)$ for residence
time in the set of all possible bound states is given by a superposition of
exponentials, which correspond to the individual states

$\psi(t)=\sum^{\infty}_{n=0}\varphi_{n}e^{-\lambda_{n}t}.$ (2.3)

Only if the time scales set by the $\lambda_{n}$ parameters are sufficiently
separated, i.e. the decay of the ground state (=dissociation) is much slower
than the decay of all other internal states, Eq. (2.3) can be replaced by Eq.
(2.1) and the Markovian description is justified. However, if there is no
sufficient separation of the involved time scales, the effective description
of all internal states as one single state cannot be Markovian any more.
Furthermore, as already mentioned in the introduction, a non-Markovian
deviation could be caused by a complex environment, for instance in biological
cellular systems.

To proceed we consider the probability density function (PDF) $p(r,t|\ast)$,
which yields the probability to find the particle at a distance equal to $r$,
given that it was initially in the bound state. The diffusion (or
Smoluchowski) equation in 2D [4, 3] governs the time evolution of
$p(r,t|\ast)$

$\frac{\partial}{\partial t}p(r,t|\ast)=D\bigg{(}\frac{\partial^{2}}{\partial
r^{2}}p(r,t|\ast)+\frac{1}{r}\frac{\partial}{\partial
r}p(r,t|\ast)\bigg{)},\quad r\geq a.$ (2.4)

The initial condition is

$p(r,t=0|\ast)=0,\quad r>a.$ (2.5)

The PDF we are interested in is only defined for $r\geq a>0$ and one has to
impose a BC for $r=a$ specifying the behavior at the encounter distance. The
(Markovian) backreaction BC incorporates association and dissociation by
relating the diffusional flux at the surface of the ”interaction disc” and the
probability $S(t|\ast)$ that the initially bound particle is unbound at time
$t>0$ as follows [2, 3, 9]

$\displaystyle 2\pi aD\frac{\partial}{\partial
r}p(r,t|\ast)|_{r=a}=\kappa_{a}p(a,t|\ast)-\kappa_{d}[1-S(t|\ast)].$ (2.6)

Note that in 2D the diffusional flux is given by

$\displaystyle-J(a|\ast)=2\pi aD\frac{\partial}{\partial
r}p(r,t|\ast)|_{r=a},$ (2.7)

and the probability $S(t|\ast)$ is defined by

$\displaystyle S(t|\ast)$ $\displaystyle=$ $\displaystyle
2\pi\int^{\infty}_{a}p(r,t|\ast)rdr.$ (2.8)

Integrating the diffusion equation over space yields

$J(a|\ast)=\frac{\partial}{\partial t}S(t|\ast),$ (2.9)

which translates in the Laplace domain to

$\tilde{J}(a,s|\ast)=s\tilde{S}(s|\ast),$ (2.10)

where we have used that $S(0|\ast)$ vanishes. Hence, the BC Eq. (2.6) becomes
in the Laplace domain [4]

$\tilde{J}(a,s|\ast)=\frac{\kappa_{d}-s\kappa_{a}\tilde{p}(a,s|\ast)}{s+\kappa_{d}}.$
(2.11)

For our analysis, we will make use of the following central identity [4] that
relates in the Laplace domain the probability that the particle is not bound
to the PDF $\tilde{p}_{\text{ref}}(a,s|a)$, which corresponds to nonreactive
diffusion and solves the diffusion equation subject to a reflecting BC,
$\tilde{J}(a,s|a)=0$,

$s\tilde{S}(s|\ast)=\frac{\kappa_{d}}{\kappa_{d}+s[1+\kappa_{a}\tilde{p}_{\text{ref}}(a,s|a)]}.$
(2.12)

A salient feature of Eq. (2.12) is that it is valid in any dimension, thereby
allowing to find an expression for $\tilde{S}$ once
$\tilde{p}_{\text{ref}}(a,s|a)$ is known. In 2D, one has [12]

$\tilde{p}_{\text{ref}}(a,s|a)=\frac{1}{2\pi
D}\bigg{[}I_{0}(qa)K_{0}(qa)+K_{0}(qa)K_{0}(qa)\frac{I_{1}(qa)}{K_{1}(qa)}\bigg{]},$
(2.13)

where $I_{0}(z),I_{1}(z),K_{z}(x),K_{1}(z)$ refer to the modified Bessel
function of first and second kind and zero and first order, respectively [1,
Sect. 9.6]. Using $I_{0}(z)K_{1}(z)+I_{1}(z)K_{0}(z)=z^{-1}$ and combining
Eqs. (2.12) and (2.13) we arrive at

$s\tilde{S}(s|\ast)=\frac{\kappa_{d}qK_{1}(qa)}{qK_{1}(qa)(s+\kappa_{d})+\kappa_{a}/(2\pi
aD)sK_{0}(qa)}.$ (2.14)

To find an expression for the probability that the particle is unbound in the
time domain, we could invert the Laplace transform explicitly by calculating a
Bromwich contour integral. However, in the Markovian backreaction BC case, an
expression has already been derived by using alternatively the Green’s
function $p(r,t|t_{0})$ of the diffusion equation subject to a backreaction
BC. The obtained result is [12]

$S(t|*)=1-2\frac{\kappa_{a}\kappa_{d}}{\pi^{3}a^{2}D^{2}}\int^{\infty}_{0}e^{-Dtx^{2}}\frac{1}{\alpha^{2}(x)+\beta^{2}(x)}\frac{1}{x}dx,$
(2.15)

where

$\displaystyle\alpha(x)$ $\displaystyle=$
$\displaystyle(x^{2}-\kappa_{D})J_{1}(xa)+hxJ_{0}(xa),$ (2.16)
$\displaystyle\beta(x)$ $\displaystyle=$
$\displaystyle(x^{2}-\kappa_{D})Y_{1}(xa)+hxY_{0}(xa),$ (2.17) $\displaystyle
h$ $\displaystyle=$ $\displaystyle\frac{\kappa_{a}}{2\pi aD},$ (2.18)
$\displaystyle\kappa_{D}$ $\displaystyle=$
$\displaystyle\frac{\kappa_{d}}{D}.$ (2.19)

$J_{n}(z),Y_{n}(z)$ denote the Bessel functions of first and second kind,
respectively [1, Sect. 9.1].

We now seek the non-Markovian generalization of Eq. (2.15). To this end, we
have to find the non-Markovian analogue of the term $\kappa_{d}[1-S(t|\ast)]$
in the BC Eq. (2.6).

###  2.1 Non-Markovian backreaction BC

In Ref. [4] it was shown that the backreaction BC can also be written in terms
of the probability density of the residence time as follows

$J(a,t|\ast)=-\kappa_{a}p(a,t|\ast)+\psi(t)+\kappa_{a}\int^{t}_{0}p(a,\tau|\ast)\psi(t-\tau)d\tau.$
(2.20)

Application of the Laplace transform leads to

$\tilde{J}(a,s|\ast)=\tilde{\psi}(s)-\kappa_{a}[1-\tilde{\psi}(s)]p(a,s|\ast).$
(2.21)

Similarly, the Laplace transform of $S(t|\ast)$ may be written as [4]

$s\tilde{S}(s|\ast)=\frac{\tilde{\psi}(s)}{1+\kappa_{a}[1-\tilde{\psi}(s)]\tilde{p}_{\text{ref}}(a,s|a)}.$
(2.22)

Here, two remarks are in order. First, upon using the Markovian version of the
residence time probability density Eq. (2.2) in Eqs. (2.21), (2.22), one
recovers Eqs. (2.11) and (2.12). On the other hand, starting from Eqs. (2.11)
and (2.12), one can simply substitute $\kappa_{d}$ by
$s\tilde{\psi}(s)/[1-\tilde{\psi}(s)]$ [3] to arrive at Eqs. (2.21), (2.22).

Next, we turn our attention to the concrete form of the residence time density
$\psi(t)$. In Ref. [4], the following choice of densities was employed

$\tilde{\psi}(s)=\frac{1}{1+(s\kappa^{-1})^{\sigma}},\quad 0<\sigma\leq 1.$
(2.23)

For a discussion and motivation for this particular form, we refer to [4].
Here, we only note that for $\sigma=1$ and by identifying $\kappa$ with
$\kappa_{d}$, one recovers Eq. (2.2). However, for $\sigma<1$, $\psi(t)$
decays slower than an exponential.

##  3 Non-Markovian survival probability in 2D

In what follows, we will use Eqs. (2.22), (2.13) and (2.23) to derive an
expression for $S(t|\ast)$ for arbitrary $0<\sigma\leq 1$ and all times. We
obtain

$\tilde{S}(s|\ast)=\frac{1}{s}\frac{\kappa^{\sigma}qK_{1}(qa)}{qK_{1}(qa)(s^{\sigma}+\kappa^{\sigma})+hs^{\sigma}K_{0}(qa)},$
(3.1)

where $q=\sqrt{s/D}$. The inversion theorem for the Laplace transformation can
be applied to find the corresponding expression of $\tilde{S}(s|\ast)$ in the
time domain

$S(t|\ast)=\frac{1}{2\pi
i}\int^{\gamma+i\infty}_{\gamma-i\infty}e^{st}\,\tilde{S}(s|\ast)ds.$ (3.2)

To calculate the Bromwich contour integral, we first note that
$\tilde{S}(s|\ast)$ has a branch point at $s=0$. Therefore, we use the contour
of FIG. 1 with a branch cut along the negative real axis. We note that the
contribution from the small circle around the origin does not vanish, but
yields 1, which can be seen by the limiting forms of the modified Bessel
functions for small arguments [1, Eqs. (9.6.7)-(9.6.9)]. Therefore, we obtain

$\displaystyle\int^{\gamma+i\infty}_{\gamma-i\infty}e^{st}\,\tilde{S}(s|\ast)ds$
$\displaystyle=$ $\displaystyle 2\pi
i-\int_{\mathcal{C}_{2}}e^{st}\,\tilde{S}(s|\ast)ds-\int_{\mathcal{C}_{4}}e^{st}\,\tilde{S}(s|\ast)ds.$
(3.3)

It remains to calculate the integrals along $\mathcal{C}_{2},\mathcal{C}_{4}$.
To this end, we choose $s=Dx^{2}e^{i\pi}$ and use [7, Append. 3, Eqs. (25),
(26))]

$\displaystyle I_{n}(xe^{\pm\pi i/2})$ $\displaystyle=$ $\displaystyle e^{\pm
n\pi i/2}J_{n}(x),$ (3.4) $\displaystyle K_{n}(xe^{\pm\pi i/2})$
$\displaystyle=$ $\displaystyle\pm\frac{1}{2}\pi ie^{\mp n\pi
i/2}[-J_{n}(x)\pm iY_{n}(x)].$ (3.5)

It follows that

$\displaystyle\int_{\mathcal{C}_{2}}e^{st}\,\tilde{S}(s|\ast)ds=2\kappa^{\sigma}_{D}\int^{\infty}_{0}e^{-Dx^{2}t}\frac{[Y_{1}(xa)+iJ_{1}(xa)][\beta_{\sigma}(x)-i\alpha_{\sigma}(x)]}{\alpha_{\sigma}(x)^{2}+\beta_{\sigma}(x)^{2}}\frac{dx}{x},$
(3.6)

where $\kappa^{\sigma}_{D}=\kappa^{\sigma}/D^{\sigma}$ and we have defined

$\displaystyle\alpha_{\sigma}$ $\displaystyle=$ $\displaystyle
x^{2\sigma-1}\bigg{[}-\big{(}xJ_{1}(xa)+hJ_{0}(xa)\big{)}\cos(\pi\sigma)-\big{(}xY_{1}(xa)+hY_{0}(xa)\big{)}\sin(\pi\sigma)\bigg{]}$
(3.7) $\displaystyle-\kappa^{\sigma}_{D}J_{1}(xa),$
$\displaystyle\beta_{\sigma}$ $\displaystyle=$ $\displaystyle
x^{2\sigma-1}\bigg{[}\big{(}xJ_{1}(xa)+hJ_{0}(xa)\big{)}\sin(\pi\sigma)-\big{(}xY_{1}(xa)+hY_{0}(xa)\big{)}\cos(\pi\sigma)\bigg{]}$
(3.8) $\displaystyle-\kappa^{\sigma}_{D}Y_{1}(xa).$

We observe that $\alpha_{\sigma=1}=\alpha,\beta_{\sigma=1}=\beta$, cp. Eqs.
(2.16), (2.17).

To evaluate the integral along the contour $\mathcal{C}_{4}$ we choose
$p=Dx^{2}e^{-i\pi}$ and after an analogous calculation one finds that
$\int_{\mathcal{C}_{2}}e^{st}\,\tilde{S}(s|\ast)ds=-\bigl{(}\int_{\mathcal{C}_{4}}e^{st}\,\tilde{S}(s|\ast)ds\bigr{)}^{\star},$
where $\star$ means complex conjugation. Thus, we have

$\displaystyle
S(t|\ast)-1=-\frac{1}{\pi}\Im\biggl{(}\int_{\mathcal{C}_{2}}e^{st}\,\tilde{S}(s|\ast)ds\biggr{)},$
(3.9)

and arrive finally at the following expression for the probability to find the
initially bound particle unbound in the time domain

$\displaystyle
S(t|\ast)=1+\frac{2}{\pi}\kappa^{\sigma}_{D}\bigg{[}\frac{2h}{a\pi}\cos(\pi\sigma)Q_{1}+\sin(\pi\sigma)Q_{2}\bigg{]},$
(3.10)

where

$\displaystyle Q_{1}$ $\displaystyle=$
$\displaystyle\int^{\infty}_{0}e^{-Dtx^{2}}\frac{x^{2\sigma-3}}{\alpha^{2}_{\sigma}(x)+\beta^{2}_{\sigma}(x)}dx,$
(3.11) $\displaystyle Q_{2}$ $\displaystyle=$
$\displaystyle-\int^{\infty}_{0}e^{-Dtx^{2}}\frac{x^{2\sigma-2}\Omega(x)}{\alpha^{2}_{\sigma}(x)+\beta^{2}_{\sigma}(x)}dx,$
(3.12)

and we have introduced

$\displaystyle\Omega(x)=x\big{[}J_{1}^{2}(xa)+Y_{1}^{2}(xa)\big{]}+h\big{[}J_{0}(xa)J_{1}(xa)+Y_{0}(xa)Y_{1}(xa)\big{]}.$
(3.13)

Eq. (3.10) shows that the parameter $\sigma$ appears as a sort of mixing angle
for the two qualitatively different contributions $Q_{1},Q_{2}$. We note that
the term $Q_{2}$ is foreign to the Markovian expression Eq. (2.15).
Correspondingly, the second term vanishes for $\sigma=1$, while one recovers
the Markovian limiting case Eq. (2.15) due to the first term $Q_{1}$ (we set
$\kappa:=\kappa_{d})$.

Furthermore, we may already conclude from Eq. (3.10) that in 2D the ultimate
escape probability is unity for all $\sigma$, like in the Markovian and in the
3D non-Markovian case. This behavior, however, is quite different to the non-
Markovian 1D case, where the escape probability depends on the parameter
$\sigma$ [4].

Finally, we are interested, how the limiting value is approached. To this end,
we derive the asymptotic behavior for $t\rightarrow\infty$. Starting point is
the expression for $\tilde{S}(s|\ast)$ Eq. (3.1). Using the series expansions
of the modified Bessel functions for small arguments [1, Eqs.
(9.6.10)-(9.6.13)], we obtain

$\tilde{S}(s|\ast)=\frac{1}{s}+\frac{ha}{\kappa^{\sigma}}s^{\sigma-1}\ln\bigg{(}\frac{1}{2}e^{\gamma_{E}}qa\bigg{)}-\frac{s^{\sigma-1}}{\kappa^{\sigma}}+\ldots,$
(3.14)

where $\gamma_{E}$ denotes Euler’s number $\gamma_{E}=0.5772156649\ldots$ [1,
Eqs. (6.1.3)]. For $\sigma=1$ we obtain in the time domain

$S(t|\ast)=1-\frac{\kappa_{a}}{\kappa_{d}}\frac{1}{4\pi Dt}+\ldots,$ (3.15)

which is the known Markovian result. If $0<\sigma<1$, we use [1, Eqs.
(29.3.99), (29.3.7)]

$\displaystyle\mathcal{L}^{-1}\bigg{[}\frac{1}{s^{k}}\ln s\bigg{]}$
$\displaystyle=$ $\displaystyle\frac{t^{k-1}}{\Gamma(k)}\bigg{[}\Psi(k)-\ln
t\bigg{]},\quad k>0,$ (3.16)
$\displaystyle\mathcal{L}^{-1}\bigg{[}\frac{1}{s^{k}}\bigg{]}$
$\displaystyle=$ $\displaystyle\frac{t^{k-1}}{\Gamma(k)},\quad k>0,$ (3.17)

where $\mathcal{L}^{-1},\Gamma,\Psi$ denote the inverse Laplace transform, the
gamma function and the digamma function [1, Eqs. (6.3.1)], respectively.
Hence, we find

$S(t|\ast)=1-\frac{\kappa_{a}}{\kappa^{\sigma}}\frac{1}{4\pi
D}\frac{t^{-\sigma}}{\Gamma(1-\sigma)}\ln
t+C\frac{t^{-\sigma}}{\Gamma(1-\sigma)}+\ldots,$ (3.18)

where we introduced the constant

$C=\frac{\kappa_{a}}{\kappa^{\sigma}}\frac{1}{2\pi
D}\bigg{[}\frac{\Psi(1-\sigma)}{2}+\ln\bigg{(}\frac{1}{2}e^{\gamma_{E}}D^{-1/2}a\bigg{)}\bigg{]}-\frac{1}{\kappa^{\sigma}}.$
(3.19)

We find again that, independent from the parameter $\sigma$, the ultimate fate
is complete dissociation, although the underlying random walk is recurrent. In
contrast, in 1D the parameter $\sigma$ decides if particles escape or remain
bound at long times. In this regard, the 2D case is more similar to the (non-
Markovian) 3D case. However, it is very different to the 3D case with respect
to how the limiting value is approached. In 3D, it was shown that the decay
goes as $t^{-\sigma}$ at $\sigma<1$, but as $t^{-3/2}$ at $\sigma=1$, which
means there is a second-order phase transition. However, in 2D, even the
functional form is different: While we observe a $t^{-\sigma}\ln t$ decay at
$\sigma<1$, we find $t^{-1}$ at $\sigma=1$.

### Acknowledgments

This research was supported by the Intramural Research Program of the NIH,
National Institute of Allergy and Infectious Diseases.

We would like to thank Bastian R. Angermann and Frederick Klauschen for
stimulating discussions.

Figure 1: Integration contour used in Eq. (3.3).

## References

  * [1] M. Abramowitz and I.A. Stegun.  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover, New York, 1965.
  * [2] N. Agmon.  J. Chem. Phys., 81:2811, 1984.
  * [3] N. Agmon and A. Szabo.  J. Chem. Phys., 92:5270, 1990.
  * [4] N. Agmon and G.H. Weiss.  J. Chem. Phys., 91:6937, 1989.
  * [5] Pines E. Agmon, N. and D. Huppert.  J. Chem. Phys., 88:5631, 1988.
  * [6] I. Bethani, S.S. Skanland, I. Dikic, and A. Acker-Palmer.  EMBO J., 29:2677, 2010.
  * [7] H.S. Carslaw and J.C. Jaeger.  Conduction of Heat in Solids.  Clarendon Press, New York, 1986.
  * [8] F.C. Collins and G.E. Kimball.  J. Colloid Sci., 4:425, 1949.
  * [9] H. Kim and K.J. Shin.  Phys. Rev. Lett., 82:1578, 1999.
  * [10] E. Pines and D. Huppert.  Chem. Phys. Lett., 126:88, 1986.
  * [11] E. Pines, D. Huppert, and N. Agmon.  J. Chem. Phys., 88:5620, 1988.
  * [12] T. Prüstel and M. Meier-Schellersheim.  J. Chem. Phys., 137:054104, 2012.
  * [13] S. A. Rice.  Diffusion Limited Reactions.  Elsevier, New York, 1985.
  * [14] D. Toussaint and F. Wilczek.  J. Chem. Phys., 78:2642, 1983.
  * [15] M. von Smoluchowski.  Z. Phys. Chem., 92:129, 1917.


# Water Preservation in Soan River Basin Using Deep Learning Techniques

Sadaqat ur Rehman
Department of Electronic Engineering
Tsinghua University
Beijing, 10084
<EMAIL_ADDRESS>
&Zhongliang Yang
Department of Electronic Engineering
Tsinghua University
Beijing, 10084
<EMAIL_ADDRESS>
&Muhammad Shahid
State Key Laboratory of Hydroscience and Engineering
Tsinghua University
Beijing, 10084
&Nan Wei
National Key Laboratory for Software Technology
Nanjing University
Nanjing, 210023
&Yongfeng Huang
Department of Electronic Engineering
Tsinghua University
Beijing, 10084
<EMAIL_ADDRESS>
&Muhammad Waqas
Department of Electronic Engineering
Tsinghua University
Beijing, 10084
<EMAIL_ADDRESS>
&Shanshan Tu
Beijing Key Laboratory of Trusted Computing
Beijing University of Technology,
Beijing, 100124
<EMAIL_ADDRESS>
&Obaid ur Rehman
Department of Electrical Engineering
Sarhad University of Science and IT
Pakistan, 25000
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

Water supplies are crucial for the development of living beings. However,
change in the hydrological process i.e. climate and land usage are the key
issues. Sustaining water level and accurate estimating for dynamic conditions
is a critical job for hydrologists, but predicting hydrological extremes is an
open issue. In this paper, we proposed two deep learning techniques and three
machine learning algorithms to predict stream flow, given the present climate
conditions. The results showed that the Recurrent Neural Network (RNN) or Long
Short-term Memory (LSTM), an artificial neural network based method,
outperform other conventional and machine-learning algorithms for predicting
stream flow. Furthermore, we analyzed that stream flow is directly affected by
precipitation, land usage, and temperature. These indexes are critical, which
can be used by hydrologists to identify the potential for stream flow. We make
the dataset publicly available (https://github.com/sadaqat007/Dataset) so that
others should be able to replicate and build upon the results published.

_K_ eywords Deep learning $\cdot$ Machine learning $\cdot$ Streamflow
estimation $\cdot$ Surface water hydrology $\cdot$ River/stream

##  1 Introduction

Figure 1: The overall framework of the proposed model, (a) CNN, (b) LSTM.
These two models are used to extract the feature vectors of amorphous
streamflow data and represent them to the feature space, which is further use
for predictive analysis

Water is an important element for all human beings as well as all living
things on the earth. Most of the water for domestic, industrial and
agriculture use is obtained from the streamflow. The researchers concluded
that streamflow has changed across the world [1, 2].

Streamflow is the main component of the hydrological cycle and it is essential
for the stability of our ecosystem. Streamflow has a vital impact on the
quality of water, living beings, and habitat in the stream [3]. Streamflow
controls the kinds of organisms that can survive in the stream, for example,
some organisms are susceptible to fast-flowing areas, while others need
stagnant pools. Moreover, it also controls the quantity of residue and
sediment carried by the stream. Large, rapidly flowing rivers are less
affected by pollution, whereas narrow small streams have less capability to
dilute and reduce wastes [4].

Automating and management of streamflow have the potential to significantly
reduce the withdrawals of water for irrigation purposes, as well as to improve
the quality of agriculture. The amount of water moving off the watershed into
the stream channel is highly related to the flow of a stream. It is affected
by the variation in weather, under different conditions such as rainstorm and
dry periods. Moreover, different seasons also fluctuate stream flow level, for
example, when evaporation is high and the ground is fertile thus removing
water from the ground due to which streamflow level decreases and vice versa.
Another serious problem of streamflow depletion is water withdrawals for
irrigation purposes, for example, industrial water withdrawals [5, 6].

Figure 2: The Location map of the Soan river basin and Distribution of Hydro-
meteorological station

Streamflow estimation and forecasting are challenging and essential for better
provision and management of water resources [7]. The streamflow prediction can
be useful for hydropower generation. Hydrologists have developed different
models and methods to estimate the streamflow [8]. Recently with the
development of the computer, programming attempts have been made to estimate
streamflow. The artificial neural network (ANN) has been applied successfully
for rainfall-runoff modeling, discharge prediction, sediment quantification,
and inflow forecasting [9, 10, 11]

In this paper, we applied deep learning techniques for high-level semantic
understanding for simulating the daily streamflow. Specifically, two deep
learning models: CNN and LSTM are adapted. This work is of great significance
and has aroused a wide spread of concern [12]. In the current era, deep
learning techniques have made significant development in fields such as
artificial intelligence [13, 14, 15] and natural language processing [16, 17].
These advancements of deep learning techniques can help the merging of
different technologies and inventions in other fields. Many applications are
witness that deep learning techniques have benchmark ability in feature
extractions and predictions [13, 18], which is free from conventional handy
feature extraction methods but carries out self-learning through data. This
can greatly leverage hydrological extremes, and alleviate the problem of
prediction for stream flow. In a nutshell, we don’t need to implement a
complex mathematical model for information learning and decision making, but
the deep learning models itself can automatically extract useful information
from the hydrological condition records by self-learning, and then conduct
streamflow prediction based on these information. This makes our model
efficient and robust compared to hand-crafted models. The complete overview of
the proposed framework is illustrated in Fig. 1

Thus, we selected a tropical basin, the Soan River basin, which is the left
bank tributary of the Indus River basin, for streamflow assessment. The
selected area is very crucial as it has high geographical importance. It is
the chief drinking water source to Islamabad city, Pakistan, and the key
source of irrigation water for Potohar plateau. In the past few decades, it
undergoes a significant climatic and land use [19, 20]. To our best knowledge,
we took the first effort to commence the trend analysis of hydrological
variables using deep learning algorithms for predicting streamflow change in
the Soan river basin. consequently, in this paper, we present and compare
results from two different deep learning models, such as CNN and RNN, and some
statistical machine learning models to assess the streamflow of the Soan River
basin. The main aims of this work are evaluating the performance of different
deep learning and machine learning techniques in predicting the stream flow in
Soan river basin (Punjab, Pakistan) based on the existing hydrological
variables, and to decide an efficient and robust model for this task.

##  2 Dataset Evaluation

###  2.1 Soan river basin.

In this paper, our main focus is on the Soan river basin. The people of
Islamabad (capital territory of Pakistan) obtain $95\%$ drinking water from a
simly dam, which is the chief reservoir for the Soan river basin. Soan river
initiating from Murree Mountains and runs through the Chirah hydrological
station, a major branch of Indus River and a significant hydrological unit of
Potohar plateau Pakistan. Figure 2 shows the visual description of the whole
site. After that, it is joined by different branches and at Dhoke Pathan
hydrological station it falls into the Indus River. Sub-basin1 is the area
above the Chirah hydrological station called Chirah sub-basin, and sub-basin 2
is the area above Dhoke Pathan hydrological station called Dhoke Pathan sub-
basin. The Soan basin covered area is $6842km^{2}$ having elevation in the
range of $265m-2274m$. Statistical data obtained from Pak govt. for the period
of $1983-2012$ shows that the annual precipitation of the basin is $1465mm$
and the mean annual temperature range from $8^{\circ}C$ to $18^{\circ}C$. The
Soan basin has both gentle and steep slope, stream flow in the basin mostly
comes from monsoon rainfall. About $60\%$ of the population associated is
rural and cultivation is the main source of economy. During last few years,
the population is greatly affected by the unpredictable variation in stream
flow of soan river basin, therefore, the implementation of a highly
predictable performance learning system is the need of the day to assess the
variation of streamflow for agriculture purpose [21, 22].

###  2.2 Data description.

Pakistan meteorological department (PMD) and Soil and Water Conservation
Research Institute (SAWCRI) provides the climatological data for this work.
The key hydroclimatic parameters, such as temperature and precipitation, were
obtained from six different meteorological stations. Water and Power
Development Authority (WAPDA) and Capital Development Authority (CDA) provide
the daily discharge data of two hydrological stations. Mwangi [23] et al.,
elaborate the discharge of Soan basin in detail and computed the runoff for
the basin, sub-basin1 and sub-basin2 using a conventional method. We obtain
statistical data for thirty years (1983-2012). The survey of Pakistan provides
the topographic information of the study area. Punjab bureau of statistics and
the Department of Agriculture (Marketing Wing) provided the agriculture
production information. The complete data descriptions of hydro-meteorological
variables are summarized in Table 1.

In this work, the predictors are precipitation and temperature, whereas
streamflow level is the response variable. These parameters are selected based
on a wealth of case studies. Specifically, streamflow level is highly
sensitive to precipitation. For similar reasons, as in Ficchi [24] et al., we
also include only two meteorological variables, such as precipitation and
temperature to provide a challenging scenario for learning algorithms to
leverage streamflow.

Table 1:  Hydrological Characteristics of the Soan basin and its sub-basins, 1983-2012. Area of Study |  Area/km2 |  Precipitation (mm/year) |  Potential Evapotranspiration (mm/year) |  Runoff (mm/year) |  Evaporation (mm/year)
---|---|---|---|---|---
Soan Basin |  6675 |  972.1 |  1382.5 |  241.8 |  730.3
Sub-basin 1 |  335 |  1649.4 |  934.5 |  695.8 |  953.6
Sub-basin 2 |  6340 |  988.3 |  1394.1 |  190.5 |  797.8

##  3 Models Design and Analysis.

In this study, we propose two different learning algorithms to extract
features from hydroclimatic variables and predict streamflow levels.
Specifically, CNN and RNN. As we discussed before, the key aim of selecting
these different types of learning algorithms is to validate streamflow
prediction on a variety of algorithms. The overall system frameworks are shown
in Figure 1. Detail descriptions and mathematical models for the individual
algorithm are as follow.

###  3.1 CNN-based prediction model.

Convolutional Neural Networks are the class of deep learning and is very
similar to ordinary Neural Networks. They are made up of neurons that have
learnable weights and biases. Each neuron receives some inputs, performs a dot
product and optionally follows it with a non-linearity. The whole network
expresses a single differentiable score function from the raw image pixels on
one end to class scores at the other [25, 26]. For the CNN-based prediction
model, a convolution neural network is introduced to extract features. The
input to the model is the daily temperature and rainfall information for the
upstream cities, and the output is the predicted river flow. The architecture
of the proposed CNN for streamflow prediction is as follows:

It consists of a convolutional layer with three various sizes of convolutional
kernels, followed by an average pooling layer and a fully-connected (fc) layer
with a supervision signal called softmax loss.

The function of the convolutional layer is to extract features map from the
input data matrix and different sizes convolutional kernels can learn various
low-level features. The pooling layer sufficed as a downsampling layer to
boost up system efficiency and robustness [13, 27]. Moreover, all these
features are fused in the fully-connected layer, which further passes them to
a supervision signal called softmax classifier for streamflow prediction. The
purpose of the softmax classifier is to compute the similarity measurement
between input feature vectors and finally display the streamflow probability.
Detail description of parameter tuning is given in the “Experimental Setting"
section.

For individual input data matrix, we describe it with a matrix
$X\in\mathbb{R}^{N\times D}$, as illustrated in Eq. 1. Whereas the $i$-th row
suggests the $i$-th day, and $D$ represents the dimension of the input
feature. Mathematically it can be represented as:

$X=\left[\begin{array}[]{c}x_{1}\\\ x_{2}\\\ \vdots\\\ x_{N}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}a_{1,1}&a_{1,2}&\cdots&a_{1,D}\\\
a_{2,1}&a_{2,2}&\cdots&a_{2,D}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{N,1}&a_{N,2}&\cdots&a_{N,D}\\\ \end{array}\right]$ (1)

Generally, let $X_{i:j}$ refer to the matrix which feature vectors from the
$i$-th day to the $j$-th day, mathematically:

$X_{i:j}=\left[\begin{array}[]{c}x_{i}\\\ x_{i+1}\\\ \vdots\\\ x_{j}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}a_{i,1}&a_{i,2}&\cdots&a_{i,D}\\\
a_{i+1,1}&a_{i+1,2}&\cdots&a_{i+1,D}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{j,1}&a_{j,2}&\cdots&a_{j,D}\\\ \end{array}\right]$ (2)

Each convolution layer is made up of several convolution kernels of various
sizes, and every size further contains many convolution kernels. The width of
the convolution kernel is aligned with the width of the input data matrix. For
example, if the height of $k$-th convolution kernel is $H$, then it can be
described as $W^{k}\in\mathbb{R}^{H\times D}$, mathematically:

$W^{k}=\left[\begin{array}[]{cccc}w^{k}_{1,1}&w^{k}_{1,2}&\cdots&w^{k}_{1,D}\\\
w^{k}_{2,1}&w^{k}_{2,2}&\cdots&w^{k}_{2,D}\\\ \vdots&\vdots&\ddots&\vdots\\\
w^{k}_{H,1}&w^{k}_{H,2}&\cdots&w^{k}_{H,D}\\\ \end{array}\right]$ (3)

The function of convolution operation is the features extraction for local
region elements of the input data matrix. For instance, when $w^{k}_{1,1}$ and
$a_{1,1}$ co-occur, then the feature $c^{k}_{1}$ extracted from $X_{1:H}$ by
the convolutional kernel can be expressed as:

$c^{k}_{1}=f(\sum_{i=1}^{H}\sum_{j=1}^{D}w^{k}_{i,j}\cdot
a_{i,j}+b^{k}_{i,j}),$ (4)

weight $w^{k}_{i,j}$ represents the significance of the $j$-th value in the
$i$-th row vector, $b^{k}_{i,j}$ denote the bias term and $f$ is a nonlinear
function, as described in [13]. The $ReLu$ function is described as nonlinear
function, which is mathematically illutrated as:

$z=ReLu(i)=max(0,i).$ (5)

The convolution process is obtained by doing the convolution between the
kernels $W^{k}$ and input data matrix $X$ with a certain step of $T_{c}$, and
computes the features of individual local region. In conclusion, the feature
extracted by the convolution kernel $W^{k}$ can be expressed as:

$C^{k}=[c^{k}_{1},c^{k}_{2},...,c^{k}_{\frac{N-H+1}{T_{c}}}]^{\top}.$ (6)

The pooling layer served to lessen the amount of neural network parameter and
maintain the overall data distribution, which eventually prevent the system
from over-fitting problem and enhance system robustness [13, 27]. Both the
pooling and convolution operations have similar function, with a difference
that pooling operatikon computes the maximum or average value of the local
area. Here, after obtaining the feature $C^{k}$ from convolution operation we
apply a max pooling operation. Let, $H_{p}$ represent the height of a pooling
kernel and $T_{p}$ represent its step size, then mathematically we compute the
output as:

$M^{k}=[m^{k}_{1},m^{k}_{2},...,m^{k}_{N_{p}}]^{\top},$ (7)

where

$m^{k}_{i}=max(c^{k}_{i},c^{k}_{i+1},...,c^{k}_{i+H_{p}-1}).$ (8)

$N_{p}=\frac{\frac{N-H+1}{T_{c}}-H_{p}+1}{T_{p}},$ (9)

The aforementioned process described that a single convolution kernel $W^{k}$
produces single feature map $M^{k}$. Finally, after the successful completion
of convolution and pooling operations, the entire extracted features are
combined in an end-to-end fashion to achieve a single feature vector of the
whole input data matrix, which can be mathematically expressed as:

$F^{\top}=[F^{\top}_{1};F^{\top}_{2};...;F^{\top}_{l}],$ (10)

where $F_{i}=M^{i}$, $l$ represents the total number of features.

Finally, for high-level features extraction fully-connected layer is used.
Let, $W_{F}$ represent the weight matrix, in order to calculate the weighted
sum of the individual element to get the final feature representation of the
input data matrix $S$, mathematically:

$y=W_{F}\cdot F+b_{f},$ (11)

where $W_{F}$ represents the learned weight matrix and $b_{f}$ represents the
bias term, however, the weight matrix $W_{F}$ highlights the significance of
every individual feature.

The parameters of the network are updated using backpropagation algorithm in
the training process, and the loss function is computed through Smooth-L1 loss
called Huber loss:

$loss(\hat{y},y)=\frac{1}{n}\sum_{i}z_{i},$ (12)

where $z_{i}$ can be calculated as:

$z_{i}=\begin{cases}0.5(\hat{y}-y)^{2},&\text{if}~{}|\hat{y}-y|<1\\\
|\hat{y}-y|-0.5,&\text{otherwise}.\end{cases}$ (13)

The main choice of selecting Smooth-$L_{1}$ loss is to diminish data
sensibility, as outliers greatly affect $L_{2}$ loss. In this paper, we used
the minimization of the $LOSS$ function over a batch size number of samples
for CNN training. The optimization algorithm is Ada-delta with $1.0$ initial
learning rate.

###  3.2 LSTM-based prediction model

Recurrent Neural Network (RNN) is an artificial neural network architecture,
which is very suitable for sequential signals modeling. Mathematically, we
described the RNN model in Eq. $(8)$, which can deal with randomly size
sequential signals. Though, the system is not able to handle long-range
problem efficiently due to the gradient vanishing problem [28]. LSTM model
[29], a variant of RNN efficiently solve this problem through rich designed
unit nodes. The key architecture refinement of LSTM is the hidden layer unit,
which contained four fundamental parts: a cell, an input gate, an output gate
and a forget gate. The basic LSTRM mathematical model is illustrated as
follows:

$\displaystyle\left\\{\begin{array}[]{l}I_{t}=\sigma(W_{i}\cdot[h_{t-1},x_{t}]+b_{i}),\\\
F_{t}=\sigma(W_{f}\cdot[h_{t-1},x_{t}]+b_{f}),\\\ C_{t}=F_{t}\cdot
C_{t-1}+I_{t}\cdot\tanh(W_{c}\cdot[h_{t-1},x_{t}]+b_{c}),\\\
O_{t}=\sigma(W_{o}\cdot[h_{t-1},x_{t}]+b_{o}),\\\
h_{t}=O_{t}\cdot\tanh(C_{t}).\\\ \end{array}\right.$ (19)

where $I_{t}$ represents the input gate, which controls the quantity of novel
information stored in the memory cell. The function of forget gate $F_{t}$ is
to discard redundant stored information from the memory cell. Hence, memory
cell $C_{t}$ is a summation of upcoming and previous information modulated by
the input gate and the forget gate, respectively. The output gate $O_{t}$
allows the memory cell to have an effect on the current hidden state and
output or block its influence. Specifically, we represent the transfer
function of LSTM units by $f_{LSTM}(\ast)$. Notice that, while calculating the
output at time step $t$, the information used is based on the input vector at
time step $t$, also includes the information stored in the cells at the
previous time step $t-1$. Hence, the output at a given time steps $t$ can be
described as:

$y_{t}=f_{LSTM}(x_{t}\mid x_{1},x_{2},...,x_{t-1}).$ (20)

As we have mentioned before, the input matrix $X$ can be regarded as a
sequential signal and the $i$-th row $x_{i}$ can be viewed as the signal at
time step $i$. In order to model this task on LSTM, we input the $i$-th row of
$X$ at time step $i$. For individual sentence $S$, we can describe it with a
matrix $S\in\mathbb{R}^{N\times D}$ as described in Eq. 1.

Overall, an RNN contains several network layers, each of which with multiple
LSTM units. Let, $n_{j}$ represent the number of LSTM units $U^{j}$ of $j$-th
hidden layer, thus the units of $j$-th layer can be described as:

$U^{j}=\\{u^{j}_{1},u^{j}_{2},\cdots,u^{j}_{n_{j}}\\}.$ (21)

In LSTM first hidden layer, the input of each unit $u^{1}_{i}$ at time step
$t$ is the weighted sum of the elements in $x_{t}$, mathematically as:

$u^{1}_{i,t}=\sum^{d}_{k=1}w^{1}_{i,k}\cdot a_{j,k}+b^{1}_{i,t},$ (22)

where $w^{1}_{.}$ represents the learned weights and $b^{1}.$ represents the
biases. Hence, the output value of $u^{1}_{i}$ at time step $t$ is:

$o^{1}_{i,t}=f_{LSTM}(u^{1}_{i,t})=f_{LSTM}(\sum_{k=1}^{d}w^{1}_{i,k}\cdot
a_{j,k}+b^{1}_{i,t}).$ (23)

The vector $Out^{j}_{t}$ is used to represent the output of the $j$-th hidden
layer at time step $t$, and each element in $Out^{j}_{t}$ represents the
output value of each unit in the $j$-th hidden layer at time step $t$,
mathematically:

$Out^{j}_{t}=f_{LSTM}(x_{t})=[o^{j}_{1,t},o^{j}_{2,t},...,o^{j}_{n_{j},t}].$
(24)

Literature showed that increasing the number of neural network layers,
increasing the ability of system to extract features [30]. Therefore, in this
paper, the network is stack with multiple layers of LSTM units. The transfer
matrix is used to connect the adjacent hidden layers. For instance, the
transfer matrix between $l$-th layer and $(l+1)$-th layer can be described as
a matrix $W^{l}\in\mathbb{R}^{{n_{l}}\times{n_{l+1}}}$:

$W^{l}=\left[\begin{array}[]{cccc}w^{l}_{1,1}&w^{l}_{1,2}&\cdots&w^{l}_{1,n_{l+1}}\\\
w^{l}_{2,1}&w^{l}_{2,2}&\cdots&w^{l}_{2,n_{l+1}}\\\
\vdots&\vdots&\ddots&\vdots\\\
w^{l}_{n_{l},1}&w^{l}_{n_{l},2}&\cdots&w^{l}_{n_{l},n_{l+1}}\\\
\end{array}\right]$ (25)

Then the input of each unit $u^{l}_{i}$ in the $l$-th hidden layer at time
step $t$ is the weighted sum of the output values of the units in the previous
layer, mathematically:

$u^{l}_{i,t}=Out^{l-1}_{t}\cdot W^{l,m}=\sum_{k=1}^{n_{l-1}}w^{l}_{i,k}\cdot
o^{l-1}_{k,t}+b^{l}_{i,t}.$ (26)

Furthermore, the output of the $l$-th layer at time interval $t$ is:

$\displaystyle
Out^{l}_{t}=f_{LSTM}(Out^{l}_{t})=[o^{l}_{1,t},o^{l}_{2,t},...,o^{l}_{n_{l},t}],$
(27) $\displaystyle
o^{l}_{i,t}=f_{LSTM}(u^{l}_{i,t})=f_{LSTM}(\sum_{k=1}^{n_{l+1}}w^{l}_{i,k}\cdot
o^{l}_{k,t}+b^{l}_{i,t}).$

The aforementioned Eq. 27 illustrated that the output at time step $t$ is not
only based on the input vector of the current time $x_{t}$, but also the
information stored in the cells at the previous $(t-1)$ intervals. Thus,
$l$-th hidden layer output at time interval $t$ can be obtained as a summation
of all former $t$ intervals, such as information fusion of former $t$
intervals $\\{x_{1},x_{2},...,x_{t}\\}$. Moreover, to compute the probability
distribution of the $(t+1)$ day, a softmax layer is added after all the hidden
layers of the LSTM. Particularly, we introduce the Prediction Weight (PW) as
matrix $W_{P}\in\mathbb{R}^{{n_{l}}\times N}$, mathematically as:

$W_{P}=\left[\begin{array}[]{cccc}w^{p}_{1,1}&w^{p}_{1,2}&\cdots&w^{p}_{1,N}\\\
w^{p}_{2,1}&w^{p}_{2,2}&\cdots&w^{p}_{2,N}\\\ \vdots&\vdots&\ddots&\vdots\\\
w^{p}_{n_{3},1}&w^{p}_{n_{3},2}&\cdots&w^{p}_{n_{l},N}\\\ \end{array}\right]$
(28)

where $N$ represents the total number of days in the dictionary $D$. Then the
learned matrix $W_{P}$ is used to compute the score for individual day in the
dictionary $D$ for a time span of 1983-2012. Mathematically, we can define as:

$y_{i}=\sum^{n_{l}}_{k=1}w^{p}_{k,i}\cdot o^{l}_{i,t}+b^{p}_{i,t},$ (29)

where $W_{P}$ represents the learned weight and $b^{p}$ represents the bias
term, the values in weight matrix $W_{P}$ reflect the significance of
individual feature in $o^{l}$. Then we can describe the dimension of output
vector $y$ as $N$.

Both CNN and LSTM based prediction models use the same Loss function and
training algorithm.

##  4 Experimental Results and Discussion

Figure 3: Precipitation and stream flow statistics of the Soan sub-basin1 for
the year 1983-2012 Figure 4: Precipitation and streamflow statistics of the
Soan sub-basin2 for the year1983-2012

###  4.1 Experiment Setting

In this paper, we performed all the experiments on two deep learning models,
such as CNN and LSTM, and three machine learning algorithms, such as RF, GBR,
and LR. For the CNN model, the input data is fed to a convolution layer with
three different kernel height: $3,5,7$. The width of kernels is equal to the
feature number. Each kernel initially generates 100 channels. After the
convolution operation, max pooling is applied to each channel. Then the
outputs are concatenated together to the last layer, which is a fully-
connected layer with a dimension of $300\times 1$. The output of which is the
predicted flow.

For the LSTM model, we use a bidirectional LSTM with cell number equal to the
features number and 300 hidden sizes. Moreover, the max-pooling layer and
fully-connected layer follow the same architecture as CNN.

As for the traditional machine learning methods, such as Logistic Regression
(LR) [31, 32], Gradient Boosting Tree (GBR) [33] and Random Forest (RF) [34],
we used the default setting from package scikit-learn.

###  4.2 Dataset Split

We divided the dataset into three subsets, such as training set, validation
set, and testing set, with a proportion of $7:1:2$. The final result is
obtained using a 10 fold cross validation with different random dataset
division. Normalization is applied to each feature in the preprocessing.

Mini-batch is used in the training process of CNN and LSTM, with 90 batch
size. The optimization algorithm is Ada-delta with $1.0$ initial learning
rate.

###  4.3 Evaluation Method

The index we adapt to evaluate results is the relative error, which is defined
as:

$re=\frac{1}{n}\sum_{i}\frac{|predict-real|}{real}.$ (30)

Notice that we have normalized the raw data in preprocessing, so before
calculating relative error, the output will be recovered to the initial
distribution.

###  4.4 Hydro-climatic Trends in Soan Basin

Fig. 3, 4 shows the hydro-climatic trends in sub-basin 1 and sub-basin 2 (all
units in mm). It can be observed that precipitation and streamflow of both
sub-basin have a strong correlation. The catchment area of the sub-basin 1 is
smaller as compared to the sub-basin 2, therefore the streamflow of the sub-
basin 1 is more as compared to sub-basin 2. The sub-basin 1 has hilly areas
and due to high mountains and slope differences, most of the precipitation is
converted into streamflow and also has small water loses. On the other hand,
sub-basin 2 has less precipitation as compared to sub-basin 1, which
eventually causes less streamflow in sub-basin 2. There is more water loss in
sub-basin 2 and has high agriculture and water consumption. Moreover, the
water storages and small check dams in sub-basin 2 can cause a reduction in
the streamflow. While, on the other hand, sub-basin 1 has no small check dams
or retention, which can stop the water and topography is favorable for high
streamflow and less water loss.

###  4.5 Predictive performance of different models.

Supervised learning algorithms are the major benchmark paradigms found in many
applications, such as classification, detection and hydrological modeling [35,
36, 37, 38]. Supervised learning algorithms contain three fundamental parts:
an input object (vector), a processing unit (hidden layers) and an output
value (supervision signal). These learning algorithms produce an inferred
function by evaluating the training data, which can be further employed for
representing new examples. An optimal setting allows these algorithms to
describe the class labels for hidden examples, which help the learning
algorithms to generalize for unseen conditions from the training data. Let,
given a set of $\aleph$ training examples $(x_{i},y_{j}),...,(x_{N})$ such
that $x_{i}$ represent the feature vector of $i^{th}$ example and $y_{j}$ is
its class label, a supervised learning algorithm search for a function $g:X\to
Y$, where $X$ and $Y$ is the input and output space respectively. The function
$g$ is an element of the hypothesis space of possible functions
$\mathfrak{G}$. Under certain circumstances $g$ is represented as a scoring
function $f:X\times Y\to{\mathbb{R}}$, such that, it return the value of $y$
that gives the maximum score: $g(x)={\underset{y}{\arg\max}}f(x,y)$ . Suppose,
$F$ represents the scoring functions.

Table 2: Relative Error of different Deep learning and Machine learning methods. Model | Mean Relative Error(%) | Standard Deviation
---|---|---
LR | 28.0377 | 1.9769
GBR | 27.4681 | 1.8108
RF | 31.5304 | 1.9430
CNN | 26.8145 | 1.6070
LSTM | 25.0984 | 1.9182

Figure 5: Prediction results using LSTM. The horizontal axis shows the real
stream flow data, whereas the vertical axis shows the predicted results
obtained through LSTM

Figure 6: Prediction results using CNN. The horizontal axis shows the real
streamflow data, whereas the vertical axis shows the predicted results
obtained through CNN

Figure 7: Prediction error changes with days

Even though, $G$ and $F$ can be any space of functions, many learning
algorithms such as CNN and LSTM are probabilistic models where $g$ takes the
form of a joint probability model $g(x,y)=P(x,y)$, and logistic regression is
a conditional probability model where $f$ takes the form of a conditional
probability model $g(x)=P(y|x)$.

In this paper, we evaluated five different learning models for measuring
predictive performance – CNN, LSTM, RF, GBR, and LR - to predict the stream
flow based on hydrological variables. These different learning models were
chosen to test the performance of a wide variety of learning algorithms on
streamflow prediction problem. The inclusion of both statistical models, such
as LR, GBR, and RF, and deep learning models, such as LSTM and CNN, help to
ensure that the predictive performance is tested on a wide variety of learning
algorithms. Statistical models are used in a variety of fields, such as
medical disease diagnosis, machine learning and social sciences for
prediction. Moreover, these models itself simply model the probability of
output in terms of input and do not perform statistical classification. CNN
and LSTM models can provide robust predictions, due to the transformations of
the input space in the inner layers [39]. The outcome of each algorithm was
summarized using Mean Relative error (MRE) of successful prediction (error
rate %), as shown in Table 2. Moreover, the advantage of deep learning models
(CNN and LSTM) over statistical models is that deep learning models continue
to improve efficiency as the number of datasets grows. Whereas, statistical
models plateau at a certain level of performance when you add more examples
and training data to the network.

Table 2 shows the parameter settings and numerical results obtained using
different algorithms including the proposed algorithms. It is clear from the
results that the LSTM surpasses other benchmark algorithms in solving
streamflow prediction. It demonstrates that the Random Forest (RF) approach
achieved the highest MRE (31.53% ) on the streamflow prediction dataset.
Followed by Improved Linear Regression (LR) with mean relative error (MRE) of
28.03%. The Gradient Boosted Regression (GBR) algorithm and CNN have similar
MRE of 27.46% and 26.81% respectively (the overall difference was no more than
0.6%). The lowest MSE is achieved by the LSTM method, which is 25.09%.
Furthermore, LSTM also shows high predictive accuracy despite having a small
error, as shown in Fig. 6. It shows the real stream flow plotted against the
predicted flow by LSTM with a $45{\circ}$ reference line. Followed by CNN,
which also shows a good predictive accuracy as shown in Fig. 6. Hence, we
concluded that the LSTM algorithm with the lowest MRE achieves the highest
prediction rate, which needs less training time to obtain an optimum solution.

The accuracy of forecasting traffic is related to the forecast duration. For
example, according to the data of the first five days or the first 10 days or
the first month, the accuracy of the next day’s traffic forecast is different.
Therefore, the abscissa of this graph is how many days before we forecast the
next day. The ordinate is the accuracy of the prediction. The main purpose is
to find the best forecast duration. Based on the data in this Fig. 7, we find
that when the predicted number of days is 7, the predicted error is the
smallest. So our model was finally determined to predict river flow of the
next day based on data from the previous 7 days.

##  5 Conclusions

Two deep learning techniques and three machine learning algorithms are
proposed to predict stream flow, given the present climate conditions. These
different learning models were chosen to evaluate the performance of a wide
variety of learning algorithms on streamflow prediction problem. The inclusion
of both statistical models, such as LR, GBR, and RF, and deep learning models,
such as LSTM and CNN, help to ensure that the predictive performance is tested
on a wide variety of learning algorithms. The results showed that the LSTM, an
artificial neural network based method, outperforms other conventional and
machine-learning algorithms for predicting stream flow.

To our best knowledge, this is the first work to introduce five different deep
learning and machine learning algorithms for streamflow prediction problem.
Furthermore, we believe that the proposed predictive algorithms serve as a
reference tool for researchers and engineers to help the design and
implementation of improving evaluation protocols.

## Acknowledgments

This work was partially supported by the National Natural Science Foundation
of China (No. 61801008, U1705261, U1405254, U1536115, U1536207), National Key
R&D Program of China (No. 2018YFB0803600, 2016YFB0801301), Beijing Natural
Science Foundation National (No. L172049), Scientific Research Common Program
of Beijing Municipal Commission of Education (No. KM201910005025).

## References

  * [1] David Labat, Yves Goddéris, Jean Luc Probst, and Jean Loup Guyot.  Evidence for global runoff increase related to climate warming.  Advances in Water Resources, 27(6):631–642, 2004.
  * [2] Nicola Gedney, PM Cox, RA Betts, O Boucher, C Huntingford, and PA Stott.  Detection of a direct carbon dioxide effect in continental river runoff records.  Nature, 439(7078):835, 2006.
  * [3] TR Petty and P Dhingra.  Streamflow hydrology estimate using machine learning (shem).  JAWRA Journal of the American Water Resources Association, 54(1):55–68, 2018.
  * [4] Gang Zhao, Huilin Gao, and Lan Cuo.  Effects of urbanization and climate change on peak flows over the san antonio river basin, texas.  Journal of Hydrometeorology, 17(9):2371–2389, 2016.
  * [5] Zhi Li and Jiming Jin.  Evaluating climate change impacts on streamflow variability based on a multisite multivariate gcm downscaling method in the jing river of china.  Hydrology and Earth System Sciences, 21(11):5531–5546, 2017.
  * [6] Bibi S Naz, Shih-Chieh Kao, Moetasim Ashfaq, Huilin Gao, Deeksha Rastogi, and Sudershan Gangrade.  Effects of climate change on streamflow extremes and implications for reservoir inflow in the united states.  Journal of Hydrology, 556:359–370, 2018.
  * [7] Manish Kumar Goyal, CSP Ojha, RD Singh, PK Swamee, et al.  Application of artificial neural network, fuzzy logic and decision tree algorithms for modelling of streamflow at kasol in india.  Water Science and Technology, 68(12):2521–2526, 2013.
  * [8] Amir Mosavi, Pinar Ozturk, Shahab Shamshirband, Hai Thanh Nguyen, and Kwok-wing Chau.  Flood prediction using machine learning.  Literature Review, Engineering Applications of Computational Fluid Mechanics, 2018.
  * [9] Caihong Hu, Qiang Wu, Hui Li, Shengqi Jian, Nan Li, and Zhengzheng Lou.  Deep learning with a long short-term memory networks approach for rainfall-runoff simulation.  Water, 10(11):1543, 2018.
  * [10] Dongwoo Jang, Hyoseon Park, and Gyewoon Choi.  Estimation of leakage ratio using principal component analysis and artificial neural network in water distribution systems.  Sustainability, 10(3):750, 2018.
  * [11] Chen Liang, Hongqing Li, Mingjun Lei, and Qingyun Du.  Dongting lake water level forecast and its relationship with the three gorges dam based on a long short-term memory network.  Water, 10(10):1389, 2018.
  * [12] Sills J. et al.  Nextgen voices: Unique identities[j].  Science, 363(6428):702–702, 2019.
  * [13] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [14] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [15] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
  * [16] Nal Kalchbrenner, Edward Grefenstette, and Phil Blunsom.  A convolutional neural network for modelling sentences.  arXiv preprint arXiv:1404.2188, 2014.
  * [17] Yoon Kim.  Convolutional neural networks for sentence classification.  arXiv preprint arXiv:1408.5882, 2014.
  * [18] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In European conference on computer vision, pages 818–833. Springer, 2014.
  * [19] Winston Yu, Yi-Chen Yang, Andre Savitsky, Donald Alford, Casey Brown, James Wescoat, Dario Debowicz, and Sherman Robinson.  The Indus basin of Pakistan: The impacts of climate risks on water and agriculture.  The World Bank, 2013.
  * [20] Muhammad Shahid, Hamza Farooq Gabriel, Amjad Nabi, Sajjad Haider, AS Khan, and AMS Shah.  Evaluation of development and land use change effects on rainfall-runoff and runoff-sediment relations of catchment area of simly lake pakistan.  Life Science Journal, 11(7s), 2014.
  * [21] A Ashfaq, M Ashraf, and A Bahzad.  Spatial and temporal assessment of groundwater behaviour in the soan basin of pakistan.  University of Engineering and Technology Taxila. Technical Journal, 19(1):12, 2014.
  * [22] Muhammad Shahid, Zhentao Cong, and Danwu Zhang.  Understanding the impacts of climate change and human activities on streamflow: a case study of the soan river basin, pakistan.  Theoretical and Applied Climatology, 134(1-2):205–219, 2018.
  * [23] Hosea M Mwangi, Stefan Julich, Sopan D Patil, Morag A McDonald, and Karl-Heinz Feger.  Relative contribution of land use change and climate variability on discharge of upper mara river, kenya.  Journal of Hydrology: Regional Studies, 5:244–260, 2016.
  * [24] Andrea Ficchi, Luciano Raso, Pierre-Olivier Malaterre, David Dorchies, Maxime Jay-Allemand, Francesca Pianosi, Peter-Jules van Overloop, and Guillaume Thirel.  Short term reservoirs operation on the seine river: Performance analysis of tree-based model predictive control.  2014\.
  * [25] Sadaqat Rehman, Shanshan Tu, Obaid Rehman, Yongfeng Huang, Chathura Magurawalage, Chin-Chen Chang, et al.  Optimization of cnn through novel training strategy for visual classification problems.  Entropy, 20(4):290, 2018.
  * [26] Sadaqat ur Rehman, Shanshan Tu, Yongfeng Huang, Guojie Liu, et al.  Csfl: A novel unsupervised convolution neural network approach for visual pattern classification.  AI Communications, 30(5):311–324, 2017.
  * [27] Y-Lan Boureau, Jean Ponce, and Yann LeCun.  A theoretical analysis of feature pooling in visual recognition.  In Proceedings of the 27th international conference on machine learning (ICML-10), pages 111–118, 2010.
  * [28] Sepp Hochreiter.  The vanishing gradient problem during learning recurrent neural nets and problem solutions.  International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems, 6(02):107–116, 1998.
  * [29] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [30] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks (2013).  arXiv preprint arXiv:1311.2901, 2013.
  * [31] Johan AK Suykens and Joos Vandewalle.  Least squares support vector machine classifiers.  Neural processing letters, 9(3):293–300, 1999.
  * [32] Frank E Harrell.  Ordinal logistic regression.  In Regression modeling strategies, pages 311–325. Springer, 2015\.
  * [33] Jerome H Friedman.  Stochastic gradient boosting.  Computational Statistics & Data Analysis, 38(4):367–378, 2002\.
  * [34] Andy Liaw, Matthew Wiener, et al.  Classification and regression by randomforest.  R news, 2(3):18–22, 2002.
  * [35] Roshanak Nateghi, Seth D Guikema, and Steven M Quiring.  Comparison and validation of statistical methods for predicting power outage durations in the event of hurricanes.  Risk Analysis: An International Journal, 31(12):1897–1906, 2011\.
  * [36] Roshanak Nateghi, Jeremy D Bricker, Seth D Guikema, and Akane Bessho.  Statistical analysis of the effectiveness of seawalls and coastal forests in mitigating tsunami impacts in iwate and miyagi prefectures.  PloS one, 11(8):e0158375, 2016.
  * [37] B Grizzetti, A Pistocchi, C Liquete, A Udias, F Bouraoui, and W Van De Bund.  Human pressures and ecological status of european rivers.  Scientific reports, 7(1):205, 2017.
  * [38] Kazuya Nishina, Mirai Watanabe, Masami K Koshikawa, Takejiro Takamatsu, Yu Morino, Tatsuya Nagashima, Kunika Soma, and Seiji Hayashi.  Varying sensitivity of mountainous streamwater base-flow no 3- concentrations to n deposition in the northern suburbs of tokyo.  Scientific reports, 7(1):7701, 2017.
  * [39] Stefano Lancini, Marco Lazzari, Alberto Masera, and Paolo Salvaneschi.  Diagnosing ancient monuments with expert software.  Structural engineering international, 7(4):288–291, 1997.