A novel model concept for modelling the leaching of natural toxins: results for the case of ptaquiloside

Abstract

Phytotoxins are a large class of highly diverse emerging environmental contaminants that have been detected at high concentrations in plants, water and soils. This study presents a novel modelling approach for assessing the fate of plant toxins in the soil–plant–atmosphere continuum, developed for the specific case of ptaquiloside (PTA), a carcinogenic phytotoxin produced by Pteridium aquilinum. The mechanistic model DAISY has been adapted for reproducing phytotoxin dynamics in plants, covering processes such as toxin generation in the canopy, wash off by precipitation and toxin recovery in the canopy after depletion events. Transport of the toxin in the soil was simulated by the advection–dispersion equation assuming weak sorption and degradation for two Danish soils. The model simulates realistic toxin contents in the plant during the growing season, where the actual PTA content is dynamic and a function of the biomass. An average of 48% of the PTA produced in the canopy is washed off by precipitation, with loads in the soil often in the order of mg m⁻² and up to a maximum of 13 mg m⁻² in a single rain event. Degradation in the soil removes 99.9% of the total PTA input to the soil, while only 0.1% leaches into the soil. The median annual flux-averaged predicted environmental concentrations during single events are often in the order of μg L⁻¹, reaching up to 60 μg L⁻¹ for the worst-case scenario. The simulated results for both degradation and wash off are of the same order of magnitude as the published data. Based on the results, we conclude that DAISY, with the newly implemented processes, is a useful tool for understanding, describing and predicting the fate of PTA in the soil. Further work comparing the model results with real data is needed for the calibration and validation of the model.

---

Environmental significance

The identification of a wide number of phytotoxins as potential aquatic micropollutants has highlighted the need for environmental models to assess the fate of phytotoxins in the environment. In our study, we have implemented relevant processes in the mechanistic model DAISY to account for toxin generation in plants, wash off after precipitation and recovery after depletion. We tested the new modelling concept for the case of ptaquiloside, one of the three terpenes produced by bracken fern (Pteridium aquilinum) with carcinogenic properties. The results show that PTA can leach in relatively high concentrations in the order of μg L⁻¹, which is comparable to those found in pesticides.

1. Introduction

Phytotoxins are defined as toxic plant secondary metabolites that play an important role in the defense of plants against herbivores, pathogens, and other plants, increasing the fitness of plants in the environment.^1–5 In mid Europe more than 500 phytotoxins have been categorized as potential aquatic micropollutants based on their high mobility, high toxicity, continuous production and environmental loads, which in some cases are in the order of kilograms per hectare.^4,6–10 Most of the studies dealing with natural toxins have focused on estimating physicochemical properties, plant concentrations, degradation in different matrices and leaching. Such studies deal with several chemical classes such as terpenes and terpenoids,^6,11,12 flavonoids^9,13,14 and alkaloids.^10,15,16 Phytotoxins have been detected in soils,^9,10,17 groundwater^9,18,19 and surface water.^10,14,20 Moreover, release mechanisms from the source plant such as root exudation^16,21 or wash off by precipitation have been identified.^20,22,23 The reported environmental concentrations are comparable to concentrations of industrial pollutants and fall in the range of ng L⁻¹ to μg L⁻¹.

Bracken fern is one of the most abundant vascular plants in the world, present in all continents but Antarctica.³ It is highly invasive and usually present as a predominant species in forest glades and open land.^3,25 Bracken is also known to produce a wide variety of well-known toxic substances, of which terpenoids are the main contributors to the high toxicity of bracken.^25,26 Terpenoids in bracken are linked to several chronic diseases in animals and humans such as enzootic hematuria, blindness and tumors.^27–30 Ptaquiloside (PTA) is one of the sesquiterpene glycosides produced by bracken. It was first identified in 1985 as the main active compound responsible for the carcinogenicity of bracken.³¹ Ptaquiloside is highly hydrophilic, water soluble and weakly sorbing, as implied by its low octanol–water partitioning coefficient (K_OW) of 0.23.¹⁰ Ptaquiloside may also bind to soil organic matter in the form of a bound residue. Ptaquiloside is degraded rapidly in water and soil due to microbial degradation and hydrolysis under equilibrium conditions.^32–34 Degradation of PTA under dynamic conditions, i.e. while being transported, might present different rates compared to that under static conditions. Hydrolysis rates are highly pH dependent, with the lowest degradation rates at neutral to slightly alkaline pH.³⁵ Ptaquiloside is produced in amounts up to 1% of plant dry weight (DW),³⁵ with estimated maximum contents (90^th percentile) reaching up to 20 kg ha⁻¹ in England.³⁷ The toxin content in bracken has been suggested to depend on the bracken subspecies (geographic region), growth stage, environmental factors such as altitude or temperature,^36–38 and previous management history, such as browsing of animals or cutting by humans.^39,40 Ptaquiloside has been detected in soils, surface water and groundwater in areas with bracken coverage, as well as in milk and meat in animals browsing on ferns.^17,20,41,42 Ptaquiloside has been detected in throughfall water from the canopy during precipitation events and also after applying water aerosols to bracken pinnae in the lab, with concentrations in the order of μg L⁻¹.^20,43 These findings have demonstrated that PTA can be washed passively at high concentrations from bracken leaves by simple contact with water. One hypothesis is that release by precipitation might be related to the presence of touch-sensitive glandular trichomes on the surface of the plant, which are known structures producing a wide array of natural products including terpenes.^37,44–46

Environmental models describing the soil–plant–atmosphere continuum are powerful tools for assessing the environmental fate of organic pollutants.⁴⁷ There are a large number of well validated models for pesticides such as HYDRUS,⁴⁸ MACRO,⁴⁹ PEARL,⁵⁰ PELMO⁵¹ or DAISY,⁵² but none of the existing models have been adapted to simulate the fate of natural toxins. Modelling the fate of natural toxins presents several challenges compared with traditional organic pollutants such as pesticides. The most significant challenges are the actual description and quantification of toxin generation in plants, the toxin release mechanisms from the plant and lack of environmental distribution data. There is at present no environmental model that can integrate plant growth with toxin generation within the plant. To the best of our knowledge, Ramwell et al. (2010) made the first and only attempt to model the fate of PTA in the environment. In their study, PTA was added to soils by spraying, identical to how pesticides are sprayed on top of the soils, and it was this application of PTA (release) that was identified as the variable with the largest uncertainty.

In the present study, we present a new modelling concept to simulate the fate of phytotoxins at the pedon-scale (square meter), applied to the specific case of bracken fern and PTA under Danish conditions. The agroecosystem model DAISY was used, as it offers the advantage of a relatively detailed description of water balance, solute transport and plant growth. This model is the first step for the description of surface processes and identification of key variables needed for simulating the fate of natural toxins in the environment.

We modified DAISY to account for toxin dynamics in the plant based on published information, with the goal of improving the suitability of the model for the specific case of natural toxins. The key processes in the newly implemented functions comprise (i) toxin generation in the canopy as a function of the plant's biomass, (ii) wash off of PTA from the plant's leaves by rainwater, and (iii) PTA recovery in the canopy after a depletion event. Sorption, degradation and transport in soil is modelled using the existing submodules in DAISY, parameterized for PTA. In the first step we present an overview of DAISY, with special focus on the newly implemented processes. Thereafter, the model is used to assess the fate of PTA in two soil types with widely differing soil hydraulic properties and known areas in Denmark with the presence of bracken. Finally, we discuss the model performance and seasonal PTA dynamics in the field.

2. Materials and methods

2.1. The DAISY model

DAISY is a physically based 1D/2D model describing the soil–plant–atmosphere system.^53,54 The model has been developed and applied in various studies concerning plant growth and solute transport, with a focus on nitrogen^55,56 and pesticide dynamics and leaching.^52–57 The main processes in DAISY are water and heat flow, nitrogen and general solute transport, carbon dynamics, plant growth and soil management. Water and solute transport in DAISY is based on the Richards equation⁵⁸ and the advection–dispersion equation,⁵⁹ respectively. For a more detailed description of the model's water flow and solute transport,^53,60 for nutrient cycling⁶¹ and plant production,⁵⁴ the reader is referred to the literature. In the next paragraphs we will only present the new processes implemented for PTA and some key processes related to the transport of PTA in the soil.

A great advantage of DAISY compared to other models is a well described plant growth module, which simulates daily plant production, phenological development and canopy development. These processes allow us to link toxin generation with biomass production and toxin wash off from the canopy by precipitation. The net production of biomass is determined from the gross canopy photosynthesis, the partitioning of nutrients between different parts of the plant and respiration rates. The phenological stage of the plant is quantified by the development stage (DS), which in DAISY has a value of 0 at emergence and 2 at maturity.⁵⁶ The canopy structure is defined by the leaf area index (LAI), which corresponds to the green leaf area per soil surface area [m² m⁻²].⁵⁴ The LAI is an important factor for growth, determining the capacity for photosynthesis, and the water interception on the canopy, which is a very important factor since it drives processes such as the wash off of toxins from the canopy.

In DAISY, a chemical present in the biomass–soil–water system can either be degraded on the plant canopy, be degraded in the soil, sorb to the soil, be taken up by plant roots or leach into deeper soil layers. DAISY describes degradation of chemicals in soil water with first order kinetics. The actual degradation rate is calculated starting with a reference degradation rate (20 °C) which is further scaled based on different abiotic factors. The calculation of the actual degradation rate ξ [T⁻¹] is described as:

ξ = −[k* f_T(T) f_h(h) f_Z(Z)] (1)

where k* [T⁻¹] is the base degradation rate under standard conditions; T [°C] is the temperature; h [L] is the soil water pressure head and Z [L] is the soil depth. The f coefficients viz. f_T, f_h, and f_Z are modifier functions for temperature, soil water pressure head and soil depth, respectively.⁵³

2.2. Implementation of the toxin generation function

A linear function for generating toxins in the canopy as a function of plant biomass has been implemented in DAISY. This was based on biomass and PTA data derived from the reports of Rasmussen and Hansen (2004)⁶² and Rasmussen et al. (2003)⁴³ for three Danish locations and it is described by:

A′(DS) = af_DS (3)

where M_MAX [M L⁻²] is the potential toxin content in the canopy per unit area; B_CAN, [M L⁻²] is the above-ground biomass, calculated using DAISY's plant growth model, and a [—] is a constant. Constant a determines the relationship between the toxin content and the biomass, i.e. mass of toxin per kg of biomass (Fig. 1a). A′(DS) is the a coefficient scaled based on the plant DS factor (f_DS), which is specified by the user. This function (eqn (3)) is used for controlling toxin generation at the beginning of plant growth and toxin decay at the latest stages of plant development.

Fig. 1 Linear relationship for defining the potential maximum canopy PTA content (M_MAX, eqn (3)) as a function of biomass (a) and illustration of toxin recovery after a wash off event at different recovery rates (r, eqn (5)) (b). MIN and MAX represent the lower and upper limit values for M_MAX and r.

Fig. 1a shows examples of this linear relationship. However, the actual toxin canopy content might be lower than M_MAX due to toxin wash off by precipitation. The amount of toxin washed off from bracken leaves is calculated as a function of the water dripping off the canopy (eqn (4)):

where W [M L⁻²] is the amount of toxin washed off from bracken by precipitation; M_t [M L⁻²] is the actual toxin content; f_WC [—] is the canopy wash off coefficient; D [L³ L⁻²] is the amount of water drip-off and I_C [L³ L⁻²] is the water interception capacity of the canopy. f_WC is a coefficient ranging from 0 to 1, where 1 means that the highest fraction of the toxin content in the canopy is washed off and 0 represents no wash-off at all. Please note that wash off is based on the amount of drip off water, and thus it will not be possible for all the toxin to be washed off under realistic conditions for precipitation and wash off coefficients values. Moreover, if the volume of intercepted water in the canopy does not surpass the canopy storage capacity, water will not drip off the canopy and hence there will be no PTA wash off. Eqn (4) is similar to the canopy wash off function in the FOCUS guidelines, with the difference of calculating the water intercepted in the canopy for a dynamic LAI. Similarly, toxin wash off from litter is calculated in the same way as for the canopy (eqn (4)).

The second process is the recovery of the toxin in the canopy after a wash off event (eqn (5)). The rationale of this process is that there is a fraction of available toxin for wash off on the plant surface that suffers a depletion in content during a rain event. For a hypothetical depletion event and assuming a constant crop biomass, i.e. no biomass production, toxin is recovered in the canopy based on the difference between the actual canopy content (M, [M L⁻²]) and the potential content (M_MAX):

where r [T⁻¹] is the recovery rate parameter, describing how fast the toxin content is recovered in the canopy after a release event. Fig. 1b illustrates how the toxin content in the canopy is replenished after a hypothetical wash off event. Combining eqn (2) and (5), we can calculate the actual toxin generation rate in the canopy:

The actual toxin content in the canopy on day t is calculated as the toxin content the previous day (t₀), plus the toxin generation (M_t, eqn (4)), minus the toxin washed off that day (W_t, eqn (3)).

2.3. Model parameterization

In this section we present the parameterization of the DAISY model for this study. Firstly, we introduce the parameterization of the plant growth module for bracken. Secondly, we describe how the uncertainty analysis of the model parameters for the newly implemented processes is performed. Lastly, we provide a brief description of the scenarios used in the simulations.

2.3.1. Bracken parameterization. The bracken module was parameterized to reproduce the maximum reported biomass of bracken observed in Denmark (550 g m⁻²).¹⁷ The LAI was defined to have a seasonal maximum value ranging between 3 and 4, based on published values for bracken.⁶³ The canopy interception capacity was set to 0.467 mm, which is the maximum storage of water per unit of the LAI.⁶³ This storage capacity is proportional to the LAI and hence is scaled up for a growing LAI during the growing season. The parameterization of the model focused on generating a realistic biomass and LAI, which are the two variables that most strongly influence toxin generation and rain wash off.

The parameterization of the plant module in DAISY aims at simulating the phenology and nutrient partitioning of bracken, and it is based on a bracken growth model reported by Pakeman et al. (1994).⁶⁴ The date of emergenence is determined by temperature with a value of 100 °C days, which corresponds to the sum of daily air temperature of days with an average temperature higher than 4 °C.⁶⁴ Partitioning of nutrients is parameterized so bracken, at early development stages (DS < 0.6), grows rapidly, helped by the translocation of carbohydrates from the rhizome pool reserves towards the young fronds. When bracken reaches maturity (DS 1), the nutrient partitioning shifts towards replenishing the carbohydrate reserves in rhizomes for the next year.

The senescence of bracken in DAISY is regulated by the average air temperature, taking place on the first day during the growing season with an average temperature lower than 2 °C, usually in October or early November depending on the year. Death of above-ground bracken biomass is simulated by a management action of cutting the plant at a 40 cm height and all biomass above this height is deposited onto the soil surface as litter. The remaining part of the plant (<40 cm) stays dormant until the next growing season. The decision of not cutting the entire plant was due to limitations for achieving annual growth from rhizome reserves in DAISY with no above ground biomass left. The litter layer is formed by the cut plant residues and its characteristics are determined by the weight of the deposited plant residues and the specific area coefficient (m³ kg L⁻¹), which is defined by the user and represents the area covered per litter mass. Water and PTA dynamics in the litter layer are described by a simple mass balance approach, defined by the water storage capacity (L kg⁻¹) and specific leaf area coefficients of the litter.⁶⁵

2.3.2. Parameter selection and uncertainty analysis. The three newly implemented parameters along with the degradation rate of PTA in the soil were selected for uncertainty analysis of the model: toxin content relationship constant (a, eqn (2)), canopy wash-off coefficient (f_WC, eqn (3)), PTA recovery rate (r, eqn (4)), and degradation in the soil (k*, eqn (1)). Table 1 lists the four parameters and their limit values. The parameter limits, except for r, have been determined based on the literature where only studies of bracken fern (var. aquilinum) in northern Europe (Denmark, England and Ireland) have been considered.^11,17,20

Table 1 Limit values for the four parameters used in the simulations

Param.	Min	Max	Unit
PTA constant (a)	0.0182	0.4727	—
Wash off (fWC)	0.15	0.76	%
Recovery (r)	4.2 × 10^−4	1	h^−1
Degrad. rate (k*)	0.0239	0.463	h^−1

---

The maximum potential PTA content in the canopy is 260 mg m⁻² as determined from the slope of the toxin relationship with biomass (a), with values defined in a previous survey of Danish bracken populations by Rasmussen et al. (2003). Although a study on an Irish population has found higher PTA contents,⁴⁰ we decided to proceed only with data from Danish bracken populations. The canopy wash off coefficient (f_WC) is in the range estimated in a UK study of PTA stream concentrations measured during precipitation events.²⁰ In this study, it was estimated that between 0.15 and 0.76% of the total PTA canopy content was washed off during precipitation events.²⁰ This range of values was estimated by calculating the fraction of total PTA that was washed off. Therefore, the coefficient gives a hint of the fraction of PTA that is available for wash off at the moment of precipitation. Due to the lack of direct measurements for wash off from litter, we assumed an equal value for litter wash off as for the canopy.¹⁸ For the recovery rate (r), due to the lack of observations with appropriate temporal resolution, values for full recovery ranging from 1 hour to 100 days were selected. These two cases were selected to represent extreme situations where there is instant recovery of toxin content (1 hour) and no significant recovery of canopy content between precipitation events (100 days). Toxin washed by precipitation from the canopy and litter are the only sources of PTA input to the soil in the current model setup. We do not consider root exudation in our model due to lack of scientific evidence of this process taking place in bracken.

Degradation of chemicals in the soil profile is calculated by applying depth, temperature and water content factors defined in FOCUS (FOrum for the Co-ordination of pesticide fate models and their USe) scenarios. FOCUS provides guidelines for environmental fate studies dealing with persistence and degradation kinetics of pesticides under the European Union registration.⁶⁶ The values for PTA degradation rates in soil are taken from a study dealing with PTA degradation in different soils from Denmark, where we considered the A-horizons only for calculating the base values (k*, eqn (1)).¹¹ The degradation rates in the soil profile were calculated applying the base values, followed by modifiers accounting for soil depth, temperature and water content recommended in the FOCUS guidelines.⁶⁷ Degradation of PTA is assumed to not take place in the living canopy, and thus takes place only when it is released by precipitation or after senescence in the litter. Sorption of PTA in the soil was calculated with a Freundlich isotherm, using a Freundlich constant (K_F) and exponential coefficient (n) of 0.17 and 0.86, respectively.¹¹ For tackling the uncertainty in the input parameters, a total of 100 simulations for each of the two soils were performed, with combinations of the values for each of the four parameters randomly selected by a Latin Hypercube Sampling Strategy (LHS), assuming a uniform distribution (no prior knowledge).

2.3.3. Simulated scenarios. We modeled the transport of PTA from the canopy and into a 2 m soil profile. We assumed that water and solute transport occurs along the vertical axis (depth), e.g. no runoff. We tested the model for two Danish locations, Humleore and Præstø, both located in Zealand. At these two sites, bracken is dominant in open areas with a dense canopy surpassing 2 meters of height in some spots (Rasmussen et al., 2004; author, personal observation). The weather for the entire simulation period is from the weather station of the University of Copenhagen in Taastrup, located 20 km west of Copenhagen (Zealand, Denmark). The climate is mild and temperate, with significant rainfall and the average temperature above the freezing point during the winter. The mean temperature of the weather input for the period 1985–2005 is 8.2 °C, with an average precipitation of 913 mm. Precipitation is distributed throughout the year, with maximum precipitation during the summer period.

The time period of the simulation was from 1980 to 2005, where the first five years were used for “warming up” the model, thus the period used for the calculations accounts for a total of 21 years. The warming up period is needed to reach stability in plant production, i.e. for building up root reserves and for building up PTA pools in the soil.

Præstø and Humleore present soil profiles of contrasting hydraulic properties, being categorized by the USDA as sandy and sandy loam, respectively. Præstø has a soil profile formed on unconsolidated sand deposits, with high hydraulic conductivity. On the other hand, Humleore has sandy loam with a high content of silt and clay in deeper horizons, with preferential transport via macropores. In Humleore, macropores are defined at depths between 0 and 100 cm, with a diameter of 4 mm and a density of 20 macropores per square meter. The physicochemical characteristics of the two soils used in the model input are shown in Table 2.

Table 2 Soil profile and texture for Præstø (Rasmussen et al., 2003) and Humleore (author, unpublished)

	Horiz.	Depth (cm)	Clay (%)	Silt (%)	Sand (%)	pH
Præstø (sandy)	A	32	6	1	93	3.3
	E	53	1	1	98	3.4
	B	118	2	0	98	3.8
	C	200	0	1	99	4.2
Humleore (sandy loam)	A	20	11	12	77	3.7
	AE	36	10	10	80	3.7
	B	70	18	26	56	5.7
	C	200	24	20	56	4.5

---

The soil profile is unsaturated and we assumed a fixed water table and a water pressure head equal to 0 cm at the bottom of the soil (2 m). This decision followed recommendations for standard parameterization of soils in Denmark in cases where no literature data are available.⁶⁸ The groundwater does not have an influence on the solute mass balance, acting simply as a sink of PTA.

3. Results and discussion

3.1. Bracken growth

DAISY simulations result in a median dry biomass production of 519 g DW m⁻² and up to 590 g DW m⁻². The variation between years in the simulated biomass is determined by the soil water content and amount of sunlight. The LAI is proportional to the biomass, reaching a seasonal maximum of 3.4 ± 0.3 at the end of the season. The simulated date of emergence and senescence takes place during April and November, respectively. The simulated biomass dynamics for a year is similar to what is found for bracken in the literature, i.e. fast initial growth, maximum biomass values during summer, nutrient partition within the plant's life cycle and the LAI.^17,24,64 Simulated maximum bracken biomass agrees with measured biomass values from three different Danish sites (Fig. 2a).⁶² The parameterization of the plant module is not calibrated and some parameters need to be adjusted to better describe the date of emergence, biomass production in the early stages of plant growth and date of senescence.

Fig. 2 Simulated results for biomass, PTA content in the canopy, wash off and leaching concentrations are shown in the figure. On the left side is shown the surface processes implemented in the model, the median biomass (a) and PTA in the canopy (b) for all the years and simulations (n = 2100). Symbols in the (a) and (b) subplots are the real data for bracken biomass and PTA content in the canopy, respectively. Each symbol represent a different Danish location (data derived from the reports of Rasmussen et al., 2003 and Rasmussen and Hansen, 2004). On the right side of the figure is shown the mass of PTA washed off by precipitation (c) and the PTA concentration at a 100 cm depth (d) for a specific year in Humleore (n = 100). The lines represent the median value, while the shaded area represents the upper (90^th) and lower (10^th) percentiles.

3.2. Toxin generation and wash off

The simulated maximum production of PTA content in the canopy has a median value of 110 mg PTA m⁻², in accordance with results reported by Rasmussen and Hansen (2004) (Fig. 2b, symbols). The content of PTA in the canopy is positively correlated with biomass until maturity, when the toxin content decreases before senescence, as seen in previous field studies (Fig. 2b). Precipitation has a big influence on the PTA contents in the canopy due to the wash off effect. The masses of PTA being washed off are typically in the order of several mg m⁻² per rain event (Fig. 2c), the same order of magnitude as that in previous findings.^18,43 The release of PTA from the bracken canopy is highly variable between years and amounts to 82 ± 48 mg m⁻² per year, with the difference between years due to variability in the amount and timing of precipitation. The highest toxin wash off is seen for a precipitation event of 60 mm in June, with a fully developed canopy and toxin contents with peak values, leading to a median wash off of 15 mg m⁻² in a single event and up to 36 mg m⁻² in the worst-case scenario.

After senescence of bracken, the amount of toxin left in the canopy is transferred to the litter. Ptaquiloside present in the litter can be released onto the soil and also, in contrast with the PTA in the canopy, undergoes degradation. The annual PTA washed off from the litter by precipitation is 14 ± 6 mg m⁻² (n = 2100). The depletion of the PTA content in the litter takes place relatively fast after senescence (<4 weeks) due to both degradation and wash off. During wet periods, litter becomes water saturated which leads to a constant release of PTA from the litter to the soil surface until the content is depleted completely.

In reality, bracken litter tends to form a loose horizontal layer 20–30 cm above the soil surface.⁶⁹ Therefore, the simulated water content of the litter might differ significantly to what is found in reality. Moreover, the timing of the first frost might be an important factor contributing to the release of PTA from bracken. After repeated cycles of freezing and thawing, the permeability of the cell's membrane increases leading to the release of its contents.⁷⁰ This principle has been used for extraction of plant proteins and it is relevant in the case of release of PTA, it being highly hydrophilic. If this frost takes place early, it might release substantial amounts of toxin from the living plants to the soil surface. This process has not been studied yet and hence it is not included in the model.

The wash off coefficient used in the model is a parameter that has been estimated just in one location although not directly measured. Therefore, the magnitude of this coefficient is uncertain. In reality, this coefficient might be variable at different development stages and depend on PTA content (g PTA per g DW), rather than on the total content (mg PTA per m²). If this is the case, the wash off coefficient will have a maximum value at the beginning of the season and decrease throughout the season. Moreover, the release of PTA from the litter might be a very different process than canopy wash off. The litter wash off coefficient in DAISY can be set to a different value or even develop the description of this process in the model further. However, data are needed to either calibrate the coefficient or describe a more complex model with dynamic values. For now, we assume a similar mechanism for the litter and the canopy, following a parsimonious approach. Regarding the production of toxins in bracken, in reality there are many biotic and abiotic factors contributing to production rates such as development stage, physical stress in the form of drought, and attack by herbivores or cutting, and it may also depend on the local variety of bracken.^17,39

3.3. Toxin fate at the pedon-scale

An example of toxin dynamics in the bracken canopy and soil for a full year is shown in Fig. 3. The canopy PTA content in the 100 scenarios presents a median seasonal maximum of 130 mg m⁻², reaching a total maximum up to 245 mg m⁻² for the scenario with the highest PTA generation. The amount of toxin being washed off increases through the season for precipitation events with the same intensity due to the growing LAI (Fig. 3b). Moreover, after a strong precipitation event, and hence high PTA wash off, PTA canopy content dynamics are quite variable with peaks and valleys. This type of dynamics can also be seen in the real measurements shown in Fig. 2b.

Fig. 3 Simulated processes influencing the fate of ptaquiloside in Humleore during the year 1994. The subplots illustrate simulated plant biomass and PTA content in the canopy (PTA can (a), PTA mass washed off the canopy by precipitation (b), degradation in the soil (c), and PTA leaching at a 100 cm soil depth (d). The lines represent the median value, while the upper and lower limits of the shaded area correspond to the 90^th and 10^th percentile, respectively (n = 100).

The PTA transferred from bracken to the soil is degraded rapidly during the summer period (Fig. 3c). This might be due to two factors, the temperature increasing the degradation rates in the soil and the long residence time, leading to a high degradation of the PTA entering the soil. In the case of sandy soil, which has no macropores, leaching of PTA takes place once the soil profile becomes water saturated during autumn. In the case of sandy loam soil, when local saturation takes place in part of the soil during high precipitation events, macropores get activated and PTA is transported quickly to deeper soil layers. This preferential flow pattern can be observed for wash off events when PTA concentrations in the leachate increase abruptly (e.g., Fig. 2d and 3d). The leaching events taking place in autumn and winter (i.e. from December until April) are explained by the transport of PTA stored in deep layers of the soil by the percolating water when the soil becomes water saturated (Fig. 3d).

Table 3 provides simulated data for the annual toxin mass balance at the field scale, considering a 100 cm depth as the bottom boundary. The standard deviation is calculated from the 100 simulations performed for each soil type. Bracken generates an annual canopy PTA mass of 177 ± 85 mg m⁻². The difference in PTA generation between the two sites is due to the different biomass generated, mainly due to water stress in the case of the sandy soil. Canopy and litter wash off accounts for on average 47% and 8% of the total PTA generated in the canopy, respectively. These two sources release to the soil an average of 54% of the total PTA mass generated in the canopy, while the remaining 46% is degraded in the litter and never released to the soil. Soil degradation is the most important process determining the fate of PTA in the soil, removing 99.9% from the total PTA input to the soil. Only 0.1% of the total PTA input to the soil ends up leaching the soil profile.

Table 3 Mass balance of ptaquiloside at the pedon-scale, showing the median annual mass (MED), standard deviations (SD) and percentage of each process or pool of the total. Values are calculated as total mass of PTA divided by the number of years (21 years, n = 100)

	Sandy			Loamy
	MED (mg m^−2)	SD (mg m^−2)	%	MED (mg m^−2)	SD (mg m^−2)	%
Total generation	173.1	83.1	100.0	180.3	86.3	100.0
Canopy wash	82.6	48.3	47.7	82.1	48.1	45.5
Litter degradation	68.7	31.7	39.7	75.1	34.6	41.6
Litter wash	12.9	6.0	7.5	14.5	6.7	8.1
Total input soil	92.7	51.8	100.0	93.9	52.3	100.0
Soil degradation	92.6	51.8	99.9	93.9	52.2	99.9
Leach macropores	0.00	0.00	0.00	0.03	0.02	0.03
Leach matrix	0.09	0.22	0.10	0.10	0.12	0.10
Total leaching	0.09	0.22	0.10	0.13	0.15	0.13

---

3.4. Leaching dynamics

The PTA concentrations in the leachate have been calculated as annual flux-averaged predicted environmental concentrations (PECs) for each of the hundred simulations, i.e. the total mass of toxin leached during the entire period (mg m⁻²) divided by the total volume of water that percolated the 100 cm deep layer (L m⁻²) and the number of years (21). Moreover, each year has been divided into two seasons: the growing season and autumn/winter. The growing season comprises the period between the day of emergence and senescence of bracken, while autumn/winter is the remaining period from senescence until day of emergence the next season. The decision of dividing the year into these periods was made to determine if there are systematic differences between the two periods.

Toxin concentrations at 100 cm are generally higher during the winter/autumn period than in the growing season, and in the sandy loam soil compared with the sandy soil (Fig. 4). The median PECs during the growing season period are 4 and 140 ng L⁻¹ for the sandy and sandy loam soil, respectively, while for the autumn/winter period, the median concentrations are 30 and 120 ng L⁻¹ for the two soils. The maximum PEC in the worst case simulations, i.e. a, r and f_WC parameters close to maximum and k* close to minimum, reaches up to 3000 ng L⁻¹ in both soils and periods. Note that these are annual averages, but single events in extreme scenarios may produce PTA concentrations up to 60 μg L⁻¹.

Fig. 4 Annual PTA leaching concentrations for both soils at a 100 cm depth for the entire simulation period (1985–2005), divided between growing season (GS) and autumn/winter period (AUT/WIN). Annual flux-averaged concentrations calculated as the total mass of PTA leached in a year divided by the total percolated water for the same period. The red line inside the box represents the median; the box edges represent the 25^th and 75^th percentiles; and the error bars represent the 10^th and 90^th percentiles (21 years, n = 2100). The red crosses are outliers which fall outside the range defined by the 10^th and 90^th percentiles.

The higher PEC in the leachate in sandy loam soil can be explained by the presence of preferential flow during intense precipitation events. Macropore transport, which only accounts for 1.5% of the total percolating water, is responsible for 23% of the total PTA leaching in the soil profile. The reason is that PTA transported in the macropores suffers no degradation in the soil profile, as it by-passes the biologically active layer of the soil. However in the case of sandy soil, PTA concentrations have a higher variability than in sandy loam soil (Fig. 4). This higher variability is explained by the fact that transport takes place by matrix flow only, leading to a longer residence time compared with sandy loam soil and thus to a higher variation in the degradation.

During the autumn/winter period, the concentrations below the root zone (100 cm) are more homogeneous in all scenarios and soils, indicated by the lower range in concentrations in the boxplot (10^th to 90^th percentile) (Fig. 4). This is due to different precipitation patterns, lower input to the soil exclusively from litter, and low degradation rates because of lower soil temperatures. The peak concentrations (90^th perc.) during autumn/winter are lower than those in the growing season for the sandy loam soil compared with the sandy soil where concentrations are higher during the growing season. This shows different dynamics in both soils, as in the sandy loam soil the toxin present in the soil will be transported preferentially through macropores and will not suffer from degradation once the topsoil has reached water saturation at the start of autumn (Fig. 3d). Based on the results, we identify the period right after senescence as the time with a higher probability of significant PTA loads reaching the biologically inactive layer of the soil, hence leading to higher toxin concentrations in the leachate. This type of dynamics for PTA leaching agrees with observations in previous studies of PTA^7,18,40 or pesticides with similar physicochemical characteristics and for temperate climates.⁷¹ This means that such areas with extensive bracken coverage, for instance England, Ireland and New Zealand, are at risk for PTA contamination of the groundwater. Moreover, the model results indicate that, on average, there is a constant leaching of PTA with base flow concentrations in the order of ng L⁻¹. These base flow PTA concentrations, of the same order of magnitude, have been previously measured in streams running through bracken infested areas in England.²⁰

Exudation of PTA from rhizomes is an unsettled question. If PTA is also found to be released from bracken rhizomes, the PEC will be higher, as the residence time of PTA in the soil is expected to be lower than PTA released from the above-ground bracken biomass. Furthermore, PTA is only one out of three major toxic terpenes produced by bracken, which can be produced in comparable amounts to those found for PTA.¹² Moreover, another factor to consider is that the PTA content in rhizomes is at its maximum during autumn and winter. Our study also highlights the importance of optimum sampling design in studies dealing with PTA leaching. Therefore, high frequency sampling for plants, soil and water may be required to obtain the peak concentrations as the leaching pulses take place over short periods of time.

4. Conclusion

This is the first modeling approach for simulating the fate of natural toxins in the soil–plant–atmosphere continuum. The model reproduces bracken biomass and canopy PTA content dynamics and agrees with previous published experimental studies. PTA can be released from the canopy in great amounts compared with pesticides for instance. Out of the total amount of PTA produced in the canopy, often in the order of kg ha⁻¹, 48% is released from the plant by precipitation and input to the soil in loads in the order of mg m⁻². The perfect conditions for an extreme toxin wash off event are when intense precipitation takes place in July or August, with a fully developed bracken canopy and maximum PTA content in the canopy.

In the soil, degradation dominates the fate of PTA, eliminating 99.9% of the total input to the soil. The time with the expected highest PEC in the leachate is during the autumn/winter period, where the median PECs for both soils were in the range of ng L⁻¹ and peak concentrations were in the order of μg L⁻¹ and hence several orders of magnitude above the maximum tolerable concentrations.

The new DAISY model serves as a valuable tool for predicting the fate of natural toxins in the field. The flexibility of the plant and toxin generation module offers the possibility of calibrating the model for different cases of plants and toxins of interest. A formal comparison of PTA monitoring data with model results is needed for the calibration of the model. We identify that further research should concentrate on the production and wash off coefficient of PTA at different development stages, as well as on the identification of other possible sources of PTA not considered in the model, such as root exudation.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors are grateful to Merete Styczen and Lars H. Rasmussen for their useful comments and advice to make this study possible. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 722493 (NaToxAq).

References

1. Inderjit, R. M. Callaway and J. M. Vivanco, Can plant biochemistry contribute to understanding of invasion ecology?, Trends Plant Sci., 2006, 11, 574–580,  DOI:10.1016/j.tplants.2006.10.004.
2. I. T. Baldwin, Finally, Proof of Weapons of Mass Destruction, Sci. Signaling, 2003, 2003, pe42,  DOI:10.1126/stke.2003.203.pe42.
3. C. Page, The taxonomy and phytogeography of bracken—a review, Bot. J. Linn. Soc., 1976, 73, 1–34.
4. A. Mithöfer and M. E. Maffei, in Plant Toxins, ed. P. Gopalakrishnakone, C. Carlini and R. Ligabue-Braun, Springer, Dordrecht, 2017, ch. 1, pp. 3–24,  DOI:10.1007/978-94-007-6464-4_21.
5. T. Isah, Stress and defense responses in plant secondary metabolites production, Biol. Res., 2019, 52, 39,  DOI:10.1186/s40659-019-0246-3.
6. K. K. Jessing, N. Cedergreen, J. Jensen and H. C. B. Hansen, Degradation and ecotoxicity of the biomedical drug artemisinin in soil, Environ. Toxicol. Chem., 2009, 28, 701–710,  DOI:10.1897/08-153r.1.
7. C. Ramwell, W. van Beinum, A. Rowbotham, H. Parry, S. Parsons, W. Luo and G. Evans, Ptaquiloside & other bracken toxins: A preliminary risk assessment, Food and Environment Research Agency, Sand Hutton, York, UK, 2010.
8. K. Refsgaard, N. Bjarnholt, B. L. Møller, M. M. Saddik and H. C. Bruun Hansen, Dissipation of cyanogenic glucosides and cyanide in soil amended with white clover (Trifolium repens L.), Soil Biol. Biochem., 2010, 42, 1108–1113,  DOI:10.1016/j.soilbio.2010.03.008.
9. C. C. Hoerger, F. E. Wettstein, H. J. Bachmann, K. Hungerbuhler and T. D. Bucheli, Occurrence and mass balance of isoflavones on an experimental grassland field, Environ. Sci. Technol., 2011, 45, 6752–6760,  DOI:10.1021/es200567b.
10. J. R. Hama and B. W. Strobel, Pyrrolizidine alkaloids quantified in soil and water using UPLC-MS/MS, RSC Adv., 2019, 9, 30350–30357,  10.1039/c9ra05301h.
11. L. H. Rasmussen, H. C. Hansen and D. Lauren, Sorption, degradation and mobility of ptaquiloside, a carcinogenic bracken (Pteridium sp.) constituent, in the soil environment, Chemosphere, 2005, 58, 823–835,  DOI:10.1016/j.chemosphere.2004.08.088.
12. V. Kisielius, D. N. Lindqvist, M. B. Thygesen, M. Rodamer, H. C. B. Hansen and L. H. Rasmussen, Fast LC-MS quantification of ptesculentoside, caudatoside, ptaquiloside and corresponding pterosins in bracken ferns, J. Chromatogr. B: Anal. Technol. Biomed. Life Sci., 2020, 121966,  DOI:10.1016/j.jchromb.2019.121966.
13. L. J. Shaw and J. E. Hooker, The fate and toxicity of the flavonoids naringenin and formononetin in soil, Soil Biol. Biochem., 2008, 40, 528–536,  DOI:10.1016/j.soilbio.2007.09.021.
14. D. W. Kolpin, C. C. Hoerger, M. T. Meyer, F. E. Wettstein, L. E. Hubbard and T. D. Bucheli, Phytoestrogens and Mycotoxins in Iowa Streams: An Examination of Underinvestigated Compounds in Agricultural Basins, J. Environ. Qual., 2010, 39, 2089,  DOI:10.2134/jeq2010.0121.
15. C. C. Hoerger, J. Schenzel, B. W. Strobel and T. D. Bucheli, Analysis of selected phytotoxins and mycotoxins in environmental samples, Anal. Bioanal. Chem., 2009, 395, 1261–1289,  DOI:10.1007/s00216-009-3088-y.
16. D. Selmar, C. Wittke, I. Beck-von Wolffersdorff, B. Klier, L. Lewerenz, M. Kleinwachter and M. Nowak, Transfer of pyrrolizidine alkaloids between living plants: A disregarded source of contaminations, Environ. Pollut., 2019, 248, 456–461,  DOI:10.1016/j.envpol.2019.02.026.
17. L. H. Rasmussen, S. Kroghsbo, J. C. Frisvad and H. C. B. Hansen, Occurrence of the carcinogenic bracken constituent ptaquiloside in fronds, topsoils and organic soil layers in Denmark, Chemosphere, 2003, 51(2), 117–127,  DOI:10.1016/s0045-6535(02)00694-x.
18. F. Clauson-Kaas, P. H. Jensen, O. S. Jacobsen, R. K. Juhler and H. C. B. Hansen, The naturally occurring carcinogen ptaquiloside is present in groundwater below bracken vegetation, Environ. Toxicol. Chem., 2014, 33, 1030–1034,  DOI:10.1002/etc.2533.
19. K.-K. Tung, C. K. Chan, Y. Zhao, K. K. J. Chan, G. Liu, N. M. Pavlovic and W. Chan, Occurrence and Environmental Stability of Aristolochic Acids in Groundwater Collected from Serbia: Links to Human Exposure and Balkan Endemic Nephropathy, Environ. Sci. Technol., 2019, 54, 1554–1561,  DOI:10.1021/acs.est.9b05337.
20. F. Clauson-Kaas, C. Ramwell, H. C. B. Hansen and B. W. Strobel, Ptaquiloside from bracken in stream water at base flow and during storm events, Water Res., 2016, 106, 155–162,  DOI:10.1016/j.watres.2016.09.049.
21. D. Selmar, A. Radwan, T. Hijazin, S. Abouzeid, M. Yahyazadeh, L. Lewerenz, M. Kleinwächter and M. Nowak, Horizontal Natural Product Transfer: Intriguing Insights into a Newly Discovered Phenomenon, J. Agric. Food Chem., 2019, 67, 8740–8745,  DOI:10.1021/acs.jafc.9b03619.
22. P. W. Hyder, E. L. Fredrickson, R. E. Estell and M. E. Lucero, Transport of phenolic compounds from leaf surface of creosotebush and tarbush to soil surface by precipitation, J. Chem. Ecol., 2002, 28, 2475–2482,  DOI:10.1023/a:1021432018512.
23. K. K. Sjøholm, B. W. Strobel and N. Cedergreen, Artemisia annua, CRC Press, 2017, pp. 131–154.
24. R. H. Marrs and A. S. Watt, Biological Flora of the British Isles: Pteridium aquilinum (L.) Kuhn, J. Ecol., 2006, 94, 1272–1321,  DOI:10.1111/j.1365-2745.2006.01177.x.
25. K. Yamada, M. Ojika and H. Kigoshi, Ptaquiloside, the major toxin of bracken, and related terpene glycosides: chemistry, biology and ecology, Nat. Prod. Rep., 2007, 24, 798–813,  10.1039/b614160a.
26. A. S. Prakash, T. N. Pereira, B. L. Smith, G. Shaw and A. A. Seawright, Mechanism of Bracken Fern Carcinogenesis: Evidence for H-Ras Activation Via Initial Adenine Alkylation by Ptaquiloside, J. Urol., 1998, 159, 600,  DOI:10.1016/s0022-5347(01)63876-x.
27. R. M. Gil Da Costa, M. M. S. M. Bastos, P. A. Oliveira and C. Lopes, Bracken-associated human and animal health hazards: Chemical, biological and pathological evidence, J. Hazard. Mater., 2012, 203, 1–12,  DOI:10.1016/j.jhazmat.2011.12.046.
28. D. M. Potter and M. S. Baird, Carcinogenic effects of ptaquiloside in bracken fern and related compounds, Br. J. Cancer, 2000, 83, 914–920,  DOI:10.1054/bjoc.2000.1368.
29. P. C. D. R. Aranha, L. H. Rasmussen, G. A. Wolf-Jäckel, H. M. E. Jensen, H. C. B. Hansen and C. Friis, Fate of ptaquiloside—A bracken fern toxin—In cattle, PLoS One, 2019, 14, e0218628,  DOI:10.1371/journal.pone.0218628.
30. P. J. O'Connor, M. E. Alonso-Amelot, S. A. Roberts and A. C. Povey, The role of bracken fern illudanes in bracken fern-induced toxicities, Mutat. Res., Rev. Mutat. Res., 2019, 782, 108276,  DOI:10.1016/j.mrrev.2019.05.001.
31. H. Mori, S. Sugie, I. Hirono, K. Yamada, H. Niwa and M. Ojika, Genotoxicity of ptaquiloside, a bracken carcinogen, in the hepatocyte primary culture/DNA-repair test, Mutat. Res. Lett., 1985, 143, 75–78,  DOI:10.1016/0165-7992(85)90108-3.
32. P. Engel, K. K. Brandt, L. H. Rasmussen, R. G. Ovesen and J. Sørensen, Microbial degradation and impact of bracken toxin ptaquiloside on microbial communities in soil, Chemosphere, 2007, 67, 202–209,  DOI:10.1016/j.chemosphere.2006.08.025.
33. R. G. Ovesen, L. H. Rasmussen and H. C. B. Hansen, Degradation kinetics of ptaquiloside in soil and soil solution, Environ. Toxicol. Chem., 2008, 27, 252–259,  DOI:10.1897/07-324r.1.
34. K. B. Ayala-Luis, P. B. Hansen, L. H. Rasmussen and H. C. B. Hansen, Kinetics of ptaquiloside hydrolysis in aqueous solution, Environ. Toxicol. Chem., 2006, 25, 2623–2629,  DOI:10.1897/05-695r.1.
35. L. Rasmussen, D. Lauren, B. Smith and H. Hansen, Variation in ptaquiloside content in bracken (Pteridium esculentum (Forst. f) Cockayne) in New Zealand, N. Z. Vet. J., 2008, 56, 304–309,  DOI:10.1080/00480169.2008.36851.
36. C. Zaccone, I. Cavoski, R. Costi, G. Sarais, P. Caboni, A. Traversa and T. M. Miano, Ptaquiloside in Pteridium aquilinum subsp. aquilinum and corresponding soils from the South of Italy: influence of physical and chemical features of soils on its occurrence, Sci. Total Environ., 2014, 496, 365–372,  DOI:10.1016/j.scitotenv.2014.07.046.
37. A. M. Rodriguez, M. G. Derita, S. A. Borkosky, C. Socolsky, A. Bardón and M. A. Hernández, Bioactive farina of Notholaena sulphurea (Pteridaceae): Morphology and histochemistry of glandular trichomes, Flora, 2018, 240, 144–151,  DOI:10.1016/j.flora.2018.01.008.
38. T. P. Atkinson, Seasonal and altitudinal variation in Pteridium aquilinum (L.) Kuhn: frond and stand types, New Phytol., 1989, 113, 359–365,  DOI:10.1111/j.1469-8137.1989.tb02414.x.
39. L. H. Rasmussen, E. Donnelly, B. W. Strobel, P. E. Holm and H. C. B. Hansen, Land management of bracken needs to account for bracken carcinogens – A case study from Britain, J. Environ. Manage., 2015, 151, 258–266,  DOI:10.1016/j.jenvman.2014.12.052.
40. C. O'Driscoll, C. Ramwell, B. Harhen, L. Morrison, F. Clauson-Kaas, H. C. Hansen, G. Campbell, J. Sheahan, B. Misstear and L. Xiao, Ptaquiloside in Irish Bracken Ferns and Receiving Waters, with Implications for Land Managers, Molecules, 2016, 21, 543,  DOI:10.3390/molecules21050543.
41. M. T. Fletcher, K. G. Reichmann, I. J. Brock, R. A. McKenzie and B. J. Blaney, Residue Potential of Norsesquiterpene Glycosides in Tissues of Cattle Fed Austral Bracken (Pteridium esculentum), J. Agric. Food Chem., 2011, 59, 8518–8523,  DOI:10.1021/jf201342t.
42. A. Virgilio, A. Sinisi, V. Russo, S. Gerardo, A. Santoro, A. Galeone, O. Taglialatela-Scafati and F. Roperto, Ptaquiloside, the major carcinogen of bracken fern, in the pooled raw milk of healthy sheep and goats: an underestimated, global concern of food safety, J. Agric. Food Chem., 2015, 63, 4886–4892,  DOI:10.1021/acs.jafc.5b01937.
43. L. H. Rasmussen, L. S. Jensen and H. C. B. Hansen, Distribution of the carcinogenic terpene ptaquiloside in bracken fronds, rhizomes (Pteridium aquilinum), and litter in Denmark, J. Chem. Ecol., 2003, 29, 771–778,  DOI:10.1023/a:1022885006742.
44. B. M. Lange and A. Ahkami, Metabolic engineering of plant monoterpenes, sesquiterpenes and diterpenes-current status and future opportunities, Plant Biotech. J., 2013, 11, 169–196,  DOI:10.1111/pbi.12022.
45. A. Huchelmann, M. Boutry and C. Hachez, Plant Glandular Trichomes: Natural Cell Factories of High Biotechnological Interest, Plant Physiol., 2017, 175, 6–22,  DOI:10.1104/pp.17.00727.
46. E. Wollenweber and H. Schneider, Lipophilic exudates of Pteridaceae – chemistry and chemotaxonomy, Biochem. Syst. Ecol., 2000, 28, 751–777,  DOI:10.1016/s0305-1978(99)00118-0.
47. N. Glæsner, E. Diamantopoulos, J. Magid, C. Kjaergaard and H. H. Gerke, Modeling solute mass exchange between pore regions in slurry-injected soil columns during intermittent irrigation, Vadose Zone J., 2018, 17, 1–10,  DOI:10.2136/vzj2018.01.0006.
48. A. Boivin, J. Šimůnek, M. Schiavon and M. T. Van Genuchten, Comparison of Pesticide Transport Processes in Three Tile-Drained Field Soils Using HYDRUS-2D, Vadose Zone J., 2006, 5, 838–849,  DOI:10.2136/vzj2005.0089.
49. G. Alavi, J. Dusek, T. Vogel, R. Green and C. Ray, Evaluation of Dual-Permeability Models for Chemical Leaching Assessment to Assist Pesticide Regulation in Hawaii, Vadose Zone J., 2007, 6, 735–745.
50. J. J. T. I. Boesten, Effects of aged sorption on pesticide leaching to groundwater simulated with PEARL, Sci. Total Environ., 2017, 576, 498–507,  DOI:10.1016/j.scitotenv.2016.10.099.
51. H. Labite, N. M. Holden, K. G. Richards, G. Kramers, A. Premrov, C. E. Coxon and E. Cummins, Comparison of pesticide leaching potential to groundwater under EU FOCUS and site specific conditions, Sci. Total Environ., 2013, 463, 432–441,  DOI:10.1016/j.scitotenv.2013.06.050.
52. S. B. Rasmussen, P. Abrahamsen, M. H. Nielsen, P. E. Holm and S. Hansen, Effects of single rainfall events on leaching of glyphosate and bentazone on two different soil types, using the DAISY model, Vadose Zone J., 2015, 14(11), 1–15,  DOI:10.2136/vzj2014.11.0164.
53. P. Abrahamsen and S. Hansen, Daisy: an open soil–crop–atmosphere system model, Environ. Model. Software, 2000, 15, 313–330,  DOI:10.1016/S1364-8152(00)00003-7.
54. S. Hansen, Daisy, a flexible soil–plant–atmosphere system model, Report. Dept. Agric., 2002.
55. K. Manevski, C. D. Børgesen, M. N. Andersen and I. S. Kristensen, Reduced nitrogen leaching by intercropping maize with red fescue on sandy soils in North Europe: a combined field and modeling study, Plant Soil, 2015, 388, 67–85,  DOI:10.1007/s11104-014-2311-6.
56. C. D. Børgesen and J. E. Olesen, A probabilistic assessment of climate change impacts on yield and nitrogen leaching from winter wheat in Denmark, Nat. Hazards Earth Syst. Sci., 2011, 11, 2541.
57. S. Hansen, C. T. Petersen, M. Mollerup, P. Abrahamsen, B. Gjettermann, M. H. Nielsen, M. E. Styczen, R. Poulsen, J. K. Lørup and K. Yamagata, Flerdimensional modellering af vandstrømning og stoftransport i de øverste 1–2 m af jorden i systemer med markdræn, 2012.
58. E. Diamantopoulos, S. C. Iden and W. Durner, Inverse modeling of dynamic nonequilibrium in water flow with an effective approach, Water Resour. Res., 2012, 48(3) DOI:10.1029/2011wr010717.
59. D. Hillel, Introduction to soil physics, Academic Press, 2013.
60. S. Hansen, H. Jensen, N. Nielsen and H. Svendsen, DAISY: Soil Plant Atmosphere System Model, NPO Report No. A 10, The National Agency for Environmental Protection, 1990, pp. 272–281.
61. S. Bruun, B. T. Christensen, E. M. Hansen, J. Magid and L. S. Jensen, Calibration and validation of the soil organic matter dynamics of the Daisy model with data from the Askov long-term experiments, Soil Biol. Biochem., 2003, 35, 67–76,  DOI:10.1016/s0038-0717(02)00237-7.
62. L. Rasmussen and H. Hansen, Growth of Bracken in Denmark and the Content of Ptaquiloside in Fronds, Poisonous Plants and Related Toxins, 2004, 354.
63. J. I. Pitman, Rainfall interception by bracken in open habitats—Relations between leaf area, canopy storage and drainage rate, J. Hydrol., 1989, 105(3–4), 317–334,  DOI:10.1016/0022-1694(89)90111-x.
64. R. J. Pakeman, R. H. Marrs and P. J. Jacob, A Model of Bracken (Pteridium aquilinum) Growth and the Effects of Control Strategies and Changing Climate, J. Appl. Ecol., 1994, 31, 145,  DOI:10.2307/2404607.
65. S. Hansen, P. Abrahamsen, C. T. Petersen and M. Styczen, Daisy: Model use, calibration, and validation, Trans. ASABE, 2012, 55, 1317–1333.
66. FOCUS, FOCUS groundwater scenarios in the EU review of active substances, Report of the FOCUS Groundwater Scenarios Workgroup, EC Document Reference Sanco/321/2000 rev.2, 2000, p. 202.
67. FOCUS, Generic guidance for FOCUS surface water scenarios. Version 1.4, 2015.
68. M. Styczen, C. T. Petersen, C. B. Koch and B. Gjettermann, Macroscopic Evidence of Sources of Particles for Facilitated Transport during Intensive Rain, Vadose Zone J., 2011, 10, 1151–1161,  DOI:10.2136/vzj2010.0124.
69. A. S. Watt, Contributions to the Ecology of Bracken (Pteridium aquilinum) VII. Bracken and Litter I. The Origin of Rings, New Phytol., 1956, 55, 369–381.
70. M. Shehadul Islam, A. Aryasomayajula and P. Selvaganapathy, A Review on Macroscale and Microscale Cell Lysis Methods, Micromachines, 2017, 8, 83,  DOI:10.3390/mi8030083.
71. A. E. Rosenbom, P. Olsen, F. Plauborg, R. Grant, R. K. Juhler, W. Brüsch and J. Kjær, Pesticide leaching through sandy and loamy fields – Long-term lessons learnt from the Danish Pesticide Leaching Assessment Programme, Environ. Pollut., 2015, 201, 75–90,  DOI:10.1016/j.envpol.2015.03.002.

---