# Biology-Informed inverse problems for insect pests detection using pheromone sensors

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, Paul-André Calatayud, Philippe Lucas, Simon Labarthe

► To cite this version:

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, et al.. Biology-Informed inverse problems for insect pests detection using pheromone sensors. 2024. [hal-04572831](https://hal.inrae.fr/hal-04572831)

HAL Id: hal-04572831

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

Preprint submitted on 12 May 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Biology-Informed inverse problems for insect pests detection using pheromone sensors

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, Paul-André Calatayud, Philippe Lucas, Simon Labarthe

► To cite this version:

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, et al.. Biology-Informed inverse problems for insect pests detection using pheromone sensors. 2025. [hal-04572831v3](https://hal.inrae.fr/hal-04572831v3)

HAL Id: hal-04572831

https://hal.inrae.fr/hal-04572831v3

Preprint submitted on 28 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Biology-Informed inverse problems for insect
pests detection using pheromone sensors

Thibault Malou¹,§, Nicolas Parisey², Katarzyna Adamczyk-Chauvat¹, Elisabeta Vergu ¹,†, Beatrice Laroche¹,³, Paul-Andre Calatayud⁴,⁵, Philippe Lucas⁶, and Simon Labarthe¹,⁷,⁸,¶

¹Université Paris-Saclay, INRAE, MaIAGE, 78350,
Jouy-en-Josas, France

²INRAE, Institute of Genetics, Environment and Plant
Protection (IGEPP—Joint Research Unit 1349), Le Rheu,
France

³Université Paris-Saclay, Inria, 91120, Palaiseau, France

⁴Institut de Diversité, Ecologie et Evolution du Vivant
(IDEEV), Université Paris-Saclay, CNRS, IRD, UMR
Evolution, Génome, Comportement et Ecologie (EGCE),
Gif-sur-Yvette, France

⁵International Center of Insect Physiology and Ecology (icipe),
Nairobi, Kenya

⁶Institute of Ecology and Environmental Sciences, INRAE,
Versailles, France,

7 Univ. Bordeaux, INRAE, BIOGECO, 33610, Cestas, France

⁸Univ. Bordeaux, Inria, INRAE, 33400, Talence, France,

<footnote>§First author: thibault.malou@inrae.fr</footnote>

<footnote>†deceased during the development of this study.</footnote>

<footnote>¶corresponding author: simon.labarthe@inrae.fr</footnote>
Biology-Informed inverse problems
for insect pests detection using
pheromone sensors

Thibault Malou¹,§, Nicolas Parisey², Katarzyna Adamczyk-Chauvat¹,
Elisabeta Vergu ¹,†, Beatrice Laroche¹,³, Paul-Andre Calatayud⁴,⁵,
Philippe Lucas⁶, and Simon Labarthe¹,⁷,⁸,¶

¹Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas,
France

²INRAE, Institute of Genetics, Environment and Plant Protection
(IGEPP—Joint Research Unit 1349), Le Rheu, France

³Université Paris-Saclay, Inria, 91120, Palaiseau, France

⁴Institut de Diversité, Ecologie et Evolution du Vivant (IDEEV),
Université Paris-Saclay, CNRS, IRD, UMR Evolution, Génome,
Comportement et Ecologie (EGCE), Gif-sur-Yvette, France

⁵International Center of Insect Physiology and Ecology (icipe),
Nairobi, Kenya

⁶Institute of Ecology and Environmental Sciences, INRAE, Versailles,
France,

⁷ Univ. Bordeaux, INRAE, BIOGECO, 33610, Cestas, France

⁸Univ. Bordeaux, Inria, INRAE, 33400, Talence, France,

<footnote>§First author: thibault.malou@inrae.fr</footnote>

<footnote>†deceased during the development of this study.</footnote>

<footnote>¶corresponding author: simon.labarthe@inrae.fr</footnote>
Abstract

Most insects have the ability to modify the odor landscape in order to communicate with their conspecies during key phases of their life cycle such as reproduction. They release pheromones in their nearby environment, volatile compounds that are detected by insects of the same species with exceptional specificity and sensitivity. Efficient pheromone detection is then an interesting lever for insect pest management in a precision agroecological culture context. A precise and early detection of pests using pheromone sensors offers a strategy for pest management before infestation.

In this paper, we develop a biology-informed inverse problem framework that leverages temporal signals from a pheromone sensor network to build insect presence maps. Prior biological knowledge is introduced in the inverse problem by the mean of a specific penalty, using population dynamics PDE residuals. We benchmark the biology-informed penalty with other regularization terms such as Tikhonov, LASSO or composite penalties in a simplified toy model. We use classical comparison criteria, such as target reconstruction error, or Jaccard distance on pest presence-absence. But we also use more task-specific criteria such as the number of informative sensors during inference. Finally, the inverse problem is solved in a realistic context of pest infestation in an agricultural landscape by the fall armyworm (Spodoptera frugiperda).

**keywords:** Inverse problem, data assimilation, biology-informed estimation, PDE, pest detection, pheromone.
Abstract

Most insects have the ability to modify the odor landscape in order to
communicate with their conspecies during key phases of their life cycle such as
reproduction. They release pheromones in their nearby environment, volatile
compounds that are detected by insects of the same species with exceptional
specificity and sensitivity. Efficient pheromone detection is then an interesting
lever for insect pest management in a precision agroecological culture context.
A precise and early detection of pests using pheromone sensors offers a strategy
for pest management before infestation.

In this paper, we develop a biology-informed inverse problem framework
that leverages temporal signals from a pheromone sensor network to build in-
sect presence maps. Prior biological knowledge is introduced in the inverse
problem by the mean of a specific penalty, using population dynamics PDE
residuals. We benchmark the biological-informed penalty with other regular-
ization terms such as Tikhonov, LASSO or composite penalties in a simplified
toy model. We use classical comparison criteria, such as target reconstruction
error, or Jaccard distance on pest presence-absence. But we also use more task-
specific criteria such as the number of informative sensors during inference.
Finally, the inverse problem is solved in a realistic context of pest infestation
in an agricultural landscape by the fall armyworm (Spodoptera frugiperda).

**keywords:** Inverse problem, data assimilation, biology-informed estimation, PDE, pest detection, pheromone.
# 1. Introduction

Insect pests represent a major threat for agricultural systems, causing direct or indirect damages to crops (Oerke, 2006), with an even increasing burden in a warming climate (Deutsch et al., 2018). World trade, travel and climate changes increase the arrival of alien species that challenge existing pest control management for native insects (EFSA PLH Panel et al., 2018). Intensified insecticide use is not sustainable, due notably to the negative impact of pesticides on biodiversity, environment and human health. Alternative solutions must be developed: a promising perspective is the development of precision agriculture which aims in particular to leverage sensors and information systems to reduce the environmental impact of agricultural practices (Gebbers et al., 2010). In that perspective, early detection of pests before infestation settlement can help designing low-pesticide crop management and epidemiological surveillance policies.

5

10

Insect pheromones are a good target for early detection of pests (Conchou et al., 2019). Indeed, insects pheromones are specie-specific, guaranteeing a good specificity of any pheromone-based detection. Pheromone specificity has long been leveraged for insect detection with the mean of pheromone traps: pheromone diffusion modifies the odor landscape, attracting the targeted insects toward traps where they can be counted (McNeil, 1991). However, if pheromone traps allow for presence/absence detections, they do not allow to accurately locate the insects in their natural habitat. Recently, new advances in sensor developments paved the way towards field-based tools to explore *in situ* pheromone plumes. These sensors have different operating principles: there are insect antenna-based biosensors (Martinez et al., 2014; Myrick et al., 2011; Pawson et al., 2020), olfactory receptor-based biosensors (Khadka et al., 2019), OBP-based biosensors (Lu et al., 2014), and e-nose technologies based on conducting polymers, metal oxide semiconductor and other devices (Moitra et al., 2016; Negri et al., 2008; Steffens et al., 2014; Wehrenfennig et al., 2012). The use of pheromone sensors could allow to directly measure the odor landscape and track back pheromone plumes towards emitters, in order to detect pests at an early stage of invasion.

15

20

25

From a mathematical point of view, the problem of locating the emission source of a molecule knowing measures of the emission plume belongs to the class of inverse problems (Isakov, 2017), the direct problem of which is a plume propagation model. Atmospheric dispersion, plume propagation, and plume-related inverse problems have a large literature, in particular for environmental applications, such as pollution detection (Lushi et al., 2010; Stockie, 2011). Inverse problems, and more specifically data assimilation problems, are most often formulated as an optimization problem aiming to minimize the discrepancies between the direct model outputs and in-field observations by fitting the model parameters (Isakov, 2017). In the specific case where insects are the emitters, locating the emission source amounts to searching for continuous density maps of pests, which induces a high dimensional parameter space to be explored during the inverse problem. Variational approaches are well suited for such high-dimensional optimization problem, for fast computation of gradients of the objective function (Bocquet, 2014). Inverse problems are most of the time ill-posed and must be regularized by additional terms in the loss function that constrain the parameter exploration, e.g. with LASSO regularization term (Adam et al., 2016). Other regularization can be searched for, that could take into account prior biological knowledge incorporated into models, like in physic or biology-informed machine learning approaches (Yazdani et al., 2020), through weak or strong constraints (Trémolet, 2007).

30

35

40

In conventional agricultural practices, farmers often use pesticide prophylactically: knowing the potential occurrence of an infestation, they preventively apply pesticide before a pest is observed. If we interpret this strategy in the context of a classification problem, it is as if they were

45
# 1. Introduction

Insect pests represent a major threat for agricultural systems, causing direct or indirect damages to crops (Oerke, 2006), with an even increasing burden in a warming climate (Deutsch et al., 2018). World trade, travel and climate changes increase the arrival of alien species that challenge existing pest control management for native insects (EFSA PLH Panel et al., 2018). Intensified insecticide use is not sustainable, due notably to the negative impact of pesticides on biodiversity, environment and human health. Alternative solutions must be developed: a promising perspective is the development of precision agriculture which aims in particular to leverage sensors and information systems to reduce the environmental impact of agricultural practices (Gebbers et al., 2010). In that perspective, early detection of pests before infestation settlement can help designing low-pesticide crop management and epidemiological surveillance policies.

Insect pheromones are a good target for early detection of pests (Conchou et al., 2019). Indeed, insects pheromones are species-specific, guaranteeing a good specificity of any pheromone-based detection. Pheromone specificity has long been leveraged for insect detection with the mean of pheromone traps: pheromone diffusion modifies the odor landscape, attracting the targeted insects toward traps where they can be counted (McNeil, 1991). However, if pheromone traps allow for presence/absence detections, they do not allow to accurately locate the insects in their natural habitat. Recently, new advances in sensor developments paved the way towards field-based tools to explore *in situ* pheromone plumes. These sensors have different operating principles: there are insect antenna-based biosensors (Martinez et al., 2014; Myrick et al., 2011; Pawson et al., 2020), olfactory receptor-based biosensors (Khadka et al., 2019), OBP-based biosensors (Lu et al., 2014), and e-nose technologies based on conducting polymers, metal oxide semiconductor and other devices (Moitra et al., 2016; Negri et al., 2008; Steffens et al., 2014; Wehrenfennig et al., 2012). The use of pheromone sensors could allow to directly measure the odor landscape and track back pheromone plumes towards emitters, in order to detect pests at an early stage of invasion.

From a mathematical point of view, the problem of locating the emission source of a molecule knowing measures of the emission plume belongs to the class of inverse problems (Isakov, 2017), the direct problem of which is a plume propagation model. Atmospheric dispersion, plume propagation, and plume-related inverse problems have a large literature, in particular for environmental applications, such as pollution detection (Lushi et al., 2010; Stockie, 2011). Inverse problems, and more specifically data assimilation problems, are most often formulated as an optimization problem aiming to minimize the discrepancies between the direct model outputs and in-field observations by fitting the model parameters (Isakov, 2017). In the specific case where insects are the emitters, locating the emission source amounts to searching for continuous density maps of pests, which induces a high dimensional parameter space to be explored during the inverse problem. Variational approaches are well suited for such high-dimensional optimization problem, for fast computation of gradients of the objective function (Bocquet, 2014). Inverse problems are most of the time ill-posed and must be regularized by additional terms in the loss function that constrain the parameter exploration, e.g. with LASSO regularization term (Adam et al., 2016). Other regularization can be searched for, that could take into account prior biological knowledge incorporated into models, like in physic or biology-informed machine learning approaches (Yazdani et al., 2020), through weak or strong constraints (Trémolet, 2007).

In conventional agricultural practices, farmers often use pesticide prophylactically: knowing the potential occurrence of an infestation, and lacking efficient ways to monitor pest populations, they apply pesticides as a preventive measure, even before a pest is observed. Consequently,
classifying all their fields as "infected" whatever their actual state of infestation, thus minimizing their miss rate, i.e. a non detected true presence of pest. In so doing, they also maximize their false discovery rate. In order to minimize pesticide use, we would like to minimize the false discovery rate, so that pesticides are only used where and when a detected pest is actually present in the field. But, according to a classical trade-off in classification problems, reducing the false discovery rate comes with an increase of the miss rate, so that the detection problem could be reformulated as follows: find the minimal false discovery rate that keeps the miss rate under an acceptable risk, from the point of view of crop protection.

The main objective of this paper is to introduce and test a data assimilation framework, namely the biology-informed data assimilation, which seeks to combine a pheromone propagation model and pheromone sensors data, like in classical data assimilation, but also biological information. In particular, we specifically aim to introduce biology-informed regularization terms, on which the biology-informed data assimilation is based, in order to assess how adding biological priors improves the accuracy of the data assimilation. Several comparison criteria will be used. In particular, we will introduce a specific criteria designed to compare the number of active sensors involved during the inference process.

## 2. Material and methods

A sketch of the global methodological framework can be found in Fig. 1. The method contains three blocks that will be detailed hereafter:

1. the direct problem that aims to model the pheromone propagation in the atmosphere over an agricultural landscape, taking into account environmental data such as land use or wind field, see Sec. 2.1;

2. the pheromone sensors that provide data of the observed variable and an observation operator that maps the state variable of the direct model to the observed variable, see Sec. 2.3.1;

3. the corresponding inverse problem that aims to infer pheromone emissions by the insects and enables the construction of pests risk maps, knowing the data and the pheromone propagation model, see Sec. 2.2.

The proposed framework is illustrated with a toy case and a more realistic case, set up in an agricultural landscape, see Sec. 2.3.2. In the different test cases, the datasets are synthetic and constructed from a spatio-temporal sampling of pheromone concentration and a realistic model of the observation process.

### 2.1. 2D reaction-convection-diffusion model for the pheromone propagation

One of the most classical model to simulate propagation of chemical species in the atmosphere is a 3D Chemical-Transport Model (CTM), that has been extensively reviewed (Leelössy et al., 2014; Sportisse, 2007). The CTM are obtained from the reactive Navier-Stokes equations by using the dilution hypothesis. This hypothesis is valid when the chemical species of interest is present in very small amounts. In this case, the compound is assumed to be chemically passive, which results in a decoupling between atmospheric dynamics and compound concentration dynamics. As a consequence, dynamical field such as the wind or the diffusion coefficients can be computed separately, offline, and given as parameters of the CTM. It results that this (linear and offline-coupled) CTM is much less complex than the full (non-linear) reactive Navier-Stokes equations, or than a CTM coupled online with a fine scale Navier-Stokes model (Leelössy et al.,
sensors capable of measuring the spatio-temporal distribution of pest populations could poten-
tially allow the farmers to reduce the frequency, dose and area of treatment by spraying when,
where and to the extent that the pest is detected. If we interpret this prophylactic strategy in
the context of a classification problem, it is as if they were classifying all their fields as "infected"
whatever their actual state of infestation, thus minimizing their miss rate, i.e. a non detected
true presence of pest. In so doing, they also maximize their false discovery rate. In order to mini-
mize pesticide use, we would like to minimize the false discovery rate, so that pesticides are only
used where and when a detected pest is actually present in the field. But, according to a classical
trade-off in classification problems, reducing the false discovery rate comes with an increase of
the miss rate, so that the detection problem could be reformulated as follows: find the minimal
false discovery rate that keeps the miss rate under an acceptable risk, from the point of view of
crop protection.

The main objective of this paper is to introduce and test a data assimilation framework,
namely the biology-informed data assimilation, which seeks to combine a pheromone propa-
gation model and pheromone sensors data, like in classical data assimilation, but also biological
information. In particular, we specifically aim to introduce biology-informed regularization terms,
on which the biology-informed data assimilation is based, in order to assess how adding biolog-
cal priors improves the accuracy of the data assimilation. Several comparison criteria will be used.
In particular, we will introduce a specific criteria designed to compare the number of active sen-
sors involved during the inference process.

## 2. Material and methods

A sketch of the global methodological framework can be found in Fig. 1. The method contains
three blocks that will be detailed hereafter:

1) the direct problem that aims to model the pheromone propagation in the atmosphere over an agricultural landscape, taking into account environmental data such as land use or wind field, see Sec. 2.1;

2) the pheromone sensors that provide data of the observed variable and an observation operator that maps the state variable of the direct model to the observed variable, see Sec. 2.3.1;

3) the corresponding inverse problem that aims to infer pheromone emissions by the insects and enables the construction of pests risk maps, knowing the data and the pheromone propagation model, see Sec. 2.2.

The proposed framework is illustrated with a toy case and a more realistic case, set up in an
agricultural landscape, see Sec. 2.3.2. In the different test cases, the datasets are synthetic and
constructed from a spatio-temporal sampling of pheromone concentration and a realistic model
of the observation process.

## 2.1. 2D reaction-convection-diffusion model for the pheromone propagation

One of the most classical model to simulate propagation of chemical species in the atmo-
sphere is a 3D Chemical-Transport Model (CTM), that has been extensively reviewed (Leelössy
et al., 2014; Sportisse, 2007). The CTM is obtained from the reactive Navier-Stokes equations
by using the dilution hypothesis. This hypothesis is valid when the chemical species of interest
is present in very small amounts. In this case, the compound is assumed to be chemically pas-
sive, which results in a decoupling between atmospheric dynamics and compound concentration
dynamics. As a consequence, dynamical field such as the wind or the diffusion coefficients can
# Biology-informed Data Assimilation

## Classical Data Assimilation

<figure><img src="image_2.png" /><figcaption>**Figure 1** – Graphical abstract: summary of the proposed Biology-Informed Data Assimilation method for the pest insects localization. A classical Data Assimilation (DA) framework is first built, combining a pheromone propagation model (blue rectangle), pheromone sensors (red rectangle) into an inverse problem (green rectangle). Additionally to this classical DA framework, biological information can be supplemented through specific regularization terms of the DA problem (purple rectangle). These regularization terms are the main insights of this study.</figcaption></figure>
# Biology-informed Data Assimilation

## Classical Data Assimilation

<figure><img src="image_2.png" /><figcaption>**Figure 1** – Graphical abstract: summary of the proposed Biology-Informed Data Assimilation method for the pest insects localization. A classical Data Assimilation (DA) framework is first built, combining a pheromone propagation model (blue rectangle), pheromone sensors (red rectangle) into an inverse problem (green rectangle). Additionally to this classical DA framework, biological information can be supplemented through specific regularization terms of the DA problem (purple rectangle). These regularization terms are the main insights of this study.</figcaption></figure>
2014), while being a reliable approximation: it is able to model the main physical phenomena at play during the transportation of pheromones.

In the present context, only the lower layer of the atmosphere, where insects and sensors are located, is of interest. Hence, for sake of computation load and data availability in real applications, a 2D depth-integrated CTM is used to model the total concentration of pheromone over the depth of the layer $c(x, y, t) [\text{g m}^{-2}]$ on a domain $\Omega \in \mathbb{R}^2$ and during a time window $[0; T]$. This depth-integrated CTM results in a 2D reaction-convection-diffusion model, that can be derived from the classical 3D CTM, see Sec. B.1 in the annex for details, and reads:

$$ (1) \qquad \frac{\partial c}{\partial t} - \nabla \cdot (\mathbf{K} \nabla c) + \nabla \cdot (\vec{u} c) + \tau_{\text{loss}} c = s \quad \forall (x, y) \in \Omega, \forall t \in [0; T] $$

The coefficient $\mathbf{K}(x, y, t) [\text{m}^2 \text{s}^{-1}]$ is the diffusion tensor. This tensor models the turbulent diffusion at an unresolved scale (classically, scale smaller than the scale of the numerical grid) using the K theory, see e.g. (Panofsky et al., 1984; Sportisse, 2007). As mentioned above, in CTM, the wind velocity field $\vec{u}(x, y, t) = (u(x, y, t), v(x, y, t)) [\text{m s}^{-1}]$ is given as parameter, where $(u, v)$ are the velocity field components along $x$ and $y$. For example, $\vec{u}$ can be computed using a large (meso-scale) meteorological model and then, taken as input data in the present CTM, see an example in Sec. 2.3.2 and on Fig. 3a. The loss coefficient $\tau_{\text{loss}}(x, y, t) [\text{s}^{-1}]$ models the loss of pheromone due to vertical phenomena, including pheromone settlement on the ground and vertical diffusion spreading the pheromone out of the layer of interest. Hence, $\tau_{\text{loss}}$ depends on the land cover and on the stability of the atmosphere, see estimations of this coefficient in Sec. B.1.3. Other features can also be included in this parameter, such as humidity-dependent deposition. In the following, the solution of the CTM model is also referred to as the pheromone plume.

In Eq. (1), the source term $s(x, y, t) [\text{g m}^{-2} \text{s}^{-1}]$ is the quantity of pheromone emitted by the insects located in $(x, y, t)$ per unit of area and time. We highlight that this source term $s$ is the quantity to infer to deduce the localization of the insects. Namely, this quantity can be decomposed in

$$ (2) \qquad s(x, y, t) = q(t)p(x, y, t) $$

with $q(t) [\text{g s}^{-1}]$ the pheromone emission rate, and $p(x, y, t) [\text{m}^{-2}]$ the insect density on the ground, so that there is a direct relationship between $s$ and the local presence of pests. The emission rate $q$ can incorporate biological information, such as insect-specific emission patterns.

The 2D model (1) is closed by the following initial and boundary conditions:

* a null initial condition $c(x, y, t = 0) = 0 \quad \forall (x, y) \in \Omega$,

* a null diffusive flux $\mathbf{K} \nabla c \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial \Omega$,

* a null convective influx $\vec{u} c \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial \Omega \cap \{(x, y) | \vec{u}(x, y, t) \cdot \vec{n} < 0\} \forall t \in [0; T]$ with $\vec{n}$ the outgoing normal vector.

The pheromone propagation model (1) is solved using an implicit upwind finite volume scheme on a Cartesian grid implemented in the git repository of the PheroSensor project at https://forgemia.inra.fr/pherosensor/pherosensor-toolbox and (Malou et al., 2024b) (see Sec. A). The companion code to reproduce the different computations and figures of this study is available at https://forgemia.inra.fr/pherosensor/companion-code-bi-da and (Malou et al., 2024a).

## 2.2. Variational Data Assimilation for insect localization

We aim to approximate insect risk maps by inferring the spatial distribution of emitted pheromones. This inverse problem consists in estimating the source term $s$ of the pheromone propagation
be computed separately, offline, and given as parameters of the CTM. It results that this (linear and offline-coupled) CTM is much less complex than the full (non-linear) reactive Navier-Stokes equations, or than a CTM coupled online with a fine scale Navier-Stokes model (Leelössy et al., 2014), while being a reliable approximation: it is able to model the main physical phenomena at play during the transportation of pheromones.

In the present context, only the lower layer of the atmosphere, where insects and sensors are located, is of interest. Hence, for sake of computation load and data availability in real applications, a 2D depth-integrated CTM is used to model the total concentration of pheromone over the depth of the layer $c(x, y, t)$ [g m$^{-2}$] on a domain $\Omega \in \mathbb{R}^2$ and during a time window $[0; T]$. This depth-integrated CTM results in a 2D reaction-convection-diffusion model, that can be derived from the classical 3D CTM, see Sec. A.1 in the annex for details, and reads:

$$ (1) \qquad \frac{\partial c}{\partial t} - \nabla \cdot (\mathbf{K}\nabla c) + \nabla \cdot (\vec{u}c) + \tau_{\text{loss}}c = s \quad \forall (x, y) \in \Omega, \forall t \in [0; T] $$

The coefficient $\mathbf{K}(x, y, t)$ [m$^2$ s$^{-1}$] is the diffusion tensor. This tensor models the turbulent diffusion at an unresolved scale (classically, scale smaller than the scale of the numerical grid) using the K theory, see e.g. (Panofsky et al., 1984; Sportisse, 2007). As mentioned above, in CTM, the wind velocity field $\vec{u}(x, y, t) = (u(x, y, t), v(x, y, t))$ [m s$^{-1}$] is given as parameter, where $(u, v)$ are the velocity field components along x and y. For example, $\vec{u}$ can be computed using a large (meso-scale) meteorological model and then, taken as input data in the present CTM, see an example in Sec. 2.3.2 and on Fig. 3a. The loss coefficient $\tau_{\text{loss}}(x, y, t)$ [s$^{-1}$] models the loss of pheromone due to vertical phenomena, including pheromone settlement on the ground and vertical diffusion spreading the pheromone out of the layer of interest. Hence, $\tau_{\text{loss}}$ depends on the land cover and on the stability of the atmosphere, see estimations of this coefficient in Sec. A.1.3. Other features can also be included in this parameter, such as humidity-dependent deposition. In the following, the solution of the CTM model is also referred to as the pheromone plume.

In Eq. (1), the source term $s(x, y, t)$ [g m$^{-2}$ s$^{-1}$] is the quantity of pheromone emitted by the insects located in $(x, y, t)$ per unit of area and time. We highlight that this source term $s$ is the quantity to infer to deduce the localization of the insects. Namely, this quantity can be decomposed in

$$ (2) \qquad s(x, y, t) = q(t)p(x, y, t) $$

with $q(t)$ [g s$^{-1}$] the pheromone emission rate, and $p(x, y, t)$ [m$^{-2}$] the insect density on the ground, so that there is a direct relationship between $s$ and the local presence of pests. The emission rate $q$ can incorporate biological information, such as insect-specific emission patterns.

The 2D model (1) is closed by the following initial and boundary conditions:

* a null initial condition $c(x, y, t = 0) = 0 \quad \forall (x, y) \in \Omega$, as the initial time is set before any pheromone emission in a context of primary invasion,

* a null diffusive flux $\mathbf{K}\nabla c \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial\Omega$,

* a null convective influx $\vec{u}c \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial\Omega \cap \{(x, y)|\vec{u}(x, y, t) \cdot \vec{n} < 0\} \forall t \in [0; T]$ with $\vec{n}$ the outgoing normal vector.

The pheromone propagation model (1) is solved using an implicit upwind finite volume scheme on a Cartesian grid. This scheme has been implemented in the git repository of the PheroSensor project, see (Malou et al., 2024c) and Sec. 5. In this repository, one can find convergence tests of this numerical scheme. The companion code (Malou et al., 2024a) is also available to reproduce the different computations and figures of this study.
model (1), using time-series of pheromone concentration. To this purpose, we use Variational Data Assimilation (VDA) techniques.

**2.2.1. Minimization problem.** VDA consists in finding the value of the control variable, here $s$, that enables the prediction made by the direct model, here the pheromone propagation model (1), to be close to some observations, here the dataset of pheromone concentration time-series. Therefore, the inference problem is defined as a minimization problem over a cost function $j$ (Bocquet, 2014) that can be split in an observation and a regularization term:

$$ (3) \qquad \left\{ \begin{array}{l} \text{Find } s_a(x, y, t) \text{ such that:} \\ s_a(x, y, t) = \arg \min_{s(x,y,t)} j(s) \text{ with } j(s) = j_{\text{obs}}(s) + j_{\text{reg}}(s) \end{array} \right. $$

The term $j_{\text{obs}}$ aims to minimize the gap between model predictions and data (Bocquet, 2014; Le Dimet et al., 1986; Sasaki, 1958). It reads

$$ (4) \qquad j_{\text{obs}}(s) = \|m(c(s)) - m^{\text{obs}}\|_{\mathbb{R}^{-1}}^2 $$

where $m^{\text{obs}}$ are the observations, $c(s)$ is the direct model output obtained for a given control $s$, $m(c)$ the observation operator that maps the direct model state variable towards the observation space and $\mathbb{R}$ is the covariance operator of the observation error. Note that the equivalent statistical notation of eq. (4) reads

$$ m^{\text{obs}} = m(c(s)) + \varepsilon $$

where $\varepsilon$ is a gaussian noise with zero mean and covariance $\mathbb{R}$.

The regularization term $j_{\text{reg}}$ is commonly added in order to ensure the well-posedness of the optimization problem, by improving the convexity of $j$ (Bocquet, 2014). Classical regularization terms are based on a rough background estimate of the control variable, either by a model with an insufficient accuracy or by coarse data.

The novelty and the results of this study lies on the introduction of Biology-Informed (BI) regularization terms that enable to drive the DA method toward an optimum $s_a$ consistent with the biological knowledge introduced in $j_{\text{reg}}$, leading to a Biology-Informed Data Assimilation framework. In general terms, the regularization term reads:

$$ (5) \qquad j_{\text{reg}} = \sum_{i \in \mathcal{L}} \alpha_{\text{reg},i} j_{\text{reg},i} $$

with $\mathcal{L}$ the list of considered BI regularization terms $j_{\text{reg},i}$ and $\alpha_{\text{reg},i}$ the associated weight coefficient. This allows for composite regularization terms that take account of several biological information. The definition of these BI regularization terms is introduced in the Results in Sec. 3.1.

**2.2.2. Adjoint model for gradient computation.** To solve the minimization problem (3), we use a gradient-based method, suitable for our convex and differentiable optimization framework. Thus, at each descent step, the gradient $\nabla_{s_j} j$ must be computed and reads:

$$ \nabla_{s_j} j(s) = \nabla_{s_j} j_{\text{obs}}(s) + \nabla_{s_j} j_{\text{reg}}(s) $$

Note that we denote $\nabla_s$ the gradient with respect to the control $s$ to avoid confusion with $\nabla$ the classical gradient with respect to the space variables. Since $s$ is defined in a high-dimensional space and $c(s)$ involves the computation of a PDE model, the numerical computation of this gradient is challenging, in particular for its observation term $\nabla_{s_j} j_{\text{obs}}$. Noting the adjoint operator with
## 2.2. Variational Data Assimilation for insect localization

We aim to approximate insect risk maps by inferring the spatial distribution of emitted pheromones. This inverse problem consists in estimating the source term s of the pheromone propagation model (1), using time-series of pheromone concentration. To this purpose, we use Variational Data Assimilation (VDA) techniques.

### 2.2.1. Minimization problem.
VDA consists in finding the value of the control variable, $s(x, y, t)$ for all $(x, y, t) \in \Omega \times [0; T]$, that enables the prediction made by the direct model, here the pheromone propagation model (1), to be close to some observations, here the dataset of pheromone concentration time-series. Therefore, the inference problem is defined as a minimization problem of a cost function $j$ with respect to the control variable (Bocquet, 2014) that can be split in an observation and a regularization term:

$$ (3) \qquad \left\{ \begin{array}{l} \text{Find } s_a(x, y, t) \text{ such that:} \\ s_a(x, y, t) = \arg \min_{s(x,y,t)} j(s) \text{ with } j(s) = j_{\text{obs}}(s) + j_{\text{reg}}(s) \end{array} \right. $$

The term $j_{\text{obs}}$ aims to minimize the gap between model predictions and data (Bocquet, 2014; Le Dimet et al., 1986; Sasaki, 1958). It reads

$$ (4) \qquad j_{\text{obs}}(s) = \|m(c(s)) - m^{\text{obs}}\|_{\mathbb{R}^{-1}}^2 $$

where $m^{\text{obs}}$ are the observations, $c(s)$ is the direct model output obtained for a given control $s$, $m(c)$ the observation operator that maps the direct model state variable towards the observation space and $\mathbb{R}$ is the covariance operator of the observation error. Note that the equivalent statistical notation of eq. (4) reads

$$ m^{\text{obs}} = m(c(s)) + \varepsilon $$

where $\varepsilon$ is a gaussian noise with zero mean and covariance $\mathbb{R}$.

The regularization term $j_{\text{reg}}$ is commonly added in order to ensure the well-posedness of the optimization problem, by improving the convexity of $j$ (Bocquet, 2014). Classical regularization terms are based on a rough background estimate of the control variable, either by a model with an insufficient accuracy or by coarse data.

The novelty and the results of this study lies on the introduction of Biology-Informed (BI) regularization terms that enable to drive the DA method toward an optimum $s_a$ consistent with the biological knowledge introduced in $j_{\text{reg}}$, leading to a Biology-Informed Data Assimilation framework. In general terms, the regularization term reads:

$$ (5) \qquad j_{\text{reg}} = \sum_{i \in \mathcal{L}} \alpha_{\text{reg},i} j_{\text{reg},i} $$

with $\mathcal{L}$ the list of considered BI regularization terms $j_{\text{reg},i}$ and $\alpha_{\text{reg},i}$ the associated weight coefficient. This allows for composite regularization terms that take account of several biological inputs. The definition of these BI regularization terms is introduced in the Results in Sec. 3.1.

### 2.2.2. Adjoint model for gradient computation.
To solve the minimization problem (3), we use a gradient-based method, suitable for our convex and differentiable optimization framework. Thus, at each descent step, the gradient $\nabla_{s,j}$ must be computed and reads:

$$ \nabla_{s,j}(s) = \nabla_{s} j_{\text{obs}}(s) + \nabla_{s} j_{\text{reg}}(s) $$
*, we compute this gradient by using the adjoint model of the CTM (Lions, 1971; Pudykiewicz, 1998) that reads:

$$ (6) \quad \partial_t c^* + \nabla \cdot (\mathbf{K}^T \nabla c^*) + \nabla(\vec{u}c^*) - (\nabla \cdot \vec{u})c^* - \tau_{loss}c^* = \left(\frac{dm}{dc}(c(s))\right)^* \cdot 2\mathbf{R}^{-1} (m(c(s)) - m^{obs}) $$

with $c^*$ the adjoint state.

The adjoint model (6) is closed by the following final and boundary conditions:

* a null final condition $c^*(x, y, t = T) = 0 \forall (x, y) \in \Omega$,

* a null diffusive flux $\mathbf{K}^T \nabla c^* \cdot \vec{n} = 0 \forall (x, y) \in \partial\Omega$,

* a null outgoing convective flux $\vec{u}c^* \cdot \vec{n} = 0 \forall (x, y) \in \partial\Omega \cap \{(x, y) | \vec{u}(x, y, t) \cdot \vec{n} > 0\} \forall t \in ]0; T]$ with $\vec{n}$ the outgoing normal vector.

Further details on the derivation of the adjoint model (6) can be found in Sec. B.3.1. This adjoint model (6) can be interpreted as the retro-propagation with diffusion of the error between the prediction of the observed variable and the observations $m(c(s)) - m^{obs}$ from the final time $t = T$ to the initial time $t = 0$. The output of the adjoint model is also referred to as the retro-plume in the following. Let us note that the adjoint model (6) requires the observation operator to be differentiable.

Once the adjoint state is computed, the gradient $\nabla_{s,j}$ is computed using the expression (see more details in Sec. B.3.1):

$$ (7) \qquad \nabla_{s,j}(s) = \nabla_{s,j_{reg}}(s) - c^*(s) $$

The gradient-descent algorithm is initiated with a null control $s_0 = 0$.

The adjoint model (6) is solved using an implicit downwind finite volume scheme on the same Cartesian grid as the direct problem. A special care has been put on the time discretization so that the discretized adjoint problem is effectively the adjoint of the discretized direct problem. The gradient-descent algorithm is implemented in a custom script. The whole VDA process is implemented in the git repository of the PheroSensor project (see Sec. A).

### 2.2.3. Sensor contribution to gradient evaluation.
We stress that in our context of a network of independent sensors dispatched over the landscape, we can assume that the observation error covariance $\mathbf{R}$ is block-diagonal. Moreover, let us note that the left-hand side of Eq. (6) is linear by definition of the adjoint model. In addition, we assume that the observations operator can be split into a sum of one-sensor observation operator $m(c) = \sum_i m_i(c)$ (avoiding the degenerate case of overlapping sensors). From this assumption, we can inject $m = \sum_i m_i$ in the adjoint model (6) and split the source term to get the following one-sensor adjoint model:

$$ (8) \quad \partial_t c_i^* + \nabla \cdot (\mathbf{K}^T \nabla c_i^*) + \nabla(\vec{u}c_i^*) - (\nabla \cdot \vec{u})c_i^* - \tau_{loss}c_i^* = (d_c m_i(c(s)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s)) - m_i^{obs}) $$

with $c_i^*$ the one-sensor adjoint state related to the $i$th sensor. Since the adjoint model (6) is linear, it results that the adjoint state $c^*$ can also be expressed as the sum of the $c_i^*$:

$$ (9) \qquad c^* = \sum_i c_i^* $$

It results that $c_i^*$ is the contribution of the $i$th sensor to $c^*$, and thus to the gradient of $j_{obs}$. This observation will be used later to compare the contribution of each sensor to the data assimilation. More details on the derivation can be found in the annex Sec. B.3.1. Let us note that in our case, the assumption that the observations operator can be split into a sum of one-sensor observation operator is true, see Sec. 2.3.1.
Note that we denote $\nabla_s$ the gradient with respect to the control $s$ to avoid confusion with $\nabla$ the classical gradient with respect to the space variables. Since $s$ is defined in a high-dimensional space and $c(s)$ involves the computation of a PDE model, the numerical computation of this gradient is challenging, in particular for its observation term $\nabla_{s,j_{obs}}$. Noting the adjoint operator with $*$, we compute this gradient by using the adjoint model of the CTM (Lions, 1971; Pudykiewicz, 1998) that reads:

$$ (6) \quad \partial_t c^* + \nabla \cdot (\mathbf{K}^T \nabla c^*) + \nabla(\vec{u}c^*) - (\nabla \cdot \vec{u})c^* - \tau_{loss}c^* = \left(\frac{dm}{dc}(c(s))\right)^* \cdot 2\mathbf{R}^{-1} (m(c(s)) - m^{obs}) $$
with $c^*$ the adjoint state.

The adjoint model (6) is closed by the following final and boundary conditions:

* a null final condition $c^*(x, y, t = T) = 0 \quad \forall (x, y) \in \Omega$,

* a null diffusive flux $\mathbf{K}^T \nabla c^* \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial\Omega$,

* a null outgoing convective flux $\vec{u}c^* \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial\Omega \cap \{(x, y) | \vec{u}(x, y, t) \cdot \vec{n} > 0\} \quad \forall t \in ]0; T[$ with $\vec{n}$ the outgoing normal vector.

Further details on the derivation of the adjoint model (6) can be found in Sec. A.3.1. This adjoint model (6) can be interpreted as the retro-propagation with diffusion of the error between the prediction of the observed variable and the observations $m(c(s)) - m^{obs}$ from the final time $t = T$ to the initial time $t = 0$. The output of the adjoint model is also referred to as the retro-plume in the following. Let us note that the adjoint model (6) requires the observation operator to be differentiable.

Once the adjoint state is computed, the gradient $\nabla_{s,j}$ is computed using the expression (see more details in Sec. A.3.1):

$$ (7) \qquad \nabla_{s,j}(s) = \nabla_{s,j_{reg}}(s) - c^*(s) $$

The gradient-descent algorithm is initiated with a null control $s_0 = 0$.

The adjoint model (6) is solved using an implicit downwind finite volume scheme on the same Cartesian grid as the direct problem. The whole VDA process is implemented in the git repository of the PheroSensor project, see (Malou et al., 2024c) and Sec. 5, including the solver of the adjoint model and a gradient-descent algorithm. In this repository, one can find convergence tests of the numerical scheme of the adjoint model, test of the adjoint operators as well as a validation of the estimation of the gradient by comparison with finite differences.

### 2.2.3. Sensor contribution to gradient evaluation.
We stress that in our context of a network of independent sensors dispatched over the landscape, we can assume that the observation error covariance $\mathbf{R}$ is block-diagonal. Moreover, let us note that the left-hand side of Eq. (6) is linear by definition of the adjoint model. In addition, we assume that the observations operator can be split into a sum of one-sensor observation operator $m(c) = \sum_i m_i(c)$ (avoiding the degenerate case of overlapping sensors). From this assumption, we can inject $m = \sum_i m_i$ in the adjoint model (6) and split the source term to get the following one-sensor adjoint model:

$$ (8) \quad \partial_t c_i^* + \nabla \cdot (\mathbf{K}^T \nabla c_i^*) + \nabla(\vec{u}c_i^*) - (\nabla \cdot \vec{u})c_i^* - \tau_{loss}c_i^* = (d_c m_i(c(s)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s)) - m_i^{obs}) $$

with $c_i^*$ the one-sensor adjoint state related to the $i^{th}$ sensor. Since the adjoint model (6) is linear, it results that the adjoint state $c^*$ can also be expressed as the sum of the $c_i^*$:

$$ (9) \qquad c^* = \sum_i c_i^* $$
2.2.4. Stopping criteria. As per Morozov's discrepancy principle, the optimization algorithm should be continued until the gap between the observations and the system state achieves the level of the accuracy of the observations. Hence, the gradient descent algorithm is terminated and the optimum $s_a$ is assumed to be reached when the cost $j_{obs}$ attains a value lower than the noise norm $||\epsilon||^2_{L^2(\Omega_{obs})}$.

## 2.3. Pheromone sensors and test cases

2.3.1. Observation operator of the pheromone sensors. The state variable $c$ is a function defined over $[0, T] \times \Omega$. We assume that we have at disposal $N_s$ sensors that are dispatched throughout the landscape and remain in the same position $(x_i, y_i)$ all along the time domain. We furthermore assume that, given $[t_1^{obs}, \dots, t_{N_T}^{obs}]$ a time sampling of the time domain, every sensor produces a data at time $t_k^{obs}$, $k = 1, \dots, N_t$, so that the observation space ($\mathcal{O}$) has dimension $N_s \times N_t$. Let us note $e_{i,k}$ the vector, for sensor $i$ and time $k$, of the canonical basis of mathcal{O}.

Then, we define the observation operator $c \mapsto m(c)$ mapping the state variable $c$ towards the observation space $\mathcal{O}$ with the equation

$$ (10) \qquad m(c) = \sum_{i=1}^{N_s} m_i(c) \quad \text{with} \quad m_i(c) = \sum_{k=1}^{N_t} e_{i,k} \int_{t_k^{obs}-\delta_t}^{t_k^{obs}} \int_{\Omega} \delta_{(x_i, y_i)}(x, y) c(\tau, x, y) d\tau dxdy $$

where $\delta_{(x_i, y_i)}$ is the Dirac function at the localization of the $i^{th}$ sensor and $\delta_t$ a time window. For observation $k$ of sensor $i$, this expression averages over a time window of width $\delta_t$ before the observation time the pheromone concentration at location $(x_i, y_i)$. This models a pre-concentrator that actively filters air during a time window $\delta_t$ and releases the accumulated quantity into the sensor for the measure. Let us note that the time sampling of the observations must satisfy $\delta_t < t_1^{obs}$ so that $t_1^{obs} - \delta_t > 0$. Let us note that in this case, the operator $c \mapsto m(c)$ is linear and is split in the sum of one-sensor observation operators $m_i$, as mentioned in Sec. 2.2.3.

2.3.2. Setting up test cases. All along the paper, we will consider two test cases. The first one is a simplified yet illustrative scenario, referred to as the Toy case, that will be used to test and analyze different biology-informed regularization strategies. The second one depicts a realistic and representative scenario for invasion of an alien insect in an agricultural landscape. We chose the Fall ArmyWorm (FAW), a.k.a. Spodoptera frugiperda, as the model insect considered in this study. The FAW is a moth that belongs to the family Noctuidae. This insect pest induces the highest yield loss for maize (Savary et al., 2019), but can grow on hundreds of plant hosts, causing severe damage in grasses (IPPC Secretariat, 2021). It is expected that global warming reduces the populations in the southern America area due to warmer and drier climates (IPPC Secretariat, 2021), but could favour the emergence of the species in the European Union (EFSA PLH Panel et al., 2018). It is then a suitable test case for a realistic scenario in an epidemiological surveillance framework in Europe. In this section, we will detail the set up of the environmental and landscape parameters needed for solving the pheromone propagation model (1) in the two test cases. The inference method and the additional biology-informed regularization terms are tested using synthetic dataset generated from a reference source term, hereinafter referred to as the target source term $s_t$. Further details for data generation and inference will be added in this sections.

**Toy case.** The computation spatial domain of the Toy case is a rectangular grid of size $40m \times 60m$ (i.e. 0.24ha) over a time window of $T = 20s$. It is built up upon a fictive insect that constantly releases pheromone at an emission rate of $q = 2.7pg.s^{-1}$ per insect. The insects density is
It results that $c_i^*$ is the contribution of the $i^{th}$ sensor to $c^*$, and thus to the gradient of $j_{obs}$. This observation will be used later to compare the contribution of each sensor to the data assimilation. The comparison of the contribution of each sensor then enables to determine which sensors contribute the most and which contribute the less, and to propose strategies for repositioning the sensors with least contribution. Repositioning strategies based on such comparison will be used in the realistic case of this study, see more details in Sec. 2.3.2 and the results of the repositioning strategy in Sec. 3.6

More details on the derivation can be found in the annex Sec. A.3.1. Let us note that in our case, the assumption that the observations operator can be split into a sum of one-sensor observation operator is true, see Sec. 2.3.1.

### 2.2.4. Stopping criteria.
As per Morozov's discrepancy principle, the optimization algorithm should be continued until the gap between the observations and the system state achieves the level of the accuracy of the observations. Hence, the gradient descent algorithm is terminated and the optimum $s_a$ is assumed to be reached when the cost $j_{obs}$ attains a value lower than the noise norm $||\epsilon||_{L^2(\Omega_{obs})}^2$.

## 2.3. Pheromone sensors and test cases

### 2.3.1. Observation operator of the pheromone sensors.
The state variable $c$ is a function defined over $[0, T] \times \Omega$. We assume that we have at disposal $N_s$ sensors that are dispatched throughout the landscape and remain in the same position $(x_i, y_i)$ all along the time domain. We furthermore assume that, given $[t_1^{obs}, \dots, t_{N_T}^{obs}]$ a time sampling of the time domain, every sensor produces a data at time $t_k^{obs}$, $k = 1, \dots, N_t$, so that the observation space $\mathcal{O}$ has dimension $N_s \times N_t$. Let us note $e_{i,k}$ the vector, for sensor $i$ and time $k$, of the canonical basis of $\mathcal{O}$.

Then, we define the observation operator $c \mapsto m(c)$ mapping the state variable $c$ towards the observation space $\mathcal{O}$ with the equation

$$ (10) \qquad m(c) = \sum_{i=1}^{N_s} m_i(c) \quad \text{with} \quad m_i(c) = \sum_{k=1}^{N_t} e_{i,k} \int_{t_k^{obs}-\delta_t}^{t_k^{obs}} \int_{\Omega} \delta_{(x_i, y_i)}(x, y) c(\tau, x, y) d\tau dxdy $$

where $\delta_{(x_i, y_i)}$ is the Dirac function at the localization of the $i^{th}$ sensor and $\delta_t$ a time window. For observation $k$ of sensor $i$, this expression averages over a time window of width $\delta_t$ before the observation time the pheromone concentration at location $(x_i, y_i)$. This models a pre-concentrator that actively filters air during a time window $\delta_t$ and releases the accumulated quantity into the sensor for the measure. Let us note that the time sampling of the observations must satisfy $\delta_t < t_1^{obs}$ so that $t_1^{obs} - \delta_t > 0$. Let us note that in this case, the operator $c \mapsto m(c)$ is linear and is split in the sum of one-sensor observation operators $m_i$, as mentioned in Sec. 2.2.3.

### 2.3.2. Setting up test cases.
All along the paper, we will consider two test cases. The first one is a simplified yet illustrative scenario, referred to as the Toy case, that will be used to test and analyze different biology-informed regularization strategies. The second one depicts a realistic and representative scenario for invasion of an alien insect in an agricultural landscape. We chose the fall armyworm (FAW), a.k.a. *Spodoptera frugiperda*, as the model insect considered in this study. The FAW is a moth that belongs to the family Noctuidae. This insect pest induces the highest yield loss for maize (Savary et al., 2019), but can grow on hundreds of plant hosts, causing severe damage in grasses (IPPC Secretariat, 2021). It is expected that global warming reduces the populations in the southern America area due to warmer and drier climates (IPPC Secretariat, 2021), but could favour the emergence of the species in the European Union (EFSA PLH Panel
constant $p_0^{\text{toy}} = 0.2m^{-2}$ in a fixed area of the domain located in a circular support $\Omega_{\text{supp}} = S$ centred in $(x_s, y_s) = (10, 10)$ and of radius $R = 5m$, see red circle in Fig. 2, which corresponds to approximatively 16 individuals. The wind field $u$ is constructed from a constant wind speed $|\bar{u}| = 8m.s^{-1}$ and a spatially-homogeneous time-dependent wind direction that linearly varies from $\theta = \frac{\pi}{4}$ at $t = 0$ to $\theta = \frac{\pi}{2}$ at $t = T$. The loss coefficient $\tau_{\text{loss}}$ is computed to approximate realistic vertical fluxes in a 3D simplified CTM model (see more details in Sec. B.1), resulting in a constant value of $\tau_{\text{loss}} = 0.01s^{-1}$. The diffusion tensor $K$ is assumed to be isotropic and homogeneous over time and space. As the diffusion value $K = 50m^2s^{-1}$ is a typical value to model an unstable atmosphere (BEA Fisher et al., 1981), we use the value $K = 10m^2s^{-1}$, which corresponds to a moderately unstable atmosphere.

$N_s = 35$ synthetic sensors are randomly placed each at a different position in the landscapes with a uniform distribution, see Fig. 2. Moreover, the accumulation time step is set to $\delta_t = 2s$. Noise is subsequently added. We take $m_{\text{obs}}(t_i, x_j, y_j) = \max(0, m(c(s_t); t_i, x_j, y_j) + \epsilon_{i,j})$, in order to avoid negative values, with $c(s_t)$ the results of the direct model (1) computed from the source term $s_t$ and $\epsilon_{i,j} \sim N(0, \sigma_{i,j}^2)$ a Gaussian noise. Since pre-concentration uncertainties scale with pheromone density, we scale the standard deviation with the model observations: $\sigma_{i,j} = 0.01m(c(s_t); t_i, x_j, y_j)$.

A favorable habitat, used in the Tikhonov regularization, see Sec. 3.1.1, is defined as a square with a width of 20m, overlapping the source support $S$, see zone delimited by the dotted red line on Fig. 2.

<figure><img src="image_3.png" /><figcaption>**Figure 2** – Map of the target source term $s_t$, circular support (constant within the red circle), the favorable habitat (square zone delimited by the red dotted line) and the sensor location (colored markers) for the Toy test case.</figcaption></figure>

**Realistic set up: the FAW case.** The FAW case is based on a real landscape in central France (Adamczyk-Chauvat et al., 2024), see Fig. 3a for the computational domain (dashed lines). The same homogeneous diffusion tensor ($K = 10m^2s^{-1}$) than the Toy case is used. The wind field $\bar{u}$ is recovered from outputs of the French weather forecast model Arome, which is the regional numerical weather forecast model from MétéoFrance (French weather forecast service, data and documentation available at https://donneespubliques.meteofrance.fr/?fond=produit&id_produit=131&id_rubrique=51) corresponding to the actual wind field in area of interest at noon
et al., 2018). It is then a suitable test case for a realistic scenario in an epidemiological surveillance framework in Europe. In this section, we will detail the set up of the environmental and landscape parameters needed for solving the pheromone propagation model (1) in the two test cases. The inference method and the additional biology-informed regularization terms are tested using synthetic dataset generated from a reference source term, hereinafter referred to as the target source term $s_t$. Further details for data generation and inference will be added in this section.

**Toy case.** The computation spatial domain of the Toy case is a rectangular grid of size 40m × 60m (i.e. 0.24ha) over a time window of $T = 20$s. It is built up upon a fictive insect that constantly releases pheromone at an emission rate of $q = 2.7\text{pg}\cdot\text{s}^{-1}$ per insect. The insects density is constant $p_0^{\text{toy}} = 0.2\text{m}^{-2}$ in a fixed area of the domain located in a circular support $\Omega_{\text{supp}} = S$ centred in $(x_s, y_s) = (10, 10)$ and of radius $R = 5\text{m}$, see red circle in Fig. 2, which corresponds to approximatively 16 individuals. The wind field $u$ is constructed from a constant wind speed $|\vec{u}| = 8\text{m}\cdot\text{s}^{-1}$ and a spatially-homogeneous time-dependent wind direction that linearly varies from $\theta = \frac{\pi}{4}$ at $t = 0$ to $\theta = \frac{\pi}{2}$ at $t = T$. The loss coefficient $\tau_{\text{loss}}$ is computed to approximate realistic vertical fluxes in a 3D simplified CTM model (see more details in Sec. A.1), resulting in a constant value of $\tau_{\text{loss}} = 0.01\text{s}^{-1}$. The diffusion tensor $\mathbf{K}$ is assumed to be isotropic and homogeneous over time and space. As the diffusion value $\mathbf{K} = 50\text{m}^2\text{s}^{-1}$ is a typical value to model an unstable atmosphere (BEA Fisher et al., 1981), we use the value $\mathbf{K} = 10\text{m}^2\text{s}^{-1}$, which corresponds to a moderately unstable atmosphere. A favorable habitat, used in the Tikhonov regularization, see Sec. 3.1.1, is defined as a square with a width of 20m, overlapping the source support $S$, see zone delimited by the dotted red line on Fig. 2.

$N_s = 35$ synthetic sensors are randomly placed each at a different position in the landscapes with a uniform distribution, see Fig. 2. Moreover, the accumulation time step is set to $\delta_t = 2$s. Noise is subsequently added. We take $m_{\text{obs}}(t_i, x_j, y_j) = \max(0, m(c(s_t); t_i, x_j, y_j) + \epsilon_{i,j})$, in order to avoid negative values, with $c(s_t)$ the results of the direct model (1) computed from the source term $s_t$ and $\epsilon_{i,j} \sim N(0, \sigma_{i,j}^2)$ a Gaussian noise. Since pre-concentration uncertainties scale with pheromone density, we scale the standard deviation with the model observations: $\sigma_{i,j} = 0.01m(c(s_t); t_i, x_j, y_j)$.

To discretize the spatio-temporal domain, fixed space and time steps are used, resp. $\Delta x = \Delta y = 0.5\text{m}$ and $\Delta t = 0.1\text{s}$. This leads to a source term $s(x, y, t)$ to infer with 1920000 degree of freedom.

**Realistic set up: the FAW case.** The FAW case is based on a real landscape in central France (Adamczyk-Chauvat et al., 2024), see Fig. 3a for the computational domain (dashed lines). The same homogeneous diffusion tensor ($\mathbf{K} = 10\text{m}^2\text{s}^{-1}$) as in the Toy case is used. The wind field $\vec{u}$ is recovered from outputs of the French weather forecast model Arome, which is the regional numerical weather forecast model from MétéoFrance (French weather forecast service, data and documentation available at https://donneespubliques.meteofrance.fr/?fond=produit&id produit=131&id_rubrique=51) corresponding to the actual wind field in area of interest at noon the 2022/11/16. As the spatial sampling of the Arome model is about $0.01^\circ$ ($\approx 1.3\text{km}$), see Fig. 3a, the wind field is projected on the computation grid by linear interpolation, and kept stationary over the computation time window.

As the loss coefficient $\tau_{\text{loss}}$ depends on the land cover, the land use heterogeneities in the landscape are taken into account in the FAW case using land cover data, see Fig. 3b. The different values of $\tau_{\text{loss}}$ depending on the land cover as recapitulated in Tab. 1, see Sec. A.1 for computations. Note that in the Toy case, the constant value 0.01 corresponds to bare soil.
the 2022/11/16. As the spatial sampling of the Arome model is about 0.01° (≈ 1.3km), see Fig. 3a, the wind field is projected on the computation grid by linear interpolation, and kept stationary over the computation time window.

As the loss coefficient $\tau_{loss}$ depends on the land cover, the land use heterogeneities in the landscape are taken into account in the FAW case using land cover data, see Fig. 3b. The different values of $\tau_{loss}$ depending on the land cover as recapitulated in Tab. 1, see Sec. B.1 for computations. Note that in the Toy case, the constant value 0.01 corresponds to bare soil.

<table><caption>**Table 1** – Loss coefficient $\tau_{loss}$ values depending on the land use. Sum. crop = Summer crop, wint. cereals = winter cereals.</caption><thead><tr><th>land use</th><th>sum. crop</th><th>maize</th><th>wint. cereals</th><th>grassland</th><th>forest</th><th>bare soil</th><th>urban</th><th>road</th></tr></thead><tbody><tr><td>τ<sub>loss</sub> (s<sup>-1</sup>)</td><td>0.0097</td><td>0.011</td><td>0.0105</td><td>0.0099</td><td>0.0095</td><td>0.01</td><td>0.009</td><td>0.009</td></tr></tbody></table>

FAW has a specific pattern of pheromone emission. During the mating period, the FAW females search for a suitable location and stay there until they lay their eggs. Each night, they emit a sexual pheromone that attracts the males. Pheromone emission is stopped just at the beginning of the mating process, so that at the population level, the pheromone emission exponentially decreases while females meet a mating partner and stop emitting. This pattern is repeated every night. We start the time window for pheromone recording 4h after nightfall, when emissions are at their peak (Tumlinson et al., 1986).

Recall that $s = pq$ with $p$ the density of insects and $q$ the quantity of pheromone emitted per time unit per insect. We use a constant emission rate per insect $q = 2.7 \text{pg} \cdot \text{s}^{-1}$. The reference density $p_{ref}^{FAW}$ of emitter insects is chosen to satisfy a realistic population dynamics model, see later Sec. 3.7, and is defined as

$$ (11) \qquad p_{ref}^{FAW}(t, x, y) = p_{0}^{FAW}(x, y) \exp \left( - \int_{0}^{t} \gamma(s) ds \right) $$

with $\gamma(t)$ the mating rate function representing the proportion of female that encounter a male at a given time and stop emitting (see Fig. 3c) and $p_0^{FAW}$ the initial density of females. In that regard, $t \mapsto \exp(-\int_0^t \gamma(s) ds)$ is the function representing the proportion of females still available for mating and that keep emitting at a given time. Initially, we suppose that approximately 60 females are uniformly distributed within a circular domain of radius $10m$. This corresponds to a constant initial density of $\bar{p}_0 = 0.2m^{-2}$ within $S = \{(x, y) \in \Omega \mid (x - x_s)^2 + (y - y_s)^2 < R^2\}$ a circular support centred in $(x_s, y_s)$ and of radius $R = 10m$, see the red circular domain on Fig. 3b:

$$ (12) \qquad p_0^{FAW}(x, y) = \bar{p}_0 1_S(x, y). $$

The FAW females usually emit when laying on plants favorable for the larva. In the FAW case landscape, we assume that maize is the most favorable plant. Therefore, the support $S$ is located within a domain region covered by maize, see land cover in dark orange on Fig. 3b. Moreover, we assume that the mating rate function $\gamma$ is a logistic function, see blue curve on Fig. 3c:

$$ (13) \qquad \gamma(t) = \frac{K_{mat}}{1 + \exp(-\lambda(t - \tau))} $$

where $K_{mat}$ is the maximal mating rate, $\tau$ the midpoint time value and $\lambda$ the logistic growth rate (see Fig. 3c). Let us note that $\tau$ is related to the average time a male necessitates to reach a female while tracking back the pheromone source. Based on biological expertise, $\tau$ is taken such that 50% of the females have mated with a male within the 30 first minutes. Hence, $p_{ref}^{FAW}(t = 30 \text{min}, x, y) = \frac{1}{2} p_0^{FAW}(x, y)$. Using the mating maximal rate $K_{mat} = 0.1$ and a logistic growth rate $\lambda \approx 8.3 \times 10^{-3}$, it results that $\tau = 127.94 \text{min}$.
<figure><img src="image_3.png" /><figcaption>**Figure 2** – Map of the target source term $s_t$ circular support (constant within the red circle), the favorable habitat (square zone delimited by the red dotted line) and the sensor location (colored markers) for the Toy test case.</figcaption></figure>

<table><caption>**Table 1** – Loss coefficient $τ_{loss}$ values depending on the land use. Sum. crop = Summer crop, wint. cereals = winter cereals.</caption><thead><tr><th>land use</th><th>sum. crop</th><th>maize</th><th>wint. cereals</th><th>grassland</th><th>forest</th><th>bare soil</th><th>urban</th><th>road</th></tr></thead><tbody><tr><td>τ<sub>loss</sub> (s<sup>-1</sup>)</td><td>0.0097</td><td>0.011</td><td>0.0105</td><td>0.0099</td><td>0.0095</td><td>0.01</td><td>0.009</td><td>0.009</td></tr></tbody></table>

FAW has a specific pattern of pheromone emission. During the mating period, the FAW females search for a suitable location and stay there until they lay their eggs. Each night, they emit a sexual pheromone that attracts the males. Pheromone emission is stopped just at the beginning of the mating process, so that at the population level, the pheromone emission exponentially decreases while females meet a mating partner and stop emitting. This pattern is repeated every night. We start the time window for pheromone recording three hours after nightfall, when emissions are at their peak (Tumlinson et al., 1986).

Recall that $s = pq$ with $p$ the density of insects and $q$ the quantity of pheromone emitted per time unit per insect. We use a constant emission rate per insect $q = 2.7pg.s^{-1}$. The reference density $p_{ref}^{FAW}$ of emitter insects is chosen to satisfy a realistic population dynamics model, see later Sec. 3.7, and is defined as

$$ (11) \qquad p_{ref}^{FAW}(t, x, y) = p_{0}^{FAW}(x, y) \exp \left( - \int_{0}^{t} \gamma(s) ds \right) $$

with $\gamma(t)$ the mating rate function representing the proportion of female that encounter a male at a given time and stop emitting (see Fig. 3c) and $p_0^{FAW}$ the initial density of females. In that regard, $t \mapsto \exp(-\int_0^t \gamma(s) ds)$ is the function representing the proportion of females still available for mating and that keep emitting at a given time. Initially, we suppose that approximately 60 females are uniformly distributed within a circular domain of radius $10m$. This corresponds to a constant initial density of $\bar{p}_0 = 0.2m^{-2}$ within $S = \{(x, y) \in \Omega | (x - x_s)^2 + (y - y_s)^2 < R^2\}$ a circular support centred in $(x_s, y_s)$ and of radius $R = 10m$, see the red circular domain on Fig. 3b:

$$ (12) \qquad p_{0}^{FAW}(x, y) = \bar{p}_{0} \mathbb{1}_{S}(x, y). $$
<figure><img src="image_4.png" /></figure>

(a) Wind field as obtained from the weather forecast model, before interpolation on the computing grid (length and color of the arrows normalized by the wind speed). Nine vectors in total are available, at the center of the figure, at the corners and in the cardinal directions.

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

(b) Fragment of agricultural landscape, land cover and sensor positions (see sec. 3.6 and fig 13).

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

(c) FAW test case, mating rate $\gamma(t)$ as in Eq. (13) (blue) and total proportion of female still emitting pheromones $\exp(-\int_{0}^{t} \gamma(s)ds)$ (red, see Eq. (11)) vs time $t$.

**Figure 3** – Map of the distribution of the computations domains (dashed lines, black for the total domain and red for the domain of the inference) with the wind field computed by the Arome model at a given time (left) and the division of agricultural parcels and soil occupation (right) for the FAW case.

In a first time, 15 sensors are randomly dispatched with uniform distribution in the maize field of interest, see the land cover in dark orange on Fig. 3b. Then, the 6 sensors that contribute the less (6 last in the contribution order for the one-sensor criteria defined in Sec. 2.4.3) are randomly relocated inside the convex hull of the 9 remaining sensors. This relocation method is detailed in the Results in sec. 3.6 and illustrated in Fig. 13.

305

Note that the computation domain is narrowed around the favorable zone of maize field to reduce computation time, see domain delimited by the red dashed line on Fig. 3a.

The accumulation time step is set to $\delta_t = 30min$, in the order of reasonable pre-concentration time. Noise is subsequently added with the same noise model than in the Toy case: $m_{obs}(t_i, x_j, y_j) = \max(0, m(c(s_t); t_i, x_j, y_j) + \epsilon_{i,j})$ with $c(s_t)$ the results of the direct model (1) computed from the source term $s_t$ and $\epsilon_{i,j} \sim N(0, \sigma_{i,j}^2)$ a Gaussian noise with standard deviation $\sigma_{i,j} = 0.05m(c(s_t); t_i, x_j, y_j)$.

310
<figure><img src="image_4.png" /></figure>

(a) Wind field as obtained from the weather forecast model, before interpolation on the computing grid (length and color of the arrows normalized by the wind speed). Nine vectors in total are available, at the center of the figure, at the corners and in the cardinal directions.

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

(b) Fragment of agricultural landscape, land cover and sensor positions (see sec. 3.6 and fig 13).

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

(c) FAW test case, mating rate $\gamma(t)$ as in Eq. (13) (blue) and total proportion of female still emitting pheromones $\exp(-\int_0^t \gamma(s)ds)$ (red, see Eq. (11)) vs time $t$.

<caption>**Figure 3** – Map of the domain and red for the domain of the inference) with the wind field computed by the Arome model at a given time (left) and the division of agricultural parcels and soil occupation (right) for the FAW case.</caption>

The FAW females usually emit when laying on plants favorable for the larva. In the FAW case landscape, we assume that maize is the most favorable plant. Therefore, the support $S$ is located within a domain region covered by maize, see land cover in dark orange on Fig. 3b. Moreover, we assume that the mating rate function $\gamma$ is a logistic function, see blue curve on Fig. 3c:

$$ (13) \qquad \gamma(t) = \frac{K_{mat}}{1 + \exp(-\lambda(t - \tau))} $$

where $K_{mat}$ is the maximal mating rate, $\tau$ the midpoint time value and $\lambda$ the logistic growth rate (see Fig. 3c). Let us note that $\tau$ is related to the average time a male necessitates to reach a female while tracking back the pheromone source. Based on biological expertise, $\tau$ is taken such that 50% of the females have mated with a male within the 30 first minutes. Hence, $p_{ref}^{FAW}(t = 30 \text{min}, x, y) = \frac{1}{2}p_0^{FAW}(x, y)$. Using the mating maximal rate $K_{mat} = 0.1$ and a logistic growth rate $\lambda \approx 8.3 \times 10^{-3}$, it results that $\tau = 127.94 \text{min}$.
## 2.4. Benchmarking criteria

Once the optimization algorithm is terminated, the optimum $s_a$ is compared against the target $s_t$, that is used to generate the data, using the criteria presented below to evaluate the accuracy of the localization of the insects.

**2.4.1. RMSE-type criteria.** A first class of criteria is based on the Root x Mean Square Error (RxMSE) between the target $s_t$ and a value of $s$, that measures the accuracy of the data as-simulation process. The x indicates here that we average over the spatial dimension only. In this study, the errors of interest are the quadratic error on the whole domain $\Omega$ and on $\Omega_{supp}$, the support of the target $s_t$. We set the total RxMSE criteria by

$$ (14) \qquad R \times \text{MSE}_{\text{tot}}(t; s) = \sqrt{\frac{1}{|\Omega|} \int_{\Omega} (s(X, t) - s_t(X, t))^2 dX} $$

and the support RMSE criteria by

$$ (15) \qquad R \times \text{MSE}_{\text{supp}}(t; s) = \sqrt{\frac{1}{|\Omega_{\text{supp}}|} \int_{\Omega_{\text{supp}}} (s(X, t) - s_t(X, t))^2 dX} $$

While the error on the whole domain gives a global assessment of the accuracy of $s_a$, by also measuring false emissions estimated outside of the target support, the error on the target support enables to quantify the error made on the specific zone where the pheromone are emitted, and thus where the insects are located.

**2.4.2. False discovery rate and miss rate.** The presence of insects at a given time and point of the domain is assessed when the quantity of pheromone emitted is above a given threshold $\epsilon_s$. This threshold depicts the minimal quantity of pheromone emitted by a single insect over the area of a numerical mesh cell. According to this presence criteria, the subdomain where the insects are located given a source term is the superlevel set

$$ (16) \qquad \Omega_{sls}(t; s, \epsilon_s) = \{(x, y) \in \Omega \mid s(x, y, t) > \epsilon_s\} $$

This implies that the superlevel set $\Omega_{sls}(t; s_a, \epsilon_s)$ (set in blue dashed lines on Fig. 4) is the sub-

<table><thead><tr><th>prediction<br>truth</th><th>presence<br>Ω<sub>sls</sub>(s<sub>a</sub>)</th><th>absence<br>Ω&#773;<sub>sls</sub>(s<sub>a</sub>)</th></tr></thead><tbody><tr><th>presence<br>Ω<sub>sls</sub>(s<sub>t</sub>)</th><td>true presence prediction<br>Ω<sub>tp</sub> = Ω<sub>sls</sub>(s<sub>a</sub>) ∩ Ω<sub>sls</sub>(s<sub>t</sub>)</td><td>false presence prediction<br>Ω<sub>fp</sub> = Ω<sub>sls</sub>(s<sub>a</sub>) ∩ Ω&#773;<sub>sls</sub>(s<sub>t</sub>)</td></tr><tr><th>absence<br>Ω<sub>sls</sub>(s<sub>t</sub>)</th><td>false absence prediction<br>Ω<sub>fa</sub> = Ω&#773;<sub>sls</sub>(s<sub>a</sub>) ∩ Ω<sub>sls</sub>(s<sub>t</sub>)</td><td>true absence prediction<br>Ω<sub>ta</sub> = Ω<sub>sls</sub>(s<sub>a</sub>) ∩ Ω&#773;<sub>sls</sub>(s<sub>t</sub>)</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4** – Illustration of the predicted and true superlevel sets and the false absence and presence predictions</figcaption></figure>

domain where the insects are predicted to be present. On the other hand, the superlevel set $\Omega_{sls}(t; s_t, \epsilon_s)$ (set in red dashed lines on Fig. 4) is the subdomain where the insects are truly present.

These two sets will be compared using dissimilarity criteria in a same way as the Jaccard similarity coefficient or the confusion matrix in supervised statistical learning, see e.g. (Stehman, 1997). The subdomain $\Omega_{fp}(t; \epsilon_s) = \Omega_{sls}(t; s_a, \epsilon_s) \cap \overline{\Omega_{sls}(t; s_t, \epsilon_s)}$ corresponds to the subdomain
To begin with, 15 sensors are randomly dispatched with uniform distribution in the maize field of interest, see the land cover in dark orange on Fig. 3b. Then, the 6 sensors that contribute the less (6 last in the contribution order for the one-sensor criteria defined in Sec. 2.4.3) are randomly relocated inside the convex hull of the 9 remaining sensors. This relocation method is detailed in the results in Sec. 3.6 and illustrated in Fig. 13.

Note that the computation domain is narrowed around the 300m long and 400m wide favorable zone of maize field to reduce computation time, see domain delimited by the red dashed line on Fig. 3a. The accumulation time step is set to $\delta_t = 30min$, in the order of reasonable pre-concentration time. Noise is subsequently added with the same noise model than in the Toy case: $m_{obs}(t_i, x_j, y_j) = max(0, m(c(s_t); t_i, x_j, y_j) + \epsilon_{i,j})$ with $c(s_t)$ the results of the direct model (1) computed from the source term $s_t$ and $\epsilon_{i,j} \sim N(0, \sigma_{i,j}^2)$ a Gaussian noise with standard deviation $\sigma_{i,j} = 0.05m(c(s_t); t_i, x_j, y_j)$.

To discretized the spatio-temporal domain of this case, fixed space and time steps are used, resp. $\Delta x = \Delta y = 2m$ and $\Delta t = 30s$. This leads to a source term $s(x, y, t)$ to infer with 10800000 degree of freedom.

## 2.4. Benchmarking criteria

Once the optimization algorithm is terminated, the optimum $s_a$ is compared against the target $s_t$, that is used to generate the data, using the criteria presented below to evaluate the accuracy of the localization of the insects.

### 2.4.1. RMSE-type criteria.
A first class of criteria is based on the Root x Mean Square Error (RxMSE) between the target $s_t$ and a value of $s$, that measures the accuracy of the data assimilation process. The x indicates here that we average over the spatial dimension only. In this study, the errors of interest are the quadratic error on the whole domain $\Omega$ and on $\Omega_{supp}$, the support of the target $s_t$. We set the total RxMSE criteria by

$$ (14) \qquad R \times \text{MSE}_{\text{tot}}(t; s) = \sqrt{\frac{1}{|\Omega|} \int_{\Omega} (s(X, t) - s_t(X, t))^2 dX} $$

and the support RMSE criteria by

$$ (15) \qquad R \times \text{MSE}_{\text{supp}}(t; s) = \sqrt{\frac{1}{|\Omega_{\text{supp}}|} \int_{\Omega_{\text{supp}}} (s(X, t) - s_t(X, t))^2 dX} $$

While the error on the whole domain gives a global assessment of the accuracy of $s_a$, by also measuring false emissions estimated outside of the target support, the error on the target support enables to quantify the error made on the specific zone where the pheromone are emitted, and thus where the insects are located.

### 2.4.2. False discovery rate and miss rate.
The presence of insects at a given time and point of the domain is assessed when the quantity of pheromone emitted is above a given threshold $\epsilon_s$. This threshold depicts the minimal quantity of pheromone emitted by a single insect over the area of a numerical mesh cell. According to this presence criteria, the subdomain where the insects are located given a source term is the superlevel set

$$ (16) \qquad \Omega_{s/s}(t; s, \epsilon_s) = \{(x, y) \in \Omega \mid s(x, y, t) > \epsilon_s\} $$

This implies that the superlevel set $\Omega_{s/s}(t; s_a, \epsilon_s)$ (set in blue dashed lines on Fig. 4) is the subdomain where the insects are predicted to be present. On the other hand, the superlevel set $\Omega_{s/s}(t; s_t, \epsilon_s)$ (set in red dashed lines on Fig. 4) is the subdomain where the insects are truly present.
where the insects are wrongly predicted to be present (false presence prediction, set in red dotted lines and blue dashed lines on Fig. 4) and $\Omega_{fa}(t; \epsilon_s) = \overline{\Omega_{s|s}(t; s_a, \epsilon_s)} \cap \Omega_{s|s}(t; s_t, \epsilon_s)$ the subdominant where the insects are wrongly predicted to be absent (false absence prediction, set in red dashed lines and blue dotted lines on Fig. 4). Noting $|\Omega|$ the area of the set $\Omega$, we then define $\rho_{fa}(t)$ the false absence prediction rate (or miss rate, i.e. the proportion of insect true presence that is predicted as absence):

340

$$ (17) \qquad \rho_{fa}(t) = \frac{|\Omega_{fa}(t)|}{|\Omega_{s|s}(t; s_t, \epsilon_s)|} = \frac{|\Omega_{fa}(t)|}{|\Omega_{fa}(t) \cup \Omega_{tp}(t)|} $$

and $\rho_{fp}(t)$ the false presence prediction rate (or false discovery rate, i.e. the proportion of predicted presence that is incorrect):

345

$$ (18) \qquad \rho_{fp}(t) = \frac{|\Omega_{fp}(t)|}{|\Omega_{s|s}(t; s_a, \epsilon_s)|} = \frac{|\Omega_{fp}(t)|}{|\Omega_{fp}(t) \cup \Omega_{tp}(t)|} $$

In a context where we aim to prevent and control the infestation, this miss rate must be kept the smallest possible, or even 0. On the other hand, in a context of precision agriculture, the false discovery rate will lead to useless application of pesticide. Thus, the objectif is to have the smallest possible false discovery rate while keeping the miss rate close to 0. In the following, we use the threshold value $\epsilon_s = 0.065$.

350

### 2.4.3. One-sensor adjoint state-based criteria.
The purpose of this criteria is to assess how the regularization strategies influence indirectly the contribution and importance of the sensors and their data, and if some strategies require less sensors to achieve the inference with a similar accuracy on $s_a$.

Let $(s_k)_{k \in \mathbb{N}}$ denotes the sequence of source term given the trajectory taken by the optimization algorithm in the control space. The gradient descent algorithm iterations can then be expressed as $\forall k > 0$:

355

$$ (19) \qquad s_{k+1} = s_k - \eta_k \nabla_{s} j(s_k) $$

with $(\eta_k)_{k \in \mathbb{N}}$ the sequence of step sizes. From this expression of one iteration, the optimum $s_a$ can then be expressed as function of the sequence $(s_k)_{k \in \mathbb{N}}$:

$$ (20) \qquad s_a = s_0 - \sum_k \eta_k \nabla_{s} j(s_k) $$

Using the expression (7) of $\nabla_{s} j$ and the splitting of the adjoint state $c^*$ into the one-sensor adjoint state $c_i^*$ as in Eq. (9), the optimum $s_a$ can then be expressed as:

360

$$ \begin{aligned} s_a &= s_0 - \sum_k \eta_k (\nabla_{s} j_{reg}(s_k) - c^*(s_k)) \\ &= s_0 - \sum_k (\eta_k \nabla_{s} j_{reg}(s_k)) + \sum_i \sum_k (\eta_k c_i^*(s_k)) \\ &= s_0 - \sum_k (\eta_k \nabla_{s} j_{reg}(s_k)) + \sum_i c_{i,tot}^* \end{aligned} $$

where

$$ (21) \qquad c_{i,tot}^*(x, y, t) = \sum_k \eta_k c_i^*(s_k; x, y, t) $$

This illustrates that $c_{i,tot}^*$ enables to quantify the contribution of the $i^{th}$ sensor to the computation of the optimal $s_a$. The comparison of the regularization strategies with regard of the sensors contribution will be done by looking at $||c_{i,tot}^*||_{L^\infty(\Omega \times [0, T])}$, the maximal contribution of
<table><thead><tr><th>prediction<br>truth</th><th>presence<br>Ω<sub>sls</sub>(s<sub>a</sub>)</th><th>absence<br>Ω&#773;<sub>sls</sub>(s<sub>a</sub>)</th></tr></thead><tbody><tr><th>presence<br>Ω<sub>sls</sub>(s<sub>t</sub>)</th><td>true presence prediction<br>Ω<sub>tp</sub> = Ω&#773;<sub>sls</sub>(s<sub>a</sub>) ∩ Ω&#773;<sub>sls</sub>(s<sub>t</sub>)</td><td>false presence prediction<br>Ω<sub>fp</sub> = Ω<sub>sls</sub>(s<sub>a</sub>) ∩ Ω&#773;<sub>sls</sub>(s<sub>t</sub>)</td></tr><tr><th>absence<br>Ω&#773;<sub>sls</sub>(s<sub>t</sub>)</th><td>false absence prediction<br>Ω<sub>fa</sub> = Ω&#773;<sub>sls</sub>(s<sub>a</sub>) ∩ Ω<sub>sls</sub>(s<sub>t</sub>)</td><td>true absence prediction<br>Ω<sub>ta</sub> = Ω&#773;<sub>sls</sub>(s<sub>a</sub>) ∩ Ω&#773;<sub>sls</sub>(s<sub>t</sub>)</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 4 - Illustration of the predicted and true superlevel sets and the false absence and presence predictions</figcaption></figure>

These two sets will be compared using dissimilarity criteria in a same way as the Jaccard similarity coefficient or the confusion matrix in supervised statistical learning, see e.g. (Stehman, 1997). The subdomain $Ω_{fp}(t; ϵ_s) = Ω_{sls}(t; s_a, ϵ_s) ∩ Ω&#773;<sub>sls</sub>(t; s_t, ϵ_s)$ corresponds to the subdomain where the insects are wrongly predicted to be present (false presence prediction, set in red dotted lines and blue dashed lines on Fig. 4) and $Ω_{fa}(t; ϵ_s) = Ω&#773;<sub>sls</sub>(t; s_a, ϵ_s) ∩ Ω<sub>sls</sub>(t; s_t, ϵ_s)$ the subdomain where the insects are wrongly predicted to be absent (false absence prediction, set in red dashed lines and blue dotted lines on Fig. 4). Noting $|Ω|$ the area of the set $Ω$, we then define $ρ_{fa}(t)$ the false absence prediction rate (or miss rate, i.e. the proportion of insect true presence that is predicted as absence):

$$ (17) \qquad ρ_{fa}(t) = \frac{|Ω_{fa}(t)|}{|Ω_{sls}(t; s_t, ϵ_s)|} = \frac{|Ω_{fa}(t)|}{|Ω_{fa}(t) \cup Ω_{tp}(t)|} $$

and $ρ_{fp}(t)$ the false presence prediction rate (or false discovery rate, i.e. the proportion of predicted presence that is incorrect):

$$ (18) \qquad ρ_{fp}(t) = \frac{|Ω_{fp}(t)|}{|Ω_{sls}(t; s_a, ϵ_s)|} = \frac{|Ω_{fp}(t)|}{|Ω_{fp}(t) \cup Ω_{tp}(t)|} $$

In a context where we aim to prevent and control the infestation, this miss rate must be kept the smallest possible, or even 0. On the other hand, in a context of precision agriculture, the false discovery rate will lead to useless application of pesticide. Thus, the objectif is to have the smallest possible false discovery rate while keeping the miss rate close to 0. In the following, we use the threshold value $ϵ_s = 0.065$.

### 2.4.3. One-sensor adjoint state-based criteria.
The purpose of this criteria is to assess how the regularization strategies influence indirectly the contribution and importance of the sensors and their data, and if some strategies require less sensors to achieve the inference with a similar accuracy on $s_a$.

Let $(s_k)_{k \in N}$ denotes the sequence of source terms given the trajectory taken by the optimization algorithm in the control space. The gradient descent algorithm iterations can then be expressed as $\forall k > 0$:

$$ (19) \qquad s_{k+1} = s_k - η_k \nabla_s j(s_k) $$

with $(η_k)_{k \in N}$ the sequence of step sizes. From this expression of one iteration, the optimum $s_a$ can then be expressed as function of the sequence $(s_k)_{k \in N}$:

$$ (20) \qquad s_a = s_0 - \sum_k η_k \nabla_s j(s_k) $$
the $i^{th}$ sensor, and the sensors will be ordered by decreasing $||c_{i,tot}^*||_{L^\infty(\Omega \times [0, T])}$, hereinafter referred to as the contribution order. This representation allows to compare the contribution of the sensors with each other and to decipher if the maximal total contribution drops from a certain rank in the contribution order, thus indicating that the sensors located upwards in the ranking may have lower, or even negligible, impact on the estimation of $s_a$. Let us note that the quantity $c_{i,tot}^*$ can be directly computed using a trajectory-integrated one-sensor adjoint model in eq. (69) and does not require to compute $c_i^*(s_k)$ for each $k$ separately, see more details in Sec. B.3.1.

## 3. Results

The core methodological improvement of this study is the Biology-Informed Data Assimilation (BI-DA) framework. It is the reason why we present it in the results and not in the material and methods section. This framework is based on the classical data assimilation optimization problem (3) but with in addition regularization terms that are informed with biological knowledge of the insect. These BI-DA framework is tested in the two test cases defined in Sec. 2.3.2, i.e. the Toy case, characterized by small spatio-temporal scales, homogeneous landscape, time-varying but space-uniform wind field and constant emission pattern of the emitting insects, and the FAW case taking into account a real agricultural landscape, observed wind field and complex emission pattern mimicking the behavior of the insect pest targetted in this study: the FAW.

We first present different biology-informed regularization terms in Sec. 3.1. In Sec. 3.2, the pheromone propagation is illustrated in the Toy case, before computing a classical VDA without regularization in Sec. 3.3, to obtain an optimal pheromone emission $s_a$ of reference. In Sec. 3.4, the BI-DA is extensively benchmarked and compared to regularization-free VDA in the Toy case. Finally, BI-DA is used to track emission source in the FAW case. In Sec. 3.5, the pheromone propagation is illustrated in this realistic setting. In Sec. 3.6, a classical VDA used as reference is performed without regularization, raising the need for a better sensor placement. In Sec. 3.7, BI-DA is performed to infer the localization of the insects.

### 3.1. Biology-Informed DA: adding biological priors in pheromone emission inference

As mentioned in Sec. 2.2, regularization terms based on additional coarse estimation of the control are often considered to improve the convexity of $j$ (Kaltenbacher et al., 2008; Tarantola, 2005). In the classical VDA framework, the optimization problem (3) is strongly constrained by the pheromone propagation model (1): the pheromone distribution $c$ is explicitly obtained from the direct problem, given the source $s$. Regularization terms can be seen as additional weak constraints, i.e. constraints introduced in the minimization problem as penalties. Inference problems that rely only on optimization problems weakly-constrained by a dynamical model are already treated in VDA (Trémolet, 2007) and machine learning contexts (Raissi et al., 2019).

We specifically introduce regularization terms based on biological models, and thus term the resulting VDA problem as BI-DA. Let us assume that $s$ satisfies a biological model that can be written $M(s) = \sigma$, where $\sigma$ is a known biological prior. Then, one can consider a regularization term of the form

$$ (22) \qquad j_{reg}(s) = \|M(s) - \sigma\|_{L^2}^2 $$

The regularization term is then designed to minimize the residual of the model $M(s) = \sigma$.

The model's operator $M$ can be any kind of operator, including partial differential operator. For instance, recall that the emission source can always be written $s = qp$ as in Eq. (2), and thus, the model $M(s) = \sigma$ can be resulting from a population dynamics model on the insect density $p$. In the VDA minimization problem (3), we work in a differentiable optimization framework.
Using the expression (7) of $\nabla_{sj}$ and the splitting of the adjoint state $c^*$ into the one-sensor adjoint state $c_i^*$ as in Eq. (9), the optimum $s_a$ can then be expressed as:

$$ 
\begin{aligned}
s_a &= s_0 - \sum_k \eta_k (\nabla_{sj} \mathrm{reg}(s_k) - c^*(s_k)) \\
&= s_0 - \sum_k (\eta_k \nabla_{sj} \mathrm{reg}(s_k)) + \sum_i \sum_k (\eta_k c_i^*(s_k)) \\
&= s_0 - \sum_k (\eta_k \nabla_{sj} \mathrm{reg}(s_k)) + \sum_i c_{i,tot}^*
\end{aligned}
$$

where

$$ (21) \qquad c_{i,tot}^*(x, y, t) = \sum_k \eta_k c_i^*(s_k; x, y, t) $$

This illustrates that $c_{i,tot}^*$ enables to quantify the contribution of the $i^{th}$ sensor to the computation of the optimal $s_a$. The comparison of the regularization strategies with regard of the sensors contribution will be done by looking at $||c_{i,tot}^*||_{L^\infty(\Omega \times [0, T])}$, the maximal contribution of the $i^{th}$ sensor, and the sensors will be ordered by decreasing $||c_{i,tot}^*||_{L^\infty(\Omega \times [0, T])}$, hereinafter referred to as the contribution order. This representation allows to compare the contribution of the sensors with each other and to decipher if the maximal total contribution drops from a certain rank in the contribution order, thus indicating that the sensors located upwards in the ranking may have lower, or even negligible, impact on the estimation of $s_a$. Let us note that the quantity $c_{i,tot}^*$ can be directly computed using a trajectory-integrated one-sensor adjoint model in eq. (69) and does not require to compute $c_i^*(s_k)$ for each $k$ separately, see more details in Sec. A.3.1.

## 3. Results

The core methodological improvement of this study is the Biology-Informed Data Assimilation (BI-DA) framework. It is the reason why we present it in the results and not in the material and methods section. This framework is based on the classical data assimilation optimization problem (3) but with in addition regularization terms that are informed with biological knowledge of the insect. These BI-DA framework is tested in the two test cases defined in Sec. 2.3.2, i.e. the Toy case, characterized by small spatio-temporal scales, homogeneous landscape, time-varying but space-uniform wind field and constant emission pattern of the emitting insects, and the FAW case taking into account a real agricultural landscape, observed wind field and complex emission pattern mimicking the behavior of the insect pest targetted in this study: the FAW.

We first present different biology-informed regularization terms in Sec. 3.1. In Sec. 3.2, the pheromone propagation is illustrated in the Toy case, before computing a classical VDA without regularization in Sec. 3.3, to obtain an optimal pheromone emission $s_a$ of reference. In Sec. 3.4, the BI-DA is extensively benchmarked and compared to regularization-free VDA in the Toy case. Finally, BI-DA is used to track emission source in the FAW case. In Sec. 3.5, the pheromone propagation is illustrated in this realistic setting. In Sec. 3.6, a classical VDA used as reference is performed without regularization, raising the need for a better sensor placement. In Sec. 3.7, BI-DA is performed to infer the localization of the insects.

### 3.1. Biology-Informed DA: adding biological priors in pheromone emission inference

As mentioned in Sec. 2.2, regularization terms based on additional coarse estimation of the control are often considered to improve the convexity of $j$ (Kaltenbacher et al., 2008; Tarantola, 2005). In the classical VDA framework, the optimization problem (3) is strongly constrained by the pheromone propagation model (1): the pheromone distribution $c$ is explicitly obtained from
Thus, the loss function $j(s)$ must be differentiable, which implies that $j_{reg}(s)$ and the underlying operator $s \mapsto M(s)$ have to be differentiable as well. However, this framework can be extended by working in a non-smooth and subdifferentiable optimization framework (Nesterov, 2018), e.g. by using the proximal-gradient algorithm (Parikh et al., 2014). This allows us to consider subdifferentiable $j_{reg}(s)$ and $M$, such as a LASSO regularization term, i.e. of the form $j_{reg}(s) = \|M(s) - \sigma\|_L^1$ or alike, see more details in Sec. 3.1.2. Recall from Eq. (5) that a weighted sums of multiple regularization terms can be considered in order to take into account multiple biological models.

We now introduce several biology-informed regularization terms that can be used to regularize the optimization problem (3). First, classical regularization terms will be introduced, such as Tikhonov (a.k.a. ridge), LASSO and group-LASSO regularization. For these regularization terms, we will pay a specific attention to the biological information we can add through these terms. Then, a population dynamic-informed $j_{reg}$ will be presented, in order to explicitly model biological mechanisms.

**3.1.1. Biology-based Tikhonov regularization term.** The most classical $j_{reg}$ is the Tikhonov regularization term (Kaltenbacher et al., 2008):

$$ (23) \qquad j_{reg, \tau}(s) = \|s - s_b\|_{C^{-1}}^2 $$

with $s_b$ a background value of $s$ and $C$ the background error covariance operator. Classically, $s_b$ is given by a prior estimate of the true state, e.g. based on additional data that are not accurate enough and/or too sparse or based on a rough physical estimate, see e.g. (Bocquet, 2014). From a biological point of view, a Biology-Informed Tikhonov $j_{reg}$ ($M(s) = s$ and $\sigma = s_b$) enables to take into account prior knowledge on the insects life habits, for instance their preferred habitat and insect life cycle. In this case, the term $s_b$ describes the expected spatial distribution and density of the insects with respect to the land use or the pheromone emission pattern.

**3.1.2. LASSO-type regularization term and corresponding biological input.** LASSO-type regularization terms are often used to ensure sparsity and perform variable selection, see (Tibshirani, 1996). In the VDA experiments based only on the data, see Secs. 3.3 and 3.6, some spurious sources are detected, including near several sensors even though they are far away from the target source. This residual detection, particularly visible near final time (see Fig. 6), can be due to numerical error, to noise or to diffusive phenomena in the adjoint model. When computing the superlevel set in eq. (16) for locating insects from the inferred quantity of pheromone emitted, these spurious pheromone sources can lead to false presence prediction.

**LASSO regularization.** In our context, we assume that the data are collected at the beginning of the infestation, so that we can expect that only a few insects are present. Hence, in order to select important local pheromones sources and to mitigate false presence predictions, a LASSO penalty is considered:

$$ (24) \qquad j_{reg, LASSO}(s) = \|s\|_1 $$

Recall that the optimization problem (3) with this $j_{reg}$ is no longer differentiable. A proximal gradient algorithm, namely the ISTA method, will be used to solve it, see more details in Annex B.3.3. Note that the ISTA is a Iterative Soft-Threshold Algorithm, meaning in particular that at each optimization iterations it thresholds any value under $\alpha_{reg, LASSO}\eta$, where $\alpha_{reg, LASSO}$ is the regularization weight of $j_{reg, LASSO}$ and $\eta$ is the descent step of the gradient-descent algorithm. Hence,
the direct problem, given the source s. Regularization terms can be seen as additional weak con- straints, i.e. constraints introduced in the minimization problem as penalties. Inference problems that rely only on optimization problems weakly-constrained by a dynamical model are already treated in VDA (Trémolet, 2007) and machine learning contexts (Raissi et al., 2019).

We specifically introduce regularization terms based on biological models, and thus term the resulting VDA problem as BI-DA. Let us assume that s satisfies a biological model that can be written $M(s) = \sigma$, where $\sigma$ is a known biological prior. Then, one can consider a regularization term of the form

$$ (22) \qquad j_{\text{reg}}(s) = \|M(s) - \sigma\|_{L^2}^2 $$

The regularization term is then designed to minimize the residual of the model $M(s) = \sigma$.

The model's operator $M$ can be any kind of operator, including partial differential operator. For instance, recall that the emission source can always be written $s = qp$ as in Eq. (2), and thus, the model $M(s) = \sigma$ can be resulting from a population dynamics model on the insect density $p$. In the VDA minimization problem (3), we work in a differentiable optimization framework. Thus, the loss function $j(s)$ must be differentiable, which implies that $j_{\text{reg}}(s)$ and the underlying operator $s \mapsto M(s)$ have to be differentiable as well. However, this framework can be extended by working in a non-smooth and subdifferentiable optimization framework (Nesterov, 2018), e.g. by using the proximal-gradient algorithm (Parikh et al., 2014). This allows us to consider subdifferentiable $j_{\text{reg}}(s)$ and $M$, such as a LASSO regularization term, i.e. of the form $j_{\text{reg}}(s) = \|M(s) - \sigma\|_{L^1}$ or alike, see more details in Sec. 3.1.2. Recall from Eq. (5) that a weighted sums of multiple regularization terms can be considered in order to take into account multiple biological models.

We now introduce several biology-informed regularization terms that can be used to regular- ize the optimization problem (3). First, classical regularization terms will be introduced, such as Tikhonov (a.k.a. ridge), LASSO and group-LASSO regularization. For these regularization terms, we will pay a specific attention to the biological information we can add through these terms. Then, a population dynamic-informed $j_{\text{reg}}$ will be presented, in order to explicitly model biological mechanisms.

### 3.1.1. Biology-based Tikhonov regularization term.
The most classical $j_{\text{reg}}$ is the Tikhonov regular- ization term (Kaltenbacher et al., 2008):

$$ (23) \qquad j_{\text{reg}, \tau}(s) = \|s - s_b\|_{C^{-1}}^2 $$

with $s_b$ a background value of s and C the background error covariance operator. Classically, $s_b$ is given by a prior estimate of the true state, e.g. based on additional data that are not accurate enough and/or too sparse or based on a rough physical estimate, see e.g. (Bocquet, 2014). From a biological point of view, a Biology-Informed Tikhonov $j_{\text{reg}}$ ($M(s) = s$ and $\sigma = s_b$) enables to take into account prior knowledge on the insects life habits, for instance their preferred habitat and insect life cycle. In this case, the term $s_b$ describes the expected spatial distribution and density of the insects with respect to the land use or the pheromone emission pattern.

### 3.1.2. LASSO-type regularization term and corresponding biological input.
LASSO-type regulariza- tion terms are often used to ensure sparsity and perform variable selection, see (Tibshirani, 1996). In the VDA experiments based only on the data, see Secs. 3.3 and 3.6, some spurious sources are detected, including near several sensors even though they are far away from the target source. This residual detection, particularly visible near final time (see Fig. 6), can be due to numerical er- ror, to noise or to diffusive phenomena in the adjoint model. When computing the superlevel set
assuming that we have a prior knowledge of a detection threshold under which pheromone emission can be discarded, a magnitude of $\alpha_{reg,LASSO}$ can be computed from $\eta$ in order to shrink the estimate of $s_a$ under this threshold at each optimization iterations.

**Group-LASSO regularization.** It is common that the insects do not move during the pheromone emission time period. It is e.g. the case for FAW, see Sec. 2.3.2. Therefore, the total quantity of pheromone emitted during this emission period can not be lower than the total amount of pheromone emitted by one insect on this time window. Conversely, if the total amount of emitted pheromone is lower, it implies that these pheromone sources are spurious and can be discarded. The following group-LASSO term is adapted to eliminate spatial positions with lowest accumulated pheromone emission:

$$ (25) \qquad j_{reg,gLASSO}(s) = \|s\|_{L^1(\Omega; L^2(0,T))} = \int_{\Omega} \left( \int_{0}^{T} s(t, x, y)^2 dt \right)^{1/2} dxdy $$

Similarly to the LASSO regularization, a proximal gradient algorithm will be used and more details are provided in B.3.3. For the group-LASSO regularization, this algorithm thresholds at each optimization iterations $s(t, x, y)$ for all $t$ if $\|s(\cdot, x, y)\|_{L^2(0,T)}$ is under $\alpha_{reg,gLASSO}\eta$, where $\alpha_{reg,gLASSO}$ is the regularization weight of $j_{reg,gLASSO}$ and $\eta$ is the descent step of the gradient-descent algorithm.

**3.1.3. Population dynamic-informed regularization term.** Population dynamics models that model the spatio-temporal evolution of the density of insects $p$ can also be used to derive biological models for the quantity of pheromone emitted $s$. These population dynamics models can be written in a generic way as follows:

$$ (26) \qquad \partial_t p + \nabla \left( \sum_i F_i(p) \right) - \sum_i R_i(p) = 0 $$

with $F_i$ and $R_i$ are resp. flux and reaction terms that model the behavior of the insect specie of interest. Recall that $s(x, y, t) = p(x, y, t)q(t) \Leftrightarrow p(x, y, t) = \frac{s(x, y, t)}{q(t)}$ with $q(t) > 0$ the quantity of pheromone emitted per insect, which is assumed to be a biological knowledge. A population dynamics PDE equation of the form (26) can be expressed as the following PDE equation, in which $s$ is the state variable:

$$ \frac{1}{q} \partial_t s - \frac{\partial_t q}{q^2} s + \nabla \left( \sum_i \mathcal{F}_i(s) \right) + \sum_i \mathcal{R}_i(s) = 0 $$

with $F_i(p = \frac{s}{q}) = \mathcal{F}_i(s)$ and $R_i(p = \frac{s}{q}) = \mathcal{R}_i(s)$. From this PDE equation, the population dynamics-informed regularization term is the following:

$$ (27) \qquad j_{reg,PD}(s) = \alpha_{reg,PD} \|\partial_t s - \frac{\partial_t q}{q} s + q \nabla \left( \sum_i \mathcal{F}_i(s) \right) - q \sum_i \mathcal{R}_i(s)\|_{L^2} $$

Several potential terms $F_i$ and $R_i$ (and the associated $\mathcal{F}_i$ and $\mathcal{R}_i$) are possible and have different biological meaning. The flux terms $F_i$ can be used to model the motion of the insects in a random direction or in a specific direction (e.g. toward the preferred habitat). The reaction terms can be used to model birth-death processes or the fact that certain insects stop emitting pheromones. An extensive list of additional possible terms is added in the supplementary material Sec. B.2, including advection, diffusion, quorum sensing, Keller-Segel term to model pheromone sensing, and additional reaction terms such as linear birth-death process, logistic growth, Allee effect, or non-local terms. Furthermore, the corresponding gradients are explicitly computed in the supplementary Sec. B.3.2, so that the reader has all the necessary elements to
in eq. (16) for locating insects from the inferred quantity of pheromone emitted, these spurious pheromone sources can lead to false presence prediction.

**LASSO regularization.** In our context, we assume that the data are collected at the beginning of the infestation, so that we can expect that only a few insects are present. Hence, in order to select important local pheromones sources and to mitigate false presence predictions, a LASSO penalty is considered:

$$ (24) \qquad j_{reg,LASSO}(s) = \|s\|_1 $$

Recall that the optimization problem (3) with this $j_{reg}$ is no longer differentiable. A proximal gradient algorithm, namely the ISTA method, will be used to solve it, see more details in Annex A.3.3. Note that the ISTA is a Iterative Soft-Threshold Algorithm, meaning in particular that at each optimization iterations it thresholds any value under $\alpha_{reg,LASSO}\eta$, where $\alpha_{reg,LASSO}$ is the regularization weight of $j_{reg,LASSO}$ and $\eta$ is the descent step of the gradient-descent algorithm. Hence, assuming that we have a prior knowledge of a detection threshold under which pheromone emission can be discarded, a magnitude of $\alpha_{reg,LASSO}$ can be computed from $\eta$ in order to shrink the estimate of $s_a$ under this threshold at each optimization iterations.

**Group-LASSO regularization.** It is common that the insects do not move during the pheromone emission time period. It is e.g. the case for FAW, see Sec. 2.3.2. Therefore, the total quantity of pheromone emitted during this emission period can not be lower than the total amount of pheromone emitted by one insect on this time window. Conversely, if the total amount of emitted pheromone is lower, it implies that these pheromone sources are spurious and can be discarded. The following group-LASSO term is adapted to eliminate spatial positions with lowest accumulated pheromone emission:

$$ (25) \qquad j_{reg,gLASSO}(s) = \|s\|_{L^1(\Omega; L^2(0,T))} = \int_{\Omega} \left( \int_{0}^{T} s(t, x, y)^2 dt \right)^{1/2} dxdy $$

Similarly to the LASSO regularization, a proximal gradient algorithm will be used and more details are provided in A.3.3. For the group-LASSO regularization, this algorithm thresholds at each optimization iterations $s(t, x, y)$ for all $t$ if $\|s(\cdot, x, y)\|_{L^2(0,T)}$ is under $\alpha_{reg,gLASSO}\eta$, where $\alpha_{reg,gLASSO}$ is the regularization weight of $j_{reg,gLASSO}$ and $\eta$ is the descent step of the gradient-descent algorithm.

### 3.1.3. Population dynamic-informed regularization term.
Population dynamics models that model the spatio-temporal evolution of the density of insects $p$ can also be used to derive biological models for the quantity of pheromone emitted $s$. These population dynamics models can be written in a generic way as follows:

$$ (26) \qquad \partial_t p + \nabla \left( \sum_i F_i(p) \right) - \sum_i R_i(p) = 0 $$

with $F_i$ and $R_i$ are resp. flux and reaction terms that model the behavior of the insect specie of interest. Recall that $s(x, y, t) = p(x, y, t)q(t) \Leftrightarrow p(x, y, t) = \frac{s(x, y, t)}{q(t)}$ with $q(t) > 0$ the quantity of pheromone emitted per insect, which is assumed to be a biological knowledge. A population dynamics PDE equation of the form (26) can be expressed as the following PDE equation, in which $s$ is the state variable:

$$ \frac{1}{q} \partial_t s - \frac{\partial_t q}{q^2} s + \nabla \left( \sum_i \mathcal{F}_i(s) \right) + \sum_i \mathcal{R}_i(s) = 0 $$
implement the population dynamic-informed penalty fitting with the behavior of other insects of interest.

For instance, the females of the model insect that we consider, FAW, start emitting pheromones once found an emission site. During the emission time lapse, the insect does not move in order to efficiently attract a mating partner and stops emitting when it finds a mating partner and the mating starts. Furthermore, a proportion of female encounter a male at a given time following a mating rate function, see Sec. 2.3.2. Thus, the dynamics of the population of the insects emitting pheromone is modelled using only a simple linear birth-death process:

$$ (28) \qquad R_{bd}(p; x, y, t) = -\gamma(x, y, t)p(x, y, t) \Leftrightarrow R_{bd}(s; x, y, t) = -\frac{\gamma(x, y, t)}{q(t)}s(x, y, t) $$

with $\gamma$ [$s^{-1}$] the birth (if $\gamma > 0$) or death (if $\gamma < 0$) rate. In general manners, this term models a linear birth-death process where a proportion $\gamma(t)$ of the population dies or reproduces at each time $t$.

### 3.1.4. Choices of hyperparameter and hot start.
As mentioned previously, we will use the different regularization terms, which includes using them simultaneously. In this regard, recall that the overall regularization term can be expressed as the weighted sum of BI regularization terms (5) and $\mathcal{L}$ is the list of the regularization terms that we take into account. In the present study, we set $\mathcal{L} = \{ 'T', 'LASSO', 'gLASSO', 'PD' \}$ where T stands for Tikhonov, PD for population dynamics and gLASSO for group-LASSO. In the following, several combination of $j_{reg,i}$ will be tested and analysed, starting from only one $j_{reg,i}$ ($\alpha_j = 0 \forall j \neq i$) and up to three $j_{reg,i}$. The composite cases 'T + LASSO', 'T + PD' and 'LASSO + PD' designate the cases when, additionally to $\alpha_{gLASSO} = 0$, we have respectively $\alpha_{PD} = 0$, $\alpha_{LASSO} = 0$ and $\alpha_{T} = 0$. Finally, the case 'T + LASSO + PD' is when only $\alpha_{gLASSO} = 0$ and is designated as 'all reg' in the Toy case. Similarly, the case 'T + gLASSO + PD' is when only $\alpha_{LASSO} = 0$ and is designated as 'all reg' in the FAW case. Moreover, in the following, 'no reg' refers to the case without regularization terms ($\alpha_i = 0 \forall i \in \mathcal{L}$).

The values of $\alpha_i$ are chosen such that they are the largest possible while satisfying the Morozov's discrepancy principle (see Sec. 2.2.4) in a reasonable computation time. In the following, in order to facilitate the solving of optimization problem, the regularization terms are progressively added. This process is also known as a hot start. In a first time, we optimize without regularization. Then, we add one by one the regularization terms initiating the optimization with the optimum obtained previously. Similarly, in order to reach the highest values possible, the weight coefficients $\alpha_i$ are increased step by step, restarting each time with hot start and repeating until the Morozov's discrepancy principle is no longer satisfied.

## 3.2. Numerical illustration of pheromone propagation in the Toy case

We first present numerical results for the pheromone propagation model (1) in the Toy case, as set up in Sec. 2.3.2, in order to illustrate the pheromone plumes computed by the model (see Fig. 5). As expected, the convection is dominant: the pheromones mainly propagate in the direction of the wind, with little additional diffusive effect. The pheromone crosses the domain within a couple of seconds during which a front propagates until reaching a quasi stationary plume. As expected, time variations of the wind field, that is north-east at the beginning and north at the end, induces deformation of the plume.
with $F_i(p = \frac{s}{q}) = F_i(s)$ and $R_i(p = \frac{s}{q}) = R_i(s)$. From this PDE equation, the population dynamics-informed regularization term is the following:

$$ (27) \qquad j_{reg,PD}(s) = \alpha_{reg,PD} \|\partial_t s - \frac{\partial_t q}{q} s + q \nabla \left( \sum_i F_i(s) \right) - q \sum_i R_i(s)\|_{L^2}^2 $$

We note that this residual uses the PDE equation expressed in strong form. We take into consideration no boundary effects, but if needed, boundary conditions could be added to the minimization problem. For regular $F_i$ and $R_i$ operators, this residual lies in $L^2(0, T; \Omega)$.

Several potential terms $F_i$ and $R_i$ (and the associated $F_i$ and $R_i$) are possible and have different biological meaning. The flux terms $F_i$ can be used to model the motion of the insects in a random direction or in a specific direction (e.g. toward the preferred habitat). The reaction terms can be used to model birth-death processes or the fact that certain insects stop emitting pheromones. An extensive list of additional possible terms is added in the supplementary material Sec. A.2, including advection, diffusion, quorum sensing, Keller-Segel term to model pheromone sensing, and additional reaction terms such as linear birth-death process, logistic growth, Allee effect, or non-local terms. Furthermore, the corresponding gradients are explicitly computed in the supplementary Sec. A.3.2, so that the reader has all the necessary elements to implement the population dynamic-informed penalty fitting with the behavior of other insects of interest.

For instance, the females of the model insect that we consider, FAW, start emitting pheromones once found an emission site. During the emission time lapse, the insect does not move in order to efficiently attract a mating partner and stops emitting when it finds a mating partner and the mating starts. Furthermore, a proportion of female encounter a male at a given time following a mating rate function, see Sec. 2.3.2. Thus, the dynamics of the population of the insects emitting pheromone is modelled using only a simple linear birth-death process:

$$ (28) \qquad R_{bd}(p; x, y, t) = -\gamma(x, y, t)p(x, y, t) \Leftrightarrow R_{bd}(s; x, y, t) = -\frac{\gamma(x, y, t)}{q(t)}s(x, y, t) $$

with $\gamma$ [$s^{-1}$] the birth (if $\gamma > 0$) or death (if $\gamma < 0$) rate. In general manners, this term models a linear birth-death process where a proportion $\gamma(t)$ of the population dies or reproduces at each time $t$.

3.1.4. *Choices of hyperparameter and hot start.* As mentioned previously, we will use the different regularization terms, which includes using them simultaneously. In this regard, recall that the overall regularization term can be expressed as the weighted sum of BI regularization terms (5) and $\mathcal{L}$ is the list of the regularization terms that we take into account. In the present study, we set $\mathcal{L} = \{ 'T', 'LASSO', 'gLASSO', 'PD' \}$ where $T$ stands for Tikhonov, $PD$ for population dynamics and $gLASSO$ for group-LASSO. In the following, several combination of $j_{reg,i}$ will be tested and analysed, starting from only one $j_{reg,i}$ ($\alpha_j = 0 \forall j \neq i$) and up to three $j_{reg,i}$. The composite cases 'T + LASSO', 'T + PD' and 'LASSO + PD' designate the cases when, additionally to $\alpha_{gLASSO} = 0$, we have respectively $\alpha_{PD} = 0$, $\alpha_{LASSO} = 0$ and $\alpha_{T} = 0$. Finally, the case 'T + LASSO + PD' is when only $\alpha_{gLASSO} = 0$ and is designated as 'all reg' in the Toy case. Similarly, the case 'T + gLASS + PD' is when only $\alpha_{LASSO} = 0$ and is designated as 'all reg' in the FAW case. Moreover, in the following, 'no reg' refers to the case without regularization terms ($\alpha_i = 0 \forall i \in \mathcal{L}$).

The values of $\alpha_i$ are chosen such that they are the largest possible while satisfying the Morozov's discrepancy principle (see Sec. 2.2.4) in a reasonable computation time. In the following, in order to facilitate the solving of optimization problem, the regularization terms are progressively
<figure><img src="image_8.png" /><figcaption>**Figure 5** – Map of the concentration in pheromone $c(s_t)$ ($\text{pg.m}^{-2}$) computed from the depth-integrated CTM model (1) with the reference source term $s_t$ for the Toy case presented in Sec. 2.3.2 at several times; in addition iso-lines of $c$ are indicated with black dashed lines, and the reference source $s_t$ support is displayed in white at final time.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6** – Classical VDA without regularization for the Toy case. The optimal quantity of pheromone emitted $s_a$ is obtained by inference without regularization term for the Toy case at $T = 0\text{s}$ (left panel) and $T = 20\text{s}$ (right panel). The support of the target quantity of emitted pheromone $s_t$ (see Fig. 2) is delimited by the dashed lines. Some spurious detection spots can be seen, as sources (light dot) or sinks (dark blue dots) in the Toy case at $T = 20\text{s}$.</figcaption></figure>
added. This process is also known as a hot start. As a first step, we optimize without regular-
ization. Then, we add one by one the regularization terms initiating the optimization with the
optimum obtained previously. Similarly, in order to reach the highest values possible, the weight
coefficients αᵢ are increased step by step, restarting each time with hot start and repeating until
the Morozov's discrepancy principle is no longer satisfied.

3.2. Numerical illustration of pheromone propagation in the Toy case

<figure><img src="image_8.png" /><figcaption>**Figure 5** – Map of the concentration in pheromone $c(s_t)$ ($\text{pg.m}^{-2}$) computed from the depth-integrated CTM model (1) with the reference source term $s_t$ for the Toy case presented in Sec. 2.3.2 at several times; in addition iso-lines of $c$ are indicated with black dashed lines, and the reference source $s_t$ support is displayed in white at final time.</figcaption></figure>

We first present numerical results for the pheromone propagation model (1) in the Toy case,
as set up in Sec. 2.3.2, in order to illustrate the pheromone plumes computed by the model
(see Fig. 5). As expected, the convection is dominant: the pheromones mainly propagate in the
direction of the wind, with little additional diffusive effect. The pheromone crosses the domain
within a couple of seconds during which a front propagates until reaching a quasi stationary
plume. As expected, time variations of the wind field, that is north-east at the beginning and
north at the end, induces deformation of the plume.

**3.3. Classical VDA without regularization used as reference in the Toy case**

We now illustrate the inference of emission sources in the Toy case without regularization
term in Fig. 6. We can observe that the optimum sₐ exhibits positive values, corresponding to
pheromone emission sources, in areas spreading out upwind of the sensor locations (see Fig. 2
for sensor locations. Some sensor locations are also visible in Fig. 6, indicated by spots). This de-
tection is referred to as the retro-plume. Indeed, the optimum sₐ is computed by descending the
gradient of the cost function (7), which is evaluated using the adjoint model (6). This equation
describes the retro-propagation with diffusion of the discrepancy between predictions and ob-
servations m(c(s)) – m[s] from final time t = T to initial time t = 0. Consequently, the gradient
descent algorithm mostly adjusts the initial prior s₀ in diffuse areas (attributable to the diffusion
term in the adjoint model) upwind (attributable to the advection term) of the sensor locations
(attributable to the source term).
### 3.3. Classical VDA without regularization used as reference in the Toy case

We now illustrate the inference of emission sources in the Toy case without regularization term in Fig. 6. We can observe that the optimum $s_a$ exhibits positive values, corresponding to pheromone emission sources, in areas spreading out upwind of the sensor locations (see Fig. 2 for sensor locations. Some sensor locations are also visible in Fig. 6, indicated by spots). This detection is referred to as the retro-plume. Indeed, the optimum $s_a$ is computed by descending the gradient of the cost function (7), which is evaluated using the adjoint model (6). This equation describes the retro-propagation with diffusion of the discrepancy between predictions and observations $m(c(s)) - m^{obs}$ from final time $t = T$ to initial time $t = 0$. Consequently, the gradient descent algorithm mostly adjusts the initial prior $s_0$ in diffuse areas (attributable to the diffusion term in the adjoint model) upwind (attributable to the advection term) of the sensor locations (attributable to the source term).

Another consequence of the retro-propagation feature can be seen in Fig. 6: if at initial time, the estimated source $s_a$ covers a large area, at final time $t = 20s$, the retro-plume is mostly confined into a small vicinity of several active sensors. Indeed, near final time, as the retro-plume is only advected for a little time interval, the propagation of the adjoint state is mostly due to diffusion.

We also note that the reconstructed source $s_a$ exhibits negative non-physical values, that can be identified by dark blue dots near some sensors (see e.g. sensor locations $X = (17, 18.5)$ and $X = (14.5, 17.0)$ at time $t = 20s$, in Fig. 6, right panel). This phenomena is not due to numerical errors: it is a feature that directly comes from the noise in the data. The VDA intends to reproduce the data, including the noise even when the noise is not consistent with physical mechanisms. Therefore, during the retro-propagation of the noisy data, the optimal source $s_a$ includes pheromone sinks, i.e. negative values, in order to approximate the noisy target, even if it is not physical.

### 3.4. Benchmark of the different regularization terms for BI-DA in the Toy case

We now aim to benchmark the different regularization terms introduced in Sec. 3.1 for BI-DA. We use the same supervised framework as in Sec. 3.3: a target source term $s_t$ is used to produce the synthetic data, and are compared to the estimated source $s_a$. Recall that $s_t = qp_{ref}^{toy}$ with $q$ the constant emission rate, and $p_{ref}^{toy}$ the pest insect density, which remains constant over time on a circular support in the Toy case, see Sec. 2.3.2. This lead to the following population dynamics model:

$$ (29) \qquad \partial_t p_{ref}^{toy}(t, x, y) = 0 $$

which corresponds to the population dynamic-informed regularization term

$$ (30) \qquad J_{reg,PD}^{toy}(s) = \|\partial_t s\|_{L^2}^2. $$

We also recall that a favorable zone for the emitting insect is defined in a square $\Omega_{square}$ that includes the support $p_{ref}^{toy}$, see Fig. 2. Moreover, from the typical density of insects and the per capita pheromone emission rate given by the literature, we compute a surface emission rate $\bar{s}_b^{toy} = 0.54 \text{pg} \cdot s^{-1} m^{-2}$ that is used as biology-informed background estimate

$$ s_b^{toy}(t, x, y) = \bar{s}_b^{toy} \mathbb{1}_{\Omega_{square}}(x, y), $$

that will be taken into consideration through the Tikhonov regularization term (23), see Sec. 3.1.1.
<figure><img src="image_9.png" /><figcaption>**Figure 6** – **Classical VDA without regularization for the Toy case.** The optimal quantity of pheromone emitted $s_a$ is obtained by inference without regularization term for the Toy case at $T = 0$s (left panel) and $T = 20$s (right panel). The support of the target quantity of emitted pheromone $s_t$ (see Fig. 2) is delimited by the dashed lines. Some spurious detection spots can be seen, as sources (light dot) or sinks (dark blue dots) in the Toy case at $T = 20$s.</figcaption></figure>

Another consequence of the retro-propagation feature can be seen in Fig. 6: if at initial time, the estimated source $s_a$ covers a large area, at final time $t = 20$s, the retro-plume is mostly confined into a small vicinity of several active sensors. Indeed, near final time, as the retro-plume is only advected for a little time interval, the propagation of the adjoint state is mostly due to diffusion.

We also note that the reconstructed source $s_a$ exhibits negative non-physical values, that can be identified by dark blue dots near some sensors (see e.g. sensor locations $X = (17, 18.5)$ and $X = (14.5, 17.0)$ at time $t = 20$s, in Fig. 6, right panel). This phenomena is not due to numerical errors: it is a feature that directly comes from the noise in the data. The VDA intends to reproduce the data, including the noise even when the noise is not consistent with physical mechanisms. Therefore, during the retro-propagation of the noisy data, the optimal source $s_a$ includes pheromone sinks, i.e. negative values, in order to approximate the noisy target, even if it is not physical. If needed, the optimization problem could be solved under positive constraints in order to avoid these drawbacks.

### 3.4. Benchmark of the different regularization terms for BI-DA in the Toy case

We now aim to benchmark the different regularization terms introduced in Sec. 3.1 for BI-DA. We use the same supervised framework as in Sec. 3.3: a target source term $s_t$ is used to produce the synthetic data, and are compared to the estimated source $s_a$. Recall that $s_t = qp_{ref}^{toy}$ with $q$ the constant emission rate, and $p_{ref}^{toy}$ the pest insect density, which remains constant over time on a circular support in the Toy case, see Sec. 2.3.2. This lead to the following population
We include in the benchmark the Tikhonov (T), LASSO, Population Dynamic (PD) regularization, and all combination of them, in order to understand how each regularization term operates and impacts the BI-DA. The optimum $s_a$ is then compared with the target source $s_t$ using the different criteria introduced in Sec. 2.4, i.e. total and support RMSE criteria (see resp. Eqs. (14) and (15)), miss and false discovery rates (see resp. Eqs. (17) and (18)), and the one-sensor adjoint state-based criteria (see Eq. (21)).

### 3.4.1. Visual comparison of estimated pheromone emissions $s_a$.
The estimated sources $s_a$ are displayed for each regularization term, including composite terms, in Fig. 7 at initial ($t = 0$s) and final times ($t = 20$s), for one single (top panel) and composite (bottom panel) regularization terms (see Sec. 3.1), together with the case without regularization ('no reg') taken as reference inference experiment.

When only a single regularization term is considered ('T', 'PD', 'LASSO' on top panel of Fig. 7), the estimate $s_a$ covers the main part of the target source $s_t$ at $t = 0.0$s, as well as when no regularization is added ('no reg'). At final time, when fewer data are available for estimation, single regularization terms do not significantly improve, or even worsen, the estimate of $s_a$ made in the reference case ('no reg'), except with the population dynamic-informed regularization term ('PD') that pulls the estimation towards constant solutions. As expected, the Tikhonov regularization ('T' on top panel of Fig. 7) tends to redistribute the emissions estimate on the square zone favorable for pests, making it visible. This includes enhancing spurious sources that are located in the favorable square but out of the circular target, and diminishing the sources in the circular target to compensate the spurious sources. At final time, the redistribution of the sources on the favorable square does not lead to a significant improvement compared to the reference case ('no reg'). Let us note that the Tikhonov regularization term also diminishes the estimate out of the favorable square. For the LASSO regularization ('LASSO' on top panel of Fig. 7), the retroplume is less diffuse, near both the initial and final time, compared to the reference ('no reg') case since small values have been thresholded by the regularization, which also corresponds to an expected feature. With the population dynamic-informed regularization ('PD' on top panel of Fig. 7), the emission estimate is slightly higher inside the target support compared to the other single regularization terms ('T' and 'LASSO') and to the reference case ('no reg'). Moreover, near the final time, the estimate $s_a$ covers the main part of the target source $s_t$.

When several regularization terms are considered (on bottom panel of Fig. 7), as long as they involve population dynamics regularization (namely 'PD + LASSO', 'PD + T' and 'all reg'), the results are visually comparable to the one obtained with the population dynamic-informed regularization ('PD'). We can note that spurious sources are estimated far from the target source at final time with the combination of population dynamics and LASSO ('PD + LASSO' on the bottom panel of Fig. 7). When the combination of population dynamics and Tikhonov is used ('PD + T' on the bottom panel of Fig. 7), we can see that the population dynamics regularization largely counterbalance the Tikhonov regularization, since few visual differences can be noticed with the case with only population dynamics ('PD'). When population dynamics, Tikhonov and LASSO terms are combined ('all reg' on the bottom panel of Fig. 7), they correctly interact and slightly improve the results at both initial and final times compared to the case of population dynamics regularization alone ('PD'), even though spurious sources in the favorable zone can be seen due to the Tikhonov regularization. On the contrary, interactions between Tikhonov and LASSO regularizations ('T + LASSO' on the bottom panel of Fig. 7) does not seem beneficial, since they hardly improve the results obtained from the reference case ('no reg').
dynamics model:

$$ (29) \qquad \partial_t p_{ref}^{toy}(t, x, y) = 0 $$

which corresponds to the population dynamic-informed regularization term

$$ (30) \qquad J_{reg,PD}^{toy}(s) = \|\partial_t s\|_{L^2}^2. $$

We also recall that a favorable zone for the emitting insect is defined in a square $\Omega_{square}$ that includes the support $p_{ref}^{toy}$, see Fig. 2. Moreover, from the typical density of insects and the per capita pheromone emission rate given by the literature, we compute a surface emission rate $\bar{s}_b^{toy} = 0.54pg \cdot s^{-1}m^{-2}$ that is used as biology-informed background estimate

$$ s_b^{toy}(t, x, y) = \bar{s}_b^{toy} \mathbb{1}_{\Omega_{square}}(x, y), $$

that will be taken into consideration through the Tikhonov regularization term (23), see Sec. 3.1.1.

We include in the benchmark the Tikhonov (T), LASSO, Population Dynamic (PD) regularization, and all combination of them, in order to understand how each regularization term operates and impacts the BI-DA. The optimum $s_a$ is then compared with the target source $s_t$ using the different criteria introduced in Sec. 2.4, i.e. total and support RMSE criteria (see resp. Eqs. (14) and (15)), miss and false discovery rates (see resp. Eqs. (17) and (18)), and the one-sensor adjoint state-based criteria (see Eq. (21)).

### 3.4.1. Visual comparison of estimated pheromone emissions $s_a$.
The estimated sources $s_a$ are displayed for each regularization term, including composite terms, in Fig. 7 at initial ($t = 0s$) and final times ($t = 20s$), for one single (top panel) and composite (bottom panel) regularization terms (see Sec. 3.1), together with the case without regularization ('no reg') taken as reference inference experiment.

At $t = 0s$, only slight discrepancies can be observed between the different regularization terms, either in single (top panel of Fig. 7) or composite ones (bottom panel), except for Tikhonov regularization ('T') that over predicts source emission on the favorable zone. Note that with LASSO regularization ('LASSO'), the retro-plume is less diffuse compared to the reference ('no reg') case since small values have been thresholded by the regularization.

At final time $t = 20s$, when fewer data are available for estimation, single regularization terms do not significantly improve the estimate of $s_a$ made in the reference case ('no reg'), except with the population dynamic-informed regularization term ('PD') which substantially increases prediction accuracy. In composite regularization terms, involving population dynamics is also key to improve predictions (see 'PD + LASSO', 'PD + T' and 'all reg' compared to 'T + LASSO'). Visually, coupling population dynamics with LASSO ('PD + LASSO') slightly improves the target estimate compared to population dynamics alone, while coupling it with a Tikhonov term seems to add false positive predictions outside the target in the favorable zone.

In conclusion, adding prior information on population behavior through the population-dynamics term seems to improve the reconstruction of target emissions $s_t$, but it is not clear if coupling it with other regularization terms is beneficial or not: quantitative assessment must be conducted.

### 3.4.2. RMSE criteria.
In order to have a more quantitative assessment of the discrepancies between $s_a$ and $s_t$ in the different cases (single and composite regularization terms), we introduced in Sec. 2.4.1 $RxMSE_{tot}$ (14), the error computed on the whole computational domain, and $RxMSE_{supp}$ (15), the error computed on the support of the target source only. The latter
<figure><img src="image_10.png" /><figcaption>**Figure 7** – Optimal sources estimated with different regularization terms. Optimal quantity of pheromone emitted $s_a$ obtained by inference with several regularization strategies for the Toy case. Upper panel, from left to right: with no regularization term (no reg), with the population dynamic-informed regularization term (30) (PD), with the LASSO regularization term (24) (LASSO) and with the Tikhonov regularization term (T) (23). Lower panel, from left to right: with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO), with a combination of Tikhonov and LASSO regularization terms (T+LASSO), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg). In both panels, the optimal sources are indicated at $t = 0.0$s (upper row) and $t = 20.0$s (lower row).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7** – Optimal sources estimated with different regularization terms. Optimal quantity of pheromone emitted $s_a$ obtained by inference with several regularization strategies for the Toy case. Upper panel, from left to right: with no regularization term (no reg), with the population dynamic-informed regularization term (30) (PD), with the LASSO regularization term (24) (LASSO) and with the Tikhonov regularization term (T) (23). Lower panel, from left to right: with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO), with a combination of Tikhonov and LASSO regularization terms (T+LASSO), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg). In both panels, the optimal sources are indicated at $t = 0.0$s (upper row) and $t = 19.9$s (lower row).</figcaption></figure>
3.4.2. RMSE criteria. In order to have a more quantitative assessment of the discrepancies between $s_a$ and $s_t$ in the different cases (single and composite regularization terms), we introduced in Sec. 2.4.1 $RxMSE_{tot}$ (14), the error computed on the whole computational domain, and $RxMSE_{supp}$ (15), the error computed on the support of the target source only. The latter enables to assess the inference accuracy on the specific zone where the insects are located. We display in Fig. 8 the time evolution of $RxMSE_{tot}$ (top panel) and $RxMSE_{supp}$ (bottom panel) for the different (single and composite) regularization terms and compared to the reference case without regularization ('no reg').

600

We first note that all the regularization terms provide sensibly equal estimation error inside the target in a large part of the time domain, until $t \approx 17.5$s, see the bottom panel of Fig. 8. This indicates that the main differences come from reconstruction outside the target support and at the end of the time window, when data are lacking.

605

The Tikhonov regularization (Fig. 8, upper panel, green curve) gives higher $RxMSE_{tot}$ than the reference estimate (blue curve), except at the extreme end of the time domain. These discrepancies seems due to the additional emissions added on the favorable square. It is confirmed by the $RxMSE_{supp}$ criteria, which is sensibly the same with Tikhonov regularization (Fig. 8, lower panel, green curve) and without regularization (blue curve), except for the last time points where the additional emissions introduced by the Tikhonov regularization improves a little the estimate on the target support compared to the reference case.

610

The LASSO regularization (Fig. 8, upper panel, pink curve) allows a sensible reduction of $RxMSE_{tot}$ compared to the reference case (blue curve), all along the time domain except at the end. By the soft-thresholding feature it provides, the LASSO regularization thresholds the small values toward zero, specifically outside the target support. However, on the target support, the LASSO regularization (Fig. 8, lower panel, pink curve) does not significantly improve the $RxMSE_{supp}$ criteria compared to the reference case (blue curve). At the end of the time domain, the LASSO regularization keeps on thresholding, even within the target support where insects are emitting, thereby impairing the accuracy of the estimate.

615

Both outside and inside the target support, the population dynamics regularization (Fig. 8, both panels, orange curve) improves the estimation at the end of the time domain compared to the reference case (blue curve). This indicates that the population dynamics adds additional priors that counterbalance the lack of data from the sensor at the end of the time series. Specifically, it allows obtaining a solution that is consistent with known population dynamics, which enables the propagation of well-reconstructed pheromone emission to zones and times where the estimate of the pheromone emission is poor due to the location of the sensors and the intrinsic weaknesses of the retro-propagation process.

620

The use of several regularizations, in particular the combination of population dynamics and Tikhonov regularizations (Fig. 8, brown curves) and the combination of all three regularizations (Fig. 8, red curve), take benefit of the enhancement provided by the LASSO regularization outside the target support and by the population dynamics regularization at the end of the time window to improve the prediction accuracy all along the time interval.

625

630

635

3.4.3. Presence prediction, false discovery and miss rates. In a practical view point, the important information for cultural practice is not the distribution of pheromone emission, but the presence of insect pests that we aim to control. Therefore, we introduced in Sec. 2.4.1 a prediction of the insects presence derived from the pheromone emission map and based on the super level set of Eq. (16) of minimal emission value, see Fig. 4.

640
enables to assess the inference accuracy on the specific zone where the insects are located. We display in Fig. 8 the time evolution of $RxMSE_{tot}$ (top panel) and $RxMSE_{supp}$ (bottom panel) for the different (single and composite) regularization terms in comparison with the reference case without regularization ('no reg').

We first note that all the regularization terms provide sensibly equal estimation error inside the target support in a large part of the time domain, until $t \approx 17.5$s (see the bottom panel of Fig. 8), while estimation error on the whole domain or near final time is more widespread. This indicates that the main differences come from reconstruction outside the target support and at the end of the time window, when data are lacking.

Again, we note that adding population dynamics to the regularization maintains low error values at final time (see 'PD', 'PD + LASSO', 'PD + T', 'all reg', compared to the other terms, bottom and top panels), as the population dynamic model enables to propagate information from the initial time to the final time. However, coupling it with LASSO and Tikhonov reduces the error both inside the target and on the whole domain (see 'all reg', red curves, Fig. 8). Compared to 'PD' term alone (orange curves), this enhancement could come from complementary effects and efficient interactions between the different regularization.

In conclusion, coupling all the regularization terms outperforms all the other combinations of regularization terms for prediction accuracy all along the time window, both on the target and the whole domain. Population dynamics is mandatory to preserve low error levels near final time, when data are lacking.

### 3.4.3. Presence prediction, false discovery and miss rates.
In a practical view point, the important information for cultural practice is not the distribution of pheromone emission (which reflects quantitative distribution of insect population), but the presence of the pests that we aim to control. Therefore, we introduced in Sec. 2.4.1 a prediction of the insects presence derived from the pheromone emission map and based on the super level set of Eq. (16) of minimal emission value, see Fig. 4.

We represent in the upper panel of Fig. 9 the presence maps of insects computed with the different regularization terms and given by the target $s_t$, at three different times ($t = 0.0s, t = 10.0s$ and $t = 20.0s$). Pests are predicted as present inside the delineated color-coded zones. As already anticipated in the visual comparison (see Sec. 3.4.1), the Tikhonov regularization overestimates the presence of pests, while the 'LASSO', 'T + LASSO' and 'no reg' strongly underestimate pest presence at final time.

Cases involving population dynamics regularization ('PD', 'PD + LASSO', 'PD + T', 'all reg') keep covering the target zone all along the time domain by propagating with the population dynamics regularization the presence predicted at initial time. We can however notice a slight overestimation at final time with the case combining all three regularizations ('all reg'), probably due to the Tikhonov term as the boundary of the square favorable zone are slightly visible (Fig. 9, upper panel, $T = 20.0s$, red line). We can also notice a slight underestimation at the final time with the case 'PD + LASSO', probably due to the effect of the LASSO regularization. We also observe spurious presence predictions with the single population dynamics regularization, the single LASSO regularization and the combination of both (Fig. 9, upper panel, resp. 'PD', 'LASSO' and 'PD + LASSO') far from the target, likely because of the noise.

To get a better view of the amount of insects that are detected and missed during the time window, we use false discovery rate $\rho_{fp}$ (18), i.e. the proportion of area where the insects are predicted to be present but are actually absent, and the miss rate $\rho_{fa}$ (17), i.e. the proportion of
<figure><img src="image_11.png" /><figcaption>**Figure 8** – Evolution of total and support RMSE of $s_a$ against $s_t$. Time evolution of the root space-mean square errors between the optimum $s_a$ obtained by inference with several regularization strategies and the target $s_t$, on the whole domain ($R \times \text{MSE}_{\text{tot}}$, Eq. (14), upper panel) and on the target support ($R \times \text{MSE}_{\text{supp}}$, Eq. (15), lower panel) in the Toy case: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (30) (PD, orange), with the LASSO regularization term (24) (LASSO, pink), with the group LASSO regularization term (25) (group LASSO, yellow), with the Tikhonov regularization term (23) (T, green), with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO, brown), with a combination of Tikhonov and LASSO regularization terms (T+LASSO, purple), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red).</figcaption></figure>

To get a better view of the amount of insects that are detected and missed during the time window, we use false discovery rate $\rho_{fp}$ (18), i.e. the proportion of area where the insects are predicted to be present but are actually absent, and the miss rate $\rho_{fa}$ (17), i.e. the proportion of area where the insects are actually present but the prediction is missed. We recall that the objective here is to reduce the false discovery rate, equivalent in our context to a useless application of pesticide in the corresponding zone, while keeping the miss rate as small as possible, or even null, to control the infestation.

We represent in the upper panel of Fig. 9 the presence maps of insects computed with different regularization terms and given by the target $s_t$, at three different times ($t = 0.0s, t = 10.0s$ and $t = 20.0s$). Pests are present inside the delineated color-coded zones. As already anticipated in the visual comparison (see Sec. 3.4.1), the Tikhonov regularization overestimates the presence of pests, while the LASSO regularization, combination of LASSO and Tikhonov regularization ($'T + LASSO'$) and the reference case ('no reg') strongly underestimate pest presence at final time.

Cases involving population dynamics regularization ('PD', 'PD + LASSO', 'PD + T', 'all reg') keep covering the target zone all along the time domain. The presence predicted with $j_{\text{reg,PD}}$ at final time benefits from the propagation by the population dynamics regularization of the presence predicted at initial time. We can however notice a slight overestimation at final time with the case combining all three regularizations ('all reg'), probably due to the Tikhonov term
<figure><img src="image_11.png" /><figcaption>**Figure 8** – Evolution of total and support RMSE of $s_a$ against $s_t$. Time evolution of the root space-mean square errors between the optimum $s_a$ obtained by inference with several regularization strategies and the target $s_t$, on the whole domain ($R \times \text{MSE}_{\text{tot}}$, Eq. (14), upper panel) and on the target support ($R \times \text{MSE}_{\text{supp}}$, Eq. (15), lower panel) in the Toy case: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (30) (PD, orange), with the LASSO regularization term (24) (LASSO, pink), with the group LASSO regularization term (25) (group LASSO, yellow), with the Tikhonov regularization term (23) (T, green), with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO, brown), with a combination of Tikhonov and LASSO regularization terms (T+LASSO, purple), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red).</figcaption></figure>

area where the insects are actually present but the prediction is missed. We recall that the objective here is to reduce the false discovery rate, equivalent in our context to a useless application of pesticide in the corresponding zone, while keeping the miss rate as small as possible, or even null, to control the infestation.

Miss rates $\rho_{fa}$ are comparable for all regularization at the beginning of the simulation, but show a strong miss rate increase at final times for penalties that do not involve '$PD$' or '$T$' (Fig. 9, lower panel, '$no reg$' and '$LASSO$', respectively blue and pink plain curves). When the Tikhonov regularization is involved, the miss rate value is low, but to the price of high values of false discovery rate ($'T'$, $'T+LASSO$ and $'PD+T'$, resp. green, purple abd grey plain and dashed curves). Adding the population dynamics term, the miss rate is controlled all along the time course of the simulation (Fig. 9, lower panel, '$all reg$', '$PD$', '$PD+LASSO$', '$PD+T$', respectively red, orange, brown and gray plain curves). The combination of all three regularizations ($'all reg$' , red curves) seems to be a good trade-off, the miss and false discovery rates are very stable, $\rho_{fa}$ remains zero on the whole time window while $\rho_{fp}$ is among the smallest one on the main part of the time window. With this combination, the three terms interact in an efficient way: the LASSO term seems to keep $\rho_{fp}$ small, the Tikhonov term reduces the $\rho_{fa}$ while the population dynamics term mitigates the rates evolution by propagating information from the initial time to the end of the time domain, when the data are scarce.
<figure><img src="image_12.png" /><figcaption>Figure 9 – Presence prediction and time evolution of false discovery and miss rates.</figcaption></figure>

Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (30) (PD, orange), with the LASSO regularization term (24) (LASSO, pink), with the group LASSO regularization term (25) (group LASSO, yellow), with the Tikhonov regularization term (T, green), with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO, brown), with a combination of Tikhonov and LASSO regularization terms (T+LASSO, purple), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). Top panel: map of presence prediction at three different time. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16). The true presence (i.e. the target support) is indicated in black.

Lower panel. Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), Eq. (18).
<figure><img src="image_12.png" /><figcaption>**Figure 9** – Presence prediction and time evolution of false discovery and miss rates.
Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (30) (PD, orange), with the LASSO regularization term (24) (LASSO, pink), with the Tikhonov regularization term (T, green), with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO, brown), with a combination of Tikhonov and LASSO regularization terms (T+LASSO, purple), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red).
Top panel: map of presence prediction at three different time. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16). The true presence (i.e. the target support) is indicated in black. Lower panel: Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), Eq. (18).</figcaption></figure>
as the boundary of the square favorable zone are slightly visible (Fig. 9, upper panel, $T = 20.0$s, red line). We also observe spurious presence predictions with the single population dynamics regularization, the single LASSO regularization and the combination of both (Fig. 9, upper panel, resp. 'PD', 'LASSO' and 'PD + LASSO') far from the target, likely because of the noise.

Miss rates $\rho_{fa}$ are comparable for all regularization at the beginning of the simulation, the differences being more visible in the false discovery rates $\rho_{fp}$. In all cases, both $\rho_{fp}$ and $\rho_{fa}$ are quite stable throughout the time domain until $t \approx 17.5$s, time at which $\rho_{fa}$ without regularization (reference case 'no reg') strongly increases (Fig. 9, lower panel, blue curves) and where discrepancies occur. As expected, with the Tikhonov regularization (Fig. 9, lower panel, green curves), $\rho_{fa}$ stays close to zero, but the $\rho_{fp}$ remains high. On the contrary, with the LASSO regularization (Fig. 9, lower panel, pink curves), $\rho_{fp}$ is the lowest but $\rho_{fa}$ is among the highest and strongly increases after $t \approx 15$s, which is earlier compared to the reference case (blue curve). With the combination of Tikhonov and LASSO regularizations (Fig. 9, lower panel, purple curves), $\rho_{fa}$ also strongly increases but after $t \approx 19$s, which is later compared to the reference case (blue curve). As soon as population dynamics regularization is involved (Fig. 9, lower panel, orange, brown, grey and red curves), $\rho_{fa}$ is rather low and kept under the one obtained in the reference case ('no reg', blue curves) especially at final time, except with the combination of population dynamics and LASSO regularizations (brown curves) that misses some part of the true presence of insects. In particular, with the combination of all three regularizations ('all reg', red curves), the miss and false discovery rates are very stable, $\rho_{fa}$ remains zero on the whole time window while $\rho_{fp}$ is the third smallest one on the main part of the time window. With this combination, the three terms interact in an efficient way: the LASSO term seems to keep $\rho_{fp}$ small, the Tikhonov term reduces the $\rho_{fa}$ while the population dynamics term smooths the rates evolution at the end of the time domain, when the data are scarce.

#### 3.4.4. Estimating active sensors with the one-sensor adjoint state criteria.
The visual comparison of optimal pheromone emission $s_a$ (see Sec. 3.4.1) indicates that sensors near the emission source seem to have greater impact on the retro-plume than sensors located further away. In order to better quantify this feature, we introduced in Sec. 2.4.3 the maximal total contribution $||c_{i,tot}^*||_\infty$ of the $i^{th}$ sensor, which is the total contribution of the data obtained with the $i^{th}$ given sensor to the decrease of the cost function along the optimization algorithm. This metric enables to understand how a specific sensor impact the estimation of $s_a$, and to rank the sensors accordingly for each regularizations, see Fig. 10. It allows in particular to understand if a regularization is able to supplement or amplify the information collected by the sensors, enabling to use less sensors to achieve a given level of accuracy in the estimate of $s_a$. Note that the sensors are not involved in the regularization terms, so that the impact of the BI-DA regularization terms on the total contribution $||c_{i,tot}^*||_\infty$ is indirect: the regularization changes the trajectory taken by the optimization algorithm in the control space, hence modifying the points where the one-sensor adjoint state $c_i^*$ is evaluated along the optimization path.

We can observe that in the reference case (Fig. 10, 'no reg' and blue line), the 4 first sensors are particularly impactful. Then, a plateau value is attained until the 31st sensor, where a drop is observed, indicating that the last sensors could possibly be removed without impairing the estimate of $s_a$. The LASSO regularization (Fig. 10, pink line) tends to amplify the contribution of all the sensors and mitigates the drops of contribution at the beginning and the end of the contribution order. Indeed, the LASSO regularization tends to threshold the small value of $s$ and the contribution of the sensors are amplified to compensate this thresholding. With Tikhonov and
In conclusion, 'PD' term is very interesting for controlling the miss rate on the whole time window and can be combined to 'T' term for enhanced reduction of $\rho_{fa}$, while coupling them with 'LASSO' reduces the false discovery rate.

### 3.4.4. Estimating active sensors with the one-sensor adjoint state criteria.
The visual comparison of optimal pheromone emission $s_a$ (see Sec. 3.4.1) indicates that sensors near the emission source seem to have greater impact on the retro-plume than sensors located further away. More generally, the number of sensors needed for an accurate prediction is an hyperparameter of the experimental design that could be evaluated. In order to obtain a proxy of this number, we developed a metric enabling to detect and number sensors that most influence the optimization process. To this end, we introduced in Sec. 2.4.3 the maximal total contribution $||c_{i,tot}^*||_\infty$, which is the total contribution of the data obtained with the $i^{th}$ sensor to the decrease of the cost function along the optimization algorithm. This metric enables to understand how a specific sensor impacts the estimation of $s_a$, and to rank the sensors accordingly for each regularizations, see Fig. 10. It allows in particular to understand if a regularization is able to supplement or amplify the information collected by the sensors. Note that the sensors are not involved in the regularization terms, so that the impact of the BI-DA regularization terms on $||c_{i,tot}^*||_\infty$ is indirect: the regularization changes the trajectory taken by the optimization algorithm in the control space, hence modifying the points where the one-sensor adjoint state $c_i^*$ is evaluated along the optimization path.

We can observe that in the reference case (Fig. 10, 'no reg', blue line), the 7 first sensors are particularly impactful. Then, a plateau value is attained until the 32nd sensor, where a drop is observed, indicating that the remaining sensors do not contribute significantly to the estimate of $s_a$. With the Tikhonov regularization, the contribution of the sensors located in the favorable zone is amplified, increasing the impact of the first sensors compared to 'no reg'. In comparison, the population dynamics regularization (Fig. 10, orange line) concentrates most of the information in the first 3 sensors: it enhances their contribution while the remaining sensors reach the plateau value. The combination of 'PD' and 'T' (grey line) shows a similar pattern than Tikhonov. Conversely, the 'LASSO' regularization, and all composite penalties containing 'LASSO' (Fig. 10, pink, brown, purple and red lines) tend to amplify the contribution of all the sensors and mitigates the drops of contribution at the beginning and the end of the contribution order.

In conclusion, by providing prior information, population dynamics and Tikhonov regularizations concentrate the information needed for inference in a limited number of sensors while the LASSO term makes use of all the available sensors during the estimation process. When the three regularizations are combined, sensor contribution is distributed among the whole sensor network, like in the 'LASSO' case. Hence, 'PD' or 'T' regularizations could be used with a reduced number of sensors, while any additional sensor will be leveraged by the 'LASSO' penalty term for thresholding.

## 3.5. Numerical illustration of pheromone propagation in the FAW case

The direct model (1) is now computed on the FAW case, as described in Sec. 2.3.2. The FAW case strongly differs from the Toy case as it relies on a realistic landscape, with real application time and space scales, taking into account a pest, the FAW, harbouring a characteristic emission pattern, see Sec. 2.3.2 for details. Namely, only the females emit pheromones to attract males, but stop emitting at the start of mating. The emission model hence describes the population dynamics of female and male pairing.
population dynamics regularizations (Fig. 10, green and orange lines), the decrease in the contri-
bution is greater for the ten first sensors compared to the reference case (blue line), which shows
that the first sensors in the contribution order concentrate most of the information. Moreover,
the maximal total contribution is almost always lower than the one obtained with the reference
case (blue line). One can conclude that these regularizations bring information that supplement
the signal carried by the sensors data.

In composite regularizations, while combining Tikhonov and LASSO or population dynam-
ics and Tikhonov regularizations (Fig. 10, purple and grey lines) give a similar pattern than the
reference case (blue line), combining population dynamics and LASSO regularizations (Fig. 10,
brown line) provides results close to the single LASSO one (pink line). The combination of all
three regularizations (Fig. 10, red line) globally enhances the overall contribution of each sensor.
Moreover, the drops of contribution at the beginning and the end of the ranking are also signif-
icantly mitigated compared to any other composite regularization. Overall, all the sensors have
very homogeneous contributions: the trend of contribution is similar to that with the LASSO
regularization alone, but amplified both in terms of values and variations.

<figure><img src="image_13.png" /><figcaption>**Figure 10 - Ranking sensor contribution in BI-DA.** Maximal total contribution $||c_{i,tot}^*||_\infty$ vs sensor rank when sorted by contribution order for several regularization strategies in the Toy case: with no regularization (no reg, blue), population dynamic-informed (PD, orange), LASSO (LASSO, pink), Tikhonov (T, green), combination of population dynamic-informed and LASSO (PD+LASSO, brown), combination of Tikhonov and LASSO (T+LASSO, purple), combination of population dynamic-informed and Tikhonov (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). The quantity $||c_{i,tot}^*||_\infty$ indicates the maximal (with respect to time and space) total contribution of sensor *i* all along the optimization path.</figcaption></figure>

**3.5. Numerical illustration of pheromone propagation in the FAW case**

The direct model (1) is now computed on the FAW case, as described in Sec. 2.3.2. The FAW
case strongly differs from the Toy case as it relies on a realistic landscape, with real application
time and space scales, taking into account a pest, the FAW, harbouring a characteristic emission
pattern, see Sec. 2.3.2 for details. Namely, only the females emit pheromones to attract male,
but stop emitting at the start of mating. The emission model hence describes the population
dynamics of female and male pairing.
<figure><img src="image_13.png" /><figcaption>**Figure 10** – Ranking sensor contribution in BI-DA. Maximal total contribution $||c_{i,tot}^*||_\infty$ vs sensor rank when sorted by contribution order for several regularization strategies in the Toy case: with no regularization ('no reg', blue), population dynamic-informed ('PD', orange), LASSO ('LASSO', pink), Tikhonov ('T', green), combination of population dynamic-informed and LASSO ('PD + LASSO', brown), combination of Tikhonov and LASSO ('T+LASSO', purple), combination of population dynamic-informed and Tikhonov ('PD + T', grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms ('all reg', red). The quantity $||c_{i,tot}^*||_\infty$ indicates the maximal (with respect to time and space) total contribution of sensor *i* all along the optimization path.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11** – Map of the concentration in pheromone $c(s_t)$ ($\text{pg.m}^{-2}$) computed from the depth-integrated CTM model (1) with the reference source term $s_t$ for the FAW case presented in Sec. 2.3.2, in logarithmic scale and at several times. In addition iso-lines of $c$ are indicated with black dashed lines. Furthermore, the agricultural parcels are displayed with grey dashed lines together with the reference source $s_t$ support in white at $T = 110.5\text{min}$ (lower right panel).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11** – Map of the concentration in pheromone $c(s_t)$ (pg.m$^{-2}$) computed from the depth-integrated CTM model (1) with the reference source term $s_t$ for the FAW case presented in Sec. 2.3.2, in logarithmic scale and at several times. In addition iso-lines of $c$ are indicated with black dashed lines. Furthermore, the agricultural parcels are displayed with grey dashed lines together with the reference source $s_t$ support in white at $T = 110.5min$ (lower right panel).</figcaption></figure>

Direct model outputs are displayed in Fig. 11. As in the Toy case, convection is also largely predominant in the FAW case: the pheromones are also transported along the wind field, with little additional diffusive effect. Note that the wind field is taken from real weather broadcasting data. The pheromones spread across the domain within a couple of minutes until reaching a quasi-steady state plume within about ten minutes. After that, the decrease in time of $s$ which reflects the emission pattern of FAW that stops releasing pheromones at the start of mating (see Fig. 3c and Sec. 2.3.2), implies that $c$ decrease everywhere until $t \approx 110min$. Finally, $s$ becomes negligible and $c$ does not variate significantly any longer. The landscape composition (see the grey dashed lines on the bottom right panel of Fig. 11), which is taken into account through the loss rate $\tau_{loss}$, does not seem to influence significantly $c$.

### 3.6. Classical VDA and sensor positioning improvement in the FAW case

A first classical VDA without any regularization term is first conducted in the FAW case with sensors randomly dispatched in the field of interest, see left panel in Fig. 12. We highlight that this estimation problem is by far more difficult to solve than the Toy case for several reasons: first, there is a change in the spatio-temporal scales by one order of magnitude (from tens of meters and seconds in the Toy case to hundreds of meters and hours in the FAW case); next, the pheromone emission pattern mimics the biological behavior of the FAW and is not stationary in the FAW case unlike in the Toy case; finally, the sensor density is much higher in the Toy case (~ $10^{-2}$ sensor/m$^2$) than in the FAW case (~ $10^{-4}$ sensor/m$^2$). In this context, classical VDA fails to detect the position of the emitting insects (left panel in Fig. 12). The VDA algorithm does not sufficiently retro-propagate the discrepancies between data and predictions: during the optimization, the VDA could be stuck in a flat zone of the objective function due to insufficient signal of the sensors.

To partially counterbalance this effect, we repositioned several sensors by leveraging the one-sensor adjoint state-based criteria (see Sec. 2.4.3) for the classical VDA. Namely, we computed
Direct model outputs are displayed in Fig. 11. As in the Toy case, convection is also largely predominant in the FAW case: the pheromones are also transported along the wind field, with little additional diffusive effect. Note that the wind field is taken from real weather broadcasting data. The pheromones spread across the domain within a couple of minutes until reaching a quasi-steady state plume within about ten minutes. After that, the decrease in time of the source $s$ which reflects the emission pattern of FAW that stops releasing pheromones at the start of mating (see Fig. 3c and Sec. 2.3.2) implies that the concentration $c$ decreases everywhere until $t \approx 110min$. Finally, $s$ becomes negligible and $c$ does not variate significantly any longer. The landscape composition (see the grey dashed lines on the bottom right panel of Fig. 11), which is taken into account through the loss rate $\tau_{loss}$, does not seem to influence significantly $c$.

### 3.6. Classical VDA and sensor positioning improvement in the FAW case

A first classical VDA without any regularization term is first conducted in the FAW case with sensors randomly dispatched in the field of interest, see left panel in Fig. 12. We highlight that this estimation problem is by far more difficult to solve than the Toy case for several reasons: first, there is a change in the spatio-temporal scales by one order of magnitude (from tens of meters and seconds in the Toy case to hundreds of meters and hours in the FAW case); next, the pheromone emission pattern mimics the biological behavior of the FAW and is no longer stationary like in the Toy case; finally, the sensor density is much higher in the Toy case ($\sim 10^{-2}$ sensor/m$^2$) than in the FAW case ($\sim 10^{-4}$ sensor/m$^2$). In this context, classical VDA fails to detect the position of the emitting insects (left panel in Fig. 12). The VDA algorithm does not sufficiently retro-propagate the discrepancies between data and predictions: during the optimization, the VDA could be stuck in a flat zone of the objective function due to insufficient signal of the sensors.

To partially counterbalance this effect, we repositioned several sensors by leveraging the one-sensor adjoint state-based criteria (see Sec. 2.4.3) for the classical VDA. Namely, we computed $c_{i,tot}^*$ in eq. (21) resulting from this first VDA without regularization, and selected the 2/5 of the sensors with the lowest contribution to the inference of $s$. These sensors were randomly repositioned within the convex envelop of the 3/5 sensors that were kept fixed, with uniform distribution (see top right panel of Fig. 13). We highlight the fact that no knowledge of the actual position of the target pheromone emission (i.e. of the emitting insects) has been used during repositioning, but only a criteria derived from the signal received by the sensors: in that sense, the repositioning method is unsupervised. On the top right panel of Fig. 13, we can observe that the sensors that have been repositioned were those located far or upwind from the targeted emitting insects, indicating that the criteria $c_{i,tot}^*$ is efficient in selecting sensors that would have been a priori identified as potentially uninformative in a supervised context where the target is known.

After sensor repositioning, the reconstructed $s_a$ by classical VDA now overlaps the target position at $T = 0min$, see right panel of Fig. 12. Furthermore, the VDA without regularization now enables to detect the insects in most of the target $s_t$ support at initial time (blue curves on top left panel of Fig. 13), unlike before relocation where no insects were correctly located (green curves). These observations were reflected in the miss and false discovery curves, that were equal to 1 before repositioning, and dropped to about 0.2 and 0.6 respectively after sensor displacement (13 lower panel).

We note that, at $T = 30min$, insects presence prediction vanishes (see Fig. 13 upper panel), and the miss and false discovery rates go up to 1, indicating inoperative presence prediction (see
<figure><img src="image_15.png" /><figcaption>**Figure 12** – VDA without regularization for the FAW case. The optimal quantity of pheromones emitted $s_a$ is obtained by inference without regularization term for the FAW case at $T = 0\,min$ (upper panels) and $T = 30\,min$ (lower panels). The support of the target quantity of emitted pheromone $s_t$ (see Fig. 3b) is delimited by the dashed lines. The left panel represents the source reconstruction with the initial random sensor distribution (uniform law over the maize field), whereas the right panel represents the same estimation after sensor re-positioning (see Sec. 3.6)</figcaption></figure>

$c_{i,tot}^*$ in eq. (21) resulting from this first VDA without regularization, and selected the 2/5 of the sensors with the lowest contribution to the inference of $s$. These sensors were randomly repositioned within the convex envelop of the 3/5 sensors that were kept fixed, with uniform distribution (see top right panel of Fig. 13). We highlight the fact that no knowledge of the actual position of the target pheromone emission (i.e. of the emitting insects) has been used during repositioning, but only a criteria derived from the signal received by the sensors: in that sense, the repositioning method is unsupervised. On the top right panel of Fig. 13, we can observe that the sensors that have been repositioned were those located far or upwind from the targeted emitting insects, indicating that the criteria $c_{i,tot}^*$ is efficient in selecting sensors that would have
<figure><img src="image_15.png" /><figcaption>**Figure 12** – VDA without regularization for the FAW case. The optimal quantity of pheromones emitted $s_a$ is obtained by inference without regularization term for the FAW case at $T = 0\min$ (upper panels) and $T = 30\min$ (lower panels). The support of the target quantity of emitted pheromone $s_t$ (see Fig. 3b) is delimited by the dashed lines. The left panel represents the source reconstruction with the initial random sensor distribution (uniform law over the maize field), whereas the right panel represents the same estimation after sensor re-positioning (see Sec. 3.6).</figcaption></figure>

Fig. 13). This feature is completely different than the drop of presence prediction accuracy observed at final time in the Toy case (see Sec. 3.4 and Fig. 9), that was inherent to boundary effects near final time and due to the retro-propagation feature of the adjoint model. Here, the collapse of discovery accuracy comes from the emission pattern of the FAW and from the accumulation time window of sensors rather than from a limitation of the data assimilation method.

## 3.7. Pheromone emission inference and pest localization with BI-DA in the FAW case

We now test the BI-DA framework in the FAW case. We first characterize the biology-informed regularization terms and benchmark them in this realistic framework.
<figure><img src="image_16.png" /><figcaption>**Figure 13** – Presence prediction and time evolution of false discovery and miss rates in the FAW case before and after sensor repositioning. The prediction of insects presence is obtained by classical VDA without regularization before and after sensors repositioning. Before repositioning, the sensors have been randomly dispatched in the maize field. After computation of $c_{i,tot}^*$, the two-fifth of the sensors with lower contribution to the source reconstruction were randomly re-allocated in the convex envelop of the three-fifth of the sensors with higher contribution. *Upper panel:* presence predictions at $T = 0$min and $T = 30$min before (green line) and after (blue line) sensor repositioning. The position of the sensors before and after relocation is indicated in the right figure. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16) (see Sec. 2.4.2). The true presence (i.e. the target support) is indicated in black. *Lower panel.* Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), in eq. (18).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 13** – Presence prediction and time evolution of false discovery and miss rates in the FAW case before and after sensor repositioning. The prediction of insects presence is obtained by classical VDA without regularization before and after sensors repositioning. Before repositioning, the sensors have been randomly dispatched in the maize field. After computation of $c_{i,tot}^*$, the two-fifth of the sensors with lower contribution to the source reconstruction were randomly re-allocated in the convex envelop of the three-fifth of the sensors with higher contribution. *Upper panel:* presence predictions at $T = 0$min and $T = 30$min before (green line) and after (blue line) sensor repositioning. The position of the sensors before and after relocation is indicated in the right figure. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16) (see Sec. 2.4.2). The true presence (i.e. the target support) is indicated in black. *Lower panel.* Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), in eq. (18).</figcaption></figure>
been a priori identified as potentially uninformative in a supervised context where the target is known.

After sensor repositioning, the reconstructed $s_a$ by classical VDA now overlaps the target position at $T = 0$min, see right panel of Fig. 12. Furthermore, the VDA without regularization now enables to detect the insects in most of the target $s_t$ support at initial time (blue curves on top left panel of Fig. 13), unlike before relocation where no insect were correctly located (green curves). These observations were reflected in the miss and false discovery curves, that were equal to 1 before repositioning, and dropped to about 0.2 and 0.6 respectively after sensor displacement (13 lower panel).

We note that, at $T = 30$min, insects presence prediction vanishes, and the miss and false discovery rates go up to 1, indicating inoperative presence prediction. This feature is completely different than the drop of presence prediction accuracy observed at final time in the Toy case (see Sec. 3.4 and Fig. 9), that was inherent to boundary effects near final time and due to the retro-propagation feature of the adjoint model. Here, the collapse of discovery accuracy comes from the emission pattern of the FAW, rather than from a limitation of the mathematical method. Actually, as FAW females stop emitting the sexual pheromone at the beginning of the mating, the quantity of pheromone measured decreases, impairing insect prediction although the insects are still present, and, what's more, fertilized. For a given level of pheromone, the VDA, and more specifically the adjoint model, retro-propagates the quantity of pheromone measured by the sensor: when the signal is too low, the retro-plume does not reach the target, and the insects are predicted outside the emission zone, see the results at $t = 30$min on the lower right panel of Fig. 12 and on the upper middle panel of Fig. 13.

### 3.7. Pheromone emission inference and pest localization with BI-DA in the FAW case

We now test the BI-DA framework in the FAW case. We first characterize the biology-informed population dynamics regularization term, that must reflect the emission pattern of the FAW. Let us note $p_{ref}^{FAW}(t, x, y)$ the proportion of females that have not mated yet, and therefore keep emitting pheromones. The dynamics of the $p_{ref}^{FAW}$ population is the following:

$$ (31) \qquad \partial_t p_{ref}^{FAW}(t, x, y) = -\gamma(t) p_{ref}^{FAW}(t, x, y) $$

with $\gamma$ the mating rate at a given time that is represented by the logistic function (13). Let us note that $p_{ref}^{FAW}$ defined in eq. (11) is an analytic solution of the population dynamics model (31). Thus, in the FAW case, the following population dynamic-informed regularization term will be considered:

$$ (32) \qquad j_{reg,PD}^{FAW}(s) = \|\partial_t s + \gamma s\|_{L^2}^2 $$

As mentioned in Sec. 2.3.2, the FAW females usually emit pheromone on plants favorable for the larvae, which in this case is assumed to be maize. However, maize covers most of the domain considered and so, it does not provide relevant information. From the biological knowledge of the insects, one can deduce an expected density of insects and a quantity of pheromone emitted per female insect, and thus a local expected value of pheromone emission $\bar{s}_{b,loc}^{FAW}$. Without any additional information, this value will be used on the whole domain at the initial time. This leads to the following Tikhonov regularization:

$$ (33) \qquad s_b^{FAW}(x, y) = \bar{s}_{b,loc}^{FAW} \quad \text{so that} \quad j_{reg,T}^{FAW}(s) = \|s(t=0) - s_b^{FAW}\|_{L^2(\Omega)}^2 $$

Finally, we note that the emission pattern of FAW is ought to be particularly impactful on the classical LASSO term. The soft-threshold will preponderantly mitigate the source reconstruction
3.7.1. Biological priors in the FAW case. First, we define the population dynamics regularization term, that must reflect the emission pattern of the FAW. Let us note $p_{ref}^{FAW}(t, x, y)$ the proportion of females that have not mated yet, and therefore keep emitting pheromones. The dynamics of the $p_{ref}^{FAW}$ population is the following:

$$ (31) \qquad \partial_t p_{ref}^{FAW}(t, x, y) = -\gamma(t) p_{ref}^{FAW}(t, x, y) $$

with $\gamma$ the mating rate at a given time that is represented by the logistic function (13). Let us note that $p_{ref}^{FAW}$ defined in eq. (11) is an analytic solution of the population dynamics model (31). Thus, in the FAW case, the following population dynamic-informed regularization term will be considered:

$$ (32) \qquad j_{reg,PD}^{FAW}(s) = \|\partial_t s + \gamma s\|_{L^2}^2 $$

As mentioned in Sec. 2.3.2, the FAW females usually emit pheromone on plants favorable for the larvae, which in this case is assumed to be maize. However, maize covers most of the domain considered and so, it does not provide relevant information. From the biological knowledge of the insects, one can deduce an expected density of insects and a quantity of pheromone emitted per female insect, and thus a local expected value of pheromone emission $\bar{s}_{b,loc}^{FAW}$. Without any additional information, this value will be used on the whole domain at the initial time. This leads to the following Tikhonov regularization:

$$ (33) \qquad s_b^{FAW}(x, y) = \bar{s}_{b,loc}^{FAW} \quad \text{so that} \quad j_{reg,T}^{FAW}(s) = \|s(t=0) - s_b^{FAW}\|_{L^2(\Omega)}^2 $$

Finally, we note that the emission pattern of FAW is ought to be particularly impactful on the classical LASSO term. The soft-threshold will preponderantly mitigate the source reconstruction when the emission is low, accelerating miss rate at the end of the simulation when the number of emitting females decreases. Conversely, we can expect that the LASSO penalty will not filter out spurious reconstructions at the beginning of the simulation when levels of emitted pheromones are still high since these spurious zones could represent higher quantities than the shrinking threshold. To counterbalance these effects, we decided to switch the LASSO penalty for a group-LASSO penalty by filtering out spatial areas that received a total amount of pheromones, during the whole time domain, under the shrinking threshold (see eq. (25)), so that

$$ j_{reg,gLASSO}^{FAW}(s) = \int_{\Omega} \left( \int_{0}^{T} s(t, x, y)^2 \right)^{1/2} dxdy. $$

A short comparison of LASSO and group-LASSO on the Toy case is made in the supplementary Sec. A.4.

3.7.2. Visual comparison of estimated pheromone emissions $s_a$. We display the map of the reconstructed pheromone emission $s_a$ obtained by BI-DA in fig. 14 for different regularization terms, namely population dynamic-informed ('PD'), group-LASSO, a combination of Tikhonov and group-LASSO ('T + group LASSO'), and a combination of all three terms ('all reg'), that are compared with the reference case without regularization term ('no reg').

We can observe that $s_a$ obtained with the population dynamics regularization (top middle panel on Fig. 14) slightly enhances the one obtained in the reference case. The group LASSO regularization (top right panel on Fig. 14) thresholds pheromone emissions at location where the total emission over time is too small. The combination of 'T + group LASSO' regularizations (bottom left panel on Fig. 14) enables to reconstruct emissions in a large area, since the Tikhonov regularization pulls the emissions on the whole domain toward the right order of magnitude, while the group-LASSO regularization counterbalances that by filtering out emission in areas with low
<figure><img src="image_17.png" /><figcaption>Figure 14 – Optimal quantity of pheromone emitted $s_a$ obtained by inference with several regularization strategies for the FAW case. Upper panel, from left to right: no regularization term (no reg), population dynamic-informed (PD), group-LASSO. Lower panel, from left to right: combination of Tikhonov and group-LASSO (T+group LASSO), combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms (all reg).</figcaption></figure>

when the emission is low, accelerating miss rate at the end of the simulation when the number of emitting females decreases. Conversely, we can expect that the LASSO penalty will not filter out spurious reconstructions at the beginning of the simulation when levels of emitted pheromones are still high since these spurious zones could represent higher quantities than the shrinking threshold. To counterbalance these effects, we decided to switch the LASSO penalty for a group-LASSO penalty by filtering out spatial areas that received a total amount of pheromones, during the whole time domain, under the shrinking threshold (see eq. (25)), so that

$$j_{reg,gLASSO}^{FAW}(s) = \int_{\Omega} \left( \int_{0}^{T} s(t, x, y)^2 \right)^{1/2} dx dy.$$

A short comparison of LASSO and group-LASSO on the Toy case is made in the supplementary Sec. B.4.
<figure><img src="image_17.png" /><figcaption>**Figure 14** – Optimal quantity of pheromone emitted $s_a$ obtained by inference with several regularization strategies for the FAW case. Upper panel, from left to right: no regularization term (no reg), population dynamic-informed (PD), group-LASSO. Lower panel, from left to right: combination of Tikhonov and group-LASSO (T+group LASSO), combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms (all reg).</figcaption></figure>

total emission. Finally, the combination of the three regularizations (bottom middle panel on Fig. 14) proposes a mix between the population dynamics regularization and the combination of group-LASSO and Tikhonov regularizations.

### 3.7.3. RMSE criteria between the optimum $s_a$ and the target $s_t$. 
The RMSE<sub>tot</sub> and RMSE<sub>supp</sub> curves show comparable error magnitude and trends for the different regularization terms.

The 'group LASSO' and '$T + group LASSO$' curves (pink and purple curves, Fig. 15), like the 'no reg' case (blue curve), present bumps around $t = 30\text{min}$ and $t = 60\text{min}$, due to the effect of the pre-concentration window in the observation operator in eq. (10). When adding the population dynamics regularization, this pattern is smoothed out, like in 'PD' (orange) and 'all reg' (red) curves. The '$T + group LASSO$' curve of RMSE<sub>supp</sub> present a bump at $t = 0\text{s}$ due to the Tikhonov regularization term, that is also smoothed on the curve of the 'all reg' case.
<figure><img src="image_18.png" /><figcaption>**Figure 15** – Root space-mean square errors between the optimum *s*<sub>a</sub> and *s*<sub>t</sub> for several regularization strategies. We display the *RxMSE*<sub>tot</sub> introduced in Eq. (14) (upper panel) and the *RxMSE*<sub>supp</sub> introduced in Eq. (15) (lower panel), versus *t* for the FAW case: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (32) (PD, orange), with the group-LASSO regularization term (25) (gLASSO, pink), with the Tikhonov regularization term and group-LASSO (T+group LASSO, brown) and with a combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms (all reg, red).</figcaption></figure>

3.7.1. Visual comparison of estimated pheromone emissions $s_a$. We display the map of the reconstructed pheromone emission $s_a$ obtained by BI-DA in fig. 14 for different regularization terms, namely population dynamic-informed ('PD'), group-LASSO, a combination of Tikhonov and group-LASSO ('T+groupLASSO'), and a combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms ('all reg'), that are compared with the reference case without regularization term ('no reg').

We can observe that $s_a$ obtained with the population dynamics regularization (top middle panel on Fig. 14) slightly enhances the one obtained in the reference case. The group LASSO regularization (top right panel on Fig. 14) thresholds pheromone emissions at location where the total emission over time is too small, and amplifies the remaining pheromone emission to compensate, including spurious emission far from the target zone. The combination of Tikhonov and group LASSO regularizations (bottom left panel on Fig. 14) enables to reconstruct emissions in a large area, since the Tikhonov regularization pulls the emissions on the whole domain toward the high order of magnitude, while the group-LASSO regularization counterbalances that by filtering out emission in areas with low total emission. Finally, the combination of all three regularizations (bottom middle panel on Fig. 14) proposes a mix between the population dynamics regularization and the combination of group-LASSO and Tikhonov regularizations.

3.7.2. RMSE criteria between the optimum $s_a$ and the target $s_t$. The $RMSE_{tot}$ and $RMSE_{supp}$ give comparable insights for the different regularizations. The population dynamics regularization (orange curve on Fig. 15) enables to smooth the curve, especially around $t = 30\text{min}$ and $t = 60\text{min}$ where discontinuities are visible in the reference case (blue curve) due to the pre-concentration effect of the observation operator (10). With the group-LASSO regularization (pink curve on Fig. 15), as consequence of the thresholding and the resulting compensation observed previously,
<figure><img src="image_18.png" /><figcaption>**Figure 15** – Root space-mean square errors between the optimum *s*<sub>a</sub> and *s*<sub>t</sub> for several regularization strategies. We display the *RxMSE*<sub>tot</sub> introduced in Eq. (14) (upper panel) and the *RxMSE*<sub>supp</sub> introduced in Eq. (15) (lower panel), versus *t* for the FAW case: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (32) (PD, orange), with the group-LASSO regularization term (25) (pink), with the Tikhonov regularization term and group-LASSO (T+group LASSO, brown) and with a combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms (all reg, red).</figcaption></figure>

Hence, the different regularization terms have low impact on the overall error, but the prior knowledge included by the population dynamics term filters out discontinuities in the observation operator due to pre-concentration.

3.7.4. Prediction of insect presence. We now focus on insects presence predictions provided by the different regularization terms together with the corresponding miss rate ρ<sub>fa</sub> and false discovery rate ρ<sub>fp</sub> (respectively upper and lower panels of Fig. 16). Without any regularization ('no reg', blue curve, Fig. 16), false discovery and miss rates stay stable until t = 30min, when pheromone emission vanishes. As expected, adding 'group LASSO', with or without Tikhonov (pink and purple curves), reduces false discoveries. Meanwhile, 'T + group LASSO' regularization increases the miss rate while 'group LASSO' alone induces a slight reduction. Adding population dynamics regularization (orange curves) allows to reduce target misses compared to 'no reg', especially during the first 20 minutes, to the price of a higher ρ<sub>fp</sub>. As observed previously on the RMSE criteria, the population dynamics regularization enables to smooth the curves, especially around t = 30min where a discontinuity is visible in the reference case (blue curve). The three regularizations interact between each other in an efficient way (Fig. 16, 'all reg', red curves) since the Tikhonov regularization enables to achieve one of the lowest ρ<sub>fa</sub> at t = 0min, the population dynamics regularization enables to propagate this good initial insect presence prediction and the group-LASSO regularization enables to threshold spurious pheromone emissions and reduces ρ<sub>fp</sub>.

Consequently, combining the different regularizations represents a good trade-off to reduce both false discovery and miss rates, which is the final operational objective for minimizing pesticide use while keeping the amount of undetected pests low.