An Extensive Review of Leaching Models for Forecasting and Integrated Management of Surface and Groundwater Quality

1. Introduction

Surface and ground water quality management is an issue of high priority in the EU. Nowadays, the scarcity of soft water of high quality is a great concern for the human wellbeing and the future economy. Nitrates Directive 91/676/EEC introduced actions for preventing nitrates pollution of water bodies (surface & groundwater) for agricultural purposes, laying emphasis on promoting good farming practices [1]. The Water Directive (98/83/EC), followed up by EU Directive 2015/1787 with the amended Annexes II & III introduced water quality intended for human consumption including water monitoring on regular basis, quality standards, organoleptic and microbiological quality) [2-3].

Nonetheless, an integral statutory water management policy was lacking at that certain period. That gap was bridged in 2000 by setting in the Water Framework Directive (WFD), which brought in the front stage new concepts of environmental approaches with the ‘Good ecological status’ (abbr. GES) and the concept of integrated management at the spatial unity of ‘River Basin’ [4]. WFD incorporates mechanisms to prioritize standards as minimum level of chemical quality, regarding the detection of hazardous chemical agents. Within the well-predefined national river basin network, specific protection zones were to be designated in which more strict objectives should be applied on local basis if required. WFD incorporates a dynamic list of priority substances in EU region which outline the major risks and concurrently set all necessary cost-effective measures to abate all anthropogenic impacts and sustain the good water status of both surface and groundwater [4].

As complementary legislation to the afore-mentioned WFD, Directive 2008/105/EC introduced Environmental Quality Standards (abbr. EQS), in water policy field [5]. The Decision No 2455/2001/EC prioritizes thirty-three (33) chemical substances in soft waters [6]. Directive 2013/39/EU amended the WFD and EQS requirements and adopted an updated, enlarged list of priority substances on published Annexes [7]. Parameter monitoring and sampling is undertaken all year round.

European water quality legislation is based on United Nation’s Integrated Methodological Framework and is expected the deeper involvement of stakeholders in the implementation properly established river basin management plans towards water sustainability pathways. The former entails the issuing of Integrated Water Management Plans on local and regional basis.

Innovative ideas for the utilization of water resources are promoted through funded European programs as a priority area (Horizon 2020, Interreg etc.), based on legislation revision upon reuse of treated wastewater (industrial and domestic) for a variety of purposes e.g. irrigation, peri-urban green and agriculture watering, water use for industrial needs etc. [8-9].

In EU countries, numerous agricultural activities are running, many of which are breeding farms i.e. cattle farming, piggeries, poultries etc. and certain administrative regions are covered with cultivated plots, mostly energy crops (i.e. sunflowers, sugar beets), cereals (i.e. wheat, maize, barley, oats), cotton, potatoes and hydroponic green houses. All the afore given crops are irrigated by river waterbodies, pumping station abstracting water from local aquifers and water collection basins, via carefully engineered complicated irrigation network. Top soil fertilization, threatens soft water quality from increased salinity and nitrates pollution. Furthermore, agents rich in nitrogen (N), via their fate in the soil and inevitable propagation in rock beds causes ion mobility, mineralization etc. Aquifers’ quality change since e.g. new nitrogen pool enters the nitrogen cycle of the local ecosystem (Figure 1), unbalancing water-soil-nutrients equilibrium. Yet, humic part of the soil, microbiota, favors the formation of nitrous oxides a strong Greenhouse Gases (GHG) with a great impact of the climate change worldwide.

Water/pollute flow through porous media demonstrate a high complexity embodied via mathematical simulation models. Complexity increase explainability to the detriment of processing difficulties to be encountered. Ontologies describe in full detail, well-defined interrelated phenomena, suitably bound in conceptual structures over fluids through porous media which are applied in artificial intelligence and affect new infiltration model generation [10].

Nutrients’ leachate prediction models are estimating not only the residual nitrate part applied on the cropland but also all potential seeping losses of nitrate pool into varied soil profile layers, combined with local weather characteristics and farming practices. Artificial Intelligence (AΙ) and machine learning models are gaining an even higher importance [11] as an alternate reliable method for handling complicated nonlinear hydrological modeling with good countrywide scale predictions. In numerous cases, lack of sufficient and proper data, necessary to attain hydrological model calibration, combined with many parameters input for demanding, profound modelling, increase drastically the processing burden and undermine prediction accuracy and reliability. Obstacles are overcome by using models based on AI and machine learning techniques presented in many publications which reveal the future trends regarding surface and ground water quality decision making [12-15].

The scope of this paper was to conduct a systematic review of leachate models useful for proactive and rehab actions to safeguard surface and subsurface soft water which become even more scarce. Well-tested solute leaching models standalone or as a part of an integral modelling package provide rapid site-specific estimates the leaching potential of chemical agents, mostly nitrates, below root zone of crops and the impact of leaching towards groundwater. Most of the models examined were process-based or conceptual approach. Nonetheless, empirical prediction models, though rather simplistic and therefore not preferable, demonstrate certain advantages since they do not require extensive calibration database information which in many cases is unavailable not to mention stochastic approach and the involving of AI.

2. Methodology

In order to isolate the suitable models adopted and agent’s fate methodology, publications derived from Science Direct®, Web of Science™ and Google Scholar huge scientific data bases. Specific keywords, i.e. ‘model’, ‘nutrients’ fate’, ‘porous medium’, ‘infiltration’, ‘percolation’, ‘nitrogen leachate’ were tested, combined or stand alone, to track down and compile the proper material which it was carefully scrutinized to extract the relevant cited content.

All findings were classified according to the emphasis drawn by the authors of the paper regarding the field of interest, though it was more or less expected a cognitive overlapping. For instance, model packages widely used, predict nutrients fate, and soil-water, soil-air balance and thereof types of crops to be planted in the near future, as well as quantitative approach of the suitable fertilize use. A basic classification was made among nutrients, ion metals and other chemical agents’ leaching process. The porous medium which stands as a hindrance to the pollutant flow justifies new reasons for further categorization. Temporal and spatial scale information, soil use along with climatic data are of great importance for any assessment. All model packages incorporate significant components which are to be calibrated by timeseries data and previous relative research, so as more reliable results to be obtained. Components could be classified as physically based or conceptually distributed models, theoretical or mechanistic (i.e. process-based or conceptual), empirical, hybrid, machine learning etc. In many model-packages prediction applied entails certain different dominant approaches governing each component. As regards the solute, simulating models predict either groundwater flow or solute transport, boundary transition zones which employ numerical models with specific conditions e.g. where salt intrusion in aquifers is to be monitored. Further categories are established based on the final prediction technique employed viz finite difference against finite element techniques [16].

3. Materials and Methods

The majority of the models applied are mechanistic models (process based), which describe satisfactorily dynamic processes and exploit accurate data used for model initialization. Mechanistic models are categorized as progressive conversion models where particles size remains unchanged during the leaching process without any penetration resistance to fluids/agents, shrinking core model with constant particle size, shrinking core model with shrinking particle size etc. [17].

The way of approach for each case selected, determine conceptual, empirical, physically based, hybrid models, steady (snapshot like pictured) or dynamic (ongoing) simulation models. Another significant categorization depends on the extent of area appliance i.e. pasture, cultivated plot, river catchment/basin, estuarine, coastal waters, aquifers, regional, countrywide etc. and dimensional fate 1-D, 2-D, spatial analysis employed e.g. distributed of lumped and huge heap of data collected and tested from data bases or models with minimum data requirements [18-19]. Groundwater modeling, when numerical models are employed, entails certain approaches e.g. finite difference [20], finite volume [21], finite element and element-free [22] methods.

Chemical agents e.g. herbicides, complex ions, organic pollutants concentration, saline and fate, runoff, drainage, infiltration, percolation, seepage phenomena etc. determine water quality deterioration, monitoring network, for integrated strategy and remedial actions. In complex soil layers, anion retention capacity varies and surely soil nitrates flux does not correspond directly to the water soil infiltration flux and it is rather a function of the organic matter content. However, major categories located and analyzed in published papers of modelling are given below:

3.1. Process-based or conceptual models

Process based model employ diagrams and interlinks among parallel processing and flowcharts of mass transfer in between different zones (entities) to facilitate complex physical and biochemical interactions of soil, crops and water bodies [18]. Soil particles redox potential, electrical conductivity etc. are ions transfer diffusion phenomena and leaching controllers [17].

Organic chemical screening is done using leaching models, which account for the mobility of the chemicals in the soil. Despite not being able to forecast concentrations, they are still straightforward and only need information that can be found in reports from soil surveys and scientific publications. They can also be readily coupled with GIS to produce an effective evaluation tool for the control of non-point sources in agricultural watersheds [23-24].

Simple models with semi-empirical process descriptions of the lumped conceptual type which utilize time series of averaged accumulated water e.g. on catchment level within predefined blocks [25] are presented next, such as: ANSWERS [26], for watersheds where sub-catchment are divided into hydraulically-interconnected elements (grids), a model structure approach known as ‘fully distributed’ that demands heaps of data and requires high computational burden, CREAMS [27-28] field scale model for chemical runoff and soil erosion, GLEAMS ([29-30], developed to assess agricultural chemicals effects on water quality, SWRRB [31], a modified CREAMS, AGNPS [32] and AquiMod [33]. In such models there has been a flux between segregated water body blocks following physical principles i.e. surface flow, groundwater recharge, infiltration, percolation, subsurface flow, water table recharge etc. [25].

Conceptually, leaching encompasses both the top layer of soil and the area beneath a plant's roots depicting the soil-ground water system. It could be considered that a soluble part of pesticide mass is instantly mobilized at the surface causing surface water infiltration into the rhizosphere.

In pastries, soil type contributes to 12% when referring to (N) leaching and a smaller contribution to (N₂) fixation. Still, soil profile and type remain the most significant factors involved in % (N) loss [34]. The agro-ecosystem model, EcoMod with followed up versions was developed as biophysical, pastoral simulation model [35]. Plants’ growth rate and carbon flux on the soil were simulated. Model modules regarding water and nutrients’ presence, incorporate calculations for water and solute circulation, nutrients’ fate, leaching, adsorption, ammonium nitrification, gaseous (N₂) losses, other geochemical phenomena i.e. gradual immobilization, mineralization etc. Fertilization and irrigation practices are included as well. Climate conditions were based upon data selected by means of stochastically created 99-year series climate documents by employing a modified version [36] of a daily climate model and incorporating data derived from the ‘Stochastic Climate Library’. By utilizing EcoMod, better regional management is applied via the feeding animal type for pastry to be selected i.e. cattle or sheep and the adequate use of nitrogen-based fertilizers to be applied, for optimal results over sustaining the adequate pasture, without all negative leaching affects caused by nitrates pollution over the catchment of the zone of experimental interest (New Zealand). Nitrogen leaching appeared to be more intense during winter time. Furthermore, the application of nitrogen-based fertilizers, though increased quantities of monthly nitrogen leaching and fate, still did not alter the pattern of % nitrogen loss i.e. the cyclic pattern of nitrogen fixation [34].

A biogeochemical model, based on DeNitrification-DeComposition processing (DNDC) was developed by [37-38] for simulating dynamics of both carbon and nitrogen forms in soils for agro-ecosystems. DNDC simulates crop yield, carbon sequestration in soil and nitrogen leaching when irrigated. Furthermore, it estimates greenhouse gas (GHG) emissions from terrestrial ecosystems and carbon sequestration [39-40]. A basic component of the model simulates crop growth and decomposition by processing data of soil conditions i.e. moisture, redox potential, temperature etc. and another component covers carbon and nitrogen transformation though aerobic nitrification, anoxic denitrification even anaerobic fermentation by means of biotic micro-population in soil. Besides, biogeochemical interactions are taken into consideration and nitrogen final releases in atmosphere or even leaching in deeper soil sequestration after hydrological processing, by means of nitrogen adsorption phenomena, plants’ uptake in rhizosphere, microbial assimilation, nitrification/denitrification etc. [41-42]. Zhang et al., (2021) [43] applied DNDC biogeochemical model, to quantify denitrification-decomposition bioprocessing, as a mechanism, a part of overall nitrogen management of Greenhouse vegetable fields (see Table 1).

The DayCent model [44] simulated nitrogen cycle in soils for various ecosystems (e.g., cropland and forest). It simulates N₂O emissions and mineral nitrogen leaching under different types of landscape exploiting a variety of meteorological data, soil properties, irrigation practices, crops particularities. More specifically, nitrogen leaching presents remarkable spatial ana temporal variations, known as hot spots and hot moments (HSHM). Nitrogen leaching and loss within HSHM, could be easily tracked and managed by means of applying the proper model to achieve suitable fertilizing management, local climate factors, soil profile properties and the local topography [45].

TETIS [46-47] is a process-based, hydrological, conceptual, distributed model, that simulates aquifers’ recharge and soil-nitrogen leaching enabling the study of the long-term groundwater quantified response to nitrogen leaching. It incorporates precipitations data from meteorological time series, different irrigation and fertilization practices and potential evapotranspiration. The latter is obtained by means of Hargreaves-Samani equation use [48].

The nitrogen cycle is dependent on fertilization and atmospheric deposition. Corine land uses, maps and ‘Pedotransfer’ functions were involved [47], [49]. They were employed in combination with functions exploiting raw soil data upgraded into useful information tools via predictive functions of certain soil properties. 50 years meteorological data were combined with irrigation and fertilizing practices as well as the digital elevation model from the Geographical Survey of the site of interest. All selected data were the input of TETIS model. Finally, the quantification and sensitivity of nitrogen recharge and nitrogen leaching corresponding to a variety of weather conditions were estimated on certain temporal bases in Valencian region (location in Spain) [50].

STICS, a one-dimensional, daily-step, soil–crop, dynamic model [51-52]. It simulates the quantities of daily drained water, (N) leaching in soil-water phases, (N) plant uptake and fixation, (N) fluxes and dynamic interrelation between hybrid tested crops and soil type under various climatic parameters e.g. air temperature, rainfall density, solar radiation, ‘Penman’ evapotranspiration, C:N ratio, pH, calcium carbonate equivalent action etc., which finally affects crops’ biomass production [53]. Pedoclimatic characteristics (i.e. soil water content, nitrogen abundance in soil, bulk density, soil infiltration rate) and soil exploitation practices are used as initial inputs for the running of the simulation. Nitrate is transferred by convection mixing phenomena. A calibration stage is mandatory before the appliance of new cultivar parameters.

Grain legume-based cropping systems require a comprehensive reevaluation, by taking advantage of bioprocesses e.g. symbiotic biota for fixation and mitigation of leaching (N), nitrogen cycle cropping systems reengineering, by adopting crop succession, improved concepts as regards nitrogen-based fertilizers appliance, irrigation practices and techniques. The solute transfer between subsoil profiles is conformed with tipping bucket concept i.e. downwards mitigation after surpassing soil’s profile water capacity and the prevention of aquifers’ nitrates pollution, focusing on soil-nitrogen balance to avoid (N) leaching, in accordance with Direct. 91/676/EC concerning agricultural activities [53].

Model of standard version 9.2 of STICS [54], is a soil-crop model which computes changes in biomass and yield in soil organic carbon (carbon flux), predicts nitrate leaching, soil water and nitrogen fate etc. by interpreting certain dependent variables and parameters e.g. weather conditions, cropping practices, [52], plant growth, carbon as (CO₂) and nitrogen fluxes [55] and solar radiation.

The latest version of STICS in based on monitoring the changes of carbon pool as expressed by the Gross Primary Productivity (GPP), dependent on autotrophic respiration, Net Primary Productivity (NPP) via plants’ growth on a predefined time basis and finally the complex component of Ecosystem Respiration (RECO) which is comprised from autotrophic plant biomass, plant nitrogen, heterotrophic residual mineralization of the organic matter and Net Ecosystem Exchange (NEE), which equals the summation of GPP and RECO. NEE represents the balance between carbon fixation via photosynthesis and respiration, that releases it to environment. Carbon flux variations are being recorded and underwent comparisons upon crops rotation [55].

AqYield model [56-57] is a simple dynamic model which requires only a few inputs (e.g. soil properties, daily climate features, dates of sowing, harvest, irrigation, tillage irrigation and soil tillage depth for crop management). The model estimates sufficiently well drainage, water flows for several crops, major nitrogen flows in the soil, plant (N) uptake and mineralization. It estimates soil-water content dynamics, water drainage and evapotranspiration. It is tested and evaluated for the precision of the predictions accomplished. Although fairly simply, it is considered to be accurate in terms of (N) leachate prediction in corps rotation as much as other more sophisticated similar models e.g. STICS. AqYield predicts satisfactorily (N) flows (fate) in the soil–plant system, including mineralization, plant uptake and mitigation on daily basis. Cultivated crops, climatic data and soil characteristics have a great gravity in shaping the final predictions.

Great bibliographical research was conducted by Baveye, (2023) [58] to emphasize the steps required for a model development to predict soil carbon dynamics, interrelated to nitrogen cycle, mostly derived from fertilizers. According to his argument, current ecosystem-scale models are introducing limitations when applied to ecosystems due to lack of the proper interdisciplinarity and “bottom-up” approach, making use of microscale soil processes. Furthermore, interesting for stipulating empirical models, applicable at large spatial scales is still valuable when upscaling efforts is the case. Thus, modelists should stay focused on micro-scale of soil microorganisms at first and progressively should try to move up regarding the spatial scale, incorporating characteristics of the region of interest. Already kwon or future empirical model are providing valuable services for the final upscaling target.

Topsoil Carbon undergoes various transformations mainly as Dissolved Organic Carbon (DOC) and Particulate Organic Carbon (POC). The latter has a vivid role in numerous flocculation and adsorption phenomena affecting soil physicochemical properties. DOC represents the soluble organic part and holds a crucial role in soil transformation via micro-organisms’ uptake availability, the high soil mobility/leaching, biochemical transformation by means of metal ion/organic chemical binding and toxicity phenomena through biosorption. [59], developed and introduced a modified TRIPLEX-DOC model focused on leaching, sorption, DOC and biodegradation phenomena [60] to simulate DOC dynamics in monsoon forest ecosystems since it does not incorporate itself a component simulating responses of DOC concentrations and fluxes to the atmospheric N deposition. By means of TRIPLEX-DOC model improvements, a more profound biogeochemical modelling was incorporated to predict more accurately DOC fluxes in forest ecosystems in the monsoon regions which differentiate them in terms of vegetation and NPP. Furthermore, the modified model predicts more accurately regional carbon cycle in seasonal patterns and assess more safely soil DOC fluxes in local ecosystems.

Nitrate leaching occurs in croplands where synthetic and animal biowaste-nitrogen based fertilizers are applied, which affect vulnerable aquifers (e.g. karsts) from nitrates’ pollution. Soil and Water Assessment Tool (SWAT), [61] developed by the US-Department of Agriculture (USDA) and Agricultural Research Service (ARS), is widely exploited to simulate the hydrological cycle, plant growth, chemical leaching in numerous watersheds and cultivated soils. SWAT(CATCH) is a distributed, semi-empirical [62], physically based, hydrologic model (DPBHM), which incorporates the Variable Storage Coefficient (VSC) method. The latter is a kinematic wave method which means that streamflow equation is simplified by assuming that the friction slope is approximately equal to the slope of the e.g. watershed channel. The proper use of SWAT simulation corrects initial assessment of certain indices e.g. Dendritic Connectivity Index (DCI) a measure of longitudinal connectivity, of aquatic ecosystem status which outlines that temporal resolution of river water routing (time steps) plays a dominant role to the final assessments. It is noteworthy that time steps are dependent on actual watershed size and configuration [63]. SWAT employs the Vegetative Filter Strip MODel (VFSMOD) which is an empirical model [64] that calculates sediment reduction. It simulates crop yields and nitrate leaching under a variety of nutrient and irrigation management practices, as a well-tested prediction tool to improve management practices upon safeguarding a high groundwater quality towards agricultural sustainability. SWAT is combined easily with Sequential Uncertainty Fitting (SUFI-2), Uncertainty Procedures (SWAT-CUP 2012) and the Nash–Sutcliffe model efficiency (NSE) validation test [65].

SWAT is being modified by coupling with shuffled frog leaping algorithms and a farm-level economic model and cost estimator (FEM). The model outcome (MOSFLA), demonstrates powerful convergence and optimization ability [66].

SWAT tool predicts and quantifies nitrate leaching pathways and prevent ground water nitrate pollution which threatens drinking water resources and unbalances ecosystems. Nitrate pollution assessment change fertilizing practices, alter the crops and prioritize new sustainable agricultural methodologies and management practices in agribusiness. The crucial nitrate pollution drivers are the type of the soil, climate variability, crops irrigation and fertilization practices over cropland [65].

The coupled SWAT (version 2012/Rev664) and MODFLOW model [67-69] performs better when complex surface-groundwater interactions analysis and explicit modeling of the groundwater system and surface water is the case [65], in comparison with field-scale cropping system models such as Decision Support System for Agrotechnology Transfer (DSSAT) [70], HYDRUS 1-D [71-72], Root Zone Water Quality Model (RZWQM) [73-74] and Leaching Estimation and Chemistry Model LEACHM [75]. All the aforementioned cannot accept as a complementary component a hydrologic model so as to drive to better predictions.

Modified croplands and watershed SWAT model coupled with the groundwater model and stream–aquifer interaction hydrologic MODFLOW compose the improved SWATMOD model [76], which demonstrate simulation capability of surface water, stream-aquifer interactions and groundwater [16].

HYDRUS-1D model was employed and calibrated before use against nitrates concentration in soil-water across all soil layers and the measured soil moisture. It validated independently and finally started to simulate N-NO₃ leaching under certain conditions: e.g. controlled fertilization rate within prearranged irrigation depth or vice versa i.e. regulated irrigation depth combined with a steady fertilization rate. Concurrently, were assumed different shifts in water table depth. Experimental site-specific soil physical and hydraulic properties were selected as well as data from an automated meteorological station i.e. net radiation (Rn), barometric pressure (BP), temp. (T), wind speed (W), rel. humidity (RH) and precipitation at 60 min intervals [72]. It is commonly employed in landfills hydrological evaluation [77].

The Penman-Monteith equation was utilized in conjunction with the meteo data obtained to determine the reference evapotranspiration [72]. HYDRUS-1D employs Richards’ equation [78] to simulate the unsaturated water flow [79]. Mass balance of nitrification - first-order reaction rate constants - and denitrification along with N-NH₄⁺ volatilization rate was taken into consideration. Soil-water characteristics curve was depicted by using van Genuchten-Mualem equation [80] and pedotransfer function model [49,72] coefficients and parameters related to the hydraulic transport and transformation (see Figure 2). Solute transport in heterogeneous waste rocks could be simulated by HYDRUS-1D since it supports multiple porosity components [81].

Dynamic leaching flux and soil water storage including dissolved CO₂ and N₂O conc. can be simulated with good results by implementing HYDRUS-3D finite element model [82-83]. It entails domain, finite element mesh, defined boundary conditions and the input of topsoil vertical layers based on Digital Elevation Model (DEM) surface and topsoil thickness information. It estimates soil and bedrocks hydraulic properties incl. CO₂ and N₂O storage along with their leaching flow rate. Appliance of HYDRUS-3D over topsoil reveals the critical influencing, leaching factors (i.e. temperature, fertilizing, precipitations), of diluted CO₂ and N₂O (the latter a strong GHG) in water and thus promotes the deeper understanding of nitrogen cycle mechanism on terrestrial ecosystems [83].

Efforts were made to estimate (N) leaching quantification of irrigated maize by testing different fertilizers, regulating, through pumping and irrigation, the water table and adopting improved irrigation management practices. During these efforts, conventional fertilization practices were tested as a reference scenario and a field experimental trial was conducted for consecutive years. Subsoil zones were identified in which the leaching fluctuates within a hydrologic year [72].

Land use investigation combined with climate change was analyzed by implementing a global vegetation model (LPJ-GUESS) [84]. Trade-off between crops yield and nitrogen soil leaching were drawn in certain scenarios of an increase in wheat and maize yield for various regions. Integrated assessment model MESSAGE was employed simulating fertilizing appliance, making use of drivers such as (Representative Concentration Pathways 8.5) [85]. Distribution of nitrogen application rate information was derived by using (LPJ-GUESS) model after being fed with data selected from the Common Agricultural Policy Regionalized Impact (CAPRI) model at a terrain grid level up to 1km². Climate and land use changes projection, assess nitrogen leaching, CO₂ mitigation and crops yield [86].

A new approach of monitoring nitrates leaching of land fertilization is obtained by the NIT-DRAIN model [87]. It simulates nitrates conc. at a subsurface drainage network's outflow. It is a drainage discharge simulation model instead of monitoring/observation model. Water and soluble nitrates are transferring through the soil profile which is separated into three successive, interconnected, subsurface stratifications, segregated by drain and mid-drain layers. Thus, there is a distinguishable water-nitrates fast transfer, throughout the macro-porous upper layer, located above the pipe network. Low and residual nitrate transfers are accomplished through the lower conceptual soil stratifications. Nitrogen load to the ground is divided into two discrete phases. Firstly, the seasonal applied nitrogen quantity which entails all predicted nitrogen transformations at the cropland scale and secondly, the remaining off seasonal nitrogen pool which undergoes extended leaching phenomena when the winter season starts. Both prediction data are required. Model also requires updated variable input at the commencement of a new hydrological cycle. It involves seven nitrates leaching and water flow parameters and employs two different transfer functions as long as nitrogen fate lasts. Model calibration based on Generalized Reduced Gradient (GRG), non-linear regression algorithm and the hierarchical scheme for validation, by making use of the split-sample test [88].

GOSSYM a mass balance mechanistic two-dimensional (2D) gridded soil model that simulates water, carbon and nitrogen interactions in soil, plant root zone and crop response to climate variables and water irrigation. It incorporates many routines as components such as daily weather information, soil profile characteristics and parameters, plant growth data etc. Modification is presented by coupling 2DSOIL mechanistic 2D finite element soil process model. The latter, simulates underground processes i.e. mineralization, of organic matter immobilization, nitrification and denitrification, [89]. GOSSYM is also to be used coupled with COMAX expert system developed ad hoc. The inference engine incorporated is very useful for cultivation practices, model results interpretation and irrigation and fertilizing regulation [18, 90]. GOSSYM-COMAX is widely validated by being tested under different environmental conditions and crop practices.

The new climate regime over certain regions in Europe incurs even more extreme rainfall incidents with catastrophic ecological results. The abstraction of the humic topsoil tends to be a severe negative outcome. Therefore, is of crucial importance the use of integrated biogeochemical-ecosystem models and the estimation of intra-/inter-annual gross primary productivity (GPP) variations e.g. (Biome-BGC) [91-93] to monitor ecosystems response against extreme events. Interactions of ten hydrological and ecological processes were taken into consideration, incl. runoff, canopy evaporation, soil-water storage, soil evaporation/respiration, transpiration, Net Primary Productivity (NPP), net ecosystem exchange, nitrogen mineralization and leaching. Once more, nitrogen propagation incurs nitrate pollution and is considered to be of high priority concern when the Good Environmental Status ‘GES’ or even the Good Ecological Status is to be assessed. The model focuses on carbon (C) and nitrogen (N) fluxes on a preset time scale (i.e. daily, monthly, annual). The model parameterization is achieved by making use of site characteristics data (e.g., soil texture/relief/depth, terrain elevation), meteorological data (e.g. precipitation, humidity, temperature, etc.) and eco-parameters (e.g. canopy growth and expansion, photosynthetic rate, carbon to nitrogen ratio of leafage, lignin fraction of the dead wood) [94].

As regards the part of the daily soil-water balance, the one-dimensional soil-plant-atmosphere DAISY, Danish model simulation (version 4.01) [95] was employed to cover the water and N balance in agro-ecosystems for crop cultivation [96]. The model comprises three main modules, i.e. a bioclimate, a vegetation, and a soil component. Soil-water retention and transport were simulated by employing DAISY model and Using the Makkink equation, evapotranspiration was determined from the daily global radiation and mean temperature [97]. Variables for nitrate leaching i.e. nitrates percolation entail crop classification into groups according to (N) uptake performance, residual (N), contribution to leachate and (N) mineralization i.e. transformation into inorganic form (NH₄⁺) [98].

Nitrate leaching predictions and nitrogen conc. in monitoring sites in Denmark, incl. various parameters i.e. soil type, climatic data, crop consequence, nitrogen fertilizers appliance rate, soil winter capping etc., were conducted making use of multivariable experiments. A linear regression was applied to identify and assess the influence gravity of a leaching year (hydrological year). A certain number of sites was selected and DAISY model (ver. 4.01) [95], (Figure 3) was implemented to simulate water percolation and balance at the depth of the rhizosphere. DAISY model is employed for quantitative description of water flow into unsaturated soil zone and makes use of certain soil quality characteristics i.e. water retention, unsaturated hydraulic conductivity which earlier were fitted to Mualem-van Genuchten equations [80, 99] and RETC optimization computer programming code [100]. In the specific case DAISY model simulates the water balance regarding different soil-crop productions. The model incorporates Richard’s equation [78] to calculate soil water flux. Different soil sites and differentiated climatic data deduce new estimated values of nitrates leaching, therefore proper cropping strategies should be followed for soil nutrients balancing [101].

‘MIKE SHE’ is a complex model package which describes distributed and physically-based water/solute flow in catchment areas. It incorporates snowmelt process along with evapotranspiration and gives numerical solutions for overland (2D), unsaturated flow (1D), channel flow (1D) and saturated flow (3D). The combination of ‘Daisy’ and ‘MIKE SHE’ models enable water and nitrates simulation transport and overland flow on a catchment level, excluding intense rainfall events i.e. Hortonian incident [102].

Strict cropland’s, nitrogen-based fertilization regulations in (kg N ha^-1), in Europe, compel farmers to resort to alternative cultivating options i.e. catch crops, to regulate post harvesting nitrogen conc. of soil-rhizosphere nutrients balance. Some cultivations were tested over time e.g. winter rye, spring barley, fodder radish, to achieve lower nitrogen leaching results. Simulation of the cultivation change, during winter and spring, by using the biophysical Agricultural Production Systems sIMulator (APSIM), simulation model [103] (Figure 4), is parameterized over mineralization of catch crops. Nitrate leaching was based on weighed percolation drainage within periods of interest. APSIM comprises certain modules i.e. SurfaceOM, SoilN, SoilWat, Canopy model and Cultivar plant type. Each one of the modules simulates different fields of interest e.g. water movement, N and C soil dynamic, plants’ solar radiation competition etc. The SoilWat is a module of APSIM. It has been developed as a cascading water balance model in which, on a daily basis, water balance is estimated by simulating runoff, saturated/unsaturated water flow, solute ground and underground movement interrelated with saturated/unsaturated flow, plant transpiration and soil evaporation. Soil Nitrogen module simulates nitrification, denitrification, urea hydrolysis and nitrogen mineralization [104]. Nitrogen originated mainly by manure appliance over cropland and inorganic based fertilizers [105].

APSIM’s appliance entity (i.e. plant-rhizosphere-soil-atmosphere), functions on daily time frame [103]. The prediction assures reliable crop management strategies, all year long, for the sake of a balanced ecosystem, minimizing all negative environmental impacts, as regards crop yield, water, nitrogen and carbon mitigation dynamics [105]. Moreover, this model has been developed for sugarcane simulations, tested and confirmed at the paddock scale. since it incorporates modules for specific crops [106-107]. Thus, it simulates management practices, water uptake and growth of sugarcane [104].

Intermediate percolation calculations conducted by using the model EVACROP v.3.0, updated of EVACROP v.1.5 [108-109]. Nitrogen leaching affects water quality and soil acidification. A bio-geophysical, process-based, multi-component ecosystem model (CoupModel) [110] was employed to estimate carbon and nitrogen fluxes, water and heat balance in eco-terrestrial cycle. It facilitates the assessment of the local ecosystem (N) balance e.g. manure appliance over tillage serving as biofertilizer. Nitrogen transformation includes (N) plant uptake and mineralization. It runs on a daily time frame and on a minimal squared meter unit grid. Model running requires certain parameters i.e. weather (precipitation) data, solar radiation, temperature, relative humidity, wind speed, precipitation and dry deposition as regards (N) information.

Carbon and nitrogen storage incorporates abiotic/inorganic i.e. ammonium and nitrates in soil and biotic part e.g. plant tissues and external characteristics such as leaf area and solar radiation expose which regulate carbon assimilation. The model regards the soil part as a series of constituent, interrelated layers, where water, heat, nitrogen and carbon interact in a dynamic way. Darcy’s law and the generalized Richards (1931) equations [100] for unsaturated conditions are governing soil water flux.

Soil evaporation phenomena were calculated by using the Penman-Monteith equation [111]. Outflows of the ground water table determined by a linear empirical drainage function. Soluble Organic Matter (SOM), (N) mineralization, plant’s water uptake is dependent on moisture and soil temperature. Other dissolved forms of (C) and (N) comprise nutrients’ pool in the roots. Microbial biomass of the seeded roots (Mycorrhiza), regulates stoichiometry of C/N ration and Carbon Use Efficiency (CUE). Any changes in the nitrifying biomass of a specific soil layer, affects eminently the soil nitrification flux and therefore the nitrification rate and the soil-water-ammonium conc., since Dissolved Organic Nitrogen (DON) and ammonium (N-NH₄) are in chemical equilibrium in both soil-water and soil-solid phase [112].

The Nitrogen Loss and Environmental Assessment Package (NLEAP) is a powerful mechanistic model that provides rapid site-specific estimates of N-NO₃ leaching potential below root zone of crops and the impact of nitrate leaching into groundwater [113]. Under various combinations of soil profile characteristics, weather data, and farming practices, the afore given model is used to estimate the residual nitrates and potential nitrate leaching losses from cropland soils into subsurface layers [114]. It runs in three temporal modes i.e. annual, monthly and event-by-event for certain water-Nitrogen cases and potential leaching.

The updated and improved NLEAP-GIS, (version 4.2), coupled with Geographic Information System (GIS) spatial soil database, as a complementary component, increases the capability of assessing (N) losses towards risky landscape and combination of different cropping systems and evaluates more accurately management practices over nitrogen transform and mitigation [115-116].

Updated version model includes infiltration and transport of nitrates and soil water; carbon and nitrogen terrestrial cycle, soil-surface and soil-profile propagation and transformations, water, nitrate/ammonium surface runoff, nitrate leaching from the root zone and ammonium uptake from the seeded plants, denitrification losses and ammonia volatilization. Each individual process is prescribed by detailed algorithms and equations [116].

Model inputs such as daily meteorological data, incl. maximum and minimum temperatures, daily rainfall, relative humidity and evaporation are prerequisite to run various scenarios enriched with soil profile information and fertilizing events. During winter wheat and summer maize season, specific (N) application rate is recommended, so as to be achievable, concurrently, higher (N) use efficiency and grain yield resulting in lower groundwater nitrate pollution risk. Nitrate leaching indicates the presence of nitrogen pool which might affect local aquifers identifying regional soil profile characteristics. Model’s outcome is being validating by using older measured nitrate leaching reports. Results in many cases are comparable with other models’ outcome e.g. DNDC model [117] on regional scale [118] (see Table 1).

The large-scale and spatial variation simulation in nitrate leaching in the region of experimental interest (location in China), the rate of fertilizer N application was positively correlated with rainfall density, however, it was negatively correlated with N use efficiency [119].

Climate change and interrelation of anthropogenic pressure to waterbodies’ catchments/basins are of great importance. Therefore, evaluation and simulation of Water Resources Management Scenarios for adopted, resilient plans over coastal agricultural watersheds were conducted. It was achieved by the implementation of an Integrated Modelling System (IMS) that was set up to cope with the complexity of croplands in vicinity with estuarine in the Med zone [120]. IMS incorporates certain components of reservoir operation (UTHRL) [121-122] surface hydrology (UTHBAL), named as (R-UTHBAL) built in R statistical language [123], groundwater hydrology (MODFLOW) [67], crop growth/nitrate leaching (REPIC) which is a coupled Environmental Policy Integrated Climate (EPIC) [124-125] model with R-ArcGIS bridge application (ESRI), nitrate pollution flow, multi-species, transport model for simulation of dispersion, advection and chemical reactions of contaminants in groundwater systems) (MT3DMS) [126] and the prediction model for seawater intrusion cases (SEAWAT) [127].

When it comes to groundwater dynamics, aquifer interactions with the stream system, surface-water and nutrient fluxes at the watershed outlet, and water and nutrient leaching from the surface to the aquifer, the Integrated Surface and Subsurface Model (ISSM) [128] produces good prediction results. It consists of the in-stream water quality model (QUAL2E), the groundwater models MODFLOW and MT3DMS, and the hydrological model SWAT.

PATRICAL is a conceptual large-scale (area grid 1×1 km²) temporal and spatial distribution model, addressing water quality and balance with remarkable time projection prediction. It enables an overall aspect of nitrate pollution in extensive regions and facilitates countermeasures for aquifers and wetland recovery to meet EU WFD requirements. It entails calibration, validation and future scenarios. When ground water bodies are the case, it incorporates tens of lumped models [129-130]. Many other models and components of minor importance and more information is figuring in Table 1.

3.2. Empirical Models Or Statistical Models

There have been cases where conventional models are considered to be inefficient when the case is soils with cracks and high heterogeneity in terms of soil texture, particles shape and orientation, burrows etc. Thus, diffusion coefficient is difficult to be adopted for extensive tortuosity. In such cases many experimental data are needed or suitable empirical correlations to be employed [17].

Nonetheless, empirical models are lacking satisfactory results when extrapolation is needed in the future. Despite this is still a useful tool in the decision making in farming or grazing [18]. The NLES5 model [98] is partially modified of an earlier NLES4 model version [131] since it has almost the same structure and share partially datasets for calibration. It considered to be an empirical prediction model of quantifying nitrate leaching on annual basis, considering the interest zone as the root zone i.e. 1 meter depth. The model accounts for incurred effects of added nitrogen (N), under certain crop sequences, seasonal crop cover, top soil types, and local weather conditions. Soil (N) leaching and mitigation prediction, caused by fertilizing over croplands, improved integrated water basin management and promoted the good quality preservation of surface water and aquifers as the main goals of model development.

Quick Plant RZWQM2 sub-model, is a simplified empirical plant module, developed to simulate the water and nitrate fate for the crops, therefore to compute on daily basis, water and nitrogen uptake of the cultivated soil (see also 4.2 below). Meteorological data collection is essential to the final outcome. It incorporates Green-Ampt infiltration equation, Richards' equation [78] for water redistribution and a plethora of submodules such as: Quick Plant, soil water and heat transfer, N balance, generic plant growth, evapotranspiration (PET), soil equilibrium chemistry, pesticide and finally management module. Soil water retention curves are acquired by using the modified Brooks-Corey equation.

Empirical prediction models demonstrate certain advantages since they do not require extensive calibration database information which in many cases is unavailable [132]. Any future strategy must address the critical issue of relatively accurate nitrate leaching prediction across common cropping systems, soil types, and climate conditions in specific geographic areas, as well as the impact of seasonal vegetation characteristics on nitrate leaching, which is dependent on (N) fertilization practices. A carefully defined exponential function was fitted to the selected data. Function incorporated parameters that were estimated by non-linear regression analysis. The prompted leaching curve yields the marginal nitrate leaching rate which corresponds to the recommended applicable (N) rate.

For the prediction of soil bulk density after tillage an empirical equation used in Water Erosion Prediction Project (WEPP) model. The equation is based on the (EPIC) model [124,133]. Vegetative Filter Strip MODel (VFSMOD) as a coupled component of the SWAT model (see also 3.1) is also an empirical model that calculates sediment reduction [64].

3.3 Deterministic Models

The deterministic models are classified into static - time is not considered as a variable – and dynamic model [18]. A further subclassification is lumped or distributed models [134]. DRAINMOD is a deterministic hydrological model employed for agricultural subsurface drainage for nutrient transport, groundwater salinity problems. Variants of the model are applied for nitrogen subsoil fate DRAINMOD-NII, salinity penetration in soil cropping DRAINMOD-S and phosphorus concentration prediction and dynamic DRAINMOD-P [16, 19,135-138]. The model is well tested over river basin in Belgium for nitrates pollution prediction [139].

MIKE-11 is a 1-D hydrodynamic model in rivers, which simulates dynamic water movement in channels and rivers. It has been applied as water quality model in England and India. It employs time series of depth, pollutants and flow data along with water quality parameters [140].

Artificial Neural Network (ANN) science field combines bidirectional long short-term memory (BiLSTM) and adaptive neurofuzzy inference system (ANFIS) to predict water quality in different groundwater by employing single exponential smoothing (SES) as a preprocessing method to adjust the weight of the dataset input, (categorized as training and testing data) [141]. The internal operations of a neural network are deterministic, not stochastic, after training is finished.

3.4 Stochastic models

When data is presented using stochastic modeling, specific degrees of randomness or unpredictability are taken into account when predicting results. Using random variables, it predicts the likelihood of different outcomes under various scenarios.

SIMCAT it considered to be a model based on deterministic, stochastic, and Monte Carlo approach. It is used for the prediction of various biochemical indices and ions’ leachate i.e. BOD, DO, and Cl^-, NH₄⁺ NO₃^- [19,142].

Monte Carlo as a specific branch of stochastic modelling, is a broad class of computational techniques that rely on periodic random sampling to obtain numerical results. It uses randomness to solve problems that might be deterministic in principle [143].

Vadose zone flow and transport model (VADOFT) is a 1-D finite-element prediction code applied on chemical transport, flow/fate in the unsaturated zone. It employs parameters such as pressure, water content, and hydraulic conductivity to solve the flow equations. The code when equipped with Monte Carlo pre- and post-processor enables the run of multi-parameter scenarios several hundred times, and provide stochastic (probabilistic) outputs [144].

3.5 Artificial Intelligence and Machine Learning models

Artificial Intelligence (AI) is getting involved, the latest years, with forecasting groundwater quality modelling methods e.g. Artificial Neural Network (ANN), Evolutionary Algorithm (EA) etc. Given its capacity to handle enormous volumes of data, it offers an alternate method for handling complicated nonlinear hydrological modeling [145-146]. The proper water quality parameter selection is from every aspect crucial in order to attain the suitable AI models training. In every strategy actual, in situ, datasets are to be directly compared with those of AI method testing [147]. AI is advantageous over non-AI models by reducing the time spent for data sampling and AI proneness to rapid nonlinearity identification patterns as regards input and output data [12].

AI models combined with metaheuristic optimization techniques demonstrate higher reliability in terms of capturing the nonlinearity of water quality parameters [147]. General categorizations of AI methodology for soft water quality estimation include artificial neural network (ANN) modelling, fuzzy logic (FL) based models, support vector machines (SVM) models, machine learning (ML) and hybrid models.

A Machine Learning (ML) model utilizes potentially either deterministic or stochastic methods upon different sectors to obtain the proper results via prediction methodology. ML models such as Adaptive Neuro-Fuzzy Inference System (ANFIS), Deep Learning (DL) [148], evolutionary computing (EC), ensemble learning (EN), hybrid-modeling (HM) or even support vector machine (SVM), support vector regression (SVR), were introduced and developed for the water quality parameters prediction and thus groundwater quality classification for irrigation purposes [11-12,147-150]. ML models improve rapidly ground water level forecasting which is regarded as critical for any water management planning [13].

ANN exploits relationship between input and output datasets [147], is produced via neuronsis a computational model that excels at content addressable memory, machine learning, pattern recognition and optimization [151-152]. ANN uses several algorithms inter back-propagation neural network (BPNN), Levenberg-Marquardt back-propagation (LMBP) [12-13], Feed forward neural networks (FFNN) subcategorized into single layer perceptron (SLP) and Multi- layer perceptron (MLP) [147]. They are preferrable to predict spatial distribution parameters in aquifer remediation [12]. ANN modelling was applied in numerous cases over various regions e.g. in India, [153-154], for groundwater quality prediction. Deep Neural Network (DNN) within ANN modelling enables accurate estimate of nitrates and heavy metals i.e. cadmium, and chromium [155].

Because of their quicker optimization and capability for generalization, EA and SVM demonstrate higher predictive performance comparable with ANN and ANFIS [156-157], albeit limited studies were completed regarding groundwater quality forecast so as to enhance confidence about the proper effectiveness for complex hydrogeological systems simulation. A major drawback is considered to be the long training time (see Table 2). Nonetheless, in the near future, ANN methodology shall be improved and ANFIS technique, which combines ANN and fuzzy inference systems advantages shall play a major role. Furthermore, ANNs, types e.g. back-propagation neural network (BPNN), Levenberg-Marquardt back-propagation (LMBP) and multilayer perceptron (MLP) algorithms are the most preferrable in AI groundwater modelling due to their accuracy [12]. ML models in many cases are analyzing landfill leachate quantity/quality. Thus, e.g. hybrid Artificial Intelligence Model (AIM) based on grey wolf metaheuristic optimization algorithm and two stage extreme learning machine (ELM-GWO) is adopted to predict landfill leachate quality parameters in Iran. Other single-stage AIMs to be mentioned (MARS, MLPANN, ELM etc.), two-stage AIMs (MLPANN-GWO and ELMGWO etc.) [14].

FL based models make use of adaptive neuro fuzzy inference systems (ANFIS) which is an integration of fuzzy inference system (FIS) and adaptive neural networks. Hanoon et al., (2021) [147] and Haggerty et al., (2023) [148], presented analytical tables of AI models and adopted learning methods to support groundwater quality forecast. Ensemble fuzzy models combine fuzzy logic methods with ensemble learning for better performance outcome. They have been applied exclusively to improving DRASTIC method performance [148,158].

Recurrent Neural Networks (RNNs), integrated with geographic information system (GIS) technology enables scientists to predict accurately groundwater quality indices and cope with health risk management [159].

3.6 Physics-Based Models

Physics- based models are mostly one-dimensional leaching models i.e. RZWQM [160-161], WAVE [162-164] and DAISY [165] which simulate root zone processing. Numerous models are hybrid i.e. process- and physics-based model e.g. MODFLOW/MODFLOW-MT3D (see 3.1).

Daisy simulated crop production and nitrogen and water balance in the root zone [102]. It includes modules and soil-water dynamic based on Richards’ equation [78]. It incorporates nitrogen transformation i.e. mineralization/immobilization, nitrification and denitrification and agricultural management practices.

Numerous RZWQM components were developed to improve simulation accuracy regarding crop root zone. Therefore, RZWQM contains a soil water and a heat transfer module [166-167], a generic plant growth module [168], a nitrogen balance module [169], a soil equilibrium chemistry module [170], an evapotranspiration (PET) module [171], a management module [172] and a pesticide module [173].

The modeling system MIKE SHE describes the distribution of solutes and water flow in a catchment using physical principles. This entails the numerical solutions of the coupled partial differential equations for the processes of evapotranspiration and snowmelt for overland (2-D), channel (1-D), unsaturated (1-D), and saturated (3-D) flow. A full modeling system for simulating the movement of water and nitrate within a watershed can be created by integrating Daisy with MIKE SHE [102].

Phosphorus (P) is an important nutrient and though less abundant in surface water compared to nitrogen, still within a constant proportion N:P, e.g. 14.7:1 in the sea column [174]. Ratio unbalancing brings about eutrophication phenomena in aquatic ecosystems [175]. Therefore, the study of phosphorus fate mechanisms and selective infiltration routes in terrestrial and aquatic ecosystems is of high priority. Pferdmenges et al., (2020) [176], introduced a comprehensive review of phosphorus-soil interaction models, classified into multiple categories, inter alia, temporal and spatial scale, surface /subsurface transport, matrix/macropore transport, mobile/immobile transition in dual porosity models.

4. Metal Ions Leaching and Solute Transport

4.1 Soil Medium

Geochemical PHREEQC code [177] performs sufficient numerical simulations incorporating mixing transport models both in aqueous and gaseous phase, taking into account geochemical reactions, such as precipitation, dissolution, complexation phenomena involving metal oxides and microbial activity. The reactive transport model is incorporating reaction kinetics for chemical processes including ion-association and specific ion interaction theory for solute activities calculations [176]. In a publication was employed to calculate Saturation Indices (SI) [178]. The simulation code was subsidized by predicting/estimating tools of model parameters such as Parameter ESTimation (PEST) software [179]. Sorption phenomena simulated by employing two-layer model introduced by [180]. Organic matter sorption was modeled by implementing WHAM model an acronym of Windermere Humic Aqueous Model, [181] which is based upon humic acids behavior according to Lewis’ theory for electron donators. Mine tailings and heavy metal leaching e.g. Pb, were encountered by activating reclamation measurements to mitigate negative leaching effects. Mining slurry combined with manure leads to metal sorption reducing any leaching processing. A new biogeochemical model was developed taking into account certain phenomena i.e. kinetically-controlled dissolution/precipitation, adsorption, water gas exchanges, surface complexation reactions, microbial respiration and growth, dissolution and precipitation. The amended model predicts more accurately Pb reactivity and soil tailings [182].

Heavily contaminated soil (brownfield) by Potentially Toxic Elements (PTE)s are a permanent threat to the local ecology via toxicity and bioaccumulation. Leaching, and runoffs phenomena are putting in peril the human wellbeing. Threats are getting even bigger considering pH change of the soil cap. Therefore, chemically burdened sites have to be under constant monitoring in terms of heavy metals quantification and ion mobility checking indices, to assure that surface watersheds and aquifers are not threatened. It is acknowledged the lack of detailed knowledge of minerals which exert control over toxic element leaching. Thus, is of great importance the investigation of soil mineralogy in PTEs release and leaching prediction [183]. The geochemical modeling program PHREEQC (version 3.1) was deployed to calculate Saturation Indices (SI), which indicate whether a specific mineral is likely to dissolve or precipitate in groundwater, in a case study in Pakistan of monitoring groundwater area suffering from Arsenic contamination [178]. Geochemical simulations were conducted making use of PHREEQC software [177], MINTEQ’s two thermodynamic databases [184] and Lawrence Livermore National Laboratory (LLNL).

Sampling investigation conducted from a brownfield site served for fertilizers production. Jarosite rich in Pb, hematite, and gypsum were the most abundant mineralogical phases/stratifications, with zinc sulfate, kintoreite and anglesite. Pb and Zn phases were identified as the dominant ones [183]. pH-dependent leaching tests were applied incl. HNO₃ and NaOH as the pH control agents in combination with geochemical models. Testing target was the reveal of leaching mechanisms and contaminants solubility in a pH ranged from 1 up to 12. All information gained advances scientific knowledge and facilitates future remediation strategies and promotes integrated sustainable management and strategy design over brownfields [183].

Barren piles of mining undergo biooxidation giving rise to metal sulfides oxidation and finally acidification phenomena which motivate ions mobility, well kwon in scientific community as Acid Mine Drainage (AMD). AMD leachates’ mitigation exert great pressure to bedrocks, surface waters and local aquifers since metal ions’ active plume put in risk the overall ecological quality of the region in vicinity to the mines.

Bio geochemical AMD processing could be regulated by applying in coal mine tailings controlled oxidative agents i.e. ozone and hydrogen peroxide, so as biooxidation to be accelerated many times as much and therefore to mitigate, in relative short term, mobilization/leaching problems by achieving the extinction of AMD processing causative sources. Selected agents demonstrate high oxidation potential and no harmful residues are to be formed when they undergo decomposition during treatment [185]. Subsurface propagation of the oxidative agents employed and the active plume formed were monitored by using well set up mathematical models [186].

A conservation equation was used in a porous medium, performing variable saturation along with convection i.e. diffusion and dispersion coefficients of liquid gas phases. Brinkman equations were adopted for porous medium flow. Finite element simulations were achieved by using COMSOL Multiphysics® for the final solution. The permeability was determined using the well tested Kozeny–Carman model [187] and the dispersion by coupling Richards’ equation [100] and van Genuchten [80] retention model. Propagation models applied in the subsurface of coal mine tailing piles are based on the kinetic rate estimated of the agents’ consumption and bacterial activity acting as catalyst. Abiotic pyrite oxidation is correlated with the surrounding redox potential (E_h) [188].

Hydrogen (H₂) is a versatile and carbon–neutral energy fuel produced by various methods, inter alia, electrolysis, gasification, Steam Methane Reforming (SMR), etc. and very promising for future energy demand cover. Nonetheless, certain conventional storage barriers render geo-storage alternative as an attractive option since it assures high storage capacity, comparably low-cost investment & exploitation viability, enhanced safety, high energy storage density etc. Salt underground caverns, which serve as natural reservoirs, seem to be a fine selection incorporating many significant advantages. Suitable geo formations of such kind are located worldwide facilitating safer energy storage planning.

Li et al., (2020b) [189] introduced an upwards dissolution rate model, simulating H₂ geo-storage in salt caverns. It entails finite volume and structural dynamic elastic mesh methodology via C++ programming. Furthermore, Wang et al., (2021) [190] developed a new prediction and optimization leaching parameter model of Solution Mining Under Gas (SMUG). Certain parameters i.e. nitrogen volume gas-brine interface depth, leaching rate, water and nitrogen injection pressure and leaching rate possess crucial role for the final accurate prediction outcome.

AbuAisha et al. (2019), [191] modelled H₂ migration to surrounding rock domain my applying Darcian percolation and Fickian diffusion concepts. The model coupled thermodynamic and hydraulic driven transport mechanisms and drew the conclusion that van Genuchten’s model parameters were overestimated and H₂ mass loss towards rocky formations was proved to be insignificant given certain pre assumptions [192].

4.2 Modified Soil Medium

PHREEQC computer code program [193], simulates leaching curves of experimental modified soil columns (see also 4.1). Sewage sludge amended to experimental calcareous soil The same soil specimen tested was enriched with injected heavy metals and nanoparticles of ZnO and MgO and zeolite functioning as absorbents. Simulation program inputs comprise four modules i.e. SOLUTION, EQUILIBRIUM PHASES, EXCHANGE, and SURFACE The columns of the enriched soil incubated for 7 days and leaching experiments that lasted 18 days were conducted by employing two fixed leaching solutions made of CaCl₂ and diethylenetriaminepentaacetic acid (DTPA) in predefined conc. PHREEQC simulation program employs organic surface complexes. Heavy metals’ Log K_s extracted from SHM and NICADONNAN databases [181]. Modified soil columns responses were studied against leaching solution and more specifically the leaching tailing of cadmium (Cd), nickel (Ni), copper (Cu) and zinc (Zn). SOLUTION module received combinations of mixed CaCl₂ and DTPA solution of determined conc. equilibrated with solid phases received from MINTEQ. The latter is a software program which simulates heavy metal transfer on account of the preprepared mixture solutions and solid phase interaction. Heavy metal propagation was simulated by means of absorption/complexation of MgO, ZnO, and zeolite serving all as reactive sites. Ion exchange sites simulated by the EXCHANGE module and heavy metal transport by TRANSPORT component. Monitoring cell (individual grid) and leaching time were properly defined for the running simulation [194].

As a result, PHREEQC displayed a respectable capacity to simulate the leaching of heavy metals through modified soil columns well prepared to undergo lab scale testing, though in earlier years [195] there have been experimental mismatches against PHREEQC predictions. Input data extracted from data bases are referring to pure solid phases which is not our case in the conducted experiments. Therefore, mismatches are to great extent justified since leaching rate by simulation was highly dependent on the Saturation Index (SI) and the solid phase type. Use of different soil type entails deeper study of surface reactions and ion binding.

In bibliography, PHREEQC code simulates the experimental results from heavy metals leachates derived from various waste forms and sewage sludge [196-198]. It is suggested that the cation exchange reaction is the predominant mechanism of Cd transport in the soil profile. PHREEQC predicts well heavy metal simulation Pb leaching in modified soil taking into consideration ion exchange mechanisms and surface complexation reactions [199].

Electric Power Research Institute (EPRI), an independent non-profit energy research, development, and deployment organization adopted MINTEQ approach model to develop FASTCHEM geochemical hydrodynamic solute transport model for modelling coal ash leachate dispersion and advection. It demonstrates useful appliances over heavily polluted soils by fly ash dispersion derived from fossil fuel powerplants operation [200].

RZWQM2 (1-D) point scale model was tested over cultivated soils and soils enriched with slurry mostly sandy and sandy loam quantities. In all cases studied, predictions were made for soil mineralization, ammonia gas and N₂O emissions and nitrates leaching regarding each individual soil type testing. The scope of the experimental effort was to set up integrated strategies over nitrogen gases abatement derived from livestock farming activities and the prevention of nitrates leaching phenomena on the ground [74].

The application of agri-residuals over cultivated soil and slopping farms is a common practice to assure higher crop yield, to minimize erosion phenomena and to increase water soil infiltration and retention along with the improved evaporation characteristics. Wheat straw mulching or rapeseed-oil residue improves rainfall infiltration outcome. Experiments conducted by employing classic Horton, Kostiakov and Philip infiltration models combined with evaporation models by Rose and Gardner denoted that water retention results regarding mixed soils with agri-additives were significantly improved. Effective water-retention techniques may be developed based on how the direct straw integration types correspond to soil infiltration and evaporation. Experimental research demonstrates that mixed straw mulching improves soil water response by increasing infiltration, improving sediment yield and by incurring runoff reduction. The former is especially helpful for implementing an integrated ecological management over agricultural regions, arid/semi-arid with limited water resources, like China [201-202].

4.3 Cement Leaching Medium

Zielina et al., (2022) [203] proposed a mathematical model to predict heavy metal leaching through fresh mortar lined porous cement (Portland cement, PC). Leaching is taking place into water at dynamic conditions, taking into account sorption phenomena apart from the already applied in literature i.e. dissolution, diffusive and advective transport. The safe prediction of heavy metals concentration at the ends of drinking water supply pipeline sections will meet all necessary actions for the suitable rehabilitation by cementation in accordance with the European active legislation [204], conc. water quality intended for human consumption [203].

(PC)s basic constituents are oxides, mainly CaO, SiO₂, Al₂O₃, Fe₂O₃ etc. Certain heavy metals i.e. Cr, As, Cd, Pb are inevitably part of the final cement production, derived from various sources, which potentially affect negatively, when accumulated, due to their high toxicity to the human health and ecosystems. (PC) underwent coupling during hydration and leaching tailing kinetics of the four above given metals and metalloids were tested over different (PC) specimens [205]. Earlier research conducted of iron leaching from steel or cast-iron model, dependent on physicochemical parameters [206], and copper leaching via hydrodynamic processes incl. the presence of biofilm [207].

Heavy metal kinetics equation is based on Elovich empirical equation. It was selected and tested, slightly modified, since it can reflect adsorption and desorption processes of heterogeneous chemical reactions. Constants of the equation reflect the diffusion rate of heavy metals from solid to liquid and the rapidness of the diffusion rate. Furthermore, the selected Freundlich dynamic equation simulates adsorption and desorption kinetics of ions and heavy metals in soil, and deduces correlation coefficient of adsorption process. Last Finally, the 2^nd order kinetics equation was also preferrable over others of the same kind since it describes satisfactorily the reaction rate and ions concentration in equilibrium state and the parabolic equation expression depicts the ion release process, controlled by multiple diffusion mechanisms [205].

4.4 Landfill Capping Layers

Numerical modeling techniques predict efficiently leachate generation [208] especially determining parameters of leachate collection, drainage systems etc. [209]. Hydraulic Evaluation of Landfill Performance (HELP) [210] is a deterministic, widely used, quasi-two-dimensional hydrologic model. It combines (1-D) soil physical and hydrological processes, both in saturated and unsaturated vertical flow and towards lateral direction (i.e. lateral drainage). Model run requirements are certain data to be collected i.e. saturated hydraulic conductivity, soil water retention characteristics, actual meteorological data, solar radiation, leaf area index, evapotranspiration, surface runoff and the interflow [211]. The model simulates water vertical flow through up to four types of landfill layers [77].

Other models evaluating hydrologically landfills are UNSAT-H [212], HYDRUS-1D [213], Finite Element subsurface FLOW system program (FEFLOW) [214]. It also incorporates evaporation, infiltration, runoffs and lateral drainage [215]. Input data for the model running are climatic (i.e. temperature, wet precipitation, solar radiation), soil characteristics (i.e. porosity, hydraulic conductivity, field capacity, wilting point), layers arrangement (i.e. surface snowmelt, surface runoff, interception or rainfall by vegetation, evaporation), vegetation. The model comprise surface (i.e. snowmelt, rainfall, interception of rainfall by vegetation, surface runoff, evaporation) and subsurface processing (i.e. soil-water evaporation, vertical drainage, plant transpiration, lateral drainage, liner leakage) [77]. The soil profile was divided into smaller profiles in order to facilitate computation.

Leachate prediction of landfill capping soil profile can be achieved by employing Robust Integrated Artificial Intelligence combined with the use of (GreyWolf) Metaheuristic Optimization Algorithm and the support of Extreme Learning Machine [14], (see also 3.5) for predicting landfill leachate quality (COD and BOD₅) and groundwater quality (turbidity and EC) at the Saravan landfill, in Rasht region, Iran.

Models and countries that were implemented in the past, according to the relative literature are presented in Table 3.

'Fertigation' is the technique of supplying dissolved fertiliser to crops through an irrigation system. When combined with an efficient irrigation system both nutrients and water can be manipulated and managed to obtain the maximum possible yield of marketable production from a given quantity of these inputs.

5. Discussion

Great research conducted the last years, focused on proactive and rehab actions to safeguard surface and subsurface soft water quality and its uses. Med river types are suffering, all year round, due to climate change (prolong drought periods, rare sudden cataclysmic rainfalls) from low waterflows and moderate to low water quality level from bio-physicochemical criteria standpoint. Intensively cultivated areas are adding huge quantities of fertilizers, mostly biowaste from breeding farms which gives rise to nutrients conc. in surface water basins and watercourses. Aquifers progressively are being nitrified on account of the nitrogen-based fertilizers’ surplus, applied on the topsoil that turn into soluble nitrates and end up to the underlying water table. Apart from nitrates limitation conc. where water is intended for human consumption, pumping up and use for irrigation purposes incur unbalance to the humic part of the topsoil and affects negatively corps yield due to increased salinity.

The alternative use of soil additives enhances crop yield, mitigates nitrate pollution, facilitates remediation methods, improve anion exchange, soil water retention and soil texture characteristics. Soil additives have their own mechanisms i.e. promote fertilizers release in a controlled way to meet metabolic needs of plants [283-284], natural derived or artificially manufactured adsorbents/absorbents, e.g. biochars, zeolites peat, hydrogels etc. [285-286]. There have been literature reports over decision support tools e.g. Multi-criteria Decision Analysis (MCDA) to have been employed in order to couple LCA with nitrate infiltration models nitrates soil management [237,287].

Therefore, well tested nitrate leaching models standalone or as a part of an integral modelling package provide rapid site-specific estimates of N-NO₃ leaching potential below root zone of crops and the impact of nitrate leaching into groundwater. Models developed within sophisticated integrated model packaging are estimating not only the residual nitrate part but also all potential nitrate pool and leaching losses from small plots up to immense catchments of great acreage. They take into consideration cropland soils into subsoil layers, under different combinations of meteorological data, soil profile properties and farming practices.

Strict cropland’s, nitrogen-based fertilization regulations in (kg N ha^-1), in Europe, compel farmers to resort to alternative cultivating options and lie emphasis on cultivation models’ prediction that change bio-fertilizing practices and support the terrestrial nitrogen balanced cycle. (GIS) component offer spatial soil database and increases predictability in demanding terrains which demonstrate individual soil characteristics.

Mechanistic models, to a great extent, are employed in scientific literature to predict quantity and quality of waste rock drainage and the incurred environmental risk. Nonetheless they demonstrate great uncertainty due to apparent heterogeneity (e.g. texture, hydraulic, geochemical etc.) the waste deposits which entails a large number of scale-dependent parameters to support the model equations. To address such cases combined deterministic and stochastic modeling approach is considered to be crucial to deal with layers’ heterogeneity [80].

Many researchers promote the development of nested catchment models, since regional applied scale models are inevitably dependent to smaller scale models such as those that provide valuable data to boundary conditions. Furthermore, Intergrated hydrologic models, i.e. models that eliminate boundaries definition between surface and subsurface, are rarely applied on regional scale due to evident computational cost which affects negatively the model calibration and incurs high predictive uncertainty. Moreover, estimation of baseflow and each water budget component usually requires additional post-processing steps [245].

Considering the effectiveness of ANN applications and AI offspring are marginal advantageous in predicting hydrochemical and hydrogeological parameters, and most specifically, extensively applied to groundwater modelling and prediction. ANN water prediction tools, could minimize water quality monitoring stations and propose alternative ways to estimate groundwater level not relying upon geoelectric characteristics [141]. Machine Learning models make use of water quality indices for the final water quality assessment [150].

Μeta-research algorithms e.g. Latent Dirichlet Allocation (LDA) are employed lately to scrutinize relevant literature topics and boost quantitative analysis in nutrients water soil transport processes. Thus, are valuable towards a better identification of the trends and to fill knowledge gaps [251].

6. Conclusions

Nitrogen leachate prediction models yield satisfactory results as a multi-component tool that co-processes data acquired from site characteristics (e.g., soil texture and depth, terrain elevation, meteorological data (e.g. precipitation, temperature, etc.) in long time series. Microbial communities of the soil and eco-physiological parameters (e.g. canopy light level, maximum photosynthetic activity) are taking into deeper consideration even more, as critical paths to nutrients uptake in the seeded plants’ root zone which regulate nitrogen surplus to be leached to deeper soil substrates. Soil water simulation modelling are valuable tools for addressing soil nitrogen pollution and contributing soil nutrient decision strategies as regards crop and fertigation management e.g. APSIM in Australia. Nitrate – water simulation models in soil profiles require a great number of input data some of which with high accuracy, measured at the field scale, a situation not always feasible. Thus, dynamic approach model systems were developed to overcome complex processes resulting from numerous interactions.

The vast majority of the models examined were process-based or conceptual models. Nonetheless, empirical prediction models, though rather simplistic and therefore not preferable, demonstrate certain advantages since they do not require extensive calibration database information which in many cases is unavailable. Artificial Intelligence and machine learning techniques extent even more prediction capability offering a great potential of accurate water quality prediction without time series use acquired from a long period measurement and databases. Hybrid structural chemical agents’ prediction models facilitate researchers to encompass even more profound sophisticated tools regarding climatic change scenarios, greenhouse gases in atmosphere, combined with soil surface temperature and proceed to climatic analyses, synthetic weather predictions and support more accurately crop assessment yields, counter measures decision making against extreme weather conditions, secure crop yield and deeper understanding of crop physiology.

Intense anthropogenic pressure and climate change with prolonged dry spells, incur saline wedge penetration phenomena to the coastal low water table resulting in water quality devaluation and soft water consumption limitation. Coastline water management involves Integrated Modelling Systems which are coupling successfully surface and groundwater hydrology, reservoir operation, agronomic/crop planning, and nitrates leaching to address aquifer’s nitrates transport and seawater intrusion and resilient adaptation plans.

Stochastic modeling is widely used as a fundamental approach tool to embed uncertainty in long-term model-based decisions. Therefore, such a methodology should operate like a regulator against policy and decision makers who should take into serious consideration stochastic models results.