Groundwater Modeling in Support of Water Resources Management and Planning under Complex Climate , Regulatory , and Economic Stresses

Groundwater is an important resource that meets part or all of the water demand in many developed basins. Since it is an integral part of the hydrologic cycle, management of groundwater resources must consider not only the management of surface flows but also the variability in climate. In addition, agricultural and urban activities both affect the availability of water resources and are affected by it. Arguably, the Central Valley of the State of California, USA, can be considered a basin where all stresses that can possibly affect the management of groundwater resources seem to have come together: a vibrant economy that depends on water, a relatively dry climate, a disparity between water demand and availability both in time and space, heavily managed stream flows that are susceptible to water quality issues and sea level rise, degradation of aquifer conditions due to over-pumping, and degradation of the environment with multiple species becoming endangered. Over the past fifteen years, the California Department of Water Resources has developed and maintained the Integrated Water Flow Model (IWFM) to aid in groundwater management and planning under complex, and often competing, requirements. This paper will describe features of IWFM as a generic modeling tool, and showcase several of its innovative applications within California.


Introduction
Groundwater is a key component of water resources management and planning in developed basins.It is generally in constant interaction with the surface flow system and acts as a buffer in replenishing the depleted surface water resources as well as meeting the unmet portion of water demand after surface water resources are utilized.Groundwater can be the only source of water to meet these demands in semi-arid or arid climates.Therefore, any water resources planning studies carried out for a developed basin must include the management of both surface and subsurface water resources while addressing the dynamic interaction between them.
Water resources management and planning in a basin are generally further complicated by the changing climate, competing interests among farmers, cities' and environmental requirements, as well as a regulatory environment that tries to balance flood control with water supply availability while honoring water rights [1,2].Additionally, maintaining the surface and subsurface water quality at or above a level that is suitable for agriculture, human consumption, and the environment tends to be an integral water resources planning component in many basins [3].Groundwater management and planning are generally affected by all these factors and a proper management plan must address many of these issues.
Modelers have developed different tools to address one or more of these issues in groundwater resources management and planning studies.Pure groundwater models simulate the flow in the aquifer but all stresses on the groundwater system must be calculated outside the model [4][5][6].The dynamic interdependence between groundwater and the rest of the hydrologic cycle as well as the regulatory environment is not simulated dynamically and must be pre-defined.
On the other hand, integrated hydrologic models incorporate the simulation of groundwater and most or all of the flow processes from the rest of the hydrologic cycle [7][8][9][10][11].They simulate groundwater flow dynamics and other flow process at varying degrees of complexity using different approaches.For example, HydroGeoSphere [7] solves the three-dimensional Richards equation for variably saturated flow in the aquifer.The location of the water table is found based on the simulated pressure field and soil moisture content.Surface water flow is simulated using a diffusion wave approximation of the two-dimensional Saint Venant equation.Surface and subsurface flows are linked to each other either by using the assumption that the hydraulic head is continuous or through a conductance term between the two domains.The Penn State Integrated Hydrologic Model (PIHM) uses a semi-discrete finite volume method to convert flow systems represented by partial differential equations (combination of saturated and unsaturated subsurface flow, one-dimensional Saint Venant equations for the channel flow and two-dimensional Saint Venant equations for the overland flow) into ordinary differential equations and fully couple them to the flow processes represented by ordinary differential equations (interception, snowmelt, and evapotranspiration) [8].MODFLOW solves the three-dimensional saturated groundwater equation using the finite difference approach [9].It has many optional components and packages that simulate different flow terms of the hydrologic cycle.The interaction between the surface and subsurface flows is achieved through a Cauchy-type boundary condition.GSFLOW integrates the Precipitation Runoff Modeling System (PRMS) rainfall runoff model with the MODFLOW groundwater model [10].Evapotranspiration, surface runoff, and deep percolation are simulated using empirical functions.Watershed Environmental Hydrology Model (WEHY) uses stochastically-averaged, physically-based conservation equations to represent the surface and subsurface flows, and their interactions [11].
Groundwater management and planning requires understanding the sources and sinks acting on the groundwater system.In developed basins, these sources and sinks are shaped by agricultural and urban development and the regulatory environment dictating the use of water resources as much as by natural forces such as precipitation and evapotranspiration.For instance, the amount of pumping is affected by the climate, agricultural crop characteristics, farm water management parameters, and availability of stream flows which are partially dictated by the operation of the upstream reservoirs and regulations dictating surface water use and in-stream flows.Therefore, to perform a comprehensive groundwater resources analysis and planning study, it is generally necessary to use more than one type of model and simulate not just the groundwater flow dynamics but also the entire hydrologic system, reservoir operations, agricultural water demands, etc.Additionally, these models must effectively and dynamically exchange information (i.e., they must be linked to each other) since each model will depend on some information another model produces.
The State of California Department of Water Resources (CADWR) has been developing and maintaining the Integrated Water Flow Model (IWFM) for the last fifteen years to make possible the planning of the groundwater resources under complex requirements [12].IWFM is an integrated hydrologic model that is both descriptive (routes the water through the hydrologic cycle when the stresses such as pumping and stream flow diversions are defined) and prescriptive (calculates the stresses dynamically based on the calculated water demands).Although the core component of IWFM is the simulation of groundwater flow in complex multi-layered aquifers, it strives to dynamically simulate the sources and sinks to the groundwater.Therefore, it also simulates a large portion of the hydrologic cycle; the agricultural and urban water demands given soil, crop, farm water management, population and per capita water use characteristics; water supply in terms of pumping and stream diversions to meet these demands; and the effect of pumping and diversions on the surface and subsurface water resources.Additionally, IWFM is designed in a modular way such that these modules (e.g., groundwater module, root zone module) are not IWFM-centric; IWFM in its entirety or its individual modules can be linked to other types of models such as reservoir systems analysis models and agricultural economics models quickly at different scales.
In this article, simulation components and the software design concepts used in the development of IWFM will be described.Several case studies that include IWFM's application to California's Central Valley, California Central Valley Groundwater-Surface Water Simulation Model (C2VSim), will be discussed to highlight the strengths of IWFM in effective groundwater resources analysis, planning, and management under complex climate, regulatory, and economic stresses.

General Description of IWFM
IWFM consists of several simulation components: groundwater, streams, lakes, unsaturated zone, and land surface and root zone flow components.All components are optional except the groundwater component.This section briefly describes each of these simulation components and how they are linked to each other.A more detailed description can be found in [13].

Groundwater
Simulation of groundwater flow is at the heart of IWFM (Figure 1).Three-dimensional, transient groundwater flow is simulated in complex, multi-layered aquifers.The aquifer layers can be a combination of confined, unconfined, and leaky layers separated by aquitards or aquicludes.IWFM can handle aquifer layers that pinch out and disappear at sections of the groundwater basin.Through the simulation period confined layers can become unconfined, or unconfined layers confined, as the simulated groundwater heads fluctuate in response to the stresses.
Three-dimensional groundwater flow is represented by the depth-integrated conservation equation at each aquifer layer interacting with the layers above and below them through leakage terms: where h = h(x,y,t) = groundwater head in an aquifer layer, (L); K = K(x,y) = saturated hydraulic conductivity, (L/T); D = D(x,y,t) = saturated thickness of the aquifer layer, (L); S s = S s (x, y, t) = storativity which is equal to the storage coefficient for a confined aquifer and specific yield for an unconfined aquifer, (dimensionless); h u = h u (x, y, t) = groundwater head at the aquifer layer that is above the layer being considered, (L); L u = L u (x, y) = leakage coefficient between the layer considered and layer above, (1/T); I u = indicator depending on the aquifer layer, equals 0 for the top-most aquifer layer and 1 for other layers, (dimensionless); h d = h d (x, y, t) = groundwater head at the aquifer layer that is below the layer being considered, (L); L d = L d (x, y) = leakage coefficient between the layer considered and layer below, (1/T); I d = indicator depending on the aquifer layer, equals 0 for the bottom-most aquifer layer and 1 for other layers, (dimensionless); q o = q o (x, y, t) = non-head-dependent sources and sinks, (L/T); q(h) = head-dependent sources and sinks, (L/T); ∇ = del operator, (1/L); x = horizontal x-coordinate, (L); y = horizontal y-coordinate, (L); and t = time, (T).Equation ( 1) is written for each aquifer layer simulated.Effects of pumping, distributed and point recharge, tile drains, and subsidence are simulated through the o q and q(h) terms that appear in Equation (1) (more details on the simulation of these processes can be found in [13]).The resulting set of equations is discretized using the Galerkin finite element method [14] and the implicit finite difference scheme in space and time, respectively.The resulting set of non-linear algebraic equations is then solved iteratively using the Newton-Raphson method [15].
Use of the finite element method allows effective representation of irregular model boundaries, line features such as streams and fault lines, and internal boundaries such as water management districts and city boundaries.It also allows varying the mesh resolution throughout the model domain, using a fine grid where necessary and a coarser grid elsewhere, without having to resort to child meshes as is common in finite difference approaches [16].IWFM allows a mixture of linear triangular and bilinear quadrilateral cells in the computational finite element mesh.The use of quadrilateral cells allows for smaller numerical problems, particularly when the root zone and land surface component (discussed later in this article) of IWFM is activated to dynamically compute groundwater pumping and recharge as well as land surface and root zone flows.
The groundwater component of IWFM accepts all the standard boundary conditions: specified head, specified flow, and general head boundary conditions.In addition, IWFM offers a "small watershed boundary condition" (Figure 1) that can be used to simulate, in a simple way, the dynamic surface and subsurface inflows into the model boundary from the adjacent watersheds.

Streams
Complex, dendritic stream channel networks are described by connecting the relevant nodes of the finite element grid, and the stream flows can be simulated using either the "instantaneous flow" Equation ( 1) is written for each aquifer layer simulated.Effects of pumping, distributed and point recharge, tile drains, and subsidence are simulated through the q o and q(h) terms that appear in Equation (1) (more details on the simulation of these processes can be found in [13]).The resulting set of equations is discretized using the Galerkin finite element method [14] and the implicit finite difference scheme in space and time, respectively.The resulting set of non-linear algebraic equations is then solved iteratively using the Newton-Raphson method [15].
Use of the finite element method allows effective representation of irregular model boundaries, line features such as streams and fault lines, and internal boundaries such as water management districts and city boundaries.It also allows varying the mesh resolution throughout the model domain, using a fine grid where necessary and a coarser grid elsewhere, without having to resort to child meshes as is common in finite difference approaches [16].IWFM allows a mixture of linear triangular and bilinear quadrilateral cells in the computational finite element mesh.The use of quadrilateral cells allows for smaller numerical problems, particularly when the root zone and land surface component (discussed later in this article) of IWFM is activated to dynamically compute groundwater pumping and recharge as well as land surface and root zone flows.
The groundwater component of IWFM accepts all the standard boundary conditions: specified head, specified flow, and general head boundary conditions.In addition, IWFM offers a "small watershed boundary condition" (Figure 1) that can be used to simulate, in a simple way, the dynamic surface and subsurface inflows into the model boundary from the adjacent watersheds.

Streams
Complex, dendritic stream channel networks are described by connecting the relevant nodes of the finite element grid, and the stream flows can be simulated using either the "instantaneous flow" approach or the kinematic wave approach.In both approaches, the following conservation equation is solved: where Q s = stream flow, (L 3 /T); A = stream flow cross-sectional area, (L 2 ); L s = stream channel length, (L); Q sin = all inflows into the stream along its length, L s , (L 3 /T); Q sout = all outflows from the stream along its length, L s , (L 3 /T); Q sint = stream-aquifer flows along its length, L s , (L 3 /T); and x s = spatial coordinate along the length of the stream, (L).
The Q sin term in Equation ( 2) consists of all inflows into the stream that IWFM can simulate: inflows from upstream channels, rainfall runoff, agricultural and urban return flows, inflows from small watersheds that are adjacent to the model boundary, inflows from tile drains, inflows from river bypass systems, inflows from upstream lakes, and inflows pre-specified by the user as time series input data.The Q sout term includes outflow to the downstream channel, surface water diversions to meet any water demands, water taken out of the stream into bypass channels, and evapotranspiration due to riparian vegetation.All of these inflow and outflow terms are optional depending on the specific application of IWFM.
In the instantaneous flow approach, it is assumed that a flow pulse travels instantaneously within the simulated stream network.This approach is applicable when the simulation time step is known to be larger than the characteristic time scale of the stream flows within the modeled basin.For instance, if it is known that stream flow travels from the upstream end to the downstream end of the stream network within the model domain within several days and a monthly simulation time step is selected for modeling, then the instantaneous flow approach can be safely used.Using this approach renders the left-hand side of Equation (2) zero: In the kinematic wave approach, the slope of the energy line is represented by the stream channel slope [17].The implementation of kinematic wave approach in IWFM is discussed in [13] in detail.This approach can be used when the stream flow does not travel through the simulated stream network within a single simulation time step requiring representation of the change in stream storage.
Stream-aquifer interaction is calculated at each stream node and the corresponding groundwater node as a head-dependent boundary condition for both systems.To convert stream flows to stream stage in order to calculate stream-aquifer interaction, the relationship between stream flow, Q s , and the stream stage is defined through user-specified rating tables specified at each stream node for the instantaneous flow approach.For the kinematic wave approach, this relationship is defined by Manning's equation [17].
Equation ( 2) is discretized using the fully implicit finite difference scheme with backward differencing [17].The resulting set of non-linear algebraic equations is combined with those obtained from the discretization of the groundwater equation and solved simultaneously.

Lakes
Lakes are represented in IWFM using the following conservation equation: where S k = lake storage, (L 3 ); Q kin = all inflows into lake, (L 3 /T); Q kout = all outflows from the lake, (L 3 /T); and Q kint = lake-aquifer flows, (L 3 /T).
Water 2016, 8, 592 The Q kin term in Equation (4) includes precipitation, inflows from streams either through streams flowing directly into the lake or diversions from streams into the lake, rainfall runoff, agricultural and urban return flows, and inflows from upstream lakes.The Q kout term includes lake surface evaporation, and lake outflow if the lake elevation exceeds a pre-specified maximum lake elevation.
Similar to stream-aquifer interaction, lake-aquifer interaction is simulated as a head-dependent boundary condition.
Lakes are identified by specifying the corresponding finite element cells.Based on the finite element nodes associated with these cells and the ground surface elevation specified at these nodes, IWFM internally generates a storage versus lake-surface-elevation rating table with the maximum storage calculated at the user-specified maximum lake elevation, which can be a time series input.
Equation ( 4) is discretized in time using a fully implicit method for each lake and the resulting set of algebraic equations is solved, along with equations representing groundwater and stream flows simultaneously, using the Newton-Raphson method.
The many optional inflows into the lakes that can be represented by IWFM make the lake component quite versatile.It has been used in the past to simulate natural and man-made lakes, managed wetlands, ponding operations in rice fields, and spreading basins used in managed aquifer recharge operations.

Land Surface and Root Zone
The land surface and root zone component of IWFM is one of the most important components and makes IWFM a powerful tool in groundwater management studies.The land-use-based irrigation water demands, agricultural return flow generation, and soil and irrigation method dependent groundwater recharge that affect groundwater management in a basin are all simulated in this component.A brief description of this component is given here; more information can be found in [18,19].
For given precipitation and irrigation rates, the land surface and root zone component simulates rainfall runoff, return flow of irrigation water, infiltration of irrigation and precipitation, actual evapotranspiration, vertical movement of moisture within the root zone, and change in soil moisture storage within the root zone for a variety of land use types including agricultural crops, managed wetlands, urban areas, native vegetation, and riparian vegetation.Rainfall runoff and irrigation return flow become inflows to streams or lakes, depending on the model setup, while moisture leaving the root zone vertically at its lower boundary eventually becomes recharge to the groundwater.The land surface and root zone component is also capable of calculating irrigation water demand in case irrigation amounts are not known, such as in a planning study or when historical measurements are not available.
The following conservation equation is used to represent soil moisture dynamics within the land surface and root zone (Figure 2): where θ = soil moisture content, (L/L); Z = rooting depth, (L); P = precipitation, (L/T); R P = rainfall runoff, (L/T); A w = irrigation, (L/T); R f,ini = initial irrigation return flow, (L/T); U = re-used portion of the initial return flow, R f,ini , (L/T); R f = net irrigation return flow, (L/T); D r = outflow due to drainage of rice or managed wetland ponds, (L/T); D p = deep percolation, i.e., moisture outflow at the lower boundary of the root zone that is destined to become groundwater recharge, (L/T); and ET = actual evapotranspiration, (L/T).In Equation ( 5) it is assumed that the soil moisture movement within the root zone is predominantly one-dimensional in the vertical and its horizontal movement is negligible.Precipitation, P, is specified as the time series input for each finite element cell.A modified version of the SCS curve number method is used to convert the original event-based methodology to a time-continuous methodology for the calculation of rainfall runoff, R P [18,20,21].Irrigation, A w , is either specified as time series data for agricultural crops or calculated dynamically (explained later in the article) based on soil, climate, crop, and farm water management parameters.Initial irrigation return flow, R f,ini , and the re-used portion of the initial return flow, U, are specified as time series fractions of the irrigation, A w .Pond drainage, D r , is calculated based on user-specified time series ponding depths; it is greater than zero if the pond depths are decreasing in time and zero otherwise.D r is only calculated for rice fields and managed wetlands and is always zero for other land use types.
It is also worth noting that soil moisture storage, Z in Equation ( 5), for ponded land use types includes both moisture storage within the soil and the pond depth.
In calculating deep percolation, D p , the Darcy equation is used with the assumption that the vertical head gradient within the root zone is unity (i.e., the soil moisture is assumed to be uniformly distributed across the root zone): where K ur = unsaturated vertical hydraulic conductivity of the root zone, (L/T); h r = total head (sum of pressure and elevation heads) within the root zone, (L), and z = vertical coordinate, (L).Based on user choice, K ur is calculated by either using the van Genuchten-Mualem approach [22,23]: or Campbell's approach [24]: where K sr = the saturated vertical hydraulic conductivity of the root zone, (L/T);  T = the total porosity of the root zone, (L/L); and  = the pore size distribution index (dimensionless).Precipitation, P, is specified as the time series input for each finite element cell.A modified version of the SCS curve number method is used to convert the original event-based methodology to a time-continuous methodology for the calculation of rainfall runoff, R P [18,20,21].Irrigation, A w , is either specified as time series data for agricultural crops or calculated dynamically (explained later in the article) based on soil, climate, crop, and farm water management parameters.Initial irrigation return flow, R f,ini , and the re-used portion of the initial return flow, U, are specified as time series fractions of the irrigation, A w .Pond drainage, D r , is calculated based on user-specified time series ponding depths; it is greater than zero if the pond depths are decreasing in time and zero otherwise.D r is only calculated for rice fields and managed wetlands and is always zero for other land use types.It is also worth noting that soil moisture storage, θZ in Equation ( 5), for ponded land use types includes both moisture storage within the soil and the pond depth.
In calculating deep percolation, D p , the Darcy equation is used with the assumption that the vertical head gradient within the root zone is unity (i.e., the soil moisture is assumed to be uniformly distributed across the root zone): where K ur = unsaturated vertical hydraulic conductivity of the root zone, (L/T); h r = total head (sum of pressure and elevation heads) within the root zone, (L), and z = vertical coordinate, (L).Based on user choice, K ur is calculated by either using the van Genuchten-Mualem approach [22,23]: or Campbell's approach [24]: where K sr = the saturated vertical hydraulic conductivity of the root zone, (L/T); θ T = the total porosity of the root zone, (L/L); and λ = the pore size distribution index (dimensionless).Finally, actual evapotranspiration, ET, is calculated as a function of potential evapotranspiration, ET pot , specified as time series input data and the soil moisture content, θ.A detailed discussion on this approach can be found in [25].Based on this calculation, ET is equal to ET pot when θ is greater than half of the field capacity, θ f , and decreases linearly with respect to θ when θ is less than half of θ f but greater than the wilting point, θ wp .When θ is below θ wp , ET is assumed to be zero.The relationship between ET and θ is shown in Figure 3.

Water Demand and Automatic Supply Adjustment
The land surface and root zone component described previously requires the irrigation water, A w , as input data.However, these data are not always available.For instance, in California, groundwater pumping, which can make up a substantial part of A w , is not measured or reported historically.Similarly, in a water resources management planning study with future climate, land use distribution, and farm water management practices, A w will not be available and needs to be estimated.To estimate A w in these situations, IWFM conveniently calculates irrigation water demand, which can be estimated to be equal to A w as long as there is enough groundwater storage for pumping and enough stream flow to divert.For the calculation of irrigation water demand, IWFM uses an irrigation scheduling type of approach that is described in [18] and [25], and combines it with the land surface and root zone conservation Equation (5):  Equation ( 5) is discretized in time using the implicit scheme and solved for each land use type (agricultural crop, managed wetland, urban outdoors, native vegetation, or riparian vegetation) at each finite element grid cell.The land surface and root zone component is linked to the groundwater component through pumping as part of the irrigation water, A w , and through deep percolation, D p .It is also linked to the stream flows through stream diversions as part of A w , rainfall runoff, R P , and irrigation return flow, R f .Because of this interdependence between the land surface and root zone flow processes and the stream and groundwater flows, Equation ( 5) is solved and the deep percolation, D p , rainfall runoff, R P , and irrigation return flow, R f , are calculated anew at each Newton-Raphson iteration for the solution of the combined stream, lake, and groundwater flow equations.

Unsaturated Zone
In basins with groundwater heads much lower than the ground surface elevations, it is sometimes necessary to simulate the vertical movement of the deep percolation, D p , through an unsaturated zone and to represent the time delay for D p to recharge the groundwater.This is achieved in IWFM by activating the optional unsaturated zone component.The unsaturated zone can be divided into as many layers as necessary and the following conservation equation is used for each layer: where D u = thickness of the unsaturated layer, (L); θ u = the moisture content of the unsaturated zone layer that is assumed to be uniform through the entire thickness of the layer, (L/L); q uin = moisture inflow into the unsaturated layer from the layer above, (L/T); and q uout = moisture outflow from the unsaturated layer into the layer below, (L/T).
For the topmost unsaturated layer, q uin is equal to the deep percolation, D p , from the root zone whereas for the lowest layer, q uout represents the recharge to the groundwater, which is called the "net deep percolation", netD p , in IWFM.q uout is calculated using either Equation ( 8) or (9) depending on the user's choice.Equation ( 10) is discretized using the implicit scheme and solved for each finite grid cell at each Newton-Raphson iteration for the solution of the combined stream, lake, and groundwater equations.As the groundwater table rises and falls, the thicknesses of the unsaturated layers may change.For instance, as the groundwater table rises, the thickness of the lowest unsaturated layer decreases until it becomes zero.When this happens, IWFM no longer simulates the moisture content within the lowest layer.As the groundwater table continues to rise, the thickness of the second lowest unsaturated layer starts decreasing until it becomes zero and IWFM stops simulating the moisture content within this layer.Similarly, as the groundwater table falls, the thicknesses of the unsaturated layers increase one at a time and as the layer thicknesses become non-zero IWFM again starts simulating the moisture content in these layers.

Water Demand and Automatic Supply Adjustment
The land surface and root zone component described previously requires the irrigation water, A w , as input data.However, these data are not always available.For instance, in California, groundwater pumping, which can make up a substantial part of A w , is not measured or reported historically.Similarly, in a water resources management planning study with future climate, land use distribution, and farm water management practices, A w will not be available and needs to be estimated.To estimate A w in these situations, IWFM conveniently calculates irrigation water demand, which can be estimated to be equal to A w as long as there is enough groundwater storage for pumping and enough stream flow to divert.
For the calculation of irrigation water demand, IWFM uses an irrigation scheduling type of approach that is described in [18,25], and combines it with the land surface and root zone conservation Equation ( 5): where A * w = irrigation water demand, (L/T); θ trg = irrigation target soil moisture content, (L/L); θ min = irrigation trigger moisture content, (L/L); D ptrg = deep percolation at irrigation target moisture content, (L/T); ET trg = actual evapotranspiration at irrigation target moisture content, (L/T); f R f,ini = fraction of irrigation amount that becomes initial return flow, (dimensionless); f U = fraction of irrigation amount that is captured and re-used after becoming initial return flow, (dimensionless); and ∆t = simulation time step length, (T).
θ trg , θ min , f R f,ini , and f U are all time series input parameters, and they can be used to represent many different irrigation practices.For instance, regulated deficit irrigation practices can be simulated by decreasing θ trg ; frequent irrigation application can be simulated by decreasing θ trg , increasing θ min , or both; increased irrigation amounts to promote leaching can be simulated by increasing θ trg ; efficiency of different irrigation methods (e.g., farrow irrigation versus drip irrigation) can be simulated using different f R f,ini values.
Once A * w is calculated, IWFM allows automatic adjustment of pumping and stream diversions to meet A w to the extent that there is enough groundwater storage for pumping and stream flow to support required diversions.Maximum groundwater pumping and diversion rates can be enforced during the automatic adjustment process to represent pump or diversion canal conveyance capacities.Ultimately, the sum of the adjusted pumping and diversions represent the A w term that appears in Equation (5).In cases where there is not enough groundwater storage or pumping capacity to support the required pumping, and not enough stream flow or conveyance capacity to support required diversions, A w will be less than A * w and all flow terms represented in Equation ( 5) will be scaled down accordingly.

Software Design of IWFM
IWFM is written in the Fortran 2008 programming language and uses the latest software design procedures through object-oriented programming.Simulation components of IWFM are extensible without affecting their previous implementations.For instance, a new stream flow routing approach with different assumptions (e.g., dynamic routing) can be developed and added to the already available stream flow routing approaches.Each simulation component of IWFM is developed such that they can quickly be turned into stand-alone modeling software.One example of this is the land surface and root zone component, which is also available as a stand-alone modeling tool named the IWFM Demand Calculator (IDC) [18,19].It has been used as an integral part of IWFM [26,27], as a stand-alone tool in developing agricultural water management plans and integrated regional water management plans for several water districts in California [28,29] and Idaho [30], and as a pre-processing tool to develop groundwater pumping rates for groundwater modeling studies carried out with other models [31,32].
To address groundwater and surface water management issues, complex models such as reservoir system analysis models or agricultural economics models are generally required along with an integrated hydrologic model such as IWFM.To facilitate easy linkage to such models, IWFM code is written in such a way that it can be compiled into an Application Programming Interface (API) that can be linked to other modeling software at the code level.This approach eliminates the need to redevelop complex algorithms within IWFM that are already available in other modelling software.Similarly, an API makes IWFM algorithms available to other modeling software without the need to redevelop these algorithms within the other software.

California Central Valley Groundwater-Surface Water Simulation Model (C2VSim)
California's Central Valley (referred to as the Valley hereafter) is a 51,000 km 2 basin that is heavily developed and has a highly engineered water distribution network (Figure 4a).It is one of the most agriculturally productive regions in the world.In 2009, the total value of agricultural production of the Valley was more than $21 billion, or 62% of the state total [33].The total area of the Valley's farms amounts to 1% of U.S. farmland, but it produces 25% of the food the U.S. population consumes, including 40% of fruits, nuts, and other table foods [34].Additionally, nearly 8 million people live in the Valley [35].
Much of the Valley's agricultural production occurs in a semi-arid to arid environment, and relies on large quantities of irrigation water.Water demand and supply are mismatched both in space and time: a large portion of the water demand occurs in the southern part of the Valley during the summer, but the northern part of the Valley receives most of the precipitation during the winter as rain on the Valley floor and snow on the surrounding mountains.Two large surface water storage and conveyance systems, the Central Valley Project (CVP) and the State Water Project (SWP), store, regulate, and move water from the northern part of the Valley to the southern part.The Sacramento and San Joaquin Rivers Delta (the Delta) is in the hub of this water resources management operation.The rainfall and snowmelt stored in the northern reservoirs in the winter and spring and released throughout the summer flow to the Delta, and are then transferred to the south through CVP and SWP canals.The combined CVP and SWP systems include 54 reservoirs with a total capacity of about 20 km 3 , 16 power plants, and more than 1900 km of canals, tunnels, and pipelines.On average the system provides around 11.5 km 3 of water annually, of which 62% is used for agricultural irrigation, 25% for municipal and industrial purposes, and the rest is dedicated to fish, wildlife, and their habitat [36,37].Besides providing water supply where and when it is needed, CVP and SWP systems are also used for flood protection, power generation, recreation, fish and wildlife enhancement, and water quality protection in the Delta.
Water 2016, 8, 12 12 of 22 (a) The C2VSim model described above has a coarse mesh and its intended use is to understand the implications of proposed water resources management projects, climate change, etc. on the Valley overall rather than their impacts at a local scale.A coarse mesh also allows a fast computer run time, which is particularly valuable when C2VSim is linked to other types of models such as reservoir system analysis models and agricultural economics models.Several examples of such studies are described in the next section.Another version of C2VSim with a finer mesh (C2VSim-FG) that includes 32,537 cells with an average area of 1.6 km 2 and 4529 stream nodes is also being developed to analyze the local-scale water resources projects.
Historically, groundwater pumping is neither measured nor reliably reported in California.To estimate the historical groundwater pumping rates, C2VSim utilizes the automatic supply Despite this heavily engineered surface water management system, groundwater is still a major source of water in the Valley.In an average year, it is estimated that groundwater meets approximately 40% of the water demand in the Valley, and as much as 60% in a dry year [38].In some parts of the Valley, it is the only source of water.In 2007, two-thirds of the municipal and industrial water was supplied by groundwater.Groundwater use in the Valley is approximately 13% of total U.S. groundwater use, making the Valley aquifers the second largest source of groundwater in the USA after the High Plains aquifer [39,40].
The California Central Valley Groundwater-Surface Water Simulation Model (C2VSim) is an application of IWFM to simulate the groundwater and surface water flow dynamics within the Valley, the interaction between these flow systems, agricultural and urban water demands, and the effect of water resources development on the surface and subsurface flows [26,41].The extent of the Valley is represented with a finite element grid that includes 1392 cells with an average cell size of 37 km 2 (Figure 4b).The top of the aquifer is defined as the land surface and the bottom is defined as the impermeable basement rocks or, where present, the base of the fresh water.It is represented in three vertical layers with varying thicknesses across the model domain, including an additional confining layer in the southwestern part of the Valley.The number of aquifer layers and their thickness were defined based on the available screening depths from about 40,000 well completion reports in the absence of enough hydro-geologic studies to parameterize the entire aquifer.Three layers were considered sufficient to represent the vertical distribution of groundwater pumping while maintaining fast computer run times.On average, 40% of pumping occurs in the top layer and 60% in the second layer.The third layer represents the deep aquifer with no pumping that generally provides upward flow into the second layer.The stream network is represented by 449 stream nodes grouped into 75 reaches and stream diversions are simulated at 246 locations.The historical C2VSim model contains measured boundary stream inflows in terms of reservoir releases, precipitation, potential evapotranspiration, surface water diversions, and land use and crop areas from October 1921 through September 2009.Subsurface boundary inflows into the model are represented by using the small-stream watershed boundary condition of IWFM for 210 ungauged small watersheds that are adjacent to the model boundary (Figure 4b).The small-stream watershed boundary condition of IWFM is used to dynamically compute groundwater boundary inflows as a function of precipitation, land surface, and aquifer storage parameters within small watersheds and is described in more detail in [13].
The C2VSim model described above has a coarse mesh and its intended use is to understand the implications of proposed water resources management projects, climate change, etc. on the Valley overall rather than their impacts at a local scale.A coarse mesh also allows a fast computer run time, which is particularly valuable when C2VSim is linked to other types of models such as reservoir system analysis models and agricultural economics models.Several examples of such studies are described in the next section.Another version of C2VSim with a finer mesh (C2VSim-FG) that includes 32,537 cells with an average area of 1.6 km 2 and 4529 stream nodes is also being developed to analyze the local-scale water resources projects.
Historically, groundwater pumping is neither measured nor reliably reported in California.To estimate the historical groundwater pumping rates, C2VSim utilizes the automatic supply adjustment feature of IWFM.Agricultural water demand is calculated by the land surface and root zone component as a function of climate, soil, crop and farm water management parameters, while urban demand is calculated as a function of population and per capita water use.Then, C2VSim estimates the historical pumping as the calculated water demand less the historically measured stream diversions to meet these demands.

Results and Discussion
There is a growing body of studies that use IWFM for groundwater resources analysis, management, and planning under complex requirements [27,[42][43][44][45][46][47][48][49].In this section, some of these studies carried out by C2VSim that highlight the features IWFM that are instrumental in effective groundwater resources analysis and management are described.

Effect of Long-Term Droughts on Groundwater in the Valley under Constant Land Use
Miller et al. [42] used C2VSim to study the effects of long-term droughts with varying intensities on the groundwater resources of the Valley.Their study stemmed from the observations of signs of paleoclimate droughts lasting up to 200 years imprinted on sediments and tree rings, as well as future climate change predictions for reduced snowpack and increased winter runoff.The goals of their study were to quantify the impacts of long-term droughts on water storage, the ability of groundwater to offset the adverse effects of the reduced stream flows, and the time it took for the system to recover, if at all, after the drought ended.
10-, 20-, 30-and 60-year droughts were represented by 30% (slight drought), 50% (moderate drought), and 70% (severe drought) reductions in mean annual precipitation, with 10-year spin-up and 30-year recovery periods before and after the droughts, respectively.The droughts were constructed by randomly selecting years from the historical record going back to 1922 that satisfied the desired reduction in the mean annual precipitation and appending these years together.The years with reduced precipitation also included inflows into the reservoirs that provide boundary surface inflows into the Valley.These reservoirs are operated based on predefined operating rules and the downstream diversions are dictated by complex water rights.In order to estimate the reservoir releases and the downstream diversions under the hypothetical droughts, CalSim, a reservoir systems analysis model for the CVP and SWP systems [50], was run using the hypothetical precipitation and reservoir inflow patterns.The resulting reservoir releases and diversion amounts were used as input data to C2VSim.To isolate the response of the system to the reduced flows, the land use and crop distribution were kept constant at 2003 levels.
With predefined precipitation, land use and crop distribution, stream boundary inflows, and stream diversions, C2VSim was able to calculate crop water demands under the hypothetical drought scenarios and, by using its automatic supply adjustment feature, estimate the groundwater pumping as the difference between the demand and the stream diversions.With the estimated groundwater pumping, it was then able to route the water through the entire surface-subsurface flow system and predict the impact of the drought on the groundwater storage, streams, and stream-aquifer interactions.
The drought scenarios imposed, on average, 39%, 22%, and 13% reductions in stream diversions across the Valley for the severe, moderate, and slight droughts, respectively.In return, groundwater pumping was increased to meet the water demand.Miller et al. [42] reported that Valley-wide average groundwater elevations fell by 46 m, 33 m, and 17 m for severe, moderate, and slight 60-year droughts, respectively, in response to increased groundwater pumping and reduced recharge due to precipitation and stream flows.Figure 5 shows the response of the groundwater system at different parts of the Valley in terms of the average change in groundwater heads for the slight and severe 60-year droughts, including the 10-year spin-up and 30-year recovery periods.The figure shows that the very southern part of the Valley, Tulare Basin, is impacted the most: groundwater levels fall, on average, 34 m and 80 m, for slight and severe droughts, respectively.Additionally, although the rest of the Valley recovers somewhat after 30 years with average precipitation and stream inflows, Tulare Basin does not.Figure 5 also shows that a new equilibrium with lower groundwater elevations compared to pre-drought levels may be reached across the Valley.
Miller et al. [42] also analyzed the impact of long-term droughts in groundwater recharge, stream-aquifer interaction, and aquifer storage.In a traditional groundwater modeling study, the increased pumping in response to decreased precipitation and stream diversions is estimated outside the groundwater model as a preprocessing step.However, in such an approach the complex and non-linear feedback mechanisms between the land surface and the groundwater cannot be predicted or effectively simulated.On the other hand, these feedback mechanisms are properly simulated in IWFM through the full linkage of the demand calculations to the climatic forcings (i.e., precipitation and evapotranspiration) and to the water supply management to meet these demands.In summary, the dynamic calculation of agricultural water demand as a function of climatic forcings, automatic supply adjustment feature to dynamically adjust pumping to meet any unmet demands after stream diversions, and the dynamic calculation of groundwater recharge in response to reduced precipitation and diversions as well as increased pumping allowed Miller et al. [42] to effectively carry out their analysis.and evapotranspiration) and to the water supply management to meet these demands.In summary, the dynamic calculation of agricultural water demand as a function of climatic forcings, automatic supply adjustment feature to dynamically adjust pumping to meet any unmet demands after stream diversions, and the dynamic calculation of groundwater recharge in response to reduced precipitation and diversions as well as increased pumping allowed Miller et al. [42] to effectively carry out their analysis.

Effect of Long-Term Droughts on Groundwater in the Valley under Varying Land Use
Miller et al. [42] assumed a constant land use and crop distribution over the long-term droughts considered to isolate the absolute impact of the drought on the groundwater resources in the Valley.However, in a real drought, farmers try to adapt to the reduced flows by taking the economic constraints into consideration and changing the crops they plant.Dale et al. [43] took the study performed by Miller et al. [42] one step further to estimate farmers' crop choices in the face of prolonged droughts and the impact, if any, of the varying crop distribution on the groundwater resources of the Valley.
To represent the farmers' crop choices, they emulated the Central Valley Production Model, (CVPM) in C2VSim.CVPM is an economic model for irrigated agricultural production that uses an optimization technique known as Positive Mathematical Programming to estimate the changes in crop area as a function of surface water availability, cost of water, and cost of energy, along with legal, physical, and economic constraints [51,52].
To emulate CVPM in C2VSim, Dale et al. [43] developed crop adaptation equations based on a multinomial logit regression analysis of crop areas generated by 1000 CVPM runs using a Monte Carlo analysis with randomly generated surface water supply availability and depth to groundwater values.The resulting logit function representation of CVPM predicted the share of a given crop in a given C2VSim subregion (Figure 4b) as a function of the percent reduction in surface water supply with respect to base conditions and the average depth to groundwater in that subregion.These logit functions were coded into C2VSim to dynamically calculate the crop areas at the beginning of each cropping season.
In this study, depth to groundwater was used as a surrogate for the cost of pumping groundwater.For average depth to groundwater to properly represent the pumping cost, it was calculated as a pumping-weighted average for each subregion.In this approach, groundwater heads at each node at each aquifer layer were weighted with respect to effective pumping on that node when calculating the average depth to groundwater.This way, groundwater heads in undeveloped

Effect of Long-Term Droughts on Groundwater in the Valley under Varying Land Use
Miller et al. [42] assumed a constant land use and crop distribution over the long-term droughts considered to isolate the absolute impact of the drought on the groundwater resources in the Valley.However, in a real drought, farmers try to adapt to the reduced flows by taking the economic constraints into consideration and changing the crops they plant.Dale et al. [43] took the study performed by Miller et al. [42] one step further to estimate farmers' crop choices in the face of prolonged droughts and the impact, if any, of the varying crop distribution on the groundwater resources of the Valley.
To represent the farmers' crop choices, they emulated the Central Valley Production Model, (CVPM) in C2VSim.CVPM is an economic model for irrigated agricultural production that uses an optimization technique known as Positive Mathematical Programming to estimate the changes in crop area as a function of surface water availability, cost of water, and cost of energy, along with legal, physical, and economic constraints [51,52].
To emulate CVPM in C2VSim, Dale et al. [43] developed crop adaptation equations based on a multinomial logit regression analysis of crop areas generated by 1000 CVPM runs using a Monte Carlo analysis with randomly generated surface water supply availability and depth to groundwater values.The resulting logit function representation of CVPM predicted the share of a given crop in a given C2VSim subregion (Figure 4b) as a function of the percent reduction in surface water supply with respect to base conditions and the average depth to groundwater in that subregion.These logit functions were coded into C2VSim to dynamically calculate the crop areas at the beginning of each cropping season.
In this study, depth to groundwater was used as a surrogate for the cost of pumping groundwater.For average depth to groundwater to properly represent the pumping cost, it was calculated as a pumping-weighted average for each subregion.In this approach, groundwater heads at each node at each aquifer layer were weighted with respect to effective pumping on that node when calculating the average depth to groundwater.This way, groundwater heads in undeveloped areas with only native vegetation did not affect the calculated averages and the groundwater heads that the farmers were most exposed to during pumping were properly represented.
At the beginning of each cropping season (generally in October), C2VSim calculated weighted-average depth to groundwater values for each subregion.Using these values along with the pre-specified percent reductions in surface water availability for each subregion, the crop adaptation functions emulating CVPM calculated crop areas for the next year in each subregion.C2VSim then used these dynamically calculated crop areas to calculate irrigation water demand and the resulting pumping requirement.The resulting tool was a hydro-economic model with a distributed integrated Water 2016, 8, 592 16 of 22 hydrologic model that was capable of simulating the complex interdependency between climate, water supply availability, agricultural economics, the impact of farmers' crop choices on the water resources, and the co-evolution of hydrology and agricultural development under any climate conditions.
The C2VSim-CVPM hydro-economic model was then used to assess the impacts of the hypothetical drought conditions generated by Miller et al. [42] on the Valley's groundwater resources when farmers adapted to the drought.Figure 6 shows the total irrigated area in the Valley for the 60-year severe drought.During the 10-year spin-up period, the Valley supports about 2.4 million hectares of irrigated land.After 60 years of severe drought with a 39% decrease in the surface water availability, the irrigated land area drops to about 2.1 million hectares, only about a 10% reduction.After the drought, irrigated land increases only to about 2.2 million hectares.

conditions.
The C2VSim-CVPM hydro-economic model was then used to assess the impacts of the hypothetical drought conditions generated by Miller et al. [42] on the Valley's groundwater resources when farmers adapted to the drought.Figure 6 shows the total irrigated area in the Valley for the 60year severe drought.During the 10-year spin-up period, the Valley supports about 2.4 million hectares of irrigated land.After 60 years of severe drought with a 39% decrease in the surface water availability, the irrigated land area drops to about 2.1 million hectares, only about a 10% reduction.After the drought, irrigated land increases only to about 2.2 million hectares.
Dale et al. [43] explained the resilience of the agricultural production in the face of severe longterm droughts with the abundance of groundwater storage in the Valley and the natural recharge into the groundwater from streams that limits large drops in the groundwater elevations.They reported an 84% increase in pumping, from 7.4 km 3 /year before drought to 13.6 km 3 /year during the drought.On average, simulated groundwater elevations in the Valley fell about 27 m, while groundwater seepage into the streams dropped from about 3.1 km 3 /year pre-drought to about 2 km 3 /year during the drought.They also reported that the amount of groundwater storage lost to subsidence was about 2.3 km 3 at the end of the 60-year drought.Dale et al. [43] explained the resilience of the agricultural production in the face of severe long-term droughts with the abundance of groundwater storage in the Valley and the natural recharge into the groundwater from streams that limits large drops in the groundwater elevations.They reported an 84% increase in pumping, from 7.4 km 3 /year before drought to 13.6 km 3 /year during the drought.On average, simulated groundwater elevations in the Valley fell about 27 m, while groundwater seepage into the streams dropped from about 3.1 km 3 /year pre-drought to about 2 km 3 /year during the drought.They also reported that the amount of groundwater storage lost to subsidence was about 2.3 km 3 at the end of the 60-year drought.
A comparison of the studies performed by Dale et al. [43] and Miller et al. [42] is revealing: farmers' crop choices based on economic constraints may offset some of the impacts of the drought on the groundwater system (27 m average decrease in the Valley-wide groundwater elevations in Dale et al. [43] with varying crop areas versus 46 m in Miller et al. [42] with fixed crop areas for the 60-year severe drought).
The Dale et al. [43] study also showed that, although an increase in groundwater pumping allowed farmers to maintain the area of irrigated lands, the cost to the aquifer due to subsidence, and to the environment and wildlife habitat due to decreased stream flows, can be undesired consequences.Additionally, they found that as the depth to groundwater increased (raising water costs), the area of low-value crops decreased and the area of both fallow and high-value crops such as tree crops and vines increased.The increase in the area of tree crops and vines causes "hardening" of the water demand and can have consequences for water management in future droughts since, unlike for annual crops, the farmers cannot fallow their lands and wait out the drought.
IWFM is both a prescriptive (groundwater pumping and stream diversions can be calculated dynamically to meet the calculated water demands) and a descriptive (once the pumping and diversions are calculated, the movement of water through the combined land surface, root zone, and groundwater systems can be simulated) tool.Dale et al. [43] took advantage of this feature of IWFM as well as its object-oriented software design to quickly modify the source code to incorporate the crop adaptation equations in their study.

Impact of 2012-2016 Drought on California's Agricultural Economy
Dale et al.'s study [43] required Monte Carlo simulations and many runs of CVPM to develop the crop adaptation equations that were used to emulate CVPM within C2VSim.Although this is a plausible approach, it has two drawbacks: i) each time CVPM is updated or modified, the Monte Carlo simulations have to be carried out again to update the crop adaptation equations, and ii) crop adaptation equations only represent a simplification of CVPM.
To avoid generating crop adaptation equations, Howitt et al. [48] linked C2VSim and the Statewide Agricultural Production (SWAP) model [53], an enhanced version of CVPM, to quantify the economic impact of the 2012-2016 droughts on California's agricultural industry during 2014.Their study also tried to estimate the impact if the drought lasted another two years through 2016.
In the linked C2VSim-SWAP hydro-economic model, similar to the study performed by Dale et al. [43], C2VSim calculated a pumping-weighted-average depth to groundwater at each subregion at the beginning of each cropping season.This information was then passed to SWAP, which, along with the percent reduction in surface water availability and cost of energy, calculated the optimum crop areas in each subregion.The crop areas were then used by C2VSim to calculate water demands, pumping and diversions to meet these demands, and all the hydrologic flow processes within the Valley for that cropping season.This process was repeated for the simulation period of three years.The groundwater simulations produced dynamic depth to groundwater values, which directly translated to the increased cost of replacing reduced surface water diversions with groundwater.These increased water costs played an important role in farmers' decisions on crop choices or fallowing.
Howitt et al. [48] concluded that the 2014 drought would increase the groundwater's share in meeting the agricultural water demand to 53%, compared to about 30% in an average year.This pointed out the importance of groundwater in offsetting the impact of drought on agriculture and of developing groundwater sustainability measures to maintain the long term health of the agricultural economy.They predicted that increased pumping cost farmers an additional $447 million while about 166,000 hectares were fallowed, leading to $800 million in lost farm revenues.Most of the fallowing occurred for low-value crops such as irrigated pasture, corn, and dry beans as the scarce water was directed to high-value crops such as vegetables, non-tree fruits, and permanent crops.Similar to the findings of Miller et al. [42] and Dale et al. [43], they concluded that Tulare Lake Basin was impacted the most due to its heavy reliance on groundwater.Assuming that the drought would continue for another two years through 2016, they predicted that the groundwater levels would continue to fall, escalating pumping costs and land subsidence, drying up of existing wells, and increasing fallowing, potentially costing Valley farmers an additional $1 billion a year.
The object-oriented software design of IWFM was instrumental in linking C2VSim to SWAP directly for this study.As with the previous studies described in this paper, the ability of IWFM to dynamically calculate water demands and allocate water supplies through pumping and diversions played a crucial role in the success of the study.

Linking a Groundwater Simulation Model to a Reservoir System Analysis Model
In developed basins where upstream reservoirs are operated to meet multiple objectives such as flood control, water supply storage, wildlife habitat protection, etc. based on pre-defined rules, it is almost impossible to perform groundwater resources management and planning separate from the operations of these reservoirs.In such basins, it is generally necessary to link groundwater simulation models with reservoir system analysis models.
Up until now, the approach to linking these two types of models has either been to over-simplify the aquifer system (lumped parameter bucket-type aquifer with a single aquifer layer) and the groundwater flow dynamics (linearization of the groundwater flow equations) [54][55][56][57], or to develop surrogate models such as linear response functions, artificial neural networks, etc. based on calibrated groundwater simulation models [58][59][60][61].The main goal of these approaches is to represent the groundwater flow dynamics as a manageable number of constraints in the optimization problem solved by the reservoir system analysis model.However, these approaches suffer from either oversimplification of the aquifer system so that the simulated groundwater flows do not represent reality, or the development of the surrogate models based on calibrated groundwater simulation models introduce extra steps and assumptions into the analysis of the linked groundwater and reservoir systems.
To avoid these drawbacks, Dogrul et al. [45] directly linked the groundwater component of C2VSim to CalSim, the reservoir systems analysis model for the CVP and SWP operations in California (Figure 4a).In this approach, the groundwater system was represented through 1392 finite-element cells and three aquifer layers without additional simplifications or assumptions.Subsidence and the effect of pumping from three separate aquifer layers on vertical flow generation were also represented.Stream-aquifer interaction was simulated through the iterative solution of the groundwater flow equations and the optimization problem for the reservoir system analysis model.In the linked model, the groundwater component of C2VSim estimated groundwater heads and the stream-aquifer interaction while CalSim estimated water-rights-based diversions and groundwater pumping to meet water demands as well as reservoir releases to provide the required diversions and instream flows for environmental and water quality requirements.Additionally, C2VSim's land surface and root zone component was used as a stand-alone model to simulate water demands and land surface and root zone flow processes that were used as input time series data to the linked model.Dogrul et al. [45] tested the linked model by comparing the results of historical simulations to a hypothetical scenario where historical precipitation was reduced by 25%.They compared the results from many angles: land surface flows, stream flows, groundwater flows and storage, subsidence, diversions and pumping, stream-aquifer interaction, horizontal and vertical groundwater head distribution, reservoir storage and releases, and the amount of surface water exports to south of Delta (Figure 4) through CVP and SWP canals to meet the water demands in the southern half of the Valley and Southern California.Their approach allowed a wholesome analysis of the codependent groundwater-surface water systems with upstream reservoir operations providing flood control, water supply, water rights compliance, and instream flow controls for environmental and legal purposes under a hypothetical climate change scenario.
In this study, IWFM's software design that allows easy conversion of its individual hydrologic simulation components into stand-alone modeling tools was instrumental in the linked models' capabilities.The stand-alone land surface and root zone component (named IWFM Demand Calculator, or IDC) was used as a pre-processor to the linked model.The groundwater component was converted into a dynamic link library (DLL) that acted as a library of groundwater simulation functions to the reservoir system analysis model, CalSim.CalSim was able to interact with the groundwater DLL and dynamically invoke appropriate groundwater simulation and data exchange functions when necessary.This allowed a seamless interaction and data exchange protocol between the reservoir system analysis and groundwater simulation models.

Conclusions
Groundwater is an integral component of the hydrologic cycle and water resources management in developed basins.There is interdependence between groundwater resources management, climate, surface water resources availability and management, economics, as well as legal and social constraints on water management.Therefore, in a developed basin, groundwater resources management is generally complex and requires attention to the co-evolution of groundwater resources with climate, surface water resource management, and economics.Although there are many groundwater simulation tools available, some in the form of integrated hydrologic models, most are not designed to handle all aspects of water resources management in a basin.
A major goal of CADWR's mission is planning for the water resources management in California to satisfy often competing needs of agricultural, urban and the environmental needs under changing climate, legal and social settings.IWFM was developed as a simulation tool to help scientists, modelers and managers effectively analyze groundwater resources in a basin under complex requirements.It is designed both as a prescriptive and descriptive tool.IWFM is prescriptive because it can calculate stresses on the groundwater and streams dynamically in terms of pumping and diversions, respectively, to meet dynamically computed water demands.It is at the same time descriptive because once these stresses are calculated, IWFM can route the water through the hydrologic system to describe where and how fast the water is moving.This feature of IWFM allows the simulation of the spatio-temporal evolution of the groundwater in response to changing climate, land use distribution, and surface water management.
Recognizing that a wholesome analysis and planning of groundwater management in a developed basin requires other simulation tools such as reservoir system analysis and agro-economic models, the IWFM code is developed in such a way that its hydrologic simulation components can quickly be turned into stand-alone simulation tools and both IWFM as a whole or these individual simulation components can be quickly and seamlessly linked to different types of models.This approach effectively taps into the power of already developed models addressing a different piece of the water management question rather than trying to re-create these tools within the IWFM code, which has generally been the standard approach.
This paper briefly describes the simulation methods used in IWFM and provides several applications highlighting its effectiveness in groundwater resources analysis, management, and planning through linkages to other types of models such as agro-economic models and reservoir system analysis models.Although the development of IWFM has been driven by the needs of California, it is a generic tool that can be used in any basin, as the issues California water managers face are not unique.IWFM source code and documentation are freely available to interested groups and can be downloaded from http://baydeltaoffice.water.ca.gov/modeling/hydrology/IWFM/.

Figure 2 .
Figure 2. Schematic representation of land surface and root zone flow processes simulated by IWFM.

Figure 2 .
Figure 2. Schematic representation of land surface and root zone flow processes simulated by IWFM.
w A = irrigation water demand, (L/T);  trg = irrigation target soil moisture content, (L/L);  min = irrigation trigger moisture content, (L/L); D ptrg = deep percolation at irrigation target moisture content, (L/T); ET trg = actual evapotranspiration at irrigation target moisture content, (L/T); R f ,ini f = fraction of irrigation amount that becomes initial return flow, (dimensionless); f U = fraction of irrigation amount that is captured and re-used after becoming initial return flow, (dimensionless); and t = simulation time step length, (T).

Figure 4 .
Figure 4. California's Central Valley (a) major reservoirs and conveyance canals of the State Water Project (SWP) and Central Valley Project (CVP) and (b) California Central Valley Groundwater-Surface Water Simulation Model (C2VSim).

Figure 5 .
Figure 5. Change in average groundwater heads across the Valley for (a) slight 60-year drought and (b) severe 60-year drought (source: Miller et al. [48])

Figure 6 .Figure 6 .
Figure 6.Total irrigated area before, during, and after a 60-year severe drought in the Valley (source: Dale et al. [49].