# Dynamic models for bird population-A parameter-varying partial differential equation identification approach

Régis Ouvrard, Guillaume Mercère, Thierry Poinot, Frédéric Jiguet, Lauriane Mouysset

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

## To cite this version:

Régis Ouvrard, Guillaume Mercère, Thierry Poinot, Frédéric Jiguet, Lauriane Mouysset. Dynamic models for bird population-A parameter-varying partial differential equation identification approach. Control Engineering Practice, 2019, 91, pp.104091. 10.1016/j.conengprac.2019.07.009. hal-02427943

HAL Id: hal-02427943

https://hal.science/hal-02427943v1

Submitted on 12 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Dynamic models for bird population – A parameter-varying partial differential equation identification approach

Régis Ouvrard<sup>a,*</sup>, Guillaume Mercère<sup>a</sup>, Thierry Poinot<sup>a</sup>, Frédéric Jiguet<sup>b</sup>, Lauriane Mouysset<sup>c</sup>

<sup>a</sup>Université de Poitiers, Laboratoire d'Informatique et d'Automatique pour les Systèmes, 2 rue Pierre Brousse, TSA 41105, 86073 Poitiers cedex 9, France

<sup>b</sup>Sorbonne Universités, Muséum National d'Histoire Naturelle, Centre d'Écologie et des Sciences de la Conservation, 43 rue Buffon, 75005 Paris, France

<sup>c</sup>Université de Bordeaux, Groupe de Recherche en Économie Théorique et Appliquée, avenue Léon Duguit, 33608 Pessac cedex, France

**Abstract**

In order to study the global decline of biodiversity, accurate models of animal population dynamics are required. In this paper, this challenging problem of biodiversity decline analysis is tackled by modelling the dynamics of bird populations. More specifically, a new data-driven modelling of bird population dynamics is suggested which resorts to a parameter-varying partial differential equation (PDE) model, the Galerkin method and the proper orthogonal decomposition. The parameter-varying formulation allows us to introduce in models prior information such as temperatures or landscape patterns. Hence, such models can be used to study the impact on biodiversity of the global warming or the agricultural intensification, for instance. In order to deal with the specific conditions of ecological applications, a specific attention is paid to the initialization as well as the implementation of an iterative identification procedure based on 3D partial moments and a Levenberg-Marquardt algorithm. These tools are tested on the data-driven modelling of the population of European Stonechat *Saxicola torquatus*, a European common bird species, by using data from the national French Breeding Bird Survey and the CORINE Land Cover. The perspectives are to model specific community of birds in order to evaluate the effect of global changes in population trends and to develop tools to help decision makers take into account biodiversity goals into public policies.

**Keywords:** Biodiversity, Ecology, Galerkin method, Parameter-varying, Partial differential equation, Population dynamics, Proper orthogonal decomposition, System identification

**1. Introduction**

The global decline of biodiversity is an undeniable fact today (Hallmann et al., 2017; IPBES, 2018a; WWF, 2018). As mentioned in the assessment reports of IPBES (2018b,c,d), in Americas, Central Asia or Europe, such a biodiversity loss is mainly due to the land degradation and, for instance, to an unsustainable agricultural. These landscape modifications have a strong impact on animal population dynamics. As a direct consequence of this correlation between land degradation and animal population dynamics, it is nowadays essential to have accurate models of animal population dynamics in order to analyse the biodiversity decline as well as to predict future trends and influence public policies. Different population dynamic models have been

developed in literature such as in (Siriwardena et al., 1998; Chamberlain et al., 2000; Mouysset et al., 2016). All these models, based on curve fitting approaches or ordinary differential equations (ODEs), describe overall population trends in an area or a country. At best, these models characterize the temporal dynamics of each patch which composes the considered area. Nowadays, new numerical tools allow us to model population dynamics with partial differential equations (PDEs) in order to characterize temporal and spatial dynamics simultaneously. To deal with this objective, automatic control community and more specifically system identification teams should investigate models based on PDEs.

Birds are known as the best taxon to measure the global biodiversity decline (Inger et al., 2014; Stanton et al., 2018). The developments introduced in

<footnote>*Corresponding author</footnote>
this paper are thus motivated by the necessary improve-
ments of the population dynamic models of birds in
order to understand past dynamics related to global
changes such as agricultural changes. The main per-
spective is to develop tools to help decision makers
take into account biodiversity goals into public poli-
cies. Thus, accurate models in terms of time and space
are needed. In (Mouysset, 2012; Mouysset et al., 2012,
2016), dynamic models for bird populations in France
have been estimated from large sets of common bird
species and accurate agricultural patterns. These mod-
els have been used to evaluate different agricultural
policy scenarios which promote targets of biodiversity
up to 2050. This is an example of tools which al-
low researchers to analyse the direct impacts of eco-
nomic scenarios of the Common Agricultural Policy
from European Union. The dynamic models estimated
in (Mouysset, 2012; Mouysset et al., 2012, 2016) are
based on the well-known logistic equation which is
often used to model a population growth in ecology
(Wu, 1996; Okubo and Levin, 2001). More specifi-
cally, these models are deduced from a nonlinear ODE,
the Verhulst's model, which characterizes the tempo-
ral dynamic. In the present paper, PDE models are
introduced to improve the modelling of the popula-
tion by adding spatial variables. Notice that, even if
no effective tools for data-driven parameter estimation
have been developed in the literature applied to ecol-
ogy so far (to the authors' knowledge), PDEs have been
used to represent population dynamics for many years
(Skellam, 1951; Levin, 1976; Shigesada et al., 1979;
Okubo, 1986; Holmes et al., 1994; Okubo and Levin,
2001). Okubo and Levin (2001) survey a wide variety
of PDEs to model dynamic behaviours of populations.
It is a way to describe spatial variations of a popula-
tion: for instance, the diffusion equation can be seen
as an extension of the random walk model of an indi-
vidual, the advection equation can represent the effect
of the non-homogeneities of environment producing an
oriented movement of animals, the logistic growth can
be used to introduce a temporal dynamic to characterize
increase or decrease of a population.

In order to introduce the landscape pattern in the
models to evaluate, for instance, the environmental foot-
print of agricultural intensification, linear parameter-
varying (LPV) PDEs are suggested. Although LPV
models have been widely studied for two decades (Toth,
2010), the LPV PDE model applications are rarely con-
sidered in literature. The following references can be
cited (Belforte et al., 2005; Schorsch et al., 2013, 2014;
Farah et al., 2016; Pham et al., 2018). In the sequel, the
term parameter-varying is favoured because the consid-

ered PDE models can be linear or nonlinear.

System identification consists in estimating the pa-
rameters of a dynamic model by adjusting these parame-
ters to minimize the discrepancy between measured data
and system behaviour. Today, the system identification
theory is widely developed for linear time-invariant sys-
tems (Ljung, 1999), i.e. for systems whose behaviour
can be described by ODEs. System identification for
distributed-parameter systems, i.e. for systems whose
behaviour can be described by PDEs, is more com-
plicated. Measurements and spatial distribution of the
sensors are critical points to achieve accurate results
(Uciński, 2005). In this paper, the considered data-sets
are given by the national French Breeding Bird Survey
described in (Jiguet et al., 2012). The sampling time is
limited to one data by breeding season, i.e. one sam-
ple per year at each location, during thirteen years. This
sampling limitation enforces us to develop new iden-
tification methods to estimate PDE parameters. More
specifically, an iterative approach based on time-space
separation (Li and Qi, 2010) is suggested. In order to
bypass some difficulties encountered when PDEs come
into play, this iterative approach consists in combining
the Galerkin method (Polis et al., 1973) and the proper
orthogonal decomposition (Newman, 1996a,b), respec-
tively. This iterative procedure requires an initialisa-
tion which is given by an estimator based on 3D partial
moments, an extension of partial moments described in
(Ouvrard and Trigeassou, 2011). These new tools are
validated in simulation before being used on real data
from the national French Breeding Bird Survey.

The present paper is a preliminary work which aims
at validating the approach by specialists of system iden-
tification. As mentioned in perspectives in the conclu-
sion section, a future work will involve specific commu-
nities of birds to evaluate the effect of global changes in
population declines.

The paper is organized as follows. Section 2 de-
scribes the parameter-varying PDE model while the
POD-Galerkin method is presented in Section 3 for a
general formulation of the PDE. The identification pro-
cedure is introduced in Section 4. Section 5 is dedicated
to the validation of the identification approach with sim-
ulated data and a simplified PDE. An application to the
population of European Stonechat *Saxicola torquatus* is
shown in Section 6 with parameters varying with the
landscape pattern in France. Section 7 concludes this
paper.
## 2. Parameter-varying PDE model

The goal is to present an identification procedure for a general formulation of a parameter-varying PDE model, i.e. with diffusion, advection and logistic growth, and with parameter variations given by polynomial functions. Thus, the following parameter-varying PDE model is considered

$$ \frac{\partial u(x, y, t)}{\partial t} = D_x(H) \frac{\partial^2 u(x, y, t)}{\partial x^2} + D_y(H) \frac{\partial^2 u(x, y, t)}{\partial y^2} - w_x(H) \frac{\partial u(x, y, t)}{\partial x} - w_y(H) \frac{\partial u(x, y, t)}{\partial y} + \beta_1(H)u(x, y, t) - \beta_2(H)u^2(x, y, t), \quad (1) $$

where x and y are spatial coordinates, t is time, u(x, y, t) is a population, H(x, y) is the so-called scheduling variable (Toth, 2010) while the varying parameters are given by

$$ \begin{aligned} D_{\bullet}(H) &= \sum_{n=0}^{n_{D_{\bullet}}} D_{\bullet, n} H^n(x, y), \\ w_{\bullet}(H) &= \sum_{n=0}^{n_{w_{\bullet}}} w_{\bullet, n} H^n(x, y), \quad (2) \\ \beta_{\bullet}(H) &= \sum_{n=0}^{n_{\beta_{\bullet}}} \beta_{\bullet, n} H^n(x, y), \end{aligned} $$

with $n_{D_{\bullet}}$, $n_{w_{\bullet}}$ and $n_{\beta_{\bullet}}$, the polynomial orders, and $D_{\bullet, n}$, $w_{\bullet, n}$ and $\beta_{\bullet, n}$, the parameters to be estimated.

The diffusion coefficients $D_x$ and $D_y$ stand for the dispersal rate (Holmes et al., 1994). The advection term relies on $w_x$ and $w_y$, the drift velocities, which characterize oriented motions due to an external stimuli or an attractive area, for instance. $\beta_1$ and $\beta_2$ are parameters which yield to a logistic growth of the population. $\beta_1$ is the intrinsic growth rate and $\beta_2/\beta_1$ represents the carrying capacity. The scheduling variable $H(x, y)$ allows us to take into account, for example, the land use to analyse agricultural changes to animal populations or the temperature to study effects of global warming.

Notice that most of the PDE models described in (Fisher, 1937; Levin, 1976; Holmes et al., 1994; Wu, 1996; Okubo and Levin, 2001) are special cases of the parameter-varying PDE model (1).

The main issue is now to estimate the parameters of polynomials (2) from data-sets.

## 3. POD-Galerkin method

The parameter estimation of the PDE (1) requires to use an optimization algorithm, i.e. an iterative procedure which minimizes a cost function based on the

measured data and the PDE simulation. But the simulation of a parameter-varying PDE such as (1) is a difficult issue. Therefore, it is preferable to calculate an approximate solution of the PDE by using the so-called POD-Galerkin method (Newman, 1996a,b) and the time-space separation. More specifically, a Galerkin method is considered and based on function basis which is built with a proper orthogonal decomposition (POD).

### 3.1. Time-space separation basic idea

The Hilbert space theory (Debnath and Mikusinski, 2005) proves that any arbitrary signal $f$ belonging to a Hilbert space like $L_2(D)$ (the space of Lebesgue-measurable square-integrable functions on a domain $D$, i.e. $\int_D |f(x)|^2 dx < \infty$), with $\langle \bullet, \bullet \rangle$ its inner product, can be decomposed into a Fourier expansion of the form

$$ f = \sum_{n=1}^{\infty} \langle f, \phi_n \rangle \phi_n, \quad (3) $$

once an orthonormal basis $\{\phi_n(x)\}_{n=1}^{\infty}$ of $L_2(D)$ has been selected. More specifically, for signals like $u(x, y, t)$ belonging to the Hilbert space $L_2(D)$, for $D = \{x \in [0, X] \text{ and } y \in [0, Y]\}$, with the inner product

$$ \langle f, g \rangle = \int_{D} \int_{D} f(x, y, t)g(x, y, t)dx dy, \quad (4) $$

the signal $u(x, y, t)$ satisfies

$$ u(x, y, t) = \sum_{n=1}^{\infty} a_n(t) \phi_n(x, y), \quad (5) $$

where $a_n(t) = \langle u, \phi_n \rangle$ are called the Fourier coefficients of this expansion. Eq. (5) is the basic idea to approximate the spatio-temporal signals involved in the population dynamics via the well-known time-space separation approach (Li and Qi, 2010).

### 3.2. Proper orthogonal decomposition

In this paper, in order to benefit from available datasets, it was decided to use the proper orthogonal decomposition to generate data-based orthonormal basis functions $\{\phi_n(x, y)\}_{n=1}^{\infty}$. As explained hereafter, such a technique

* does not require specific prior knowledge on the way the system is excited (which is the case, e.g., when basis functions like cosine and sine functions are introduced),

* helps the user to select the truncation order by looking at the magnitude of singular or eigenvalues generated by the POD.
Also called Karhunen-Loeve decomposition or principal component analysis (PCA) (Liang et al., 2002), the POD consists in generating an empirical covariance data matrix, then determining the orthonormal basis {$\phi_n(x, y)$}$_{n=1}^N$ by extracting the main vectors spanning its column space by using an eigenvalue or singular value decomposition (Newman, 1996a,b). The empirical sampled covariance matrix used by the POD method is more precisely defined by

$$ \mathbf{K} = \begin{bmatrix} \alpha_{11} & \cdots & \alpha_{1M} \\ \vdots & \ddots & \vdots \\ \alpha_{M1} & \cdots & \alpha_{MM} \end{bmatrix}, \quad (6) $$

where the components $\alpha_{kl}$ are determined by considering the inner product between time samples as follows<fnref n="1" />

$$ \alpha_{kl} = \frac{1}{M} \int \int_{\mathcal{D}} \delta(x, y, k \, dt) \delta(x, y, \ell \, dt) \, dx dy, \quad (7) $$

and from the so-called M zero-mean snapshots

$$ \delta(x, y, k \, dt) = u(x, y, k \, dt) - \bar{u}(x, y), \quad (8) $$

with $\bar{u}(x, y) = \text{mean}_t(u(x, y, t))$, the time average of the population.

**Remark 1.** For the sake of conciseness, in the sequel, the double integration $\int \int_{\mathcal{D}} \cdot dx dy$ introduced with the inner product will be denoted $\int \cdot$ and, for instance, Eq. (7) is now denoted

$$ \alpha_{kl} = \frac{1}{M} \int \delta(x, y, k \, dt) \delta(x, y, \ell \, dt). \quad (9) $$

Once the matrix **K** is generated, the determination of the basis functions {$\phi_n(x, y)$}$_{n=1}^\infty$ is performed by computing the eigenvalue decomposition of the real symmetric matrix **K**, feature which guarantees that the eigenvalues are all real and the eigenvectors orthogonal to each other. More specifically,

$$ \mathbf{K} = \mathbf{Q} \mathbf{R} \mathbf{Q}^{\top}, \quad (10) $$

with $\mathbf{Q} \in \mathbb{R}^{M \times M}$ with $\mathbf{Q}\mathbf{Q}^{\top} = \mathbf{I}_M$, a $M \times M$ identity matrix, while $\mathbf{R} \in \mathbb{R}^{M \times M}$ contains the eigenvalues of **K** assumed to be ordered in a decreasing order. By using the direct link between the POD and the PCA, it is clear that the eigenvectors corresponding to the largest eigenvalues are the main modes of the covariance matrix **K**. Said differently, by selecting the basis functions

{$\phi_n(x, y)$}$_{n=1}^N$ from the first N eigenvectors of **K** with the largest eigenvalues (i.e., the largest energy among all the available components), it is guaranteed that

$$ \sum_{n=1}^{N} a_n(t) \phi_n(x, y) \quad (11) $$

is the best approximation of $u(x, y, t)$ (in the least-squares sense) given {$\phi_n(x, y)$}$_{n=1}^N$ (Liang et al., 2002). Finally, the basis function set {$\phi_n(x, y)$}$_{n=1}^N$ is determined as a linear combination of the data snapshots

$$ \phi_n(x, y) = \sum_{k=1}^{M} \mathbf{Q}_k^{(n)} \delta(x, y, k \, dt), \quad (12) $$

where $\mathbf{Q}_k^{(n)}$ is the k-th component of the n-th eigenvector of **K**. These computations are performed one time as an initial step of the iterative identification procedure.

### 3.3. Galerkin method

As previously mentioned, the Galerkin method operates a time-space separation which transforms the PDE problem into a system of ODEs. Some details of the calculations are given in Appendix A. It consists in finding an approximate solution of the PDE defined by

$$ \hat{u}(x, y, t) = \sum_{n=1}^{N} a_n(t) \phi_n(x, y) + \bar{u}(x, y), \quad (13) $$

where {$\phi_n(x, y)$}$_{n=1}^N$ is the known function basis determined by POD, and {$a_n(t)$}$_{n=1}^N$ is the time coefficient set to be estimated.

<footnote n="1">¹dt is the time sampling period.</footnote>
The substitution of (13) into (1) leads to

$$
\begin{align}
\sum_{n=1}^{N} \dot{a}_n(t) \phi_n(x, y) = & \sum_{n=1}^{N} a_n(t) \left( D_x(H) \frac{\partial^2 \phi_n(x, y)}{\partial x^2} + D_y(H) \frac{\partial^2 \phi_n(x, y)}{\partial y^2} \right. \nonumber \\
& \left. - w_x(H) \frac{\partial \phi_n(x, y)}{\partial x} - w_y(H) \frac{\partial \phi_n(x, y)}{\partial y} \right. \nonumber \\
& \left. + \beta_1(H) \phi_n(x, y) \right) \nonumber \\
& - 2\beta_2(H) \sum_{n=1}^{N} a_n(t) \phi_n(x, y) \bar{u}(x, y) \tag{14} \\
& - \beta_2(H) \sum_{n=1}^{N} a_n^2(t) \phi_n^2(x, y) \nonumber \\
& - 2\beta_2(H) \sum_{n=1}^{N-1} \sum_{m=n+1}^{N} a_n(t) a_m(t) \phi_n(x, y) \phi_m(x, y) \nonumber \\
& + D_x(H) \frac{\partial^2 \bar{u}(x, y)}{\partial x^2} + D_y(H) \frac{\partial^2 \bar{u}(x, y)}{\partial y^2} \nonumber \\
& - w_x(H) \frac{\partial \bar{u}(x, y)}{\partial x} - w_y(H) \frac{\partial \bar{u}(x, y)}{\partial y} \nonumber \\
& \left. + \beta_1(H) \bar{u}(x, y) - \beta_2(H) \bar{u}^2(x, y) \right) \nonumber
\end{align}
$$

By applying the inner product $\langle \text{Eq.}(14), \phi_\ell \rangle$, the or-
thonormal property of $\{\phi_n(x, y)\}_{n=1}^N$ yields the following
system of ODEs for $\ell = 1, \dots, N$

$$
\begin{align}
\dot{a}_\ell(t) = & \sum_{n=1}^{N} a_n(t) \int \left\{ D_x(H) \frac{\partial^2 \phi_n}{\partial x^2} + D_y(H) \frac{\partial^2 \phi_n}{\partial y^2} \right. \nonumber \\
& \left. - w_x(H) \frac{\partial \phi_n}{\partial x} - w_y(H) \frac{\partial \phi_n}{\partial y} + \beta_1(H) \phi_n \right\} \phi_\ell \nonumber \\
& - 2 \sum_{n=1}^{N} a_n(t) \int \beta_2(H) \phi_n \phi_\ell \bar{u} \nonumber \\
& - \sum_{n=1}^{N} a_n^2(t) \int \beta_2(H) \phi_n^2 \phi_\ell \nonumber \\
& - 2 \sum_{n=1}^{N-1} \sum_{m=n+1}^{N} a_n(t) a_m(t) \int \beta_2(H) \phi_n \phi_m \phi_\ell \nonumber \\
& + \int \left\{ D_x(H) \frac{\partial^2 \bar{u}}{\partial x^2} + D_y(H) \frac{\partial^2 \bar{u}}{\partial y^2} - w_x(H) \frac{\partial \bar{u}}{\partial x} \right. \nonumber \\
& \left. - w_y(H) \frac{\partial \bar{u}}{\partial y} + \beta_1(H) \bar{u} - \beta_2(H) \bar{u}^2 \right\} \phi_\ell \tag{15}
\end{align}
$$

For an easier implementation, this system of ODEs can
be rewritten with a matrix formulation

$$
\dot{\mathbf{a}}(t) = (\mathbf{\Gamma}_1 + \beta_{10} \mathbf{I}_N - 2\mathbf{\Gamma}_2) \mathbf{a}(t) \\
- (\mathbf{I}_N \otimes \mathbf{a}^{\top}(t)) \mathbf{\Lambda} \mathbf{a}(t) + \mathbf{b}_1, \quad (16)
$$

with

⊗, Kronecker product,

$$
\mathbf{a}(t) = \left[ \begin{array}{c} a_1(t) \quad \dots \quad a_N(t) \end{array} \right]^\top,
\mathbf{I}_N, N \times N \text{ identity matrix,}
$$

$$
\Gamma_1 = \begin{bmatrix} \int \mathcal{F}_1(\phi_1)\phi_1 \dots \int \mathcal{F}_1(\phi_N)\phi_1 \\ \vdots & \ddots & \vdots \\ \int \mathcal{F}_1(\phi_1)\phi_N \dots \int \mathcal{F}_1(\phi_N)\phi_N \end{bmatrix},
$$
$$
\Gamma_2 = \begin{bmatrix} \int \beta_2(H)\phi_1\phi_1\bar{u} \dots \int \beta_2(H)\phi_N\phi_1\bar{u} \\ \vdots & \ddots & \vdots \\ \int \beta_2(H)\phi_1\phi_N\bar{u} \dots \int \beta_2(H)\phi_N\phi_N\bar{u} \end{bmatrix}, \quad (17)
$$

$$
\mathbf{\Lambda} = \int \beta_2(H) \mathbf{\Phi} \otimes (\mathbf{\Phi} \otimes \mathbf{\Phi}^{\top}),
$$

$$
\mathbf{\Phi} = \left[ \begin{array}{c} \phi_1(x, y) \quad \dots \quad \phi_N(x, y) \end{array} \right]^\top,
$$

$$
\mathbf{b}_1 = \left[ \int \mathcal{F}_2(\bar{u})\phi_1 \dots \int \mathcal{F}_2(\bar{u})\phi_N \right]^\top,
$$

$$
\begin{align*}
\mathcal{F}_1(\bullet) &= D_x(H) \frac{\partial^2 \bullet}{\partial x^2} + D_y(H) \frac{\partial^2 \bullet}{\partial y^2} - w_x(H) \frac{\partial \bullet}{\partial x} \\
&\quad - w_y(H) \frac{\partial \bullet}{\partial y} + \sum_{\nu=1}^{n_{\beta_1}} \beta_{1\nu} H^{\nu}(x, y) \bullet,
\end{align*}
$$

$$
\begin{align*}
\mathcal{F}_2(\bullet) = D_x(H) \frac{\partial^2 \bullet}{\partial x^2} + D_y(H) \frac{\partial^2 \bullet}{\partial y^2} - w_x(H) \frac{\partial \bullet}{\partial x} \\
&\quad - w_y(H) \frac{\partial \bullet}{\partial y} + \beta_1(H) \bullet - \beta_2(H) \bullet^2.
\end{align*}
$$

The initial conditions of the ODEs system (16) is ob-
tained by making the initial signal error orthogonal to
the basis functions $\phi_\ell(x, y)$, i.e.

$$
\langle u(\cdot, \cdot, 0) - \hat{u}(\cdot, \cdot, 0), \phi_\ell \rangle_{\ell=1,\dots,N} = 0, \quad (18)
$$

or

$$
a_{\ell}(0) = \int (u(x, y, 0) - \bar{u}(x, y)) \phi_{\ell}(x, y), \quad (19)
$$

for $\ell = 1, \dots, N$.

Therefore, given parameters, the simulation of the
system of ODEs (16) from the initial conditions (19)
provides the approximate solution (13). The latter can
be used in the following iterative identification proce-
dure.

4. Iterative identification procedure

4.1. Parameter estimation

The parameter vector $\boldsymbol{\theta} \in \mathbb{R}^{N_{\theta}}$ given by

$$
\boldsymbol{\theta} = \begin{bmatrix} D_{x_0} \dots D_{x_{n_{D_x}}} & D_{y_0} \dots D_{y_{n_{D_y}}} & w_{x_0} \dots w_{x_{n_{w_x}}} \\ w_{y_0} \dots w_{y_{n_{w_y}}} & \beta_{1_0} \dots \beta_{1_{\beta_1}} & \beta_{2_0} \dots \beta_{2_{\beta_2}} \end{bmatrix}^{\top} \quad (20)
$$

is estimated from a measured population<fnref n="2" />
${\{u(i, j, k)\}_{i=0,...,N_s; \ j=0,...,N_y; \ k=0,...,N_t}}$ with an iterative

<footnote n="2">²For simplicity, $u(i, j, k)$ denotes sampled data $u(i \, dx, j \, dy, k \, dt)$ with corresponding sampling distances $dx, dy$ and $dt$, respectively.</footnote>
procedure by minimizing the quadratic criterion

$$J(\theta_{iter}) = \sum_{k=0}^{N_t} \sum_{j=0}^{N_y} \sum_{i=0}^{N_x} (u(i, j, k) - \hat{u}(i, j, k))^2, \quad (21)$$

where $\hat{u}(i, j, k)$ is the approximate solution (13) of the PDE model (1), and $\theta_{iter}$ is the parameter vector at the iter-th iteration. This procedure repeats iteratively the following two steps from an initial vector $\theta_0$ and until convergence

* approximate a solution $\hat{u}(i, j, k)$ by the Galerkin method, Eq. (13),

* compute a new parameter vector $\theta_{iter+1}$ with a Levenberg-Marquardt algorithm which minimizes the criterion (21).

## 4.2. Levenberg-Marquardt algorithm

Because the considered PDE model (1) is non-linear in parameters, the Levenberg-Marquardt algorithm is used to identify its parameters (Björck, 1996; Nocedal and Wright, 2006). This model does not have exogenous input, and the unique excitation is the initial condition $u(x, y, 0)$. Moreover, the ecological dataset often has a small number of time samples. Thus, in order to limit the approximation error of the optimization algorithm, the true parametric sensitivity functions (Knudsen, 1994) are considered with initial conditions deduced from (19).

At the iter-th iteration, the new parameter vector is given by

$$\theta_{iter+1} = \theta_{iter} - (\mathbf{J}''_{\theta\theta} + \mu \mathbf{I}_{N_\theta})^{-1} \mathbf{J}'_{\theta}, \quad (22)$$

with $\mu$, a tuning parameter, $\mathbf{I}_{N_\theta}$, an identity matrix with appropriate dimensions, and the gradient and the Hessian respectively defined by

$$\mathbf{J}'_{\theta} = -2 \sum_{k=0}^{N_t} \sum_{j=0}^{N_y} \sum_{i=0}^{N_x} (u(i, j, k) - \hat{u}(i, j, k)) \Xi(k) \Phi(i, j),$$

$$\mathbf{J}''_{\theta\theta} = 2 \sum_{k=0}^{N_t} \sum_{j=0}^{N_y} \sum_{i=0}^{N_x} \Xi(k) \Phi(i, j) \Phi^{\top}(i, j) \Xi^{\top}(k), \quad (23)$$

where $\hat{u}(i, j, k)$ is the approximate solution (13), $\Phi(x, y)$ is defined in (17), and $\Xi(t)$ is the $N_\theta \times N$ time sensitivity

function matrix given by

$$\Xi(t) = \begin{bmatrix} \sigma_{D_{x0}}^{a}(t) & \cdots & \sigma_{D_{xN_p}}^{a}(t) & \sigma_{D_{y0}}^{a}(t) \\ \cdots & \sigma_{D_{xN_p}}^{a}(t) & \sigma_{w_{x0}}^{a}(t) & \cdots & \sigma_{w_{xN_x}}^{a}(t) \\ \sigma_{w_{y0}}^{a}(t) & \cdots & \sigma_{w_{yN_y}}^{a}(t) & \sigma_{\beta_{10}}^{a}(t) & \cdots \\ \sigma_{\beta_{1\beta_1}}^{a}(t) & \sigma_{\beta_{20}}^{a}(t) & \cdots & \sigma_{\beta_{2\beta_2}}^{a}(t) \end{bmatrix}^{\top}, \quad (24) \\ \sigma_{D_{\bullet n}}^{a}(t) = \begin{bmatrix} \sigma_{D_{\bullet n}}^{a1}(t) & \cdots & \sigma_{D_{\bullet n}}^{aN}(t) \end{bmatrix}^{\top}, \\ \sigma_{w_{\bullet n}}^{a}(t) = \begin{bmatrix} \sigma_{w_{\bullet n}}^{a1}(t) & \cdots & \sigma_{w_{\bullet n}}^{aN}(t) \end{bmatrix}^{\top}, \\ \sigma_{\beta_{\bullet n}}^{a}(t) = \begin{bmatrix} \sigma_{\beta_{\bullet n}}^{a1}(t) & \cdots & \sigma_{\beta_{\bullet n}}^{aN}(t) \end{bmatrix}^{\top}, \\ \sigma_{\circ}^{a_m}(t) = \frac{\partial a_m(t)}{\partial \circ}.$$

The parametric sensitivity functions are simulated by considering the following matrix formulations deduced from derivatives of (16) with respect to each parameter

$$\begin{align} \dot{\sigma}_{D_{\bullet n}}^{a}(t) &= (\Gamma_1 + \beta_{10} \mathbf{I}_N - 2\Gamma_2) \sigma_{D_{\bullet n}}^{a}(t) \nonumber \\ &\quad + \Gamma_{D_{\bullet n}} \mathbf{a}(t) - (\mathbf{I}_N \otimes \sigma_{D_{\bullet n}}^{a\top}(t)) \mathbf{\Lambda} \mathbf{a}(t) \nonumber \\ &\quad - (\mathbf{I}_N \otimes \mathbf{a}^{\top}(t)) \mathbf{\Lambda} \sigma_{D_{\bullet n}}^{a}(t) + \mathbf{b}_{D_{\bullet n}}, \nonumber \\ \dot{\sigma}_{w_{\bullet n}}^{a}(t) &= (\Gamma_1 + \beta_{10} \mathbf{I}_N - 2\Gamma_2) \sigma_{w_{\bullet n}}^{a}(t) \nonumber \\ &\quad + \Gamma_{w_{\bullet n}} \mathbf{a}(t) - (\mathbf{I}_N \otimes \sigma_{w_{\bullet n}}^{a\top}(t)) \mathbf{\Lambda} \mathbf{a}(t) \nonumber \\ &\quad - (\mathbf{I}_N \otimes \mathbf{a}^{\top}(t)) \mathbf{\Lambda} \sigma_{w_{\bullet n}}^{a}(t) + \mathbf{b}_{w_{\bullet n}}, \nonumber \\ \dot{\sigma}_{\beta_{1n}}^{a}(t) &= (\Gamma_1 + \beta_{10} \mathbf{I}_N - 2\Gamma_2) \sigma_{\beta_{1n}}^{a}(t) \nonumber \\ &\quad + \Gamma_{\beta_{1n}} \mathbf{a}(t) - (\mathbf{I}_N \otimes \sigma_{\beta_{1n}}^{a\top}(t)) \mathbf{\Lambda} \mathbf{a}(t) \nonumber \\ &\quad - (\mathbf{I}_N \otimes \mathbf{a}^{\top}(t)) \mathbf{\Lambda} \sigma_{\beta_{1n}}^{a}(t) + \mathbf{b}_{\beta_{1n}}, \nonumber \\ \dot{\sigma}_{\beta_{2n}}^{a}(t) &= (\Gamma_1 + \beta_{10} \mathbf{I}_N - 2\Gamma_2) \sigma_{\beta_{2n}}^{a}(t) \nonumber \\ &\quad - 2\Gamma_{\beta_{2n}} \mathbf{a}(t) - (\mathbf{I}_N \otimes \sigma_{\beta_{2n}}^{a\top}(t)) \mathbf{\Lambda} \mathbf{a}(t) \nonumber \\ &\quad - (\mathbf{I}_N \otimes \mathbf{a}^{\top}(t)) \mathbf{\Lambda} \sigma_{\beta_{2n}}^{a}(t) \nonumber \\ &\quad - (\mathbf{I}_N \otimes \mathbf{a}^{\top}(t)) \mathbf{\Lambda} \sigma_{\beta_{2n}}^{a}(t) - \mathbf{b}_{\beta_{2n}}, \end{align} \quad (25)$$
for $n = 0, \dots, n_{D_*}, n_{w_*}$ or $n_{\beta_*}$ and where

$$
\begin{align*}
\Gamma_{D_{\bullet n}} &= \left[ \begin{array}{cccc}
\int \mathcal{F}_{D_{\bullet n}}(\phi_1)\phi_1 & \dots & \int \mathcal{F}_{D_{\bullet n}}(\phi_N)\phi_1 \\
\vdots & \ddots & \vdots \\
\int \mathcal{F}_{D_{\bullet n}}(\phi_1)\phi_N & \dots & \int \mathcal{F}_{D_{\bullet n}}(\phi_N)\phi_N
\end{array} \right], \\
\Gamma_{w_{\bullet n}} &= \left[ \begin{array}{cccc}
\int \mathcal{F}_{w_{\bullet n}}(\phi_1)\phi_1 & \dots & \int \mathcal{F}_{w_{\bullet n}}(\phi_N)\phi_1 \\
\vdots & \ddots & \vdots \\
\int \mathcal{F}_{w_{\bullet n}}(\phi_1)\phi_N & \dots & \int \mathcal{F}_{w_{\bullet n}}(\phi_N)\phi_N
\end{array} \right], \\
\Gamma_{\beta_{10}} &= \mathbf{I}_N,
\end{align*}
$$

$$
\Gamma_{\beta_{1n}} = \begin{bmatrix} \int H^n \phi_1 \phi_1 \dots \int H^n \phi_N \phi_1 \\ \vdots \vdots \vdots \vdots \\ \int H^n \phi_1 \phi_N \dots \int H^n \phi_N \phi_N \end{bmatrix},
$$

$$
\begin{align}
\Gamma_{\beta_{2n}} &= \begin{bmatrix} \int H^n \phi_1 \phi_1 \bar{u} \dots \int H^n \phi_N \phi_1 \bar{u} \\ \vdots \vdots \vdots \vdots \\ \int H^n \phi_1 \phi_N \bar{u} \dots \int H^n \phi_N \phi_N \bar{u} \end{bmatrix}, \nonumber \\
\Lambda_{\beta_{2n}} &= \int H^n \Phi \otimes (\Phi \otimes \Phi^\top) \nonumber \\
\mathbf{b}_{D_{\bullet n}} &= \left[ \int \mathcal{F}_{D_{\bullet n}}(\bar{u})\phi_1 \dots \int \mathcal{F}_{D_{\bullet n}}(\bar{u})\phi_N \right]^\top, \nonumber \\
\mathbf{b}_{w_{\bullet n}} &= \left[ \int \mathcal{F}_{w_{\bullet n}}(\bar{u})\phi_1 \dots \int \mathcal{F}_{w_{\bullet n}}(\bar{u})\phi_N \right]^\top, \nonumber \\
\mathbf{b}_{\beta_{1n}} &= \left[ \int H^n \bar{u} \phi_1 \dots \int H^n \bar{u} \phi_N \right]^\top, \nonumber \\
\mathbf{b}_{\beta_{2n}} &= \left[ \int H^n \bar{u}^2 \phi_1 \dots \int H^n \bar{u}^2 \phi_N \right]^\top, \nonumber \\
\mathcal{F}_{D_{\bullet n}}(\circ) &= H^n(x, y) \frac{\partial^2 \circ}{\partial \bullet^2}, \quad \mathcal{F}_{w_{\bullet n}}(\circ) = H^n(x, y) \frac{\partial \circ}{\partial \bullet}, \tag{26}
\end{align}
$$

and with the initial conditions

$$
\begin{align}
\sigma_{D_{\bullet n}}^{a}(0) &= \Gamma_{D_{\bullet n}} \mathbf{a}(0) + \mathbf{b}_{D_{\bullet n}}, \nonumber \\
\sigma_{w_{\bullet n}}^{a}(0) &= \Gamma_{w_{\bullet n}} \mathbf{a}(0) + \mathbf{b}_{w_{\bullet n}}, \nonumber \\
\sigma_{\beta_{1n}}^{a}(0) &= \Gamma_{\beta_{1n}} \mathbf{a}(0) + \mathbf{b}_{\beta_{1n}}, \tag{27} \\
\sigma_{\beta_{2n}}^{a}(0) &= -2\Gamma_{2} \mathbf{a}(0) \nonumber \\
&\quad - (\mathbf{I}_N \otimes \mathbf{a}^{\top}(0)) \Lambda_{\beta_{2n}} \mathbf{a}(0) - \mathbf{b}_{\beta_{2n}}. \nonumber
\end{align}
$$

4.3. Initial estimate based on 3D partial moments

The aim of this section is to obtain an initial estimate
$\theta_0$ of the parameters for the initialisation of the itera-
tive identification procedure given in (22). For the PDE
model, the choice of the initialisation is particularly cru-
cial to tend toward the global optimum of the crite-
rion. It can be done from prior knowledge. In ecolog-
ical applications, such a prior knowledge is often miss-
ing. Therefore, a least-squares method is suggested. It
is based on partial moments developed in (Trigeassou,
1987; Ouvrard and Trigeassou, 2011) for systems repre-
sented by ODEs and extended to the PDE models there-
after.

Define the 3D partial moment with orders $\ell$, $m$ and $n$ of an arbitrary signal $f(x, y, t)$ belonging to a Hilbert space by

$$
\begin{equation}
\begin{split}
\mathcal{M}_{\ell,m,n}^{f}(X, Y, T) &= \int_{0}^{X} \int_{0}^{Y} \int_{0}^{T} x^{\ell} y^{m} t^{n} f(x, y, t) \, dt \, dy \, dx, \\
&\text{with } \ell, m \text{ and } n \in \{\mathbb{N}, -\}^{3}.
\end{split}
\tag{28}
\end{equation}
$$

The idea of the partial moment-based approach is to
apply to (1) with $n_{D_{\bullet}} = 0$, $n_{w_{\bullet}} = 0$ and $n_{\beta_{\bullet}} = 0$ the
following integration

$$
\int_0^X \int_0^{x_1} x_2^2 \int_0^Y \int_0^{y_1} y_2^2 \int_0^T t \quad (1) \quad dt \, dy_2 \, dy_1 \, dx_2 \, dx_1
$$

in such a way that all partial derivatives disappear. Thus,
a formulation with only partial moments and parameters
is obtained. Indeed, by considering the integration by
parts

$$
\int_0^{\tau_0} \tau_1 \frac{\partial f}{\partial \tau} d\tau_1 = \tau_0 f(\tau_0) - \int_0^{\tau_0} f(\tau_1) d\tau_1 \quad (30)
$$

or

$$
\int_0^{\tau_0} \int_0^{\tau_1} \tau_2 \frac{\partial^2 f}{\partial \tau^2} d\tau_2 d\tau_1 = \tau_0^2 f(\tau_0) \\
+ 2\tau_0 \int_0^{\tau_0} f(\tau_1) d\tau_1 - 6 \int_0^{\tau_0} \tau_1 f(\tau_1) d\tau_1, \tag{31}
$$

and the Cauchy formula of repeated integrations

$$
\int_0^t \int_0^{\tau_1} \dots \int_0^{\tau_{n-1}} f(\tau_n) d\tau_n \dots d\tau_2 d\tau_1 = \frac{1}{(n-1)!} \int_0^t (t-\tau_1)^{n-1} f(\tau_1) d\tau_1, \quad (32)
$$

Eq. (29) can be reformulated with 3D partial moments
as follows

$$
\begin{equation}
\begin{split}
\mathcal{U}_1(X, Y, T) ={}& D_{x_0} \mathcal{U}_2(X, Y, T) + D_{y_0} \mathcal{U}_3(X, Y, T) \\
& - w_{x_0} \mathcal{U}_4(X, Y, T) - w_{y_0} \mathcal{U}_5(X, Y, T) \\
& + \beta_{10} \mathcal{U}_6(X, Y, T) + \beta_{20} \mathcal{U}_7(X, Y, T),
\end{split}
\tag{33}
\end{equation}
$$

<footnote>$^3$If $\ell$, $m$ or $n$ is the symbol $-$, then there is no integration with respect to the corresponding variable. For instance, $M_{\ell,m,-}^f(X, Y, T) = \int_0^X \int_0^Y x^\ell y^m f(x, y, T) \, dy \, dx$</footnote>
where

$$
\begin{align*}
\mathcal{U}_1 &= XYT \mathcal{M}_{2,2,-}^{u} - XT \mathcal{M}_{2,3,-}^{u} - YTM_{2,3,-}^{u} \\
&\quad + T \mathcal{M}_{3,3,-}^{u} - XYM_{2,2,0}^{u} + XM_{2,3,0}^{u} \\
&\quad + YM_{3,2,0}^{u} - \mathcal{M}_{3,3,0}^{u},
\end{align*}
$$

$$
\mathcal{U}_2 = X^2 Y \mathcal{M}_{-2,1}^{u} - X^2 \mathcal{M}_{-3,1}^{u} + 2XY \mathcal{M}_{0,2,1}^{u} \\ - 2X \mathcal{M}_{0,3,1}^{u} - 6Y \mathcal{M}_{1,2,1}^{u} + 6 \mathcal{M}_{1,3,1}^{u},
$$

$$
\mathcal{U}_3 = XY^2 \mathcal{M}_{2,-,1}^{u} - Y^2 \mathcal{M}_{3,-,1}^{u} + 2XY \mathcal{M}_{2,0,1}^{u} \\ - 2Y \mathcal{M}_{3,0,1}^{u} - 6X \mathcal{M}_{2,1,1}^{u} + 6 \mathcal{M}_{3,1,1}^{u},
$$

$$
\mathcal{U}_4 = 3Y \mathcal{M}_{2,2,1}^{u} - 3 \mathcal{M}_{2,3,1}^{u} - 2XY \mathcal{M}_{1,2,1}^{u} + 2X \mathcal{M}_{1,3,1}^{u},
$$

$$
\mathcal{U}_5 = 3X \mathcal{M}_{2,2,1}^{u} - 3 \mathcal{M}_{3,2,1}^{u} - 2XY \mathcal{M}_{2,1,1}^{u} + 2Y \mathcal{M}_{3,1,1}^{u},
$$

$$
\mathcal{U}_6 = XY \mathcal{M}_{2,2,1}^{u} - Y \mathcal{M}_{3,2,1}^{u} - X \mathcal{M}_{2,3,1}^{u} + \mathcal{M}_{3,3,1}^{u},
$$

$$
\mathcal{U}_7 = XY \mathcal{M}_{2,2,1}^{u2} - Y \mathcal{M}_{3,2,1}^{u2} - X \mathcal{M}_{2,3,1}^{u2} + \mathcal{M}_{3,3,1}^{u2}.
$$

Hence, (33) can be rewritten in a linear regression
form

$$
\mathcal{U}_1(X, Y, T) = \varphi^\top(X, Y, T) \theta_{LS}, \quad (34)
$$

with

$$
\varphi(X, Y, T) = \begin{bmatrix} \mathcal{U}_2(X, Y, T) & \mathcal{U}_3(X, Y, T) \\ -\mathcal{U}_4(X, Y, T) & -\mathcal{U}_5(X, Y, T) \\ \mathcal{U}_6(X, Y, T) & \mathcal{U}_7(X, Y, T) \end{bmatrix}^\top, \quad (35)
$$

and a least-squares estimate of parameters is given by

$$
\theta_{LS} = \left( \sum_{k=0}^{N_t} \sum_{j=0}^{N_y} \sum_{i=0}^{N_x} \varphi(i, j, k) \varphi^{\top}(i, j, k) \right)^{-1} \\
\cdot \sum_{k=0}^{N_t} \sum_{j=0}^{N_y} \sum_{i=0}^{N_x} \varphi(i, j, k) \mathcal{U}_1(i, j, k). \qquad (36)
$$

In that way, the vector $\theta_{LS}$ gives the initial parameters
of a PDE model with $n_D = 0$, $n_{w\cdot} = 0$ and $n_{\beta\cdot} = 0$.
To complete the initial vector (20), small values for the
missing parameters corresponding to $n_D > 0$, $n_{w\cdot} > 0$
and $n_{\beta\cdot} > 0$ are considered.

5. Simulation

The goal of this simulated test is to validate the pro-
posed tools for applications in ecology, more specifi-
cally for the study of bird population evolution based on
data collected with the monitoring program presented in
Jiguet et al. (2012). Section 6 presents such an applica-
tion with real data.

5.1. Simulated system

Consider a population described by a parameter-
varying formulation of the Fisher PDE (Fisher, 1937)

given by

$$
\frac{\partial u(x, y, t)}{\partial t} = D(H) \left( \frac{\partial^2 u(x, y, t)}{\partial x^2} + \frac{\partial^2 u(x, y, t)}{\partial y^2} \right) + \beta_1(H)u(x, y, t) - \beta_2(H)u^2(x, y, t), \quad (37)
$$

with

$$
\begin{equation}
\begin{aligned}
D(H) &= D_0 + D_1 H(x, y), \\
\beta_{\bullet}(H) &= \beta_{\bullet 0} + \beta_{\bullet 1} H(x, y).
\end{aligned}
\tag{38}
\end{equation}
$$

For the identification procedure, the developments are
those presented in the previous section without parame-
ters w• and with n<sub>D</sub> = 1 and n<sub>β•</sub> = 1.

5.2. *Simulation protocol*

An ecological application often presents some speci-
ficities: a small number of time samples, a spatial do-
main given by a map and a kriging preprocessing of
data. It is proposed to test and to validate the tools
in these specific conditions with 21 time samples, the
France map and noise-free data since the kriging pre-
processing smooths the noise.

The data-set is provided by a simulation of the PDE
(37) using a finite element method. A triangular mesh
with a maximum mesh edge length of 20 is generated
for the spatial domain. The geometry of the spatial do-
main is presented in Figure 1. In order to limit the effect
of the boundary conditions, the simulation domain is a
wide square (blue square in Figure 1) and, to build the
data-set, only the data that belongs to the France map
are considered. In other words, all grid points outside
the France map are simulated with the finite element
method, but they are ignored in the final data-set. The
rectangles specify the value of the scheduling variable:
H(x, y) = 0 for the coordinates inside the bottom rect-
angle, and H(x, y) = 1 inside the top one.

Let us select the true parameters D₀ = 2000, D₁ =
-200, β₁₀ = -0.04, β₁₁ = 0.006, β₂₀ = -0.01 and
β₂₁ = -0.02, and the sampling distances dt = 1, dx = 2
and dy = 2, respectively. The units are omitted in this
simulation application. Neumann boundary conditions
are specified and initial conditions are defined as fol-
lows

$$
u(x, y, 0) = \sum_{n=1}^{10} e^{-(x-x_n)^2+(y-y_n)^2}/2\sigma_n^2, \quad (39)
$$

i.e. a sum of Gaussian centered on (x_n, y_n) ∈ [100, 900]
with variance σ_n^2 ∈ [3000, 6000].

By analogy with the ecological application of Sec-
tion 6, the scheduling variable *H(x, y)* has been chosen
to represent the agricultural intensification: the bottom
<figure><img src="image_3.png" /><figcaption>Figure 1: Simulation domain (blue square), data domain (France map) and scheduling variable (blue rectangles).</figcaption></figure>

rectangle is an area unfavourable to biodiversity and the top rectangle is favourable. Therefore, the diffusion linked to the parameter *D* is larger in the bottom rectangle than the top one, the decreasing rate ($\beta_1$) is faster and the carrying capacity ($\beta_2/\beta_1$) is smaller. In other words, the decreasing trend of the bird population is larger in the bottom rectangle than in the top rectangle.

The data-set size is $N_x = 500$, $N_y = 500$ and $N_t = 20$, respectively. Figure 2 presents the population $u(x, y, t)$ for some time samples.

<figure><img src="image_4.png" /><figcaption>Figure 2: Data $u(x, y, t)$ for $t = 0, 1, 2, 5, 10$ and 20.</figcaption></figure>

Notice that, hereafter, the estimation data-set presented in Figure 2 will be used for the estimation step. For the estimated model validation, another data-set is generated in the same way but with different values of ($x_n, y_n$) and $\sigma_n^2$ for the initial condition (39).

### 5.3. Selection of the design parameters

The iterative identification procedure introduced in Section 4 requires to choose two design parameters: *M*, the number of snapshots in Eq. (6) and *N*, the number of basis functions $\phi_n(x, y)$ in Eq. (13). For *M*, in the proposed context of a small number of time samples, it is preferable to consider all the snapshots, i.e. $M = 21$ for the above simulation protocol. The design parameter *N* must be large enough to consider the largest eigenvalues $\lambda_n$ and to represent correctly the system behaviour, and a limited value gives a parsimonious model. The POD yields the normalized eigenvalues $\lambda_n$ shown in Figure 3. By considering increasing values of *N* and by comparing Levenberg-Marquardt algorithm results in terms of fitting defined by

$$ FIT = 100 \left( 1 - \frac{\|u - \hat{u}\|_2}{\|u - \text{mean}(u)\|_2} \right), \quad (40) $$

table 1 presents the fittings obtained for the estimation and validation data-sets. A number *N* = 7 gives a sufficiently good fitting and a parsimonious model. In other words, the eigenvalues from $\lambda_1 = 0.866$ to $\lambda_7 = 6.6e-8$ are sufficient to represent the system behaviour. The first seven orthonormal basis functions $\phi_n(x, y)$ are presented in Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 3: Normalized eigenvalue $\lambda_n$ in a decreasing order.</figcaption></figure>

### 5.4. Comparison with a constant-parameter PDE model

Here, the goal is to evaluate the benefit of the parameter-varying PDE model compared to a constant-parameter PDE model. By considering data-sets presented in Section 5.2, two PDE models defined by (37) are estimated: one with parameters varying with $n_D = 1$,
<table><caption>Table 1: FIT (%) vs N</caption><thead><tr><th>N</th><th>Results with estimation data-set</th><th>Results with validation data-set</th></tr></thead><tbody><tr><td>3</td><td>61.80</td><td>54.83</td></tr><tr><td>4</td><td>83.31</td><td>80.98</td></tr><tr><td>5</td><td>94.13</td><td>93.41</td></tr><tr><td>6</td><td>99.31</td><td>99.24</td></tr><tr><td>7</td><td>99.44</td><td>99.43</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: The seven considered basis functions $\phi_n(x, y)$ with the corresponding eigenvalue $\lambda_n$.</figcaption></figure>

$n_{\beta_1} = 1, n_{\beta_2} = 1$, and one with constant parameters, i.e. with $n_D = 0, n_{\beta_1} = 0, n_{\beta_2} = 0$, respectively. The design parameters are $M = 21$ and $N = 7$.

Table 2 shows the obtained fittings and Figure 5 presents the spatial error at some time samples. It is clear that the parameter-varying PDE model gives a better representation of the system than the constant-parameter PDE model. The spatial error of the constant-parameter PDE model in Figure 5 presents large amplitudes between top and bottom rectangles.

<table><caption>Table 2: FIT (%) for PDE models with parameters varying and constant parameters</caption><thead><tr><th>PDE model with</th><th>Estimation data</th><th>Validation data</th></tr></thead><tbody><tr><td>parameters varying</td><td>99.44</td><td>99.43</td></tr><tr><td>constant parameters</td><td>93.53</td><td>92.79</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 5: Spatial error (%) defined by $(u(x, y, t) - \hat{u}(x, y, t)) / \text{mean}_{x,y}(u(x, y, t))$ for any x, any y and $t = 1, 10, 20$ with the validation data-set.</figcaption></figure>

### 5.5. Numerical implementation of integrals

Here, the goal is to evaluate the impact of the numerical implementation of integrals in the convergence of the Levenberg-Marquardt algorithm, and more specifically in the context of a small number of time samples. In other words, the numerical implementation of integrals in (16-17) and (25-26) causes some approximations. How do these approximations affect the convergence?

In order to observe the effect of the numerical implementation of integrals, a new data-set is used. Consider the so-called POD-Galerkin data-set where $u$ is generated with the POD-Galerkin solution (13) with the true
parameters and the basis functions set {$\phi_n(x, y)$}$_{n=1}$ obtained from the data-set simulated with the finite element method of Section 5.2. Thus, the POD-Galerkin data-set must lead to a zero error on the estimated parameters, i.e. the Levenberg-Marquardt algorithm converges towards the exact true parameters.

All integrals are implemented with the trapezoidal rule. The design parameter *N* is 7. Two different sampling times are used in the simulation of the ODE system (16) and two different snapshot numbers *M* in Eq. (6) are considered: *dt* = 1 and *M* = 21, for the first case, and *dt* = 0.5 and *M* = 41, for the second one.

The relative errors between true and estimated parameters are given in Table 3. It can be concluded that the estimation error tends towards zero for a decreasing sampling time. For applications with a small number of time samples such as *dt* = 1 and *M* = 21, the approximation due to the numerical implementation of integrals is acceptable (less than 1 %, even if it is not negligible).

<table><caption>Table 3: Relative errors (%) after convergence of the Levenberg-Marquardt algorithm with the POD-Galerkin data-set</caption><thead><tr><th></th><th>dt = 1 and M = 21</th><th>dt = 0.5 and M = 41</th></tr></thead><tbody><tr><td>D<sub>0</sub></td><td>0.04</td><td>0.001</td></tr><tr><td>D<sub>1</sub></td><td>0.68</td><td>0.018</td></tr><tr><td>β<sub>10</sub></td><td>0.02</td><td>0.001</td></tr><tr><td>β<sub>11</sub></td><td>0.26</td><td>0.015</td></tr><tr><td>β<sub>20</sub></td><td>0.74</td><td>0.029</td></tr><tr><td>β<sub>21</sub></td><td>0.72</td><td>0.032</td></tr></tbody></table>

## 5.6. Fitting in terms of trend

As mentioned in the introduction, in ecological applications, the commonly used tools are based on ODE models which describe the global trend of the population in time domain (Wu, 1996; Okubo and Levin, 2001; Mouysset, 2012; Mouysset et al., 2012, 2016). Even if our PDE models represent the population in both spatial domain and time domain, the fitting in terms of global trend is now considered. The normalized temporal trend is defined by

$$f_t(t) = \sum_{j=0}^{N_v} \sum_{i=0}^{N_x} f(i(dx, j.dy, t)) / \sum_{j=0}^{N_v} \sum_{i=0}^{N_x} f(i(dx, j.dy, 0)). \quad (41)$$

The considered data-set is the one presented in Section 5.2 and the design parameters are *M* = 21 and *N* = 7, respectively. Figure 6 presents the true temporal trend of the population with the one obtained by simulating the parameter-varying PDE model estimated with the

Levenberg-Marquardt algorithm. Notice that the estimated PDE model represents almost perfectly the global trend. Once again, in order to evaluate the benefit of the parameters varying, Figure 7 presents the same curves but with the constant-parameter PDE model. In this case, the error is more significant which justifies the use of a parameter-varying formulation.

<figure><img src="image_8.png" /><figcaption>Figure 6: Trends calculated with the true data (blue circles) and the estimated parameter-varying PDE model simulation (green).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Trends calculated with the true data (blue circles) and the estimated constant-parameter PDE model simulation (red).</figcaption></figure>

## 5.7. Validation with a richer scheduling variable

Here, the goal is to test the proposed tool with a more complicated scheduling variable than the one used in the previous sub-sections. Consider the simulation protocol introduced in Section 5.2 with the simulation domain and the scheduling variable presented in Figure
8. Here, the four inner squares specify the value of the scheduling variable: $H(x, y) = 1$ for the coordinates inside the right top square, $H(x, y) = 0.5$ inside the right bottom one, $H(x, y) = 0.25$ inside the left bottom one and $H(x, y) = 0$ inside the left top one. With this new scheduling variable, the selected true parameters are $D_0 = 2000$, $D_1 = -250$, $\beta_{10} = -0.04$, $\beta_{11} = 0.007$, $\beta_{20} = -0.015$ and $\beta_{21} = -0.03$.

<figure><img src="image_10.png" /><figcaption>Figure 8: Simulation domain (outer blue square), data domain (France map) and scheduling variable (inner blue squares) for the global implementation validation.</figcaption></figure>

In order to test the convergence, the Levenberg-Marquardt algorithm is initialized with a parameter vector $\theta_0$ at 80% of the true values. The design parameters are $M = 21$ and $N = 6$, respectively. The convergence is shown in Figure 9.

After convergence, the estimated parameters are close to true values and the obtained fitting is 99.94 %. Table 4 shows the obtained relative errors on the parameters.

It is important to note that the estimation error is due to both the numerical implementation of integrals and the residual error between the simulations of $u$ and $\hat{u}$. This residual error cannot be null because, in the present simulation example, $u$ is given by a finite element method which is a more accurate way to simulate the PDE (37) than the Galerkin method used to calculate $\hat{u}$. This explains why the relative error is larger than the values shown in Table 3. Nevertheless, this error, with a maximum of 11.34 %, is acceptable taking into account the difficult context of the PDE model identification with only 21 time samples.

<figure><img src="image_11.png" /><figcaption>Figure 9: Convergence of the Levenberg-Marquardt algorithm: true parameters (red line) and estimated parameters with the finite-element data-set (blue triangles).</figcaption></figure>

# 6. European Stonechat *Saxicola torquatus* population

## 6.1. Bird data

The considered bird data is given by the national French Breeding Bird Survey (Jiguet et al., 2012). This is a volunteer-based standardized monitoring based on birdwatchers which count spring birds annually during 5 minutes exactly on 10 fix points within a randomly selected square. Figure 10 shows the spatial distribution of the monitored squares. Therefore, this database provides the information related to the bird abundances across the whole country. Abundance values for each species are available for the period 2002-2014. In this paper, only one species is considered, the European Stonechat *Saxicola torquatus* drawn at the bottom right of Figure 11. A spatial interpolation of these abundance data is performed to obtain relative abundance values for each possible 2 × 2 km square in the country
<table><caption>Table 4: Relative errors (%) after convergence of the Levenberg-Marquardt algorithm</caption><thead><tr><th></th><th>Relative error</th></tr></thead><tbody><tr><td>D<sub>0</sub></td><td>0.28</td></tr><tr><td>D<sub>1</sub></td><td>0.01</td></tr><tr><td>β<sub>10</sub></td><td>0.89</td></tr><tr><td>β<sub>11</sub></td><td>6.54</td></tr><tr><td>β<sub>20</sub></td><td>11.34</td></tr><tr><td>β<sub>21</sub></td><td>5.73</td></tr></tbody></table>

(i.e. about 136 000 squares) using kriging models based on spatial autocorrelation and the exponential function (Doxa et al., 2010; Mouysset et al., 2012). The considered kriged data-set for 2002-2014 is presented in Figure 11.

<figure><img src="image_12.png" /><figcaption>Figure 10: Spatial distribution of the squares monitored in continental France by the Breeding Bird Survey scheme.</figcaption></figure>

## 6.2. Land use index

In order to introduce the landscape pattern given by the CORINE Land Cover<fnref n="4" /> in the parameter-varying PDE model, the scheduling variable $H(x, y)$ stands for a land use diversity index inspired by the ecological Shannon index (Hill, 1973) defined on a disc $C$ of centre $(x, y)$ by

$$H(x, y) = - \sum_{i=1}^{R} p_i \ln p_i, \quad (42)$$

with $p_i$, the proportion of the $i$-th land use and $R$, the number of land uses on $C$ (radius 1.5 km) deduced from

<figure><img src="image_13.png" /><figcaption>Figure 11: European Stonechat *Saxicola torquatus* ©Katia LIPOVOI-LPO (bottom right) and kriged abundance of this specie in France for 2002-2014.</figcaption></figure>

the CORINE Land Cover. Figure 12 shows two examples of computation. Thereby, the scheduling variable is small if the landscape diversity is small (Figure 12 right), for instance in area with an intensive agriculture or a wide urban density. On the contrary, the scheduling variable is large for region with varied biotopes (Figure 12 left). The considered land use index has been calculated from CORINE Land Cover 2006 and is plotted in Figure 13.

<footnote n="4">⁴https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-corine</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 12: Examples of computation of the land use index. CORINE Land Cover frontiers (black), disc C (blue) and centre (red cross) with H(456, 400) = 2.71 (left) and H(458, 676) = 0.27 (right).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Land use index H from CORINE Land Cover 2006.</figcaption></figure>

### 6.3. Considered model

The parameter-varying PDE model has been chosen with the following structure

$$ \frac{\partial u(x, y, t)}{\partial t} = D(H) \left( \frac{\partial^2 u(x, y, t)}{\partial x^2} + \frac{\partial^2 u(x, y, t)}{\partial y^2} \right) - w_x(H) \frac{\partial u(x, y, t)}{\partial x} - w_y(H) \frac{\partial u(x, y, t)}{\partial y} + \beta_1(H)u(x, y, t) - \beta_2(H)u^2(x, y, t), \quad (43) $$

where

$$ \begin{aligned} D(H) &= D_0 + D_1 H(x, y), \\ w_{\bullet}(H) &= w_{\bullet 0} + w_{\bullet 1} H(x, y), \\ \beta_{\bullet}(H) &= \beta_{\bullet 0} + \beta_{\bullet 1} H(x, y). \end{aligned} \quad (44) $$

Therefore, the parameter vector to be estimated is

$$ \boldsymbol{\theta} = [D_0 \ D_1 \ w_{x_0} \ w_{x_1} \ w_{y_0} \ w_{y_1} \ \beta_{10} \ \beta_{11} \ \beta_{20} \ \beta_{21}]^\mathrm{T}. \quad (45) $$

### 6.4. Identification results

From the 3D partial moments and the least-squares estimate (36), the following estimated parameters are

obtained: $D_0 = 113.638$, $w_{x_0} = 4.718$, $w_{y_0} = 2.823$, $\beta_{10} = -0.124$ and $\beta_{20} = -0.036$. Then small values for the missing parameters are selected in the initialization vector as follows

$$ \boldsymbol{\theta}_0 = [113.638 \ -0.1 \ 4.718 \ -0.1 \ 2.823 \\ -0.1 \ -0.124 \ 0.01 \ -0.036 \ 0.001]^\mathrm{T}. \quad (46) $$

From this initialization, the Levenberg-Marquardt algorithm yields the estimated parameter vector

$$ \boldsymbol{\theta}_{22} = [37.617 \ 67.858 \ -35.243 \ 54.480 \ 19.641 \\ -14.391 \ 0.697 \ -0.824 \ 0.228 \ -0.255]^\mathrm{T}. \quad (47) $$

The simulation of this estimated model is shown in Figure 14.

<figure><img src="image_16.png" /><figcaption>Figure 14: Estimated abundance given by the parameter-varying PDE model (43) with parameters (47).</figcaption></figure>
### 6.5. Fitting in terms of trend

Once again, even if our PDE model represents the population in both spatial domain and time domain, the fitting in terms of global trend is considered by using the normalized temporal trend defined by (41). More specifically, the goal is to compare the estimated PDE model with the ODE logistic model (Mouysset et al. (2016)) defined by

$$u_t((k+1) dt) = u_t(k dt)+u_t(k dt)R_s\left(1-\frac{u_t(k dt)}{K}\right), (48)$$

where $u_t(t)$ is the data-set with $t \in [2002, 2014]$ calculated from (41) and plotted (blue circles) in Figure 15, $R_s$ and $K$ are parameters estimated by least-squares from $u_t(t)$. As shown in Figure 15, the population trend given by the estimated PDE model (green dashed curve) approximates the global population trend like the ODE logistic model (magenta dashed curve).

<figure><img src="image_17.png" /><figcaption>Figure 15: Measured population trend $u_t(t)$ (blue circles), ODE logistic model simulation (magenta dashed line) and estimated PDE model trend simulated with (47) and calculated from (41) (green dashed line).</figcaption></figure>

## 7. Conclusion

In this paper, a specific attention has been paid to the estimation of parameters of a parameter-varying partial differential equation model for the bird population dynamic. In order to reach this goal, a new identification method based on the Galerkin method and the proper orthogonal decomposition has been presented. This method is implemented with a Levenberg-Marquardt algorithm initialized by a least-squares estimate based on 3D partial moments.

The proposed tools have been validated in simulation in the following conditions:

* a parametric estimation of nonlinear PDE models with parameters varying,

* considered systems without exogenous input and with a unique excitation provided by initial conditions,

* and a small number of time samples.

Finally, the algorithms are applied to European Stonechat population with data collected in the national French Breeding Bird Survey. The scheduling variable of the parameter-varying formulation is the landscape pattern given by the CORINE Land Cover data.

The future works consists in applying tools to a community of birds specialized of agricultural area, considering a Shannon index in the parameter-varying formulation which characterize the agricultural practice and developing tools to help take biodiversity goals into public policies into account for Common Agricultural Policy from European Union, for instance.

## References

* Belforte, G., Dabbene, F., Gay, P., 2005. LPV approximation of distributed parameter systems in environmental modelling. Environmental Modelling and Software 20, 1063–1070.

* Björck, A., 1996. Numerical methods for least squares problems. SIAM.

* Chamberlain, D., Fuller, R., Bunce, R., Duckworth, J., Shrubb, M., 2000. Changes in the abundance of farmland birds in relation to the timing of agricultural intensification in England and Wales. Journal of Applied Ecology 37, 771–788.

* Debnath, L., Mikusinski, P., 2005. Introduction to Hilbert spaces with applications. Academic Press.

* Doxa, A., Bas, Y., Paracchini, M. L., Pointereau, P., Terres, J. M., Jiguet, F., 2010. Low-intensity agriculture increases farmland bird abundances in France. Journal of Applied Ecology 47, 1348–1356.

* Farah, M., Mercère, G., Ouvrard, R., Poinot, T., 2016. Combining least-squares and gradient-based algorithms for the identification of a co-current flow heat exchanger. Int. J. Control, 1–13.

* Fisher, R. A., 1937. The wave of advance of advantageous genes. Annals of Eugenics 7(4), 355–369.

* Hallmann, C. A., Sorg, M., Jongejans, E., Siepel, H., Hofland, N., Schwan, H., Stenmans, W., Müller, A., Sumser, H., Hörren, T., Goulson, D., de Kroon, H., 2017. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12(10).

* Hill, M. O., 1973. Diversity and evenness: a unifying notation and its consequences. Ecology 54, 427–432.

* Holmes, E. E., Lewis, M. A., Banks, J. E., Veit, R. R., 1994. Partial differential equations in ecology: spatial interactions and population dynamics. Ecology 75(1), 17–29.

* Inger, R., Gregory, R., Duffy, J. P., Stott, I., Vorisek, P., Gaston, K. J., 2014. Common European birds are declining rapidly while less abundant species? numbers are rising. Ecology Letters.

* IPBES (Ed.), March 2018a. IPBES-6 Plenary. Medellin, Colombia.

* IPBES, 2018b. Summary for policymakers of the assessment report on land degradation and restoration of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. R. Sc-holes, L. Montanarella, A. Brainich, N. Barger, B. ten Brink,
M. Cantele, B. Erasmus, J. Fisher, T. Gardner, T. G. Holland, F. Kohler, J. S. Kotiaho, G. Von Maltitz, G. Nangendo, R. Pandit, J. Parrotta, M. D. Potts, S. Prince, M. Sankaran and L. Willemen (Eds.). IPBES secretariat, Bonn, Germany.

IPBES, 2018c. Summary for policymakers of the regional assessment report on biodiversity and ecosystem services for Europe and Central Asia of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. M. Fischer, M. Rounsevell, A. Torre-Marin Rando, A. Mader, A. Church, M. Elbakidze, V. Elias, T. Hahn, P.A. Harrison, J. Hauck, B. Martín-López, I. Ring, C. Sandström, I. Sousa Pinto, P. Visconti, N.E. Zimmermann and M. Christie (Eds.), IPBES secretariat, Bonn, Germany.

IPBES, 2018d. Summary for policymakers of the regional assessment report on biodiversity and ecosystem services for the Americas of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. J. Rice, C.S. Seixas, M.E. Zaccagnini, M. Bedoya-Gaitán, N. Valderrama, C.B. Anderson, M.T.K. Arroyo, M. Bustamante, J. Cavender-Bares, A. Diaz-de-Leon, S. Fennessy, J. R. García Márquez, K. Garcia, E.H. Helmer, B. Herrera, B. Klatt, J.P. Ometo, V. Rodríguez Osuna, F.R. Scarano, S. Schill and J. S. Farinaci (Eds.), IPBES secretariat, Bonn, Germany.

Jiguet, F., Devictor, V., Julliard, R., Couvet, D., October 2012. French citizens monitoring ordinary birds provide tools for conservation and ecological sciences. Acta Oecologica 44, 58–66.

Knudsen, M., 1994. A sensitivity approach for estimation of physical parameters. In: 10th IFAC Symposium on System Identification. Vol. 2. pp. 231–237.

Levin, S., 1976. Population dynamic models in heterogeneous environments. Annual review of ecology and systematics 7, 287–310.

Li, H., Qi, C., 2010. Modeling of distributed parameter systems for applications - A synthesized review from time-space separation. Journal of Process Control 20 (2010) 20, 891–901.

Liang, Y., Lee, H., Lim, S., Lin, W., Lee, K., Wu, C., 2002. Proper orthogonal decomposition and its applications. Part I: theory. Journal of Sound and Vibration 252, 527–544.

Ljung, L., 1999. System identification - Theory for the User. Prentice-Hall, Upper Saddle River, N.J., 2nd edition.

Mouysset, L., 2012. Les politiques publiques au défi de la biodiversité : modèles et scénarios bio-économiques pour une agriculture durable. Ph.D. thesis, Muséum National d'Histoire Naturelle.

Mouysset, L., Doyen, L., Jiguet, F., 2012. Different policy scenarios to promote various targets of biodiversity. Ecological Indicators 14, 209–221.

Mouysset, L., Miglianico, M., Makowski, D., Jiguet, F., Doyen, L., 2016. Selection of dynamic models for bird populations in farms. Environmental Modeling & Assessment 21 (3), 407–418.

Newman, A., 1996a. Model reduction via the Karhunen-Loeve expansion, Part I: An exposition. Technical report, T.R. 96-32, University of Maryland, ISR. MD, USA.

Newman, A., 1996b. Model reduction via the Karhunen-Loeve expansion, Part II: Some elementary examples. Technical report, T.R. 96-32, University of Maryland, ISR. MD, USA.

Nocedal, J., Wright, S., 2006. Numerical optimization, 2nd Edition. Springer.

Okubo, A., 1986. Dynamical aspects of animal grouping: Swarms, schools, flocks, and herds. Advances in Biophysics 22, 1–94.

Okubo, A., Levin, S., 2001. Diffusion and ecological problems: modern perspectives, Second Edition. Springer.

Ouvrard, R., Trigeassou, J. C., 2011. On embedded FIR filter models for identifying continuous-time and discrete-time transfer functions: the RPM approach. Int. J. Control 84 (3), 616–632.

Pham, D., Mercère, G., Ouvrard, R., Poinot, T., 2018. Heat equation parameter estimation based on the POD-Galerkin approach. In: 18th IFAC Symposium on System Identification.

Polis, M., Goodson, R., Wozny, M., 1973. On parameter identification for distributed systems using Galerkin's criterion. Automatica 9, 53–64.

Schorsch, J., Gilson, M., Laurain, V., Garnier, H., 2013. Identification of LPV partial differential equation models. In: 52nd IEEE Conference on Decision and Control. Florence, Italy.

Schorsch, J., Laurain, V., Gilson, M., Garnier, H., 2014. Refined IV-based method for LPV partial differential equation model identification. In: 13th European Control Conference. Strasbourg, France.

Shigesada, N., Kawasaki, K., Teramoto, E., 1979. Spatial segregation of interacting species. Journal of Theoretical Biology 79-1, 83–99.

Siriwardena, G., Baillie, S., Buckland, S., Fewster, R., Marchant, J., Wilson, J., 1998. Trends in the abundance of farmland birds: a quantitative comparison of smoothed Common Birds Census indices. Journal of Applied Ecology 35, 24–43.

Skellam, J. G., 1951. Random dispersal in theoretical populations. Biometrika 38, 196–218.

Stanton, R. L., Morrissey, C. A., Clark, R. G., 2018. Analysis of trends and agricultural drivers of farmland bird declines in North America: A review. Agriculture, Ecosystems and Environment 254, 244–254.

Toth, R., 2010. Modeling and identification of linear parameter-varying systems. Springer.

Trigeassou, J. C., 1987. Contribution à l'extension de la méthode des moments en automatique. application à l'identification des systèmes linéaires. Thèse d'etat, Université de Poitiers, France.

Uciński, D., 2005. Optimal measurement methods for distributed parameter system identification. CRC Press.

Wu, J., 1996. Theory and applications of partial functional differential equations. Vol. 119. Springer, New York.

WWF, 2018. Living planet report - 2018: Aiming higher. Grooten, M. and Almond, R.E.A. (Eds.), WWF, Gland, Switzerland.

## Appendix A. Galerkin method calculations

The Galerkin method (Section 3.3) aims at determining the time coefficients ${a_n(t)}_{n=1}^N$ of (13) from the known function basis ${\phi_n(x, y)}_{n=1}^N$. With the substitution of the approximate solution (13) in (1), in the left-hand-side of the corresponding equation, the time derivative becomes

$$ \begin{aligned} \frac{\partial \hat{u}(x, y, t)}{\partial t} &= \sum_{n=1}^{N} \frac{\partial a_n(t)}{\partial t} \phi_n(x, y), \\ &= \sum_{n=1}^{N} \dot{a}_n(t) \phi_n(x, y), \end{aligned} \quad (A.1) $$
while, on the right-hand-side, the spatial derivatives become

$$ 
\begin{aligned}
\frac{\partial^2 \hat{u}(x, y, t)}{\partial x^2} &= \sum_{n=1}^{N} a_n(t) \frac{\partial^2 \phi_n(x, y)}{\partial x^2} + \frac{\partial^2 \bar{u}(x, y)}{\partial x^2}, \\
\frac{\partial^2 \hat{u}(x, y, t)}{\partial y^2} &= \sum_{n=1}^{N} a_n(t) \frac{\partial^2 \phi_n(x, y)}{\partial y^2} + \frac{\partial^2 \bar{u}(x, y)}{\partial y^2}, \\
\frac{\partial \hat{u}(x, y, t)}{\partial x} &= \sum_{n=1}^{N} a_n(t) \frac{\partial \phi_n(x, y)}{\partial x} + \frac{\partial \bar{u}(x, y)}{\partial x}, \\
\frac{\partial \hat{u}(x, y, t)}{\partial y} &= \sum_{n=1}^{N} a_n(t) \frac{\partial \phi_n(x, y)}{\partial y} + \frac{\partial \bar{u}(x, y)}{\partial y}.
\end{aligned}
\quad (A.2) $$

Lastly, the quadratic term in (1) yields

$$ 
\begin{aligned}
\hat{u}^2(x, y, t) = & \left( \sum_{n=1}^{N} a_n(t) \phi_n(x, y) \right)^2 \\
& + 2 \sum_{n=1}^{N} a_n(t) \phi_n(x, y) \bar{u}(x, y) + \bar{u}^2(x, y),
\end{aligned}
\quad (A.3) $$

and

$$ 
\begin{aligned}
\left( \sum_{n=1}^{N} a_n(t) \phi_n(x, y) \right)^2 = & \sum_{n=1}^{N} a_n^2(t) \phi_n^2(x, y) \\
& + 2 \sum_{n=1}^{N-1} \sum_{m=n+1}^{N} a_n(t) a_m(t) \phi_n(x, y) \phi_m(x, y).
\end{aligned}
\quad (A.4) $$

Therefore, the substitution of (13) in (1) leads to Equation (14).

By applying the inner product $\langle Eq.(14), \phi_l \rangle = \int Eq.(14) \phi_l$, the orthonormal property of $\{\phi_n(x, y)\}_{n=1}^N$ given by

$$ \langle \phi_n, \phi_l \rangle = \int \phi_n \phi_l = \begin{cases} 1 & \text{if } n = l \\ 0 & \text{if } n \neq l \end{cases} \quad (A.5) $$

leads to the following simplification on the left-hand-side

$$ \int \sum_{n=1}^{N} \dot{a}_n(t) \phi_n \phi_l = \dot{a}_l(t) \quad (A.6) $$

and to Equation (15).

Remind that

$$ \beta_1(H) = \sum_{\nu=0}^{n_{\beta_1}} \beta_{1_\nu} H^\nu(x, y). \quad (A.7) $$

Therefore, on the right-hand-side of (15), the following simplification can be done

$$ 
\begin{aligned}
\sum_{n=1}^{N} a_n(t) \int \beta_1(H) \phi_n \phi_l &= \beta_{10} a_l(t) \\
&+ \sum_{n=1}^{N} a_n(t) \int \sum_{\nu=1}^{n_{\beta_1}} \beta_{1_\nu} H^\nu(x, y) \phi_n \phi_l.
\end{aligned}
\quad (A.8) $$

Hence, (15) can be rewritten as follows

$$ 
\begin{aligned}
\dot{a}_l(t) = & \sum_{n=1}^{N} a_n(t) \int \left\{ D_x(H) \frac{\partial^2 \phi_n}{\partial x^2} + D_y(H) \frac{\partial^2 \phi_n}{\partial y^2} - w_x(H) \frac{\partial \phi_n}{\partial x} - w_y(H) \frac{\partial \phi_n}{\partial y} + \sum_{\nu=1}^{n_{\beta_1}} \beta_{1_\nu} H^\nu(x, y) \phi_n \right\} \phi_l \\
& - 2 \sum_{n=1}^{N} a_n(t) \int \beta_2(H) \phi_n \phi_l \bar{u} \\
& - \sum_{n=1}^{N} a_n^2(t) \int \beta_2(H) \phi_n^2 \phi_l \\
& - 2 \sum_{n=1}^{N-1} \sum_{m=n+1}^{N} a_n(t) a_m(t) \int \beta_2(H) \phi_n \phi_m \phi_l \\
& + \int \left\{ D_x(H) \frac{\partial^2 \bar{u}}{\partial x^2} + D_y(H) \frac{\partial^2 \bar{u}}{\partial y^2} - w_x(H) \frac{\partial \bar{u}}{\partial x} - w_y(H) \frac{\partial \bar{u}}{\partial y} + \beta_1(H) \bar{u} - \beta_2(H) \bar{u}^2 \right\} \phi_l
\end{aligned}
\quad (A.9) $$

for $l = 1, \dots, N$.

By considering the vector $a(t)$, the matrix formulation (16) is obtained where

* $\Gamma_1, \Gamma_2$ and $b_1$ are directly deduced from (A.9),

* the simplification (A.8) leads to the term $\beta_{10} \mathbf{I}_N \mathbf{a}(t)$ on the right-hand-side of (16) and to the definition of $\mathcal{F}_1(\bullet)$ which presents the particular sum $\sum_{\nu=1}^{n_{\beta_1}} \beta_{1_\nu} H^\nu(x, y) \bullet$ in (17),

* all terms $\int \bullet \phi_n \phi_m \phi_l$ due to the square non-linearity in the PDE definition are included in the term $(\mathbf{I}_N \otimes \mathbf{a}^\top(t)) \Lambda \mathbf{a}(t)$.