Artificial drainage has been known to be an important water management practice for farming of the most productive soils of the Midwest  . Draining water from the soil profile is an important hydrologic component in most agricultural soils  . Artificial drainage may be provided by installing drainage ditches or drain tubes. These systems are usually installed in irrigated arid and semi arid lands to control water logging and salinity. The successful performance of a drainage system depends on optimal design of drain depth and drain space  . The need to have guidelines for drainage design and water management for different soils and climates has driven both the experimental field research and computer modeling  . Computer-based simulation models can predict subsurface drain flows, water table fluctuations, and crop yields in a greater variety of conditions than what is feasible through monitoring, which allows timely decisions to be made about complex problems when field data are both difficult and expensive to obtain   . System dynamics is a methodology developed by Forrester to analyses dynamic behavior of complex systems containing biological, economic, social, technological and political elements, aided by computer  . In this study, the system dynamics technique was used to simulate the performance of a drainage system under wheat crop in a clay soil. The system dynamics tool, Stella is used to provide a fully integrated simulation system to conceptualize, document, simulate and analyze subsurface drainage water system. The objective of this paper is to develop a system dynamic model to predict drain discharge behavior for different drain systems (depth and spacing), under wheat crop in a clay soil in Mashtul Pilot Area (MPA), Egypt. The results can be used to design the drainage system geometry for better water management on both field and catchment scales.
2. Mashtul Pilot Area as a Case Study
The field work (sampling and measurements) was carried out in MPA. MPA was constructed in 1980 in south-eastern part of the Nile Delta  . It is situated 90 km northeast of Cairo in a rather flat area (Figure 1). The area is approximately 260 feddans. Table 1 showed irrigation schedule & amount and crop productivity in MPA, 2005. The layout and design of this pilot area were already planned by the DRI of the National Water Research Center (NWRC), as shown in Figure 1.
Figure 1. Layout of Mashtul pilot area in the Nile Delta.
Table 1. Irrigation schedule & amount and crop productivity in MPA, 2005.
The southern and western boundaries are formed by the Mahmoudia Drain and its branch; the northern and eastern are bound by tertiary irrigation canals. It is characterized by a deep clay top layer and a sandy aquifer. The clay layer, which is approximately 6.0 m thick, contains about 35% silt and 65% clay. Irrigation water is delivered by gravity to the tertiary canals and lifted approximately 0.5 m to field level by pumps.
The area is drained through a subsurface drainage system that consists of parallel PVC lateral drains, which discharge into buried concrete collector drains through a manhole. The design of the subsurface drainage system was made according to the standard criteria of the Egyptian Public Authority for Drainage Projects (EPADP). Two types of criteria can be distinguished; namely, the agricultural and the technical ones. The agricultural criteria are an average depth of the groundwater table midway between the drains of 1.0 m and an average drainage rate of 1.0 mm/day to permit sufficient leaching. The technical criteria are a design discharge rate for the determination of drain pipe capacity of 4 mm/day for rice areas and 3 mm/day for other crops, a safety factor of 25% in the design of the collector drains to account for sedimentation and misalignment and change in diameters and maximum depth of 1.5 m for laterals and 2.5 m for collectors. The area was divided into eighteen drainage units with different drain depths and spacing. The units were cultivated with a single crop for each unit during each cropping season. Berseem (Egyptian clover) and wheat were cultivated as winter crops and cotton, rice, and maize as summer crops  .
MPA was controlled under the existing farming conditions, aiming to apply the data measurement program for two years, from June 2005 to May 2007, and consequently two winter and two summer seasons. The area was surveyed with a Global Positioning System (GPS) to locate its boundaries, and with a Geographic Information System (GIS) for sampling locations, as shown in Figure 2. The sampling locations included the manholes for field lateral drainage water sampling and collectors’ outlets for collector drainage water sampling. The measurement program was applied in the area for the four seasons as follows;
Determination crop pattern, fertilizer amount, and time required to apply each crop for each unit.
Collected drainage water samples before cultivation, before and after applying fertilizers, and periodically every 10 days.
3. Method and Materials
3.1. Model Development: Conceptual Model
In order to quantifying analyze the subsurface drainage system in MPA, a dynamic modeling of the system is made (Figure 3). Stella software is used to develop
Figure 2. Sampling locations in MPA.
Figure 3. Overview of the main processes included in the conceptual model to calculate the drainage and related water table management.
the model. Furthermore, the model calculates the future drainage water flow. Solving the water and solute transport equations requires two soil water relations, namely the soil water content-water potential relation and the soil water potential-hydraulic conductivity relation. They were taken according to  . Values of several coefficients of van Genuchten equation such as “n”, “m” constants were calculated based on the relationship between “n”, “m” and the pore size distribution index, “l” (e.g. n = l + 1, m = l/n). The later (l) is available in literature and varies according to the soil textural class. “l”, values are given in the soil data base and can be edited by the users. Bubbling pressure, soil water content at saturation, field capacity and wilting point, and saturated hydraulic conductivity were given in data base for different soil types and can be edited by users should filed data become available. The data base information was collected from different sources worldwide and reference has been made to the sources.
3.2. Modelling System
The conceptual model was applied on the study site to model the water balance predictions of the subsurface drainage which is employed to simulate the performance of drainage and related water management systems.
The conceptual model uses dynamic modeling to quantify subsurface drainage, deep seepage, infiltration, and evapotranspiration. Subsurface drainage flux is calculated based on the assumption that mainly lateral water movement occurs in the saturated region. The flux is determined by the water table elevation at the mid plane of the drains and the water level in the drains. The Hooghoudt’s steady-state equation  is used when the depth of water on the surface is less than the surface storage (S1):
where qH is the water flux [L T−1], mdr is the midpoint water table height above the drain [L], Ke is the effective lateral hydraulic conductivity [L T−1] and Ldr is the horizontal distance between drains [L]. de is the equivalent depth from the drain to the restrictive layer [L] and is specified in the model to correct for convergence flux near the drains. It can be defined as a function of Ldr, the drainage tube radio rdr [L] and d [L], the real depth from the drain to the restrictive layer . crf is the ratio of the average flux between the drains to the flux midway between the drains and is controlled by the shape of the water table profile. crf is approximated as 1.0, which implies that Equation (1) corresponds to the ellipse equation that is often used to determine drain spacing  . Because the soil is not a homogeneous media, Ke is defined as:
where N is the number of soil (horizontal) layers, (KH)i is the horizontal hydraulic conductivity of the i-th horizontal layer and Ei is the thickness of the i-th layer [L]. This equation should be used observing the convention that soil layers are numbered progressively from up to down. If the water table is located in soil layer P (with 1 ≤ P ≤ N), then dP is the distance between the midpoint water table elevation and the bottom of layer P.
In the event that the depth of water on the surface exceeds S1 (the surface storage), the assumption of a curved water table completely below the soil surface fails (Bouwer and van Schilfgaarde, 1963) and in this case, the flux is calculated as:
where zpn is the ponded depth [L], zdr is the distance from the soil surface to drain [L] and cpn [--] is a constant defined for a given depth of soil profile and drain-size, -depth and -spacing (Skaggs, 1981). The deep vertical seepage flux is computed using Darcy’s law to calculate the flux through the restrictive layer:
where qv is the deep vertical seepage flux [L T−1], KV is the effective vertical hydraulic conductivity of the restrictive layer [L T−1], h1 is the average distance from the bottom of the restrictive layer to the water table [L], h2 is the hydraulic head in the groundwater aquifer referenced to the bottom of the restrictive layer [L], and Erest is the thickness of the restrictive layer [L]. Infiltration is computed through a modified Green-Ampt procedure. The potential evapotranspiration (PET) is estimated using the Thornthwaite equation  .
4. Results and Discussion
Despite the efforts of several alternative approaches to manage intangible factors, none has been sufficient to fully incorporate relationships between variables, delays and feedback, all of which characterize the behavior of intangible resources. So, Managers continue taking decisions only based (or support) on their experience, knowledge that constitute their mental models   . Therefore, there is a need to explore new tools to represent the complex relationships found in systems. One promising option is system dynamics which is a feedback-based, object-oriented approach. Although system dynamics is not a novel approach, it offers a new way of modeling for future dynamics of complex systems. According to Simonovic and Fahmy  , system dynamics is based on a theory of system structure and a set of tools for representing complex systems and analyzing their dynamic behavior. The most important feature of system dynamics is that it helps to elucidate the endogenous structure of the system under consideration, and demonstrate how different elements of a system actually relate to each other. This facilitates experimentation as relations within the system are changed to reflect different decisions   . Agricultural systems and their environmental effects, like many other environmental problems, constitute complex systems, which study requires systemic approaches capable of explicitly managing the temporal dimension, sustainability conditions, uncertainty and externalities   . Therefore, the system dynamic is good approach to model this system.
4.1. Causal Loop Diagram
Causal loop diagram is an important tool for representing the feedback structure of systems. A causal diagram consists of variables connected by arrows denoting the causal influences among the variables. A feedback loop is a succession of causes and effects such that a change in a given variable travels around the loop and comes back to affect the same variable. If an initial increase in a variable in a feedback loop eventually results in an increasing effect on the same variable, then, the feedback loop is identified as a “reinforcing or positive” feedback loop. If an initial increase in a variable eventually results in a decreasing effect on the same variable, then the feedback loop is identified as a negative, counteracting or balancing’ loop   .
The causal loop diagram in this study has shown in Figure 4. The first negative feedback loop represents the evapotranspiration effect: the larger the evapotranspiration, the less the soil water content and soil moisture stress “ket”, which in turn decreases evapotranspiration. The second feedback loop represents the interaction between evapotranspiration and upward flux: the larger the evapotranspiration, the larger the upward flux, then the larger the soil water content and ket, which in turn increases evapotranspiration. The third feedback loop represents the interaction between soil water storage and percolation: the larger the storage, the larger the hydraulic conductivity, then the larger the percolation, which in turn decreases soil water storage. The fourth feedback loop represents interaction between water table and soil water storage: an increase in the percolation increases water table and upward flux and soil water content. In the fifth feedback loop as water table rises by deep percolation, the depth of water above the drain increases which increases the drain discharge, and in turn decreases the water table.
4.2. Model Calibration
The model was calibrated by considering the time series of observed subsurface flux at the outlet of the drainage system (Figure 5). The calibration was carried out through a trial and error approach by which the parameters were varied manually on the basis of the assessment of model performance. Model performance evaluation was based on visual inspection of calibration plots (hydrographs of observed and predicted time series, hydrographs of observed and predicted cumulative time series, scatter plots of observed and predicted time series, etc.) and on the statistical assessment of model performance indexes such as; Mean Absolute Error (MAE), Relative Root Mean Square Error (RRMSE), Model Efficiency (EF), Coefficient of Residual Mass (CRM), Coefficient of Determination (CD) and, Goodness of Fit (R2). The characteristic of the different statistical criteria is given in Table 2.
Figure 4. Causal loop diagram.
Figure 5. Measured and simulated subsurface drainage water under wheat crop in a MPA using Stella dynamic modelling.
Table 2. The characteristic of the different statistical criteria.
The measured data and simulated results were compared in terms of subsurface drainage discharge (Figure 6). The statistical analysis results of water table are shown in Table 3. The statistical analysis implies a good fit between measured and simulated values. The comparative study reveals that the model performs well, reliable and accurate for predicting subsurface drainage water in clay soils.
Figure 6. Measured and simulated subsurface drainage water in a MPA.
Table 3. Statistical analysis results of water table in MPA.
MAE: mean absolute error; RRMSE: relative root mean square error; CD: coefficient of determination; EF: model efficiency; CRM: coefficient of residual mass; R2: goodness of fit.
A dynamic model was developed to predict drain discharge behavior for different drain systems (depth and spacing). The developed model was used to simulate the daily drainage water at the midpoint of drain spacing in a clay soil. The measured data and simulated results were compared in terms of subsurface drainage discharge. The statistical analysis implies a good fit between measured and simulated values. The comparative study reveals that the model performs well, reliable and accurate for predicting subsurface drainage water in clay soils. Results (predictions) can be used to design the drainage system geometry for better water management on both field and catchment scales. Results indicated that, the model can be used as a decision support tool to help policy makers in long term strategic management for irrigation projects. The developed model can potentially help in setting guidelines for using subsurface drainage water in agricultural sector in Egypt.