# Use of the Hayami diffusive wave equation to model the relationship infected–recoveries–deaths of Covid-19 pandemic

Roger Moussa, Samer Majdalani

► To cite this version:

Roger Moussa, Samer Majdalani. Use of the Hayami diffusive wave equation to model the relationship infected–recoveries–deaths of Covid-19 pandemic. Epidemiology and Infection, 2021, 149, 10.1017/s0950268821001011 hal-03299735

HAL Id: hal-03299735

https://hal.inrae.fr/hal-03299735

Submitted on 26 Jul 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

## Original Paper

**Cite this article:** Moussa R, Majdalani S (2021). Use of the Hayami diffusive wave equation to model the relationship infected-recoveries-deaths of Covid-19 pandemic. *Epidemiology and Infection* **149**, e138, 1–13. https://doi.org/10.1017/S0950268821001011

Received: 28 June 2020

Revised: 18 February 2021

Accepted: 17 March 2021

**Key words:**

Coronavirus Covid-19; diffusive wave equation; epidemiology; Hayami kernel function; hydraulic; hydrology; SIRD compartmental models; unit hydrograph

**Author for correspondence:**

Roger Moussa, E-mail: roger.moussa@inrae.fr

# Use of the Hayami diffusive wave equation to model the relationship infected-recoveries-deaths of Covid-19 pandemic

Roger Moussa¹ and Samer Majdalani²

¹LISAH, Univ. Montpellier, INRAE, IRD, Montpellier SupAgro, Montpellier, France and ²HSM, CNRS, IRD, Univ. Montpellier, Montpellier, France

## Abstract

Susceptible *S*-Infected *I*-Recovered *R*-Death *D* (*SIRD*) compartmental models are often used for modelling of infectious diseases. On the basis of the analogy between *SIRD* and compartmental models in hydrology, this study makes mathematical formulations developed in hydrology available for modelling in epidemiology. We adapt the Hayami model solution of the diffusive wave equation generally used in hydrological modelling to compartmental *I*–*R*–*D* models in epidemiology by simulating the relationships between the number of infectious *I*(*t*), the number of recoveries *R*(*t*) and the number of deaths *D*(*t*). The Hayami model is easy-to-use, robust and parsimonious. We compare the empirical one-parameter exponential model usually used in *SIRD* models to the two-parameter Hayami model. Applications were implemented on the recent Covid-19 pandemic. The application on data from 24 countries shows that both models give comparable performances for modelling the *I*–*D* relationship. However, for modelling the *I*–*R* relationship and the active cases, the exponential model gives fair performances whereas the Hayami model substantially improves the model performances. The Hayami model also presents the advantage that its parameters can be easily estimated from the analysis of the data distributions of *I*(*t*), *R*(*t*) and *D*(*t*). The Hayami model is parsimonious with only two parameters which are useful to compare the temporal evolution of recoveries and deaths in different countries based on different contamination rates and recoveries strategies. This study highlights the interest of knowledge transfer between different scientific disciplines in order to model different processes.

## Introduction

Since the pioneer study of [1] in epidemiology, *SIRD* (Susceptible *S*-Infected *I*-Recovered *R*-Death *D*) compartmental models are often used for modelling of infectious diseases, and are now used for the recent Covid-19 pandemic in China [2–4]. *SIRD* models can be used at the patch scale or over a large area [5–10], and can be coupled to probabilistic approaches [11, 12] and various processes such as population mobility [13–16]. When used at the scale of a patch, or of a country when the latter is considered as one homogeneous patch, *SIRD* models can be viewed similar to lumped hydrological models representing the links between several water cycle compartments [17] with rain as input (equivalent to the input *S*) and river flow, infiltration and evaporation as outputs (equivalent to the outputs *R* and *D*). Different hydrological modelling approaches were developed on the basis of an interaction between two categories of ‘functions’ called ‘production function’ and ‘transfer function’ [17, 18]. The production function separates the rain into runoff, infiltration and evaporation, and therefore distributes the mass balance between different compartments. Based on mass conservation, the transfer function is considered as a filter which transforms an input signal into an output signal, and calculates the time distribution of the production function outputs. In epidemiology, we can define by analogy a ‘production function’ to represent the relationship *S*–*I* which calculates the part of *S* contributing to *I*, and a ‘transfer function’ to represent the relationships *I*–*R*–*D*.

For the Covid-19 pandemic, the daily number of infected cases *I*(*t*), recoveries *R*(*t*) and deaths *D*(*t*) are freely available on websites of governmental institutions [19–21]. From graphics, we observe that the relationships *I*–*R* and *I*–*D* have similar shapes to those obtained from hydrologic transfer functions. First, mass conservation is verified because the total number of *I*(*t*) is equal to the total numbers of *R*(*t*) and *D*(*t*). Second, the output signals *R*(*t*) and *D*(*t*) can be derived from the input signal *I*(*t*) using a mathematical formulation identical to the transfer function describing the physical advection-dispersion processes: a lag time translation (for the advection process) and an attenuation of the peak (for the dispersion process). Although on the basis of analogical properties, transfer functions have been successfully used in other disciplines such as astronomy, geophysics, soil mechanics, meteorology, oceanography, traffic simulation or biological flows [22, 23], to our knowledge, no applications of the transfer

© The Author(s), 2021. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.

CAMBRIDGE
UNIVERSITY PRESS
functions were conducted for modelling I-R-D relationships in compartmental models in epidemiology.

This study aims to adapt transfer functions used in hydrology for modelling I-R-D relationships in epidemiology. I-R-D models are generally described by partial differential equations [24–26] where parsimonious first-order kinetics remains to be largely used because of their simplicity, robustness and their low number of parameters [27, 28]. In hydrology, under the hypotheses that the system is linear and time invariant, simplified versions of the transfer function were developed like the ‘unit hydrograph’, that is equivalent to a convolution with a kernel function representing the response to a Dirac input [29]. Furthermore, a linear reservoir is a one-parameter kernel function (or unit hydrograph) with an exponential decrease equivalent to first-order kinetics used in epidemiology (dR/dt = γI or dD/dt = δI with γ and δ as parameters). However, one-parameter unit hydrograph is not adapted to represent processes such as flood routing in channels where at least two-parameter unit hydrographs are needed to describe both advection and dispersion processes.

In epidemiology, first-order kinetics remains to be the simplest model [2, 16], and a question arises whether one-parameter first-order kinetics is sufficient to model at the country scale the I-R-D relationships for Covid-19, and whether simulations could be improved by using other mathematical formulations such as the two-parameter kernel functions. Various unit hydrographs have been developed [29], but one is of interest for possible applications in epidemiology: the Hayami model [30] used as a unit hydrograph resolution of the diffusive wave equation with parameters having a physical interpretation [31–34]. For applications in epidemiology, questions arise on the adaptability of the Hayami model to simulate the relationships I-R-D, the interpretation of model parameters and how could these parameters be used to compare different case studies.

This study attempts to make the Hayami solution of the diffusive wave equation available for modelling the I-R-D relationships in epidemiology. Applications were implemented on the actual 2020 Covid-19 pandemic. First, we present the I-R-D model, the Hayami model, its adaptation for compartmental models in epidemiology, and an application example to simulate R and D (given I) in the case of China where the pandemic is ending. Then, applications were implemented on data from 24 countries all round the world, with different levels of evolution of the pandemic, contamination rates and recoveries strategies. Finally, we show how the model parameters can be used to compare recoveries and deaths in different countries, and discuss the usefulness and the limitations of the model. The Supplementary material presents the theory of the unit hydrograph, and shows a comparison between the observed and calculated R(t), D(t) an A(t) for 24 countries.

## The I-R-D model

The I-R-D model developed herein is based on the unit hydrograph theory in hydrology adapted for applications in epidemiology. Inputs are the observed number of daily infected cases $I_o(t)$ and the observed mortality ratio $\mu$, and outputs are the calculated daily recovered cases $R_c(t)$ and the calculated daily death cases $D_c(t)$. The performances of the I-R-D model are obtained by comparing the observed daily recovered cases $R_o(t)$ to the calculated $R_c(t)$, and the observed daily death cases $D_o(t)$ to the calculated $D_c(t)$.

The distribution $I_o(t)$ is separated into two series: $I_R(t)$ representing the part of $I_o(t)$ that would recover, and $I_D(t)$ the part of $I_o(t)$ that would die (Fig. 1). We suppose a simple proportionality

relationship between $I_o(t)$ and each of $I_R(t)$ and $I_D(t)$ according to:

$$I_o(t) = I_R(t) + I_D(t) \quad (1)$$

$$I_R(t) = (1-\mu)I_o(t) \quad \text{and} \quad I_D(t) = \mu I_o(t) \quad (2)$$

The I-R-D model is divided into two models: I-R and I-D. On the basis of the unit hydrograph theory [17, 29] (see details in the Supplementary material) we have

$$\begin{array}{l} \text{I-R model: } R_c(t) = I_R(t) * u_R(t) \quad \text{and} \\ \text{I-D model: } D_c(t) = I_D(t) * u_D(t) \end{array} \quad (3)$$

where the symbol '*' represents the convolution relation, and $u_R(t)$ and $u_D(t)$ [$T^{-1}$] are mathematical kernel functions with an integral equal to 1. For a Dirac (unit) input at $t=0$, the output has the same mathematical equation as $u_R(t)$ (or $u_D(t)$). In hydrology, the functions $u_R(t)$ and $u_D(t)$ are called unit hydrographs [29]. Section ‘Adapting the Hayami kernel function to epidemiology’ shows how these functions can be calculated.

The model also enables us to compare the observed cumulated number of recoveries $R_{to}(t)$ to the calculated $R_{tc}(t)$, the observed cumulated number of deaths $D_{to}(t)$ to the calculated $D_{tc}(t)$, and the observed daily active cases $A_o(t)$ to the calculated $A_c(t)$ with the below equation:

$$\begin{array}{l} A_o(t) = A_o(t-1) + I_o(t) - R_o(t) - D_o(t) \quad \text{and} \\ A_c(t) = A_c(t-1) + I_c(t) - R_c(t) - D_c(t) \end{array} \quad (4)$$

The performances of the I-R model, the I-D model and the active cases are calculated using the Nash–Sutcliffe efficiency (NSE) criteria [35] traditionally used in hydrology, with the below equations, respectively:

$$\begin{align} \text{NSE}_R &= 1 - \frac{\sum_{t=1}^{n} [R_o(t) - R_c(t)]^2}{\sum_{t=1}^{n} [R_o(t) - R_m]^2}; \notag \\ \text{NSE}_D &= 1 - \frac{\sum_{t=1}^{n} [D_o(t) - D_c(t)]^2}{\sum_{t=1}^{n} [D_o(t) - D_m]^2}; \tag{5} \\ \text{NSE}_A &= 1 - \frac{\sum_{t=1}^{n} [A_o(t) - A_c(t)]^2}{\sum_{t=1}^{n} [A_o(t) - A_m]^2} \notag \end{align}$$

with $R_m$ the mean value of $R_o(t)$, $D_m$ the mean value of $D_o(t)$, $A_m$ the mean value of $A_o(t)$ and $n$ the number of time steps. The criteria NSE is inferior to 1 with an optimum corresponding to NSE = 1. NSE = 0 corresponds to using the mean value of observations as a benchmark predictor and is regularly used as a benchmark to compare models [36, 37]. NSE < 0 indicates that the model is a worse predictor than the mean of observations. However the NSE is an indicator of performance which is sensitive to peak values. Increasingly an alternative metric, the Kling–Gupta efficiency (KGE) [38] is used instead with the below equations, respectively

$$KGE_R = 1 - \sqrt{(r_R - 1)^2 + \left(\frac{\sigma_{R,c}}{\sigma_{R,o}} - 1\right)^2 + \left(\frac{\mu_{R,c}}{\mu_{R,o}} - 1\right)^2}; \quad (6a)$$
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** I-R-D model structure. The input $I_o(t)$ is the observed number of daily infected cases divided into $I_p(t)$ and $I_D(t)$ proportional to the mortality ratio $\mu$. $t$ is the time expressed in days. The I-R model calculates the daily number of recoveries $R_c(t)$, with a performance $NSE_R$. The I-D model calculates the daily number of deaths $D_c(t)$, with a performance $NSE_D$. The active cases $A_c(t)$ are calculated with a performance $NSE_A$.</figcaption></figure>

where $r_R$ is the linear correlation between $R_o(t)$ and $R_c(t)$, $\sigma_{R,o}$ is the standard deviation of $R_o(t)$, $\sigma_{R,c}$ is the standard deviation of $R_c(t)$, $\mu_{R,o}$ is the mean of $R_o(t)$ and $\mu_{R,c}$ is the mean of $R_c(t)$:

$$KGE_D = 1 - \sqrt{(r_D - 1)^2 + \left(\frac{\sigma_{D,c}}{\sigma_{D,o}} - 1\right)^2 + \left(\frac{\mu_{D,c}}{\mu_{D,o}} - 1\right)^2}; \quad (6b)$$

where $r_D$ is the linear correlation between $D_o(t)$ and $D_c(t)$, $\sigma_{D,o}$ is the standard deviation of $D_o(t)$, $\sigma_{D,c}$ is the standard deviation of $D_c(t)$, $\mu_{D,o}$ is the mean of $D_o(t)$ and $\mu_{D,c}$ is the mean of $D_c(t)$:

$$KGE_A = 1 - \sqrt{(r_A - 1)^2 + \left(\frac{\sigma_{A,c}}{\sigma_{A,o}} - 1\right)^2 + \left(\frac{\mu_{A,c}}{\mu_{A,o}} - 1\right)^2}; \quad (6c)$$

where $r_A$ is the linear correlation between $A_o(t)$ and $A_c(t)$, $\sigma_{A,o}$ is the standard deviation of $A_o(t)$, $\sigma_{A,c}$ is the standard deviation of $A_c(t)$, $\mu_{A,o}$ is the mean of $A_o(t)$ and $\mu_{A,c}$ is the mean of $A_c(t)$.

Knoben et al. [39] analysed the relationships between the NSE and the KGE, and concluded that NSE and KGE values cannot be directly compared because their relationship is non-unique. They also showed that using the mean observations as a benchmark does not result in KGE = 0 but KGE = -0.41. Thus KGE > -0.41 indicate that a model improves upon the mean observation benchmark.

In the following, we note the input $i(t)$ for $I_R(t)$ or $I_D(t)$ and the output $o(t)$ for $R_c(t)$ or $D_c(t)$. Although $i(t)$ and $o(t)$ are continuous functions in time in hydrological models, herein $i(t)$ and $o(t)$ are discrete functions representing the distribution of the number of cases per day, with $t$ an integer representing the day number. In hydrology, different approaches can be used to calculate $u(t)$ [29].

The exponential model, similar to the solution of the first-order kinetics used in compartmental models in epidemiology, is a unit hydrograph considered as a reference model as below:

$$u(t) = \frac{1}{k} e^{-t/k} \quad (7)$$

where $k$ [T] is a parameter representing linear reservoir retention.

## Adapting the Hayami kernel function to epidemiology

Saint-Venant [40] formulated the system of partial differential equations (continuity and momentum) to describe one-

dimensional, gradually-varied, unsteady flow in rivers. The solution of the Saint-Venant equations has given rise to a number of numerical methods because no analytical solution is available. In practical applications of flood routing in natural channels, the acceleration terms in the Saint-Venant equations can be neglected, and the system is reduced to one parabolic equation, the diffusive wave equation [31, 32]:

$$\frac{\partial Q}{\partial t} + C(Q) \frac{\partial Q}{\partial x} - D(Q) \frac{\partial^2 Q}{\partial x^2} = 0 \quad (8)$$

where $x$ [L] is the length along the channel, $t$ [T] is the time and the celerity $C(Q)$ [$LT^{-1}$] and the diffusivity $D(Q)$ [$L^2T^{-1}$] are functions of the discharge $Q$ [$L^3T^{-1}$]. Let $I(t)$ and $O(t)$ be, respectively, the upstream inflow and the downstream outflow. In the particular case of a semi-infinite channel, no physical downstream boundary condition exists, and $C(Q)$ and $D(Q)$ constant, the diffusive wave equation has as solution the analytical Hayami (1951) equation [30]:

$$O(t) = \int_{0}^{t} I(t - \tau)u(\tau)d\tau = I(t) * u(t) \quad (9)$$

where $u(t)$ is the Hayami kernel function:

$$u(t) = \left(\frac{\theta z}{\pi}\right)^{1/2} \frac{e^{[z(2-\theta/t-t/\theta)]}}{t^{3/2}} \quad (10)$$

with two parameters $\theta$ [T] the lag time and $z$ [dimensionless] a shape parameter. When $z$ tends to zero, the Hayami equation has comparable shape as the exponential model equation (7). For hydraulic applications, equation (10) is generally written as a function of two physically based parameters, changing the two-parameter ($\theta$, $z$) into ($C$, $D$), with $C = L/\theta$ ($L$ being the length) and $D = L^2/(4\theta z)$. The use of the two-parameter ($C$, $D$) is more adapted for hydraulic application because both parameters have physical significance. The use of the two-parameter ($\theta$, $z$) is more adapted for hydrologic applications on catchments because $\theta$ represents the lag time between $i(t)$ and $o(t)$, and can be thus estimated. However, $z$ is empirical and needs to be calibrated. For applications in epidemiology, the parameter $\theta$ represents the lag time between either $I_R(t)$ and $R_o(t)$, or $I_D(t)$ and $D_o(t)$, and hence can be easily interpreted, whereas the parameter $z$ has no significance for the relationships I-R and I-D. This is why equation (10) needs to be rewritten according to two parameters that can be more easily interpreted in epidemiology. For that, $i(t)$ and $o(t)$ can be considered as distributions characterised respectively as below (Fig. 2):

* $G_I$ and $G_O$: the centres of gravity.

* $T_I$ and $T_O$: the abscissae of $G_I$ and $G_O$ which represent the means of $i(t)$ and $o(t)$.

* $s_I$ and $s_O$: the standard deviations of $i(t)$ and $o(t)$, and hence $s_I^2$ and $s_O^2$ are the variances.

The values of $T_I$ and $T_O$ depend on the time origin ($t = 0$) whereas $s_I$ and $s_O$ do not depend on time origin. Following Moussa [31, 32], we have

$$\theta = T_O - T_I \quad \text{and} \quad z = \frac{(T_O - T_I)^2}{2(s_O^2 - s_I^2)} \quad (11)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of the distributions of the input $i(t)$ and the output $o(t)$. $G_I$ and $G_O$ are, respectively, the centres of gravity of $i(t)$ and $o(t)$. $T_I$ and $T_O$ are the abscissae of $G_I$ and $G_O$ representing the means of $i(t)$ and $o(t)$. $s_I$ and $s_O$ are the standard deviations of $i(t)$ and $o(t)$. $\theta$ is the time delay between the two centres of gravity.</figcaption></figure>

The parameter $\theta$ represents the time delay between the two centres of gravity $G_I$ and $G_O$. Let $\tau$ be

$$ \tau = \sqrt{s_O^2 - s_I^2} \qquad (12) $$

The new parameter $\tau$ has a time dimension, and is proportional to the square root of the difference of the variance of $o(t)$ and $i(t)$. Combining equations (11) and (12), we have the relationships:

$$ z = \frac{\theta^2}{2\tau^2} \quad \text{and} \quad \tau = \frac{\theta}{\sqrt{2z}} $$

Substituting equation (13) in (10) we obtain the expression of the Hayami unit hydrograph function of $(\theta, \tau)$ instead of $(\theta, z)$:

$$ u(t) = \left( \frac{\theta^3}{2\pi\tau^2} \right)^{1/2} \frac{e^{\left[ \left( \frac{\theta^2}{2\tau^2} \right)^2 \left( 2-\theta/t-t/\theta \right) \right]}}{t^{3/2}} \quad (14) $$

The Hayami $u(t)$ given by equation (14) is applicable for the continuous function data. But, the Covid-19 pandemic problem dealing with $I-R-D$ and active variables are observed in discrete time interval of 1 day. This necessitates the conversion of the continuous function $u(t)$ to discrete 1-day time interval. To obtain at the discrete $u(t)$ of the Hayami model at $\Delta t$ time interval (with $\Delta t = 1$ day herein), we use the discrete $u(t)$ of the Hayami model using a similar approach as the method described by [41] (see details in the Supplementary material). Another solution consists of subdividing the daily time step into finer numerical time steps under the hypothesis of uniform distribution of data. The finer the time step, the closer one gets to a continuous function of time (see the sensitivity analysis for $\Delta t = 10$ min, 1 h, 3 h and 1 day in Section ‘Sensitivity analysis of the $I-R-D$ Hayami model to the time step of calculation’). The $I-R-D$ model proposed herein was developed in Matlab® and integrated in the MHYDAS hydrological model (MHYDAS-IRD) [42].

The two parameters $(\theta, \tau)$ are related to the means and variances of the two distributions $i(t)$ and $o(t)$ by equations (10) and (11). The advantage of equation (14) using $(\theta, \tau)$ instead of $(\theta, z)$ is that both parameters $(\theta, \tau)$ can be estimated from data

analysis without any calibration. However, the calibration of $(\theta, \tau)$ can be also undertaken in order to improve the model performances.

Figure 3a shows the shape of the output distribution for different values of $\theta$. A small $\theta$ value gives a sharp output distribution with a small time delay between the input and the output, whereas a large $\theta$ value gives a damped output distribution with a large time delay between the input and the output. Figure 3b shows the output distribution for different values of $\tau$. A small $\tau$ value gives a damped and more asymmetric output distribution with a small time delay between the input and the output, whereas a large $\tau$ value gives a sharp and more symmetric output distribution with a large time delay between the input and the output. Figures 3c and d show the cumulative distribution function (CDF) of Figures 3a and b: they represent the time evolution of the total number of $o(t)$ when $i(t)$ is a Dirac delta function.

In the applications, we compare the reference one-parameter exponential model (equation 7) with $k_R$ the parameter of $u_R(t)$ and $k_D$ the parameter of $u_D(t)$, to the two-parameter Hayami model (equation 14) with $(\theta_R, \tau_R)$ the parameters of $u_R(t)$ and $(\theta_D, \tau_D)$ the parameters of $u_D(t)$. The mortality ratio $\mu$ is supposed to be known from the observed data. For both models, the parameters of $u_R(t)$ are calculated separately from those of $u_D(t)$. The performances are $NSE_R$ and $KGE_R$ for $I-R$, $NSE_D$ and $KGE_D$ for $I-D$ and $NSE_A$ and $KGE_A$ for the active cases (equations (5) and (6)). The parameters are calibrated using an iterative automatic trial-and-error method maximising the NSE criteria. This method can be used for both exponential and Hayami models, for datasets $(I_o(t), R_o(t)$ and $D_o(t))$ that are partial (pandemic is evolving) or complete (pandemic is ending).

## Applications on the Covid-19 pandemic

### Data

Applications were implemented on the actual Covid-19 pandemic which occurred in Wuhan, Hubei Province, China, in December 2019. Data are available from Worldometer [20] Covid-19 data used by Johns Hopkins CSSE [19], governmental institutions and many others [21]. Daily data for $I_o(t), R_o(t)$ and $D_o(t)$ are available from 1 January until 19 May 2020 for 24 countries: Australia
<figure><img src="image_5.png" /><figcaption>Fig. 3. Examples of the Hayami unit hydrograph $u(t)$. (a) For $\tau = 10$ days and different values of $\theta$. (b) For $\theta = 20$ days and different values of $\tau$. (c) CDF of (a) and (d) CDF of (b).</figcaption></figure>

(denoted AU), Austria (AT), Belgium (BE), China (CN), Cuba (CU), Czechia (CZ), Denmark (DK), France (FR), Germany (DE), Iceland (IS), Iran (IR), Italy (IT), Japan (JP), Malaysia (MY), New Zealand (NZ), Romania (RO), Slovakia (SK), South Korea (KR), Spain (ES), Switzerland (CH), Thailand (TH), Turkey (TU), USA (US) and the whole world (World). The latter case enables us to study an overall average trend in the whole world.

Figure 4 shows the characteristics of the data for all 24 countries, where for each country $I_t$ is the total number of infected cases, $R_t$ the total number of recoveries and $D_t$ the total number of deaths. Figure 4a shows the large range of variability of $I_t$ with $2 \times 10^3$ (IS, NZ, SK) $< I_t < 5 \times 10^6$ (World). We define also the index $(R_t + D_t)/I_t$ representing the evolution of the pandemic at the date of data availability 19 May 2020: the index is close to 1 when the pandemic is ending (e.g. CN, IS, NZ) and lowest than 0.5 (e.g. BE, US, World) when the pandemic is still actively evolving. Figure 4b shows the large range of variability of the mortality ratio $\mu = D_t/I_t$ with 0.005 (IS) $< \mu < 0.16$ (BE). Probably different data collection strategies have been carried out among the different countries. However, the I-R-D model was applied separately to each country, and therefore the calibrated sets of parameters can be considered as descriptors of the I-R and I-D relationships for each country.

We compare the exponential and the Hayami I-R-D models to simulate the relationships I-R, I-D and the active cases. As we observe noisy signals for $I_o(t)$, $R_o(t)$ and $D_o(t)$, simulations can be improved by smoothing these signals. For both the exponential and the Hayami models, we also compare four smoothing strategies: without smoothing, 3-days moving average, 5-days moving average, and 7-days moving average. A sensitivity analysis was also undertaken to study the impact of the time step of calculation on the calibrated parameters and the model performances, by

subdividing the daily time step into finer numerical time steps ($\Delta t = 10$ min, 1 h, 3 h and 1 day). For the I-R model, we note $NSE_R$, $NSE_{R3}$, $NSE_{R5}$ and $NSE_{R7}$, the values of NSE respectively for the four smoothing strategies. For the I-D model, we note $NSE_D$, $NSE_{D3}$, $NSE_{D5}$ and $NSE_{D7}$. For the active cases, we note $NSE_A$, $NSE_{A3}$, $NSE_{A5}$ and $NSE_{A7}$. For the I-R model, we note $KGE_R$, $KGE_{R3}$, $KGE_{R5}$ and $KGE_{R7}$, the values of KGE respectively for the four smoothing strategies. For the I-D model, we note $KGE_D$, $KGE_{D3}$, $KGE_{D5}$ and $KGE_{D7}$. For the active cases, we note $KGE_A$, $KGE_{A3}$, $KGE_{A5}$ and $KGE_{A7}$.

First, we present an example of application and the sensitivity analysis on China where the pandemic is ending. Then, we show the results for all 24 countries, for both the exponential and Hayami models, and for the three smoothing strategies. The Supplementary material shows the comparison between the observed $\{R_o(t), D_o(t), A_o(t)\}$ and the calculated $\{R_c(t), D_c(t), A_c(t)\}$ for all countries.

## Application on China

First, we study in detail the case of China where the pandemic was ending by the mid of April 2020. The total number of infected $I$ is 82 123 cases, recoveries $R$ is 78 042 cases and deaths $D$ is 3324 cases and the mortality ratio is $\mu = 0.041$. Figures 5 and 6 show the observed $\{R_o(t), D_o(t), A_o(t)\}$ and the calculated $\{R_c(t), D_c(t), A_c(t)\}$ for respectively the exponential and the Hayami models after calibration of the parameters using data with 3-days moving average. For the I-R model, the calibrated parameters of $u_R(t)$ are for the exponential model $k_R = 26.5$ days with fair performance $NSE_{R3} = 0.648$ (Fig. 5a), and for the Hayami model $\theta_R = 21.8$ days and $\tau_R = 10.3$ days (very close to those calculated from equation (10) $\theta_R = 20.3$ days and $\tau_R = 9.5$) with excellent performance
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Data characteristics for the 24 countries, where for each country $I_t$ is the total number of infected cases, $R_t$ the total number of recoveries and $D_t$ the total number of deaths. Data are available from 1 January to 19 May 2020. (a) Pandemic evolution index $(R_t + D_t)/I_t$ function of $I_t$ and (b) mortality ratio: $\mu = D_t/I_t$.</figcaption></figure>

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

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

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

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** I-R-D exponential model application on Covid-19 in China (with 3-days smoothing average of data) using the calibrated parameters. (a) Comparison of the observed and calculated recoveries; (b) comparison of the observed and calculated deaths; (c) comparison of the observed and calculated active cases; (d) comparison of the cumulated observed infected cases, the observed and calculated recoveries, and the observed and calculated deaths. Data are available from 1 January 2020 to 19 May 2020.</figcaption></figure>

$NSE_{R3} = 0.99$ (Fig. 6a). The calibrated values of the lag time $k_r$ and $\theta_R$ are comparable to cure rate for infections (17–20 days) obtained in [4]. For the I-D model, the calibrated parameters of $u_D(t)$ are for the exponential model $k_D = 9.3$ days with very good performance $NSE_{D3} = 0.93$ (Fig. 5b), and for the Hayami model $\theta_D = 10.1$ days and $\tau_D = 11.1$ days with excellent performance $NSE_{D3} = 0.95$ (Fig. 6b). When comparing the observed and calculated active

cases $A(t)$, we observe that the exponential model gives good results $NSE_{A3} = 0.82$ (Fig. 5c) whereas the Hayami model gives much better results $NSE_{A3} = 0.95$ (Fig. 6c). Finally, when comparing the cumulated observed and calculated number of recoveries, and the cumulated observed and calculated number of deaths, we observe also fair performance of the exponential model (Fig. 5d) and excellent performance of the Hayami model (Fig. 6d).
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** I-R-D Hayami model application on Covid-19 in China (with 3-days smoothing average of data) using the calibrated parameters. (a) Comparison of the observed and calculated recoveries; (b) comparison of the observed and calculated deaths; (c) comparison of the observed and calculated active cases; (d) comparison of the cumulated observed infected cases, the observed and calculated recoveries, and the observed and calculated deaths. Data are available from 1 January 2020 to 19 May 2020.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7.** Comparison of the unit hydrographs calibrated for the exponential and Hayami models for the Covid-19 in China: (a) $u_R(t)$; (b) $u_D(t)$; (c) CDF of (a) and (d) CDF of (b).</figcaption></figure>
<table><caption>**Table 1.** I–R–D Hayami model application on Covid-19 in China for different time steps analysis by subdividing the daily time step into finer time steps (Δt = 10 min, 1 h, 3 h and 1 day) and four smoothing strategies (without smoothing, 3-days moving average, 5-days moving average and 7-days moving average): the calibrated parameters of the I–R model (θ<sub>R</sub> and τ<sub>R</sub>) and the corresponding criteria functions (NSE<sub>R</sub> and KGE<sub>R</sub>), the calibrated parameters of the I–D model (θ<sub>D</sub> and τ<sub>D</sub>) and the corresponding criteria functions (NSE<sub>D</sub> and KGE<sub>D</sub>), and the criteria functions corresponding to the actives cases (NSE<sub>A</sub> and KGE<sub>A</sub>).</caption><thead><tr><th rowspan="2">Δt</th><th colspan="4">I-R model</th><th colspan="4">I-D model</th><th colspan="2">Active cases</th></tr><tr><th>θ<sub>R</sub><br>(days)</th><th>τ<sub>R</sub><br>(days)</th><th>NSE<sub>R</sub></th><th>KGE<sub>R</sub></th><th>θ<sub>D</sub><br>(days)</th><th>τ<sub>D</sub><br>(days)</th><th>NSE<sub>D</sub></th><th>KGE<sub>D</sub></th><th>NSE<sub>A</sub></th><th>KGE<sub>A</sub></th></tr></thead><tbody><tr><td colspan="11" style="text-align:center; font-weight:bold;">Without smoothing</td></tr><tr><td>10 min</td><td>21.6</td><td>10.5</td><td>0.973</td><td>0.606</td><td>10.2</td><td>11.0</td><td>0.919</td><td>0.967</td><td>0.997</td><td>0.739</td></tr><tr><td>1 h</td><td>21.6</td><td>10.5</td><td>0.971</td><td>0.604</td><td>10.2</td><td>11.0</td><td>0.918</td><td>0.967</td><td>0.997</td><td>0.739</td></tr><tr><td>3 h</td><td>21.7</td><td>10.4</td><td>0.969</td><td>0.604</td><td>10.1</td><td>11.0</td><td>0.918</td><td>0.967</td><td>0.997</td><td>0.739</td></tr><tr><td>1 day</td><td>21.8</td><td>10.3</td><td>0.966</td><td>0.603</td><td>10.1</td><td>11.1</td><td>0.915</td><td>0.966</td><td>0.997</td><td>0.739</td></tr><tr><td colspan="11" style="text-align:center; font-weight:bold;">Smoothing using 3-days moving average</td></tr><tr><td>10 min</td><td>21.5</td><td>10.5</td><td>0.986</td><td>0.614</td><td>10.3</td><td>11.0</td><td>0.956</td><td>0.982</td><td>0.997</td><td>0.739</td></tr><tr><td>1 h</td><td>21.6</td><td>10.5</td><td>0.986</td><td>0.614</td><td>10.2</td><td>11.0</td><td>0.955</td><td>0.982</td><td>0.997</td><td>0.739</td></tr><tr><td>3 h</td><td>21.6</td><td>10.4</td><td>0.985</td><td>0.612</td><td>10.2</td><td>11.1</td><td>0.953</td><td>0.982</td><td>0.997</td><td>0.739</td></tr><tr><td>1 day</td><td>21.7</td><td>10.3</td><td>0.984</td><td>0.611</td><td>10.2</td><td>11.1</td><td>0.953</td><td>0.982</td><td>0.997</td><td>0.739</td></tr><tr><td colspan="11" style="text-align:center; font-weight:bold;">Smoothing using 5-days moving average</td></tr><tr><td>10 min</td><td>21.5</td><td>10.5</td><td>0.990</td><td>0.618</td><td>10.3</td><td>11.1</td><td>0.971</td><td>0.988</td><td>0.998</td><td>0.740</td></tr><tr><td>1 h</td><td>21.5</td><td>10.5</td><td>0.990</td><td>0.617</td><td>10.3</td><td>11.1</td><td>0.970</td><td>0.987</td><td>0.998</td><td>0.740</td></tr><tr><td>3 h</td><td>21.6</td><td>10.4</td><td>0.989</td><td>0.615</td><td>10.2</td><td>11.2</td><td>0.969</td><td>0.986</td><td>0.998</td><td>0.740</td></tr><tr><td>1 day</td><td>21.6</td><td>10.4</td><td>0.988</td><td>0.615</td><td>10.2</td><td>11.2</td><td>0.968</td><td>0.985</td><td>0.998</td><td>0.740</td></tr><tr><td colspan="11" style="text-align:center; font-weight:bold;">Smoothing using 7-days moving average</td></tr><tr><td>10 min</td><td>21.5</td><td>10.5</td><td>0.993</td><td>0.621</td><td>10.2</td><td>11.1</td><td>0.980</td><td>0.988</td><td>0.998</td><td>0.740</td></tr><tr><td>1 h</td><td>21.6</td><td>10.5</td><td>0.992</td><td>0.619</td><td>10.2</td><td>11.2</td><td>0.980</td><td>0.988</td><td>0.998</td><td>0.740</td></tr><tr><td>3 h</td><td>21.6</td><td>10.5</td><td>0.992</td><td>0.619</td><td>10.3</td><td>11.2</td><td>0.979</td><td>0.987</td><td>0.998</td><td>0.740</td></tr><tr><td>1 day</td><td>21.6</td><td>10.4</td><td>0.992</td><td>0.619</td><td>10.3</td><td>11.2</td><td>0.978</td><td>0.987</td><td>0.998</td><td>0.740</td></tr></tbody></table>

These results show that the one-parameter exponential model is sufficient and give comparable results to the Hayami model only for modelling the relationship I–D. However, for modelling the relationship I–R and the active cases, the Hayami model highly improves the performances in comparison with the exponential model, presenting also the advantage that its parameters can be easily estimated from the analysis of the data distributions of I(t) and R(t). Figure 7 shows the $u_R(t)$ and $u_D(t)$ obtained with the calibrated parameters for China. For $u_R(t)$ (Fig. 7a and the corresponding CDF in Fig. 7c), we observe different behaviours for the two models. The exponential model cannot reproduce a unit hydrograph that rises to a maximum and then falls down, which explains the fair performance of the exponential model. However, the Hayami model succeeds because it enables us to take into account a lag time translation (comparable to the advection process) and an attenuation of the peak (comparable to the dispersion process). For $u_D(t)$ (Fig. 7b and the corresponding CDF in Fig. 7d), we observe very comparable results for both models.

## Sensitivity analysis of the I-R-D Hayami model to the time step of calculation

This section presents a sensitivity analysis of the Hayami I–R–D model on China. Compared to applications in hydrology where the discharge is a continuous function of time, the major

difference for epidemiological applications of the Hayami I–R–D model is that in epidemiology the data are a discrete time function representing the daily number of cases. However, all unit hydrograph approaches, including the Hayami model, can be easily adapted for applications on discrete or continuous data as discussed in Section ‘Adapting the Hayami kernel function to epidemiology’ [41, 43]. Around the reference data time step Δt = 1 day, we conduct a sensitivity analysis by subdividing the daily time step into finer time steps (Δt = 10 min, 1 h and 3 h) under the assumption of a uniform distribution of the number of cases. The finer the time step Δt, the closer one gets to a continuous function of time. For the data used, we also compare four smoothing strategies: without smoothing, 3-days moving average, 5-days moving average and 7-days moving average.

Table 1 shows the values of the calibrated parameters of the I–R model (θ<sub>R</sub> and τ<sub>R</sub>) and the corresponding criteria functions (NSE<sub>R</sub> and KGE<sub>R</sub>), the calibrated parameters of the I–D model (θ<sub>D</sub> and τ<sub>D</sub>) and the corresponding criteria functions (NSE<sub>D</sub> and KGE<sub>D</sub>), and the criteria functions corresponding to the actives cases (NSE<sub>A</sub> and KGE<sub>A</sub>). For the fine time steps (10 min, 1 h and 3 h), the calibrated parameters and the corresponding performance criteria remain very close in comparison with those obtained with the reference time step of 1 day. This is due to the fact that the use of a fine time step only allows tending to continuous time functions (i.e. Δt = 10 min) but using the same
<figure><img src="image_13.png" /><figcaption>**Fig. 8.** For the I–R model (a and d), the I–D model (b and e) and the active cases (c and f), comparison of the exponential model (denoted Exp) and the Hayami model (Hay) for four different smoothing strategies: without smoothing, 3-days moving average (3d), 5-days moving average (5d), and 7-days moving average (7d). The values of the NSE (a, b and c) and KGE (d, e and f) performance criteria are classified by ascending order for the 24 countries.</figcaption></figure>

data as for $\Delta t = 1$ day. However, we observe an improvement in model performances and a slight change in the values of the set parameters when the original data are smoothed ($21.5 < \theta_R < 21.8$ days and $10.3 < \tau_R < 10.5$ days). The improvement of model performances comparing the data without smoothing to 7-days moving average smoothed data are as follows: $NSE_R = 0.966$ and $0.992$, $KGE_R = 0.603$ and $0.619$, $NSE_D = 0.915$ and $0.978$, $KGE_R = 0.966$ and $0.987$, $NSE_A = 0.997$ and $0.998$, $KGE_A = 0.739$ and $0.740$. As the model results are more sensitive to smoothing data then to the use of finer time steps, in the following we limit the applications to compare the performance of the models using the daily time step for the four smoothing strategies (without smoothing, 3-days, 5-days moving and 7-days moving average).

## Application on 24 countries

This section aims to show how the methodology can be applied, and how the parameters could be used to analyse the variability among different countries. For a part of the 24 countries, the pandemic is still evolving, and the series $I_o(t)$, $R_o(t)$ and $D_o(t)$ are still incomplete. Hence, the parameters of both the exponential and the Hayami models calibrated on the available data, may slightly change when the pandemic ends. For each country, we obtain generally similar values of the calibrated parameters for the four smoothing strategies (without smoothing, 3-days moving average, 5-days moving average and 7-days moving average) but the performances (NSE and KGE) may vary drastically with the smoothing strategy.

Figure 8 shows a comparison of the NSE (Figs 8a–c) and KGE (Figs 8d–f) of the exponential and the Hayami models for the four smoothing strategies for all studied countries. For the I–R

model, the 7-days moving average improves the performances NSE (Fig. 8a) of both the exponential and the Hayami models in comparison with the two remaining smoothing strategies. Moreover, the Hayami model performs much better than the exponential model: for the 5-days moving average, $NSE_{R5} > 0.90$ for four countries with the exponential model, and for 20 countries with the Hayami model (Fig. 9a). For the I–D model, Figure 8b shows that for both the exponential and the Hayami models, the 3-days, 5-days and 7-days moving average give similar results but better than without smoothing. The Hayami model gives comparable and slightly better results than the exponential model: for the 5-days moving average, $NSE_{D5} > 0.90$ for 12 countries with the exponential model, and for 16 countries with the Hayami model (Fig. 9b). For the active cases, Figure 8c shows that the three smoothing strategies give similar results for both the exponential and the Hayami models, because the active cases is a cumulated function, and consequently is less sensitive to noisy data. However, the Hayami model gives better results than the exponential model: for the 5-days moving average, $NSE_{A5} > 0.95$ for nine countries with the exponential model, and for all 24 countries with the Hayami model (Fig. 9c).

Similar results are obtained with the KGE. The Hayami model performs much better than the exponential model: for the 5-days moving average, $KGE_{R5} > 0.90$ for two countries with the exponential model, and for 16 countries with the Hayami model (Fig. 9d); $KGE_{D5} > 0.90$ for six countries with the exponential model, and for 12 countries with the Hayami model (Fig. 9e); $KGE_{A5} > 0.90$ for seven countries with the exponential model, and for 14 countries with the Hayami model (Fig. 9e).

Finally, Figure 10 shows the Hayami calibrated parameters for the 24 countries for both the I–R model ($\theta_R$ and $\tau_R$) and the I–D model ($\theta_D$ and $\tau_D$). The parameters range are: 11 days
<figure><img src="image_14.png" /><figcaption>**Fig. 9.** Comparison of the performances NSE (a, b and c) and KGE (d, e and f) of the exponential (denoted Exp) and the Hayami models (smoothing data with 5-days moving average) for: the I–R model (a and d), the I–D model (b and e), and the active cases (c and f).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 10.** Hayami calibrated parameters for the 24 countries for: (a) I–R model ($\theta_R$ and $\tau_R$) and (b) I–D model ($\theta_D$ and $\tau_D$).</figcaption></figure>

(IR) $< \theta_R < 89$ days (US), 1.5 days (IR) $< \tau_R < 99$ days (US), 4 days (DK, MY) $< \theta_D < 30$ days (KR), 1 day (NZ) $< \tau_D < 27$ days (KR). High values of $\theta_R$ (or $\theta_D$) correspond to a long time interval between infection and recoveries (or deaths). High values of $\tau_R$ (or $\tau_D$) correspond to a high difference of the variance between $I_o(t)$ and $R_o(t)$ (or $I_o(t)$ and $D_o(t)$) as given in equation (9) and shown in Figure 2. Some particular cases can be examined such as Iran which has the shortest recovery period of $\theta_R = 11$ days with small $\tau_R = 1$ day in comparison with results of other

countries. A small value of $\theta_R$ can be due either to a late date of detection of infected cases and consequently a shorter recovery period, or to the recovery criterion with, for example a single negative test instead of two for recovered patients, which shortens the contamination period. The small value of $\tau_R$ is an indicator that both $I(t)$ and $R(t)$ have the same dynamic with low diffusion. Conversely, France and the US data are characterised by a large $\tau_R$ (respectively, 82 and 99). This is mainly due to the fact that the $R(t)$ curves were rising and incomplete by mid-May 2020.
<figure><img src="image_16.png" /><figcaption>**Fig. 11.** Comparison between the observed $R_0(t)$ and the calculated $R_c(t)$ recoveries, and between the observed $A_0(t)$ and the calculated $A_c(t)$ active cases, for 24 countries using the Hayami I-R-D model (smoothing data with 5-days moving average).</figcaption></figure>
Consequently, the calibration of the I–R–D model on the NSE criterion favouring the high values of the observation can induce a poor estimation of the parameters. This is not the case where the pandemic is over (e.g. Austria, China, Germany, Iceland, New Zealand, etc.), where small values of $τ_R$ (between 6 and 10 days) are calibrated.

The four parameters ($θ_R$, $τ_R$, $θ_D$ and $τ_D$) are descriptors of the relationships between the distributions $I(t)$, $R(t)$ and $D(t)$. Consequently, they depend on various factors such as the measurement strategy of each country, the health policy, the population density, the presence and date of closure of the local airport, etc. Given the heterogeneities of countries, and the time evolution of health policies, it could be hard to explain as well. However, these parameters remain useful indicators on data acquisition, recoveries strategies and pandemic evolution. They are also useful to compare and classify countries and regions.

Finally, Figure 11 shows the overall very good performances of the I–R–D Hayami model (with 5-days moving average data smoothing) for the 24 countries, comparing the observed $R_0(t)$ and the calculated $R_c(t)$ recoveries, and the observed $A_0(t)$ and the calculated $A_c(t)$ active cases. Calibrated parameters and model performances will probably change when additional data will be available. However, these applications should be seen as first tests that show that the Hayami solution of the diffusive wave equation, and more generally the unit hydrograph theory, can be easily adapted to compartmental I–R–D models in epidemiology. In hydrology, the input $i(t)$ and the output $o(t)$ signals are continuous functions, whereas $I_0(t)$, $R_0(t)$ and $D_0(t)$ functions in epidemiology are discrete functions. Despite discontinuities in data, and the different methods among countries to define the observed $I_0(t)$, $R_0(t)$ and $D_0(t)$, the I–R–D model provides simulations mostly qualified as excellent when slightly smoothing the noisy data. The Hayami model is easy-to-use and parsimonious with only two parameters for each of $u_R(t)$ and $u_D(t)$. One main advantage of the Hayami model is that both parameters ($θ, t$) can be estimated using equation (11) from the analysis of $I_0(t)$, $R_0(t)$ and $D_0(t)$ when the pandemic ends, avoiding the calibration procedure. The parameters, either calibrated or estimated, are useful simple describers to compare the temporal evolution of recoveries and deaths in different countries.

## Conclusion

On the basis of the analogy between SIRD and compartmental models in hydrology, this study makes mathematical formulations developed in hydrology available for modelling in epidemiology. We adapt the ‘transfer functions’ generally used in hydrological modelling to compartmental I–R–D models in epidemiology in order to simulate the relationships between the number of infectious $I(t)$, the number of recovered $R(t)$, the number of death cases $D(t)$ and the number of active cases $A(t)$. Simplified approaches of the transfer functions such as the unit hydrograph are easy-to-use and parsimonious with a low number of parameters. We compare the one-parameter exponential model usually used in SIRD epidemiologic model to the two-parameter physically based Hayami model solution of the diffusive wave equation. Applications were implemented on the recent Covid-19 pandemic.

The application on 24 countries enables us to compare the performances of the two models. The exponential model gives very good performances for modelling the relationship I–D, but fair performances for modelling I–R and the number of active cases.

For I–R, the Hayami model improves significantly the performances with excellent performances for all variables. The Hayami model presents also the advantage that its parameters can be easily estimated from the analysis of the data distributions of $I(t)$, $R(t)$ and $D(t)$.

The Hayami model is parsimonious with only two parameters which are useful simple describers to compare the temporal evolution of recoveries and deaths in different countries with different contamination rates and strategies for recoveries.

These first results illustrate the interest of adapting mathematical formulations developed in a physical discipline like hydrology for applications in epidemiology. This allows epidemiology to benefit from the numerous advances in hydrology, and provides epidemiological modellers simple and easy-to-use parsimonious tools that have been evaluated in the literature, and could possibly make a modest contribution to the complex modelling exercise in epidemiology.

**Supplementary material.** The supplementary material for this article can be found at https://doi.org/10.1017/S0950268821001011

**Acknowledgements.** The authors thank Claude Bocquillon at University of Montpellier, Jean-Claude Davidian at Montpellier SupAgro, Claude Casellas at University of Montpellier and Jean-Paul Lhomme at IRD for fruitful discussions and for their comments on this study. The authors thank the Editor Eduardo Massad and the reviewers for their constructive comments.

**Financial support.** This is a purely theoretical work.

**Conflict of interest.** None.

**Ethical standards.** Not applicable (this is a purely theoretical work with no human subject involved).

**Data availability statement.** Datasets for this research are available on the Worldometer's website https://www.worldometers.info/coronavirus/.

## References

1. Kermack WO and McKendrick AG (1927) Contributions to the mathematical theory of epidemics. *Proceedings of the Royal Society A* **115**, 700–721.

2. Massad E et al. (2020) Two complementary model-based methods for calculating the risk of international spreading of a novel virus from the outbreak epicentre. The case of COVID-19. *Epidemiology and Infection* **148**, e109, 1–6.

3. Wu JT, Leung K and Leung GM (2020) Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. *Lancet (London, England)* **395**, 689–697.

4. Zhou F et al. (2020) Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study. *Lancet (London, England)* **395**, 1054–1062.

5. Sattenspiel L and Dietz K (1995) A structured epidemic model incorporating geographic mobility among regions. *Mathematical Biosciences* **128**, 1–2.

6. Fromont E, Pontier D and Langlais M (2003) Disease propagation in connected host population with density-dependant dynamics: the case of the feline leukemia virus. *Journal of Theoretical Biology* **223**, 465.

7. Auger P et al. (2008) The Ross–Macdonald model in a patchy environment. *Mathematical Biosciences* **216**, 123–131.

8. Finger F et al. (2014) Cholera in the Lake Kivu region (DRC): integrating remote sensing and spatially explicit epidemiological modeling. *Water Resources Research* **50**, 5624–5637.

9. Greenhalgh S and Day T (2017) Time-varying and state-dependent recovery rates in epidemiological models. *Infectious Disease Modelling* **2**, 419–430.
10. Vaidya NK et al. (2015) Modelling the epidemic spread of an H1N1 influ- enza outbreak in a rural university town. *Epidemiology and Infection* **143**, 1610–1620.

11. Kadowaki H et al. (2018) The risk of rabies spread in Japan: a mathem- atical modelling assessment. *Epidemiology and Infection* **146**, 1245–1252.

12. Aswi A et al. (2019) Bayesian spatial and spatio-temporal approaches to modelling dengue fever: a systematic review. *Epidemiology and Infection* **147**, e33, 1–14.

13. Bertuzzo E et al. (2008) On the space-time evolution of a cholera epi- demic. *Water Resources Research* **44**, W01424.

14. Barbosa H et al. (2018) Human mobility: models and applications. *Physics Reports* **734**, 1–74.

15. Soriano-Paños D et al. (2020) Vector-borne epidemics driven by human mobility. *Physical Review Research* **2**, 013312. doi: 10.1103/ PhysRevResearch.2.013312.

16. Massad E, Coutinho FAB and Wilder-Smith A (2019). Modelling an optimum vaccination strategy against ZIKA virus for outbreak use. *Epidemiology and Infection* **147**, e196, 1–8.

17. Beven KJ (2000) *Rainfall-Runoff Modelling: The Primer*. Chichester, New York: John Wiley & Sons.

18. Moussa R, Chahinian N and Bocquillon C (2007) Distributed hydro- logical modelling of a Mediterranean mountainous catchment – model construction and multi-site validation. *Journal of Hydrology* **337**, 35–51.

19. Johns Hopkins Center for Health Security (2020) Daily updates on the emerging novel Coronavirus COVID-19 Global Cases by the Center for Systems Science and Engineering (CSSE) at the Johns Hopkins Center for Health Security. Available at https://gisanddata.maps.arcgis.com/ apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6.

20. **Worldometers** (2020) Worldometer's Covid-19 data. Available at https:// www.worldometers.info/coronavirus/.

21. Dong E, Du H and Gardner L (2020) An interactive web-based dash- board to track COVID-19 in real time. Correspondence. *The Lancet Infectious Diseases* **20**, 533–534. doi: https://doi.org/10.1016/S1473-3099 (20)30120-1.

22. Singh VP (2002) Is hydrology kinematic? *Hydrological Processes* **16**, 667–716.

23. Chevron B and Moussa R (2016) Determinants of modelling choices for 1-D free-surface flow and morphodynamics in hydrology and hydraulics: a review. *Hydrology and Earth System Sciences* **20**, 3799–3830.

24. Feng Z, Xu D and Zhao H (2007) Epidemiological models with non- exponentially distributed disease stages and applications to disease control. *Bulletin of Mathematical Biology* **69**, 1511–1536.

25. Fowler AC and Hollingsworth TD (2015) Simple approximations for epi- demics with exponential and fixed infectious periods. *Bulletin of Mathematical Biology* **77**, 1539–1555.

26. Krylova O and Earn DJD (2013) Effects of the infectious period distribu- tion on predicted transitions in childhood disease dynamics. *Journal of the*

Royal Society Interface **10**, 1–15. doi: https://doi.org/10.1098/ rsif.2013.0098.

27. Murray JD (2007) *Mathematical Biology*, 3rd Edn., Vol. 1. New York: Springer.

28. Gallos LK and Fefferman NH (2015) The effect of disease-induced mor- tality on structural network properties. *PLoS ONE* **10**, e0136704.

29. Singh PK, Mishra SK and Jain MK (2014) A review of the synthetic unit hydrograph: from the empirical UH to advanced geomorphological meth- ods. *Hydrological Sciences Journal* **59**, 239–261.

30. Hayami S (1951) On the propagation of flood waves. *Bulletins – Disaster Prevention Research Institute, Kyoto University, Kyoto*, pp. 1–16.

31. Moussa R (1996) Analytical Hayami solution for the diffusive wave flood routing problem with lateral inflow. *Hydrological Processes* **10**, 1209–1227.

32. Moussa R and Bocquillon C (1996) Criteria for the choice of flood rout- ing methods in natural channels. *Journal of Hydrology* **186**, 1–30.

33. Moussa R and Bocquillon C (1996) Algorithms for solving the diffusive wave flood routing equation. *Hydrological Processes* **10**, 105–124.

34. Moussa R and Majdalani S (2019) Evaluating lateral flow in an experi- mental channel using the diffusive wave inverse problem. *Advances in Water Resources* **127**, 120–133.

35. Nash JE and Sutcliffe JV (1970) River flow forecasting through concep- tual models I: a discussion of principles. *Journal of Hydrology* **10**, 282–290.

36. Schaefli B and Gupta HV (2007) Do Nash values have value? *Hydrological Processes* **21**, 2075–2080.

37. Moussa R (2010) When monstrosity can be beautiful while normality can be ugly: assessing the performance of event-based flood models. *Hydrological Sciences Journal* **55**, 1074–1084.

38. Gupta HV et al. (2009) Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological mod- eling. *Journal of Hydrology* **377**, 80–91.

39. Knoben WJM, Freer JE and Woods RA (2019) Technical note: inherent benchmark or not? Comparing Nash–Sutcliffe and Kling–Gupta efficiency scores, *Hydrology and Earth System Sciences* **23**, 4323–4331.

40. Saint-Venant B (1871) Théorie du movement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leurs lits. *Comptes Rendus des Séminaires de l'Académie des Sciences* **73**, 147–154, 237–240.

41. Todini E and Bossi A (1986) PAB (parabolic and backwater) an unconditionally stable flood routing scheme particularly suited for real time forecasting and control. *Journal of Hydraulic Research, IAHR* **24**, 405–424.

42. Moussa R, Voltz M and Andrieux P (2002) Effects of the spatial organization of agricultural management on the hydrological behaviour of a farmed catchment during flood events. *Hydrological Processes* **16**, 393–412.

43. O'Connor KM (1976) A discrete linear cascade model for hydrology. *Journal of Hydrology* **29**, 203–242.