**4. Water and solute transport using HYDRUS software: case study in Croatia**

The example will show the procedure of modelling water flow and nitrogen species transport in soil vadose zone. In this example, urea and NPK fertilizers were added to soil profile. The study site is located in the eastern Croatia in the intensive agricultural production area. The more information about the study can be found in [45]. The project main goal was to evaluate the influence of high fertilizer load, mostly nitrogen based, on the groundwater resources. Here, the results of modelling study on one selected site will be presented from year 2014. HYDRUS‐1D [37] was used to simulate water flow and nitrate transport and its ability to reproduce observed water and nitrate outflows (collected by lysimeters) was assessed. Simulations are carried out in one‐dimensional domain using measured soil hydraulic parameters, climatic daily data, crop growth parameters, and fertilizer application rates as the input for model. Evapotranspiration rates were calculated using CROPWAT model [46] (based on the Penman‐Monteith approach). The soil type at the site was classified as Haplic Gleysol Calcaric Eutric Siltic (Horizons: Ap‐Bg‐Cr‐Cg), from which saturated hydraulic conductivity, retention curve (at different pressure head values), and saturated water content were meas‐ ured. Optimization of necessary remaining hydraulic parameters (*θr*, *α*, and *n*) was performed using RETC software [17] by fitting the measured data. **Figure 7** shows the graph of the water retention curve data (circle) and the fitted water retention equation (line) for two upper soil layers. The fitting of the van Genuchten equation to data was confirmed by high *R*<sup>2</sup> values (>0.97).

**3.5. Available vadose zone models**

150 Groundwater - Contaminant and Resource Management

partially saturated conditions.

reactive solute transport in soil [30].

agricultural management systems [41].

series [40].

[43].

**Croatia**

in agro‐ecosystems [44].

A large number of programs are available for the description of water flow and solute transport in soil vadose zone. Below is the list and brief description of the most widely used software for modelling water flow and solute transport, which are applicable to unsaturated and

HYDRUS‐1D software includes one‐dimensional finite‐element model for simulating the

HYDRUS (2D/3D) is a software package for simulating water, heat, and solute movement in

MACRO is a one‐dimensional, process‐oriented, dual‐permeability model for water flow and

TOUGH ('transport of unsaturated groundwater and heat') suite of software codes are multidimensional numerical models for simulating the coupled transport of water, vapour,

SWAP (soil, water, atmosphere, and plant) simulates transport of water, solutes and heat in unsaturated/saturated soils at field scale level, during growing seasons and for long‐term time

The RZWQM is an integrated physical, biological, and chemical process model that simulates plant growth and movement of water, nutrients, and pesticides in run‐off and percolate within

Animo is a detailed process‐oriented simulation model for the evaluation of nitrate leaching

LEACHM (the leaching estimation and chemistry model) refers to a suite of simulation models describing the water and chemical regime in unsaturated or partially saturated soil profiles

Daisy is a dynamic model for the simulation of water and nitrogen dynamics and crop growth

The example will show the procedure of modelling water flow and nitrogen species transport in soil vadose zone. In this example, urea and NPK fertilizers were added to soil profile. The study site is located in the eastern Croatia in the intensive agricultural production area. The more information about the study can be found in [45]. The project main goal was to evaluate the influence of high fertilizer load, mostly nitrogen based, on the groundwater resources. Here, the results of modelling study on one selected site will be presented from year 2014.

**4. Water and solute transport using HYDRUS software: case study in**

to groundwater, N and P loads on surface waters and greenhouse gas emission [42].

movement of water, heat, and multiple solutes in variably saturated media [37].

two‐ and three‐dimensional variably saturated media [38].

non‐condensable gas, and heat in porous and fractured media [39].

**Figure 7.** Soil water retention curves for (a) 0–25 soil layer, and (b) 25–50 cm soil layer derived from RETC software.

After acquiring all the necessary input data, we can proceed to modelling with HYDRUS. In the main processes window (**Figure 8** left), water flow and solute transport are selected with additional root water uptake and root growth option since Barley (*Hordeum vulgare* L.) was

**Figure 8.** Main modelling process selection in HYDRUS‐1D (left) and input of soil hydraulic parameters (right).

**Figure 9.** Snapshot of soil hydraulic (left) and solute transport (right) selection window in HYDRUS‐1D.

grown at the site. After selecting geometry (length of the profile i.e. 50 cm) and time information (duration of the simulation, that is, 365 days) hydraulic parameters (*θr, θs, α, n, Ks, l)* for selected model are inserted (**Figure 8** right).

In the soil hydraulic model window (**Figure 9** left), there are multiple option to choose from, varying from single porosity to dual porosity/permeability model with different approaches. In this example, van Genuchten‐Mualem equation was used without considering the effect of soil hysteresis. Initial and boundary conditions were chosen. Atmospheric boundary condition (which includes precipitation and evapotranspiration) was selected at the top and seepage face at the bottom boundary (to mimic lysimeter plate). Root water uptake was simulated using Feddes [20] model. As for water flow, it is necessary to select solute transport parameters (**Figure 9** right). The number of solute are set to three since the simulation consider nitrification chain (urea > ammonium > nitrate). The first‐order reaction term representing nitrification of urea to ammonium was 0.38 per day [47]. The first‐order reaction term representing nitrifica‐ tion of ammonium to nitrate was 0.2 per day [47]. The first‐order reaction term for the volatilization of ammonium to ammonia was 0.0552 per day [48]. The distribution coefficient for ammonium (*Kd*) is assumed to be 3.5 cm3 g-1 [47]. The unrestricted passive root uptake of urea, ammonium and nitrate was assumed. For solute boundary condition, concentration flux and zero concentration gradient were selected for upper and lower boundary condition, respectively.

After inserting all the necessary data (derived from field observation and laboratory meas‐ urements of addition optimization) and running the model, pre‐processing (input parameters, **Figure 10** left) and post‐processing (results, **Figure 10** right) window are displayed in the HYDRUS.


**Figure 10.** Snapshot of pre‐processing (input data) and post‐processing (results) window in HYDRUS‐1D.

grown at the site. After selecting geometry (length of the profile i.e. 50 cm) and time information (duration of the simulation, that is, 365 days) hydraulic parameters (*θr, θs, α, n, Ks, l)* for selected

**Figure 9.** Snapshot of soil hydraulic (left) and solute transport (right) selection window in HYDRUS‐1D.

In the soil hydraulic model window (**Figure 9** left), there are multiple option to choose from, varying from single porosity to dual porosity/permeability model with different approaches. In this example, van Genuchten‐Mualem equation was used without considering the effect of soil hysteresis. Initial and boundary conditions were chosen. Atmospheric boundary condition (which includes precipitation and evapotranspiration) was selected at the top and seepage face at the bottom boundary (to mimic lysimeter plate). Root water uptake was simulated using Feddes [20] model. As for water flow, it is necessary to select solute transport parameters (**Figure 9** right). The number of solute are set to three since the simulation consider nitrification chain (urea > ammonium > nitrate). The first‐order reaction term representing nitrification of urea to ammonium was 0.38 per day [47]. The first‐order reaction term representing nitrifica‐ tion of ammonium to nitrate was 0.2 per day [47]. The first‐order reaction term for the volatilization of ammonium to ammonia was 0.0552 per day [48]. The distribution coefficient

urea, ammonium and nitrate was assumed. For solute boundary condition, concentration flux and zero concentration gradient were selected for upper and lower boundary condition,

After inserting all the necessary data (derived from field observation and laboratory meas‐ urements of addition optimization) and running the model, pre‐processing (input parameters, **Figure 10** left) and post‐processing (results, **Figure 10** right) window are displayed in the

g-1 [47]. The unrestricted passive root uptake of

model are inserted (**Figure 8** right).

152 Groundwater - Contaminant and Resource Management

for ammonium (*Kd*) is assumed to be 3.5 cm3

respectively.

HYDRUS.

Cumulative water outflow, which was measured in lysimeter at site 4 during 2014 and simulated ones using HYDRUS‐1D are presented on the **Figure 11**. The amount of water outflow from lysimeter was mainly the result of high precipitation events. The largest outflows were collected during wet part of the year, while during main crop vegetation season the outflows were very small due to large crop water demand (large transpiration intensity). Very high value of *R*<sup>2</sup> (0.96) between the measured and simulated values indicates that the HYDRUS model was capable to reproduce field data with high efficiency. Simulated values of nitrate outflow reflect the water flow pattern (**Figure 12**) which was measured in the same lysimeter (site 4) during 2014. Although the model derived larger cumulative nitrate outflows compared to the measured ones, the *R*<sup>2</sup> value of 0.72 indicate good model ability to simulate nitrogen transformation (from urea and ammonium to nitrate) and transport. The larger simulated nitrate values could be due to denitrification process that might occur in soil which was not considered in the modelling. From **Figure 12**, it can be seen that the main nitrate leaching occurs at the end of the simulation period after barley harvesting (day 164) and during wet period (autumn/winter). Our results indicate that numerical models can be very helpful for estimating water flow and nitrate dynamics under field conditions. From this example, it can be seen the possible negative influence of the nitrogen fertilizer application and their potential of leaching below root zone. One of the possible numerical models usages is their application in crop water demand (irrigation) and fertilization optimization or pesticide management in agriculture which can eventually lead to the protection of environment.

**Figure 11.** Observed cumulative water outflows from lysimeter (site 4) in 2014 and simulated ones using HYDRUS‐1D.

**Figure 12.** Observed cumulative nitrate outflows from lysimeter (site 4) in 2014 and simulated ones using HYDRUS‐ 1D.

#### **5. Trace metals mobility in soil‐plant system**

In some recent reports by European Commission [49], contamination by trace metal elements (Cd, Cu, Pb, Zn, Ni, Cr) and some nutrients (N, P) was defined as one of the main pressures to environmental resources in Europe. Contamination of arable soil resources by trace metal elements due to different anthropogenic activities is increasing rapidly and continuously in the last decades at the global level and often with detrimental ecological scenarios [50]. Trace metal elements represent relatively wide group of essential and some of the most toxic elements for biota which occur generally at very low levels in the 'clean' pedospheres. For instance, the total contents of trace elements in non‐contaminated mineral top soils range from 1 to 100 ppm [51], but in contaminated soils, their concentrations may be higher by up to several orders of magnitude.

Trace metal behaviour in the ecosystems, their chemical forms (species), mobility and risk of inclusion into the food chain, greatly depend on the environmental conditions. For example, the contamination of soil by trace metals is followed by a cascade of reactions with soil surfaces and the concentration of metals in the soil solution is controlled by a number of inter‐related processes: oxidation/reduction, precipitation/dissolution, adsorption/desorption, inorganic and organic complex formation [52]. Which processes will predominate depends on the physical, chemical, and biological properties of the (non)contaminated soil, as well as the important environmental factors (e.g. moisture, temperature, aeration). The most dominant parameters controlling soil trace metals chemistry and availability to plants are pH and soil OM. At lower soil pH values, hydrogen ions are adsorbed to soil particles, which increases the positive charge on inorganic and organic soil components, resulting in weaker adsorption of metal cations and their increased mobility in soil. SOM has a two‐sided role in metal mobility in soil: (i) particulate SOM is retaining trace metals through the formation of metal‐SOM complexes, thus decreasing their mobility [53], and (ii) dissolved organic matter (DOM) can increase trace metal mobility due to formation of metal‐DOM complexes [54] that are sub‐ stantially less bounded to soil particles than a free metal ions.

Particular metal forms in soil fractions (soluble, solid, liquid) are possible to detect by adequate analytical technique [55] or predict by some of computational models (e.g. next section). However, inside the cultivated soil (vs. soil without plants), even at very small scales (up to several mm), the physically‐chemical and/or biological characteristics may differ substantially, and thus consequently influence trace metal biogeochemistry. It is especially pronounced at the soil‐root interface, or so‐called rhizosphere microarea. From the biogeochemical perspec‐ tive, the rhizosphere is very dynamic and heterogenic microarea dominantly controlled by plant roots and released different organic metabolites.

Inside of relatively numerous group of trace metals, during the last decades, copper (Cu) has been intensively studied from different scientific perspectives given on several next facts: it is essential phytonutrient at small concentrations but easily becomes phytotoxic at higher levels; its biochemistry is pH‐dependent; it is one of the most reactive trace metal elements with soil OM and over the huge reactive (mostly negatively charged surface) interface of OM, Cu strongly competes with some other positively charged metals (e.g. Cd, Zn) and thus impacts on its bioavailability. In the next section, one of the most advanced biogeochemical modelling approaches elaborates the influence of particular soil OM substances on Cu chemical speciation in relatively homogenous and realistic rhizosphere conditions.

#### **5.1. Visual MinteQ example: copper chemical speciation**

**Figure 11.** Observed cumulative water outflows from lysimeter (site 4) in 2014 and simulated ones using HYDRUS‐1D.

**Figure 12.** Observed cumulative nitrate outflows from lysimeter (site 4) in 2014 and simulated ones using HYDRUS‐

In some recent reports by European Commission [49], contamination by trace metal elements (Cd, Cu, Pb, Zn, Ni, Cr) and some nutrients (N, P) was defined as one of the main pressures to environmental resources in Europe. Contamination of arable soil resources by trace metal elements due to different anthropogenic activities is increasing rapidly and continuously in the last decades at the global level and often with detrimental ecological scenarios [50]. Trace metal elements represent relatively wide group of essential and some of the most toxic

**5. Trace metals mobility in soil‐plant system**

154 Groundwater - Contaminant and Resource Management

1D.

Chemical speciation modelling was recognized as a very useful tool for studying various types of elements (nutrients, trace metals) and their mobility in natural systems such as rhizosphere. In the next, modelling in Visual MINTEQ chemical equilibrium program [56]) was performed for three levels of Cu (non‐contaminated and contaminated system), at three pH levels (covering the most naturally‐occurring pH reactions in different rhizosphere environments) and at three levels of soil dissolved OM (DOC; corresponds for mineral to organic soil types) (**Table 2**).


1 HA‐Cu: humic acid‐complexed Cu via 1‐carboxylic and 2‐phenolic functional groups. 2 Species with the < 0.5% of total concentration are not shown.

**Table 2.** Distribution (%) of Cu species in tested soil solution, estimated by Visual MINTEQ chemical equilibrium software (NICA‐Donnan model) as affected by soil pH (5, 7, and 9), dissolved organic carbon (DOC; 10, 20, and 30 mgL-1) and different soil Cu total concentration (40, 250, and 500 mg kg-1).

Even though the total soil Cu content by itself is not an adequate measure to determine Cu mobility and phytoavailability, a strong positive correlation between total element concentra‐ tion and its bioavailable fraction is often reported, especially in contaminated soils or in trials on soils spiked by metals [57]. In this model, chosen soil total Cu concentrations were (in mg kg‐-1); 40 (correspond to most of non‐contaminated soil conditions), 250 and 500 (corresponds from medium to highly contaminated soil conditions). The mobility of trace metals in soil depends ultimately on their chemical speciation, which is actually a function of pH and the presence of inorganic and organic ligands in the soil solution. The non‐ideal competitive adsorption (NICA)‐Donnan (sub)model for the adsorption of cations onto dissolved organic matter (DOC) and software database on equilibrium constants [58] was applied in this model. Cu concentrations used for Cu speciation here were average values obtained from experimen‐ tal data from CaCl<sup>2</sup> extracts, and for other elements in soil solution from saturated soil water extracts [53]. Ionic strength was set to be calculated and temperature of 25°C for all the calculations. Other software default settings were not modified.

From the data presented in **Table 2**, it can be seen that the majority of Cu in the modelled solutions was founded as a Cu‐DOC complex, suggesting that DOC might be the main factor affecting Cu ultimate fate after its release from the soil solid phase into the rhizosphere solution. Increase in total Cu concentration led to an increase in free Cu2+ ion in the modelled solutions, but only at lower concentrations of DOC and at acid pH (**Table 2**). Even though in every investigated scenario, the majority of Cu in the modelled solution was found to be complexed with DOC (**Table 2**), an increased soil pH caused the reduction of the percentage of free Cu ion in the solution. It is known that free metal ion is considered the most mobile species in the soil, thus data are implicating a decreased metal mobility in soil with an increase in soil pH. Furthermore, with an increase in soil pH Cu shifted from carboxylic groups in humic acids (HA1‐Cu (6)) to phenolic functional groups (HA2‐ Cu (6)), which is in agreement with some previous models that phenolic groups might be more important for metal (e.g. Cd, Cu) binding under higher pH of surrounding media [55].
