B Usremax8 for Heterogeneous Multilayer Systems

1. Model Background

For optimizing management of complex heterogeneous systems, one would rather use US/ REMAX8 [16] than US/WELLSD. This is the basic version of the Utah State response matrix model. To develop influence coefficients, it uses code modified from MODFLOW, a modular finite-difference groundwater flow simulation model [4], and STR, a related stream routing module [17]. The physical system data needed by US/REMAXD can be input in the same format as is used by MODFLOW and STR. Internally, US/REMAX0 also uses a portion of PLUM AN, a decision support system for optimal groundwater contaminant plume management [18], and other code.

The optimization model formulation capabilities are similar to those of US/WELLSD (Table 2). For steady state, the generic objective is to minimize

where Wj is the weight assigned to pumping in cell j, dimensionless or [$T/L3].US/REMAXB can employ constraints 1-3 of US/WELLS° for multiple layers. Similar to the US/WELLSD constraint 4, US/REMAXB can force total extraction to exceed, equal, or be less than total injection. Again, via the sign on the weighting coefficients, one can perform maximization. One can also achieve multiobjective optimization by the weighting method. Whereas in US/ WELLS15 the same weight must be applied to all extraction wells in a time step (and a different weight can be used for injection wells, but the same must be applied to all such wells in a particular time step), in US/REMAXB each well can employ a different weight.

2. Application and Results

Introduction. For illustration, we discuss addressing a contaminant plume in a representative study area. First, the study area is described and the results of continuing current management are predicted, using MODFLOW+STR for flow simulation and MOC [19] for transport simulation. Then an approach to developing an optimal strategy is discussed, the S/O model is applied, and an optimal strategy is computed. Next, the system response to implementing the optimal strategy is verified using MODFLOW+STR and MOC. Finally, slight variations in the management goal or situation are assumed and new optimal strategies are developed. Computed optimal strategies are compared. Suguino [20] first addressed this study area using PLUM AN. Some of the discussion below follows his development.

Study Area Description and Situation. The area (Figure 3) measures about 4.3 km by 4.3 km. It is bounded on the north by a large saltwater body; on the south, east, and northwest by impermeable material; and on the west by a lake. A river transects the area from south to north. Aquifer parameters of this example study area were obtained from ranges reported by Todd [21].

For the unconfined upper layer (layer 1), parameters are as follows. Hydraulic conductivity:

1st zone: 45 m/day (coarse sand) from lake to contaminant spill area (columns 1-36 and 57-58)

2nd zone: 30 m/day (medium sand) in irrigated area (columns 51-56) 3rd zone: 450 m/day (fine gravel) in contaminant spill area (columns 37-50).

Specific yield:

1st zone: 0.27 (coarse sand) 2nd zone: 0.28 (medium sand) 3rd zone: 0.25 (fine gravel)

Figure 3 Finite-difference grid for the area addressable with US/REMAX8.

Recharge by deep percolation and/or irrigation:

1.167 x 10~8 m/sec in nonirrigated area 1.928 x 10~8 m/sec in irrigated area

In the confined lower layer (layer 2):

Transmissivity: 0.1564 m2/sec Saturated thickness: 30.0 m Storage coefficient: 0.0001

Finite-difference models are to be used in this study. This requires system discretization. The resulting block-centered cell grid (Figure 3) has 58 columns and 39 rows. Cell side lengths range from 3 to 400 m. Because MOC will be used for transport simulation near the plume, cells of uniform size are specified for that region. The resulting 17 row by 20 column region (subsystem) near the plume has square cells of 15.2 m (50 ft) side length.

A conservative (nonreactive) contaminant is assumed to be spilled in the top aquifer layer (layer 1) of cell (22, 18) or (lls, 3S). (The subscript "s" after a cell row or column index indicates that the cell is in the subsystem.) This cell is treated as a continuous source during the management period.

Initially, pumping for water supply occurs in two cells between the plume and the river. One well is in layer 1 of (23, 15) or (12s, 15s). The other well is in layer 2 of (18, 18) or (7S, 18s). There is immediate concern about the potential for contamination reaching the supply well in layer 1.

Nonoptimal System Response Determination (Step 1). Before one attempts to develop an optimal strategy, one usually demonstrates the need for such a strategy. This requires predicting system response if no optimal strategy is implemented. Frequently, simulation models are used for this action. Here, MODFLOW+STR computes the potentiometric surface that will result from assumed steady-state conditions (Figure 4).

Because of the gradient, the contaminant will tend to migrate toward the supply wells. MOC is used to quantify the migration resulting in the subsystem from the steady flow. Figure 5 shows the 210 ppb contour expected to result 60 days after contamination begins. Furthermore, concentration in the cell containing the drinking well (12s, 15s) reaches 317 ppb 8 months after the spill. We assume that this concentration level exceeds the health advisory for human consumption and that developing a plume capture strategy is desirable.

Management Goals Specification and S/O Model Formulation for Scenario 1 (Step 2). The assumed goal is to minimize the steady pumping (extraction and injection) needed to capture the plume. Plume capture will presumably be achieved when hydraulic gradients, just outside the plume boundary, all point toward the plume interior. We also want the head at extraction wells not to drop too far (to avoid reducing saturated thickness by more than about 10%) or the head at injection wells not to rise to the ground surface. These criteria identify the example problem termed Scenario 1.

The S/O model formulation for this scenario is shown below. The model computes the pumping strategy that minimizes the value of the objective function, subject to the stated constraints and bounds. Locations of potential injection and extraction wells to be considered by the model are shown in Figure. 5. Figure 6 identifies head difference (gradient) control cell pairs and shows the direction that will be imposed on the hydraulic gradient by any computed optimal strategy. These are placed to enclose the plume projected to exist by day 60. A modeler can select potential well locations on the basis of practical experience. For example, the closer the injection wells are to the head gradient control locations, the less pumping is needed to

Figure 5 Subsystem discretization, potential well locations for Scenario 1, and 210 ppm contour, 60 days after contamination begins.

Figure 4 Nonoptimal (unmanaged) steady-state Potentiometrie surface contour map for the study area of Scenario 1 (meters above MSL).

o iI

Figure 4 Nonoptimal (unmanaged) steady-state Potentiometrie surface contour map for the study area of Scenario 1 (meters above MSL).

Legend:

[♦1 unmanaged well ( 1 ) layer 1 (21 layer 2 rn potential managed injection tuell [Ô] potential managed extraction tuell

| contaminant source Q 210 ppb concentration contour meter 0 IS 30 45

ES23

Figure 5 Subsystem discretization, potential well locations for Scenario 1, and 210 ppm contour, 60 days after contamination begins.

10 il 12

16 17

X

X

X

X

X

X

X

X

X

1

1

1

1

1

1

1

^

X

• (2)

I

I

I

I

I

Is

X

0

X

0

0

X

ô

0

4

X

o

X

(1)

t

t

Î

t

t

t

ÎN

X

1

1

1

1

1

I

I

X

X

X

X

X

X

X

X

X

25 27 J

29 31 33 35

16 17 t8

20 21 22

Legend;

( I ) layer 1 (2) lager 2 ri potential managed ■ Injection ujell (upper bound on head)

[Ô] potential managed cHtractian well Oomer bound on head)

| contaminant source «— gradient constraint 0 15 30 45 m

16 17 19 21 23

25 27 J

29 31 33 35

Figure 6 Head-difference constraint locations applied within the S/O model in Scenario 1.

satisfy the head-difference constraint. Thus, the modeler might want the model to consider pumping sites near the location where heads need most to be affected.

The model objective is to minimize the value of Equation (3), using weights of 1, subject to

Gà > 0.01, hi > 15.0, hè < 25.0, for ô = 1, . for ê = 1, . for e = 1, .

where G-0 is the difference in head between a pair of cells, the first located farther from the plume. A positive value denotes a higher head farther from the plume, [L]; he is the hydraulic head just outside the casing of pumping well e located in the center of a pumping cell, [L]; o is the index denoting pair of cells head-difference (gradient) control pair; and e is the index denoting pumping well at the center of cell j or k. Here j = 6 and k = 25.

Note that identifying the location of potential extraction and injection wells for the model (Figure 5) does not mean that the model will choose to pump at those locations. Via the optimization process, the model might choose to pump at only a few of the potential sites. The computed strategy will require less total pumping than any other strategy possible for the specified potential well locations and imposed bounds and constraints. Furthermore, since this is a steady-state problem, steady-state system response to implementing the strategy computed by the model will satisfy all those bounds and constraints. This is verified in the next step.

Optimal Strategy Computation and Verification for Scenario 1 (Step 3). The optimal strategy computed for Scenario 1 is shown in Table 3. Because the model is minimizing pumping only for plume containment in layer 1, no extraction is shown for layer 2. The original unmanaged pumping does continue from original supply wells in both layers (Figure 3) but is not included in Table 3 because the model is not optimizing that pumping.

Table 3 Pumping Results for the Sample Scenarios

Scenario g(extr)

Constraints

1st Layer 2nd Layer

Gradient constraint on heads located on the same layer, head constraint on injection and extraction well.

Added pumping constraint: total sum of extraction = total sum of injection.

Gradient constraint on heads located on the same and on different layers, head constraint on injection and extraction wells.

Figure 7 shows the locations of wells that will pump, according to the optimal strategy. It also shows the head-difference constraints [Equation (4)] that will be tight. Tight constraints are those that are satisfied exactly. The other gradient constraints are also satisfied, but the model had no difficulty in doing so. These latter head-difference constraints are "loose" (there is more than 0.01 m difference between the heads at the two cells coupled by an arrow in Figure 6 but not shown at all in Figure 7). No heads are against their bounds. Therefore neither Equation (5) nor (6) is tight.

It is appropriate to verify that the computed strategy accomplished its goal of plume capture. MODFLOW+STR can be used to demonstrate how quickly the optimal steady pumping s

X

X

1

1

1

/

X

121

I

I

J.

j

X

O

O

O

X

(11

X

20 21 22

Legend:

Ml actiue managed injection well

[ÔH actiue managed entraction well

I | contaminant source tight gradient constraint meter 0 15 30 45

Figure 7 Location of optimal pumping wells and tight head-difference (gradient) constraints for Scenario 1.

strategy will cause the desired gradients to occur. Transient simulation demonstrated that the gradient constraints would be satisfied 30 days after implementing the optimal pumping strategy (Figure 8). Figures 9 and 10 show the ultimate steady-state surface resulting from strategy implementation. Clearly, a groundwater divide has been formed between the plume and the supply well.

MOC is used to predict the pollutant transport that would result from strategy implementation. No contaminant moved past the injection wells.

Theoretical verification of the optimality of the computed strategy is beyond the scope of this document. However, many texts on operations research and linear programming assure the optimality of solutions to models having a linear objective function and constraints.

Alternative Scenarios.

Scenario 2. This scenario differs from the previous in the addition of a constraint forcing total injection to equal total extraction around the plume. Again, pumping from the two supply wells is not included in the total.

aquifer optimal all gradients check of contain, pumping constraints cont.conc.

begins achieved

0 60 90 240 days

Figure 8 Time scale of Scenario 1.

Figure 9 Subsystem Potentiometrie surface resulting from implementing the optimal pumping strategy for Scenario 1 (meters above MSL).

Well

Ips

gpm

llfel)

Ips

gpm

1

1.36

21.56

6

-3.09

- 49.02

2

0.62

9.89

7

5.02

79.58

3

-0.60

-9.49

8

4.46

70.67

A

-9.69

- >53.55

9

7.88

124.96

5

O.BS

13X6

10

—100.00

-IS84J82

Figure 10 Subsystem potentiometric surface resulting after 6 months of optimal pumping for Scenario 1 (meters above MSL).

Results in Table 3 show an increase in extraction and a decrease in injection. Total pumping needed for plume containment increased slightly (1.4%). This illustrates the phenomenon—increasing the number or restrictiveness of constraints does not improve the value of the objective function.

Although total pumping increased, one less extraction well is used in this strategy than in the previous (Table 4). The same number of gradient constraints are tight, but the locations of the tight gradient constraints differ slightly.

Scenario 3. This scenario demonstrates what might happen if involved managers have conflicting goals. It differs slightly from Scenario 1. In addition to controlling the plume, the agency wishes to extract more from layer 2 for water supply. Three new potential extraction wells are located in cells (19, 25), (20, 25), and (21, 25), as if along a nearby road. Pumping is not permitted to change at the two initial supply wells.

Table 4 Numbers of Managed Wells that Will Pump Under the Optimal Strategies for the Tested Scenarios

Scenario

S(extr.) 1st Layer 2nd Layer

K inj.)

(g + b) total

1

3 —

6

9

2

2 —

6

8

3

4 3

14

21

As a result, the objective function is altered to maximize new extraction from layer 2 while still minimizing the pumping in layer 1 needed to capture the plume. This is achieved by assigning a negative sign to extraction from the supply wells, and minimizing:

Since minimizing a negative number is the same as maximizing a positive number, minimizing negative extraction in layer 2 means maximizing that extraction.

Also added are new constraints imposed on vertical flow in cells (21, 21) and (22, 21). There, the head in the lower layer is forced to exceed that in the upper layer by 0.01 m, preventing the downward migration of contaminant.

Figure 11 shows the resulting optimal injection and extraction well locations and tight gradient constraints. The optimal pumping strategy includes seven extraction wells and 14 injection wells. Although extraction of polluted water decreases, injection increases with respect to Scenario 1 (Table 3). Extraction of water for public supply increases by 31% above the un-managed rate.

Although the gradient constraints are all satisfied by the optimal strategy, subsequent simulation demonstrated that the vertical gradient is reversed in some plume-containing cells in which the gradient was unconstrained. This illustrates that one must be careful in placing head or gradient control in appropriate locations. In practice, another optimization would be performed, using additional vertical head-difference or gradient constraints. Processing Considerations. It is useful to consider the resources required to address optimization problems. First, the total computer time needed to solve an optimization problem is of

1 2 3 4 5 6 7 8 9 10 1112 B 1415 1617 18 19 20

X

X

X

X

X

X

I

I

1

1

1

1

/

X

• (2)

i

1

I

1

I

Iv

o

X

O

X

+4

o

+

O

X

X

X

(1)

X

X

20 21 22

Legend:

nn actiue managed injection tuell fp] actiue managed eKtraction tuell

PI contaminant source

(heads on same layer) « tight grad. constr. (heads on diff. layer)

meter 0 15 30 45

Figure 11 Location of optimal pumping wells and tight head-difference (gradient) constraints for Scenario 3.

interest. Table 5 illustrates the time needed to address Scenario 1. Included stages use either the discussed simulation models or the PLUMAN code on a 386 PC running at 33 MHz and having 4 MB RAM. Time required for US/REMAX8 is comparable to that of PLUMAN, since it uses many of the same solution procedures.

Clearly, the stage of computing influence coefficients, arranging the optimization model, and calculating an optimal strategy is the most computationally intensive. For this scenario and stage, two steps can be distinguished. The first involves computing influence coefficients. The second is model organization and optimization problem solution.

Here, the step of generating influence coefficients requires by far the most time. This results because this act essentially involves one simulation of a modified MODFLOW+STR per potential pumping location. Since there are 31 potential pumping locations, 31 simulations are performed to develop the influence coefficients needed for the response matrix. The more decision variables (potential pumping rates), the more computer time involved in this step.

The step involving model formulation and calculation of the optimal strategy is fairly short. The time needed to perform the optimization is a function of the number of decision variables (potential pumping rates) and state variables (heads or gradients that must be constrained within the optimization model). The larger these numbers, the more time required.

Second, the size of the optimization problem being solved is of interest. For example, the special versions of US/WELLSd and US/REMAXB that are released in shortcourses are limited in the number of nonzero values they can have in the optimization formulation. (Even optimization algorithms that are not part of water management models are commonly limited either in the number of nonzeros or in the number of rows and columns in their constraint equations.)

By way of explanation, there is one row in the response matrix per head or gradient constraint equation per time step of constraint. There is one column in the matrix per decision variable. For a steady-state problem, total matrix size is the product of the number of control locations and the number of decision variables. The matrix contains one nonzero coefficient for each potential pumping location-head control location pair (per time step of active constraint).

For the steady-state Scenario 1, there are 31 x (22 + 31), or 1643, nonzeros due to influence coefficients. There are also 31 nonzeros due to the weighting coefficients (even if they are 1 in value) assigned to decision variables in the objective function. Thus, the optimization model formulation for Scenario 1 employs almost 1700 nonzeros. (That of Scenario A using US/WELLSd includes 919 nonzeros.) This number can be reduced significantly by considering injection in only every other cell on the plume periphery rather than in each cell. For example, if only 12 injection wells were considered, the number of nonzeros would be about 18 x (22 + 18) + 18, or 738. In addition to reducing problem size, this would significantly reduce computational time.

Table 5 Computer Time Required to Perform Each Activity for Scenario 1

Step

Software used

Time (min)

1

MODFLOW+STR (compute nonoptimal head)

5.0

2

MOC (predict solute transport in a nonoptimal potent surface)

35.0

3

PLUMAN (compute influence coefficients, formulate management model

150.0

and determine optimal pumping strategy)

4

MODFLOW+STR (compute transient head response to optimal pumping)

1.3

5

MOC (compute head and solute transport response to optimal pumping)

8.0

Reducing the number of nonzeros below 1000 is important because that is the upper limit on problem size in the inexpensive "special" versions of US/REMAXB and US/WELLSD. If problem size increases beyond that, software price increases dramatically. The full professional versions of the software can address problems of virtually unlimited size.

Was this article helpful?

0 0

Post a comment