$R \times \text{MSE}_{\text{tot}}$ is slightly higher and $R \times \text{MSE}_{\text{supp}}$ slightly lower. With the combination of Tikhonov and group-LASSO regularizations (purple curve on Fig. 15), the trend is similar. However, due to the additional Tikhonov regularization that aims to pull the emissions on the whole domain towards the right order of magnitude at the initial time, $\text{RMSE}_{\text{supp}}$ is lower and $\text{RMSE}_{\text{tot}}$ higher, due to additional spurious emissions, at only $t = 0\text{min}$. The combination of all three regularizations (red curve on Fig. 15) enables to harness the advantages offered by each of them: the Tikhonov regularization enables to recover the right order magnitude at $t = 0$, while the group-LASSO regularization filters out the spurious emissions (including the one introduced by the Tikhonov regularization) and the population dynamics regularization propagates the emissions reconstructed at $t = 0\text{s}$ through the population dynamics model (which smooths the $R \times \text{MSE}$ curve). The different regularization terms achieve small errors for $t > 60\text{min}$, when the insects stop to emit, indicating that the drop of pheromone levels is correctly reconstructed.

### 3.7.3. Prediction of insect presence.
We now focus on insects presence predictions provided by the different regularization terms together with the corresponding miss rate $\rho_{fa}$ and false discovery rate $\rho_{fp}$, see respectively the upper and lower panels of Fig. 16. The population dynamics regularization (orange curves on Fig. 16) provides among the lowest $\rho_{fa}$, especially during the first 20 minutes, but a higher $\rho_{fp}$. As observed previously on the RMSE criteria, the population dynamics regularization enables to smooth the curves, especially around $t = 30\text{min}$ where a discontinuity is visible in the reference case (blue curve). This smoothing feature enables also to achieve lower $\rho_{fa}$ and $\rho_{fp}$ just after $t = 30\text{min}$ compared to the reference case (blue curve). With the group-LASSO regularisation (pink curves on Fig. 16), $\rho_{fa}$ and $\rho_{fp}$ stay stable and does not vary significantly except at the $t = 30\text{min}$ discontinuity as in the reference case without regularization (blue curve). However, due to the thresholding feature of the group-LASSO regularization, $\rho_{fa}$ and $\rho_{fp}$ are respectively higher and lower than the one obtained without regularization (blue curve), indicating a too stringent variable selection. With the combination of Tikhonov and group-LASSO regularizations (purple curves on Fig. 16), as observed on the RMSE criteria, the presence of insects is over-predicted, $\rho_{fa}$ is very low and $\rho_{fp}$ very high at $t = 0\text{min}$ due to the Tikhonov regularization implying additional but potential spurious pheromone emission over the whole domain. After that, $\rho_{fa}$ rapidly increases while $\rho_{fp}$ rapidly drops, and the trend are similar to the trend observed with the group-LASSO regularization alone (pink curve), but still with lower values likely due to the interaction between the group-LASSO and Tikhonov regularizations. The three regularizations interact between each other in a efficient way (Fig. 16, red curves) and the combination enables to articulate the advantages offered by each regularizations: the Tikhonov regularization enables to achieve one of the lowest $\rho_{fa}$ at $t = 0\text{min}$, 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 $\rho_{fp}$. The smoothing feature of the population dynamics regularization is also visible on $\rho_{fa}$ around $t = 30\text{min}$ where $\rho_{fa}$ obtained with the reference case (blue curves) is discontinuous. Let us note that with the combination of all three regularizations, in contrast to the case involving only group-LASSO regularization (pink curves), $\rho_{fp}$, while being among the lowest in particular around $t = 20\text{min}$, exhibits a relatively high value at $t = 0\text{min}$. Furthermore, in contrast to the case involving the combination of Tikhonov and group-LASSO regularizations (purple curves), $\rho_{fa}$ takes time to decrease and become one of the lowest among the different regularizations.

### 3.7.4. Estimating active sensors with the one-sensor adjoint state criteria.
We now look at the maximal total contributions $||c_{i,tot}^*||_\infty$ of the sensors and how the regularizations influence these sensors contributions, see Sec. 2.4.3. Similar trends at those observed in the Toy case are visible
<figure><img src="image_19.png" /><figcaption>**Figure 16** – Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (32) (PD, orange), with the group-LASSO regularization term (25) (group LASSO, pink), with Tikhonov and group-LASSO regularization terms (T+group LASSO, purple), and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). Upper panel: map of presence prediction at three different time, true presence in black. Lower panel, miss rate $ρ_{fa}$ (17) (solid lines) and false discovery rate $ρ_{fp}$ (18) (dashed lines) vs t.</figcaption></figure>

3.7.5. Estimating active sensors with the one-sensor adjoint state criteria. We now look at the maximal total contributions $||c_{i,tot}^*||_\infty$ of the sensors and how the regularizations influence these sensors contributions in the FAW case, see Sec. 2.4.3. Similar trends than those observed in the Toy case are visible without regularization (blue curve on Fig. 17): a few sensors mainly contribute, then a plateau is attained and finally, a few sensors without significant contribution. In the FAW case, adding population dynamics has low impact on sensor contribution, except for the last ones. As in the Toy case, adding shrinkage through the 'groupLASSO' penalty (pink, purple and red curves) amplifies and standardizes the contributions of the sensors.
<figure><img src="image_19.png" /><figcaption>**Figure 16** – Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (32) (PD, orange), with the group-LASSO regularization term (25) (group LASSO, pink), with Tikhonov and group-LASSO regularization terms (T+group LASSO, purple), and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). Upper panel: map of presence prediction at three different time, true presence in black. Lower panel, miss rate $ρ_{fa}$ (17) (solid lines) and false discovery rate $ρ_{fp}$ (18) (dashed lines) vs $t$.</figcaption></figure>

here, especially without regularization (blue curve on Fig. 17): a few sensors mainly contribute, then a plateau is attained (less visible than in the Toy case but still observable) and finally, a few sensors (here approximately 4) without significant contribution.

In the FAW case, the addition of the population dynamics regularization (orange curve on Fig. 17) does not strongly influence $||c_{i,tot}^*||_\infty$ since the curve of $||c_{i,tot}^*||_\infty$ vs the sensor rank in the contribution order obtained with the population dynamic regularization is very similar to that obtained in the situation without regularization (blue curves). With the group-LASSO regularization (pink curve on Fig. 17), the contribution of the sensors that already contribute the most in the
<figure><img src="image_20.png" /><figcaption>**Figure 17 - Sensor importance in BI-DA.** Maximal total contribution $||c_{i,tot}^*||_\infty$ vs sensor rank when sorted by contribution order for several regularization strategies in the FAW case. It includes no regularization (no reg, blue), population dynamic-informed (PD, orange), group-LASSO (group LASSO, pink), Tikhonov and group-LASSO (T+group LASSO, purple), 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>

In conclusion, the different penalties tend to act in the same way than in the Toy case. We however note that sensor contributions are less contrasted compared to the Toy case since the less informative sensors have been relocated in the convex hull of the most informative ones.

# 4. Discussion

## 4.1. Incorporating biological knowledge with BI-DA

In this study, we introduced a methodology to infer insect pest presence maps from signals of a network of pheromone sensors dispatched over a field that must be protected against a pest. The inference method is built upon variational data assimilation techniques: the aim is to reconstruct a control variable or a hidden state (here the quantity of pheromone emitted by the insects) by combining a direct model (here the pheromone propagation model derived from the atmospheric fluid dynamics) and data (here from the pheromone sensors). This combination is done by minimizing the distance between the data and an estimate of the observed variable. Once the pheromone emission reconstructed, they are taken as a proxy of the pest presence map. The numerical experiments of classical Data Assimilation (DA) problem highlight that classical DA enables to reconstruct the quantity of pheromone emitted almost only in a vicinity of the sensors. To tackle this issue, we proposed several regularization terms specifically designed to add prior Biological Information (BI) of the insects to improve the accuracy of the optimal solution. In particular, we introduced a population dynamics regularization that explicitly model the behavior of the pest of interest, particularly for pheromone emission. Thus, the BI-DA reconstruction is pulled towards a solution which is consistent with the physiology of the pest. The different regularization terms have been rigorously tested in a toy case, using different criteria, such as RMSE or false discovery and miss rates. The population dynamics regularization reduces the
<figure><img src="image_20.png" /><figcaption>**Figure 17** – **Sensor importance in BI-DA.** Maximal total contribution $||c_{i,tot}^*||_\infty$ vs sensor rank when sorted by contribution order for several regularization strategies in the FAW case. It includes no regularization (no reg, blue), population dynamic-informed (PD, orange), group-LASSO (group LASSO, pink), Tikhonov and group-LASSO (T+group LASSO, purple), 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>

reference case (blue curve) are emphasized, while the contribution of the sensors that contribute
the less in the reference case are lowered. This is due to the fact that the group-LASSO regular-
ization thresholds all the pheromone emissions at places where the total amount of pheromone
are low. Consequently to these thresholding, there is an amplification of the contributions of the
major sensors to compensate for the thresholding of the locations that contribute less.

With the combination of Tikhonov and group-LASSO regularizations (purple curve on Fig. 17), the trend is close to the group-LASSO regularization only but with lower absolute values for $||c_{i,tot}^*||_\infty$, indicating that the Tikhonov regularization brings information at $t = 0\min$ that supplement information brought by the sensors. When adding population dynamics to Tikhonov and group-LASSO (Fig. 17, red curve), the sensor contribution ranking is comparable to the situation with combination of Tikhonov and group-LASSO regularizations (purple curve).

4. Discussion

4.1. Incorporating biological knowledge with BI-DA

In this study, we introduced a methodology to infer insect pest presence maps from signals
of a network of pheromone sensors dispatched over a field that must be protected against a
pest. The inference method is built upon variational data assimilation techniques: the aim is to
reconstruct a control variable or a hidden state (here the quantity of pheromone emitted by the
insects) by combining a direct model (here the pheromone propagation model derived from the
atmospheric fluid dynamics) and data (here from the pheromone sensors). This combination is
done by minimizing the distance between the data and an estimate of the observed variable.
Once the pheromone emission reconstructed, they are taken as a proxy of the pest presence
RMSE, especially at the end of the time window when the data are missing, and reduces the miss rates without too much enhancing the false discovery rate, compared to the reference case without regularization. But using the population dynamics term conjointly with Tikhonov and LASSO terms provide the best results regarding presence/absence prediction, pulling the miss rates towards zero with lower rate of false discovery. However, if few sensors have a strong contribution to the reconstruction with the population dynamics regularization, which allows to expect to get comparable estimates with less sensors, coupling the different regularizations increase the number of sensors that have a strong contribution to the estimate.

## 4.2. Strong vs weak constraint in biology-informed data assimilation

One of the most classical method, referred in this study as classical VDA, consists in considering the direct model as strong constraint of the minimization problem: the estimate of the observed variable is the solution of the direct model, knowing a given estimate of the control variable. On the other hand, if the hidden state to reconstruct is the state variable of the direct model, one can use methods where the direct model acts as weak constraint of the minimization problem: the minimization problem aims to minimize the distance between the data and the estimate of the observed variable computed from the hidden state, and the residual of the direct model, see e.g. Trémolet, 2007.

In the BI-DA method we propose, the pheromone propagation model and the pheromone sensor data are combined with biological information, especially model of the dynamics of the insect population. Whilst the pheromone propagation model is used as strong constraint of the minimization problem, the biological information are incorporated in the minimization problem as a weak constraint through the biology-informed regularization terms, see Sec. 3.1. Thus, the combination of data, pheromone propagation model and biological information is based on a mix of strong and weak constraints of the minimization problem.

One could propose a method where the minimization problem is strongly constrained by both the pheromone propagation model and the population dynamics model. For instance, in such situation, the control to infer would be the initial density of insects $p(t = 0, x, y)$ and the observed variable would be estimated through the coupling of the population dynamics model and the pheromone propagation model in the following way. For an estimate of $p(t = 0)$, the density of insects $p(t, x, y) \ \forall (t, x, y)$ is propagated with the population dynamics model and thereby the quantity of pheromone emitted $s(t, x, y) = q(t)p(t, x, y)$ can be computed. Then, the pheromone propagation model (given $s(t, x, y) \ \forall (t, x, y)$) is used to evaluate the concentration of pheromone $c(t, x, y) \ \forall (t, x, y)$ and the observed variable can be computed. One could go even further by coupling the population dynamics model of the female insect population emitting pheromone (depending on the density of male insects), the population dynamics model of the male population (depending on the pheromone concentration and its gradient) and the pheromone propagation model (depending on the density of female insects emitting pheromone). However, these methods based only on strong constraints require a good knowledge of the population dynamics. Actually, adding a strong constraint remains to assume that the observed data are perfectly explained by the model included in the strong constraint, with small modelling error.

On the contrary, the weak constraints resulting from the biological informations included in the BI-DA intrinsically take into account the modelling errors of the underlying biological processes. As the weak constraint tries to minimize the residual of the population dynamics model, an error model can be added to this residual to control its spatio-temporal correlation structure. Therefore, the BI-DA method does not require as much precise knowledge of the pest biology as strong constraints. It is then better suited for situations where the biological informations is
map. As the classical Data Assimilation (DA) problem seems ill-conditioned and enables to re-
construct the quantity of pheromone emitted in a vicinity of the sensors, we proposed several
regularization terms specifically designed to add prior Biological Information (BI) of the insects
to improve the accuracy of the optimal solution. In particular, we introduced a population dy-
namics regularization that explicitly model the behavior of the pest of interest, particularly for
pheromone emission. Thus, the BI-DA reconstruction is pulled towards a solution which is con-
sistent with the physiology of the pest. The different regularization terms have been rigorously
tested in a toy case, using different criteria, such as RMSE or false discovery and miss rates.
The population dynamics regularization reduces the RMSE, especially at the end of the time
window when the data are missing, and reduces the miss rates without too much enhancing
the false discovery rate, compared to the reference case without regularization. But using the
population dynamics term conjointly with Tikhonov and LASSO terms provide the best results
regarding presence/absence prediction, pulling the miss rates towards zero with lower rate of
false discovery. However, if few sensors have a strong contribution to the reconstruction with
the population dynamics regularization, which allows to expect to get comparable estimates with
less sensors, coupling the different regularizations increase the number of sensors that have a
strong contribution to the estimate.

## 4.2. Strong vs weak constraint in biology-informed data assimilation

One of the most classical method, referred in this study as classical VDA, consists in con-
sidering the direct model as strong constraint of the minimization problem: the estimate of the
observed variable is the solution of the direct model, knowing a given estimate of the control
variable. On the other hand, if the hidden state to reconstruct is the state variable of the direct
model, one can use methods where the direct model acts as weak constraint of the minimization
problem: the minimization problem aims to minimize the distance between the data and the es-
timate of the observed variable computed from the hidden state, and the residual of the direct
model, see e.g. Trémolet, 2007.

In the BI-DA method we propose, the pheromone propagation model and the pheromone
sensor data are combined with biological information, especially model of the dynamics of the
insect population. Whilst the pheromone propagation model is used as strong constraint of the
minimization problem, the biological information are incorporated in the minimization problem
as a weak constraint through the biology-informed regularization terms, see Sec. 3.1. Thus, the
combination of data, pheromone propagation model and biological information is based on a mix
of strong and weak constraints of the minimization problem.

One could propose a method where the minimization problem is strongly constrained by
both the pheromone propagation model and the population dynamics model. For instance, in
such situation, the control to infer would be the initial density of insects p(t = 0, x, y) and the
observed variable would be estimated through the coupling of the population dynamics model
and the pheromone propagation model in the following way. For an estimate of p(t = 0), the
density of insects p(t, x, y) ∀(t, x, y) is propagated with the population dynamics model and
thereby the quantity of pheromone emitted s(t, x, y) = q(t)p(t, x, y) can be computed. Then,
the pheromone propagation model (given s(t, x, y) ∀(t, x, y)) is used to evaluate the concentra-
tion of pheromone c(t, x, y) ∀(t, x, y) and the observed variable can be computed. One could
go even further by coupling the population dynamics model of the female insect population
emitting pheromone (depending on the density of male insects), the population dynamics model
of the male population (depending on the pheromone concentration and its gradient) and the
pheromone propagation model (depending on the density of female insects emitting pheromone).
only partially known or may be subject to errors or uncertainties. Moreover, strongly constrain-
ing by both the pheromone propagation model and the population dynamics model could make
the minimization problem more difficult to solve and lead to identifiability issues. Indeed, the
minimum of the cost function may not be unique if the population dynamics is non-linear, and
adding multiple strong constraints may worsen the conditioning of the cost function.

4.3. Toward an effective use of BI-DA for pheromone sensor analysis

The proposed BI-DA method has been developped with consideration of practical constraints encountered in real applications. However, several points still require improvements before an effective use in an agricultural context. These improvements must be carried out in the framework of a global benefit-risks analysis of replacing pesticide-based prophylaxis by BI-DA-guided treatments: a life cycle analysis (LCA) of the sensors should be conducted, including manufacture, transport, deployment, use, recycling and computations and compared to the LCA of pesticides. Beside carbon emissions, benefit-risks analysis should also compare the impact of the two practices on health, environment and biodiversity.

First, computation time can be redhibiting for real applications, since solving BI-DA problem presently requires tens of hours in HPC facilities. This computational load is mostly due to solving the direct and adjoint pheromone propagation models, and to the number of iterations needed by the gradient descent algorithm. To speed-up direct and adjoint models computations, adapted metamodel technics, such as POD-based methods, can be studied. In this framework, PDE computations are replaced by fast approximations obtained by the surrogate model. Moreover, such metamodels can be learned from more complex models than the 2D pheromone propagation model used in this study. For instance, a 3D pheromone propagation model that accounts for vertical variability could be considered. While such models may be computationally infeasible for real-time applications, they can be solved and used offline to train a that could be used in-line for BI-DA. Speeding-up the gradient descent algorithm can be obtained by improving the convergence rate, thus reducing the number of iterations needed for convergence, using preconditionning technics.

Second, the selectivity of pheromone sensors, which is key for precise detection, has as the counterpart that a new sensor network should be deployed for a new target species. From a modelling point of view, increasing the number of target species will not change the CTM direct model which will be identical to all the species, so that a sub-linear increase of DA computation load can be expected. However, as the new species may have its proper behavior, BI-DA regularisation terms must be adapted to the new target, inducing model development time.

Finally, we chose for this study a data assimilation framework where the whole time series
of sensor signals are obtained at once and available for inference. A more sequential frame-
work where presence predictions are updated when new data are available would be usefull for
an adaptative and sequential sensor placement and data acquisition. This framework is particu-
larly suitable for epidemio-surveillance where a dynamic screening is key to adapt the treatment
to the pest dynamics. We note that our variational approach is mathematically equivalent to
Kalman smoothing on the whole time window whenever the biology informed-regularization
terms are differentiable and quadratic. However, we note that the LASSO and group-LASSO
terms are not differentiable and that a non-linear population dynamics model would result in a
non-quadratic PD term. Furthermore, the variational approach avoid to build and store the gain
and covariance matrices needed for Kalman filtering, thus reducing the computational load.
955

However, these methods based only on strong constraints require a good knowledge of the pop-
ulation dynamics. Actually, adding a strong constraint remains to assume that the observed data
are perfectly explained by the model included in the strong constraint, with small modelling error.

960

On the contrary, the weak constraints resulting from the biological informations included in
the BI-DA intrinsically take into account the modelling errors of the underlying biological pro-
cesses. As the weak constraint tries to minimize the residual of the population dynamics model,
an error model can be added to this residual to control its spatio-temporal correlation structure.
Therefore, the BI-DA method does not require as much precise knowledge of the pest biology
as strong constraints. It is then better suited for situations where the biological informations is
only partially known or may be subject to errors or uncertainties. Moreover, strongly constrain-
ing by both the pheromone propagation model and the population dynamics model could make
the minimization problem more difficult to solve and lead to identifiability issues. Indeed, the
minimum of the cost function may not be unique if the population dynamics is non-linear, and
adding multiple strong constraints may worsen the conditioning of the cost function.

965

4.3. Uncertainty sources in the model

970

The proposed BI-DA method and the inference of pheromone sources are subject to various
source of uncertainties. These uncertainties are of several types: uncertainties related to the en-
vironmental parameters, to the pheromone propagation model, to the test cases, to the sensors
and the data. In this section, these uncertainties are briefly presented and discussed.

975

*Uncertainties related to the environmental parameters.* The pheromone propagation is modelled using the 2D CTM (1). Thus, the resulting pheromone concentration depends on environmental parameters such as the wind velocity field $\vec{u}$ and the loss coefficient $\tau_{loss}$. However, the estimation of these two environmental parameters is itself uncertain.

980

In the realistic FAW case, the wind velocity is estimated using the data provided by the Arome
model, see Sec. 2.3.2. This model have a spatial sampling of 0.01°(≈ 1.3km) and a temporal
sampling of 1h. Thus, the smaller scale are not taken into account in the wind data. On the other
hand, the spatial and temporal scale considered in realistic application are smaller: of the order of
magnitude of a couple of minutes and of tens of meters. Therefore, in this context, uncertainties
related to the small spatio-temporal scale of the wind field could be taken into account.

985

The loss coefficient $\tau_{loss}$ enables to parametrize all the vertical fluxes not directly modelled in the 2D model, see Sec. B.1.2. This coefficient depends in particular on the vegetation that absorbs a part of the pheromone that settles on the ground, and on the atmosphere stability. In this study and up to now, this coefficient is calibrated so the 2D CTM fits to a particular solution of the 3D CTM, see Sec. B.1.3. However, the calibration of $\tau_{loss}$ with respect to the land use remains rather empirical and source of uncertainties.

990

In order to understand the impact of these uncertainties on the pheromone propagation
and on the pheromone source inference, an uncertainty quantification and a sensitivity analysis
of the model with respect to these environmental parameters could be conducted, taking into
account these different scales.

995

*Uncertainties related to the pheromone propagation model.* The 2D CTM model is derived from a 3D CTM model by averaging the 3D model over an atmospheric layer depth. This approximation is source of uncertainties since potential updrafts or vertical turbulence are roughly taken into account through the loss term. Additional modelling errors can come from inaccurate model parametrization or inaccurate modelling hypothesis. These uncertainties could be tackled by
4.4. Uncertainty sources in the model

The proposed BI-DA method and the inference of pheromone sources are subject to various
source of uncertainties. These uncertainties are of several types: uncertainties related to the en-
vironmental parameters, to the pheromone propagation model, to the test cases, to the sensors
and the data.

*Uncertainties related to the environmental parameters.* The pheromone propagation is modelled using the 2D CTM (1). Thus, the resulting pheromone concentration depends on environmental parameters such as the wind velocity field $\vec{u}$ and the loss coefficient $\tau_{loss}$. However, the estimation of these two environmental parameters is itself uncertain.

In the realistic FAW case, the wind velocity is estimated using the data provided by the Arome model, see Sec. 2.3.2. This model have a spatial sampling of 0.01°(≈ 1.3km) and a temporal sampling of 1h. Thus, the smaller scale are not taken into account in the wind data. On the other hand, the spatial and temporal scale considered in realistic application are smaller: of the order of magnitude of a couple of minutes and of tens of meters. Therefore, in this context, uncertainties related to the small spatio-temporal scale of the wind field could be taken into account.

The loss coefficient $\tau_{loss}$ enables to parametrize all the vertical fluxes not directly modelled in the 2D model, see Sec. A.1.2. This coefficient depends in particular on the vegetation that absorbs a part of the pheromone that settles on the ground, and on the atmosphere stability. In this study and up to now, this coefficient is calibrated so the 2D CTM fits to a particular solution of the 3D CTM, see Sec. A.1.3. However, the calibration of $\tau_{loss}$ with respect to the land use remains rather empirical and source of uncertainties.

In order to understand the impact of these uncertainties on the pheromone propagation
and on the pheromone source inference, an uncertainty quantification and a sensitivity analysis
of the model with respect to these environmental parameters could be conducted, taking into
account these different scales.

*Uncertainties related to the pheromone propagation model.* The 2D CTM model is derived from a 3D CTM model by averaging the 3D model over an atmospheric layer depth. This approximation is source of uncertainties since potential vertical variability of the wind in the layer, updrafts or vertical turbulence are roughly taken into account through the loss term. Additional modelling errors can come from inaccurate model parametrization or inaccurate modelling hypothesis. These uncertainties could be tackled by devoting more effort to modelling pheromone propagation. This includes the use of models that takes into account the random nature of the uncertainties and of metamodels that potentially take into account the vertical variability of the atmosphere.

*Uncertainties related to the test cases.* To illustrate the performance of the BI-DA method compared to a classical DA method, the data are generated from a target pheromone emission that we aim to retrieve. In these test cases, the computation domain, the location of the target emission and of the sensors and more generally the landscape design or the meteorological information are sources of uncertainties. A systematic quantification of these uncertainties can be performed using sensitivity analysis to assess the impact of these parameters on the BI-DA problem. In particular, the importance of some simplificatory assumption can be explored.

For instance, in both Toy and FAW cases, it is assumed that all the pheromone emissions are inside the computation domain. If some pheromone emission are located out of the computation domain, one might expect that small spurious pheromone emissions will be retrieved inside the
devoting more effort to modelling pheromone propagation. This includes the use of models that takes into account the random nature of the uncertainties and of metamodels.

***Uncertainties related to the test cases.*** To illustrate the performance of the BI-DA method compared to a classical DA method, the data are generated from a target pheromone emission that we aim to retrieve. In these test cases, the computation domain, the location of the target emission and of the sensors could be taken with much worst characteristics. For instance, in both Toy and FAW cases, it is assumed that all the pheromone emissions are inside the computation domain. If some pheromone emission are located out of the computation domain, one might expect that small spurious pheromone emissions will be retrieved inside the domain by the inference to compensate influx through the domain boundaries, that are not consider in the inverse model.

However, as done in the FAW case with the sensors reallocation, see Sec. 3.6, this type of uncertainty can be bypassed following the following steps. First, we start over a large domain and perform the inference. Then, we perform the relocation process based on the one-sensor adjoint state-based criteria as done in the FAW case and reduce the computation domain (and potentially refine the mesh) to a narrower domain around the new sensors location. These steps can be repeated several times.

***Uncertainties related to the sensors and the data.*** The data used in this study are synthetic: they are generated from a target pheromone emission and an observation model that includes a pre-concentration phase, see Sec. 2.3.2. As soon as suitable sensors are available, the classical DA and new BI-DA methods should be tested in progressively less controlled environment and conditions. These tests should start with controlled laboratory experiment with pheromone spray and controlled wind condition, progressing to tests conducted in open fields with uncontrolled environment. This will enable us to test the proposed methods in real condition and improve them by taking into account the properties of the sensors and the statistics of their data.

## 4.4. Sensor placement and reallocation

In BI-DA, we try to leverage spatio-temporal data obtained through a network of sensors that measure spatio-temporal concentration maps of pheromone. The time and space sampling of these map is key for the accuracy of the inference, as illustrated in the FAW case, see Sec. 3.6. The time sampling is given by the intrinsic functioning of the sensors. Due to the small amount of pheromone in the landscape, pre-concentration chambers accumulating the signal to detect are necessary, which determine the minimal time step of the time sampling. The space sampling is given by the spatial repartition of the sensor over the landscape.

In this study, we did not investigate the sensitivity of the inverse problem to the spatial placement of the sensors. Optimal design for building sampling design minimizing the variance of the inference can be constructed in comparable spatio-temporal inference problems (Parisey et al., 2022). This optimal design indicates, given environmental and landscape characteristics, where the sensors should be placed in order to get the best possible confidence on the estimate. Optimal placement is beyond the scope of this study and will be explored in future works.

However, we introduced a specific metric measuring the maximal total contribution of a given sensor. This metric allows to order the sensors by contribution in the inference, hence distinguishing between key and accessory sensors. One could then imagine an iterative method for a dynamics placement of the sensors where the marginally contributing sensors are reallocated and moved to another place where they could enhance the accuracy of pest presence/absence estimates. We investigated how this criteria was impacted for varying BI-DA regularization terms. We found that some regularization terms supplement the information taken out of most sensors
domain by the inference to compensate influx through the domain boundaries, that are not considered in the inverse model. However, as done in the FAW case with the sensors reallocation, see Sec. 3.6, this type of uncertainty can be bypassed using the following steps. First, we start over a large domain and perform the inference. Then, we perform the relocation process based on the one-sensor adjoint state-based criteria as done in the FAW case and reduce the computation domain (and potentially refine the mesh) to a narrower domain around the new sensors location. These steps can be repeated several times.

***Uncertainties related to the sensors and the data.*** The data used in this study are synthetic: they are generated from a target pheromone emission and an observation model that includes a pre-concentration phase, see Sec. 2.3.2. As soon as suitable sensors are available, the classical DA and new BI-DA methods should be tested in progressively less controlled environment and conditions. These tests should start with controlled laboratory experiment with pheromone spray and controlled wind condition, progressing to tests conducted in open fields with uncontrolled environment. The technological developments of pheromone sensors will aim to obtain compact devices, making it possible to move the sensors in order to implement adaptative strategies such as sensor repositioning. This will enable us to test the proposed methods in real condition and improve them by taking into account the properties of the sensors and the statistics of their data. Efficient data acquisition strategies should be considered in order to mitigate data acquisition uncertainties, such as optimal experimental design or dynamical sensor repositioning strategies.

## 4.5. Sensor placement and reallocation

In BI-DA, we try to leverage spatio-temporal data obtained through a network of sensors that measure spatio-temporal concentration maps of pheromone. The time and space sampling of these map is key for the accuracy of the inference, as illustrated in the FAW case, see Sec. 3.6. The time sampling is given by the intrinsic functioning of the sensors. Due to the small amount of pheromone in the landscape, pre-concentration chambers accumulating the signal to detect are necessary, which determine the minimal time step of the time sampling. The space sampling is given by the spatial repartition of the sensor over the landscape.

In this study, we did not investigate the sensitivity of the inverse problem to the spatial placement of the sensors. Optimal design for building sampling design minimizing the variance of the inference can be constructed in comparable spatio-temporal inference problems (Parisey et al., 2022). This optimal design indicates, given environmental and landscape characteristics, where the sensors should be placed in order to get the best possible confidence on the estimate. It also provides, after hyperparameter selection, the number of sensors necessary to achieve a given accuracy.

Optimal placement is beyond the scope of this study and will be explored in future works. However, we introduced a specific metric measuring the maximal total contribution of a given sensor. This metric allows to order the sensors by contribution in the inference, hence distinguishing between key and accessory sensors. One could then imagine an iterative method for a dynamics placement of the sensors where the marginally contributing sensors are reallocated and moved to another place where they could enhance the accuracy of pest presence/absence estimates. We investigated how this criteria was impacted for varying BI-DA regularization terms. We found that some regularization terms supplement the information taken out of most sensors and only a few key sensors still significantly contribute to the inference. This is the case for population dynamics and Tikhonov regularization terms. As these regularization terms highlight the
and only a few key sensors still significantly contribute to the inference. This is the case for pop-
ulation dynamics and Tikhonov regularization terms. As these regularization terms highlight the
sensors that do not significantly contribute in the light of both the propagation of the pheromone
and the biological information, they can be used to study the placement and reallocation of sen-
sors with regards to the pheromone propagation, including the atmosphere dynamics and envi-
ronmental parameter, and to the biological knowledge of the insects, including the dynamics of
its population. On the contrary, other regularization terms amplify the information taken out of
the sensors and the contribution of every sensors. This is particularly the case for regularizations
involving a LASSO term. This stresses the fact that optimal placement strategies will depend on
the selected regularization terms.

5. Conclusion

We introduced in this study a numerical framework to solve the inverse problem for pest
detection from the assimilation of pheromone concentration data. We specifically worked on
adding prior biological knowledge on the population dynamics and the behavior of the pest
to the inference problem. Diverse regularizations to biologically inform the data assimilation
have been proposed and benchmarked. The population dynamics model enables to propagate
emissions well estimated by the data to places and times that are not covered by the assimilation
of the data set.

To reduce the pest miss rate, synonymous of a better control of the pest infestation, without
increasing too much the false discovery rate, synonymous of unnecessary pesticide application,
the following conclusions can be formulated:

(1) Some regularization terms, such as LASSO or group-LASSO, filter out spurious signals, hence enhancing the information provided by the sensors. This leads to lower false discovery rate but higher miss rate.

(2) Some regularization terms, such as the population dynamic-informed or Tikhonov regularization terms, provide prior information that can partially supplement the signal provided by the sensors, leading to lower miss rate, sometimes to the cost of a higher false discovery rate.

(3) The sensor contribution to the inference can be measured with the one-sensor adjoint state criteria. Regularization terms can decrease the contribution of some sensors and enable to efficiently discriminate sensors with high and low contributions. This implies that less sensors are required to achieve equivalent accuracy in the inference. Moreover, sensors with the lowest contribution can be moved to increase inference accuracy of emission sources. In this respect, population-dynamic or Tikhonov regularization terms are convenient to benchmark sensor position.

(4) Once sensor location is fixed, cumulating regularization terms can reduce miss rate while keeping false discovery rates in reasonable values. In particular, it is worth using a composite regularization mixing knowledge on favorable habitat and population dynamics and thresholding spurious pheromone emissions.

Acknowledgments

Funding

The authors acknowledge the support of the French National Research Agency (ANR) under
the grant ANR-20-PCPA-0007 (PheroSensor).
sensors that do not significantly contribute in the light of both the propagation of the pheromone and the biological information, they can be used to study the placement and reallocation of sensors with regards to the pheromone propagation, including the atmosphere dynamics and environmental parameter, and to the biological knowledge of the insects, including the dynamics of its population. On the contrary, other regularization terms amplify the information taken out of the sensors and the contribution of every sensors. This is particularly the case for regularizations involving a LASSO term. This stresses the fact that optimal placement strategies will depend on the selected regularization terms.

## 5. Conclusion

We introduced in this study a numerical framework to solve the inverse problem for pest detection from the assimilation of pheromone concentration data. We specifically worked on adding prior biological knowledge on the population dynamics and the behavior of the pest to the inference problem. Diverse regularizations to biologically inform the data assimilation have been proposed and benchmarked. The population dynamics model enables to propagate emissions well estimated by the data to places and times that are not covered by the assimilation of the data set.

To reduce the pest miss rate, synonymous of a better control of the pest infestation, without increasing too much the false discovery rate, synonymous of unnecessary pesticide application, the following conclusions can be formulated:

(1) Some regularization terms, such as LASSO or group-LASSO, filter out spurious signals, hence enhancing the information provided by the sensors. This leads to lower false discovery rate but higher miss rate.

(2) Some regularization terms, such as the population dynamic-informed or Tikhonov regularization terms, provide prior information that can partially supplement the signal provided by the sensors, leading to lower miss rate, sometimes to the cost of a higher false discovery rate.

(3) The sensor contribution to the inference can be measured with the one-sensor adjoint state criteria. Regularization terms can decrease the contribution of some sensors and enable to efficiently discriminate sensors with high and low contributions. This implies that less sensors are required to achieve equivalent accuracy in the inference. Moreover, sensors with the lowest contribution can be moved to increase inference accuracy of emission sources. In this respect, population-dynamic or Tikhonov regularization terms are convenient to benchmark sensor position.

(4) Once sensor location is fixed, cumulating regularization terms can reduce miss rate while keeping false discovery rates in reasonable values. In particular, it is worth using a composite regularization mixing knowledge on favorable habitat and population dynamics and thresholding spurious pheromone emissions.

## Acknowledgments

Preprint version 2 of this article has been peer-reviewed and recommended by Peer Community In Mathematical & Computational Biology (https://doi.org/10.24072/pci.mcb.100313; Tannier, 2025)
Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of
interest in relation to the content of the article.

6. References

1090

Adam L, Branda M (2016). *Sparse Optimization for Inverse Problems in Atmospheric Modelling*. *Environmental Modelling & Software* **79**, 256-266. https://doi.org/10.1016/j.envsoft.2016.02.002. (Visited on 10/09/2023).

1095

Adamczyk-Chauvat K, Kassa M, Papaïx J, Kiêu K, Stoica RS (2024). *Statistical Inference for Random T-tessellations Models*. Application to Agricultural Landscape Modeling. Annals of the Institute of Statistical Mathematics **76**, 447-479. https://doi.org/10.1007/s10463-023-00893-3. (Visited on 05/01/2024).

Bocquet M (2014). *Introduction to the Principles and Methods of Data Assimilation in the Geosciences*.

1100

Boukal DS, Berec L (2009). *Modelling Mate-Finding Allee Effects and Populations Dynamics, with Applications in Pest Control*. *Population Ecology* **51**, 445-458. https://doi.org/10.1007/s10144-009-0154-4.

Britton N (1989). *Aggregation and the Competitive Exclusion Principle*. *Journal of Theoretical Biology*. https://doi.org/10.1016/S0022-5193(89)80189-4.

1105

Conchou L, Lucas P, Meslin C, Proffit M, Staudt M, Renou M (2019). *Insect Odorscapes: From Plant Volatiles to Natural Olfactory Scenes*. *Frontiers in Physiology* **10**, 972. https://doi.org/10.3389/fphys.2019.00972. (Visited on 07/28/2023).

Deutsch CA, Tewksbury JJ, Tigchelaar M, Battisti DS, Merrill SC, Huey RB, Naylor RL (2018). *Increase in Crop Losses to Insect Pests in a Warming Climate*. *Science* **361**, 916-919. https://doi.org/10.1126/science.aat3466. (Visited on 06/26/2023).

1110

EFSA PLH Panel et al. (2018). *Pest Risk Assessment of Spodoptera Frugiperda for the European Union*. *EFSA Journal* **16**. https://doi.org/10.2903/j.efsa.2018.5351. (Visited on 06/26/2023).

Fisher BEA, MacQueen JF (1981). *A Theoretical Model for Particulate Transport from an Elevated Source in the Atmosphere*. *IMA Journal of Applied Mathematics* **27**, 359-371. https://doi.org/10.1093/imamat/27.3.359. eprint: https://academic.oup.com/imamat/article-pdf/27/3/359/1735525/27-3-359.pdf.

1115

Fisher RA (1937). *The Wave of Advance of Advantageous Genes*. *Annals of eugenics* **7**, 355-369.

Furter J, Grinfeld M (1989). *Local vs. Non-Local Interactions in Population Dynamics*. *Journal of Mathematical Biology*. https://doi.org/10.1007/BF00276081.

1120

Ganga Prasath S, Mahadevan L (2023). *Rheomergy: Collective Behaviour Mediated by Active Flow-Based Recruitment*. *Proceedings of the Royal Society A* **479**, 20220470. https://doi.org/10.1098/rspa.2022.0470.

Gebbers R, Adamchuk VI (2010). *Precision Agriculture and Food Security*. *Science* **327**, 828-831. https://doi.org/10.1126/science.1183899. (Visited on 07/27/2023).

1125

Genieys S, Volpert V, Auger P (2006). *Pattern and Waves for a Model in Population Dynamics with Nonlocal Consumption of Resources*. *Mathematical Modelling of Natural Phenomena*. https://doi.org/10.1051/mmnp:2006004.

IPPC Secretariat. (2021). *Scientific Review of the Impact of Climate Change on Plant Pests*. FAO on behalf of the IPPC Secretariat. https://doi.org/10.4060/cb4769en. (Visited on 06/26/2023).

1130
Funding

The authors acknowledge the support of the French National Research Agency (ANR) under
the grant ANR-20-PCPA-0007 (PheroSensor).

Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of
interest in relation to the content of the article.

Data, script, code, and supplementary information availability

The toolbox containing the core functions for BI-DA is available at https://forgemia.inra.fr/pherosensor/pherosensor-toolbox and (Malou et al., 2024b). The companion codes needed to reproduce the different computation and figures of this paper are available at https://forgemia.inra.fr/pherosensor/companion-code-bi-da and (Malou et al., 2024a).

References

Adam L, Branda M (2016). *Sparse Optimization for Inverse Problems in Atmospheric Modelling*. *Environmental Modelling & Software* **79**, 256–266. https://doi.org/10.1016/j.envsoft.2016.02.002.

Adamczyk-Chauvat K, Kassa M, Papaïx J, Kiêu K, Stoica RS (2024). *Statistical Inference for Random T-tessellations Models*. Application to Agricultural Landscape Modeling. *Annals of the Institute of Statistical Mathematics* **76**, 447–479. https://doi.org/10.1007/s10463-023-00893-3.

Bocquet M (2014). *Introduction to the Principles and Methods of Data Assimilation in the Geosciences*.

Boukal DS, Berec L (2009). *Modelling Mate-Finding Allee Effects and Populations Dynamics, with Applications in Pest Control*. *Population Ecology* **51**, 445–458. https://doi.org/10.1007/s10144-009-0154-4.

Britton N (1989). *Aggregation and the Competitive Exclusion Principle*. *Journal of Theoretical Biology*. https://doi.org/10.1016/S0022-5193(89)80189-4.

Conchou L, Lucas P, Meslin C, Proffit M, Staudt M, Renou M (2019). *Insect Odorscapes: From Plant Volatiles to Natural Olfactory Scenes*. *Frontiers in Physiology* **10**, 972. https://doi.org/10.3389/fphys.2019.00972.

Deutsch CA, Tewksbury JJ, Tigchelaar M, Battisti DS, Merrill SC, Huey RB, Naylor RL (2018). *Increase in Crop Losses to Insect Pests in a Warming Climate*. *Science* **361**, 916–919. https://doi.org/10.1126/science.aat3466.

EFSA PLH Panel et al. (2018). *Pest Risk Assessment of Spodoptera Frugiperda for the European Union*. *EFSA Journal* **16**. https://doi.org/10.2903/j.efsa.2018.5351.

Fisher BEA, MacQueen JF (1981). *A Theoretical Model for Particulate Transport from an Elevated Source in the Atmosphere*. *IMA Journal of Applied Mathematics* **27**, 359–371. https://doi.org/10.1093/imamat/27.3.359. eprint: https://academic.oup.com/imamat/article-pdf/27/3/359/1735525/27-3-359.pdf.

Fisher RA (1937). *The Wave of Advance of Advantageous Genes*. *Annals of eugenics* **7**, 355–369.

Furter J, Grinfeld M (1989). *Local vs. Non-Local Interactions in Population Dynamics*. *Journal of Mathematical Biology*. https://doi.org/10.1007/BF00276081.

Ganga Prasath S, Mahadevan L (2023). *Rheomergy: Collective Behaviour Mediated by Active Flow-Based Recruitment*. *Proceedings of the Royal Society A* **479**, 20220470. https://doi.org/10.1098/rspa.2022.0470.
Isakov V (2017). *Inverse Problems for Partial Differential Equations*. Vol. 127. Applied Mathematical Sciences. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-51658-5. (Visited on 10/09/2023).

1135

Kaltenbacher B, Neubauer A, Scherzer O (2008). *Iterative Regularization Methods for Nonlinear Ill-Posed Problems*. Vol. 6. Walter de Gruyter. https://doi.org/10.1515/9783110208276.

Kareiva PM (1983). *Local Movement in Herbivorous Insects: Applying a Passive Diffusion Model to Mark-Recapture Field Experiments*. *Oecologia* **57**, 322–327. JSTOR: 4216967.

Keller EF, Segel LA (1970). *Initiation of Slime Mold Aggregation Viewed as an Instability*. *Journal of Theoretical Biology* **26**, 399–415. https://doi.org/10.1016/0022-5193(70)90092-5.

1140

Keller EF, Segel LA (1971). *Model for Chemotaxis*. *Journal of Theoretical Biology* **30**, 225-234. https://doi.org/10.1016/0022-5193(71)90050-6.

Khadka R, Aydemir N, Carraher C, Hamiaux C, Colbert D, Cheema J, Malmström J, Kralicek A, Travas-Sejdic J (2019). *An Ultrasensitive Electrochemical Impedance-Based Biosensor Using Insect Odorant Receptors to Detect Odorants*. *Biosensors and Bioelectronics* **126**, 207-213. https://doi.org/10.1016/j.bios.2018.10.043. (Visited on 04/28/2024).

1145

Le Dimet FX, Talagrand O (1986). *Variational Algorithms for Analysis and Assimilation of Meteorological Observations: Theoretical Aspects*. *Tellus A: Dynamic Meteorology and Oceanography* **38**, 97-110. https://doi.org/10.3402/tellusa.v38i2.11706.

Leelössy Á, Molnár F, Izsák F, Havasi Á, Lagzi I, Mészáros R (2014). *Dispersion Modeling of Air Pollutants in the Atmosphere: A Review*. *Open Geosciences*. https://doi.org/doi:10.2478/s13533-012-0188-6.

1150

Lions JL (1971). *Optimal Control of Systems Governed by Partial Differential Equations*. Vol. 170. Springer Verlag.

Lu Y, Li H, Zhuang S, Zhang D, Zhang Q, Zhou J, Dong S, Liu Q, Wang P (2014). *Olfactory Biosensor Using Odorant-Binding Proteins from Honeybee: Ligands of Floral Odors and Pheromones Detection by Electrochemical Impedance*. *Sensors and Actuators B: Chemical* **193**, 420-427. https://doi.org/10.1016/j.snb.2013.11.045. (Visited on 04/28/2024).

1155

Lushi E, Stockie JM (2010). *An Inverse Gaussian Plume Approach for Estimating Atmospheric Pollutant Emissions from Multiple Point Sources*. *Atmospheric Environment* **44**, 1097-1107. https://doi.org/10.1016/j.atmosenv.2009.11.039. (Visited on 10/09/2023).

1160

Malou T, Labarthe S (2024a). *BI-DA-companion-code*. https://doi.org/10.5281/ZENODO.11175422. (Visited on 05/10/2024).

Malou T, Labarthe S (2024b). *Pherosensor-Toolbox*. [object Object]. https://doi.org/10.5281/ZENODO.11110047. (Visited on 05/10/2024).

1165

Martinez D, Arhidi L, Demondion E, Masson JB, Lucas P (2014). *Using Insect Electroantennogram Sensors on Autonomous Robots for Olfactory Searches*. *Journal of Visualized Experiments*, 51704. https://doi.org/10.3791/51704. (Visited on 04/28/2024).

McMahon T, Denison P (1979). *Empirical Atmospheric Deposition Parameters—a Survey*. *Atmospheric Environment* (1967) **13**, 571-585.

1170

McNeil JN (1991). *Behavioral Ecology of Pheromone-Mediated Communication in Moths and Its Importance in the Use of Pheromone Traps*. *Annual Review of Entomology* **36**, 407-430. https://doi.org/10.1146/annurev.en.36.010191.002203. (Visited on 07/29/2023).

Moitra P, Bhagat D, Pratap R, Bhattacharya S (2016). *A Novel Bio-Engineering Approach to Generate an Eminent Surface-Functionalized Template for Selective Detection of Female Sex Pheromone of Helicoverpa Armigera*. *Scientific Reports* **6**, 37355. https://doi.org/10.1038/srep37355. (Visited on 04/28/2024).

1175

Murray JD (2002). *Mathematical Biology: I. An Introduction*. Springer.
Gebbers R, Adamchuk VI (2010). *Precision Agriculture and Food Security*. Science **327**, 828-831.
https://doi.org/10.1126/science.1183899.

Genieys S, Volpert V, Auger P (2006). *Pattern and Waves for a Model in Population Dynamics with Nonlocal Consumption of Resources*. Mathematical Modelling of Natural Phenomena. https://doi.org/10.1051/mmnp:2006004.

IPPC Secretariat. (2021). *Scientific Review of the Impact of Climate Change on Plant Pests*. FAO on behalf of the IPPC Secretariat. https://doi.org/10.4060/cb4769en.

Isakov V (2017). *Inverse Problems for Partial Differential Equations*. Vol. 127. Applied Mathematical Sciences. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-51658-5.

Kaltenbacher B, Neubauer A, Scherzer O (2008). *Iterative Regularization Methods for Nonlinear Ill-Posed Problems*. Vol. 6. Walter de Gruyter. https://doi.org/10.1515/9783110208276.

Kareiva PM (1983). *Local Movement in Herbivorous Insects: Applying a Passive Diffusion Model to Mark-Recapture Field Experiments*. Oecologia **57**, 322-327. JSTOR: 4216967.

Keller EF, Segel LA (1970). *Initiation of Slime Mold Aggregation Viewed as an Instability*. Journal of Theoretical Biology **26**, 399-415. https://doi.org/10.1016/0022-5193(70)90092-5.

Keller EF, Segel LA (1971). *Model for Chemotaxis*. Journal of Theoretical Biology **30**, 225-234.
https://doi.org/10.1016/0022-5193(71)90050-6.

Khadka R, Aydemir N, Carraher C, Hamiaux C, Colbert D, Cheema J, Malmström J, Kralicek A, Travas-Sejdic J (2019). *An Ultrasensitive Electrochemical Impedance-Based Biosensor Using Insect Odorant Receptors to Detect Odorants*. Biosensors and Bioelectronics **126**, 207-213. https://doi.org/10.1016/j.bios.2018.10.043.

Le Dimet FX, Talagrand O (1986). *Variational Algorithms for Analysis and Assimilation of Meteorological Observations: Theoretical Aspects*. Tellus A: Dynamic Meteorology and Oceanography **38**, 97-110. https://doi.org/10.3402/tellusa.v38i2.11706.

Leelössy Á, Molnár F, Izsák F, Havasi Á, Lagzi I, Mészáros R (2014). *Dispersion Modeling of Air Pollutants in the Atmosphere: A Review*. Open Geosciences. https://doi.org/doi:10.2478/s13533-012-0188-6.

Lions JL (1971). *Optimal Control of Systems Governed by Partial Differential Equations*. Vol. 170. Springer Verlag.

Lu Y, Li H, Zhuang S, Zhang D, Zhang Q, Zhou J, Dong S, Liu Q, Wang P (2014). *Olfactory Biosensor Using Odorant-Binding Proteins from Honeybee: Ligands of Floral Odors and Pheromones Detection by Electrochemical Impedance*. Sensors and Actuators B: Chemical **193**, 420-427. https://doi.org/10.1016/j.snb.2013.11.045.

Lushi E, Stockie JM (2010). *An Inverse Gaussian Plume Approach for Estimating Atmospheric Pollutant Emissions from Multiple Point Sources*. Atmospheric Environment **44**, 1097-1107. https://doi.org/10.1016/j.atmosenv.2009.11.039.

Malou T, Labarthe S (2024a). *BI-DA-companion-code*. https://doi.org/10.5281/ZENODO.11175422.

Malou T, Labarthe S (2024b). *Pherosensor-Toolbox*. [object Object]. https://doi.org/10.5281/ZENODO.11110047.

Malou T, Labarthe S (2024c). *Pherosensor-toolbox: a Python package for Biology-Informed Data Assimilation*. Journal of Open Source Software **9**, 6863. https://doi.org/10.21105/joss.06863.

Martinez D, Arhidi L, Demondion E, Masson JB, Lucas P (2014). *Using Insect Electroantennogram Sensors on Autonomous Robots for Olfactory Searches*. Journal of Visualized Experiments, 51704.
https://doi.org/10.3791/51704.
1180

Myrick AJ, Baker TC (2011). Locating a Compact Odor Source Using a Four-Channel Insect Electroantennogram Sensor. Bioinspiration & Biomimetics **6**, 016002. https://doi.org/10.1088/1748-3182/6/1/016002. (Visited on 04/28/2024).

Negri RM, Bernik D (2008). Tracking the Sex Pheromone of Codling Moth against a Background of Host Volatiles with an Electronic Nose. Crop Protection **27**, 1295-1302. https://doi.org/10.1016/j.cropro.2008.02.008. (Visited on 04/28/2024).

Nesterov Y (2018). Lectures on Convex Optimization. Vol. 137. Springer.

1185

Nisbet RM, Gurney W (1982). Modelling Fluctuating Populations. John Wiley.

Oerke EC (2006). Crop Losses to Pests. The Journal of Agricultural Science **144**, 31-43. https://doi.org/10.1017/S0021859605005708. (Visited on 06/26/2023).

Okubo A (1986). Dynamical Aspects of Animal Grouping: Swarms, Schools, Flocks, and Herds. Advances in Biophysics **22**, 1-94. https://doi.org/10.1016/0065-227X(86)90003-1.

1190

Panofsky HA, Dutton JA (1984). Atmospheric Turbulence. Models and Methods for Engineering Applications. New York: Wiley.

Parikh N, Boyd S, et al. (2014). Proximal Algorithms. Foundations and trends® in Optimization **1**, 127-239.

1195

Parisey N, Leclerc M, Adamczyk-Chauvat K (2022). Optimal Spatial Monitoring of Populations Described by Reaction-Diffusion Models. Journal of Theoretical Biology **534**, 110976.

Pawson SM, Kerr JL, O'Connor BC, Lucas P, Martinez D, Allison JD, Strand TM (2020). Light-Weight Portable Electroantennography Device as a Future Field-Based Tool for Applied Chemical Ecology. Journal of Chemical Ecology **46**, 557-566. https://doi.org/10.1007/s10886-020-01190-6. (Visited on 07/28/2023).

1200

Perthame B (2006). Transport Equations in Biology. Birkhäuser Basel.

Pudykiewicz JA (1998). Application of Adjoint Tracer Transport Equations for Evaluating Source Parameters. Atmospheric Environment **32**, 3039-3050. https://doi.org/10.1016/S1352-2310(97)00480-9.

1205

Raissi M, Perdikaris P, Karniadakis G (2019). Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations. Journal of Computational Physics **378**, 686-707. https://doi.org/10.1016/j.jcp.2018.10.045.

Sasaki Y (1958). An Objective Analysis Based on the Variational Method. Journal of the Meteorological Society of Japan. Ser. II **36**, 77-88. https://doi.org/10.2151/jmsj1923.36.3\_77.

1210

Savary S, Willocquet L, Pethybridge SJ, Esker P, McRoberts N, Nelson A (2019). The Global Burden of Pathogens and Pests on Major Food Crops. Nature Ecology & Evolution **3**, 430-439. https://doi.org/10.1038/s41559-018-0793-y. (Visited on 06/26/2023).

Shigesada N (1980). Spatial Distribution of Dispersing Animals. Journal of mathematical biology **9**, 85-96.

1215

Shigesada N, Kawasaki K, Teramoto E (1979). Spatial Segregation of Interacting Species. Journal of Theoretical Biology **79**, 83-99. https://doi.org/10.1016/0022-5193(79)90258-3.

Skellam JG (1951). Random Dispersal in Theoretical Populations. Biometrika **38**, 196-218.

Sportisse B (2007). A Review of Current Issues in Air Pollution Modeling and Simulation. Computational Geosciences.

1220

Steffens C, Manzoli A, Oliveira JE, Leite FL, Correa DS, Herrmann PSP (2014). Bio-Inspired Sensor for Insect Pheromone Analysis Based on Polyaniline Functionalized AFM Cantilever Sensor. Sensors and Actuators B: Chemical **191**, 643-649. https://doi.org/10.1016/j.snb.2013.10.053. (Visited on 04/28/2024).
McMahon T, Denison P (1979). *Empirical Atmospheric Deposition Parameters—a Survey*. Atmospheric Environment (1967) **13**, 571-585. https://doi.org/10.1016/0004-6981(79)90186-0.

McNeil JN (1991). *Behavioral Ecology of Pheromone-Mediated Communication in Moths and Its Importance in the Use of Pheromone Traps*. Annual Review of Entomology **36**, 407-430. https://doi.org/10.1146/annurev.en.36.010191.002203.

Moitra P, Bhagat D, Pratap R, Bhattacharya S (2016). A Novel Bio-Engineering Approach to Generate an Eminent Surface-Functionalized Template for Selective Detection of Female Sex Pheromone of Helicoverpa Armigera. Scientific Reports **6**, 37355. https://doi.org/10.1038/srep37355.

Murray JD (2002). *Mathematical Biology: I. An Introduction*. Springer.

Myrick AJ, Baker TC (2011). Locating a Compact Odor Source Using a Four-Channel Insect Electroantennogram Sensor. Bioinspiration & Biomimetics **6**, 016002. https://doi.org/10.1088/1748-3182/6/1/016002.

Negri RM, Bernik D (2008). Tracking the Sex Pheromone of Codling Moth against a Background of Host Volatiles with an Electronic Nose. Crop Protection **27**, 1295-1302. https://doi.org/10.1016/j.cropro.2008.02.008.

Nesterov Y (2018). *Lectures on Convex Optimization*. Vol. 137. Springer.

Nisbet RM, Gurney W (1982). *Modelling Fluctuating Populations*. John Wiley.

Oerke EC (2006). Crop Losses to Pests. The Journal of Agricultural Science **144**, 31-43. https://doi.org/10.1017/S0021859605005708.

Okubo A (1986). *Dynamical Aspects of Animal Grouping: Swarms, Schools, Flocks, and Herds*. Advances in Biophysics **22**, 1-94. https://doi.org/10.1016/0065-227X(86)90003-1.

Panofsky HA, Dutton JA (1984). Atmospheric Turbulence. Models and Methods for Engineering Applications. New York: Wiley.

Parikh N, Boyd S, et al. (2014). *Proximal Algorithms*. Foundations and trends in Optimization **1**, 127-239. https://doi.org/10.1561/2400000003.

Parisey N, Leclerc M, Adamczyk-Chauvat K (2022). Optimal Spatial Monitoring of Populations Described by Reaction-Diffusion Models. Journal of Theoretical Biology **534**, 110976. https://doi.org/10.1016/j.jtbi.2021.110976.

Pawson SM, Kerr JL, O'Connor BC, Lucas P, Martinez D, Allison JD, Strand TM (2020). Light-Weight Portable Electroantennography Device as a Future Field-Based Tool for Applied Chemical Ecology. Journal of Chemical Ecology **46**, 557-566. https://doi.org/10.1007/s10886-020-01190-6.

Perthame B (2006). *Transport Equations in Biology*. Birkhäuser Basel.

Pudykiewicz JA (1998). Application of Adjoint Tracer Transport Equations for Evaluating Source Parameters. Atmospheric Environment **32**, 3039-3050. https://doi.org/10.1016/S1352-2310(97)00480-9.

Raissi M, Perdikaris P, Karniadakis G (2019). *Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations*. Journal of Computational Physics **378**, 686-707. https://doi.org/10.1016/j.jcp.2018.10.045.

Sasaki Y (1958). *An Objective Analysis Based on the Variational Method*. Journal of the Meteorological Society of Japan. Ser. II **36**, 77-88. https://doi.org/10.2151/jmsj1923.36.3\_77.

Savary S, Willocquet L, Pethybridge SJ, Esker P, McRoberts N, Nelson A (2019). *The Global Burden of Pathogens and Pests on Major Food Crops*. Nature Ecology & Evolution **3**, 430-439. https://doi.org/10.1038/s41559-018-0793-y.
1225

Stehman SV (1997). *Selecting and Interpreting Measures of Thematic Classification Accuracy*. *Remote Sensing of Environment* **62**, 77-89. https://doi.org/10.1016/S0034-4257(97)00083-7.

Stockie JM (2011). *The Mathematics of Atmospheric Dispersion Modeling*. SIAM Review **53**, 349-372. https://doi.org/10.1137/10080991X. (Visited on 10/09/2023).

1230

Sun GQ (2016). *Mathematical Modeling of Population Dynamics with Allee Effect*. *Nonlinear Dynamics*. https://doi.org/10.1007/s11071-016-2671-y.

Takahashi LT, Maidana NA, Ferreira WC, Pulino P, Yang HM (2005). *Mathematical Models for the Aedes Aegypti Dispersal Dynamics: Travelling Waves by Wing and Wind*. *Bulletin of mathematical Biology* **67**, 509-528. https://doi.org/10.1016/j.bulm.2004.08.005.

1235

Tarantola A (2005). *Inverse Problem Theory and Methods for Model Parameter Estimation*. Vol. 89. SIAM.

Tibshirani R (1996). *Regression Shrinkage and Selection via the Lasso*. *Journal of the Royal Statistical Society: Series B (Methodological)* **58**, 267-288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x. eprint: https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1996.tb02080.x.

1240

Trémolet Y (2007). *Model-Error Estimation in 4D-Var*. *Quarterly Journal of the Royal Meteorological Society* **133**, 1267-1280. https://doi.org/10.1002/qj.94.

Tumlinson JH, Mitchell ER, Teal PEA, Heath RR, Mengelkoch LJ (1986). *Sex Pheromone of Fall Armyworm, Spodoptera Frugiperda (J.E. Smith): Identification of Components Critical toattraction in the Field*. *Journal of Chemical Ecology* **12**, 1909-1926. https://doi.org/10.1007/BF01041855. (Visited on 04/28/2024).

1245

Wehrenfennig C, Schott M, Gasch T, Sauerwald T, Düring RA, Vilcinskas A, Kohl CD (2012). *Laboratory Characterization of Metal-oxide Sensors Intended for in Situ Analyses of Pheromones - SOMMSA Approach*. *physica status solidi* (a) **209**, 935-939. https://doi.org/10.1002/pssa.201100586. (Visited on 04/28/2024).

1250

Yazdani A, Lu L, Raissi M, Karniadakis GE (2020). *Systems Biology Informed Deep Learning for Inferring Parameters and Hidden Dynamics*. *PLOS Computational Biology* **16**. Ed. by Vassily Hatzimanikatis, e1007575. https://doi.org/10.1371/journal.pcbi.1007575. (Visited on 10/09/2023).

## Appendix A. Data, script, code, and supplementary information availability

1255

The toolbox containing the core functions for BI-DA is available at https://forgemia.inra.fr/pherosensor/pherosensor-toolbox and (Malou et al., 2024b). The companion codes needed to reproduce the different computation and figures of this paper are available at https://forgemia.inra.fr/pherosensor/companion-code-bi-da and (Malou et al., 2024a).

## Appendix B. Supplementary material

1260

### B.1. Derivation of the depth-integrated chemical transport model

This section aims to present the derivation of the 2D model of pheromone propagation (1). In a first part, a classical 3D model of the pheromone propagation in the atmosphere, also called 3D CTM (Chemical Transport Model), is presented. Then, the 2D model, presented in Sec. 2.1 and used in this study, is derived from this 3D model. In the last part, the calibration of the 2D model using the 3D model is presented.

1265
Shigesada N (1980). *Spatial Distribution of Dispersing Animals*. Journal of mathematical biology **9**, 85-96. https://doi.org/10.1007/BF00276037.

Shigesada N, Kawasaki K, Teramoto E (1979). *Spatial Segregation of Interacting Species*. Journal of Theoretical Biology **79**, 83-99. https://doi.org/10.1016/0022-5193(79)90258-3.

Skellam JG (1951). *Random Dispersal in Theoretical Populations*. Biometrika **38**, 196-218. https://doi.org/10.2307/2332328.

Sportisse B (2007). *A Review of Current Issues in Air Pollution Modeling and Simulation*. Computational Geosciences **11**, 159-181. https://doi.org/10.1007/s10596-006-9036-4.

Steffens C, Manzoli A, Oliveira JE, Leite FL, Correa DS, Herrmann PSP (2014). *Bio-Inspired Sensor for Insect Pheromone Analysis Based on Polyaniline Functionalized AFM Cantilever Sensor*. Sensors and Actuators B: Chemical **191**, 643-649. https://doi.org/10.1016/j.snb.2013.10.053.

Stehman SV (1997). *Selecting and Interpreting Measures of Thematic Classification Accuracy*. Remote Sensing of Environment **62**, 77-89. https://doi.org/10.1016/S0034-4257(97)00083-7.

Stockie JM (2011). *The Mathematics of Atmospheric Dispersion Modeling*. SIAM Review **53**, 349-372. https://doi.org/10.1137/10080991X.

Sun GQ (2016). *Mathematical Modeling of Population Dynamics with Allee Effect*. Nonlinear Dynamics. https://doi.org/10.1007/s11071-016-2671-y.

Takahashi LT, Maidana NA, Ferreira WC, Pulino P, Yang HM (2005). *Mathematical Models for the Aedes Aegypti Dispersal Dynamics: Travelling Waves by Wing and Wind*. Bulletin of mathematical Biology **67**, 509-528. https://doi.org/10.1016/j.bulm.2004.08.005.

Tarantola A (2005). *Inverse Problem Theory and Methods for Model Parameter Estimation*. Vol. 89. SIAM.

Tibshirani R (1996). *Regression Shrinkage and Selection via the Lasso*. Journal of the Royal Statistical Society: Series B (Methodological) **58**, 267-288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x. eprint: https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1996.tb02080.x.

Trémolet Y (2007). *Model-Error Estimation in 4D-Var*. Quarterly Journal of the Royal Meteorological Society **133**, 1267-1280. https://doi.org/10.1002/qj.94.

Tumlinson JH, Mitchell ER, Teal PEA, Heath RR, Mengelkoch LJ (1986). *Sex Pheromone of Fall Armyworm, Spodoptera Frugiperda (J.E. Smith): Identification of Components Critical to Attraction in the Field*. Journal of Chemical Ecology **12**, 1909-1926. https://doi.org/10.1007/BF01041855.

Wehrenfennig C, Schott M, Gasch T, Sauerwald T, Düring RA, Vilcinskas A, Kohl CD (2012). *Laboratory Characterization of Metal-oxide Sensors Intended for in Situ Analyses of Pheromones - SOMMSA Approach*. physica status solidi (a) **209**, 935-939. https://doi.org/10.1002/pssa.201100586.

Yazdani A, Lu L, Raissi M, Karniadakis GE (2020). *Systems Biology Informed Deep Learning for Inferring Parameters and Hidden Dynamics*. PLOS Computational Biology **16**. Ed. by Vassily Hatzi-manikatis, e1007575. https://doi.org/10.1371/journal.pcbi.1007575.

# Appendix A. Supplementary material

## A.1. Derivation of the depth-integrated chemical transport model

This section aims to present the derivation of the 2D model of pheromone propagation (1). In a first part, a classical 3D model of the pheromone propagation in the atmosphere, also called
### B.1.1. 3D chemical transport model.
In this section, $C(\vec{x}, t) (g.m^{-3})$ denotes the density of chemical (here the pheromone) at a given point $\vec{x} = (x, y, z) \in \mathbb{R}^2 \times \mathbb{R}^+$ and at a time $t \in \mathbb{R}^+$. The equation of mass conservation of $C$ can be written as the following PDE, see e.g. reviews in (Leelössy et al., 2014; Sportisse, 2007):

$$ (34) \qquad \frac{\partial C}{\partial t}(\vec{x}, t) + \nabla_{3D} \cdot (-\mathbf{K}_{3D}(\vec{x}, t)\nabla_{3D}C(\vec{x}, t) + \vec{u}_{3D}(\vec{x}, t)C(\vec{x}, t)) = S(\vec{x}, t) $$

with $-\mathbf{K}_{3D}\nabla_{3D}C + \vec{u}_{3D}C (g.m^{-3}.m.s^{-1})$ the mass flux of pheromone due to the diffusion and the advection of the pheromone in the atmosphere, and $S(\vec{x}, t) (g.m^{-3}.s^{-1})$ the source term, i.e. the quantity of pheromone emitted by the insects per seconds and per air volume in the present context. In this section, $\nabla_{3D}$ denotes the 3D gradient operator $\nabla_{3D} = (\partial_x \ \partial_y \ \partial_z)$, in order to distinguish with the 2D gradient operator used in the rest of the article.

This equation can be obtained from the reactant Navier-Stokes equations using the dilution hypothesis. This hypothesis consists in assuming that the chemical of interest is present in the atmosphere in a very small quantity and does not significantly impact the dynamic of the atmosphere. It results that the dynamics of the atmosphere and the dynamics of the pheromone are decoupled.

As classically done in the K-theory, see e.g. Panofsky et al., 1984; Sportisse, 2007, the diffusion flux due to fluctuations of the air motion caused by the turbulence of the atmosphere at non-resolved scales is closed using the Fick's law. This law relates the diffusion flux to the gradient of the concentration and leads to the diffusion term $-\mathbf{K}_{3D}\nabla_{3D}C$ with $\mathbf{K}_{3D}(\vec{x}, t)$ the diffusion tensor. In the present case, the diffusion tensor, a priori non-homogeneous, is assumed anisotropic in the $Oxy$ plane but isotropic along the $z$ axis. Thus, the diffusion tensor can be

$$ \text{written as: } \mathbf{K}_{3D} = \begin{pmatrix} K_{xx} & K_{xy} & 0 \\ K_{yx} & K_{yy} & 0 \\ 0 & 0 & K_{zz} \end{pmatrix}. $$

The advection flux of mass $\vec{u}_{3D}C$ is due to the wind and the gravity, and $\vec{u}_{3D}(\vec{x}, t)$ is the velocity of the pheromones in the atmosphere. Thus, the velocity of the pheromones can be written as $\vec{u}_{3D} = (u, v, -w_{set})$. The components $u$ and $v$ of the velocity are given by the wind velocity. Since the vertical wind velocity is assumed to be negligible and the pheromones are heavier than the air, the vertical component $-w_{set}$ is the settling velocity, i.e. the velocity at which the pheromones fall due to the gravity.

The dilution hypothesis and decoupling of the pheromone and atmosphere dynamics imply that the diffusion tensor $\mathbf{K}_{3D}$ and the wind field $\vec{u}_{3D}$ can be computed separately by a meteorological model and given as environmental parameters of the CTM model (34).

As the pheromones are settled under the effect of the gravity, they end up on the ground and a part of the pheromones are absorbed by the ground. To model this phenomenon, the following boundary condition on the vertical mass flux is commonly used:

$$ (35) \qquad (K_{zz} \partial_z C + w_{set} C)|_{z=0} = w_{dep} C|_{z=0} $$

with $w_{dep}$ the deposition velocity.

### B.1.2. Depth integration of the 3D CTM model.
The 3D CTM (34) is now integrated over a layer of the atmosphere $z \in [h_{-}; h_{+}]$ in order to obtain the 2D model used and presented in Sec. 2.1:

$$ (36) \qquad \int_{h_{-}}^{h_{+}} \frac{\partial C}{\partial t} - \nabla_{3D} \cdot (\mathbf{K}_{3D} \nabla_{3D} C) + \nabla_{3D} \cdot (\vec{u}_{3D} C) dz = \int_{h_{-}}^{h_{+}} S dz $$

In the following, the total concentration of pheromone in the layer is denoted $c(x, y, t) = \int_{h_{-}}^{h_{+}} C(x, y, z, t) dz$. The integral of the different terms are now computed in order to obtain a 2D model.
3D CTM (Chemical Transport Model), is presented. Then, the 2D model, presented in Sec. 2.1 and used in this study, is derived from this 3D model. In the last part, the calibration of the 2D model using the 3D model is presented.

A.1.1. 3D chemical transport model. In this section, $C(\vec{x}, t) (g.m^{-3})$ denotes the density of chemical (here the pheromone) at a given point $\vec{x} = (x, y, z) \in \mathbb{R}^2 \times \mathbb{R}^+$ and at a time $t \in \mathbb{R}^+$. The equation of mass conservation of $C$ can be written as the following PDE, see e.g. reviews in (Leelössy et al., 2014; Sportisse, 2007):

$$ (34) \qquad \frac{\partial C}{\partial t}(\vec{x}, t) + \nabla_{3D} \cdot (-\mathbf{K}_{3D}(\vec{x}, t)\nabla_{3D}C(\vec{x}, t) + \vec{u}_{3D}(\vec{x}, t)C(\vec{x}, t)) = S(\vec{x}, t) $$

with $-\mathbf{K}_{3D}\nabla_{3D}C + \vec{u}_{3D}C(g.m^{-3}.m.s^{-1})$ the mass flux of pheromone due to the diffusion and the advection of the pheromone in the atmosphere, and $S(\vec{x}, t) (g.m^{-3}.s^{-1})$ the source term, i.e. the quantity of pheromone emitted by the insects per seconds and per air volume in the present context. In this section, $\nabla_{3D}$ denotes the 3D gradient operator $\nabla_{3D} = (\partial_x \ \partial_y \ \partial_z)$, in order to distinguish with the 2D gradient operator used in the rest of the article.

This equation can be obtained from the reactant Navier-Stokes equations using the dilution hypothesis. This hypothesis consists in assuming that the chemical of interest is present in the atmosphere in a very small quantity and does not significantly impact the dynamic of the atmosphere. It results that the dynamics of the atmosphere and the dynamics of the pheromone are decoupled.

As classically done in the K-theory, see e.g. Panofsky et al., 1984; Sportisse, 2007, the diffusion flux due to fluctuations of the air motion caused by the turbulence of the atmosphere at non-resolved scales is closed using the Fick's law. This law relates the diffusion flux to the gradient of the concentration and leads to the diffusion term $-\mathbf{K}_{3D}\nabla_{3D}C$ with $\mathbf{K}_{3D}(\vec{x}, t)$ the diffusion tensor. In the present case, the diffusion tensor, a priori non-homogeneous, is assumed anisotropic in the $Oxy$ plane but isotropic along the $z$ axis. Thus, the diffusion tensor can be

$$ \text{written as: } \mathbf{K}_{3D} = \begin{pmatrix} K_{xx} & K_{xy} & 0 \\ K_{yx} & K_{yy} & 0 \\ 0 & 0 & K_{zz} \end{pmatrix}. $$

The advection flux of mass $\vec{u}_{3D}C$ is due to the wind and the gravity, and $\vec{u}_{3D}(\vec{x}, t)$ is the velocity of the pheromones in the atmosphere. Thus, the velocity of the pheromones can be written as $\vec{u}_{3D} = (u, v, -w_{set})$. The components $u$ and $v$ of the velocity are given by the wind velocity. Since the vertical wind velocity is assumed to be negligible and the pheromones are heavier than the air, the vertical component $-w_{set}$ is the settling velocity, i.e. the velocity at which the pheromones fall due to the gravity.

The dilution hypothesis and decoupling of the pheromone and atmosphere dynamics imply that the diffusion tensor $\mathbf{K}_{3D}$ and the wind field $\vec{u}_{3D}$ can be computed separately by a meteorological model and given as environmental parameters of the CTM model (34).

As the pheromones are settled under the effect of the gravity, they end up on the ground and a part of the pheromones are absorbed by the ground. To model this phenomenon, the following boundary condition on the vertical mass flux is commonly used:

$$ (35) \qquad (K_{zz}\partial_z C + w_{set}C)|_{z=0} = w_{dep}C|_{z=0} $$

with $w_{dep}$ the deposition velocity.
The integral of the time derivative becomes simply:

$$ (37) \qquad \int_{h_{-}}^{h_{+}} \partial_{t} C d z=\partial_{t} c $$

The integral of the advection term can be written as:

$$ \int_{h_{-}}^{h_{+}} \nabla_{3D} \cdot (\vec{u}_{3D} C) d z=\underbrace{\int_{h_{-}}^{h_{+}} \partial_{x}(u C)+\partial_{y}(v C) d z}_{=l_{1}}+\underbrace{\int_{h_{-}}^{h_{+}} \partial_{z}(-w_{\text {set }} C) d z}_{=l_{2}} $$

Assuming that the velocity components u and v are constant or near the average value over the layer of interest, the integral $l_1$ becomes:

$$ (38) \qquad l_1 = \partial_x \left( \int_{h_{-}}^{h_{+}} u C dz \right) + \partial_y \left( \int_{h_{-}}^{h_{+}} v C dz \right) = \partial_x (u c) + \partial_y (v c) $$

The integral of the diffusion term can be written as:

$$ -\int_{h_{-}}^{h_{+}} \nabla_{3D} \cdot (\mathbf{K}_{3D} \nabla_{3D} C) d z=-\underbrace{\int_{h_{-}}^{h_{+}} \partial_{x}(K_{xx} \partial_{x} C+K_{xy} \partial_{y} C)+\partial_{y}(K_{yx} \partial_{x} C+K_{yy} \partial_{y} C) d z}_{=l_{3}}-\underbrace{\int_{h_{-}}^{h_{+}} \partial_{z}(K_{zz} \partial_{z} C) d z}_{=l_{4}} $$

In the present case, the diffusion coefficients in the plane Oxy $K_{xx}$, $K_{xy}$, $K_{yx}$ and $K_{yy}$ is assumed to not depend on the depth z which is a common assumption. Thus, the integral $l_3$ becomes:

$$ (39) \qquad l_3 = \partial_x \left( \int_{h_{-}}^{h_{+}} K_{xx} \partial_x C + K_{xy} \partial_y C dz \right) + \partial_y \left( \int_{h_{-}}^{h_{+}} K_{yx} \partial_x C + K_{yy} \partial_y C dz \right) \\ = \partial_x (K_{xx} \partial_x c + K_{xy} \partial_y c) + \partial_y (K_{yx} \partial_x c + K_{yy} \partial_y c) $$

From the integrals of the advection and the diffusion terms, it remains the vertical mass fluxes $l_2 - l_4$ that become:

$$ (40) \qquad l_2 - l_4 = -\int_{h_{-}}^{h_{+}} \partial_z (K_{zz} \partial_z C + w_{\text{set}} C) dz = -[K_{zz} \partial_z C + w_{\text{set}} C]_{z=h_{-}}^{z=h_{+}} $$

In the present study, these vertical fluxes, which model the quantity of pheromones leaving the layer of interest, are closed using a 1st order term, i.e. by a reaction term weighted by a loss rate $\tau_{loss}(x, y)$:

$$ (41) \qquad l_2 - l_4 = -[K_{zz} \partial_z C + w_{\text{set}} C]_{z=h_{-}}^{z=h_{+}} \approx \tau_{\text{loss}} c $$

At last, the total amount of pheromone emitted by the insects in the layer of interest is denoted:

$$ (42) \qquad \int_{h_{-}}^{h_{+}} S(x, y, z, t) d z=s(x, y, t) $$

By introducing the expression of the integrals (37), (38), (39), (41) and (42) in the 3D CTM integrated over a layer of the atmosphere (36), we obtain the 2D CTM:

$$ (43) \qquad \frac{\partial c}{\partial t}-\nabla \cdot(\mathbf{K} \nabla c)+\nabla \cdot(\vec{u} c)+\tau_{\text {loss }} c=s $$

with the $\vec{u}$ and $\mathbf{K}$ resp. the wind field and diffusion tensor in the Oxy plane.

### B.1.3. Calibration of the loss coefficient $\tau_{loss}$. 
In the derivation of the 2D CTM from the 3D CTM, the vertical fluxes are closed by a reaction term weighted by a loss coefficient $\tau_{loss}$. This coefficient depends on the vertical fluxes. Thus, it depends on the atmosphere (through the vertical diffusion), on the pheromone's molecules (through the settling velocity) and on the vegetation (through the deposition velocity). However, the estimation of $\tau_{loss}$ is not straightforward. In this section, a method of estimation of $\tau_{loss}$ is proposed. In the present study, $\tau_{loss}$ is estimated by
A.1.2. *Depth integration of the 3D CTM model.* The 3D CTM (34) is now integrated over a layer of the atmosphere $z \in [h_{-}; h_{+}]$ in order to obtain the 2D model used and presented in Sec. 2.1:

$$ (36) \qquad \int_{h_{-}}^{h_{+}} \frac{\partial C}{\partial t} - \nabla_{3D} \cdot (\mathbf{K}_{3D} \nabla_{3D} C) + \nabla_{3D} \cdot (\vec{u}_{3D} C) \, dz = \int_{h_{-}}^{h_{+}} S \, dz $$

In the following, the total concentration of pheromone in the layer is denoted $c(x, y, t) = \int_{h_{-}}^{h_{+}} C(x, y, z, t) dz$. The integral of the different terms are now computed in order to obtain a 2D model.

The integral of the time derivative becomes simply:

$$ (37) \qquad \int_{h_{-}}^{h_{+}} \partial_t C dz = \partial_t c $$

The integral of the advection term can be written as:

$$ \int_{h_{-}}^{h_{+}} \nabla_{3D} \cdot (\vec{u}_{3D} C) dz = \underbrace{\int_{h_{-}}^{h_{+}} \partial_x (uC) + \partial_y (vC) dz}_{=l_1} + \underbrace{\int_{h_{-}}^{h_{+}} \partial_z (-w_{\text{set}} C) dz}_{=l_2} $$

Assuming that the velocity components $u$ and $v$ are constant or near the average value over the layer of interest, the integral $l_1$ becomes:

$$ (38) \qquad l_1 = \partial_x \left( \int_{h_{-}}^{h_{+}} u C dz \right) + \partial_y \left( \int_{h_{-}}^{h_{+}} v C dz \right) = \partial_x (u c) + \partial_y (v c) $$

The integral of the diffusion term can be written as:

$$ - \int_{h_{-}}^{h_{+}} \nabla_{3D} \cdot (\mathbf{K}_{3D} \nabla_{3D} C) \, dz = - \underbrace{\int_{h_{-}}^{h_{+}} \partial_x (K_{xx} \partial_x C + K_{xy} \partial_y C) + \partial_y (K_{yx} \partial_x C + K_{yy} \partial_y C) \, dz}_{=l_3} - \underbrace{\int_{h_{-}}^{h_{+}} \partial_z (K_{zz} \partial_z C) \, dz}_{=l_4} $$

In the present case, the diffusion coefficients in the plane $Oxy$ $K_{xx}$, $K_{xy}$, $K_{yx}$ and $K_{yy}$ is assumed to not depend on the depth $z$ which is a common assumption. Thus, the integral $l_3$ becomes:

$$ (39) \qquad \begin{aligned} l_3 &= \partial_x \left( \int_{h_{-}}^{h_{+}} K_{xx} \partial_x C + K_{xy} \partial_y C dz \right) + \partial_y \left( \int_{h_{-}}^{h_{+}} K_{yx} \partial_x C + K_{yy} \partial_y C dz \right) \\ &= \partial_x (K_{xx} \partial_x c + K_{xy} \partial_y c) + \partial_y (K_{yx} \partial_x c + K_{yy} \partial_y c) \end{aligned} $$

From the integrals of the advection and the diffusion terms, it remains the vertical mass fluxes $l_2 - l_4$ that become:

$$ (40) \qquad l_2 - l_4 = - \int_{h_{-}}^{h_{+}} \partial_z (K_{zz} \partial_z C + w_{\text{set}} C) \, dz = - [K_{zz} \partial_z C + w_{\text{set}} C]_{z=h_{-}}^{z=h_{+}} $$

In the present study, these vertical fluxes, which model the quantity of pheromones leaving the layer of interest, are closed using a 1st order term, i.e. by a reaction term weighted by a loss rate $\tau_{\text{loss}}(x, y)$:

$$ (41) \qquad l_2 - l_4 = - [K_{zz} \partial_z C + w_{\text{set}} C]_{z=h_{-}}^{z=h_{+}} \approx \tau_{\text{loss}} c $$

At last, the total amount of pheromone emitted by the insects in the layer of interest is denoted:

$$ (42) \qquad \int_{h_{-}}^{h_{+}} S(x, y, z, t) dz = s(x, y, t) $$

By introducing the expression of the integrals (37), (38), (39), (41) and (42) in the 3D CTM integrated over a layer of the atmosphere (36), we obtain the 2D CTM:

$$ (43) \qquad \frac{\partial c}{\partial t} - \nabla \cdot (\mathbf{K} \nabla c) + \nabla \cdot (\vec{u} c) + \tau_{\text{loss}} c = s $$
comparing the solution of the 2D CTM derived in the previous section and the integral over the layer of the Gaussian plume, the solution of the 3D CTM (34) in a particular case.

The following assumptions are made. It is assumed that the diffusion tensor is isotropic and constant ($K_{3D} = K_{ld}$) and that the wind is unidirectional and constant ($\vec{u}_{3D} = (u, 0, -w_{set})$). Moreover, it is assumed that the source term is of the form $S(\vec{x}) = Q\delta(x)\delta(y)\delta(z - H)$ with Q a constant emission rate and H the height at which the source of emission is located. The Gaussian plume is the solution of the 3D CTM under these assumptions and is of the following form, see more details and a review on the topic in (Stockie, 2011):

$$ (44) \qquad C_{gp}(r, y, z) = \frac{Q}{u} \frac{1}{\sqrt{4\pi r}} \mathcal{G}_h(r, y, z) \mathcal{G}_v(r, y, z) \mathcal{F}_{bc}(r, y, z) $$

with:

$$ \mathcal{G}_h(r, y, z) = \exp\left(-\frac{y^2}{4r}\right), $$

$$ \mathcal{G}_v(r, y, z) = \exp\left(-\frac{w_{set}(z-H)}{2K} - \frac{w_{set}^2 r}{4K^2}\right) \left(\exp\left(-\frac{(z-H)^2}{4r}\right) + \exp\left(-\frac{(z+H)^2}{4r}\right)\right), $$

$$ \mathcal{F}_{bc}(r, y, z) = \frac{4w_0\pi r}{K} \exp\left(-\frac{w_{set}(z-H)}{2K} - \frac{w_{set}^2 r}{4K^2}\right) \exp\left(\frac{w_0(z+H)}{K} + \frac{w_0^2 r}{K^2}\right) \text{erfc}\left(\frac{z+H}{2\sqrt{r}} + \frac{w_0\sqrt{r}}{K}\right), $$

the change of variable $r = \frac{w_{set}}{2}$.

In order to get realistic estimations of the loss coefficient $\tau_{loss}$, the solution of the 2D CTM (43) is empirically compared with the integral over the layer of the Gaussian plume $\int_{h-}^{h+} C_{gp}(r, y, z) dz$ under the assumptions presented above in a realistic case.

Based on the average of the wind speed of the meteorological data, the value $u = 12m.s^{-1}$ is used in this section. The settling velocity is set as $w_{set} = 2.7 \times 10^{-3}m.s^{-1}$ as used in (Stockie, 2011). It is assumed that the source of pheromone is located at $H = 2m$ and the layer of interest is $2m$ wide centred at the height of the pheromone source ($h_{\pm} = H \pm 1m$). The values of deposition velocity $w_{dep}$ with respect to the vegetation used here are found in (McMahon et al., 1979). From this realistic case, it results the value of $\tau_{loss}$ given in Tab. 1.

## B.2. List of additional population dynamics-informed regularization terms

We detail here additional regularization terms that could be used during the inference process. These regularization terms are not included in the results, since they were not included in the numerical exploration of the BI-DA problems, in order to limit the number of possibilities and to focus on the known behavior of the FallArmy Worm (FAW), our insect of interest. However, we indicate here a large list of features that could be introduced in the population dynamic-informed regularization term (27). The corresponding gradient of the terms are indicated in section B.3.2.

### B.2.1. Logarithmic barrier term for exclusion zone.
As mentioned in Sec. 3.1.1, the biology-based Tikhonov enables to take into account the favorite habitat of the insects. Nonetheless, this approach does not preclude the possibility of locating insects in areas that they will avoid.

To do so, one can consider a $j_{reg}$ based on a logarithmic barrier function, see e.g. (Nesterov, 2018). This approach, and barrier function in general, is in particular used in penalty method that replaces an inequality constraint (e.g. $h(s; x, y, t) < b \forall t \in [0; T] \forall (x, y) \in \Omega_c$) by a penalizing term (e.g. $\int_0^T \int_{\Omega_c} -\epsilon \log(b - h(s; x, y, t)) dxdydt$ with $\Omega_c$ the constraint domain $\epsilon$ a very small constant) easier to handle in optimization problems.

In the present case, this penalizing term aims to preclude locating insects in areas that they will avoid or where they will not emit pheromones, i.e. $s = 0$ or very small on the exclusion zone $\Omega_{exc}$ derived from the biological knowledge of the insect. For instance, it is known that the FAW
with the $\vec{u}$ and $K$ resp. the wind field and diffusion tensor in the $Oxy$ plane.

### A.1.3. Calibration of the loss coefficient $\tau_{loss}$.
In the derivation of the 2D CTM from the 3D CTM, the vertical fluxes are closed by a reaction term weighted by a loss coefficient $\tau_{loss}$. This coefficient depends on the vertical fluxes. Thus, it depends on the atmosphere (through the vertical diffusion), on the pheromone's molecules (through the settling velocity) and on the vegetation (through the deposition velocity). However, the estimation of $\tau_{loss}$ is not straightforward. In this section, a method of estimation of $\tau_{loss}$ is proposed. In the present study, $\tau_{loss}$ is estimated by comparing the solution of the 2D CTM derived in the previous section and the integral over the layer of the Gaussian plume, the solution of the 3D CTM (34) in a particular case.

The following assumptions are made. It is assumed that the diffusion tensor is isotropic and constant ($K_{3D} = K_{Id}$) and that the wind is unidirectional and constant ($\vec{u}_{3D} = (u, 0, -w_{set})$). Moreover, it is assumed that the source term is of the form $S(\vec{x}) = Q\delta(x)\delta(y)\delta(z - H)$ with $Q$ a constant emission rate and $H$ the height at which the source of emission is located. The Gaussian plume is the solution of the 3D CTM under these assumptions and is of the following form, see more details and a review on the topic in (Stockie, 2011):

$$ (44) \qquad C_{gp}(r, y, z) = \frac{Q}{u} \frac{1}{\sqrt{4\pi r}} \mathcal{G}_h(r, y, z) \mathcal{G}_v(r, y, z) \mathcal{F}_{bc}(r, y, z) $$

with:

$$ \mathcal{G}_h(r, y, z) = \exp\left(-\frac{y^2}{4r}\right), $$

$$ \mathcal{G}_v(r, y, z) = \exp\left(-\frac{w_{set}(z-H)}{2K} - \frac{w_{set}^2 r}{4K^2}\right) \left(\exp\left(-\frac{(z-H)^2}{4r}\right) + \exp\left(-\frac{(z+H)^2}{4r}\right)\right), $$

$$ \mathcal{F}_{bc}(r, y, z) = \frac{4w_0\pi r}{K} \exp\left(-\frac{w_{set}(z-H)}{2K} - \frac{w_{set}^2 r}{4K^2}\right) \exp\left(\frac{w_0(z+H)}{K} + \frac{w_0^2 r}{K^2}\right) \text{erfc}\left(\frac{z+H}{2\sqrt{r}} + \frac{w_0\sqrt{r}}{K}\right), $$

the change of variable $r = \frac{1}{u} \int_0^x K(\xi) d\xi$, erfc the complementary error function and $w_0 = w_{dep} - \frac{w_{set}}{2}$.

In order to get realistic estimations of the loss coefficient $\tau_{loss}$, the solution of the 2D CTM (43) $c$ is empirically compared with the integral over the layer of the Gaussian plume $\int_{h-}^{h+} C_{gp}(r, y, z) dz$ under the assumptions presented above in a realistic case.

Based on the average of the wind speed of the meteorological data, the value $u = 12m.s^{-1}$ is used in this section. The settling velocity is set as $w_{set} = 2.7 \times 10^{-3}m.s^{-1}$ as used in (Stockie, 2011). It is assumed that the source of pheromone is located at $H = 2m$ and the layer of interest is $2m$ wide centred at the height of the pheromone source ($h_{\pm} = H \pm 1m$). The values of deposition velocity $w_{dep}$ with respect to the vegetation used here are found in (McMahon et al., 1979). From this realistic case, it results the value of $\tau_{loss}$ given in Tab. 1.

## A.2. List of additional population dynamics-informed regularization terms

We detail here additional regularization terms that could be used during the inference process. These regularization terms are not included in the results, since they were not included in the numerical exploration of the BI-DA problems, in order to limit the number of possibilities and to focus on the known behavior of the FallArmy Worm (FAW), our insect of interest. However, we indicate here a large list of features that could be introduced in the population dynamic-informed regularization term (27). The corresponding gradient of the terms are indicated in section A.3.2.

### A.2.1. Logarithmic barrier term for exclusion zone.
As mentioned in Sec. 3.1.1, the biology-based Tikhonov enables to take into account the favorite habitat of the insects. Nonetheless, this approach does not preclude the possibility of locating insects in areas that they will avoid.
females emit pheromones in a vertical position on vegetation that are favorable for the growth of the larva. Hence, the FAW females will not emit pheromones in some areas that are not suitable such as roads.

Therefore, the following regularization/penalization term is considered:

$$ (45) \qquad j_{\text{reg,exc}}(s) = \int_0^T \int_{\Omega_{\text{exc}}} -\epsilon \log(\epsilon - s(x, y, t)^2) dxdydt $$

Let us note that this logarithmic barrier function is smooth (and hence differentiable) and tends to the indicator function (0 if $s = 0$ locally, $+\infty$ otherwise) when $\epsilon \to 0$.

### B.2.2. Population-dynamics models.
Recall that the quantity of pheromone emitted $s$ can be written as $s = pq$ with $q$ the quantity of pheromone emitted per insect, which is assumed to be a biological knowledge, and $p$ the density of insect. Moreover, the density $p$, and consequently $s$, follows a population dynamic model that can be written in the generic way as following:

$$ (46) \qquad \begin{aligned} & \partial_t p + \nabla \left( \sum_i F_i(p) \right) - \sum_i R_i(p) = 0 \\ & \Leftrightarrow \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 \end{aligned} $$

with $F_i$ and $R_i$ are resp. flux and reaction terms that model the behavior of the insect specie of interest, and $F_i(p = \frac{s}{q}) = \mathcal{F}_i(s)$ and $R_i(p = \frac{s}{q}) = \mathcal{R}_i(s)$.

As mentioned earlier, only a little number of terms have been introduced in this study. In this section, we aim to indicate a large list of features and resulting population dynamics terms that could be introduced in the population dynamic-informed regularization term (27).

**Diffusion term.** The most classical term modelling population dynamics is the diffusion term:

$$ (47) \qquad F_d(p; x, y, t) = -\sigma(x, y) \nabla p(x, y, t) \Leftrightarrow \mathcal{F}_d(s) = -\frac{1}{q(t)} \sigma(x, y) \nabla s(x, y, t) $$

with $\sigma [m^2 \cdot s^{-1}]$ the diffusion coefficient of the population, depending notably on the specie.

This term models the spread of the population at a macroscopic scale caused by a motion of each individual in a random direction with a probability independent of other individuals or the environment. This diffusion term can be understood as a continuous and large population limit of a discrete microscopic unbiased random walk process with a fixed probability of jump, see e.g. (Nisbet et al., 1982) or (Murray, 2002).

Simple models based on the diffusion term (47) are used to model insects population dynamics in (Kareiva, 1983) where the author provides constant values of $\sigma$ for 12 species based on experimental data.

**Advection term.** This term models the advection motion of the population toward a determined direction, see e.g. (Murray, 2002):

$$ (48) \qquad F_a(p; x, y, t) = \vec{v}_i(x, y, t) p(x, y, t) \Leftrightarrow \mathcal{F}_a(s) = \frac{1}{q(t)} \vec{v}_i(x, y, t) s(x, y, t) $$

with $\vec{v}_i [m \cdot s^{-1}]$ the velocity field of the population (assumed here independent of $p$).

For instance, one can model the attraction of the population toward the preferred habitat by choosing the velocity field $\vec{v}_i(x, y) = \nabla \phi(x, y)$, where $\phi$ denotes the iso-distance line function from this habitat. This term can also be understood as an attractive force produced by the potential function $\phi$ representing how favorable the environment is, see e.g. (Shigesada et al., 1979) and in the review (Okubo, 1986). The two following terms are specific cases of advection terms.
To do so, one can consider a $j_{reg}$ based on a logarithmic barrier function, see e.g. (Nesterov, 2018). This approach, and barrier function in general, is in particular used in penalty method that replaces an inequality constraint (e.g. $h(s; x, y, t) < b \forall t \in [0; T] \forall (x, y) \in \Omega_c$) by a penalizing term (e.g. $\int_0^T \int_{\Omega_c} -\epsilon \log(b - h(s; x, y, t))dxdydt$ with $\Omega_c$ the constraint domain $\epsilon$ a very small constant) easier to handle in optimization problems.

In the present case, this penalizing term aims to preclude locating insects in areas that they will avoid or where they will not emit pheromones, i.e. $s = 0$ or very small on the exclusion zone $\Omega_{exc}$ derived from the biological knowledge of the insect. For instance, it is known that the FAW females emit pheromones in a vertical position on vegetation that are favorable for the growth of the larva. Hence, the FAW females will not emit pheromones in some areas that are not suitable such as roads.

Therefore, the following regularization/penalization term is considered:

$$ (45) \qquad j_{reg,exc}(s) = \int_{0}^{T} \int_{\Omega_{exc}} -\epsilon \log(\epsilon - s(x, y, t)^2) dxdydt $$

Let us note that this logarithmic barrier function is smooth (and hence differentiable) and tends to the indicator function (0 if $s = 0$ locally, $+\infty$ otherwise) when $\epsilon \to 0$.

### A.2.2. Population-dynamics models.
Recall that the quantity of pheromone emitted $s$ can be written as $s = pq$ with $q$ the quantity of pheromone emitted per insect, which is assumed to be a biological knowledge, and $p$ the density of insect. Moreover, the density $p$, and consequently $s$, follows a population dynamic model that can be written in the generic way as following:

$$ (46) \qquad \begin{aligned} & \partial_t p + \nabla \left( \sum_i F_i(p) \right) - \sum_i R_i(p) = 0 \\ & \Leftrightarrow \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 \end{aligned} $$

with $F_i$ and $R_i$ are resp. flux and reaction terms that model the behavior of the insect specie of interest, and $F_i(p = \frac{s}{q}) = \mathcal{F}_i(s)$ and $R_i(p = \frac{s}{q}) = \mathcal{R}_i(s)$.

As mentioned earlier, only a little number of terms have been introduced in this study. In this section, we aim to indicate a large list of features and resulting population dynamics terms that could be introduced in the population dynamic-informed regularization term (27).

**Diffusion term.** The most classical term modelling population dynamics is the diffusion term:

$$ (47) \qquad F_d(p; x, y, t) = -\sigma(x, y) \nabla p(x, y, t) \Leftrightarrow \mathcal{F}_d(s) = -\frac{1}{q(t)} \sigma(x, y) \nabla s(x, y, t) $$

with $\sigma [m^2 \cdot s^{-1}]$ the diffusion coefficient of the population, depending notably on the specie.

This term models the spread of the population at a macroscopic scale caused by a motion of each individual in a random direction with a probability independent of other individuals or the environment. This diffusion term can be understood as a continuous and large population limit of a discrete microscopic unbiased random walk process with a fixed probability of jump, see e.g. (Nisbet et al., 1982) or (Murray, 2002).

Simple models based on the diffusion term (47) are used to model insects population dynamics in (Kareiva, 1983) where the author provides constant values of $\sigma$ for 12 species based on experimental data.
**Quorum sensing term.** This term models the advection motion of the population but toward a direction depending on the population density:

$$ (49) \qquad F_{qs}(p; x, y, t) = \vec{v}_i(p; x, y, t)p(x, y, t) $$

with $\vec{v}_i(p)$ the density-dependent velocity field of the population.

This term enables to model the advection and spread of the population where the individuals interact with each other. The most classical density-dependent velocity field is

$$ (50) \qquad \vec{v}_i(p) = -\sigma_{qs} \nabla p $$

with $\sigma_{qs} [m^2 \cdot m^2 \cdot s^{-1}]$ the diffusion coefficient per density unit. This velocity field enables to model a motion due to avoidance of overcrowded places and population pressure, see e.g. (Murray, 2002). It can also be understood as the diffusion of the population at a macroscopic scale caused by a motion of each individual in a biased random direction with a probability increasing with the density of individuals, see e.g. (Murray, 2002; Nisbet et al., 1982). This density-dependent velocity field (50) leads to the term in variable s

$$ (51) \qquad \mathcal{F}_{qs}(s) = -\frac{\sigma_{qs}}{q^2} s \nabla s $$

Avoidance of overcrowded places and of population pressure is a common behavior among insects, see e.g. (Murray, 2002; Shigesada, 1980) and references therein. For instance, (Shigesada, 1980) proposes and applies a model with the term (49) and the velocity field (50) to an ant lion population. This model is based on experiments highlighting population pressure in this population.

**Pheromone sensing term.** This term models the advection motion of the population but toward a direction depending on the concentration of pheromone:

$$ (52) \qquad F_{ps}(p; x, y, t) = \chi p(x, y, t) \nabla c(s; x, y, t) \Leftrightarrow \mathcal{F}_{ps}(s; x, y, t) = \frac{\chi}{q(t)} s(x, y, t) \nabla c(s; x, y, t) $$

with $\chi [m^2 \cdot g^{-1} \cdot m^2 \cdot s^{-1}]$ the chemotactic gain constant ($\chi \in \mathbb{R}^{+*}$).

The insects employ pheromones to attract individuals of the same sex, such as in the case of sex pheromones for FAW, or of both sexes, as is the case with aggregation pheromones for the red palm weevils. To reach the source of the pheromone, the insects moves in the direction of the gradient of concentration of pheromones. This motion induced by a chemical signal is commonly known as chemotaxis, especially in cellular biology. A Keller-Segel-type system (Keller et al., 1970, 1971) can be used to model at the same time both the population dynamics (here of the insects) and the dynamics of the chemical attractant (here the pheromone), see e.g. (Murray, 2002; Perthame, 2006).

For instant, (Ganga Prasath et al., 2023) proposes a Keller-Segel-type system to model the dynamics of swarm of bees communicating using a pheromone.

**Logistic term.** Introduced by Verhulst in the 19th century, the logistic growth is used to model a self-limiting birth-death process, see e.g. (Murray, 2002) or (Perthame, 2006) and reference therein:

$$ (53) \qquad R_l(p; x, y, t) = \frac{r}{\kappa} p(x, y, t)(\kappa - p(x, y, t)) \Leftrightarrow R_l(s; x, y, t) = \frac{1}{q(t)^2} \frac{r}{\kappa} s(x, y, t)(\kappa q(t) - s(x, y, t)) $$

with $r [s^{-1}]$ the linear reproduction rate ($r \in \mathbb{R}^{+*}$) and $\kappa [m^{-2}]$ the carrying capacity of the environment ($\kappa \in \mathbb{R}^{+*}$).

This term enables to model a saturation phenomena that arises when the maximal occupation density (the carrying capacity) is exceeded ($p > \kappa$, the death process overcome the birth process),
**Advection term.** This term models the advection motion of the population toward a determined direction, see e.g. (Murray, 2002):

$$ (48) \qquad F_a(p; x, y, t) = \vec{v}_i(x, y, t)p(x, y, t) \Leftrightarrow \mathcal{F}_a(s) = \frac{1}{q(t)} \vec{v}_i(x, y, t)s(x, y, t) $$

with $\vec{v}_i [m.s^{-1}]$ the velocity field of the population (assumed here independent of $p$).

For instance, one can model the attraction of the population toward the preferred habitat by choosing the velocity field $\vec{v}_i(x, y) = \nabla\phi(x, y)$, where $\phi$ denotes the iso-distance line function from this habitat. This term can also be understood as an attractive force produced by the potential function $\phi$ representing how favorable the environment is, see e.g. (Shigesada et al., 1979) and in the review (Okubo, 1986). The two following terms are specific cases of advection terms.

**Quorum sensing term.** This term models the advection motion of the population but toward a direction depending on the population density:

$$ (49) \qquad F_{qs}(p; x, y, t) = \vec{v}_i(p; x, y, t)p(x, y, t) $$

with $\vec{v}_i(p)$ the density-dependent velocity field of the population.

This term enables to model the advection and spread of the population where the individuals interact with each other. The most classical density-dependent velocity field is

$$ (50) \qquad \vec{v}_i(p) = -\sigma_{qs} \nabla p $$

with $\sigma_{qs} [m^2.m^2.s^{-1}]$ the diffusion coefficient per density unit. This velocity field enables to model a motion due to avoidance of overcrowded places and population pressure, see e.g. (Murray, 2002). It can also be understood as the diffusion of the population at a macroscopic scale caused by a motion of each individual in a biased random direction with a probability increasing with the density of individuals, see e.g. (Murray, 2002; Nisbet et al., 1982). This density-dependent velocity field (50) leads to the term in variable $s$

$$ (51) \qquad \mathcal{F}_{qs}(s) = -\frac{\sigma_{qs}}{q^2} s \nabla s $$

Avoidance of overcrowded places and of population pressure is a common behavior among insects, see e.g. (Murray, 2002; Shigesada, 1980) and references therein. For instance, (Shigesada, 1980) proposes and applies a model with the term (49) and the velocity field (50) to an ant lion population. This model is based on experiments highlighting population pressure in this population.

**Pheromone sensing term.** This term models the advection motion of the population but toward a direction depending on the concentration of pheromone:

$$ (52) \qquad F_{ps}(p; x, y, t) = \chi p(x, y, t) \nabla c(s; x, y, t) \Leftrightarrow \mathcal{F}_{ps}(s; x, y, t) = \frac{\chi}{q(t)} s(x, y, t) \nabla c(s; x, y, t) $$

with $\chi [m^2.g^{-1}.m^2.s^{-1}]$ the chemotactic gain constant ($\chi \in \mathbb{R}^{+*}$).

The insects employ pheromones to attract individuals of the same sex, such as in the case of sex pheromones for FAW, or of both sexes, as is the case with aggregation pheromones for the red palm weevils. To reach the source of the pheromone, the insects moves in the direction of the gradient of concentration of pheromones. This motion induced by a chemical signal is commonly known as chemotaxis, especially in cellular biology. A Keller-Segel-type system (Keller et al., 1970, 1971) can be used to model at the same time both the population dynamics (here of the insects) and the dynamics of the chemical attractant (here the pheromone), see e.g. (Murray, 2002; Perthame, 2006).
and invasion front. This saturation is for instance due to intra-specific competition for a limited resource that is consumed faster than it can renewed.

The equation obtained when considering diffusion and logistic source is the well-known Fisher-KPP equation, introduced in (RA Fisher, 1937) to model the spread of favorable gene and first used for population dynamics problems in (Skellam, 1951). A Fisher-KPP-type system is for instance used in (Takahashi et al., 2005) to model the dynamics of mosquitoes populations where the clear water (required for egg-laying) and the amount of nutrient (human blood) are limited, potentially leading to saturation phenomenon.

**Allee effect term.** The Allee effect is characterized by the decrease of the population due to a small or sparse population. The following term is an extension of the logistic term (53) which additionally models the Allee effect, see e.g. a review in (Sun, 2016) and reference therein:

$$ (54) \qquad R_{\text{Ae}}(p; x, y, t) = \frac{r}{\kappa_+\kappa_-} p(x, y, t)(\kappa_+ - p(x, y, t))(p(x, y, t) - \kappa_-) \\ \Leftrightarrow R_{\text{Ae}}(s; x, y, t) = \frac{1}{q(t)^3} \frac{r}{\kappa_+\kappa_-} s(x, y, t)(q(t)\kappa_+ - s(x, y, t))(s(x, y, t) - q(t)\kappa_-) $$

with $r [s^{-1}]$ the linear reproduction rate ($r \in \mathbb{R}^{+*}$), $\kappa_+ [m^{-2}]$ the carrying capacity of the environment ($\kappa_+ \in \mathbb{R}^{+*}$) and $\kappa_- [m^{-2}]$ the Allee threshold ($\kappa_- \in \mathbb{R}^{+*}$ and $\kappa_- < \kappa_+$) that is the critical density under which the population decrease.

This effect is e.g. due to the lack of interactions between individuals of the population. In insect populations, the Allee effect commonly occurs since lower densities reduce the chance of an individual to locate a mate for purpose of reproduction.

Since the Allee effect can lead to the extinction of the insect population, it is often used in pest control strategy, see e.g. (Boukal et al., 2009).

### B.2.3. Non-local terms.
The non-local terms are generalization of the previously presented terms. These non-local terms enable to model non-local interaction and long range diffusion, leading for instance to spatial patterns, see e.g. (Furter et al., 1989), (Genieys et al., 2006) or (Murray, 2002) and reference therein. Non-local terms arise especially when the individuals of the population communicate with each other remotely, e.g. using chemical or visual signals, or when the individuals interact with the environment not local but surrounding.

For example, one can extend the logistic term (53) into a non-local version:

$$ (55) \qquad R_{nll}(p; x, y, t) = \frac{r}{\kappa} p(x, y, t)(\kappa - k * p(x, y, t)) \\ \Leftrightarrow R_{nll}(s; x, y, t) = \frac{1}{q(t)^2} \frac{r}{\kappa} s(x, y, t)(\kappa q(t) - k * s(x, y, t)) $$

with the convolution product $k * s(x, y, t) = \int_{\Omega} k(x - x', y - y')s(x', y', t)dx'dy'$ and k the convolution kernel.

This implies that individuals are not competing for the resource locally available, as with the term (53), but for the resource available in a vicinity. In this term, the convolution product can be seen as an average weighted by the kernel k that defines how the individuals interact with the vicinity and how large is this vicinity. Let us note that the local term (53) is retrieved if k is the Dirac kernel $\delta$. Similarly, one can consider a non-local Allee effect term to model the decrease of the population due to a small or sparse population in a vicinity. For instance, (Britton, 1989) proposed a non-local version of the Fisher-KPP equation based on this non-local logistic term.
For instant, (Ganga Prasath et al., 2023) proposes a Keller-Segel-type system to model the dynamics of swarm of bees communicating using a pheromone.

**Logistic term.** Introduced by Verhulst in the 19th century, the logistic growth is used to model a self-limiting birth-death process, see e.g. (Murray, 2002) or (Perthame, 2006) and reference therein:

(53)

$$R_I(p; x, y, t) = \frac{r}{\kappa} p(x, y, t)(\kappa - p(x, y, t)) \Leftrightarrow R_I(s; x, y, t) = \frac{1}{q(t)^2} \frac{r}{\kappa} s(x, y, t)(\kappa q(t) - s(x, y, t))$$

with $r$ [$s^{-1}$] the linear reproduction rate ($r \in \mathbb{R}^{+*}$) and $\kappa$ [$m^{-2}$] the carrying capacity of the environment ($\kappa \in \mathbb{R}^{+*}$).

This term enables to model a saturation phenomena that arises when the maximal occupation density (the carrying capacity) is exceeded ($p > \kappa$, the death process overcome the birth process), and invasion front. This saturation is for instance due to intra-specific competition for a limited resource that is consumed faster than it can renewed.

The equation obtained when considering diffusion and logistic source is the well-known Fisher-KPP equation, introduced in (RA Fisher, 1937) to model the spread of favorable gene and first used for population dynamics problems in (Skellam, 1951). A Fisher-KPP-type system is for instance used in (Takahashi et al., 2005) to model the dynamics of mosquitoes populations where the clear water (required for egg-laying) and the amount of nutrient (human blood) are limited, potentially leading to saturation phenomenon.

**Allee effect term.** The Allee effect is characterized by the decrease of the population due to a small or sparse population. The following term is an extension of the logistic term (53) which additionally models the Allee effect, see e.g. a review in (Sun, 2016) and reference therein:

$$ (54) \qquad \begin{aligned} R_{Ae}(p; x, y, t) &= \frac{r}{\kappa_+ \kappa_-} p(x, y, t)(\kappa_+ - p(x, y, t))(p(x, y, t) - \kappa_-) \\ \Leftrightarrow R_{Ae}(s; x, y, t) &= \frac{1}{q(t)^3} \frac{r}{\kappa_+ \kappa_-} s(x, y, t)(q(t)\kappa_+ - s(x, y, t))(s(x, y, t) - q(t)\kappa_-) \end{aligned} $$

with $r$ [$s^{-1}$] the linear reproduction rate ($r \in \mathbb{R}^{+*}$), $\kappa_+$ [$m^{-2}$] the carrying capacity of the environment ($\kappa_+ \in \mathbb{R}^{+*}$) and $\kappa_-$ [$m^{-2}$] the Allee threshold ($\kappa_- \in \mathbb{R}^{+*}$ and $\kappa_- < \kappa_+$) that is the critical density under which the population decrease.

This effect is e.g. due to the lack of interactions between individuals of the population. In insect populations, the Allee effect commonly occurs since lower densities reduce the chance of an individual to to locate a mate for purpose of reproduction.

Since the Allee effect can lead to the extinction of the insect population, it is often used in pest control strategy, see e.g. (Boukal et al., 2009).

### A.2.3. Non-local terms.
The non-local terms are generalization of the previously presented terms. These non-local terms enable to model non-local interaction and long range diffusion, leading for instance to spatial patterns, see e.g. (Furter et al., 1989), (Genieys et al., 2006) or (Murray, 2002) and reference therein. Non-local terms arise especially when the individuals of the population communicate with each other remotely, e.g. using chemical or visual signals, or when the individuals interact with the environment not local but surrounding.

For example, one can extend the logistic term (53) into a non-local version:

$$ (55) \qquad \begin{aligned} R_{nll}(p; x, y, t) &= \frac{r}{\kappa} p(x, y, t)(\kappa - k * p(x, y, t)) \\ \Leftrightarrow R_{nll}(s; x, y, t) &= \frac{1}{q(t)^2} \frac{r}{\kappa} s(x, y, t)(\kappa q(t) - k * s(x, y, t)) \end{aligned} $$
As mentioned before, a Keller-Segel type system can be used to model both the dynamics of the population and the dynamics of the chemical attractant, in which the quantity of chemical attractant emitted depends on the population density. In special cases, the equation of the chemical attractant dynamics can be analytically solved. The concentration of the attractant can then be expressed as a convolutional product between the density of the population and a kernel describing the attractant dynamics, which leads to a non-local flux term in the population dynamics model, see e.g. (Perthame, 2006).

For example, in the present case, one can approximate the pheromone dynamics, described in the direct model (1), by a Gaussian plume $c_{gp}(s) = g_{gp} * s$ with $g_{gp}$ the Gaussian plume kernel. Injecting this approximation in the pheromone sensing term (52) leads to a non-local pheromone sensing term:

$$ (56) \qquad \mathcal{F}_{nlps}(s; x, y, t) = \frac{\chi}{q(t)} s(x, y, t) \nabla(g_{gp} * s(x, y, t)) $$

Non-local terms of the shape $R_{lr}(p) \propto k * p$ can be used to model long range diffusion processes. The range of the diffusion processes depends on the choice of the kernel $k$ and its moments, see e.g. (Murray, 2002).

### B.3. On solving the optimization problem (3)

In this section, we present all the calculation details for solving the optimization problem (3) presented in Sec. 2.2. In the present context, the cost function $j$ is optimized using gradient descent or proximal gradient methods. Recall that the gradient of $j$, or at least of its differentiable part, is given by

$$ (57) \qquad \nabla_{s}j = \nabla_{s}j_{obs} + \nabla_{s}j_{reg} $$
with $j_{reg} = \sum_{i} \alpha_{reg,i}j_{reg,i} \Rightarrow \nabla_{s}j_{reg} = \sum_{i} \alpha_{reg,i} \nabla_{s}j_{reg,i}$.

In Sec. 2.2.2, it is mentioned that $\nabla_{s}j_{obs} = -c^*$ with $c^*$ the adjoint state, solution of the adjoint model (6). In a first part, the derivation of the adjoint model and of the expression (7) are detailed. Moreover, related models used in this study, such that the one-sensor adjoint model, will be detailed.

In Secs. 3.1 and B.2, several Biology-Informed differentiable $j_{reg,i}$ are proposed. In a second part, the gradients of these $j_{reg,i}$ are calculated.

In Sec. 3.1.2, LASSO and group-LASSO regularization terms are proposed. However, these $j_{reg}$ are not differentiable. Thus, a proximal gradient algorithm is used in this study. In a third part, we present the proximal-gradient algorithm and its application to the LASSO and group-LASSO regularization terms, called ISTA algorithm.

#### B.3.1. Derivation of the adjoint model and related.
In this section, $K_T$ denotes the control space and $K_T = \mathcal{L}^2([0; T]; \mathcal{K})$ with $\mathcal{K}$ the space of functions from $\Omega \subset \mathbb{R}^2$ to $\mathbb{R}$ which regularity in space is required in order to solve the direct model (1). We assume that the spaces $K_T$ and $\mathcal{K}$ are Hilbert spaces and we have

$$ < f, g >_{K_T} = \int_{0}^{T} < f(t, \cdot), g(t, \cdot) >_{K} dt. $$

In a similar way, $V_T = \{f; f \in \mathcal{L}^2([0; T]; V) \text{ and } \partial_t f \in \mathcal{L}^2([0; T]; V)\}$ denotes the space of the state variable, solution of the direct model and we have $< f, g >_{V_T} = \int_{0}^{T} < f(t, \cdot), g(t, \cdot) >_{V} dt$. Moreover, $\mathcal{D}$ denotes the observations space.
with the convolution product $k * s(x, y, t) = \int_{\Omega} k(x - x', y - y')s(x', y', t)dx'dy'$ and k the convolution kernel.

This implies that individuals are not competing for the resource locally available, as with the term (53), but for the resource available in a vicinity. In this term, the convolution product can be seen as an average weighted by the kernel k that defines how the individuals interact with the vicinity and how large is this vicinity. Let us note that the local term (53) is retrieved if k is the Dirac kernel $\delta$. Similarly, one can consider a non-local Allee effect term to model the decrease of the population due to a small or sparse population in a vicinity. For instance, (Britton, 1989) proposed a non-local version of the Fisher-KPP equation based on this non-local logistic term.

As mentioned before, a Keller-Segel type system can be used to model both the dynamics of the population and the dynamics of the chemical attractant, in which the quantity of chemical attractant emitted depends on the population density. In special cases, the equation of the chemical attractant dynamics can be analytically solved. The concentration of the attractant can then be expressed as a convolutional product between the density of the population and a kernel describing the attractant dynamics, which leads to a non-local flux term in the population dynamics model, see e.g. (Perthame, 2006).

For example, in the present case, one can approximate the pheromone dynamics, described in the direct model (1), by a Gaussian plume $c_{gp}(s) = g_{gp} * s$ with $g_{gp}$ the Gaussian plume kernel. Injecting this approximation in the pheromone sensing term (52) leads to a non-local pheromone sensing term:

$$ (56) \qquad \mathcal{F}_{nlps}(s; x, y, t) = \frac{\chi}{q(t)} s(x, y, t) \nabla(g_{gp} * s(x, y, t)) $$

Non-local terms of the shape $R_{lr}(p) \propto k * p$ can be used to model long range diffusion processes. The range of the diffusion processes depends on the choice of the kernel k and its moments, see e.g. (Murray, 2002).

### A.3. On solving the optimization problem (3)

In this section, we present all the calculation details for solving the optimization problem (3) presented in Sec. 2.2. In the present context, the cost function j is optimized using gradient descent or proximal gradient methods. Recall that the gradient of j, or at least of its differentiable part, is given by

$$ (57) \qquad \nabla_s j = \nabla_s j_{obs} + \nabla_s j_{reg} $$

with $j_{reg} = \sum_i \alpha_{reg,i} j_{reg,i} \Rightarrow \nabla_s j_{reg} = \sum_i \alpha_{reg,i} \nabla_s j_{reg,i}$.

In Sec. 2.2.2, it is mentioned that $\nabla_s j_{obs} = -c^*$ with $c^*$ the adjoint state, solution of the adjoint model (6). In a first part, the derivation of the adjoint model and of the expression (7) are detailed. Moreover, related models used in this study, such that the one-sensor adjoint model, will be detailed.

In Secs. 3.1 and A.2, several Biology-Informed differentiable $j_{reg,i}$ are proposed. In a second part, the gradients of these $j_{reg,i}$ are calculated.

In Sec. 3.1.2, LASSO and group-LASSO regularization terms are proposed. However, these $j_{reg}$ are not differentiable. Thus, a proximal gradient algorithm is used in this study. In a third part, we present the proximal-gradient algorithm and its application to the LASSO and group-LASSO regularization terms, called ISTA algorithm.
As mentioned previously, the gradient of the cost function is $\nabla_{s}j(s) = \nabla_{s}j_{obs}(s) + \alpha_{reg} \nabla_{s}j_{reg}(s)$ with the gradient of the observation term (4) that is: $\forall \delta s \in \mathcal{K}_{T}$

$$ (58) \quad \begin{aligned} \langle \nabla_{s}j_{obs}(s), \delta s \rangle_{\mathcal{K}_{T}} &= 2 \langle \mathbf{R}^{-1}(m(c(s)) - m^{\text{obs}}), \frac{dm}{dc}(c(s)) \cdot \frac{dc}{ds}(s) \cdot \delta s \rangle_{\mathcal{D}} \\ &= 2 \left( \frac{dm}{dc}(c(s)) \right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \frac{dc}{ds}(s) \cdot \delta s \rangle_{\mathcal{V}_{T}} \\ &= 2 \left( \frac{dc}{ds}(s) \right)^{*} \cdot \left( \frac{dm}{dc}(c(s)) \right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \delta s \rangle_{\mathcal{K}_{T}} \end{aligned} $$

with $s \mapsto c(s)$ the operator that given $s$ returns the solution of the direct model (1) and $d_s c(s)$ the derivative of this operator with respect to $s$. However, the operator $(d_s c(s))^*$ may not be easy to calculate. The adjoint model provides a method to estimate the results of this operator.

**Linear tangent model.** In this section, the direct model is written as:

$$ (59) \quad \frac{\partial c}{\partial t} + \mathcal{A}(c) = \mathcal{B}(s) $$

In our context, the differentiable operators are:

$$ (60) \quad \mathcal{A}(c) = \nabla \cdot (\vec{u}c - \mathbf{K}\nabla c) + \tau_{\text{loss}}c \quad \text{and} \quad \mathcal{B}(s) = s $$

with the boundary conditions:

* 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.

Denoting $\omega = d_s c \cdot \delta s \in \mathcal{V}_T$, the linear tangent model is then:

$$ (61) \quad \frac{\partial \omega}{\partial t} + \frac{d\mathcal{A}}{dc}(c(s)) \cdot \omega = \frac{d\mathcal{B}}{ds}(s) \cdot \delta s $$

with the initial condition $\omega(t=0) = 0$ (since $c(t=0) = 0 \quad \forall x, y \forall s$).

Let us note that, in the present case, $d_s \mathcal{A} = \mathcal{A}$ since $\mathcal{A}$ is a linear operator. By applying the scalar product with a test function $z \in \mathcal{V}_T$ and integration by part, we then have:

$$ (62) \quad \left\langle \frac{\partial z}{\partial t} - \left( \frac{d\mathcal{A}}{dc}(c(s)) \right)^* \cdot z, \omega \right\rangle_{\mathcal{V}_T} - \left[ \left\langle z, \omega \right\rangle_{\mathcal{V}} \right]_{t=0}^{t=T} = - \left\langle \frac{d\mathcal{B}}{ds}(s) \cdot \delta s, z \right\rangle_{\mathcal{V}_T} \quad \forall z \in \mathcal{V}_T $$

**Adjoint model.** On the other hand, based on the expression of the gradient of $j$ (57) and of $j_{obs}$ (58), the gradient of $j$ satisfies: $\forall \delta s \in \mathcal{K}_{T}$

$$ (63) \quad \langle \nabla_{s}j, \delta s \rangle_{\mathcal{K}_{T}} = \langle \nabla_{s}j_{reg}, \delta s \rangle_{\mathcal{K}_{T}} + 2 \left\langle \frac{dm}{dc}(c(s)) \right\rangle^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \omega \rangle_{\mathcal{V}_{T}} $$

By introducing the linear tangent model (62) in the equation of gradient (63), we then have: $\forall \delta s \in \mathcal{K}_{T}$ and $\forall z \in \mathcal{V}_{T}$

(64)

$$ \langle \nabla_{s}j, \delta s \rangle_{\mathcal{K}_{T}} = \langle \nabla_{s}j_{reg}, \delta s \rangle_{\mathcal{K}_{T}} - \left\langle \frac{d\mathcal{B}}{ds}(s) \cdot \delta s, z \right\rangle_{\mathcal{V}_{T}} $$

$$ (65) \quad \begin{aligned} & - \left\langle \frac{\partial z}{\partial t} - \left( \frac{d\mathcal{A}}{dc}(c(s)) \right)^* \cdot z - 2 \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \omega \right\rangle_{\mathcal{V}_{T}} \\ & + \left[ \left\langle z, \omega \right\rangle_{\mathcal{V}} \right]_{t=0}^{t=T} \end{aligned} $$
### A.3.1. Derivation of the adjoint model and related.
In this section, $\mathcal{K}_T$ denotes the control space and $\mathcal{K}_T = \mathcal{L}^2([0; T]; \mathcal{K})$ with $\mathcal{K}$ the space of functions from $\Omega \subset \mathbb{R}^2$ to $\mathbb{R}$ which regularity in space is required in order to solve the direct model (1). We assume that the spaces $\mathcal{K}_T$ and $\mathcal{K}$ are Hilbert spaces and we have

$$< f, g >_{\mathcal{K}_T} = \int_0^T < f(t, \cdot), g(t, \cdot) >_{\mathcal{K}} dt.$$

In a similar way, $\mathcal{V}_T = \{f; f \in \mathcal{L}^2([0; T]; \mathcal{V})\}$ and $\partial_t f \in \mathcal{L}^2([0; T]; \mathcal{V})$} denotes the space of the state variable, solution of the direct model and we have $< f, g >_{\mathcal{V}_T} = \int_0^T < f(t, \cdot), g(t, \cdot) >_{\mathcal{V}} dt$. Moreover, $\mathcal{D}$ denotes the observations space.

As mentioned previously, the gradient of the cost function is $\nabla_{s}j(s) = \nabla_{s}j_{obs}(s) + \alpha_{reg} \nabla_{s}j_{reg}(s)$ with the gradient of the observation term (4) that is: $\forall \delta s \in \mathcal{K}_T$

$$ (58) \quad \begin{aligned} < \nabla_{s}j_{obs}(s), \delta s >_{\mathcal{K}_T} &= 2 < \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \frac{dm}{dc}(c(s)) \cdot \frac{dc}{ds}(s) \cdot \delta s >_{\mathcal{D}} \\ &= 2 < \left(\frac{dm}{dc}(c(s))\right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \frac{dc}{ds}(s) \cdot \delta s >_{\mathcal{V}_T} \\ &= 2 < \left(\frac{dc}{ds}(s)\right)^{*} \cdot \left(\frac{dm}{dc}(c(s))\right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \delta s >_{\mathcal{K}_T} \end{aligned} $$

with $s \mapsto c(s)$ the operator that given $s$ returns the solution of the direct model (1) and $d_s c(s)$ the derivative of this operator with respect to $s$. However, the operator $(d_s c(s))^*$ may not be easy to calculate. The adjoint model provides a method to estimate the results of this operator.

**Linear tangent model.** In this section, the direct model is written as :

$$ (59) \qquad \frac{\partial c}{\partial t} + \mathcal{A}(c) = \mathcal{B}(s) $$

In our context, the differentiable operators are:

$$ (60) \qquad \mathcal{A}(c) = \nabla \cdot (\vec{u}c - \mathbf{K}\nabla c) + \tau_{\text{loss}}c \quad \text{and} \quad \mathcal{B}(s) = s $$

with the boundary conditions:

* 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.

Denoting $\omega = d_s c \cdot \delta s \in \mathcal{V}_T$, the linear tangent model is then:

$$ (61) \qquad \frac{\partial \omega}{\partial t} + \frac{d\mathcal{A}}{dc}(c(s)) \cdot \omega = \frac{d\mathcal{B}}{ds}(s) \cdot \delta s $$

with the initial condition $\omega(t=0) = 0$ (since $c(t=0) = 0 \quad \forall x, y \forall s$).

Let us note that, in the present case, $d_s \mathcal{A} = \mathcal{A}$ since $\mathcal{A}$ is a linear operator. By applying the scalar product with a test function $z \in \mathcal{V}_T$ and integration by part, we then have:

$$ (62) \quad < \frac{\partial z}{\partial t} - \left( \frac{d\mathcal{A}}{dc}(c(s)) \right)^* \cdot z, \omega >_{\mathcal{V}_T} - \left[ < z, \omega >_{\mathcal{V}} \right]_{t=0}^{t=T} = - < \frac{d\mathcal{B}}{ds}(s) \cdot \delta s, z >_{\mathcal{V}_T} \quad \forall z \in \mathcal{V}_T $$

**Adjoint model.** On the other hand, based on the expression of the gradient of $j$ (57) and of $j_{obs}$ (58), the gradient of $j$ satisfies: $\forall \delta s \in \mathcal{K}_T$

$$ (63) \quad < \nabla_{s}j, \delta s >_{\mathcal{K}_T} = < \nabla_{s}j_{reg}, \delta s >_{\mathcal{K}_T} + 2 < \left( \frac{dm}{dc}(c(s)) \right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \omega >_{\mathcal{V}_T} $$
Thus, the adjoint model aims to find a particular test function $c^*$ of $V_T$, so-called adjoint state, that satisfies:

$$ (66) \qquad \frac{\partial c^*}{\partial t} - \left( \frac{dA}{dc}(c(s)) \right)^* \cdot c^* = 2 \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}) $$

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

For the differential operator $A$ given by Eq. (60), the adjoint operator $(d_c A)^*$ is simply $(d_c A)^* \cdot z = -\nabla \cdot (\mathbf{K}^T \nabla z) - \nabla(\vec{u}z) + (\nabla \cdot \vec{u})z + \tau_{\text{loss}}z$ and the boundary conditions of the adjoint model are:

* 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\} \forall t \in ]0; T[$ with $\vec{n}$ the outgoing normal vector.

Let us note that by definition, the adjoint model is a linear model.

**Computation of the gradient.** By introducing the adjoint state as test function in the equation (64), it returns that $\nabla_{sj}$ satisfies the following equation:

$$ (67) \qquad \langle \nabla_{sj}(s), \delta s \rangle_{\mathcal{K}_T} = \langle \nabla_{sj_{\text{reg}}}(s), \delta s \rangle_{\mathcal{K}_T} - \langle \left( \frac{dB}{ds}(s) \right)^*, c^*(c(s)), \delta s \rangle_{\mathcal{V}_T} \quad \forall \delta s \in \mathcal{K}_T $$

with $(d_s B)^* = \text{Id}$ for the differential operator $B$ given by (60). Thus, it results the expression of the gradient:

$$ (68) \qquad \nabla_{sj}(s) = \nabla_{sj_{\text{reg}}}(s) - \left( \frac{dB}{ds}(s) \right)^* \cdot c^*(c(s)) $$

To estimate the gradient at each optimization iteration, it requires to solve in order the direct model (1) and then the adjoint model (6).

**One-sensor and trajectory-integrated one-sensor adjoint models.** Recall that the one-sensor adjoint model is mostly obtained by splitting the observation operator $c \mapsto m(c)$ into the sum of one-sensor observation operator $m = \sum_i m_i$. By introducing this observation operator splitting in the adjoint model (6), the adjoint model can be written the following way:

$$ \partial_t c^* + \nabla \cdot (\mathbf{K}^T \nabla c^*) + \nabla(\vec{u}c^*) - (\nabla \cdot \vec{u})c^* - \tau_{\text{loss}}c^* = \sum_i \sum_k \left( \frac{dm_k}{dc}(c(s)) \right)^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{\text{obs}}) $$

Assuming that $(\frac{dm_k}{dc}(c(s)))^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{\text{obs}}) = 0$ for $i \neq k$ and as the left hand side is linear with respect to $c^*$, the adjoint model written (B.3.1) can be split into the one-sensor adjoint models (8). Let us note that, as mentioned in Sec. 2.2.3, the hypothesis $(\frac{dm_k}{dc}(c(s)))^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{\text{obs}}) = 0$ for $i \neq k$ is satisfied in our case if the sensors are independent $(\mathbf{R}^{-1}$ is block-diagonal) and if two sensors are not located in the same place.

Recall from Sec. 2.2.3 that the resulting one-sensor adjoint state is the contribution of a sensor to the gradient $\nabla_{sj_{\text{obs}}}$.

In this study, the sequence $(s_j)_{j \in \mathbb{N}}$ given by the trajectory taken by the optimization algorithm using different DA strategies are compared through one-sensor adjoint state-based criteria. The main one-sensor adjoint state-based criteria used, introduced in Sec. 2.4.3, is the maximal total contribution $\| \sum_j \eta_j c_i^*(s_j) \|_{L^\infty(\Omega \times ]0, T])}$, with $(\eta_j)_{j \in \mathbb{N}}$ the sequence of the step sizes and $c_i^*(s_j)$ is the solution of the one-sensor adjoint model:

$$ \partial_t c_{i,j}^* + \nabla \cdot (\mathbf{K}^T \nabla c_{i,j}^*) + \nabla(\vec{u}c_{i,j}^*) - (\nabla \cdot \vec{u})c_{i,j}^* - \tau_{\text{loss}}c_{i,j}^* = (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{\text{obs}}) $$

Thus, by multiplying with the step size and summing over the sequences the previous one-sensor adjoint model and due the linearity of the left hand side (that is always true by definition of the
By introducing the linear tangent model (62) in the equation of gradient (63), we then have:

$\forall \delta s \in \mathcal{K}_T$ and $\forall z \in \mathcal{V}_T$

$$ (64) $$

$$ < \nabla_s j, \delta s >_{\mathcal{K}_T} = < \nabla_{sjreg}, \delta s >_{\mathcal{K}_T} - < \frac{d\mathcal{B}}{ds}(s) \cdot \delta s, z >_{\mathcal{V}_T} $$

$$ (65) \qquad - < \frac{\partial z}{\partial t} - \left( \frac{d\mathcal{A}}{dc}(c(s)) \right)^* \cdot z - 2 \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \omega >_{\mathcal{V}_T} \\ \qquad \qquad + [ < z, \omega >_{\mathcal{V}} ]_{t=0}^{t=T} $$

Thus, the adjoint model aims to find a particular test function $c^*$ of $\mathcal{V}_T$, so-called adjoint state, that satisfies:

$$ (66) \qquad \frac{\partial c^*}{\partial t} - \left( \frac{d\mathcal{A}}{dc}(c(s)) \right)^* \cdot c^* = 2 \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}) $$

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

For the differential operator $\mathcal{A}$ given by Eq. (60), the adjoint operator $(d_c\mathcal{A})^*$ is simply $(d_c\mathcal{A})^*$. $z = -\nabla \cdot (\mathbf{K}^T \nabla z) - \nabla(\vec{u}z) + (\nabla \cdot \vec{u})z + \tau_{\text{loss}}z$ and the boundary conditions of the adjoint model are:

* 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.

Let us note that by definition, the adjoint model is a linear model.

**Computation of the gradient.** By introducing the adjoint state as test function in the equation (64), it returns that $\nabla_s j$ satisfies the following equation:

$$ (67) \qquad < \nabla_s j(s), \delta s >_{\mathcal{K}_T} = < \nabla_{sjreg}(s), \delta s >_{\mathcal{K}_T} - < \left( \frac{d\mathcal{B}}{ds}(s) \right)^* \cdot c^*(c(s)), \delta s >_{\mathcal{V}_T} \forall \delta s \in \mathcal{K}_T $$

with $(d_s\mathcal{B})^* = \text{Id}$ for the differential operator $\mathcal{B}$ given by (60). Thus, it results the expression of the gradient:

$$ (68) \qquad \nabla_s j(s) = \nabla_{sjreg}(s) - \left( \frac{d\mathcal{B}}{ds}(s) \right)^* \cdot c^*(c(s)) $$

To estimate the gradient at each optimization iteration, it requires to solve in order the direct model (1) and then the adjoint model (6).

**One-sensor and trajectory-integrated one-sensor adjoint models.** Recall that the one-sensor adjoint model is mostly obtained by splitting the observation operator $c \mapsto m(c)$ into the sum of one-sensor observation operator $m = \sum_i m_i$. By introducing this observation operator splitting in the adjoint model (6), the adjoint model can be written the following way:

$$ \partial_t c^* + \nabla \cdot (\mathbf{K}^T \nabla c^*) + \nabla(\vec{u}c^*) - (\nabla \cdot \vec{u})c^* - \tau_{\text{loss}}c^* = \sum_i \sum_k \left( \frac{dm_k}{dc}(c(s)) \right)^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{\text{obs}}) $$

Assuming that $(\frac{dm_k}{dc}(c(s)))^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{\text{obs}}) = 0$ for $i \neq k$ and as the left hand side is linear with respect to $c^*$, the adjoint model written (A.3.1) can be split into the one-sensor adjoint models (8). Let us note that, as mentioned in Sec. 2.2.3, the hypothesis $(\frac{dm_k}{dc}(c(s)))^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{\text{obs}}) = 0$ for $i \neq k$ is satisfied in our case if the sensors are independent $(\mathbf{R}^{-1}$ is block-diagonal) and if two sensors are not located in the same place.

Recall from Sec. 2.2.3 that the resulting one-sensor adjoint state is the contribution of a sensor to the gradient $\nabla_{sjobs}$.
adjoint model), $c_{i,tot}^* = \sum_j \eta_j c_i^*(s_j)$ is solution of the following trajectory-integrated one-sensor adjoint model:

$$ (69) \qquad \begin{aligned} & \partial_t c_{i,tot}^* + \nabla \cdot (\mathbf{K}^T \nabla c_{i,tot}^*) + \nabla(\vec{u} c_{i,tot}^*) - (\nabla \cdot \vec{u}) c_{i,tot}^* - \tau_{loss} c_{i,tot}^* \\ & = \sum_j (\eta_j (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs})) \end{aligned} $$

This trajectory-integrated one-sensor adjoint model enable to compute the criteria without having to compute $c_i^*(s_j)$ for each $s_j$.

Let us note that, if the one-sensor observation operator $c \mapsto m$ is linear (which is the case in the present study), the right hand side of the trajectory-integrated one-sensor adjoint model

$$ \text{RHS} = \sum_j (\eta_j (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{\text{obs}})) $$

$$ \text{RHS} = (d_c m_i)^* \cdot 2\mathbf{R}^{-1} \cdot (\sum_j \eta_j (m_i(c(s_j)) - m_i^{\text{obs}})). $$

Moreover, if the direct model is linear with respect to $s$ (which is the case in the present study), the right hand side of the trajectory-integrated one-sensor adjoint model even becomes

$$ \text{RHS} = (d_c m_i)^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(\sum_j \eta_j s_j)) - (\sum_j \eta_j) m_i^{\text{obs}}). $$

### B.3.2. Computation of the gradients $\nabla_{s_j} j_{reg,i}$.
In this section, we calculate the gradient of the differentiable regularization terms $j_{reg,i}$ introduced in Sec. 3.1 and in Sec. B.2. Recall that $\nabla_s$ is the gradient operator with respect to $s$ and $\nabla_{s_j} j_{reg,i}(s)$ is the function given by $\frac{d j_{reg,i}(s)}{ds} \cdot \delta s = \langle \nabla_{s_j} j_{reg,i}(s), \delta s \rangle \forall \delta s$.

**Gradient of the Tikhonov regularization term.** The gradient of the Tikhonov regularization term (23), presented in Sec. 3.1.1 is simply:

$$ \nabla_{s_j} j_{reg, \tau}(s) = 2\mathbf{C}^{-1}(s - s_b) $$

**Gradient of the logarithmic barrier regularization term.** The gradient of the logarithmic barrier regularization term (45), presented in Sec. B.2.1 is simply:

$$ \nabla_{s_j} j_{reg, \text{exc}}(s) = \frac{2\epsilon}{\epsilon - s^2} \mathbb{1}_{\Omega_{\text{exc}}} $$

**Gradient of the population dynamics-informed regularization term.** Recall that the population dynamics-informed regularization term (27), presented in Sec. 3.1.3, can be written as:

$$ \begin{aligned} j_{reg,PD}(s) &= \|\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}^2 \\ &= \langle \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), \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) \rangle \end{aligned} $$

and hence, its differentiable is

$$ \frac{d j_{reg,PD}(s)}{ds} \cdot \delta s = 2 \langle \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) , \partial_t \delta s - \frac{\partial_t q}{q} \delta s + q \nabla \left( \sum_i \frac{d \mathcal{F}_i(s)}{ds} \cdot \delta s \right) + q \sum_i \frac{d \mathcal{R}_i(s)}{ds} \cdot \delta s \rangle $$

Hence $\frac{d j_{reg,PD}(s)}{ds} \cdot \delta s = \langle \nabla_{s_j} j_{reg,PD}(s), \delta s \rangle$ with:

$$ \nabla_{s_j} j_{reg,PD}(s) = 2 \left( \partial_t - \frac{\partial_t q}{q} + q \sum_i \frac{d \nabla \mathcal{F}_i(s)}{ds} + q \sum_i \frac{d \mathcal{R}_i(s)}{ds} \right)^* \left( \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) \right) $$
In this study, the sequence $(s_j)_{j \in \mathbb{N}}$ given by the trajectory taken by the optimization algorithm using different DA strategies are compared through one-sensor adjoint state-based criteria. The main one-sensor adjoint state-based criteria used, introduced in Sec. 2.4.3, is the maximal total contribution $\| \sum_j \eta_j c_i^*(s_j) \|_{L^\infty(\Omega \times [0, T])}$, with $(\eta_j)_{j \in \mathbb{N}}$ the sequence of the step sizes and $c_i^*(s_j)$ is the solution of the one-sensor adjoint model:

$$ \partial_t c_{i,j}^* + \nabla \cdot (\mathbf{K}^T \nabla c_{i,j}^*) + \nabla(\vec{u} c_{i,j}^*) - (\nabla \cdot \vec{u}) c_{i,j}^* - \tau_{loss} c_{i,j}^* = (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs}) $$

Thus, by multiplying with the step size and summing over the sequences the previous one-sensor adjoint model and due the linearity of the left hand side (that is always true by definition of the adjoint model), $c_{i,tot}^* = \sum_j \eta_j c_i^*(s_j)$ is solution of the following trajectory-integrated one-sensor adjoint model:

$$ (69) \qquad \begin{aligned} & \partial_t c_{i,tot}^* + \nabla \cdot (\mathbf{K}^T \nabla c_{i,tot}^*) + \nabla(\vec{u} c_{i,tot}^*) - (\nabla \cdot \vec{u}) c_{i,tot}^* - \tau_{loss} c_{i,tot}^* \\ & = \sum_j (\eta_j (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs})) \end{aligned} $$

This trajectory-integrated one-sensor adjoint model enable to compute the criteria without having to compute $c_i^*(s_j)$ for each $s_j$.

Let us note that, if the one-sensor observation operator $c \mapsto m$ is linear (which is the case in the present study), the right hand side of the trajectory-integrated one-sensor adjoint model

$$ \text{RHS} = \sum_j (\eta_j (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs})) $$

$$ \text{RHS} = (d_c m_i)^* \cdot 2\mathbf{R}^{-1} \cdot (\sum_j \eta_j (m_i(c(s_j)) - m_i^{obs})). $$

Moreover, if the direct model is linear with respect to $s$ (which is the case in the present study), the right hand side of the trajectory-integrated one-sensor adjoint model even becomes

$$ \text{RHS} = (d_c m_i)^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(\sum_j \eta_j s_j)) - (\sum_j \eta_j) m_i^{obs}). $$

### A.3.2. Computation of the gradients $\nabla_s j_{reg,i}$.
In this section, we calculate the gradient of the differentiable regularization terms $j_{reg,i}$ introduced in Sec. 3.1 and in Sec. A.2. Recall that $\nabla_s$ is the gradient operator with respect to $s$ and $\nabla_s j_{reg,i}(s)$ is the function given by $\frac{d j_{reg,i}(s)}{ds} \cdot \delta s = \langle \nabla_s j_{reg,i}(s), \delta s \rangle \forall \delta s$.

**Gradient of the Tikhonov regularization term.** The gradient of the Tikhonov regularization term (23), presented in Sec. 3.1.1 is simply:

$$ \nabla_s j_{reg,T}(s) = 2\mathbf{C}^{-1}(s - s_b) $$

**Gradient of the logarithmic barrier regularization term.** The gradient of the logarithmic barrier regularization term (45), presented in Sec. A.2.1 is simply:

$$ \nabla_s j_{reg,exc}(s) = \frac{2\epsilon}{\epsilon - s^2} \mathbb{1}_{\Omega_{exc}} $$

**Gradient of the population dynamics-informed regularization term.** Recall that the population dynamics-informed regularization term (27), presented in Sec. 3.1.3, can be written as:

$$ \begin{aligned} j_{reg,PD}(s) &= \|\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}^2 \\ &= \langle \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), \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) \rangle \end{aligned} $$
Recall that the adjoint operator is a linear operator: $(F_1 + F_2)^* = F_1^* + F_2^*$

We know have to calculate the adjoint operators $(\partial_t - \frac{\partial_t q}{q})^*$, $(\frac{d\nabla F_i(s)}{ds})^*$ and $(\frac{dR_i(s)}{ds})^*$ for the $F_i$ and $R_i$ presented in Sec. 3.1.3 and in Sec. B.2.2. Let us note that no boundary conditions are imposed in the estimation of the population dynamics model. Hence, all the boundary terms are taken into account in the adjoint operators.

**Temporal dynamics terms.** Here, we calculate the operator $(\partial_t - \frac{\partial_t q}{q})^*$. By applying the scalar product with a test function $\varphi$, we have:

$$ \langle \partial_t \delta s - \frac{\partial_t q}{q} \delta s, \varphi \rangle = \langle -\partial_t \varphi - \frac{\partial_t q}{q} \varphi, \delta s \rangle + \int_{\Omega} [\varphi \delta s]_{t=0}^{t=T} dX $$

Therefore, the operator $(\partial_t - \frac{\partial_t q}{q})^*$ is the operator such that

$$ \varphi \mapsto -\partial_t \varphi - \frac{\partial_t q}{q} \varphi \quad \forall (x, y, t) \in \Omega \times ]0; T[$$

for $\varphi$ such that $\varphi = 0 \ \forall (x, y, t) \in \Omega \times \{0; T\}$

**Diffusion term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the diffusion operator $F_i = F_d$ (47). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_d(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle -\nabla \cdot (\frac{1}{q}\sigma \nabla \delta s), \varphi \rangle \\ &= \langle -\nabla \cdot (\frac{1}{q}\sigma^T \nabla \varphi), \delta s \rangle + \int_0^T \int_{\partial\Omega} \left(\frac{1}{q}\sigma^T \nabla \varphi \delta s - \frac{1}{q}\varphi \sigma \nabla \delta s\right) \cdot \vec{n} d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla F_d(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto -\nabla \cdot \left( \frac{1}{q} \sigma^T \nabla \varphi \right) \quad \forall (x, y, t) \in \Omega \times ]0; T[$$

for $\varphi$ such that $\frac{1}{q}\sigma^T \nabla \varphi \cdot \vec{n} = 0$ and $\varphi = 0 \ \forall (x, y, t) \in \partial\Omega \times ]0; T[$

**Advection term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the advection operator $F_i = F_a$ (48). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_a(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle \nabla \cdot (\frac{1}{q}\vec{v}_i \delta s), \varphi \rangle \\ &= \langle \nabla \cdot (\frac{1}{q}\vec{v}_i) \varphi - \nabla \cdot (\frac{1}{q}\vec{v}_i \varphi), \delta s \rangle + \int_0^T \int_{\partial\Omega} \frac{1}{q} \vec{v}_i \varphi \delta s \cdot \vec{n} d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla F_a(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \nabla \cdot \left( \frac{1}{q} \vec{v}_i \right) \varphi - \nabla \cdot \left( \frac{1}{q} \vec{v}_i \varphi \right) \quad \forall (x, y, t) \in \Omega \times ]0; T[$$

for $\varphi$ such that $\varphi = 0 \ \forall (x, y, t) \in \partial\Omega \times ]0; T[$
and hence, its differentiable is

$$ \frac{d j_{reg,PD}(s)}{ds} \cdot \delta s = 2 \langle \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) \cdot \partial_t \delta s - \frac{\partial_t q}{q} \delta s + q \nabla \left( \sum_i \frac{d \mathcal{F}_i(s)}{ds} \cdot \delta s \right) + q \sum_i \frac{d \mathcal{R}_i(s)}{ds} \cdot \delta s \rangle $$

Hence $\frac{d j_{reg,PD}(s)}{ds} \cdot \delta s = \langle \nabla_s j_{reg,PD}(s), \delta s \rangle$ with:

$$ \nabla_s j_{reg,PD}(s) = 2 \left( \partial_t - \frac{\partial_t q}{q} + q \sum_i \frac{d \nabla \mathcal{F}_i(s)}{ds} + q \sum_i \frac{d \mathcal{R}_i(s)}{ds} \right)^* \left( \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) \right) $$

Recall that the adjoint operator is a linear operator: $(F_1 + F_2)^* = F_1^* + F_2^*$
We know have to calculate the adjoint operators $(\partial_t - \frac{\partial_t q}{q})^*$, $(\frac{d\nabla \mathcal{F}_i(s)}{ds})^*$ and $(\frac{d\mathcal{R}_i(s)}{ds})^*$ for the $\mathcal{F}_i$ and $\mathcal{R}_i$ presented in Sec. 3.1.3 and in Sec. A.2.2. Let us note that no boundary conditions are imposed in the estimation of the population dynamics model. Hence, all the boundary terms are taken into account in the adjoint operators.

**Temporal dynamics terms.** Here, we calculate the operator $(\partial_t - \frac{\partial_t q}{q})^*$. By applying the scalar product with a test function $\varphi$, we have:

$$ \langle \partial_t \delta s - \frac{\partial_t q}{q} \delta s, \varphi \rangle = \langle -\partial_t \varphi - \frac{\partial_t q}{q} \varphi, \delta s \rangle + \int_{\Omega} [\varphi \delta s]_{t=0}^{t=T} dX $$

Therefore, the operator $(\partial_t - \frac{\partial_t q}{q})^*$ is the operator such that

$$ \varphi \mapsto -\partial_t \varphi - \frac{\partial_t q}{q} \varphi \quad \forall (x, y, t) \in \Omega \times ]0; T[ \\ \text{for } \varphi \text{ such that } \varphi = 0 \ \forall (x, y, t) \in \Omega \times \{0; T\} $$

**Diffusion term.** Here, we calculate the operator $(\frac{d\nabla \mathcal{F}_i(s)}{ds})^*$ for the diffusion operator $\mathcal{F}_i = \mathcal{F}_d$ (47). By applying the scalar product with a test function $\varphi$, we have:

$$ \langle \frac{d \nabla \mathcal{F}_d(s)}{ds} \cdot \delta s, \varphi \rangle = \langle -\nabla \cdot \left( \frac{1}{q} \sigma \nabla \delta s \right), \varphi \rangle \\ = \langle -\nabla \cdot \left( \frac{1}{q} \sigma^T \nabla \varphi \right), \delta s \rangle + \int_0^T \int_{\partial \Omega} \left( \frac{1}{q} \sigma^T \nabla \varphi \delta s - \frac{1}{q} \varphi \sigma \nabla \delta s \right) \cdot \vec{n} d\Gamma dt $$

Therefore, the operator $(\frac{d\nabla \mathcal{F}_d(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto -\nabla \cdot \left( \frac{1}{q} \sigma^T \nabla \varphi \right) \quad \forall (x, y, t) \in \Omega \times ]0; T[ \\ \text{for } \varphi \text{ such that } \frac{1}{q} \sigma^T \nabla \varphi \cdot \vec{n} = 0 \text{ and } \varphi = 0 \ \forall (x, y, t) \in \partial \Omega \times ]0; T[ $$

**Advection term.** Here, we calculate the operator $(\frac{d\nabla \mathcal{F}_i(s)}{ds})^*$ for the advection operator $\mathcal{F}_i = \mathcal{F}_a$ (48). By applying the scalar product with a test function $\varphi$, we have:

$$ \langle \frac{d \nabla \mathcal{F}_a(s)}{ds} \cdot \delta s, \varphi \rangle = \langle \nabla \cdot \left( \frac{1}{q} \vec{v}_i \delta s \right), \varphi \rangle \\ = \langle \nabla \cdot \left( \frac{1}{q} \vec{v}_i \right) \varphi - \nabla \cdot \left( \frac{1}{q} \vec{v}_i \varphi \right), \delta s \rangle + \int_0^T \int_{\partial \Omega} \frac{1}{q} \vec{v}_i \varphi \delta s \cdot \vec{n} d\Gamma dt $$
**Quorum sensing term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the quorum sensing operator $\mathcal{F}_i = \mathcal{F}_{qs}$ (51). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \left\langle \frac{d\nabla \mathcal{F}_{qs}(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle -\nabla \cdot \left( \frac{\sigma_{qs}}{q^2} (s\nabla\delta s + \nabla s\delta s) \right), \varphi \right\rangle \\ &= \left\langle -\nabla \cdot \left( \frac{\sigma_{qs}}{q^2} (s\nabla\varphi - \nabla s\varphi) \right) - \nabla \cdot \left( \frac{\sigma_{qs}}{q^2} \nabla s \right) \varphi, \delta s \right\rangle \\ &\quad + \int_0^T \int_{\partial\Omega} \frac{\sigma_{qs}}{q^2} (s\nabla\varphi\delta s - s\varphi\nabla\delta s + \nabla s\varphi\delta s) \cdot \vec{n} d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla \mathcal{F}_{qs}(s)}{ds})^*$ is the operator such that

$$ \begin{aligned} \varphi \mapsto & -\nabla \cdot \left( \frac{\sigma_{qs}}{q^2} (s\nabla\varphi - \nabla s\varphi) \right) - \nabla \cdot \left( \frac{\sigma_{qs}}{q^2} \nabla s \right) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] \\ \text{for } \varphi \text{ such that } \varphi = 0 \text{ and } \frac{\sigma_{qs}}{q^2} s\nabla\varphi \cdot \vec{n} = 0 & \quad \forall (x, y, t) \in \partial\Omega \times [0; T] \end{aligned} $$

**Non local pheromone sensing term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the non-local pheromone sensing operator $\mathcal{F}_i = \mathcal{F}_{nlps}$ (52). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \left\langle \frac{d\nabla \mathcal{F}_{nlps}(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle \nabla \cdot \left( \frac{\chi}{q} (s\nabla(g_{gp} * \delta s) + \nabla(g_{gp} * s)\delta s) \right), \varphi \right\rangle \\ &= \left\langle \bar{g}_{pg} * \nabla \cdot \left( \frac{\chi}{q} s\nabla\varphi \right) - \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s)\varphi \right) + \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s) \right) \varphi, \delta s \right\rangle \\ &\quad + \int_0^T \int_{\partial\Omega} \frac{\chi}{q} (s\varphi\nabla(g_{gp} * \delta s) - s\nabla\varphi(g_{gp} * \delta s) + \nabla(g_{gp} * s)\varphi\delta s) \cdot \vec{n} d\Gamma dt \end{aligned} $$

with $\bar{g}_{pg}$ the function such that $\bar{g}_{pg}(x, y, t) = g_{pg}(-x, -y, -t)$. Therefore, the operator $(\frac{d\nabla \mathcal{F}_{nlps}(s)}{ds})^*$ is the operator such that

$$ \begin{aligned} \varphi \mapsto & \bar{g}_{pg} * \nabla \cdot \left( \frac{\chi}{q} s \nabla \varphi \right) - \nabla \cdot \left( \frac{\chi}{q} \nabla (g_{gp} * s) \varphi \right) + \nabla \cdot \left( \frac{\chi}{q} \nabla (g_{gp} * s) \right) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] \\ \text{for } \varphi \text{ such that } \varphi = 0 \text{ and } \frac{\chi}{q} s \nabla \varphi \cdot \vec{n} = 0 & \quad \forall (x, y, t) \in \partial\Omega \times [0; T] \end{aligned} $$

**Logistic term.** Here, we calculate the operator $(\frac{d\mathcal{R}_i(s)}{ds})^*$ for the logistic operator $\mathcal{R}_i = \mathcal{R}_I$ (53). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \left\langle \frac{d\mathcal{R}_I(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \delta s, \varphi \right\rangle \\ &= \left\langle \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \varphi, \delta s \right\rangle \end{aligned} $$

Therefore, the operator $(\frac{d\mathcal{R}_I(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$

**Allee effect term.** Here, we calculate the operator $(\frac{d\mathcal{R}_I(s)}{ds})^*$ for the Allee effect operator $\mathcal{R}_i = \mathcal{R}_{Ae}$ (54). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \left\langle \frac{d\mathcal{R}_{Ae}(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2\kappa_+\kappa_-)\delta s, \varphi \right\rangle \\ &= \left\langle \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2\kappa_+\kappa_-)\varphi, \delta s \right\rangle \end{aligned} $$
Therefore, the operator $(\frac{d\nabla F_a(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \nabla \cdot \left(\frac{1}{q} \vec{v}_i\right) \varphi - \nabla \cdot \left(\frac{1}{q} \vec{v}_i \varphi\right) \quad \forall (x, y, t) \in \Omega \times [0; T] $$
for $\varphi$ such that $\varphi = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Quorum sensing term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the quorum sensing operator $F_i = F_{qs}$ (51). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_{qs}(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle -\nabla \cdot \left( \frac{\sigma_{qs}}{q^2} (s\nabla\delta s + \nabla s\delta s) \right), \varphi \rangle \\ &= \langle -\nabla \cdot \left( \frac{\sigma_{qs}}{q^2} (s\nabla\varphi - \nabla s\varphi) \right) - \nabla \cdot \left( \frac{\sigma_{qs}}{q^2} \nabla s \right) \varphi, \delta s \rangle \\ &\quad + \int_0^T \int_{\partial\Omega} \frac{\sigma_{qs}}{q^2} (s\nabla\varphi\delta s - s\varphi\nabla\delta s + \nabla s\varphi\delta s) \cdot \vec{n} d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla F_{qs}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto -\nabla \cdot \left( \frac{\sigma_{qs}}{q^2} (s \nabla \varphi - \nabla s \varphi) \right) - \nabla \cdot \left( \frac{\sigma_{qs}}{q^2} \nabla s \right) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$
for $\varphi$ such that $\varphi = 0$ and $\frac{\sigma_{qs}}{q^2} s \nabla \varphi \cdot \vec{n} = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Non local pheromone sensing term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the non-local pheromone sensing operator $F_i = F_{nlps}$ (52). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_{nlps}(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle \nabla \cdot \left( \frac{\chi}{q} (s\nabla(g_{gp} * \delta s) + \nabla(g_{gp} * s)\delta s) \right), \varphi \rangle \\ &= \langle \bar{g}_{pg} * \nabla \cdot \left( \frac{\chi}{q} s \nabla \varphi \right) - \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s)\varphi \right) + \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s) \right) \varphi, \delta s \rangle \\ &\quad + \int_0^T \int_{\partial\Omega} \frac{\chi}{q} (s\varphi\nabla(g_{gp} * \delta s) - s\nabla\varphi(g_{gp} * \delta s) + \nabla(g_{gp} * s)\varphi\delta s) \cdot \vec{n} d\Gamma dt \end{aligned} $$

with $\bar{g}_{pg}$ the function such that $\bar{g}_{pg}(x, y, t) = g_{pg}(-x, -y, -t)$. Therefore, the operator $(\frac{d\nabla F_{nlps}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \bar{g}_{pg} * \nabla \cdot \left( \frac{\chi}{q} s \nabla \varphi \right) - \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s)\varphi \right) + \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s) \right) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$
for $\varphi$ such that $\varphi = 0$ and $\frac{\chi}{q} s \nabla \varphi \cdot \vec{n} = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Logistic term.** Here, we calculate the operator $(\frac{dR_i(s)}{ds})^*$ for the logistic operator $R_i = R_I$ (53). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{dR_i(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \delta s, \varphi \rangle \\ &= \langle \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \varphi, \delta s \rangle \end{aligned} $$

Therefore, the operator $(\frac{dR_I(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$
Therefore, the operator $(\frac{dR_{Ae}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2 \kappa_+ \kappa_-) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$

### B.3.3. On the proximal gradient algorithm for the LASSO regularization term.
As mentioned in Sec. 3.1, if the regularization terms $j_{reg}$ are convex and differentiable, then the cost $j$ is convex and differentiable and the optimization problem (3) is solved using a gradient descent algorithm. However, as mentioned in Secs 3.1 and 3.1.2, the differentiable framework can be extended to a non-smooth and subdifferentiable framework to consider LASSO-type regularization terms. In this context, the gradient descent algorithm can be extended to the proximal-gradient algorithm, also called the Forward-Backward algorithm.

In this section, the proximal-gradient algorithm is briefly presented in the context of the optimization problem (3) with several $j_{reg}$ including the LASSO $j_{reg}$ (24) and group-LASSO $j_{reg}$ (25).

The proximal-gradient algorithm, see e.g. (Parikh et al., 2014) is an extension of the gradient descent algorithm used when the function to minimize can be split in $j(s) = g(s) + h(s)$ with:

* $g(s)$ a convex and differentiable function, which is in the present context the differentiable part of $j$, i.e. the observation term $j_{obs}$ and the differentiable $j_{reg}$,

* $h(s)$ a convex and subdifferentiable function, which is in the present case the LASSO or group-LASSO $j_{reg}$.

For such cost function, the optimum $s_a$ is a fixed point $s_a = \text{prox}_{\lambda h}(s_a - \lambda \nabla g(s_a))$ with $\text{prox}_f$ the proximal operator applied to the function $f$ and $\lambda > 0$ the gradient descent step size.

The proximal-gradient algorithm is then a fixed point algorithm based on the iteration loop:

$$ s_{k+1} = \text{prox}_{\lambda h}(s_k - \lambda \nabla g(s_k)) $$

In the present case, the gradient $\nabla g$ is computed using the expression (7). For the LASSO regularization term $h = \alpha_{reg,LASSO}j_{reg,LASSO}$ given by Eq. (24), the proximal operator $\text{prox}_{\lambda h}$ is

$$ \text{prox}_{\lambda h}(s) = \text{sign}(s) \max(0, |s| - \lambda \alpha_{reg,LASSO}). $$

For the group-LASSO regularization term $h = \alpha_{reg,gLASSO}j_{reg,gLASSO}$ given by Eq. (25), the proximal operator $\text{prox}_{\lambda h}$ is

$$ \text{prox}_{\lambda h}(s) = \text{sign}(s) \max \left( 0, 1 - \frac{\lambda \alpha_{reg,gLASSO}}{\|s\|_{L^2(0;T)}} \right) s $$

with $\|s\|_{L^2(0;T)} = (\int_0^T s(t)^2 dt)^{1/2}$. As highlighted by the proximal operator $\text{prox}_{\lambda h}$, this algorithm acts as an iterative soft-threshold, shrinking to 0 the value of $s$ lower than the threshold defined by $\lambda \alpha_{reg,LASSO}$ for LASSO and the value of $s$ at places where the total amount $\|s\|_{L^2(0;T)}$ is lower than the threshold $\lambda \alpha_{reg,gLASSO}$ for group-LASSO. Consequently, the adaptation of the proximal-gradient algorithm to the LASSO regularization term is commonly called Iterative Shrinkage-Threshold Algorithm (ISTA).

## B.4. Comparison of the LASSO and group-LASSO regularization terms in the Toy case

In this section, we briefly compare the LASSO (24) and group-LASSO (25) regularization terms in the Toy case presented in Sec. 2.3.2. Recall that the effect of the LASSO regularization term in the Toy case has already been highlighted in Sec. 3.4.
Allee effect term. Here, we calculate the operator $(\frac{dR_i(s)}{ds})^*$ for the Allee effect operator $R_i = R_{Ae}$ (54). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{dR_{Ae}(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2 \kappa_+ \kappa_-) \delta s, \varphi \rangle \\ &= \langle \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2 \kappa_+ \kappa_-) \varphi, \delta s \rangle \end{aligned} $$

Therefore, the operator $(\frac{dR_{Ae}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2 \kappa_+ \kappa_-) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$

### A.3.3. On the proximal gradient algorithm for the LASSO regularization term.
As mentioned in Sec. 3.1, if the regularization terms $j_{reg}$ are convex and differentiable, then the cost $j$ is convex and differentiable and the optimization problem (3) is solved using a gradient descent algorithm. However, as mentioned in Secs 3.1 and 3.1.2, the differentiable framework can be extended to a non-smooth and subdifferentiable framework to consider LASSO-type regularization terms. In this context, the gradient descent algorithm can be extended to the proximal-gradient algorithm, also called the Forward-Backward algorithm.

In this section, the proximal-gradient algorithm is briefly presented in the context of the optimization problem (3) with several $j_{reg}$ including the LASSO $j_{reg}$ (24) and group-LASSO $j_{reg}$ (25).

The proximal-gradient algorithm, see e.g. (Parikh et al., 2014) is an extension of the gradient descent algorithm used when the function to minimize can be split in $j(s) = g(s) + h(s)$ with:

* $g(s)$ a convex and differentiable function, which is in the present context the differentiable part of $j$, i.e. the observation term $j_{obs}$ and the differentiable $j_{reg}$,

* $h(s)$ a convex and subdifferentiable function, which is in the present case the LASSO or group-LASSO $j_{reg}$.

For such cost function, the optimum $s_a$ is a fixed point $s_a = \text{prox}_{\lambda h}(s_a - \lambda \nabla g(s_a))$ with $\text{prox}_f$ the proximal operator applied to the function $f$ and $\lambda > 0$ the gradient descent step size.

The proximal-gradient algorithm is then a fixed point algorithm based on the iteration loop:

$$ s_{k+1} = \text{prox}_{\lambda h}(s_k - \lambda \nabla g(s_k)) $$

In the present case, the gradient $\nabla g$ is computed using the expression (7). For the LASSO regularization term $h = \alpha_{reg,LASSO} j_{reg,LASSO}$ given by Eq. (24), the proximal operator $\text{prox}_{\lambda h}$ is

$$ \text{prox}_{\lambda h}(s) = \text{sign}(s) \max(0, |s| - \lambda \alpha_{reg,LASSO}). $$

For the group-LASSO regularization term $h = \alpha_{reg,gLASSO} j_{reg,gLASSO}$ given by Eq. (25), the proximal operator $\text{prox}_{\lambda h}$ is

$$ \text{prox}_{\lambda h}(s) = \text{sign}(s) \max \left( 0, 1 - \frac{\lambda \alpha_{reg,gLASSO}}{\|s\|_{L^2(0;T)}} \right) s $$

with $\|s\|_{L^2(0;T)} = (\int_0^T s(t)^2 dt)^{1/2}$. As highlighted by the proximal operator $\text{prox}_{\lambda h}$, this algorithm acts as an iterative soft-threshold, shrinking to 0 the value of $s$ lower than the threshold defined by $\lambda \alpha_{reg,LASSO}$ for LASSO and the value of $s$ at places where the total amount $\|s\|_{L^2(0;T)}$ is lower than the threshold $\lambda \alpha_{reg,gLASSO}$ for group-LASSO. Consequently, the adaptation of the proximal-gradient algorithm to the LASSO regularization term is commonly called *Iterative Shrinkage-Threshold Algorithm (ISTA)*.
<figure><img src="image_21.png" /><figcaption>Figure 18 - Optimal sources estimated with different regularization terms. Optimal quantity of pheromone emitted $s_a$ obtained by inference without regularization (left), with the LASSO regularization (24) (centre) and with the group-LASSO regularization (25) (right) for the Toy case. The optimal sources are indicated at $t = 0.0s$ (upper row) and $t = 20.0s$ (lower row).</figcaption></figure>

### B.4.1. Visual comparison of estimated pheromone emissions $s_a$.
Recall that the LASSO regularization ('LASSO', centre of Fig. 18) implies that the retro-plume is less diffuse at both initial ($t = 0s$, top row) and final ($t = 20s$, bottom row) times. On the other hand, the group-LASSO ('group LASSO', right of Fig. 18) implies that the retro-plume is even less diffuse at initial time but more diffuse at final time. This is due to the fact that the group-LASSO regularization aims to threshold all the pheromone emissions at location where the total emissions are too low, unlike the LASSO regularization that shrinks all pheromone emissions that are too low. To compensate these retro-plumes less diffuse at initial time but more diffuse at final time, the pheromone emissions are higher at initial time but lower at final time with the group-LASSO regularization than with the LASSO regularization. Let us note that more spurious emissions far from the target support are recovered with the group-LASSO regularization than with the LASSO regularization.

### B.4.2. RMSE criteria.
Recall that with the LASSO regularization (pink curves on Fig. 19), the total error $RxMSE_{tot}$ is sensibly lower and the error on the target support $RxMSE_{supp}$ similar to the ones obtained without regularization (blue curve) at the beginning of the time window and up to $t \approx 17.5s$. After that time, with the LASSO regularization, both $RxMSE_{tot}$ and $RxMSE_{supp}$ significantly increase and become higher than without regularization. With the group-LASSO
<figure><img src="image_21.png" /><figcaption>Figure 18 - Optimal sources estimated with different regularization terms. Optimal quantity of pheromone emitted $s_a$ obtained by inference without regularization (left), with the LASSO regularization (24) (centre) and with the group-LASSO regularization (25) (right) for the Toy case. The optimal sources are indicated at $t = 0.0s$ (upper row) and $t = 20.0s$ (lower row).</figcaption></figure>

### A.4. Comparison of the LASSO and group-LASSO regularization terms in the Toy case

In this section, we briefly compare the LASSO (24) and group-LASSO (25) regularization terms in the Toy case presented in Sec. 2.3.2. Recall that the effect of the LASSO regularization term in the Toy case has already been highlighted in Sec. 3.4.

#### A.4.1. Visual comparison of estimated pheromone emissions $s_a$.
Recall that the LASSO regularization ('LASSO', centre of Fig. 18) implies that the retro-plume is less diffuse at both initial ($t = 0s$, top row) and final ($t = 20s$, bottom row) times. On the other hand, the group-LASSO ('group LASSO', right of Fig. 18) implies that the retro-plume is even less diffuse at initial time but more diffuse at final time. This is due to the fact that the group-LASSO regularization aims to threshold all the pheromone emissions at location where the total emissions are too low, unlike the LASSO regularization that shrinks all pheromone emissions that are too low. To compensate these retro-plumes less diffuse at initial time but more diffuse at final time, the pheromone emissions are higher at initial time but lower at final time with the group-LASSO regularization than with the LASSO regularization. Let us note that more spurious emissions far from the target support are recovered with the group-LASSO regularization than with the LASSO regularization.
<figure><img src="image_22.png" /><figcaption>**Figure 19** – 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 LASSO regularization term (24) (LASSO, pink) and with the group-LASSO regularization term (25) (group LASSO, yellow)</figcaption></figure>

regularization (yellow curves on Fig. 19), the errors follow similar trend but are nearly always slightly lower. Let us note that, at the very beginning of the time window (between $t = 0$s and $t = 2.5$s), $R \times \text{MSE}_{\text{supp}}$ is higher with group-LASSO regularization than with LASSO regularization. This comes from the fact that the group-LASSO regularization thresholds all the emissions, including at the initial time, located at certain positions in the target supports.

B.4.3. *Presence prediction, false discovery and miss rates.* As already anticipated with the visual comparison, the presence prediction (upper row of Fig. 20) is less diffuse at initial time, similarly diffuse at middle time and more diffuse at final time with the group-LASSO regularization (yellow curves) than with the LASSO regularization (pink curve). Recall the with the LASSO regularization (pink curves), the miss rate $\rho_{fa}$ and false discovery rate $\rho_{fp}$ (bottom row of Fig. 20) are quite stable until $t \approx 17.5$s, time at which $\rho_{fa}$ significantly increases. With the group-LASSO regularization (yellow curves), the trend is very similar, but two little difference can be noted. First, as observed on the errors, at the initial times, $\rho_{fa}$ is slightly higher and $\rho_{fp}$ is slightly lower with the group-LASSO regularization than with the LASSO regularization. With the LASSO regularization, $\rho_{fa}$ is even null for a short amount of time, while it is always non-zero with the group-LASSO regularization. Also, one can note that $\rho_{fa}$ significantly increases a bit later with the group-LASSO regularization than with the LASSO regularization.
<figure><img src="image_22.png" /><figcaption>**Figure 19** – 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 ($RxMSE_{tot}$, Eq. (14), upper panel) and on the target support ($RxMSE_{supp}$, Eq. (15), lower panel) in the Toy case: with no regularization term (no reg, blue), with the LASSO regularization term (24) (LASSO, pink) and with the group-LASSO regularization term (25) (group LASSO, yellow)</figcaption></figure>

### A.4.2. RMSE criteria.
Recall that with the LASSO regularization (pink curves on Fig. 19), the total error $RxMSE_{tot}$ is sensibly lower and the error on the target support $RxMSE_{supp}$ similar to the ones obtained without regularization (blue curve) at the beginning of the time window and up to $t \approx 17.5$s. After that time, with the LASSO regularization, both $RxMSE_{tot}$ and $RxMSE_{supp}$ significantly increase and become higher than without regularization. With the group-LASSO regularization (yellow curves on Fig. 19), the errors follow similar trend but are nearly always slightly lower. Let us note that, at the very beginning of the time window (between $t = 0$s and $t = 2.5$s), $RxMSE_{supp}$ is higher with group-LASSO regularization than with LASSO regularization. This comes from the fact that the group-LASSO regularization thresholds all the emissions, including at the initial time, located at certain positions in the target supports.

### A.4.3. Presence prediction, false discovery and miss rates.
As already anticipated with the visual comparison, the presence prediction (upper row of Fig. 20) is less diffuse at initial time, similarly diffuse at middle time and more diffuse at final time with the group-LASSO regularization (yellow curves) than with the LASSO regularization (pink curve). Recall the with the LASSO regularization (pink curves), the miss rate $\rho_{fa}$ and false discovery rate $\rho_{fp}$ (bottom row of Fig. 20) are quite stable until $t \approx 17.5$s, time at which $\rho_{fa}$ significantly increases. With the group-LASSO regularization (yellow curves), the trend is very similar, but two little difference can be noted. First, as observed on the errors, at the initial times, $\rho_{fa}$ is slightly higher and $\rho_{fp}$ is slightly lower with the group-LASSO regularization than with the LASSO regularization. With the LASSO regularization, $\rho_{fa}$ is even null for a short amount of time, while it is always non-zero with the group-LASSO regularization. Also, one can note that $\rho_{fa}$ significantly increases a bit later with the group-LASSO regularization than with the LASSO regularization.
<figure><img src="image_23.png" /><figcaption>Figure 20 – 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 LASSO regularization term (24) (LASSO, pink) and with the group-LASSO regularization term (25) (group LASSO, yellow). 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.</figcaption></figure>

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_23.png" /><figcaption>Figure 20 – 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 LASSO regularization term (24) (LASSO, pink) and with the group-LASSO regularization term (25) (group LASSO, yellow). 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.</figcaption></figure>

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).