COMSOL 6.4 - Pesticide Transport and Reaction in Soil

Pesticide Transport and Reaction in Soil

Aldicarb is a commercial pesticide, used on a variety of crops, including cotton, sugar beet, citrus fruits, potatoes, and beans. The general population may be exposed to aldicarb primarily through the ingestion of contaminated water and foods.

This example looks at the degradation kinetics of aldicarb and its toxic by-products, investigating both the degradation time-scale as well as the spatial concentration distribution of toxic components. In the first model the chemicals are contained in a water pond, treated as a perfectly mixed system. The second model tracks the detailed distribution of chemicals in soil as the pesticide leaches out of the pond and is transported in water through the ground.

Model Definition

Aldicarb degrades by transformation to the corresponding sulfoxide and the sulfone (both of which are toxic), and is detoxified by hydrolysis to oximes and nitriles. The chain of reactions is illustrated in Figure 1. The toxicity of a chemical species is indicated by its LD50 value, signifying the median lethal dose (mg/kg) to half of a test population of rats. As indicated, both the sulfoxide and sulfone analogues of aldicarb are also relatively toxic.

Figure 1: Reaction pathways of aldicarb degradation.

Each of the j unimolecular reactions outlined above has a rate expression of the form

Note that in this example the concentration unit is mol/m3 and the rate constants are expressed in 1/day.

Perfectly Mixed System

The first model solves for the decomposition kinetics of aldicarb occurring in a water pond. The pond is treated as a closed and perfectly mixed system. The reaction mechanism illustrated in Figure 1 translates into the following mass balance equations:

•

For aldicarb (a)

•

For aldicarb sulfoxide (asx)

•

For aldicarb sulfone (asn)

•

For aldicarb oxime (ao)

•

For aldicarb sulfoxide oxime (asxo)

•

For aldicarb sulfone oxime (asno)

Solving this set of coupled ODEs provides information on the time scales of the degradation processes.

Space- and Time-Dependent System

In a more detailed model, assume that aldicarb moves from the pond into relatively dry soil. In the soil, the aldicarb decomposes according to the mechanism illustrated in Figure 1. In addition, the pesticide and its decay products are transported by convection, dispersion, sorption, and volatilization.

Geometry

Water is ponded by a ring sitting on the ground. The soil is layered and rests on rocks, with the top layer slightly less permeable than the bottom one. Water moves through the bottom of the ring into the soil. The water level in the ring is known, as is the initial distribution of pressure heads in the soil. There is no flow through the vertical walls or the surface outside of the ring.

Figure 2: Geometry of the infiltration ring and soil column.

Aldicarb moves with water from the pond into the soil at a constant concentration. In the soil, the chemicals react and also adsorb onto soil particles. Aldicarb and the aldicarb sulfone volatilize to the atmosphere. The sorption, biodegradation, and volatilization proceed in linear proportion to the aqueous concentrations. The soil is initially pristine with zero concentration of the involved chemicals. At the ground surface outside the ring there is volatilization to the atmosphere for ca and casn. Model the problem with 2D axisymmetry and track the solute transport for 10 days. The vertical axis to the left is the symmetry axis. Add infinite elements to the right axis such that the solutes can freely leave the soil column with the fluid flow.

Fluid Flow

Richards’ equation governs the saturated-unsaturated flow of water in the soil. The soil pores are connected to the atmosphere, so you can assume that pressure changes in the air do not affect the flow and use Richards’ equation for single-phase flow. Given by Ref. 1, Richards’ equation in pressure head reads

where Cm denotes specific moisture capacity (m−1); Se is the effective saturation of the soil (dimensionless); Sp is a storage coefficient (m−1); Hp is the pressure head (m), which is proportional to the dependent variable, p (Pa); t is time; K equals the hydraulic conductivity (m/s); D is the direction (typically, the z direction) that represents vertical elevation (m).

To be able to combine boundary conditions and sources with the Darcy’s Law formulation, COMSOL Multiphysics converts Richards’ equation to SI units and solves for the pressure (SI unit: Pa). Hydraulic head, H, pressure head, Hp, and elevation D are related to pressure p as

Also, the permeability κ (SI unit: 1/m2) and hydraulic conductivity K (SI unit: m/s) are related to the viscosity μ (SI unit: Pa·s) and density ρ (SI unit: kg/m3) of the fluid, and the acceleration of gravity g (SI unit: m/s2) by

In this problem, Sp = (θs − θr)/(1 m·ρg) where θs and θr denote the volume fraction of fluid at saturation and after drainage, respectively. For more details see The Richards’ Equation Interface in the Subsurface Flow Module User’s Guide.

Mass Transport

The governing equation for solute transport describes advection and dispersion of a sorbing, volatilizing, and decaying solute in variably saturated soil.

(1)

The Transport of Diluted Species in Porous Media interface implements Equation 1. It describes the time rate of change in two terms: c denotes dissolved concentration (mol/m3) and cP is the mass of adsorbed contaminant per dry unit weight of solid (mg/kg). Further, θ denotes the volume fraction of fluid (porosity), and ρb is the bulk density (kg/m3). Because ρb amounts to the dried solid mass per bulk volume of the solids and pores together, the term ρbcP gives solute mass attached to the soil as the concentration changes with time.

Solute spreading now includes mechanical dispersion in water plus molecular diffusion for water and air. These three processes appear in the liquid-gas dispersion tensor, whose entries are

In these equations, DLGii are the principal components of the liquid-gas dispersion tensor; DLGij and DLGji are the cross terms; α is the dispersivity (m) where the subscripts “1” and “2” denote longitudinal and transverse dispersivities, respectively; Dm and DG (m2/d) are molecular diffusion coefficients; and τL and τG give the tortuosity factors for liquid (water) and gas (air), respectively.

The three solutes — aldicarb, aldicarb sulfoxide, and aldicarb sulfone — have different decay terms, RLi, partition coefficients, kPi, and volatilization constants, kGi. All of the solutes adsorb to soil particles, but only two of the solutes volatilize; sulfoxide does not.

Model Data

The following table provides data for the fluid-flow model:


Variable	Unit	Description	upper layer	Lower layer
Ks	m/d	Saturated hydraulic conductivity	0.298	0.454
θs		Porosity/void fraction	0.399	0.339
θr		Residual saturation	0.001	0.001
α	m-1	alpha parameter	1.74	1.39
n		n parameter	1.38	1.60
Hp0	m	Pressure head in ring	0.01
Hp,init	m	Initial pressure head	-(z+1.2) -0.2(z+0.4)	-(z+1.2)

The inputs needed for the solute-transport model are:


Variable	Units	Description	Value
ρb	kg/m3	Bulk density	1400
kp	m3/kg	Partition coefficient	0.0001
Dm	m2/d	Coefficient of molecular diffusion	0.00374
αr	m	Longitudinal dispersivity	0.005
αz	m	Transverse dispersivity	0.001
	d-1	Decay rate in liquid	0.05
	d-1	Decay rate on solid	0.01
c0	mol/m3	Solute concentration in ring	1.0

Results and Discussion

First, review the results of the perfectly mixed reactor model. Figure 3 shows the concentration profiles of aldicarb and all of its decay products as well as the sum of the three most toxic species — aldicarb, aldicarb sulfoxide, and aldicarb sulfone (see Figure 1 for LD50 values). Only small amounts of aldicarb remain after 10 days. Considering the summed-up contributions, contamination levels clearly remain high even after several months.

Figure 3: Concentration profiles as reactions occur during a 100 day time period.

The following results come from the space- and time-dependent model setup. Figure 4 shows the Pressure Head distribution after 10 days

Figure 4: Pressure head after 10 days.

Figure 5 shows the fluid flow in soil after 0.3 days (left) and 1.0 day (right). The plots illustrate the wetting of the soil with time. As indicated by the arrows, the fluid velocities are relatively high beneath the ponded water.

Figure 5: The effective saturation (surface plot), and flow velocity (arrows) in a variably saturated soil after 0.3 days (left) and 1 day (right).

Figure 6 shows the concentration distributions of aldicarb and the equally toxic aldicarb sulfoxide after 1, 5, and 10 days of infiltration. Consistent with the evolving flow field, the main direction of transport is in the vertical direction.

Figure 6: Concentration of aldicarb (left) and aldicarb sulfoxide (right) after 1, 5, and 10 days (note the differing color ranges).

The distribution of aldicarb has clearly reached steady-state conditions after 10 days, a time frame that was also predicted by the ideal reactor model (see Figure 3). Results also show that the soil contamination is rather local with respect to the aldicarb source. The aldicarb sulfoxide, on the other hand, can be expected to affect a considerably larger soil volume for a significantly longer time.

Notes About the COMSOL Implementation

This model makes use of the feature. It performs a coordinate scaling to the selected domain such that boundary conditions on the outside of the infinite element layer are effectively applied at a very large distance. Therefore unwanted effects of artificial boundary conditions on the region of interest are suppressed. This allows to model details in an area which is actually very large or infinite.

References

1. J. Bear, Hydraulics of Groundwater, McGraw Hill, 1978.

2. M.Th. van Genuchten, “A closed-form equation for predicting the hydraulic of conductivity of unsaturated soils,” Soil Sci. Soc. Am. J., vol. 44, pp. 892–898, 1980.

Subsurface_Flow_Module/Solute_Transport/pesticide_transport

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Load the rate constants from file.

Rate constants

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pesticide_transport_parameters_1.txt.

In the text field, type Rate constants.

First consider the aldicarb decomposition kinetics in the water pond treated as a perfectly mixed system.

Global ODEs and DAEs (ge)

Global Equations 1 (ODE1)

Read in a set of the equations defining the reactions.

In the window, under > click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pesticide_transport_reactions.txt.

Locate the section. Click

In the dialog, type concentration in the text field.

In the tree, select > .

Click .

In the window for , locate the section.

Click

In the dialog, type reactionrate in the text field.

In the tree, select > .

Click .

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the text field, type range(0,1,100).

In the toolbar, click

Results

Concentration of species (100 days)

A plot showing the time dependent concentration is added automatically. Edit this plot to reproduce Figure 3 with the following steps:

In the window for , type Concentration of species (100 days) in the text field.

Click to expand the section. From the list, choose .

Locate the section.

Select the checkbox. In the associated text field, type Concentration (mol/m<sup>3</sup>).

Global 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

In the text field, type t.

From the list, choose .

The concentration for all species is added automatically. Add the expression for the sum of the three most toxic species to the table.

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
cpm_a+cpm_asn+cpm_asx	mol/m^3	Sum of most toxic species

In the toolbar, click

Now solve the time- and space-dependent transport and reaction problem in the soil.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Find the subsection. In the table, clear the checkbox for .

Click the button in the window toolbar.

In the tree, select > .

Click to expand the section. In the text field, type 3.

In the table, enter the following settings:


c_a
c_asx
c_asn

In the table, clear the checkbox for .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the table, clear the checkbox for .

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

Step 1: Time Dependent

Load the parameters defining the material properties and the geometry from file.

Global Definitions

Transport parameters

In the toolbar, click

and choose > .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pesticide_transport_parameters_2.txt.

In the text field, type Transport parameters.

Geometry 1

The modeling domain is made up of the two permeable soil layers, each of which is represented by a rectangular domain in 2D axisymmetry.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 1.5.

In the text field, type 0.9.

Locate the section. In the text field, type -1.3.

Click to expand the section. Select the checkbox.

Clear the checkbox.

In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.25

This additional layer to the right will later be used to define an . Read more about it in the Notes About the COMSOL Implementation section. Proceed with the second soil layer.

Rectangle 2 (r2)

Right-click and choose .

In the window for , locate the section.

In the text field, type 0.4.

Locate the section. In the text field, type -0.4.

To finish the model geometry, add a point on the top boundary marking the pond’s outer rim.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.25.

Click

Click the

button in the toolbar.

Now, define the .

Definitions

Infinite Element Domain 1 (ie1)

In the toolbar, click

Select Domains 3 and 4 only.

In the window for , locate the section.

From the list, choose .

Variables 1

In the toolbar, click

Load the rate expressions from file.

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pesticide_transport_variables.txt.

Richards’ Equation (dl)

Begin by specifying the properties for the bottom soil layer in the default node, then duplicate this node and modify the domain selection and properties to match the top layer.

Unsaturated Porous Medium 1

In the window, under > click .

In the window for , locate the section.

From the list, choose . In the Sp text field, type Sp_1.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the Ks text field, type Ks_1.

Locate the section. In the α text field, type alpha_1.

In the n text field, type n_1.

In the θr text field, type thetar_1.

Unsaturated Porous Medium 2

In the window, under > right-click and choose .

Select Domains 2 and 4 only.

In the window for , locate the section.

In the Sp text field, type Sp_2.

Porous Matrix 1

In the window, expand the node, then click .

In the window for , locate the section.

In the Ks text field, type Ks_2.

Locate the section. In the α text field, type alpha_2.

In the n text field, type n_2.

In the θr text field, type thetar_2.

Gravity 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

Click the button.

In the Hp text field, type -(z+1.2).

Initial Values 2

Right-click > > and choose .

Select Domains 2 and 4 only.

In the window for , locate the section.

In the Hp text field, type -(z+1.2)-0.2*(z+0.4).

Pressure Head 1

In the toolbar, click

and choose .

Select Boundary 5 only.

In the window for , locate the section.

In the Hp0 text field, type Hp0.

Pervious Layer 1

In the toolbar, click

and choose .

Select Boundaries 2 and 8 only.

In the window for , locate the section.

In the Hb text field, type -2.

In the Rb text field, type 1/5[d].

Transport of Diluted Species in Porous Media (tds)

Unsaturated Porous Medium 1

Now, set up the transport equation for an unsaturated porous medium, accounting for adsorption, dispersion and volatilization.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Liquid 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the θl text field, type dl.theta_l.

This corresponds to the liquid volume fraction.

From the list, choose .

From the dHp/dt list, choose .

In the Cm text field, type dl.Cm.

Locate the section. From the u list, choose .

The diffusion coefficients are the same for all species.

Locate the section. In the DL,ca text field, type Dl. Also type Dl in the text fields for the liquid diffusion coefficients of aldicarb sulfoxide and aldicarb sulfone.

Gas 1

Again, the diffusion coefficients are the same for all species.

In the window, click .

In the window for , locate the section.

In the DG,ca text field, type Dg. Also type Dg in the text fields for the gas diffusion coefficients of aldicarb sulfoxide and aldicarb sulfone.

Locate the section. In the kG,ca text field, type kg_a.

In the kG,casn text field, type kg_asn.

Unsaturated Porous Medium 1

In the window, click .

Adsorption 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select the checkbox.

In the cP,ca text field, type kp_a*c_a.

Select the checkbox.

In the cP,casx text field, type kp_asx*c_asx.

Select the checkbox.

In the cP,casn text field, type kp_asn*c_asn.

Unsaturated Porous Medium 1

In the window, click .

Dispersion 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Specify the α matrix as


alphar
alphaz

Reactions 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Rca text field, type dl.theta_l*(-r_1-r_3).

In the Rcasx text field, type dl.theta_l*(r_1-r_2-r_4).

In the Rcasn text field, type dl.theta_l*(r_2-r_5).

Outflow 1

In the toolbar, click

and choose .

Select Boundaries 2, 8, 12, and 13 only. This corresponds to the right and bottom boundaries.

Concentration 1

In the toolbar, click

and choose .

Select Boundary 5 only.

In the window for , locate the section.

Select the checkbox.

In the c0,ca text field, type c0.

Select the checkbox.

Volatilization 1

In the toolbar, click

and choose .

Select Boundaries 6 and 11 only.

In the window for , locate the section.

In the hc text field, type Dg/d_s.

Select the checkbox.

Materials

Some required material properties have not yet been defined. This is indicated by a small red cross at the material node. Continue as follows to add the missing properties.

Porous Material: Lower Layer

In the window, under right-click and choose > .

In the window for , type Porous Material: Lower Layer in the text field.

Locate the section. Click

Fluid 1 (pmat1.fluid1)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rhof	kg/m³	Basic

Solid 1 (pmat1.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 1-poro_1.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rhob	kg/m³	Basic

Porous Material: Upper Layer

In the window, under > right-click and choose .

In the window for , type Porous Material: Upper Layer in the text field.

Select Domains 2 and 4 only.

Solid 1 (pmat2.solid1)

In the window, expand the > > node, then click .

In the window for , locate the section.

In the θs text field, type 1-poro_2.

Mesh 1

Using a mapped mesh is a good idea for this geometry. It uses less mesh elements while keeping the accuracy compared to using a triangular mesh with the same mesh size.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Mapped 1

In the toolbar, click

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 0.02.

In the window, right-click and choose .

Study 2

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the text field, type range(0,0.1,0.9) range(1,1,10).

In the window, click .

In the window for , locate the section.

Clear the checkbox, to disable the automatic generation of default plots. Instead, choose the relevant plots from the Result Templates after computation.

In the toolbar, click