# Metamodeling of a physically based pesticide runoff model with a long short-term memory approach

Guillaume Métayer, Cécile Dagès, Marc Voltz, Jean-Stéphane Bailly

► To cite this version:

Guillaume Métayer, Cécile Dagès, Marc Voltz, Jean-Stéphane Bailly. Metamodeling of a physically based pesticide runoff model with a long short-term memory approach. Journal of Hydrology, 2025, 662 (Part A), pp.133800. 10.1016/j.jhydrol.2025.133800. hal-05192397

HAL Id: hal-05192397

https://hal.inrae.fr/hal-05192397v1

Submitted on 30 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Research papers

# Metamodeling of a physically based pesticide runoff model with a long short-term memory approach

Guillaume Métayer* ID, Cécile Dagès, Marc Voltz, Jean-Stéphane Bailly

UMR LISAH, University of Montpellier, INRAE, AgroParisTech, Institut Agro Montpellier, IRD, Montpellier, France

## ARTICLE INFO

This manuscript was handled by D. Han, Editor-in-Chief, with the assistance of Bo Guo, Associate Editor

Keywords:

Surrogate model
LSTM
Concentration
Prediction
Surface water

## ABSTRACT

Predicting the temporal variability in pesticide runoff from agricultural fields is essential for assessing the impacts of agricultural practices on surface water quality and the associated environmental risks. Richards-based models are valuable tools for this purpose, particularly in contexts characterized by high temporal variability in rainfall intensity, such as in the Mediterranean region. These methods can be applied to a wide range of operational applications at both the field and catchment scales. However, their operational capability is limited, primarily due to their high computation times. In this study, we assessed a metamodeling approach for a Richards-based model designed to simulate the hourly variability in water and pesticide runoff over a year while reducing computation times. The proposed approach is based on long short-term memory (LSTM) models and is applied to a field-scale pesticide runoff model coupling Richards and convection-dispersion equations with the uniform mixing cell concept for pesticide remobilization and an overland flow routine. The resulting metamodel accounts for meteorological conditions, compound properties, and pesticide application dates and rates. Although the LSTM architecture is not optimized, the metamodel reduces the computation time by 90% compared with that of the initial model. The metamodel shows high overall performance in simulating hourly discharges and pesticide concentrations, with Nash–Sutcliffe efficiency values of 0.99, confirming the ability of the LSTM models to handle hydrochemical time series. The results nonetheless indicate room for improvement in specific areas, such as simulating runoff occurrence or concentrations for compounds with low half-life values. Furthermore, the results encourage further consideration of adding static parameters as inputs to LSTM models to increase their generalizability. Overall, this study provides key insights for deploying methods to assess the impacts of agricultural practices on a large scale.

## 1. Introduction

The use of pesticides in agriculture can affect surface water quality and, consequently, both environmental and human health. A major pathway for the transfer of pesticides from agricultural fields to surface waters is surface runoff (Wauchope et al., 1995; Louchart et al., 2001; Reichenberger et al., 2007), which is influenced by many factors that vary in space and time (Shipitalo and Owens, 2003; Calvet et al., 2005), such as the type of compound, topography, climate, soil, and agricultural practices. Numerical models that simulate pesticide transfer through surface runoff over time are valuable tools for assessing the impact of pesticide application on surface water quality over time under varying climatic conditions and soil characteristics.

Various types of models can predict the temporal variability of pesticide runoff from an agricultural field, differing primarily in their

methods for representing the water runoff process (Siimes and Kämäri, 2003). Some models estimate daily runoff by applying the empirical Soil Conservation Service (SCS) curve-number (CN) method (e.g., the pesticide root zone model (PRZM)—Carsel et al., 1985; Groundwater Loading Effects of Agricultural Management Systems (GLEAMS)—Leonard et al., 1987; Agricultural Policy/Environmental eXtender (APEX)—Williams et al., 2006). This simplified representation of runoff processes is computationally efficient, thereby facilitating model application. Consequently, these models are used for operational applications, such as assessing environmental exposure to pesticides in water bodies as part of the regulatory evaluation of environmental risks associated with pesticide use in the European Union (EU) and the United States (US) (FOCUS, 2001; Boivin and Poulsen, 2017; Young, 2019). However, these models face limitations in specific meteorological contexts (Sreejith et al., 2024), such as Mediterranean conditions characterized by stormy

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: guillaume.metayer10@gmail.com (G. Métayer).</footnote>
precipitation with high intensity and significant temporal variability. Accurately simulating Hortonian runoff in these contexts requires sub-daily temporal resolution (Bronstert et al., 2023). Richards-based models can simulate surface runoff as soil infiltration excess at vari- able subdaily time steps. They are thus well suited for addressing these challenges. They are capable of accurately simulating pesticide runoff over multiple years (Métayer et al., 2024). Nonetheless, these modeling approaches, which rely on the numerical solution of differential equa- tions, require expertise for their implementation, numerous soil pa- rameters and a long computing time. Thus, Richards-based models are primarily used for research purposes. These methods are impractical for operational applications requiring thousands of model executions, such as the regulatory evaluation of environmental risks associated with pesticide use in the EU and US.

Metamodeling approaches, also known as surrogate modeling or
emulation modeling (Razavi et al., 2012), provide a means to reduce the
computation time of an initial physically based model while maintaining
similar predictive capability. These approaches aim to fit an approxi-
mation mathematical function to reproduce outputs from an initial
model, with reduced computation time (Kleijnen, 2009; Castelletti et al.,
2012). Thus, metamodeling relies on the generation of a training dataset
composed of simulations from the initial model (inputs/outputs of in-
terest) that cover the intended application domain and on the selection
of a mathematical approximation function. Two main metamodeling
approaches can be distinguished on the basis of the nature of the chosen
approximation function (Razavi et al., 2012). On the one hand, process-
based approaches involve simplifying the mechanisms of the initial
physically based model (e.g., Tiktak et al., 2006). As such, they are
conceptual simplifications driven by the modeler. On the other hand,
statistically based approaches replace the initial model with a mathe-
matical function that has no explicit physical meaning. Unlike process-
based metamodeling, it does not rely on a priori knowledge. The
resulting metamodel has no explicit physical meaning (Lim and Wang,
2022) but implicitly accounts for all the mechanisms represented by the
initial model (Stenemo et al., 2007). Thus, statistically based approaches
are appropriate for metamodeling physically based models when pre-
serving their full physical complexity, which is relevant for a specific
application, such as Richards-based models for runoff simulation under
Mediterranean conditions.

In the field of environmental sciences, several statistically based approaches for metamodeling have already been presented and tested in the literature. These methods can be broadly divided into two main categories: (a) methods for interpolating the response surface of the initial model from a small subset of outputs from the initial models (e.g., Gaussian processes; Marrel and Iooss, 2023) and (b) methods approximating response surfaces of the initial model from any input n-uplet (e. g., deep learning; Cohen, 2023). These methods have been applied for metamodeling Richards-based pesticide fate models to simulate time-average pesticide fate quantities. For example, Reichenberger et al. (2021) and Stenemo et al. (2007) developed metamodels of the MACRO model (Larsbo et al., 2005) to simulate the average pesticide concentration in soil at depths of 2 m over 20 years and the 80th percentile of the average concentration at a depth of 1 m over the same period, respectively. Reichenberger et al. (2021) applied trilinear interpolation in the three-dimensional space of the sorption coefficient of the compound, the compound's half-life, and the Freundlich exponent, whereas Stenemo et al. (2007) employed a perceptron neural network. To our knowledge, no metamodeling method has been tested for simulating the temporal variation in surface or groundwater contamination by pesticides.

Among statistically based approaches for approximating response
surfaces, deep learning models are complex nonlinear systems built by
extracting essential features or representations directly from the raw
input data by organizing their computational processing units, called
artificial neurons, into a hierarchy of layers (Alvi et al., 2023). Deep
learning models have generated significant enthusiasm in hydrology in

recent years and have demonstrated considerable potential for this
purpose. Specifically, long short-term memory neural networks (LSTM;
Hochreiter and Schmidhuber, 1997) have recently been identified as
particularly effective for predicting runoff time series (Kratzert et al.,
2018, 2019a; b; Gauch et al., 2021; Frame et al., 2021; Nearing et al.,
2022) and, for this reason, are currently among the most studied
modeling methods (Tripathy and Mishra, 2024). Several studies have
also demonstrated their potential to simulate temporal variations in
variables related to water quality, such as sediment loads (e.g., AlDahoul
et al., 2021) or chlorophyll-a concentrations (e.g., Abbas et al., 2023).
LSTM neural networks have primarily been assessed for the prediction of
daily flow or water quality variables on the basis of observed data at the
catchment scale. Daily predictions have been assessed over multiple
years (Kratzert et al., 2018, 2019a, 2019b; AlDahoul et al., 2021; Frame
et al., 2021; Lees et al., 2021; Yin et al., 2021; Abbas et al., 2023;
Gholizadeh et al., 2023; Liu et al., 2024). Their measured prediction
performance is very good according to the performance thresholds
recommended by Moriasi et al. (2015), with Nash–Sutcliffe efficiency
(NSE) values often exceeding 0.7. Similarly, simulations at an hourly
time step have demonstrated excellent results, with NSE values
exceeding 0.8 (Kao et al., 2020; Xiang et al., 2020). In these cases,
performance was evaluated over relatively short periods, typically
spanning from just a few hours to one day. On the basis of these findings,
we hypothesize that LSTM neural networks can effectively simulate
subdaily variations in water and pesticide runoff, positioning them as a
promising metamodeling approach for physically based pesticide runoff
models.

Thus, the main aim of the present work was to develop and assess a
metamodeling approach for a physically based field-scale model
designed to simulate the hourly variabilities in water and pesticide
runoff via an LSTM model. The specific objectives were to i) determine
whether an LSTM model can effectively metamodel a Richards-based
pesticide fate model for simulating hourly series of Hortonian runoff
discharges and associated pesticide concentrations on an annual scale
for a given Mediterranean agricultural field and ii) identify the strengths
and limits of the LSTM model for achieving this goal. In this work, the
model subjected to metamodeling estimates pesticide surface runoff by
simulating pesticide transport within the soil using the Richards equa-
tion, coupled with a convection-dispersion equation, pesticide remobi-
lization in surface runoff using the uniform mixing cell concept and
overland flow at the field scale using the diffusive wave approach. This
model is implemented in the fully spatially distributed hydrological
model MHYDAS-Pesticide (Crevoisier et al., 2021). Hereafter, MHYDAS-
Pesticide refers specifically to its field-scale implementation.

## 2. Materials and methods

The metamodeling of MHYDAS-Pesticide is based on fitting an
approximation function using a dataset consisting of MHYDAS-Pesticide
simulations (inputs/outputs). Version 1.0 of MHYDAS-Pesticide
(Crevoisier et al., 2021) was used for this study. This section first pro-
vides an overview of MHYDAS-Pesticide, followed by a description of
the approximation functions tested and the experimental design used to
generate the dataset for fitting and performance evaluation. It then de-
tails the method used to fit the approximation functions. Finally, the
approach for assessing the performance of the resulting metamodel in
simulating hourly runoff discharges and pesticide concentrations over
one year is presented, and the reduction in computation time compared
with MHYDAS-Pesticide is evaluated.

2.1. MHYDAS-pesticide description

MHYDAS-Pesticide is a fully spatially distributed physically based
model that simulates, over a continuous period, the hydrological fate of
a pesticide compound in a heterogeneous agricultural landscape. In this
work, only its simulation capabilities at the field scale have been
<table><caption>**Table 1**
Number of LSTM layers and their hidden state length for the two tested LSTM networks.</caption><thead><tr><td>LSTM layer</td><td>Hidden state length</td></tr></thead><tbody><tr><td>1</td><td>64</td></tr><tr><td>2</td><td>128</td></tr><tr><td>3</td><td>128</td></tr><tr><td>4</td><td>64</td></tr><tr><td>5</td><td>1</td></tr></tbody></table>

mobilized. A brief description of the numerical representation of the main physicochemical processes involved in this application is presented below. A detailed description of the processes represented by MHYDAS-Pesticide is available in the model manual (Crevoisier et al., 2021).

With respect to hydrological processes, the model simulates infiltration, water movement in the soil, surface runoff, and water propagation. Excess infiltration and water movement in soil are simulated by solving the Richards equation (Richards, 1931) via the Ross numerical method (Ross, 2003; Crevoisier et al., 2009). The relationships among the water content, pressure head and hydraulic conductivity are expressed by the modified Mualem–van Genuchten hydrodynamic model (Schaap and van Genuchten, 2006). Surface runoff results from the simulated infiltration excess once the maximum ponding height is exceeded. Then, surface runoff is propagated to the outlet of the field by overland flow via the diffusive wave model, a simplification of the Saint-Venant system (Saint-Venant, 1871), which is solved via the Hayami method (Moussa, 1996).

With respect to pesticide fate processes, the model simulates the vertical reactive transport of compounds in soil and transport by runoff. The former is represented by a convection-dispersion equation. An instantaneous sorption equilibrium relationship is assumed between the liquid and solid soil phase concentrations; for this study, a linear model is used. Compound degradation is represented by first-order kinetics for both the liquid and solid soil phases. The degradation rate is adjusted by a coefficient depending on the soil water content (Walker, 1974). For calculating concentrations in runoff, uniform and instantaneous mixing between runoff and the topsoil mixing cell is considered. For this work,

the latter topsoil mixing cell is set to 1 cm.

The target outputs were runoff discharge and the associated pesticide concentrations. The model time step is variable, ranging from a fraction of seconds to 5 min during rainfall events and up to an hour otherwise, adjusting according to water fluxes in the soil.

## 2.2. Tested approximation functions

The chosen approximation functions are LSTM neural networks. The general principles of these models are explained in the supplementary materials (Text S1). Two LSTM-based sequence-to-sequence networks were fitted to the outputs of MHYDAS-Pesticide: i) LSTM-runoff for the simulation of hourly runoff discharges and ii) LSTM-conc for simulating hourly pesticide concentrations.

The two networks share the same LSTM-based sequence-to-sequence architecture (Table 1). The optimal architecture for simulating a variable depends on the objective and the data to be processed (Bengio, 2012). The tested architecture was chosen on the basis of results from studies presented in the literature. Several LSTM layers were arranged in series. Salman et al. (2018) demonstrated that a multilayer LSTM network with 4 layers in series performs better than a network with only a single LSTM layer for processing meteorological data. Therefore, the networks were configured with 5 LSTM layers in series (Table 1), with the last layer modified to simulate the target variable (hourly runoff discharge or concentration). This final layer, accordingly, has a hidden state length of 1 and uses the ReLU activation function instead of the tanh function to calculate the output value and the potential update vector for $ct$ (Text S1). Since the outputs of an LSTM layer range between -1 and 1 when the tanh function is used, the use of the ReLU function allows for the calculation of higher values for the target variables. This function has been used in several examples of discharge or water quality calculations (e.g., Xiang et al., 2020; AlDahoul et al., 2021; Arsenault et al., 2023). The first 4 layers have a hidden state length of 64 or 128, which is close to those chosen by several authors (e.g., AlDahoul et al., 2021; Luppichini et al., 2022; Abbas et al., 2023). This choice is also based on the results of Arsenault et al. (2023), who reported that the performance of an LSTM network containing two LSTM layers with 128 neurons for simulating daily discharge is similar to that of an LSTM

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Diagram of the approach used for metamodel development and assessment.</figcaption></figure>
network with the same number of layers but with 512 neurons.

The two tested LSTM networks differ only in the nature of their inputs and outputs (Fig. 1):

- LSTM-runoff computes a series of hourly runoff discharges, taking as input a series of hourly rainfall data and hourly reference evapotranspiration (ET0) data.

- LSTM-conc computes a series of hourly pesticide concentrations. A series of hourly meteorological data, such as the difference between rainfall and reference evapotranspiration (ET0), a series including the date of application and the compound's half-life ($t_{1/2}$), a series of constant values for the compound's sorption coefficient ($K_D$), and a series of hourly discharges computed by LSTM-runoff, are used as inputs. An alternative strategy, which uses the same input variables but excludes the discharges computed by LSTM-runoff, was also tested. However, owing to its poor performance, it is not presented in this work, but further details can be found in Métayer (2024).

The inputs and outputs of the two networks are time series, each with a length equivalent to the number of hours in the simulation period. The simulation period was chosen to span from September 1st of year n−1 to February 1st of year n + 1, covering a total duration of 17 months or 12,432 h. The first 5 months served as a warm-up period, preceding the 12 months over which the metamodel performance was assessed.

The final metamodel is constructed by sequentially combining the LSTM-runoff and LSTM-conc networks to simulate runoff discharges and pesticide concentrations while accounting for the temporal variability of pesticide applications and the variability in compound properties ($K_D$ and $t_{1/2}$).

## 2.3. MHYDAS-pesticide simulations

The MHYDAS-Pesticide simulation dataset, which was used for fitting and assessing the metamodel, was generated to encompass the metamodel application domain, including variable meteorological conditions, pesticide application dates, and pesticide properties ($K_D$ and $t_{1/2}$).

### 2.3.1. MHYDAS-pesticide setup

To generate the MHYDAS-Pesticide 1.0 outputs, a 1.2 m-deep soil profile was considered. It was divided into 15 cells, with thicknesses increasing from 1 cm at the surface to 10 cm at the bottom of the profile. The values for hydraulic conductivity (Ks) and soil thickness for the uniform mixing of runoff water (zmix) correspond to those adjusted by Métayer et al. (2024) for a specific soil type in the same context: 0.07 cm·h⁻¹ and 1 cm, respectively. The other parameters were also set as in Métayer et al. (2024).

The simulation period corresponds to that targeted by the metamodel outputs, i.e., a 17-month period from September 1st of year n-1 to February 1st of year n + 1.

### 2.3.2. Variations in meteorological data and pesticide properties

The considered meteorological data included 5-minute rainfall time series and hourly reference evapotranspiration (ET₀) time series. These samples were obtained in two stages. First, daily series were estimated via the Multivariate Weather-State Approach Conditionally Skew-Normal Generator (WACS – Flecher et al., 2010), which is a stochastic weather series generator. The latter was calibrated on the basis of 19 years of observations from the Mediterranean Observatory of the Rural Environment and Water (Observatoire Méditerranéen de l'Environnement Rural et de l'Eau or OMERE) observatory (Molénat et al., 2018) used for the assessment of MHYDAS-Pesticide performance conducted by Métayer et al. (2024). In the second stage, the daily rainfall series was disaggregated to a 5-minute time step, and the ET0 series was disaggregated to an hourly time step. The disaggregation of rainfall was carried out via the *HyetosMinute* program (Kossieris et al., 2018), which

is based on the Bartlett-Lewis model and calibrated on the same rainfall data as WACS. The daily ET0 was disaggregated according to the hourly distribution of the theoretical global daily radiation for clear skies. This distribution was estimated for all days of the year for clear skies on the basis of astronomical variables, considering the geographical coordinates of the Roujan experimental catchment from OMERE. Fig. S1 shows the variability in the monthly rainfall and ET0 amounts obtained and their average values. The herbicide application date was determined for each simulation within a time window consistent with viticulture practices, i.e., between February 1st and July 31st (Simonovici, 2019; Chambre d'agriculture de l'Hérault, 2022). The application date was randomly selected from rain-free days. In total, 40 dates were determined for each ET0/rainfall pair. For each application, a unit dose of 1 kg·ha⁻¹ was assumed.

The considered pesticide properties were $t_{1/2}$ and $K_D$, with ranges determined on the basis of the minimum and maximum values reported in the Pesticide Properties DataBase (PPDB, 2022) for herbicide compounds authorized in France in 2022 for vine cultivation (Anses, 2023). The $t_{1/2}$ values used range from 0.1 days to 400 days; the $K_D$ values range from 0.05 l·kg⁻¹ to 500 l·kg⁻¹.

### 2.3.3. Experimental design

Two sets of MHYDAS-Pesticide simulations were generated.

The first dataset was used for metamodel fitting and validation. These data were used to assess the metamodel performance in simulating hourly runoff and concentration values across a wide range of dummy pesticide compound properties. These MHYDAS-Pesticide simulations resulted from a sampling scheme that combined 286 pairs of rainfall-ET0 records, 40 pesticide application dates per rainfall-ET0 pair, 7 levels of $t_{1/2}$ values, and 7 levels of $K_D$ values. The 7 levels of $t_{1/2}$ and $K_D$ values were chosen to divide the studied variation ranges of these pesticide parameters into 7 intervals of equal amplitude on a logarithmic scale. Each level was then combined with all levels of the other input variables (Faivre et al., 2013). The values of $t_{1/2}$ and $K_D$ for a given combination were determined by random selection within the corresponding level. In total, 560,560 simulations of MHYDAS-Pesticide were obtained. For the calculation of pesticide concentrations, 392,000 simulations (70 %) were used for metamodel fitting, and 168,560 (30 %) were used for validation. For the runoff calculation, 200 simulations—corresponding to the 200 different rainfall-ET0 pairs from the fitting dataset of the metamodel for calculating concentrations—were used for metamodel fitting, whereas 86 simulations were used for validation.

A second dataset was generated to assess the accuracy of the metamodel in simulating the interannual distribution of pesticide concentrations, using two actual herbicide compounds as examples: diuron (N'-(3,4-dichlorphenyl)-N,N-dimethylurea) and glyphosate (N-(phosphonomethyl)glycine). This dataset was generated by executing MHYDAS-Pesticide via the 286 rainfall-ET0 pairs employed for fitting and validating the metamodel, with one application date and the specific $t_{1/2}$ and $K_D$ value pair assigned to diuron and glyphosate. The $K_D$ and $t_{1/2}$ values for each compound correspond to those calibrated by Métayer et al. (2024). For diuron, $K_D$ was set to 34 l·kg⁻¹ and $t_{1/2}$ at 35 days, whereas for glyphosate, $K_D$ was set to 182 l·kg⁻¹ and $t_{1/2}$ at 52 days. Consequently, 286 realizations of MHYDAS-Pesticide were generated for each tested compound.

After computation, the outputs of MHYDAS-Pesticide were aggregated to an hourly time step, which was the target time step for the metamodel.

## 2.4. Fitting of the tested approximation functions

The parameters of the LSTM networks to be fitted to approximate the MHYDAS-Pesticide simulations are the set of weights and biases of the neural networks (Text S1). In total, 297,224 parameters were fitted for LSTM-runoff, and 297,736 were fitted for LSTM-conc.
## 2.4.1. Input data preprocessing

The input data were preprocessed to facilitate fitting and enhance predictive performance, particularly for LSTM-conc in hourly pesticide concentrations. LSTM-conc takes both temporal variables (rainfall, ETO and pesticide application date) and static variables ($K_D$ and $t_{1/2}$) as inputs, whereas the structure of an LSTM network is suited only to temporal input sequences. Therefore, adaptations had to be made to account for $K_D$ and $t_{1/2}$ as inputs. Thus, the value of $t_{1/2}$ was associated with the application date to calculate the theoretical evolution of a unit stock of pesticide, $y$, in the soil on the basis of first-order kinetics, according to equation (1):

$$y = \begin{cases} \log \left( e^{\frac{-\ln2}{t_1} t} \right) & \text{for } t > 0 \text{ and } y > -13 \\ -13 & \text{otherwise} \end{cases} \quad (1)$$

$t$ is the time elapsed since application (h), and $t_{1/2}$ is the compound half-life (h).

Consequently, the evolution of $y$ is integrated as input data into LSTM-conc, which represents the application date and the pesticide $t_{1/2}$. This preprocessing also facilitates the consideration of the time delay after application in the network calculation of concentrations. $K_D$ was integrated into LSTM-conc as a constant series of logarithmic values. The hourly rainfall and ETO data were combined into a single synthetic meteorological time series by subtracting one variable from the other.

All the input data were standardized by centering and scaling on the basis of the means and standard deviations calculated from the dataset used for fitting, as recommended in the literature (LeCun et al., 2012).

In summary:

- LSTM-runoff was fitted using 2 input vectors: an hourly rainfall series and an hourly ETO series;

- LSTM-conc was fitted using 4 input vectors: a synthetic meteorological data series, an hourly runoff discharge series, a series that integrates the application date and $t_{1/2}$, and a constant series of logarithmic values of $K_D$.

## 2.4.2. Fitting criteria

An LSTM network is fitted iteratively using subsets of the fitting dataset, referred to as batches (Kratzert et al., 2018). At each iteration, the parameters are updated to minimize a cost function via an optimization algorithm (Goodfellow et al., 2016). The fitting dataset can be traversed in its entirety multiple times. The number of complete passes through the dataset, known as epochs, is an additional hyperparameter to set, along with the batch size.

For fitting the LSTM-runoff and LSTM-conc networks, we chose the mean square error (MSE, Equation (2)) as the cost function, which was minimized using the ADAM optimization algorithm (Kingma and Ba, 2015). The target variable selected for fitting LSTM-conc was the decimal logarithm of hourly concentrations, with values below 1e-07 µg·1⁻¹ set to zero, and for LSTM-runoff, it was the hourly runoff discharge, with values below 0.01 mm·h⁻¹ also set to zero.

$$MSE = \frac{1}{n} \sum_{i=1}^{n} (y_{i}^{\text{meta}} - y_{i}^{\text{MHYDAS}})^{2} \quad (2)$$

$n$ is the number of discharge values/concentration log values to be approximated by the LSTM network in a batch, $y_i^{\text{meta}}$ is the i-th value of discharge values/concentration log values simulated by the LSTM network, and $y_i^{\text{MHYDAS}}$ is the i-th value of discharge values/concentration log values simulated by MHYDAS-Pesticide.

LSTM-conc was fitted by aligning runoff occurrences between LSTM-runoff and MHYDAS-Pesticide: mismatched discharges and the associated MHYDAS-Pesticide concentrations were treated as zero. This

approach minimized noise in the relationship between runoff and concentration, ensuring a more robust LSTM-conc fitting process.

LSTM-conc was fitted with a batch size of 256, the maximum supported by the computer resources (28 CPUs and 128 Go of RAM), resulting in 1,531 parameter update iterations per epoch. Due to the high computational cost of a single epoch (approximately 3 days), fitting was limited to 10 epochs. Despite this relatively low number compared to other studies (e.g., 500 epochs in Salman et al., 2018), the high iteration count per epoch was sufficient to ensure cost function convergence. LSTM-runoff was fitted with a smaller batch size of 16, in line with the smaller number of MHYDAS-Pesticide simulations available for fitting (200), and over 150 epochs with 12 iterations each. For both LSTM networks, the final parameters values were those that minimized the cost function after the completed epochs.

## 2.5. Assessment of the metamodel performance

The fitted metamodel was first assessed for its ability to simulate hourly runoff discharges and pesticide concentrations over a year across a wide range of compound properties, i.e., $t_{1/2}$ and $K_D$. Next, its accuracy in simulating the interannual distribution of concentrations for two herbicide compounds, diuron and glyphosate, was assessed. Finally, the metamodel was assessed for its ability to reduce the computational time of MHYDAS-Pesticide.

### 2.5.1. Prediction of hourly runoff discharge and concentration over a year

The performances of LSTM-runoff and LSTM-conc in simulating hourly runoff discharge and concentration series over a year were measured both for the MHYDAS-Pesticide realizations used for the fitting dataset and for those from the validation dataset (see Section 2.3). For concentrations, the overall performance assessment was conducted on samples of 5,000 randomly selected realizations within the validation dataset.

#### 2.5.1.1. Overall performance analysis.
The overall performance of the metamodel in reproducing hourly discharge or concentration log-value series was quantified using a statistical criterion of accuracy, namely, the NSE (Equation (3)), and a correctness criterion, namely, the bias. The latter was calculated as a percentage of the standard deviation of the MHYDAS-Pesticide values to be reproduced (Equation (4)).

$$NSE = 1 - \frac{\sum_{i=1}^{n} (y_{i}^{\text{MHYDAS}} - y_{i}^{\text{meta}})^{2}}{\sum_{i=1}^{n} (y_{i}^{\text{MHYDAS}} - \bar{y}^{\text{MHYDAS}})^{2}} \quad (3)$$

$$Bias = 100 \times \frac{\frac{1}{n} \sum_{i=1}^{n} (y_{i}^{\text{meta}} - y_{i}^{\text{MHYDAS}})}{\sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_{i}^{\text{MHYDAS}} - \bar{y}^{\text{MHYDAS}})^{2}}} \quad (4)$$

$y_i^{\text{meta}}$ is the i-th hourly discharge value/concentration log value simulated by one of the evaluated LSTM metamodels, $y_i^{\text{MHYDAS}}$ is the i-th hourly discharge value/concentration log value simulated by MHYDAS-Pesticide, $\bar{y}^{\text{MHYDAS}}$ is the average of the hourly discharge values/concentration log values simulated by MHYDAS-Pesticide, and $n$ is the total number of discharge values/concentration log values simulated by MHYDAS-Pesticide.

The NSE and bias criteria were chosen because reference values exist for assessing the performance of hydrological models (Moriasi et al., 2015). The bias was calculated via equation (4) instead of its usual expression as a percentage of the mean of the values since the mean of the hourly discharge values and the decimal logarithms of the concentrations simulated by MHYDAS-Pesticide can approach zero (Figs. S2 and S3).

In addition to the overall performance criteria, the performances of the metamodels were also analyzed over time by graphically comparing the simulated LSTM time series with those of MHYDAS-Pesticide.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Hourly discharge simulated by the metamodel as a function of the hourly discharge simulated by MHYDAS-Pesticide 1.0 for the fitting data (a) and validation data (b).</figcaption></figure>

2.5.1.2. Performance of runoff occurrence simulation. The ability of LSTM-runoff to simulate runoff occurrence influences the performance assessment of LSTM-conc for concentration simulation. The fitting of the latter was, in fact, guided to emphasize the relationship between the presence of runoff and the simulation of concentration occurrence (see Section 2.4.2). Thus, in cases where LSTM-runoff simulates runoff but MHYDAS-Pesticide does not, LSTM-conc is likely to simulate a high concentration, whereas MHYDAS-Pesticide produces no runoff and therefore a null concentration, and vice versa. These significant differences between the concentration values of LSTM-conc and MHYDAS-Pesticide may strongly affect the statistical performance criteria, even if the associated runoff values are low. In this sense, to understand and characterize the situations involved, the accuracy of LSTM-runoff in simulating runoff occurrence was specifically assessed. The analysis was conducted for two cases: when LSTM-runoff simulates runoff but MHYDAS Pesticide does not, and when LSTM-runoff does not simulate runoff, contrary to MHYDAS-Pesticide. In both cases, the accuracy of LSTM-runoff was assessed by analyzing the distributions of the relevant discharge values and corresponding rainfall values.

2.5.1.3. Variation in the performance of the concentration simulation due to differences in compound properties. The metamodel performance in simulating concentrations was evaluated for each combination of $K_D$ and $t_{1/2}$ levels used for generating the fitting and validation datasets (see section 2.3.3). This performance was assessed for concentration values corresponding to runoff simulated by both LSTM-runoff and MHYDAS-Pesticide.

For each level combination, the bias was expressed as a percentage of the standard deviation of all the MHYDAS-Pesticide values

corresponding to the values simulated by the assessed metamodels that were greater than $1e-07 \ \mu g \cdot l^{-1} \cdot h^{-1}$. Additionally, a linear regression between the values simulated by the LSTM metamodels and those simulated by MHYDAS-Pesticide was performed for each level combination, and the coefficient of determination ($R^2$) was calculated.

2.5.2. Prediction of the interannual distributions of diuron and glyphosate concentrations

The metamodel accuracy in simulating the interannual distribution of pesticide concentrations in surface waters was assessed for two examples, diuron and glyphosate. This temporal distribution is used for assessing the risk of surface water contamination by pesticides at the European scale (FOCUS, 2001), making it a variable of particular importance. For each compound and each of the 286 simulated years (see Section 2.3.1), the maximum hourly concentration and the annual average concentration were estimated for both the metamodel and MHYDAS-Pesticide. The 80th percentile and median values of the 286 maximum hourly concentrations and 286 annual average concentrations were then compared. The choice of the 80th percentile is based on conclusions drawn by the FOCUS working group regarding the update of the European method for assessing the risk of surface water contamination by pesticides (Adriaanse et al., 2020). This working group, in fact, proposes using a value between the 50th and 90th percentiles as an estimate of exposure.

2.5.3. Reduction in computation time

The computation time of the metamodel was compared to that of the initial model, MHYDAS-Pesticide. The computation time of the metamodel was defined as the sum of the computation times for LSTM-runoff

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison of hydrographs simulated by the metamodel (orange) and MHYDAS-Pesticide 1.0 (blue) for a year selected from the validation dataset.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. The percentage of hourly discharges simulated as zero by the metamodel but nonzero by MHYDAS-Pesticide 1.0 among all hourly discharges simulated as nonzero by MHYDAS-Pesticide 1.0 (a) and the percentage of hourly discharges simulated as zero by MHYDAS-Pesticide 1.0 but nonzero by the metamodel among all hourly discharges simulated as nonzero by the metamodel (b) as a function of the corresponding hourly rainfall and the ratio between the maximum 5-minute rainfall intensity within the hour and the hourly rainfall.</figcaption></figure>

and LSTM-conc. The computation times were estimated for 10 simulations, corresponding to 10 randomly selected realizations from the fitting dataset. The relative difference between the computation time of MHYDAS-Pesticide and that of the metamodel was calculated for each simulation according to equation (5):

$$ \text{Relativedifference (\%)} = 100 \times \frac{t^{\text{meta}} - t^{\text{MHYDAS}}}{t^{\text{MHYDAS}}} \qquad (5) $$

$t^{\text{meta}}$ and $t^{\text{MHYDAS}}$ are the computation times of the metamodel and MHYDAS-Pesticide, respectively.

# 3. Results

## 3.1. Prediction of hourly runoff discharge over a year

### 3.1.1. Overall performance

The overall performance of LSTM-runoff in simulating hourly runoff discharge was similar for both fitting and validation (Fig. 2). In both cases, the NSE and bias were 0.99 and -0.05 %, respectively. The metamodel did not simulate zero runoff values but instead produced noise between 0.03 and 0.05 mm.h⁻¹. To eliminate this noise, all values below 0.05 mm.h⁻¹ were assumed to be zero in the following analysis. The performance of LSTM-runoff in simulating nonzero hourly discharges for both fitting and validation exhibited NSE and bias values of

0.99 and -0.4 %, respectively. These high NSE values and low biases (Moriasi et al., 2015) reflect the very high accuracy of the metamodel in reproducing the hydrographs simulated by MHYDAS-Pesticide. Fig. 3 graphically illustrates this accuracy for a one-year simulation extracted from the validation set.

### 3.1.2. Performance of the runoff occurrence simulation

The hourly discharges simulated as nonzero by MHYDAS-Pesticide but zero by the metamodel represented less than 0.1 % of all the simulated discharges (both zero and nonzero) and nearly 5 % of all the nonzero discharges simulated by MHYDAS-Pesticide for both the fitting and validation data. These erroneous zero discharges simulated by the metamodel mostly corresponded to very low hourly discharges from MHYDAS-Pesticide (99 % of the erroneous values correspond to hourly discharges from MHYDAS-Pesticide that are less than 0.4 mm.h⁻¹) and low-intensity rainfall (99 % of the erroneous values correspond to rainfall events with rainfall intensities less than 2.8 mm.h⁻¹) (Fig. S4). The low number of concentration values from MHYDAS-Pesticide that were affected indicated that the modification of the dataset used for fitting the LSTM-conc data was marginal (see section 2.4.2).

The hourly discharges simulated as zero by MHYDAS-Pesticide but nonzero by the metamodel were more numerous. They represented nearly 0.4 % of all the simulated discharges (both zero and nonzero values) and approximately 30 % of all the nonzero discharges simulated by the metamodel for both the fitting and validation data. These cases again involved low simulated discharges (99 % of the erroneous nonzero simulated discharges by the metamodel were below 0.8 mm.h⁻¹), but the corresponding rainfall could be significant (99 % of the erroneous nonzero hourly simulated discharges by the metamodel corresponded to rainfall intensities below 5.2 mm.h⁻¹) (Fig. S5).

The discrepancies observed between the metamodel and MHYDAS-Pesticide in simulating runoff occurrence can be attributed in part to the difference in temporal resolution of the rainfall data used by the two models. MHYDAS-Pesticide simulated discharges using rainfall at a resolution of 5 min during rainfall events, whereas the metamodel used hourly data. Indeed, hourly rainfall data do not provide information on maximum rainfall intensity or 5 min rainfall data. Accordingly, there were mismatches between the two model predictions, especially for low hourly rainfall amounts (Fig. 4). (Fig. S6 illustrates these two situations with two examples). The lower the hourly rainfall values are, the greater the proportion of discrepancies between the metamodel and MHYDAS-Pesticide. Similarly, the proportion of discharges simulated as zero by the metamodel but nonzero by MHYDAS-Pesticide increased with a higher ratio of the maximum 5-minute rainfall intensity to the hourly rainfall total (Fig. 4a). In contrast, the proportion of discharges simulated as nonzero by the metamodel but zero by MHYDAS-Pesticide increased as this ratio decreased (Fig. 4b).

## 3.2. Prediction of hourly pesticide concentrations over a year

### 3.2.1. Overall performance

The overall performance of LSTM-conc remained stable between fitting and validation (Fig. 5). When considering all the simulated concentrations and including the cases of discrepancies in simulating runoff occurrence between LSTM-runoff and MHYDAS-Pesticide, the performance of LSTM-conc was moderate, with an NSE and bias of approximately 0.52 and 2.3 %, respectively (Fig. 5a and b). Excluding the cases of discrepancies in simulating runoff occurrence, the performance of LSTM-conc improved drastically, with an NSE and bias of 0.998 and -0.03 %, respectively (Fig. 5c and d). The low NSE value when accounting for discrepancies in runoff occurrence highlights the influence on the NSE calculation of concentration values associated with mismatches in runoff occurrence simulations. The differences between the concentrations simulated by the two models were indeed large for these cases, although the temporal dynamics simulated by LSTM-conc remained consistent (Fig. 6). Moreover, the concentrations involved
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Densities of hourly concentrations simulated by LSTM-conc as a function of hourly concentrations simulated by MHYDAS-Pesticide 1.0 for fitting data (a), (c), (e) and validation data (b), (d), (f). This includes cases that consider all discharges simulated by LSTM-runoff (a), (b), match discharge simulations (c), (d), and assign a value of $10^{-9} \mu g \cdot l^{-1} \cdot h^{-1}$ to concentrations simulated by LSTM-conc and MHYDAS-Pesticide 1.0 when the hourly rainfall is <5 mm.h$^{-1}$. Instead of displaying values, 2D density plots are presented to visualize point distributions. Points with zero concentrations for both LSTM-conc and MHYDAS are excluded from the density calculations.</figcaption></figure>

corresponded to low discharge values and thus low exported mass values. The sum of the exported masses corresponding to the concentrations simulated by LSTM-conc but not by MHYDAS-Pesticide due to discrepancies in runoff occurrence represented 3 % of the total masses simulated by LSTM-conc. The sum of the exported masses corresponding to the concentrations simulated by MHYDAS-Pesticide but not by LSTM-conc due to discrepancies in runoff occurrence represented less than 0.4 % of the total masses simulated by MHYDAS-Pesticide. Therefore, one way to use the metamodel while avoiding more than 99 % of the cases of incorrect runoff occurrence simulation was to consider only concentrations simulated for hourly rainfall greater than 5 mm.h$^{-1}$ (see Section 3.2.2).

### 3.1.2). Assigning zero values to concentrations corresponding to rainfall below this threshold did not significantly alter the annual contamination pattern simulated by MHYDAS-Pesticide (Fig. S7), and the corresponding NSE and bias values remained very high, namely, 0.98 and 0.06 %, respectively (Fig. 5e and f). Furthermore, LSTM-conc slightly underestimated concentrations below $1 \mu g \cdot l^{-1}$ (bias of approximately -7.5 %) but slightly overestimated higher values (bias of approximately 2 %).

### 3.2.2. Performance of the concentration simulation across different ranges of compound properties

The performance of LSTM-conc varied across different $K_D$ and $t_{1/2}$
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Temporal dynamics of concentrations simulated by the metamodel (orange) and MHYDAS-Pesticide 1.0 (blue) for a selected year from the validation dataset, using various $t_{1/2}$ and KD parameter pairs.</figcaption></figure>

levels, with lower accuracy observed for extreme $K_D$ and $t_{1/2}$ values. Specifically, the largest absolute bias values were associated with low $t_{1/2}$ levels (Fig. 7b), which are further illustrated in Fig. S8. Additionally, for low $t_{1/2}$ levels, Fig. 7b reveals a pronounced bias gradient across $K_D$ levels. The $R^2$ values remained above 0.7 across all $K_D$ and $t_{1/2}$ levels, reaching 0.99 (Fig. 7a).

These results suggest that the metamodel does not fully account for the effect of $K_D$ in the simulation of concentrations for those corresponding to low $t_{1/2}$ values. This issue can be related to the amount of data used to fit the metamodel. The nonzero concentrations simulated by MHYDAS-Pesticide associated with these low $t_{1/2}$ values were, in fact, lower than those associated with higher $t_{1/2}$ values (Fig. 7c). As a result, the metamodels were primarily fitted with nonzero concentration values corresponding to higher $t_{1/2}$ values.

### 3.3. Prediction of the interannual distributions of diuron and glyphosate concentrations

The metamodel generally overestimated the annual maximum and annual average concentrations of the two tested compounds, although these overestimations remained moderate (Fig. 8). The overestimations were stronger for glyphosate than for diuron, which is in line with the bias values estimated for the metamodel at the corresponding $t_{1/2}$ and $K_D$ levels (Fig. 7b). For glyphosate, the metamodel overestimated the 80th percentile and median values of the annual maximum concentrations by a factor of 2.7 and the annual average concentrations by a factor of 1.8. For diuron, the metamodel overestimated the 80th percentile and median values of the annual maximum concentrations by a factor of 1.95 and the annual average concentrations by a factor of 1.3.

Compared with MHYDAS-Pesticide, the metamodel exhibited greater variability in both the maximum and annual average concentrations, with more pronounced differences observed for the maximum concentrations. For glyphosate, the standard deviation of the annual maximum concentration was 24.8 µg·l⁻¹ for the metamodel and 7.6 µg·l⁻¹ for MHYDAS-Pesticide. Similarly, the standard deviations of the annual average concentrations were 6.9 µg·l⁻¹ and 3.2 µg·l⁻¹ for the metamodel and MHYDAS-Pesticide, respectively. For diuron, the standard deviations of the annual maximum concentration were 88.4 µg·l⁻¹ for the metamodel and 48.8 µg·l⁻¹ for MHYDAS-Pesticide, whereas the standard deviations of the annual average concentration were 12.5 µg·l⁻¹ and 10.0 µg·l⁻¹, respectively.

These results were consistent with those presented in Section 3.2. The overestimation of the maximum and annual average concentrations confirms the metamodel’s tendency to overestimate high concentrations. Furthermore, the greater variability of the simulated values observed for the metamodel can be explained by its tendency to underestimate lower concentrations (below 1 µg·l⁻¹) to compensate for the overestimation of higher concentrations.

### 3.4. Reduction in computation time

The developed metamodel significantly reduces the computation time of MHYDAS-Pesticide. The execution times for LSTM-runoff and LSTM-conc were similar, nearing 2.7 s. A complete metamodel, composed of LSTM-runoff and LSTM-conc, reduced the simulation time by an average of 90%.

## 4. Discussion

This study highlights the potential of an LSTM model for metamodeling a Richards-based pesticide runoff model and offers several methodological insights. The LSTM metamodel is effective in simulating the hourly variation in water contamination for a wide range of pesticide properties and dates of pesticide application. Consequently, it offers broader application potential than the metamodels proposed by Stenemo et al. (2007) and Reichenberger et al. (2021) for analyzing water contamination by pesticides from agricultural fields. The latter are limited to simulating static variables and do not account for the variability in pesticide application over time, unlike the proposed LSTM-based metamodel. Our metamodel also accounts for variations in the pesticide application rate and in the number of applications throughout the year. This capability arises from the first-order kinetics and linear sorption equilibrium assumptions in the initial MHYDAS-Pesticide model, which establish a linear relationship between the concentration and pesticide application rate. As a result, simulated concentrations can be adjusted using a multiplicative factor to accommodate different application rates. However, the metamodeling approach was developed only for a specific soil type, which limits its operational applicability. Future research should therefore focus on incorporating the effects of soil type on runoff generation and pesticide fate. For example, influential physical parameters, such as topsoil hydraulic conductivity, could be integrated as additional inputs to the metamodel.
<figure><img src="image_8.png" /><figcaption>Fig. 7. R² a) and bias b) values calculated with the concentration log values simulated by LSTM-conc and MHYDAS-Pesticide 1.0 and c) the number of concentration values for each combination of KD level and t1/2 studied.</figcaption></figure>

As expected, the developed metamodeling approach significantly reduces the computation time of MHYDAS-Pesticide. Indeed, neural networks are known for their short execution times despite requiring considerable computational effort during the fitting phase (Adikari et al., 2021). However, the reduction achieved is slightly smaller than that typically observed in metamodeling, where reductions of up to 98% have been reported (Razavi et al., 2012). This limitation can be attributed to the multilayer architecture of the tested LSTM networks. The metamodel comprises two LSTM networks, each with 5 layers and approximately 300,000 parameters. Optimizing the number of layers and the length of the hidden state, as suggested by Jin et al. (2022), could increase computational efficiency. Additionally, the metamodel and MHYDAS-Pesticide are not directly comparable due to differences in the sequences they simulate. The metamodel generates hourly discharges and concentrations at a fixed time step, whereas MHYDAS-Pesticide uses a variable time step of less than 5 min during runoff events. Fitting the LSTM networks on sequences equivalent to 17 months

at a 5-minute resolution (149,184 time steps) would be computationally prohibitive and would also constrain the operational implementation of the metamodel. Subhourly meteorological data, which are required as model inputs, are indeed often more difficult to obtain.

The overall performance in predicting runoff discharge and pesticide concentration is excellent, based on the performance thresholds set by Moriasi et al. (2015). The achieved accuracy in simulating pesticide concentrations makes the metamodel well suited for operational use. Notably, the results indicate its applicability for assessing both chronic and acute risks of surface water contamination. Nonetheless, there is room for improvement in metamodel performance in specific areas. First, the metamodel's tendency to overestimate high concentrations and underestimate low concentrations could be reduced by excluding low concentration values from the fitting phase. For example, concentrations below the detection thresholds of current analytical methods (approximately 10⁻³ µg·l⁻¹; Eaufrance, 2023) could be omitted, as the accuracy of the model cannot be verified for these unmeasurable concentrations. Second, the metamodel accuracy varies on the basis of the values of K_D and t₁/₂. It notably deteriorates for extreme parameter values, particularly for low t₁/₂ values, where only a few nonzero concentrations are simulated. This observation is common in the fitting of a statistical model (Geman et al., 1992). The metamodel accuracy could then be improved by fitting several local sub-metamodels with narrower ranges of compound properties. For example, an LSTM-conc network fitted with narrower ranges of K_D and t₁/₂, which include the values selected for glyphosate (K_D ∈ [134;500 l·kg⁻¹] and t₁/₂ ∈ [37;122 days]), significantly improves accuracy (Fig. S9). The interannual distributions of glyphosate concentrations obtained are similar to those simulated by MHYDAS-Pesticide, without overestimation (Fig. S10). Other strategies to increase the performance of LSTM-conc can also be proposed. In particular, the results indicate that the effect of variations in K_D on concentration simulation is not fully accounted for by LSTM-conc, especially for low values of t₁/₂. Thus, our choice of integrating this static parameter into the LSTM network as a constant time series does not seem to be fully satisfactory. The use of static variables as inputs to LSTM models has been discussed in the literature. Alternative approaches that could be explored include (i) preprocessing these variables (Pyo et al., 2023) in a manner similar to how we handled t₁/₂ or (ii) modifying the LSTM cells as proposed by Kratzert et al. (2019b). Finally, the metamodel also encounters limitations in accurately simulating the occurrence of hourly discharges. The results suggest that this issue is partly due to the difference in temporal resolution between the original MHYDAS-Pesticide model (up to 5 min) and the metamodel (1 h). These findings underscore the challenge of accurately predicting surface runoff at larger time steps. A potential way to improve the metamodel performance could involve incorporating an input variable related to the maximum subhourly rainfall intensity for each hour of rainfall. For example, the ratio of the maximum subhourly intensity to the total hourly rainfall could be included.

This work confirms the strength of an LSTM model in simulating subdaily variations in runoff discharge and water quality variables over a year, even with a limited number of input variables. This finding supports previous studies that demonstrated the ability of LSTM to model the rainfall-runoff relationship at a daily time step over a year (Kratzert et al., 2018, 2019a; b; AlDahoul et al., 2021; Frame et al., 2021; Lees et al., 2021; Yin et al., 2021; Abbas et al., 2023; Gholizadeh et al., 2023; Liu et al., 2024) and at hourly time steps over shorter periods (Kao et al., 2020; Xiang et al., 2020). Some authors have attempted to explain LSTM performance by interpreting its calculations in physical terms. For example, Kratzert et al. (2018) identified a correlation between variations in influential state variables (e.g., soil water content) and LSTM internal values. As proposed by Pyo et al. (2023), a similar interpretation may apply to water quality variables. When simulating pesticide concentrations, the internal values of the LSTM models could be correlated with the pesticide content in the soil. However, interpreting LSTM models generally remains challenging. Their black-box
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the distributions of annual maximum hourly concentration values simulated by the metamodel and MHYDAS-Pesticide 1.0 for glyphosate (a) and diuron (c) and comparisons of the distributions of annual average concentration values for glyphosate (b) and diuron (d).</figcaption></figure>

nature is frequently criticized, particularly when fitting datasets are derived from observed data (Tripathy and Mishra, 2024). Consequently, LSTM models may lack physical consistency and may capture incorrect relationships between hydrochemical variables. Metamodeling helps address this issue by enhancing interpretability, as all the physical mechanisms involved in generating the target variables, which the metamodel must learn, are known. Furthermore, in our approach, we incorporated a representation of a physical mechanism—namely, first-order kinetics for compound degradation—into the LSTM model to predict pesticide runoff concentrations. These approaches are part of ongoing efforts to overcome the limitations of black-box models by hybridizing them with physically based models (Zhang et al., 2022; Tripathy and Mishra, 2024).

## 5. Conclusion

This study offers insights for metamodeling a physically based pesticide fate model to simulate the hourly variability in water and pesticide runoff over a year. The results show the potential of an LSTM based model for this purpose. The proposed LSTM-based metamodel accurately simulates water and pesticide runoff by accounting for both static and temporal key factors. It enables predictions under highly variable rainfall conditions, accounting for a wide range of pesticide properties and the temporal variability of pesticide applications. This is achieved while reducing computation time by up to 90 %. Thus this LSTM-based metamodeling approach shows great potential for large-scale and computationally demanding applications, such as regional water quality assessments, that involve simulations of the area-wide distribution of influential input factors, or landscape-scale analyses, by integrating it into a fully spatially distributed model. However, achieving this goal requires further development of the approach to enhance the metamodel's genericity. In particular, further efforts are required to account for the effects of soil type variations on runoff

generation. This could be addressed by integrating additional inputs into the metamodel, such as influential soil physical parameters. In fact, this work raises the question of how to effectively integrate static physical inputs into LSTM architectures for simulating water quantity and quality dynamics. One potential solution is to encode such parameters into temporal signals tailored to their physical significance. Overall, the proposed metamodeling approach could be adapted for other water quality models, with modifications according to the specific physical processes represented by the initial physically based model.

This work also highlights the challenge of developing accurate runoff prediction models at temporal resolutions aligned with data availability. The strong overall performance of the hourly metamodel demonstrates its ability to capture the subtleties of a subhourly Richards-based model. However, its limitations in simulating runoff occurrence underlines the critical importance of incorporating subhourly rainfall variability to achieve accurate simulations across the full spectrum of runoff variability.

## CRediT authorship contribution statement

**Guillaume Métayer:** Writing – original draft, Validation, Software, Investigation, Formal analysis, Conceptualization. **Cécile Dagès:** Writing – review & editing, Validation, Supervision, Funding acquisition, Conceptualization. **Marc Voltz:** Writing – review & editing, Supervision, Funding acquisition, Conceptualization. **Jean-Stéphane Bailly:** Writing – review & editing, Supervision, Funding acquisition, Conceptualization.

## Funding

The senior author's PhD grant was funded by the French National Research Institute for Agriculture, Food and Environment (INRAE) and the Occitanie Region. This study was conducted as part of the GREET
project, supported by the Occitanie Region and the European Regional Development Fund (FEDER) (grant LR 0014594/2017-007745-01).

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

We thank the OMERE observatory (http://wwwobs-omere.org) for providing the data that supported this work. We also thank David Crevoisier, a scientist at the French National Research Institute for Agriculture, Food and Environment (INRAE), for his help in optimizing model calculations.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jhydrol.2025.133800.

## Data availability

Data will be made available on request.

## References

Abbas, A., Park, M., Baek, S.S., Cho, K.H., 2023. Deep learning-based algorithms or long-term prediction o chlorophyll-a in catchment streams. J. Hydrol. 626, 130240. https://doi.org/10.1016/j.jhydrol.2023.130240.

Adikari, K.E., Shrestha, S., Ratnayake, D.T., Budhathoki, A., Mohanasundaram, S., Dailey, M.N., 2021. Evaluation of artificial intelligence models for flood and drought forecasting in arid and tropical regions. Environ. Model. Softw. 144, 105136. https://doi.org/10.1016/j.envsoft.2021.105136.

Adriaanse, P., Boivin, A., Klein, M., Jarvis, N., Stemmer, M., Fait, G., Egsmose, M., 2020. Scientific report of EFSA on the 'repair action' of the FOCUS surface water scenarios. Ej EFSA J. 18 (6), 6119. https://doi.org/10.2903/j.efsa.2020.6119.

AlDahoul, N., Essam, Y., Kumar, P., Ahmed, A.N., Sheri, M., Seelnasr, A., Elshafe, A., 2021. Suspended sediment load prediction using long short-term memory neural network. Sci. Rep. 11 (1), 7826. https://doi.org/10.1038/s41598-021-87415-4.

Alvi, M., Batstone, D., Mbamba, C.K., Keymer, P., French, T., Ward, A., Dwyer, J., Cardell-Oliver, R., 2023. Deep learning in wastewater treatment: a critical review. Water Res. 245, 120518. https://doi.org/10.1016/j.watres.2023.120518.

Agence nationale de sécurité sanitaire, de l'alimentation, de l'environnement et du travail (Anses), Données ouvertes du catalogue E-Phy des produits phytopharmaceutiques, matières fertilisantes et supports de culture, adjuvants, produits mixtes et mélanges. 2023.

Arsenault, R., Martel, J.-L., Brunet, F., Brissette, F., Mai, J., 2023. Continuous streamflow, prediction in ungauged basins: long short-term memory neural networks clearly outperform traditional hydrological models. Hydrol. Earth Syst. Sci. 27, 139-157. https://doi.org/10.5194/j.hess-27-139-2023.

Bengio, Y. 2012. Practical recommendations for gradient-based training of deep architectures, in: Neural networks: Tricks of the trade, edited by: Montavon, G., Orr, G. B., and Müller, K.-R., Springer, 437-478, doi: 10.1007/978-3-642-35289-8_26.

Boivin, A., Poulsen, V., 2017. Environmental risk assessment of pesticides: state of the art and prospective improvement from science. Environ. Sci. Pollut. Res. 24, 6889-6894. https://doi.org/10.1007/s11356-016-8289-2.

Bronstert, A., Niehoff, D., Schiffler, G.R., 2023. Modeling infiltration and infiltration excess: the importance of fast and local processes. Hydrol. Process. 37, e14875. https://doi.org/10.1002/hyp.14875.

Calvet, R., Barriuso, E., Bedos, C., Benoit, P., Charnay, M., Coquet, Y. 2005. Les pesticides dans les sols. Conséquences agronomiques et environnementales. France Agricole, Paris.

Carsel, R.F., Mulkey, L.A., Lorber, M.N., Baskin, L.B., 1985. The Pesticide root Zone Model (PRZM): a procedure for evaluating pesticide leaching threats to groundwater. Ecol. Model. 30, 49-69. https://doi.org/10.1016/0304-3800(85)90036-5.

Castelletti, A., Galelli, S., Ratto, M., Soncini-Sessa, R., Young, P.C., 2012. A general framework for dynamic emulation modeling in environmental problems. Environ. Model. Softw. 34, 5-18. https://doi.org/10.1016/j.envsoft.2012.01.002.

Chambre d'agriculture de l'Hérault, 2022. Note d'information technique "glyphosate" Evolution des conditions d'utilisation en viticulture. Note technique, 1-2.

Cohen, M. 2023. Metamodel and bayesian approaches for dynamic systems. General Mathematics [math.GM]. Institut Polytechnique de Paris, English. NNT: 2023IPPAS003. Tel.04337924.

Crevoisier, D., Chanzy, A., Voltz, M., 2009. Evaluation of the Ross fast solution of Richards' equation in unfavourable conditions for standard finite element methods. Adv. Water Resour. 32, 936-947. https://doi.org/10.1016/j.advwatres.2009.03.008.

Crevoisier, D., Dagès, C., Fabre, J.C., Roman-Villafane, S., Voltz, M. 2021. MHYDAS-Pesticides-1.0 : Modélisation Hydrologique Distribuée des AgroSystèmes - Eau et Pesticides. Rapport dans le cadre du projet collaboratif PollDiff Capture 2018-2021 (AAP Readynov FEDER-FSE-Occitanie), INRAE Montpellier, UMR LISAH. https://hal.inrae.fr/hal-04090048v1.

Eaufrance, 2023. Naiades, données sur la qualité des eaux de surface. https://naiades.eaufrance.fr/acces-donnees#/physicochimie.

Faivre, R., Iooss, B., Mahevas, S., Makowski, D., Monod, H. 2013. Analyse de sensibilité et exploration de modèles : Application aux sciences de la nature et de l'environnement. Editions Quae, Versailles.

Flecher, C., Naveau, P., Allard, D., Brisson, N.N., 2010. A stochastic daily weather generator for skewed data. Water Resour. Res. 2010 (46), W07519. https://doi.org/10.1029/2009WR008098.

FOCUS, 2001. Focus surface water scenarios in the EU evaluation process under 91/414/EEC. Report prepared by the FOCUS Working Group on Surface Water Scenarios, EC Document Reference SANCO/4802/2001-rev.2 final.

Frame, J.M., Kratzert, F., Raney, A., Rahman, M., Salas, F.R., Nearing, G.S., 2021. Post-processing the national water model with long short-term memory networks for streamflow predictions and model diagnostics. J. Am. Water Resour. Assoc. 57 (6), 885-905. https://doi.org/10.1111/1752-1688.12964.

Gauch, M., Mai, J., Lin, J., 2021. The proper care and feeding of CAMELS: how limited training data affects streamflow prediction. Environ. Model Softw. 135, 104926. https://doi.org/10.1016/j.envsoft.2020.104926.

Geman, S., Bienenstock, E., Doursat, R., 1992. Neural networks and the bias/variance dilemma. Neural Comput. 4 (1), 1-58.

Gholizadeh, H., Zhang, Y., Frame, J., Gu, X., Green, C.T., 2023. Long short-term memory models to quantify long-term evolution of streamflow discharge and groundwater depth in Alabama. Sci. Total Environ. 901, 165884. https://doi.org/10.1016/j.scitotenv.2023.165884.

Goodfellow, I., Bengio, Y., Courville, A. 2016. Deep Learning. MIT Press, available at: http://www.deeplearningbook.org.

Hochreiter, S., Schmidhuber, J., 1997. Long short-term memory. Neural Comput. 9, 1735-1780. https://doi.org/10.1162/neco.1997.9.1735.

Jin, M., Liao, Q., Patil, S., Abdulaheem, A., Al-Shehri, D., Glatz, G., 2022. Hyperparameter tuning of artificial neural networks for well production estimation considering the uncertainty in initialized parameters. ACS Omega 7 (28), 24145-24156. https://doi.org/10.1021/acsomega.2c00498.

Kao, I., Zhou, Y., Chang, L., Chang, F., 2020. Exploring a long short-term memory based encoder-decoder framework for multi-step-ahead flood forecasting. J. Hydrol. 583, 124631. https://doi.org/10.1016/j.jhydrol.2020.124631. ISSN 0022-1694.

Kingma, D.P., Ba, J. 2015. Adam: a method for stochastic optimization. 3rd International Conference for Learning Representations, San DiegoCoRR, abs/1412.6980., doi: 10.48550/arXiv.1412.6980.

Kleijnen, J.P.C., 2009. Kriging metamodeling in simulation: a review. Eur. J. Oper. Res. 192 (3), 707-716. https://doi.org/10.1016/j.ejor.2007.10.013.

Kossieris, P., Makropoulos, C., Onof, C., Koutsoyiannis, D., 2018. A rainfall disaggregation scheme for sub-hourly time scales: coupling a Bartlett-Lewis based model with adjusting procedures. J. Hydrol. 556, 980-992. https://doi.org/10.1016/j.jhydrol.2016.07.015.

Kratzert, F., Klotz, D., Brenner, C., Schulz, K., Herrnegger, M., 2018. Rainall-runoff modeling using long short-term memory (LSTM) networks. Hydrol. Earth Syst. Sci. 22, 6005-6022. https://doi.org/10.5194/hess-22-6005-2018.

Kratzert, F., Klotz, D., Herrnegger, M., Sampson, A.K., Hochreiter, S., Nearing, G.S., 2019a. Toward improved predictions in ungauged basins: exploiting the power o machine learning. Water Resour. Res. 55, 11344-11354. https://doi.org/10.1029/2019WR026065.

Kratzert, F., Klotz, D., Shalev, G., Klambauer, G., Hochreiter, S., Nearing, G., 2019b. Towards learning universal, regional, and local hydrological behaviors via machine learning applied to large-sample datasets. Hydrol. Earth Syst. Sci. 23, 5089-5110. https://doi.org/10.5194/hess-23-5089-2019.

Larsbo, M., Roulier, S., Stenemo, F., Kasteel, R., Jarvis, N., 2005. An improved dual-permeability model of water flow and solute transport in the vadose zone. Vadose Zone J. 4, 398-406. https://doi.org/10.2136/vzj2004.0137.

LeCun, Y.A., Bottou, L., Orr, G.B., Müller, K.R., 2012. Efficient Backprop. Springer, Berlin, Heidelberg, Germany.

Lees, T., Buechel, M., Anderson, B., Slater, L., Reece, S., Coxon, G., Dadson, S.J., 2021. Benchmarking data-driven rainfall-runoff models in Great Britain: a comparison of long short-term memory (LSTM)-based models with four lumped conceptual models. Hydrol. Earth Syst. Sci. 25, 5517-5534. https://doi.org/10.5194/hess-25-5517-2021.

Leonard, R.A., Knisel, W.G., Still, D.A., 1987. Gleams: groundwater loading effects of agricultural management systems. Tran. Am. Soc. Agric. Biol. Eng. 30, 1403-1418. https://doi.org/10.13031/2013.30578.

Lim, T., Wang, K., 2022. Comparison of machine learning algorithms for emulation of a gridded hydrologic model given spatially explicit inputs. Comput. Geosci. 159, 105025. https://doi.org/10.1016/j.cageo.2021.105025. ISSN 0098-3004.

Liu, J., Bian, Y., Lawson, K., Shen, C., 2024. Probing the limit of hydrologic predictability with the transformer network. J. Hydrol. 637, 131389. https://doi.org/10.1016/j.jhydrol.2024.131389.

Louchart, X., Voltz, M., Andrieux, P., Moussa, R., 2001. Herbicide transport to surface waters at field and watershed scales in a Mediterranean vineyard area. J. Environ. Qual. 30, 982-991. https://doi.org/10.2134/jeq2001.303982x.
Luppichini, M., Barsanti, M., Giannecchini, R., Bini, M., 2022. Deep learning models to predict flood events in fast-flowing watersheds. Sci. Total Environ. 813, 151885. https://doi.org/10.1016/j.scitotenv.2021.151885. ISSN 0048-9697.

Marrel, A., Iooss, B. 2023. Probabilistic surrogate modeling by Gaussian process: a new estimation algorithm for more reliable prediction. ea-04322818v2.

Métayer, G. 2024. Simulation par approche mixte mécaniste-statistique du devenir des pesticides dans l'environnement pour le diagnostic de la contamination des sols et des eaux. PhD Thesis, Institut Agro Montpellier. 110-155.

Métayer, G., Dagès, C., Crevoisier, D., Bailly, J.-S., Voltz, M., 2024. Performance of a pesticide fate model for predicting multi-year surface runoff contamination in a Mediterranean vineyard. Sci. Total Environ. 906, 167357. https://doi.org/10.1016/j.scitotenv.2023.167357.

Molénat, J., Raclot, D., Zitouna, R., Andrieux, P., Coulouma, G., Feurer, D., Grünberger, O., Lamachère, J.-M., Bailly, J.-S., Belotti, J.-L., Ben Azzez, K., Ben Mechlia, N., Ben Younès Louati, M., Biarnès, A., Blanca, Y., Carriere, D., Chaabane, H., Dagès, C., Debabria, A., Dubreuil, A., Fabre, J.C., Fages, D., Floure, C., Garnier, F., Geniez, C., Gomez, C., Hamdi, R., Huttel, O., Jacob, F., Jenhaoui, Z., Lagacherie, P., Le Bissonnais, Y., Louati, R., Louchart, X., Mekki, I., Moussa, R., Negro, S., Pepin, Y., Prevot, L., Samouélian, A., Seidel, J.L., Trotoux, G., Troiano, S., Vinatier, F., Zante, P., Zrelli, J., Albergel, J., Voltz, M., 2018. OMERE: a long-term observatory of soil and water resources, in interaction with agricultural and land management in Mediterranean hilly catchments. Vadose Zone J. 17, 1–18. https://doi.org/10.2136/vzj2018.04.0086.

Moriasi, D.N., Gitau, M.W., Pai, N., Daggupati, P., 2015. Hydrologic and water quality models: performance measures and evaluation criteria. Trans. Am. Soc. Agric. Biol. Eng. 58, 1763–1785. https://doi.org/10.13031/trans.58.10715.

Moussa, R., 1996. Analytical Hayami solution for the diffusive wave flood routing problem with lateral inflow. Hydrol. Process. 10, 1209–1227. https://doi.org/10.1002/(SICI)1099-1085(199609)10:9.

Nearing, G.S., Klotz, D., Sampson, A.K., Kratzert, F., Gauch, M., Frame, J.M., Shalev, G., Nevo, S., 2022. Data assimilation and autoregression for using near-real-time streamflow observations in long short-term memory networks. Hydrol. Earth Syst. Sci. 26 (21), 5493–5513. https://doi.org/10.5194/hess-26-5493-2022.

PPDB, 2022. Pesticide properties database, University of Hertford-shire, https://sitem.herts.ac.uk/aeru/ppdb/en/.

Pyo, J., Pachepsky, Y., Kim, S., Abbas, A., Kim, M., Sung Kwon, Y., Ligaray, M., Cho, K. H., 2023. Long short-term memory models of water quality in inland water environments. Water Res. X 21, 100207. https://doi.org/10.1016/j.wroa.2023.100207. ISSN 2589-9147.

Razavi, S., Tolson, B.A., Burn, D.H., 2012. Review of surrogate modeling in water resources. Water Resour. Res. 48, W07401. https://doi.org/10.1029/2011WR011527.

Reichenberger, S., Bach, M., Skitschak, A., Frede, H.G., 2007. Mitigation strategies to reduce pesticide inputs into ground- and surface water and their effectiveness; a review. Sci. Total Environ. 384, 1–35. https://doi.org/10.1016/j.scitotenv.2007.04.046.

Reichenberger, S., Multsch, S., Pohlert, T., Jarvis, N., Gönczi, M., Kreuger, J. 2021. Developing a MACRO meta-model for Swedish drinking water abstraction zones, SETAC Europe 31 st annual meeting, Helsinki; poster pres.

Richards, A.L., 1931. Capillary conduction of liquids in porous mediums. Physics 1, 318–333. https://doi.org/10.1063/1.1745010.

Ross, P.J., 2003. Modeling soil water and solute transport – Fast, simplified numerical solutions. Agron. J. 95, 1352–1361. https://doi.org/10.2134/agronj2003.1352.

Salman, A.G., Heryadi, Y., Abdurahman, E., Suparta, W., 2018. Single layer and multi-layer long short-term memory (LSTM) model with intermediate variables for weather forecasting. Procedia Comput. Sci. 135, 89–98. https://doi.org/10.1016/j.procs.2018.08.153.

Saint-Venant, B. 1871. Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leurs lits, Comptes rendus de l'Académie des Sciences de Paris, 73, 147-154, 237-240.

Schaap, M.G., van Genuchten, M.T., 2006. A Modified Mualem-van Genuchten formulation for improved description of the hydraulic conductivity near saturation. Vadose Zone J. 5, 27–34. https://doi.org/10.2136/vzj2005.0005.

Shipitalo, M.J., Owens, L.B., 2003. Atrazine, deethylatrazine, and deisopropylatrazine in surface runoff from conservation tilled watersheds. Environ. Sci. Tech. 37, 944–950. https://doi.org/10.1021/es020870b.

Siimes, K., Kämäri, J., 2003. A review of available pesticide leaching models: selection of models for simulation of herbicide fate in finnish sugar beet cultivation. Boreal Environ. Res. 8, 31–51.

Simonovici, M., 2019. Enquête pratiques phytosanitaires en viticulture en 2016 Nombre de traitements et indicateurs de fréquence de traitement. Agreste Les Dossiers, 1-50.

Sreejith, K.S., Kumar, G.P., Dwarakish, G.S., 2024. A critical review of the soil conservation services - curve number method in hydrological modelling. Wetlands 44. https://doi.org/10.1007/s13157-024-01873-w.

Stenemo, F., Lindahl, A.M.L., Gärdenäs, A., Jarvis, N., 2007. Meta-modeling of the pesticide fate model MACRO for groundwater exposure assessments using artificial neural networks. J. Contam. Hydrol. 93 (1-4), 270–283. https://doi.org/10.1016/j.jconhyd.2007.03.003.

Tiktak, A., Boesten, J.J.T.I., van der Linden, A.M.A., Vancoloster, M., 2006. Mapping ground water vulnerability to pesticide leaching with a process-based metamodel of EuroPEARL. J. Environ. Qual. 35 (4), 1213–1226. https://doi.org/10.2134/jeq2005.0377.

Tripathy, K.P., Mishra, A.K., 2024. Deep learning in hydrology and water resources disciplines: concepts, methods, applications, and research directions. J. Hydrol. 628, 130458. https://doi.org/10.1016/j.jhydrol.2023.130458. ISSN 0022-1694.

Walker, A., 1974. A simulation model for prediction of herbicide persistence. J. Environ. Qual. 3, 396-401. https://doi.org/10.2134/jeq1974.00472425000300040021x.

Wauchope, R.D., Graney, R.L., Cryer, S., Eadsforth, C., Klein, A.W., Racke, K.D., 1995. Pesticide runoff: methods and interpretation of field studies. Pure Appl. Chem. 67, 2089–2108. https://doi.org/10.1351/pac199567122089.

Williams, J.R., Wang, E., Meinardus, A., Harman, W.L., Siemens, M., Atwood, J.D., 2006. APEX Users Guide. Version 2110. Temple, TX: Texas A&M University, Texas Agricultural Extension Service, Texas Agricultural Experiment Station, Blacklands Research Center.

Xiang, Z., Yan, J., Demir, I., 2020. A rainfall-runoff model with LSTM-based sequence-to-sequence learning. Water Resour. Res. 56, e2019WR025326. https://doi.org/10.1029/2019WR025326.

Yin, H., Zhang, X., Wang, F., Zhang, Y., Xia, R., Jin, J., 2021. Rainfall-runoff modeling using LSTM-based multi-state-vector sequence-to-sequence model. J. Hydrol. 598. https://doi.org/10.1016/j.jhydrol.2021.126378.

Young, D.F., 2019. U.S. Environmental protection agency model for estimating pesticides in surface water. ACS Sympos. Ser. 1308. https://doi.org/10.1021/bk-2019-1308. ch016.

Zhang, C., Xue, H., Dong, G., Jing, H., He, S., 2022. Runoff estimation based on hybrid-physics-data model. ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci., V-3-2022, 347–352. https://doi.org/10.5194/isprs-annals-V-3-2022-347-2022.