Dynamics of trace metals in a shallow coastal ecosystem: insights from the Gulf of Gabès (southern Mediterranean Sea)

1.   Introduction

1.1.   Problem statement

This paper contributes to the macroscopic modelling of crowd movements. We consider the following initial-boundary value problem for a non-local scalar conservation law that describes the evolution of the local density $\rho$ of pedestrians as a function of time $t$ and position $\mathbf{x} = (x_1, x_2)$ on a crowd evolution domain $\Omega\subset \mathbb R^2$:

Here $\boldsymbol\nu = (\nu_1, \nu_2)$ is a vector field that (with slight abuse of notation) is defined as

where $\boldsymbol\mu$ is the (normalized) fixed smooth vector field of preferred directions (e.g., given by the regularized solution of an eikonal equation), and ${\mathcal I}[\rho(t)]$ is a non-local correction term that depends on the current density distribution. This notation indicates a functional dependence, i.e., ${\mathcal I}$ depends on the function $\rho(t) : = \rho(t, \cdot)$ as a whole. The function $V: \mathbb R^+ \to \mathbb R$ represents the average speed of pedestrians as a function of their density.

We assume that pedestrians move in a space surrounded by walls, and that the vector field $\boldsymbol\nu$ points inward along the boundary $\partial\Omega$ of $\Omega$, that is

where $\mathbf{n}$ is the outward normal to $\Omega$. (Of course, we may assume $\boldsymbol\mu(\mathbf{x})\cdot \mathbf{n} (\mathbf{x})\leq 0$ for all $\mathbf{x}\in \partial\Omega$, then it is enough to ensure that also ${\mathcal I}[\rho(t)](\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\leq 0$ for all $\mathbf{x}\in \partial\Omega$.) In this case, the condition $\rho(t, \mathbf{x}) = 0$ on $\partial\Omega$ corresponds to a zero-flux condition. If, for simulation reasons, we need to consider smaller domains and to add adsorbing conditions on the part of the boundary not corresponding to walls (and where the vector field points outwards), suitable modifications of the model are needed.

If $\operatorname*{supp} (\rho_0)\subset\Omega$ and $\boldsymbol\nu(t, \mathbf{x})\cdot \mathbf{n} (\mathbf{x})\leq 0$ for all $\mathbf{x}\in \partial\Omega$ and $t > 0$, then problem (1.1) is equivalent to the Cauchy problem

Following Colombo et al. ^[1], we consider a non-local interaction term of the form

where $\eta$ is a smooth non-negative kernel with compact support such that $\iint_{ \mathbb R^2} \eta(\mathbf{x})\, \mathrm{d} \mathbf{x} = 1$ and $\varepsilon$ is a model parameter modulating the magnitude of the interaction term. The term (1.5) models how pedestrians account for other pedestrian distribution close to them to correct their path. To better account for the reaction of pedestrians to densities ahead of them, one may consider anisotropic kernels $\eta$, see e.g., ^[1] and [2,Appendix D]. To account for the presence of boundaries, and walls (or obstacles) into boundaries, unlike ^[3,4], we modify the usual convolution product as follows:

where $\rho_w: \mathbb R^2\to \mathbb R^+$ is defined as

with $R_{w} \geq R$ big enough so that (1.3) holds. Here we denote by

the diameter of a set $A\subset \mathbb R^2$ and by

the "ball" of radius $\ell$ around $\Omega$. Furthermore, we define

which is finite if $\Omega$ is bounded ($\chi$ denotes the characteristic function). See Figure 1.

The presence of high density values at the wall and obstacle locations included in (1.6) and (1.7) is intended to mimic their effect on the pedestrian dynamics. Indeed, in this way the non-local correction term (1.5) "sees" the presence of the wall and deviates the movement from the desired trajectory, thus acting as a discomfort term expressing the tendency of pedestrians to stay away from obstacles.

Remark 1. An explicit condition ensuring $\nabla (\eta\ast_w\rho(t))(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\geq 0$ for all $\mathbf{x}\in \partial\Omega$ (and thus $\boldsymbol\nu(t, \mathbf{x})\cdot \mathbf{n} (\mathbf{x})\leq 0$) is the following:

where $f(\mathbf{x})^-: = \max\{-f(\mathbf{x}), 0\} = -\min\{f(\mathbf{x}), 0\}$ denotes the negative part of a function $f$ and $\Omega^c: = \mathbb R^2\setminus \Omega$ is the complement of $\Omega$. Indeed we have

1.2.   Related works

Macroscopic models of (vehicular and pedestrian) traffic flow are based on balance laws describing the spatio-temporal evolution of averaged quantities such as density and mean velocity. In analogy to vehicular traffic models, macroscopic crowd motion models were introduced starting from the beginning of this century based on scalar conservation laws ^[5,6], gas dynamics equations ^[7,8], gradient flow methods ^[9,10] and time evolving measures ^[11].

More recently, models consisting in non-local conservation laws in two-space dimension were proposed by different authors ^{[1,3,4,12,13,14,15,16]}. Conservation laws with non-local flux functions arise in a large variety of applications, such as traffic flow ^[17,18,19], sedimentation ^[20,21], and material flows on conveyor belts ^[22]. In particular, anisotropic kernels have been introduced in traffic models to account for downstream interactions, see e.g., ^[17,19]. Physically based derivations of vision-based interaction pedestrian dynamics can be found for example in ^{[11,23,24,25]}. We emphasize that none of the previously cited models (with the obvious exception of ^[1]) fits into the framework of the present approach, which is focused on the appropriate incorporation of boundary conditions and obstacles. That said, we are well aware of a vast body of literature on empirical and numerical studies of pedestrian flows in confined areas that call for the explicit incorporation of such boundaries, see for instance recent proceedings volumes of conferences on Pedestrian and Evacuation Dynamics ^[26,27,28]. (We refer in Section 4 to some specific contributions within these collections.)

The computation of numerical solutions for non-local models is challenging due to the high non-linearity of the system and the dependence of the flux function on the convolution product. In order to overcome the computational bottleneck, high order schemes have been developed for scalar equation ^[29] and systems ^[30,31] in one space dimension. For crowd dynamics models, first order finite volume approximations based on the Lax-Friedrichs scheme have been used in ^[1,3,4], aiming at demonstrating convergence properties for existence results. Here, we will consider the finite difference weighted essentially non-oscillatory (WENO) schemes developed in ^[32] to achieve high-resolution spatial accuracy. WENO schemes ^[33,34,35] are widely employed for the simulation of complex flow fields, their main advantage being their capability to achieve arbitrarily high order formal accuracy in smooth regions while maintaining stable, non-oscillatory and sharp discontinuity transitions.

1.3.   Outline of the paper

The remainder of this paper is organized as follows. In Section 2, we deal with the well-posedness of problem (1.4), (1.5)–(1.7) under the hypothesis (1.3), so that (1.4), (1.5)–(1.7) is equivalent to (1.1), (1.5)–(1.7). Section 3 describes the WENO scheme used to compute approximate solutions and Section 4 collects three different numerical tests investigating the model characteristics. Conclusions and perspectives are elaborated in Section 5.

2.   Well-posedness

We summarize the assumptions concerning the domain $\Omega$ and the functions $V$, $\boldsymbol{\nu}$, $\boldsymbol{\mu}$ and $\eta$:

$\mathbf{(I1)}$ The domain $\Omega\subseteq \mathbb R^2$ is a non-empty bounded open set with smooth boundary $\partial\Omega$, so that the outward normal $\mathbf{n}(\mathbf{x})$ is uniquely defined for all $\mathbf{x}\in \partial\Omega$.

$\mathbf{(I2)}$ The speed function $V\in\mathbf{C^{2}}(\mathbb R; \mathbb R)$ is non-increasing, $V(0) = V_{\max}$ and $V(R) = 0$ for some constants $V_{\max}, R > 0$.

$\mathbf{(I3)}$ The vector field of preferred directions $\boldsymbol\mu\in (\mathbf{C^{2}}\cap\mathbf{W^{{1, \infty}}})(\mathbb R^2; \mathbb R^2)$ is defined such that $\mathinner{ \mathop{ {{\rm{div}}} } } \boldsymbol\mu \in (\mathbf{W^{{1, 1}}}\cap\mathbf{W^{{1, \infty}}})(\mathbb R^2; \mathbb R)$.

$\mathbf{(I4)}$ The vector field $\boldsymbol\nu$ points inward along the boundary $\partial\Omega$ of $\Omega$, i.e., $\boldsymbol\nu(t, \mathbf{x})\cdot \mathbf{n} (\mathbf{x})\leq 0$ for all $\mathbf{x}\in \partial\Omega$, $t\geq 0$.

$\mathbf{(I5)}$ The kernel function $\eta\in\mathbf{C_c^{3}}(\mathbb R^2; \mathbb R^+)$ satisfies $\iint_{ \mathbb R^2} \eta(\mathbf{x})\; \mathrm{d} \mathbf{x} = 1$.

Solutions of problem (1.4), (1.5)–(1.7), equivalently (1.1), (1.5)–(1.7), are intended in the following sense.

Definition 1. [1,Def. 2.1] For any $T > 0$ and $\rho_0\in \mathbf{L^1}(\mathbb R^2, [0, R])$ such that $\operatorname*{supp} \rho_0 \subset\Omega$, a function $\rho\in\mathbf{C^{0}}([0, T], \mathbf{L^1}(\mathbb R^2;[0, R])$ is said to be a weak entropy solution to (1.4), (1.5)–(1.7), if it is a Kružkov entropy solution to the Cauchy problem

i.e., for all $\kappa\in \mathbb R$ and all test functions $\phi\in {\mathbf{C}}_{\mathbf{c}}^\infty (]-\infty, T[\, \times \mathbb R^2; \mathbb R^+)$ there holds

Indeed, by mass conservation, we have that $\operatorname*{supp} \rho(t, \cdot) \subset\Omega$ for all $t > 0$, and therefore $\rho(t, \mathbf{x}) = 0$ for a.e. $\mathbf{x}\in\Omega^c$, see [1,Proposition 3.1]. Therefore, by abuse of notation, $\rho\in{{\mathbf{L}}^\infty } ([0, T]\times\Omega; \mathbb R)$ can be seen as a Kružkov semi-entropy solution in the sense of [37,Def. 1]: namely, from (2.2) we have

where $f(\rho): = \rho V(\rho)$, $s^+ = \max\{s, 0\}$, $s^- = \max\{-s, 0\}$, $\text{sgn}^+(s) = \text{sgn}(s^+)$ and $\text{sgn}^-(s) = -\text{sgn}(s^-)$.

If $\rho\in({{\mathbf{L}}^\infty } \cap \rm{BV})([0, T]\times\Omega; \mathbb R)$, then the classical definition introduced by Bardos, Le Roux and Nédélec ^[38] holds:

where $\text{tr}\rho$ denotes the trace of $\rho$ at the boundary $\partial\Omega$.

We refer the reader to ^[36] for a discussion on the different notions of admissible solutions for scalar multi-dimensional initial-boundary value problems and their equivalence.

Under hypotheses (I1)–(I3), the non-local term ${\mathcal I}$ defined by (1.5)-(1.6)-(1.7) satisfies the following estimates for $\varepsilon$ as in (1.5), $R_w$ as in (1.7) and $M$ as in (1.8):

Moreover, for any $r_1, r_2\in\mathbf{L^1}(\Omega; [0;R])$ there hold

see [1,Proof of Lemma 3.1].

Recalling that $f(\rho) = \rho V(\rho)$ and following [1,Theorem 2.1], we have the following well-posedness result.

Theorem 1. Let (I1)–(I4) hold and $\rho_0 \in (\mathbf{L^1}\cap \rm{BV})(\mathbb R^2;[0, R])$ with $\operatorname*{supp} \rho_0 \subset\Omega$. Then there exists a unique weak entropy solution $\rho\in\mathbf{C^{0}}(\mathbb R^+; \mathbf{L^1}(\mathbb R^2;[0, R]))$ to (1.4), (1.5)–(1.7) with $\operatorname*{supp} \rho(t, \cdot) \subset\Omega$ for $t > 0$. Moreover, $\rho$ satisfies the following estimates

where we define

Stability with respect to $\rho_0$, $v$ and $\boldsymbol\mu$ also holds from [1,Theorem 2.2].

PROOF OF THEOREM 1. Following ^[1], we have to check that we fit the required hypotheses. First of all, given any $r\in\mathbf{C^{0}}([0, T]; \mathbf{L^1}(\Omega; [0, R]))$, we verify that the scalar conservation law

with $\mathbf{w}(t, \mathbf{x}) = \boldsymbol\mu(\mathbf{x})+ {\mathcal I}[r(t)](\mathbf{x})$, satisfies the assumptions of [39,Theorem 5 and Sec. 5.4], and therefore admits a weak entropy solution $\rho\in{{\mathbf{L}}^\infty } (\mathbb R^+; \mathbf{L^1}(\Omega; [0, R]))$ (see [1,Lemma 2.1]). Setting $\varphi(t, \mathbf{x}, \rho) = f(\rho)\mathbf{w}(t, \mathbf{x})$, it is easy to check that $\varphi, \partial_\rho\varphi, \partial^2_{x_i, \rho}\varphi, \partial^2_{x_i, x_j}\varphi\in\mathbf{C^{0}}([0, T]\times \mathbb R^2\times[0, R])$, $\partial_\rho\varphi\in{{\mathbf{L}}^\infty } ([0, T]\times \mathbb R^2\times[0, R])$ and $\text{div} \ \varphi\in{{\mathbf{L}}^\infty } ([0, T]\times \mathbb R^2\times[0, R])$, thanks to (I2), (I3) and (I4). Moreover, by (2.6), we have that

which guarantees that $\rho\in\mathbf{C^{0}}(\mathbb R^+; \mathbf{L^1}(\Omega; [0, R]))$. Moreover, by (2.5),

and by (2.7)

Therefore, by [40,Theorem 2.2], we get (2.11).

Estimates (2.8)–(2.10) ensure that the same fixed point argument used in [1,Proof of Theorem 2.1] can be applied to prove the existence and uniqueness of solutions to (1.4).

3.   High-resolution numerical schemes

3.1.   Notation and semi-discrete formulation

We take $\Omega = [x_1^{\min}, x_1^{\max}]\times[x_2^{\min}, x_2^{\max}]$ and denote by $\rho:[0, \infty[\, \times\Omega\to[0, R]$ and

the solution and the flux function of problem (1.1)–(1.4). We use a uniform Cartesian grid with nodes $(x_1^i, x_2^j)$, $i = 1, \dots, M_1$ and $j = 1, \dots, M_2$ such that $x_1^i = (i-1/2)h$, $x_2^j = (j-1/2)h$, $h = (x_1^{\max}-x_1^{\min})/M_1 = (x_2^{\max}-x_2^{\min})/M_2$. This corresponds to $M_1\times M_2$ grid points $\mathbf{x}_{\mathbf{i}}: = (x_1^i, x_2^j)$, where $\mathbf{i} = (i, j)\in\mathcal{M}: = \{1, \dots, M_1\}\times\{1, \dots, M_2\}\subset {\mathbb N}^2$, and as in ^[9] we utilize two-dimensional unit vectors $\mathbf{e}_1: = (1, 0)$ and $\mathbf{e}_2 : = (0, 1)$ to address neighbouring grid points $\mathbf{x}_{\mathbf{i}+\mathbf{e}_1} = (x_1^{i+1}, x_2^j)$ and $\mathbf{x}_{\mathbf{i}+ \mathbf{e}_2} = (x_1^{i}, x_2^{j+1})$.

We define $\mathbf{u}:[0, \infty)\to \mathbb R^{M_1\times M_2}$ as a solution computed at an instant $t$ in the grid points where

In Section 3.2 we will discuss the discretization of the convolution product. Using this notation, we may approximate the solution of (1.1)–(1.4) in semi-discrete form (that is, discrete in space but continuous in time) by a system of ODEs

where $\mathcal{C}(\mathbf{u})$ represents the spatial discretization of the convective term with entries given by $\mathcal{C}(\mathbf{u}) = (\mathcal{C}(\mathbf{u})_{\mathbf{i}})_{\mathbf{i}\in\mathcal{M}}$ with

where $\smash{\hat{f}_{\mathbf{i}+\frac12\mathbf{e}_1}}$ and $\smash{\hat{f}_{\mathbf{i}+\frac12\mathbf{e}_2}}$ are the numerical fluxes, which in this paper will be a fifth-order version. To this end, we require the summands for $l = 1$ and $l = 2$ in (3.2) to be of the same order of approximation to $\partial f /\partial{x_1}$ and $\partial f /\partial {x_2}$, respectively, at $\mathbf{x} = \mathbf{x}_{\mathbf{i}}$. For upwinding and stability, a flux function $f(\rho)$ is split as follows:

Then each component is approximated separately using its own "wind direction" with respect to $\boldsymbol{e}_l$. The simple Lax-Friedrichs flux splitting

with a suitable viscosity coefficient $\alpha > 0$ is used in this paper. We herein use

in direction $\boldsymbol{e}_k$, $k = 1, 2$. The numerical fluxes are split accordingly, i.e.,

where $\mathcal{R}^{\pm}\bigl({f}_{\mathbf{i}+ (-r:r)\mathbf{e}_k}\bigr) = \mathcal{R}^{\pm}\bigl({f}_{\mathbf{i}-r\mathbf{e}_k}, \dots, {f}_{\mathbf{i}+r\mathbf{e}_k} \bigr)$ denotes $(2r-1)$th-order WENO upwind-biased reconstructions for $r = 2, 3, 4$, see ^[33,34,35].

Ghost cells are needed to compute numerical fluxes near the boundary. To handle these cases we use the boundary condition in (1.1) and set $u_{\mathbf{i}} = 0$ if $\mathbf{i} \notin\mathcal{M}$ and as the vector field $\boldsymbol\nu$ points inward along $\partial\Omega$, i.e., $\boldsymbol\nu\cdot \mathbf{n} (\mathbf{x})\leq 0$ for $\mathbf{x}\in \partial\Omega$, we set

For evacuation problems, to avoid dealing with extended domains, we have to handle a vector field $\boldsymbol\nu$ which points outward at the exit door $D\subset \partial\Omega$, i.e., $\boldsymbol\nu\cdot \mathbf{n} (\mathbf{x}) > 0$ for $\mathbf{x}\in D$. In this case, we set

For more details about the implementation of high-order finite difference WENO schemes for crowd dynamics see ^[9,32].

3.2.   Discretization of the convolution term

In order to evaluate the non-local term in (1.5), we take into account that $\nabla(\eta\ast_w \rho) = \nabla\eta \ast_w \rho$, where the convolution term $\ast_w$ is defined by (1.6). The corresponding convolutions $(\partial\eta/\partial x_1)\ast_w \rho$ and $(\partial\eta/\partial x_2)\ast_w \rho$ are calculated approximately on the discrete grid via a quadrature formula, in our cases a composite Simpson rule. Since $\mbox{supp}(\eta)\subset[-n_0h, n_0h]\times[-n_0h, n_0h]$ for $n_0\in {\mathbb N}$ large enough, any convolution product is given by

where $c_p$ and $c_q$ are the coefficients in the quadrature rule and $\mathbf{p} = (p, q)$. This formula for $\mathbf{u} = (u_{\mathbf{i}})\in \mathbb R^{M_1\times M_2}$ and for the convolution product (1.6) can be written as

where $u_{w, \mathbf{i}}$ is a discrete version of the function (1.7) defined by

Clearly, the discrete convolution (3.6) causes a computational bottleneck. This is a classical problem in scientific computing that is effectively handled by fast convolution algorithms, mainly based on Fast Fourier Transforms ^[41] (see also ^[13,14]).

3.3.   Time discretization

Finally, the semi-discrete scheme (3.1) is discretized by a third-order TVD Runge-Kutta time discretization method ^[34] that can be specified as follows. Assume that $\mathbf{u}^n$ is the vector of approximate solutions at $t = t_n$. Then the approximate values $\mathbf{u}^{n+1}$ associated with $t_{n+1} = t_n + \Delta t$ are calculated by

The combined space and time discretizations define a fully discrete scheme. Observe that the time evolution is only of third order of accuracy, while the spatial discretization is fifth-order accurate. To have fifth-order accuracy also in time, a two-step SSP Runge-Kutta method should be used.

4.   Numerical simulations

We aim at investigating the effects of the non-local operator (1.5)–(1.7) from the crowd dynamics modelling point of view. To this end, in the following numerical examples, we solve numerically (1.4)–(1.7) for $t\in[0, T]$ and $\mathbf{x}\in\Omega$ by using the high-resolution numerical scheme described in Section 3. In particular, we consider here a FD-WENO5 scheme with fifth-order of accuracy in space, see ^[32] for a discussion of high-order finite difference WENO schemes with different orders of accuracy for crowd dynamics. For each iteration, the time step $\Delta t$ in (3.8) is determined by the formula

We choose $C_{\mathrm{cfl}}$ as the largest multiple of $0.05$ that yields oscillation-free numerical solutions. In all numerical tests we have used $C_{\mathrm{cfl}} = 0.2$.

4.1.   Example 1: Comparison with Colombo and Rossi ^[3,4].

In contrast to (1.6)–(1.7), the model proposed in ^[3,4] uses the following definition of the convolution product in the non-local term (1.5)

We remark that the model in ^[3,4] also displays non-locality in the speed, but this prevents a global maximum principle from holding. Therefore, here we keep the speed dependence on the density local.

We consider the evacuation problem proposed in [4,Section 2.1], where $\Omega = R\setminus (C_1\cup C2)$ is composed by a room $R = [0, 8]\times[-2, 2]$ with two symmetric columns $C_1 = \, ]4.5, 7[\, \times\, ]0.8, 1.5[$ and $C_2 = \, ]4.5, 7[\, \times\, ]-1.5, -0.8[$ that guide the evacuation through the exit door set at $D = \{8\}\times\, ]-0.8, 0.8[$. The functions and parameters are chosen as

where $\chi_{[0, l]}(\mathbf{x})$ denotes the characteristic function. The fixed vector field $\boldsymbol\mu(\mathbf{x})$ is given by the unit vector tangent to the geodesic from $\mathbf{x}$ to the exit door, see Figure 2a. Since the non-local term defined by (4.1) does not guarantee that the resulting direction of motion points inside the domain ($\boldsymbol\nu(t, \mathbf{x})\cdot \mathbf{n} (\mathbf{x})\leq 0$ for $\mathbf{x}\in \partial\Omega$, $t\geq 0$), we need to add to $\boldsymbol\mu$ a (fixed) discomfort vector field with maximal intensity along the walls as in ^[3,4], resulting in the vector field showed in Figure 2b.

We display numerical approximations computed with FD-WENO5 scheme with $h = 1/80$ at times $T = 1, 3$ and $6$ for two kernel supports ($l = 0.45$ in Figure 3 and $l = 0.9$ in Figure 4). We observe that the non-local correction term (1.5) allows pedestrians to "see" the presence of the wall and obstacles, and to deviate the movement from the desired trajectory. For $l = 0.45$, we can observe that pedestrians can pass between the obstacles and the wall, as in the Colombo-Rossi model, however this effect is less remarkable for larger kernel supports like $l = 0.9$. Indeed, comparing Figures 3 and 4, we can see that a larger kernel support corresponds to a wider discomfort effect, impacting the velocity vector field on larger portions of the walkable domain. More generally, qualitative differences between the two models depend on the parameter choices of $R_w$ and the discomfort vector field. Nevertheless, we remark that our definition of the convolution (1.6)–(1.7) qualitatively captures the discomfort due to the presence of walls and obstacles delimiting the walking domain.

Figure 5 illustrates the effect of the parameter value $R_w$. Specifically, the scenario of the left column of Figure 5 is the same as in the left column of Figure 3, with $l = 0.45$ and $R_w = 1.5$ replaced by $R_w = 2$. On the right column of Figure 5, we change also $l$ to $l = 0.9$ and compare it to Figure 4, left. Larger values of $R_w$ increase the "discomfort" effect exerted by walls and obstacles, as it becomes apparent if one compares, for instance, the widths of the streams of pedestrians to either side and between the obstacles: they are thinner for $R_w = 2$ than for $R_w = 1.5$.

Finally, Figure 6 shows the impact of the parameters $R_w$ and $l$ on the evacuation time for both models. We can observe that the model is quite sensitive to changes in the convolution support since, as previously mentioned, this affects heavily the resulting velocity field near obstacles and walls. Furthermore, the evacuation process for a given value of $l$ is slower for $R_w = 2$ than for $R_w = 1.5$. We mention that the parameters and obstacles have been chosen to illustrate the model and the difference to ^[3,4], but do not correspond to a specific real-world scenario. There are however, abundant experimental studies providing data on pedestrian flows with obstacles that could form the basis for future comparison with simulations, see e.g., ^{[42,43,44,45]}.

4.2.   Example 2: Comparison between isotropic and anisotropic kernels (lane formation).

In this section, we consider that pedestrians have a limited vision field oriented in a given direction. We study a simple example of evacuation of a rectangular room $\Omega = [0, 8]\times[-3, 3]$, where the vector field $\boldsymbol\mu(\mathbf{x}) = (1, 0)$ is fixed constantly oriented towards the right of the domain. We investigate the influence of a conic convolution kernel on the evacuation dynamics and the pattern formation. We take the kernel function $\eta(\mathbf{x})$ given in (4.2) and cut a conic section $\eta(\mathbf{x})\chi_{\mathcal{S}(\mathbf{x}, l, \alpha, \boldsymbol{\gamma})}(\mathbf{x})$ of angle $2\alpha$ oriented to direction $\boldsymbol{\gamma}(\mathbf{x})$ which is described by the region

(see Figure 7). The section $\eta \chi_{\mathcal{S}(\mathbf{x}, l, \alpha, \boldsymbol{\gamma})}$ is smoothed by convolution with a Gaussian kernel $g(\mathbf{x}) = \exp(-(\| \mathbf{x}\|^2/2\sigma))$ with $\sigma = 5\times 10^{-4}$, then normalized and finally shifted by the distance $0.04$ to ensure that the maximum of the normalized smoothed kernel is centered in $(0, 0)$. Other parameters are taken as

We compare the dynamics given by different kernel orientations and different angle amplitudes. Besides the circular symmetric kernel $\eta$ (i.e., $\alpha = \pi$), we consider

● $\boldsymbol{\gamma}(\mathbf{x}) = (-1, 0)$, corresponding to forward interaction (by central symmetry of the convolution product), and $\alpha = \pi/2, \; \pi/4$;

● $\boldsymbol{\gamma}(\mathbf{x}) = (1, 0)$, corresponding to backward interaction, and $\alpha = \pi/4$;

● $\boldsymbol{\gamma}(\mathbf{x}) = (0, -1)$, modeling pedestrian reacting to the presence of other individuals (or obstacles) on their left, and $\alpha = \pi/4$;

● $\boldsymbol{\gamma}(\mathbf{x}) = (-1, -1)$, for forward-left interactions, and $\alpha = \pi/4$.

The corresponding convolution kernels are depicted in Figure 7. In Figures 8–10, we display numerical approximations and the corresponding vector fields of preferred directions computed with FD-WENO5 scheme with $h = 1/80$ at times $T = 0.1$, $T = 0.2$ and $T = 0.4$. We observe that symmetric and half-circular downstream interaction kernels lead to similar patterns consisting of horizontal lanes (see Figure 8), while narrower angles $\alpha$ lead to vertical patterns, both for forward and backward interactions (displayed in Figure 9). Finally, lateral interactions mostly lead to diagonal waves (Figure 10). Again, the scenario is hypothetical, but could be compared with experimental results from the literature ^{[26,27,28,46]}.

4.3.   Example 3: Room evacuation with and without obstacles

In this section, we consider the problem of evacuating people in a room through a door. In particular, we are interested in studying the impact of the presence of obstacles in front of the door on the evacuation time. This problem has already been discussed by several authors, see for example ^[4,47,48] and references therein. In these works, the authors infer that, in some cases, the location and size of the obstacles may speed up the population to the exit. The interest in this phenomenon, sometimes called "faster-is-slower effect" ^[49,50], goes back to the well-known paper by Helbing et al. ^[51]. Delays during evacuation due to other effects have been investigated, for instance, in ^[52,53]. For this example, we consider the walking domain available to pedestrians is $\Omega = R\setminus C_{i}$, $i = 1, 2, 3$, where the room $R = [0, 8]\times[-3, 3]$ contains obstacles $C_{i}$. The door $D$, the functions $v$ and $\eta$, and the parameter $\varepsilon$ and $R_w$ are the same as in (4.2). The initial condition is a linear combination of characteristic functions with values $0.9$ in $[0.5, 2]\times[-2.2, 0]$, $0.6$ in $[0.5, 2.2]\times[0, 2.2]$, $0.5$ in $[2,4]\times[-2.2, 0]$ and $0.8$ in $[2.2, 4]\times[0, 2.2]$ (see Figure 11). The obstacles and parameter $l$ for the three different scenarios considered are

The numerical solutions for the three scenarios at three different times is displayed in Figure 12, where we have used the FD-WENO5 scheme to computed the approximate solutions.

Figure 13 shows the time evolution of the mass inside the room. We observe that in case $C_1$ the evacuation time is higher in comparison to cases $C_2$ and $C_3$, where the evacuation time is lower, especially for case $C_3$. This is confirmed by Figure 12, which shows the evacuation of the populations at times $T = 6, 20$ and $32$, for the cases $C_1$, $C_2$ and $C_3$ respectively.

5.   Conclusion

In this work, we have proposed and studied a non-local macroscopic pedestrian flow model accounting for the presence of walls and obstacles limiting the walking domain. In particular, the proposed model is able to capture pedestrians' discomfort near obstacles and walls. Under suitable regularity assumptions, the model turns out to be well-posed. Moreover, we analyzed the impact of different anisotropic kernels on the formation of patterns in the solutions. High-resolution numerical schemes of WENO type allow to perform accurate simulations, bypassing the computational bottleneck given by the dependence of the flux function on integral terms. Future research should focus on multi-population models accounting for groups with different characteristics and/or destinations, and on the theoretical analysis of the observed pattern formation.

Acknowledgements

This research was supported by the INRIA Associated Team "Efficient numerical schemes for non-local transport phenomena" (NOLOCO; 2018–2020). RB is supported by Fondecyt project 1170473 and CRHIAM, project ANID/FONDAP/15130015. LMV is supported by Fondecyt project 1181511. RB, DI and LMV are also supported by CONICYT/PIA/Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal AFB170001.

Conflict of interest

The authors declare there is no conflicts of interest in this paper.