Journal of Pesticide Science
Online ISSN : 1349-0923
Print ISSN : 1348-589X
ISSN-L : 0385-1559
Society Awards 2022
Use of mathematical modeling and its inverse analysis for precise assessment of pesticide dissipation in a paddy environment
Kei Kondo
著者情報
ジャーナル オープンアクセス HTML

2022 年 47 巻 3 号 p. 146-153

詳細
Abstract

The extrapolability of the lysimeter test as a dissipation simulator in an actual paddy field was evaluated using mathematical models and their inverse analyses for predicting pesticide fate and transport processes in paddy test systems. As a source of experimental data, a four-year comparative experiment in lysimeters and paddy fields was conducted using various paddy pesticides. First, the dissipations for various active ingredients in granule pesticides under submerged applications were statistically compared using simple kinetic modeling. Second, the dissipation pathways, unobserved experimental components, and effect of the experimental setting were evaluated using a higher tier mathematical model with a novel inverse analysis protocol. Finally, owing to experimental constraints, the unobtainable parameters were extracted from the laboratory container test before being transferred to compare the outdoor experimental data under different formulation types.

Introduction

More than 50% of agricultural land comprises paddy fields,1) thus rice cultivation is one of the most important agricultural industries in Japan. Paddy pesticides are indispensable agricultural materials that prevent pest and disease crop damage during the cultivation season because of the Asian monsoon conditions, which bring higher precipitation and humidity than in the arid regions.2) On the other hand, it is widely known that pesticides used in paddy fields are remarkably prone to runoff to the open environment because of flooding conditions in the fields.3) Therefore, the surface runoff containing pesticides, mainly due to irrigation and meteorological events, is of great concern because of its potential impact on the aquatic environment in Japan.

To account for the potential environmental risk of pesticides, the Ministry of Environment has developed a simple screening model comprising a conceptual 100 km2 watershed, including main streams, tributaries, paddy fields, and non-paddy agricultural fields. In the model, the exposure is characterized as the predicted environmental concentration (PEC) and compared with the water quality standard calculated from the acceptable daily intake (ADI) or the lowest acute effect concentration (AEC). For the PEC determination of the surface runoff from paddy fields, the tier system is adopted and comprises three stages: numerical calculation (tier 1), lysimeter test (tier 2), and field test (tier 3). Regardless of the abundant studies dealing with the dissipation of paddy pesticides in lysimeters and actual paddy fields over the last three decades,46) few studies have compared the dissipation of paddy pesticides under coordinated test conditions. Hence, scientific knowledge regarding the extrapolability of the lysimeter test as a dissipation simulator for actual paddy fields is insufficient.

Modeling is a part of the scientific method and should not be separated from the experiment.7) Either statistical or mathematical models can facilitate the interpretation of the mechanisms hidden in the experimental data through abstraction, interpolation, and extrapolation. For the analysis of pesticide dissipation, a mathematical model is preferable because of its explicitly defined structure in mathematical language. Indeed, various mathematical models predicting the fate and transport of paddy pesticides in Japan have been developed and their applicability verified.813) However, the major drawback of these models is that model prediction easily loses accuracy without appropriate calibration of numerous model parameters to site-specific conditions. Unlike research purpose modeling, where input parameters are usually well prepared, regulatory modeling must provide a generic assessment with limited data and limited relevant information.14)

The present study proposes a novel modeling approach for the precise assessment of pesticide dissipation in paddy test systems using a mathematical model and inverse analysis. This approach first inversely calibrates the model parameters from minimal experimental data, and then extracts the quantitative information behind the data from the behavior of the calibrated model. A schematic diagram of this approach is shown in Fig. 1, and the three perspectives (comparison, extraction, and relation) regarding experimental data are explained in the following sections.

Fig. 1. Schematic diagram of mathematical modeling and its inverse analysis under the paddy test system.

1. Comparison of experimental data15)

1.1. Comparative experiment using lysimeters and paddy fields

To investigate the extrapolability of the lysimeter as a simulator of pesticide dissipation in actual paddy fields, a four-year (2012–2015) comparative experiment with lysimeters and paddy fields was conducted using various paddy pesticides. The lysimeter experiments were conducted under outdoor conditions using the experimental facilities of the Institute of Environmental Toxicology (IET). Each lysimeter was made of concrete with a surface area of 1 m2 (1 m×1 m) and a depth of 1 m and was covered with a roof panel, penetrating UV radiation. Two types of soil, alluvial and volcanic ash, were packed in the top 50 cm soil layer and labeled respective lysimeters as LA and LV. The flooding conditions were maintained by controlling the daily percolation rate at 1.5 cm/day. The daily evapotranspiration (ET) was calculated from the water requirement and percolation amount. Irrigation water was pumped from groundwater to maintain a paddy water depth of 5 cm.

Two well-managed experimental paddy fields containing alluvial and volcanic ash soils (hereafter labeled as FA and FV, respectively) were selected for the paddy field experiments. Both fields had 800 m2 surface areas and were located within 30 km of the IET lysimeter facilities. Daily ET amounts were estimated by the FAO Penman-Monteith method16) and the daily percolation amounts were estimated from the water balance in the experimental paddy plot.17) A summary of the key water balance components highlighted in this paper is shown in Table 1.

Table 1. Summary of key components for water balance in experimental paddy fields
2012 2013 2014 2015
Alluvial soil plot
Paddy water depth 4.4±1.19 6.4±1.60 6.3±2.39 4.0±1.20
Average daily percolation 0.43 0.14 0.09 0.58
Cumulative surface runoff 6.5 0.0 0.0 0.0
Volcanic ash soil plot
Paddy water depth 5.4±1.67 5.5±0.88 5.0±0.80 5.9±1.12
Average daily percolation 0.54 0.67 0.41 0.54
Cumulative surface runoff 6.0 0.0 0.0 0.0

Unit: cm

Over four years, a total of 14 pesticide formulation products were applied, and 20 active ingredients and 4 metabolites were analyzed. Every year, the application date and sampling periods of paddy water (PW) were coordinated, with samples being taken before application and 0.125, 1, 2, 3, 5, 7, 8, 10, 14, and 21 days after treatment (DAT). Concentrations in the water samples were determined using validated analytical methods.15)

1.2. Simple kinetic modeling and grouping analysis

Among the abovementioned experimental data, the dissipation data of 13 active ingredients in 6 granule formulation products under submerged applications were selected to evaluate the experimental performance of the lysimeter in an actual paddy field (the experimental details are found in the previous paper15)). To model the dissipation of the submerged application of the granule, the single first-order (SFO) model was coupled with another kinetic phase expressing the release from the granules, which was denoted as the SFOR model. The integral form of the SFOR model is as follows18):

  
(1)

where Ct is the pesticide concentration in the PW at time t, Cdiss is the dissolved concentration (mg/L), kr is the release rate from the granules (1/day), and ke is the decrease rate of the concentration in the PW (1/day). The goodness-of-fit of the fitted model was evaluated visually and statistically. As a statistical measure, the χ2 test was used to evaluate the agreement between the calculated and observed values.19) Next, the similarity of the dissipation was statistically compared between test plots with the same soil types using grouping analysis. In the analysis, as shown in Fig. 2, compared datasets were first fitted to a different parameter (separated) model and an entire or partial common parameter (grouped) model using the parameters (Cdiss, kr and ke) of the SFOR model. The two models were then compared by means of a one-way analysis of variance (ANOVA) F-test.20) If the null hypothesis was accepted with a p-value above 0.05, the grouped model was adopted, which meant the success of the grouping regarding commonly fitted parameter(s); otherwise, the individual model was adopted.

Fig. 2. Example of grouping analysis: (a) separated models; (b) grouped model. The data with closed form means “<LOQ (limit of quantitation)”.

Table 2 presents the results of the grouping analyses. The results between the LA and LV plots showed that 3 and 7 out of the 9 pesticides were entirely grouped, respectively. For datasets that failed in the entire grouping, significant differences (p<0.05%) were observed with respect to Cdiss. Nevertheless, the quantitative fact that 56% of the entire process and approximately 90% of the decrease phases in the dissipation between the lysimeters ensured high data reproducibility between replicates. The results from the lysimeters and paddy fields indicated that no pesticides were grouped entirely. Based on this process, it was found that the lysimeters could simulate dissipation in actual paddy fields with accuracies of 10% for Cdiss, 34% for kr, and 48% for ke. Again, the low grouping rate was remarkable for Cdiss. This difference is possibly associated with individual or combined effects of the experimental variables, such as the hydrological components, physicochemical properties of pesticides, and the effect of the experimental design of the lysimeter.

Table 2. Result of grouping analysis
Compared dataa) Alluvial soil plot Volcanic ash soil plot
C diss k r k e C diss k r k e
Lysimeter vs. Lysimeter 3 6 8 7 7 8
Lysimeter vs. Paddy field 2 9 15 3 8 9

a) The numbers of analyzed data for Lysimeter vs. Lysimeter and Lysimeter vs. Paddy field are 9 and 25 for each soil type, respectively.

2. Extraction of information from experimental data21)

2.1. PCPF-1R model

To clarify the factors affecting the experimental performance of the lysimeter simulating the actual conditions, the preferred choice was to quantify the contributions of the experimental variables using the model calibrated to the site-specific conditions. In this case, rather than a lower-tier model such as the SFOR model, a higher-tier mathematical model was more suited for the integration and abstraction of the experimental variables. Therefore, this study adopted the PCPF (pesticide concentration in paddy field) -1 model, a well-validated higher tier mathematical model for simulating the fate and transport processes of paddy pesticides.22,23) The difficulty in using the original model is that the calibration of the model-specific parameters has been mostly derived from experiments designed for model validation; therefore, users may face serious workloads during model calibration, especially for the simulation of new compounds. To avoid this problem, the governing equations of the original model were first modified to reduce the model-specific parameters and substitute more general terms. Subsequently, a new code was developed in an open software R language as the PCPF-1R model.

2.2. Novel inverse analysis protocol

The second attempt to overcome the model calibration problem was to introduce automatic calibration, which minimizes the objective function between the simulated and observed values based on the mathematical algorithm. This calibrated result ensures reproducibility and objectivity, and the quality of the results can be further refined by coupling with sensitivity and uncertainty analysis.24) Recently, the open software R has provided the package “FME” that enabled us to easily implement such inverse modeling.25) Therefore, the automatic calibration of the PCPF-1R model could be implemented by manipulating the R resources including “FME.”

Here, highlights of the developed R-based inverse analysis protocol are briefly explained. The most daunting but important task is the parameterization of the physicochemical properties of the target compound because users have to acquire reliable data manually from existing literature or databases (appropriateness depends on user subjectivity). The prior parameter uncertainty is then visualized using multiple parameter sets randomly selected by a Monte Carlo simulation (MCS) to confirm that this sufficiently covers the observed data. Subsequently, the model parameters sensitive to the objective function are filtered by the standardized rank regression coefficients (SRRCs) based on multiple regression analysis using the outputs of MCS.26,27) After the parameter identifiability check, model fitting to the observed data, using the pseudo-random search algorithm (PseudoOptim),28) is carried out to optimize the filtered parameters. In the final step, the parameter uncertainties associated with the experimental data are investigated using the Markov chain Monte Carlo (MCMC) method. In this protocol, the adaptive Metropolis algorithm is used to improve the acceptance/rejection efficiency by tuning the covariance of the jumping (proposal) distribution based on the history of the chain at a certain frequency.29) Setting the prior distributions and starting values using the result of model fitting, a sufficient number of MCMC sampling is administered with single or multiple chain(s) and an appropriate burn-in period. The calibrated parameters are expressed by posterior distributions, and the posterior parameter uncertainty is visualized through MCS.

2.3. Application example

The applicability of the developed R-based inverse analysis was verified using the dissipation data of the herbicides simetryn and molinate, obtained in 2012. The post-inverse analysis (calibrated result) for simetryn is shown in Fig. 3, which clearly indicated the reduction of the parameter uncertainty and good agreement of the calibrated simulation with the observed data through the inverse analysis. Moreover, the fluctuations of the concentrations in the FA/FV plots due to variations of the paddy water depth were well described. Next, examples of posterior distributions, given as the equilibrium sorption coefficient (Kd) for simetryn and volatilization rate constant (kL–A) for molinate, are shown in Fig. 4. While the high variation of Kd for simetryn, depending on soils, was consistent with the literature data,30) faster volatilization tendency was found in the FA/FV plots for molinate with high vapor pressure.

Fig. 3. Observed data and calibrated simulations of concentrations in paddy water for simetryn in each test plot, where the calibrated simulation is given by the parameter sets that give the highest probability within the chain, light and dark gray bands are prior and posterior minimum–maximum ranges of MCS outputs associated with site-specific parameter uncertainties (modified from Kondo et al.21) with permission from Wily, Copyright 2019).
Fig. 4. Posterior distributions of equilibrium soil sorption coefficient (Kd) for simetryn and volatilization rate constant (kL–A) for molinate (modified from Kondo et al.21) with permission from Wily, Copyright 2019).

As shown in Table 1, two surface runoff events were observed in the FA/FV plots in 2012, and a decrease in the concentrations in PW at 9 and 11 DATs was confirmed (Fig. 3). The mass balance of the calibrated simulation was analyzed to quantify the runoff amounts of the two herbicides. Figure 5 shows the mass balance of the two herbicides at 21 DAT. According to this result, the cumulative runoff amounts of simetryn and molinate were 2.8–10.7% and 0.9–1.1% of the applied mass, respectively, and a non-negligible effect was found in simetryn in the FA plot. This analysis also revealed a remarkable difference in the dissipation pathways between the lysimeter and paddy fields.

Fig. 5. Mass distributions of simetryn and molinate at 21 days after treatment (DAT) calculated from the calibrated simulation (modified from Kondo et al.21) with permission from Wily, Copyright 2019).

It was considered that the experimental setting, such as the daily percolation rate in the lysimeter, may cause the difference in the dissipation pathways from those of paddy fields; therefore, a simple case study was performed using the measure of a time weighted average concentration (TWAC).27) First, the TWACs at all test plots were calculated from the posterior MCS outputs. Table 3 presents the results for 72 hr. The TWACs at the lysimeters were underestimated in the actual paddy fields. Second, the daily percolation amounts in the water balance data from the LA/LV plot were substituted with those averaged from the FA/FV plot, and then the posterior MCSs were again run to calculate the TWAC as the modified lysimeter case. These values were equivalent to the actual paddy fields; thus, the setting of the daily percolation rate in the lysimeter experiment was found to be the key parameter for the extrapolation of the lysimeter test as the dissipation simulator of actual paddy fields.

Table 3. Summary of time weighted average concentration (TWAC) at 72 hr
TWAC (mg/L)
Paddy field Conventional lysimeter Modified lysimeter
Simetryn
Alluvial soil 0.45±0.010 0.21±0.006 0.48±0.010
Volcanic ash soil 0.41±0.010 0.19±0.003 0.35±0.005
Molinate
Alluvial soil 2.02±0.032 1.00±0.023 2.09±0.056
Volcanic ash soil 1.79±0.040 1.12±0.100 1.95±0.159

3. Relating different scales of experimental data31)

3.1. First step: inverse modeling of laboratory experiment

The discussion in the previous sections on the application of the mathematical model and its inverse analysis only focused on the pesticide concentration in PW, owing to experimental constraints. As a result, the fate and transport parameters in paddy soils relied only on literature data. The literature base parameterization, on the other hand, may cause over-calibration of the model during the inverse analysis. Recently, mathematical modeling studies that assess the persistence of pesticides in water-sediment test systems have been reported.32,33) Inspired by these studies, inverse modeling of the paddy soil test under laboratory conditions was attempted in order to extract the unobtainable parameters in outdoor experiments such as lysimeters and paddy fields.

A container test with flooded soils collected from the test plots was conducted to target four herbicides: daimuron, fentrazamide, bromobutide, and bensulfuron-methyl. Mixed analytical standards of the four herbicides were applied to replicate the test systems after two weeks of pre-incubation. At 28 DAT, each test system was separated into a water phase (WP) and soil phase (SP), and the concentration of each phase was determined by the validated methods. The analytical concentrations were subjected to the inverse analysis of the PCPF-LR model, a simple two-compartment mathematical model that is structurally compatible with the improved PCPF-1R model. The procedures for inverse analysis were based on a previously developed protocol.

The calibrated simulations using the PCPF-LR model are projected at the top of Fig. 6. It was confirmed that the model appropriately captured the decreasing trend of the observed data. From these results, the model calculated the apparent Kd (Kd, app) values, defined as the ratio of the concentration in the SP to that in the WP, to monitor the time-dependent sorption state. The Kd, app values for daimuron, shown in Fig. 7, increased with time, as previous studies pointed out18,34,35) and then reached a plateau, the level of which was within or close to the range derived from the OECD 106 batch method. Thus, the present approach is also useful for assessing the soil sorption phenomena of pesticides. Note that the cause of wide q5–q95 range in Kd, app at the FV soil was associated with relatively higher variation of the concentration in the WP comparing to the other soils due to the posterior parameter uncertainty.

Fig. 6. Plots of observed data versus calibrated simulations for daimuron in laboratory and outdoor experiments at each test soil, where q5 to q95 correspond to the 5th to 95th quantiles of MCS outputs. Closed symbols in measured data stand “<LOQ (limit of quantitation)” (reprinted from Kondo et al.31) with permission from Wily, Copyright 2020).
Fig. 7. Measured and simulated time course of apparent sorption coefficient (Kd, app), where q5 to q95 correspond to the 5th to 95th quantiles of MCS outputs. Dashed lines show the upper and lower KF values derived from the batch experiment (modified from Kondo et al.31) with permission from Wily, Copyright 2020).

3.2. Second step: parameter transfer to outdoor experiments

To utilize the calibrated results of laboratory modeling for outdoor modeling as the second step, the parameters related to the degradation in WP and SP, phase transfer between the two phases, and soil sorption were randomly sampled from the respective posterior distributions. A total of 100 parameter sets were created and transferred to the outdoor model using the PCPF-1R model. With these parameter sets, PseudoOptim was iteratively run to calibrate site-specific parameters, such as the dissolution rate from granules, photolysis rate, and instant partitioning between PW and paddy soil.

The second step was verified using the experimental data of the four herbicides mentioned in section 3.1, applied as granules in 2012 and as a flowable mode in 2013. The post-processed simulations reasonably described the dissipations of the target herbicides, depending on the formulation type (see middle and bottom of Fig. 6). To compare the dissipation characteristics of the four herbicides between the formulations, 1 : 1 comparisons were performed of the dissipation indicator (DT50), derived from the kinetic models and that of persistence (DegT50), derived from the PCPF-1R and PCPF-LR models. As indicated in Fig. 8, while significantly faster dissipation was found in the flowable mode (p-value=0.004), the degradation was closer to the 1 : 1 relationship (p-value=0.1). This result suggests that the persistence of the target herbicides in PW was not affected by the formulation type, although their overall dissipation patterns were significantly different.

Fig. 8. 1 : 1 comparison of time required for 50% dissipation (DT50) and degradation half-lives (DegT50) in paddy water (PW) under flowable application (subscript, “PW-Flowable”) and under granule application (subscript, “PW-Granule”) for outdoor experiments (modified from Kondo et al.31) with permission from Wily, Copyright 2020).

Concluding remarks

By means of mathematical modeling with inverse analysis, our interpretation of the information hidden in the experimental data was facilitated. This study successfully demonstrates the utility of this method through a case study of pesticide dissipation in paddy test systems from three perspectives. As shown in Fig. 1, laboratory modeling can be extended to the OECD 307 aerobic flooded soil test, using labeled compounds. The inversely derived parameters, including transformation processes from OECD 307 data, can be transferred to predict outdoor metabolite modeling under various formulation types. Finally, these in-field modeling results calibrated from the registration data can be exported to the modeling of the post-registration data taken from the edge-of-field36) and catchment scales.37) The iterative runs of the aforementioned experiment and modeling can help cost-optimized experimental design.38)

Traditionally, researchers have developed mathematical models for predicting pesticide fate and transport using user-friendly environments such as MS Excel®. The present study also attempted to build all modeling materials using open software R, and the main codes were archived in a repository (https://github.com/k-kondo-IET/PCPF-1R). In the upcoming open science era, the developed models, protocols, and their codes are expected to contribute to the innovation of the environmental pesticide science field.

Acknowledgments

The author is pleased to be recommended by the Pesticide Science Society of Japan for this honorable Society Award. This study was conducted as part of the Continuing Research Project in compliance with the Expenditure Plan for Public Interest Purposes of the Institute of Environmental Toxicology (Corporate code: A012055). The author sincerely thank Dr. Yasuhiro Yogo, Prof. Hisashi Miyagawa, and Prof. Masahiro Natsume of the Advisory Panel for their valuable suggestions, cooperation, and technical support for our study. The author also expresses thanks to the Japan Association for Advancement of Phyto-Regulators for their cooperation and advice during the field experiments. The author expresses thanks and respect to Prof. Hirozumi Watanabe for his teaching, especially in mathematical modeling.

References
 
© Pesticide Science Society of Japan 2022. This is an open access article distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) License (https://creativecommons.org/licenses/by-nc-nd/4.0/)

This article is licensed under a Creative Commons [Attribution-NonCommercial-NoDerivatives 4.0 International] license.
https://creativecommons.org/licenses/by-nc-nd/4.0/
feedback
Top