3.2.1  Introduction

Table 3.1 summarizes characteristics of optimization problem types and whether computed optimality is global or local. As mentioned elsewhere, a globally optimal solution is truly the best solution of the entire feasible solution space. A locally optimal solution is the best solution for part of the feasible solution space.

Below, we assume X_n is a real variable being optimized, I_n is an integer variable being optimized, and c_n , $c_n^{IP}$

3.2.2  Linear Optimization Problem

A linear programming (LP) optimization problem has only linear equations—variables exist only to the first power and there are no products of two or more variables. Equations 3.1 through 3.3 show a sample LP problem representation:

• Variables X ₁ and X ₂ are real numbers
• Linear objective function
  3.1 $$\text{maximize }{a}_1{X}_1+a_2{X}_2$$

  
• Linear constraints
  3.2 $$\begin{array}{c}{c}_{1,1}{X}_1+c_{1,2}{X}_2\le 5\\ c_{2,1}{X}_1+c_{2,2}{X}_2=3\end{array}$$

  
• Bounds
  3.3 $$X_1\ge b_1,X_2\le b_2$$

  

A sample linear objective function is the total sum of pumping. If pumping is positive in sign, this objective function is commonly maximized for water supply problems, and minimized for environmental or remediation problems.

3.2.3  Quadratic Optimization Problem

A quadratic programming (QP) optimization problem has all linear equations except for a quadratic objective function. The objective function includes one or more terms having a variable raised to the second power, or a product of two variables each to the first power. Equation 3.4 shows a sample QP problem representation:

• Variables X ₁ and X ₂ are real numbers
• Quadratic objective function
  3.4 $$\text{maximize }{a}_1{X}_1+a_2{X}_1{X}_2+a_3{X}_2$$

  
• Linear constraints (as in linear model)
• Bounds (as in linear model)

An example quadratic objective function computes the total cost of pumping (if cost includes pumping volume or rate times dynamic lift that is affected by pumping). Normally optimization minimizes this value.

3.2.4  Integer Optimization Problem

Integer programming (IP) problems involve integer variables. Sample integer variables are 0, 1, 2, 3, and so forth. Problems in which integer variables must be either 0 or 1 are termed binary integer or 0-1 programming problems. Binary (0 and 1) variables are often used to indicate that an action will not be performed, or will be performed, respectively.

Equations 3.5 through 3.7 show a sample IP problem:

• Variables I ₁ and I ₂ are nonnegative integers
• Objective function (Equation 3.5)
  3.5 $$\text{Maximize }{a}_1{I}_1+a_2{I}_2$$

  
• Constraints
  3.6 $$c_1{I}_1-c_2{I}_2<10$$

  
• Bounds
  3.7 $$I_1>0\left(\text{automatic}\right),I_2<5$$

  

A sample binary objective function represents the sum of the products of respective well installation cost times 0 (for a well that will not be installed) or 1 (for a well that will be installed). A sample nonbinary objective function includes the numbers of wells of each different type times the respective cost of each type.

3.2.5  Mixed Integer Optimization Problem

In mixed integer programming (MIP) problems, some variables are integers and some are real numbers.

Equations 3.8 and 3.9 and attendant bounds describe a sample MIP formulation:

• Variables include both real and integer numbers
  • X ₁, X ₂ are real numbers
  • I ₁, I ₂ are nonnegative integers
• Objective function
  3.8 $$\text{Minimize }{a}_1{X}_1+a_2{X}_2+a_3{I}_1+a_4{I}_2$$

  
• Constraints can be real, integer, and mixed
• As in LP for real variables
• As in IP for integer variables
• Mixture of real and integer variables as shown in Equation 3.9
  3.9 $$\begin{array}{l}{c}_{3,1}{X}_1+ c_{3,3}{I}_1>12\hfill \\ c_{4,2}{X}_2+c_{4,4}{I}_2> 8\hfill \end{array}$$

  
• Bounds
  • As in LP for real variables
  • As in IP for integer variables

A sample MIP objective function represents the total sum of the costs of groundwater pumping plus the costs of installing pumped wells. Optimization usually minimizes this sum.

3.2.6  Nonlinear Optimization Problem

Equations 3.10 and 3.11 and related bounds describe a nonlinear programming (NLP) optimization problem. An NLP problem has at least one nonlinear constraint, or an objective function of a power different than one or quadratic.

• Objective function
  • If only linear or quadratic, there must be a nonlinear constraint
  • Nonlinear
    3.10 $$\text{Maximize }{a}_1{X}_1^{0.33}+a_2{X}_1{X}_2^4$$

    
• Constraints
  • Linear (optional)
  • Nonlinear (OF or at least one constraint must be nonlinear)
    3.11 $$c_1{X}_1{X}_2^2-c_2{X}_2^3=10$$

    
• Bounds (as in linear model)

A sample NLP optimization problem is to minimize the extraction pumping needed to achieve target nonlinear groundwater contaminant concentrations by a certain time. Optimizable concentrations can also exist within an objective function. They are more nonlinear than quadratic.

3.2.7  Mixed Integer Nonlinear Optimization Problem

A mixed integer nonlinear programming (MINLP) combines MIP and NLP features. Commonly it has a MIP objective function and nonlinear constraints. A sample MINLP problem is to minimize cost of installing and operating wells to achieve groundwater contamination remediation goals by extraction pumping or pump-and-treat. Groundwater contaminant concentrations are nonlinear functions of pumping extraction.

3.3.1  Introduction

An S-O model includes representations of a physical system and an optimization problem. The physical system, optimization problem, and their representations can be either linear or nonlinear. Degrees of nonlinearity also exist for physical systems and optimization problems. For example, a quadratic programming optimization problem is only slightly nonlinear and is solved more easily than a more nonlinear problem, and sometimes yields a globally optimal solution (Chapter 4 discusses convergence of optimization methods).

Different S-O models have been developed and reported to solve all of the above optimization problem types for linear or nonlinear systems and representations. Generally, the different S-O procedures have been developed to improve computational performance for particular system-problem combinations. Some techniques require more computational, technical, or computer programming resources than others. This reference emphasizes methods practicably useful to most individuals developing pumping strategies and designs.

3.3.2  Linearity and Nonlinearity in Physical System and Simulators

3.3.2.1  Groundwater Flow

To reiterate, in equations describing linear processes, no variable can exist as a power other than one, and no products of two or more variables exist. Head response to pumping in a confined aquifer is linear. For a confined aquifer, Figure 3.1 illustrates Theis equation terms.

Note that if the potentiometric head remains in the semiconfining or confining layer, saturated thickness (and transmissivity), and storativity do not change in response to pumping. For a particular location, distance, and time, u and W(u) are independent of pumping Q. Thus drawdown s equals a linear product of Q and W(u)/4πkb—the head response to pumping is linearly proportional to pumping rate. This multiplicative property is a commonly employed characteristic of a linear system.

Figure 3.1   Theis equation and applicable physical setting.

Equation 3.12 illustrates use of the multiplicative property to compute groundwater head. Here, assume that the initial water table is horizontal and at equilibrium. Groundwater is extracted at a single well, index number ê.

3.12 $$\Delta h_{\widehat{o}}={\delta }_{\widehat{o},\widehat{e}}^h\frac{p_{\widehat{e}}}{p_{\widehat{e}}^{ut}}$$

where Δh_ô is a change in steady-state aquifer potentiometric surface elevation at observation location ô (L), ${\delta }_{ô,ê}^h$

is the influence coefficient describing effect of steady groundwater pumping at location ê on steady-state potentiometric surface elevation at location ô (L), p_ê is the pumping rate at location ê (L³/T), and $p_{ê}^{ut}$

is the magnitude of steady “unit” pumping stimulus in location ê used to generate the influence coefficient (L³/T). This does not necessarily equal 1.

Assume that a unit steady pumping extraction rate of 1 m³/min at well ê causes a drawdown of 1 m at observation point ô. In that case, ${\delta }_{ô,ê}^h$

equals (−1) and $p_{ê}^{ut}=1$

. Equation 3.12 shows that, ${\delta }_{ô,ê}^h$

if $p_{ê}^{ut}$

and are known, the change in head caused by any pumping rate can be easily computed. If pumping p_ê equals 2 m³/min, head change will equal (−1)(2)/(1) or −2. This linear response is typical of confined aquifers (or approximates behavior of unconfined aquifers where the change in transmissivity due to pumping is small by comparison with the original transmissivity).

Another characteristic of linear systems is the additive property, commonly taught as superposition and image well theory in groundwater texts (multiplicative and additive properties are employed more fully in Chapter 8). The additive property sums the effects of pumping at M^p locations on head at location ô. Equation 3.13 shows that this cumulative effect is the result of adding the individual effects of M^p pumping rates.

3.13 $$\Delta h_{\widehat{o}}={\displaystyle \sum _{\widehat{e}=1}^{M^p}{\delta }_{\widehat{o},\widehat{e}}^h\frac{p_{\widehat{e}}}{p_{\widehat{e}}^{ut}}}$$

where M^p is the total number of locations at which water is being pumped from the aquifer.

The additive property can be used to make composite influence coefficients. For example, it can describe the effect of pumping at M^p locations on the difference in head between two locations. The resulting Equation 3.14 is used in a later example.

3.14 $$\Delta {\Omega }_{\widehat{u}}={\displaystyle \sum _{\widehat{e}=1}^{M^p}{\delta }_{\widehat{u},\widehat{e}}^h}\frac{p_{\widehat{e}}}{p_{\widehat{e}}^{ut}}$$

where Ω _û is the difference in potentiometric surface elevation between locations 1 and 2 of pair û, (L). Here, since the initial steady-state potentiometric surface is horizontal, Ω _û also equals ΔΩ(û), the change in the head difference due to pumping.

Head response to stimuli in an unconfined aquifer is nonlinear because saturated thickness and transmissivity change in response to pumping. If the change is insignificant, the aquifer can be considered to be approximately linear. Historically, a general rule of thumb was that a system could be considered linear if the transmissivity change due to pumping was less than 10 percent. The smaller the initial transmissivity, the greater the head decline response to extraction.

Flow and head responses to stimuli in a river–aquifer system are represented in most numerical S models as piecewise linear functions. Figure 3.2 paraphrases how MODFLOW represents this and other head-dependent boundary flows using piecewise linear functions.

Different pieces of the functions can have different slopes. Figures 3.2a, b, and c have two, and three pieces, respectively. These functions are reasonable simplifications. At the beginning of a simulation iteration, MODFLOW evaluates head(s) in each cell and determines transmissivity and the piece of each flow equation to use for that cell. MODFLOW iterates until newly computed values match assumed values closely enough to satisfy convergence criteria. In other words, MODFLOW treats a nonlinear physical system as if it is linear, and iterates until convergence. Some S-O techniques also use an iterative procedure to address nonlinear system-problem combinations. Such methods assume linearity and iterate until they converge to a solution for a nonlinear situation. Other names for this iteration are cycling and successive linear optimization.

Figure 3.2   Sample piecewise linear head-dependent boundary flows in MODFLOW. (a) Seepage to drain, (b) river–aquifer seepage, (c) evapotranspiration. (Modified from McDonald, M.G., and A.W. Harbaugh, MODFLOW: A modular three-dimensional finite-difference ground-water flow model, USGS Techniques of Water-Resources Investigations, Washington, D.C., Book 6, Chapter A1, 1988.)

3.3.2.2  Groundwater Contamination Concentration

Groundwater concentration response to extraction stimuli is usually nonlinear, as illustrated in this groundwater contamination remediation. During the effort, we also empirically prove that S-O modeling can identify globally optimal solutions for nonlinear concentration response surfaces. The empirical proof required that we first determine what the response surfaces look like, and where the globally optimal solutions lie.

Figure 3.3 shows a plan view of initial trichloroethylene (TCE) contamination existing about 150 feet (46 m) beneath the ground surface. The management intent is to determine the optimal 4-year steady extraction rates from wells 1–3, which will cause the highest remaining concentration anywhere in the field (cmax) to be as low as possible. In Figure 3.3, the greatest initial concentration exceeds 125 ppb. The employed simulator is MT3DMS.

Figure 3.3   Initial TCE concentrations for three-well effort minimizing cmax after 4 years. (From Aly, A.H., and R.C. Peralta, Water Resour. Res., 35(8):2415–2425, 1999.)

Because the total pumping rate is fixed at 2,000 gpm, all possible strategies can be shown on a two-dimensional plane with axes for pumping at wells 1 and 2 (shaded plane in Figure 3.4). Pumping from well 3 is automatically known for any point (strategy) on the plane because it equals 2,000 −p ₁−p ₂. The feasible solution space for this optimization problem has two dimensions.

Figure 3.5 shows the cmax response surface that results from simulating all feasible pumping strategies. Note that the cmax surface is the same surface that an S-O model will navigate in searching for a globally optimal solution. The S-O model and the exhaustive simulations both yielded a globally optimal steady pumping strategy. This strategy will reduce contamination such that after 4 years, the highest concentration will be less than 1 ppb. Figure 3.5 contours show only one locality having the lowest concentrations. The contours do not reveal any other depressed areas that could contain locally optimal solutions.

Minimizing a state-variable cmax within even a single realization is a nonlinear optimization problem because cmax is a nonlinear function of pumping rates. Figure 3.5 problem is relatively straightforward because there are few choices concerning which wells will pump, and total pumping is constant.

Figure 3.4   Feasible solution space (shaded) for (P(1) + P(2) + P(3)) = 2,000 gpm.

Figure 3.5   cmax response surface after three years of pumping from three wells. (From Aly, A.H., and R.C. Peralta, Water Resour. Res., 35(8):2415–2425, 1999.)

3.3.3  Nonlinearity in Optimization Problem

A quadratic optimization problem is a special case of nonlinear optimization problem and is relatively easily solved. Nonlinear optimization problems are more challenging if the optimizer must select wells from among many candidate wells, and must add an extra significant cost to the objective function value if a well is pumped at all. The extra cost might be for construction, installation, special surcharge, or administration. Including the decision whether to construct is usually an integer or mixed integer optimization problem. Combining constrained nonlinear or objective function variables and MIP components produces MINLP problems—which are among the most challenging problems for optimizing pumping strategies.

3.3.4  Summary

The combined representations of the physical system and optimization problem determine how nonlinear and complex an optimization problem is. For large, complex optimization problems, especially MINLP and stochastic problems, reaching the globally optimal solution cannot be guaranteed. Except by testing solutions in all somewhat different parts of the feasible solution space, empirical certainty of global optimality cannot be achieved. In practice, that is acceptable. For complicated problems, getting close to global optimality using an S-O model will yield a better strategy than trial and error using an S model.

3.6.1  Multiple Module Models

Effective groundwater management sometimes involves multiple physical and nonphysical processes. An example of multiple physical processes includes groundwater extraction, seepage from leaky conveyance system and water use with attendant deep percolation, recharge, and chemical transport and reactions, and effects on sustainability. Chapters 7 through 9 discuss common S-O model simulator types, and how they can be coupled with optimizers. They cover the flow and transport processes most commonly considered in groundwater management. Discussing more processes is beyond the book scope.

Figure 3.7 illustrates multiple flow processes optimizable in an S-O model that could optimize urban and agricultural conjunctive water management without unacceptably affecting water quality. Solving the nonlinear S-O model would involve cycling techniques discussed in Chapter 4. Developing response or trade-off curves (Figure 3.8) involves multiobjective methods discussed in Chapter 6. Preparing Figure 3.8a involves employing constraints derived from a surface water quality model, within a dynamic stream–aquifer conjunctive use model. Creating Figure 3.8b involves detailed vadose zone flow and transport simulation. The curves can be included within groundwater management S-O model constraints.

Figure 3.7   Sample flow processes in a multiple module S-O model.

Figure 3.8   Sample objective response curves: (a) Maximum population supportable via conjunctive use, while protecting surface water quality, for alternative upstream inflow rates (After Ejaz, M.S., and R.C. Peralta, Adv. Water Resour., 18(2):67–75, 1995.); (b) Minimum crop yield reduction necessary to protect groundwater quality, as affected by depth to groundwater and irrigation frequency. (After Peralta, R.C. et al., Water Resour. Res., 30(11): 3183–3193, 1994.)

3.6.2  Decomposition and Hierarchical Optimization

Sometimes a management situation is so complicated and nonlinear that one cannot be certain of computing globally optimal strategies. Accepting this inability can be liberating. Acceptance frees one to look at the problem differently. Rather than solving one large complicated problem, one seeks solutions to several simpler problems that are somehow related. The adage of eating an elephant one bite at a time applies by analogy to digesting large or complex optimization problems.

Decomposition means the breaking of a large problem into smaller problems. Presumably the smaller problems are more readily solvable. The order in which the problems should be solved constitutes a hierarchy. The first subproblems are at the lowest level of the hierarchy. The final optimization to be solved is at the highest level—the master decision level. Partial or full repetitive iteration through the hierarchy might be necessary.

Figure 3.9 shows a sequence of models of different sizes used to aid critical area planning. The central image is a finite difference model grid of the Arkansas Grand Prairie. A hierarchical solution involves multiple steps, and coordinated assumptions concerning boundary conditions and flows. Information concerning S-O modeling constraints can go in the direction of the arrows, and in the opposite direction.

In the illustrated top-down approach, large-scale regional sustained groundwater yield and conjunctive use strategy would drive a developing Grand Prairie optimal strategy, which would determine pumping at individual wells. In actuality, the Grand Prairie model was developed first. Information concerning how much saturated thickness needed to be maintained in a critical subsystem was derived through daily irrigation scheduling and groundwater simulation of all pumping wells in the subsystem (assuring adequate water throughout the year). Then, developing a regional groundwater model could aid in assigning Grand Prairie boundary assumptions and coordinating conjunctive use with river system navigation and ecosystem requirements.

Optimal sustained yield and conjunctive use planning strategies for the Grand Prairie attempted to satisfy irrigation water needs (Figure 3.10). Figure 3.11 illustrates trade-offs that might be involved in selecting a strategy. Implementing one strategy would eventually cause the evolution of an acceptable equilibrium target potentiometric surface (Figure 3.12).

Figure 3.13 shows sustainable groundwater pumping rates that would maintain adequate saturated thickness throughout the area. However, water needs and groundwater pumping at that time, far exceeded those rates. So, imported surface water would be needed. Mathematical optimization or water-law based allocation simulation aids developing a transient conjunctive use strategy for the irrigation season (Table 3.2). This would use groundwater in accordance with the annual sustainable strategy.

Figure 3.9   Hierarchical sustained groundwater yield and conjunctive water use.

Figure 3.10   Upper bounds (qU) on water to be delivered to Grand Prairie cells when developing optimal strategies for a climatically average season. (From Ranjha, A.Y., Peralta, R.C. and A. Yazdanian, ICID Bulletin: Irrigation, Drainage and Flood Control, 39(1):1–12, 1990.)

Figure 3.11   Conceptual three-dimensional trade-off curve. (From Ranjha, A.Y., Peralta, R.C. and A. Yazdanian, ICID Bulletin: Irrigation, Drainage and Flood Control, 39(1):1–12, 1990.)

Figure 3.12   Target sustainable potentiometric surface in m above sea level. (From Ranjha, A.Y., Peralta, R.C. and A. Yazdanian, ICID Bulletin: Irrigation, Drainage and Flood Control, 39(1):1–12, 1990.)

Figure 3.13   Proportion of crop water use that can sustainably come from groundwater. (From Ranjha, A.Y., Peralta, R.C. and A. Yazdanian, ICID Bulletin: Irrigation, Drainage and Flood Control, 39(1):1–12, 1990.)

3.6.3  Multiple Study Area or Model Optimization

Figure 3.14 illustrates a setting for optimizing multiple groundwater simulation modeled study areas simultaneously. Because of the overlap in boundaries, the boundary conditions of both models must be interrelated and consistent. For example, one can maximize steady pumping from all wells in both areas, subject to normal constraints, plus constraints forcing both models to have the same variable head cell values and complementary flows in overlapping cells.

Figure 3.14   Plan view for matching boundary conditions when simultaneously optimizing two modeled groundwater areas.