Quoted from: Zipper, Samuel C., Mehmet Evren Soylu, Christopher J. Kucharik, and Steven P. Loheide II. "Quantifying indirect groundwater-mediated effects of urbanization on agroecosystem productivity using MODFLOW-AgroIBIS (MAGI), a complete critical zone model." Ecological Modelling 359 (2017): 201-219. https://doi.org/10.1016/j.ecolmodel.2017.06.002 

        To provide a new tool for addressing interactions and feedbacks between groundwater and ecosystem processes with the explicit representation of dynamic vegetation growth, we coupled the AgroIBIS ecosystem model to the modular groundwater flow model, MODFLOW-2005 (Harbaugh, 2005); hereafter, we refer to the coupled models as MODFLOW-AgroIBIS (MAGI). MODFLOW-2005 (Harbaugh, 2005), is a version of the MODFLOW family of finite-difference groundwater flow models, which have been widely used and validated since being introduced in the 1980s (McDonald and Harbaugh, 1988). MODFLOW has been successfully applied at scales ranging from riparian floodplains (Doble et al., 2006, Faulkner et al., 2012, Mastrocicco et al., 2014, Zell et al., 2015) to watersheds (Cho et al., 2009, Falke et al., 2011, Lam et al., 2011, Sutanudjaja et al., 2011, Sutanudjaja et al., 2014) to global (de Graaf et al., 2015, de Graaf et al., 2017).

        MODFLOW simulates three-dimensional saturated groundwater flow for a porous medium using a block-centered finite-difference approach and solving the partial-differential equation:

        where Kxx, Kyy and Kzz are values of hydraulic conductivity in the x, y, and z directions [L T−1]; h is the potentiometric head [L]; W is a volumetric flux per unit volume representing sources and/or sinks of water, e.g. pumping wells [T−1]; S is the storage coefficient accounting for elasticity of the porous material [L−1]; and t is time [T] (Harbaugh, 2005). MODFLOW is designed and built around the concept of modular programming, in which specific packages simulate different components of groundwater flow systems. Existing packages simulate processes including solute transport (Guo and Langevin, 2002, Hornberger et al., 2002, Zheng et al., 2001), groundwater management (Ahlfeld et al., 2005), aquifer subsidence (Hoffmann et al., 2003), and interactions with surface water features (Merritt and Konikow, 2000, Niswonger and Prudic, 2005, Prudic et al., 2004).

        To couple AgroIBIS and MODFLOW, we converted AgroIBIS to a MODFLOW package. Similar to the HYDRUS package for MODFLOW (Seo et al., 2007), AgroIBIS simulates all land surface and unsaturated zone processes, and calculates the average flux across the bottom of the soil profile for each MODFLOW timestep. This flux is passed to MODFLOW as a source/sink (W in Eq. (2)) to the upper layer of the groundwater flow equation. MODFLOW then solves the groundwater flow equation to calculate a new hydraulic head distribution over the model domain, which is provided to AgroIBIS as a pressure head bottom boundary condition corresponding to the height of the water table above the bottom of the AgroIBIS soil profile. Within this coupled framework, all changes in storage due to the draining and filling of pore spaces occur within the AgroIBIS soil column via the response of the pressure head profile to changes in the bottom boundary condition and corresponding changes in soil moisture; within MODFLOW, only storage changes resulting from aquifer and fluid compression/expansion are represented. Thus, all AgroIBIS cells must be discretized such that hydraulic head never drops below the bottom of the soil profile during a simulation.

        MAGI simulates myriad processes within a single modeling framework and these processes operate at timesteps ranging from minutes to years. “Fast” processes (e.g. near-surface fluxes of energy, water and momentum, photosynthesis, and respiration) operate at timesteps of minutes to hours; “medium” processes (e.g. vegetation phenology and carbon/nutrient cycling) operate at timesteps of days to weeks; and “slow” processes (e.g. agricultural management, planting decisions, and land cover transitions) operate at monthly to yearly timesteps. In MAGI, these processes are organized in a hierarchical framework allowing each process to operate internally at the proper timestep, while the coupling between AgroIBIS and MODFLOW is designed in a way such that the groundwater flow equation can be solved at a user-defined frequency as determined by the needs of a particular simulation. For example, in simulations where lateral flow is expected to change slowly and/or groundwater is expected to be relatively unimportant to ecosystem processes, the groundwater flow equation can be solved at weekly or longer timesteps, while in applications where hydraulic conductivity is very high or lateral groundwater flow is very important, MODFLOW and AgroIBIS can exchange fluxes hourly.

        One issue associated with differences in timestep length between AgroIBIS and MODFLOW is that holding a constant pressure head bottom boundary condition in AgroIBIS for an entire MODFLOW timestep can lead to overestimates of groundwater recharge or capillary rise, causing oscillatory groundwater table dynamics over time. To eliminate these issues, we implemented an intermediate prediction mechanism into the coupling between AgroIBIS and MODFLOW, shown in Fig. 1b. This prediction mechanism estimates a new AgroIBIS bottom boundary condition based on the storage coefficient of the uppermost MODFLOW layer and the calculated bottom boundary flux each time the Richards Equation is solved, and which is used for the next Richards Equation timestep, and this process is completed until a new MODFLOW timestep is reached. At this point, the cumulative bottom boundary flux from AgroIBIS over the entire MODFLOW timestep is input into the MODFLOW groundwater flow equation, and the AgroIBIS bottom boundary condition is updated based on the new MODFLOW solution. Effectively, this prediction step can be thought of as approximating the groundwater response as one-dimensional in the vertical direction during each AgroIBIS timestep, with lateral flow occurring only during each MODFLOW timestep.