Teaching Numerical Groundwater Flow Modeling with Spreadsheets: Unconfined Aquifers and Multilayered Vertical Cross-Sections

Abstract

The use of spreadsheets in numerical groundwater flow modeling, although not new, has yet to be fully exploited in educational settings. This teaching aid complements the one published earlier, broadening its scope from the numerical solution for a single-layer confined aquifer to include unconfined aquifers. It covers flow in horizontal layers and multilayered vertical cross-sections. The new spreadsheets can also handle hydraulic conductivity anisotropy, provided its principal directions align with the Cartesian axes. It can also manage irregular aquifer shapes, incorporate all three types of boundary conditions, and account for pumping and recharge in the modeling process. However, the flow regime remains at steady state. With the integration of these new features, the spreadsheet model becomes a versatile tool for simulating a great number of steady-state groundwater flow situations. This user-friendly platform enables students to swiftly experiment with diverse configurations, delve into the underlying mathematical expressions, and gain insight into the workings of more sophisticated numerical codes like MODFLOW. As a complement, and for verification purposes, Python Flopy scripts are provided for all examples in the GitHub repository where the spreadsheets can be accessed.

1 Introduction

This teaching aid is a companion to the one published by the first author in which the solution of the partial differential equation describing steady-state flow in a vertically averaged single-layer confined aquifer using a spreadsheet was presented (Gómez-Hernández 2022). It is suggested the reader looks at that paper to get familiar with the basic architecture of the workbooks since it will be the same one used in this work.

The application of spreadsheets for solving the groundwater flow equation is not a recent development. The pioneering work in this area can be attributed to Olsthoorn (1985), who, in his Ph.D. thesis (Olsthoorn 1999), reports having extensively worked on groundwater flow modeling using the Lotus 1, 2, 3 program as early as 1983. Since then, numerous research papers have been published that utilize spreadsheets for both steady-state and transient groundwater modeling (Akhter et al. 2006; Ankor and Tyler 2019; Bair and Lahm 2006; Bhattacharjya 2011; Elfeki and Bahrawi 2015; Fox 1996; Karahan and Ayvaz 2005a, b; Li and Davis 2022; Molano 2014; Niazkar and Afzali 2015; Olsthoorn 1999; Ousey 1986).

In this paper, emphasis is placed on how the workflow described in Gómez-Hernández (2022) can be extended to solve groundwater flow in more generic, complex settings. Specifically, the following two cases will be addressed: unconfined flow in a single-layer horizontal aquifer and unconfined flow in a vertical cross-section of a multilayered aquifer. In addition, every cell in the model can be anisotropic to flow as long as the principal directions of the anisotropy conductivity tensors are parallel to the Cartesian axes. Examples of the three common types of boundary conditions are included, plus the difficult-to-handle phreatic surface boundary condition in unconfined flow. In all cases, flow will be at steady-state. All the examples have been prepared using Excel spreadsheets, and they have been tested in Google spreadsheets. The main difference between the two implementations is the superior speed of Excel versus Google.

The paper starts with a review of the same aquifer analyzed in the earlier paper but considers the aquifer as unconfined. This “minor” conceptual change has a significant impact on the numerical solution of the flow problem. The paper describes how this change is addressed, and the findings are carried over to the solution of flow in a vertical cross-section of an unconfined aquifer. The paper continues showing how this latter solution can be used to analyze seepage through an earth dam and ends with a discussion about possible extensions and limitations. The paper also ends with a link to the GitHub repository where the spreadsheets and Python Flopy scripts to verify the solutions can be downloaded.

2 Two-Dimensional Unconfined Flow in a Horizontal Plane

The paper starts with the same aquifer used in the previous paper depicted in plan view at the top of Fig. 1. It is a quasi-rectangular unconfined aquifer that is crossed by a river from north to south, impermeable in all its boundaries except for the east boundary, where it is in touch with a lake to which is perfectly connected. There are three pumping wells, and there is also infiltration over the entire aquifer from the top. Aquifer hydraulic conductivity takes three different values in three zones as shown in the figure. It is assumed that the aquifer is at steady state.

Fig. 1

Aquifer sketch used for both the confined and unconfined aquifer simulations. The difference is in the cross-sections

Fig. 2

Aquifer discretization

The equation that gives the hydraulic head for a generic cell and its derivation can be looked up in Sect. 2 of Gómez-Hernández (2022), and it will not be repeated here. In the next section, the specific implementation details that make the solution for the unconfined aquifer depart from that of the confined one are explained in detail. The reader who wishes to go directly into the use of the spreadsheets can jump to the corresponding “workflow" subsections.

2.1 Implementation Details

Figure 1 provides the same plan view of the aquifer as in previous work and includes two cross-sections: the upper section represents the confined aquifer, while the lower section represents the unconfined aquifer. In the case of the confined aquifer, the saturated thickness is constant and not affected by the piezometric head. As a result, the transmissivity—calculated as the conductivity multiplied by the saturated thickness—can be determined in advance and was included as part of the input data. Conversely, for the unconfined aquifer, the saturated thickness varies depending on the piezometric head and equals the difference between the piezometric head and the aquifer bottom. This means transmissivity cannot be predetermined and included in the input data; instead, conductivity will be provided, and transmissivity computed. Figure 2 shows the discretization used for the model, the same one used in the previous teaching aid: 19 rows by 33 columns of 100 -m by 100 -m cells.

The dependency of transmissivity on the saturated thickness, which in turn, depends on the piezometric head, yields the problem highly nonlinear, and the convergence of the iterative solution unstable, particularly because during the iterations, some cells may become dry, and there is a need to account for a possible rewetting of the cells in the next iteration if the surrounding cells are not dry. Acknowledging this situation brings the following complications:

• 1. Input data must be hydraulic conductivity, not transmissivity.

• 2. At each iteration, the cell saturated thickness must be computed for each cell as a function of the piezometric head in the previous iteration.

• 3. At each iteration, transmissivities must be recomputed as the product of hydraulic conductivity and saturated thickness. This makes the solution particularly sensitive to the choice of initial heads to start the iterations because initial values very far from the solution will result in unrealistic transmissivities that will result in unrealistic head estimates that might result in a lack of convergence.

• 4. Because the saturated thickness may become zero if the hydraulic head estimate is below the cell bottom, a mechanism must be established to deactivate a cell when its saturated thickness becomes zero.

• 5. At the same time, a mechanism must be established to reactivate (or rewet) an inactive cell when the hydraulic head estimates in the surrounding cells are above the aquifer bottom of the inactive cell.

Recall that the numerical solution is iterative and that the recalculate option of Microsoft Excel is used to iterate on the cells that include circular references. We learned, the hard way, that it was not sufficient to add new spreadsheets for the calculation of saturated thickness and transmissivity. We had to enforce the order of calculation, that is, first, the saturated thickness then the transmissivity, but it was not obvious which the calculation order of the spreadsheets was. We found out that computations are made in the alphabetic order of the spreadsheet name (Williams 2001), which forced us to rename all the spreadsheets as described below.

The changes needed with respect to the confined aquifer in the previous paper are as follows:

• 1. Renaming all spreadsheets to ensure that their alphabetic order coincides with the order in which they have to be computed.

• 2. Adding a spreadsheet with the cell bottom elevations.

• 3. Adding two spreadsheets for conductivity data, one with the values of the conductivity tensor along the x axis, which will serve to compute the transmissivities $T_W$ and $T_E$, and another one with the values of the conductivity tensor along the y axis, which will serve to compute the transmissivities $T_N$ and $T_S$.

• 4. Adding a spreadsheet to compute the saturated thickness.

• 5. The transmissivity spreadsheets become calculated spreadsheets that are computed after the saturated thickness.

• 6. During the iteration process, the calculated heads may result below the cell bottom, rendering it dry, but in the next iteration, if the surrounding cell heads are above the cell bottom, it is indicative that the cell should be partially saturated. Therefore, there is a need for a mechanism that allows the drying and rewetting of any cell. For this purpose, two new spreadsheets are needed, one containing a flag indicating whether a cell should be wet as a function of the heads in the surrounding cells, and another one containing the head value that should be assigned to cells that were dry but should be wet because the surrounding cells are above the cell bottom.

• 7. Because the cells may dry out and become inactive during the iteration process, an intermediary binary spreadsheet with the active cells after the last iteration is also needed.

• 8. Infiltration due to recharge is only calculated for the wet cells.

Now let us revise the new collection of spreadsheets and the changes introduced in the definition of their formulas.

Final results (six sheets) The final results contain the same ones as for the unconfined aquifer, with the addition of a spreadsheet where the flow balance is shown for each cell. This new spreadsheet serves to debug the functioning of the model since it pinpoints the cells in which the conservation of mass may be failing.

• B_h: final results with the hydraulic heads in the aquifer [L]

• F_QNorth: flows through cell north boundaries [L$^3$ T$^{-1}$]

• F_QSouth: flows through cell south boundaries [L$^3$ T$^{-1}$]

• F_QWest: flows through cell west boundaries [L$^3$ T$^{-1}$]

• F_QEast: flows through cell east boundaries [L$^3$ T$^{-1}$]

• G_CellBal: residual flow in the active cells (it should be zero, if convergence has been achieved) [L$^3$ T$^{-1}$]

Input parameters (nine sheets) A new sheet is added with the bottom elevation, and the transmissivity sheet is replaced by the hydraulic conductivity tensor components along the principal directions, which are assumed parallel to the Cartesian axes.

• A_i: active cells

• A_Bot: bottom elevation [L]

• A_hfix: prescribed heads [L]

• A_Kx: conductivities along the x axis [L T$^{-1}$]

• A_Ky: conductivities along the y axis [L T$^{-1}$]

• A_W: pumping well extraction [L$^3$ T$^{-1}$] (negative if injection)

• A_hR: river stage [L]

• A_hB: elevation of riverbed bottom [L]

• A_R: riverbed conductance [L$^2$ T$^{-1}$]

Intermediate variables (14 sheets) There are eight more sheets than for the confined case. All 14 sheets are described next.

• C_i: active cells after the last iteration

• C_Sat: saturated thickness [L]

• D_Wet: Boolean variable indicating whether a cell should be wet or not as a function of the hydraulic heads in the surrounding cells

• D_ReWet: value of the hydraulic head to be used only for those cells that became dry but should rewet because D_Wet is True

• D_Tx: transmissivity along the x axis [L$^2$ T$^{-1}$]

• D_Ty: transmissivity along the y axis [L$^2$ T$^{-1}$]

• D_Qriv: river recharge flow [L$^3$ T$^{-1}$]

• D_QN: infiltration flow [L$^3$ T$^{-1}$]

• E_TN: transmissivity at the interface with the north cell [L$^2$ T$^{-1}$] computed as the harmonic average of adjacent y transmissivities

• E_TS: transmissivity at the interface with the south cell [L$^2$ T$^{-1}$] computed as the harmonic average of adjacent y transmissivities

• E_TW: transmissivity at the interface with the west cell [L$^2$ T$^{-1}$] computed as the harmonic average of adjacent x transmissivities

• E_TE: transmissivity at the interface with the east cell [L$^2$ T$^{-1}$] computed as the harmonic average of adjacent x transmissivities

• E_TZ: sum of E_TN, E_TS, E_TW, and E_TE. (Note that the name has been changed from TT to TZ to ensure that it is computed after the four summands have been updated.)

Each spreadsheet includes the spatially distributed values of the respective variable. The user must choose a consistent set of units to describe each of the input parameters; the intermediate results and the final results will be in consonance with that decision. In the example shown herein, meters and days have been chosen as the basic units to represent all variables.

Next we describe how the workbook is set up, indicating only the differences with respect to the setup already described in the previous paper:

Step 0. Create workbook The only difference is that there is a need for more spreadsheets as listed above and that their names have to be changed to ensure proper calculation of the intermediate variables.

Step 1. Define active cells and cells with prescribed heads No change. The active cells and the prescribed cells are shown in Figs. 3 and 4.

Fig. 3

Sheet with active cells

Fig. 4

Sheet with prescribed head cells and their values

Step 2. Input data The only difference with the previous paper is the need to introduce the cell bottom elevation and the hydraulic conductivities, which can be anisotropic. For the bottom elevation, we have chosen an undulating profile with elevations that range between 68 m and 80 m. We have also included a bump in the aquifer bottom towards the center as high as 100 m to demonstrate the drying capabilities of the new workbook. Figure 5 shows the bottom elevation for each cell.

Fig. 5

Sheet with bottom elevations. Notice the promontory in the center of the lower half that will end up with dry cells

The hydraulic conductivities followed the same zonation used for the confined case, although their values have been modified so that the calculated transmissivities are of the same order of magnitude. Conductivities are anisotropic, with 10 times higher conductivity along the x axis than the y axis. Figures 6 and 7 show the input values.

Fig. 6

Sheet with hydraulic conductivities along the x axis

Fig. 7

Sheet with hydraulic conductivities along the x axis

Wells are located in the same cells as in the previous paper; however, the pumping rates have been lowered from a total of 35,000 to 23,000 m$^3$d$^{-1}$ distributed in three wells pumping at 5,000, 10,000 and 8,000 m$^3$d$^{-1}$ as shown in Fig. 8. The reason for this reduction is that with the previous values, the aquifer would go dry at the well locations.

Fig. 8

Sheet with well pumping rates

The remaining input data for the active cells, prescribed heads, river stage, elevation of riverbed bottom, and riverbed conductance remain as in the previous paper.

Step 3. Compute intermediate variables Given the heads for the last iteration, there is a need to compute a number of intermediary variables before the next iteration.

Step 3.1. Active cells at the end of the iteration First, check if any active cell has dried out. It will occur if the calculated head is below the cell bottom elevation. This check is performed in spreadsheet C_i. As for all spreadsheets, a unique expression is used for all cells within the rectangle enclosing the aquifer. This expression for cell I5 is

It reads as follows:

• If we are initializing the calculations or restarting them for any reason, then the input information about active cells is used; else

  • If the head at the current cell is “hundef" (the default value for initially inactive cells, and the value for cells whose transmissivity has become zero because its saturated thickness is zero), then set the cell as inactive; else

    • Set the cell as active.

This expression is copied within the range I5:AO23 enclosing the aquifer. The result is a spreadsheet with the active cells after the last iteration, which are all active cells in the input spreadsheet minus those cells that have dried out. From now on, the expression for cell I5 will be reported without making reference to it. Additionally, it will not be reiterated that this expression needs to be copied to all other cells in the aquifer.

Step 3.2. Saturated thickness at the end of the iteration The saturated thickness, computed in spreadsheet C_Sat, is set to zero if the cell was defined as inactive in the input data or if the head after the last iteration is below the cell bottom; otherwise, the saturated thickness is set to the difference between the piezometric head and the cell bottom. This is achieved with the following expression

Step 3.3. Decide if a cell should be rewetted in the next iteration Sheet D_Wet is built specifically to determine whether a cell should be rewetted in the next iteration. Rewetting should happen if any of the active surrounding cells has a piezometric head above the bottom of the current cell. This is controlled with a Boolean flag that is computed with the following expression

This expression reads as follows:

• If the cell is inactive in the input data, then it should not be wet; else

  • If any of the four adjacent cells is active and its head is above the cell bottom, it should be wet in the next iteration.

This calculation is performed for all cells in the aquifer, although it will be used only for those cells that are dry in the previous iteration, as explained later.

Step 3.4. Compute the head value to be used in those cells that should rewet When a cell has to move from dry to wet, the head assigned to that cell is computed as a function of the heads in the surrounding cells. The procedure is similar to the one used in the program MODFLOW (Harbaugh 2005): a “WetFactor” is defined in spreadsheet D_ReWet, which is used to compute the rewetting head. If all surrounding cells are active and have piezometric heads that are above the cell bottom, the new head is set to the cell bottom plus the average of the differences between the heads in the surrounding cells and the cell bottom multiplied by the WetFactor; else, the surrounding heads are searched in the order top, bottom, left, right, the first one that is active and with a head value above the cell bottom is chosen, and the rewetting head is set to the cell bottom plus the difference between the adjacent cell head minus the cell bottom multiplied by the WetFactor. This description translates to the following expression

The objective of this calculation is to set the value of the rewetted cell above its bottom by a fraction of the difference between the bottom and the surrounding cell heads, allowing cells that went dry during an iteration to rewet if the surrounding cells have piezometric heads that justify such a rewetting. In this example, the WetFactor is 0.01 (specified in spreadsheet D_ReWet).

Step 3.5. Computing transmissivities After computing the saturated thickness, transmissivities must be computed as the product of the hydraulic conductivity and the saturated thickness. This is done for both the x and y components of the transmissivity tensor. For instance, the transmissivity in the x direction in D_Tx is given by

Notice that the transmissivity is set to zero if the cell, after the last iteration, is inactive (because it went dry or was inactive in the input data). The same procedure is repeated for the y direction in D_y.

Step 3.6. Compute the flow entering the aquifer from the river It uses the same equation as in the previous paper with the only difference being that it is limited to the cells that are active in the current iteration, a condition that has been computed in step 3.1 and corresponds to those cells with a value of 1 in spreadsheet C_i.

Step 3.7. Compute the flow entering the aquifer through infiltration Infiltration could be an input array as the rest of the input arrays, but we have decided to compute it as a function of a constant net infiltration given in m d$^{-1}$ and the cell size in m. Both values are specified in spreadsheet D_QN. The expression used takes into account that the infiltration, for the purpose of computing the global balance later on, should be calculated only on the active cells in the last iteration. For this reason the expression has suffered a small modification with respect to the one in the previous paper as follows

For the inactive cells, the recharge is undefined, and for the cells with a prescribed head, the recharge is set to zero; otherwise, the recharge is given by the product of the infiltration rate N1_ and the cell area Delta*Delta.

Step 3.8. Compute the interblock transmissivities The interblock transmissivities are computed as in the previous paper as the harmonic average of adjacent cell transmissivities, with the only difference being that now for the computation of the horizontal interblock transmissivities in E_TE and E_TW, the transmissivities along x are used, and for the interblock transmissivities in E_TN and E_TS, the transmissivities along y are used. Similar expressions are written for the other three interblock transmissivities. In addition, to speed up the calculations, sheet E_TZ is constructed with the sum of the four interblock transmissivities.

Step 4. Iterative solution for piezometric heads The calculation of the heads is performed in sheet B_h and takes advantage of the iterative capabilities of Excel when dealing with circular references. The equation that has to be met by the head at any given active cell was derived in the previous paper and repeated here as

$$\begin{aligned} h=\frac{T_N'h_N+T_S'h_S+T_W'h_W+T_E'h_E-W+N\varDelta ^2+Q_{Riv}}{T_N'+T_S'+T_W'+T_E'}, \end{aligned}$$

(1)

where h is the head at the cell being computed, $T_N', T_S', T_W', and T_E'$ are the interblock transmissivities to the north, south, east, and west of the cell, respectively, $h_N, h_S, h_W,$ and $h_E$ are the heads at the cells to the north, south, east, and west of the cell, respectively, W is the extraction flow at the cell, N is the net infiltration rate, $\varDelta ^2$ is the cell area (assuming square cells), and $Q_{Riv}$ is the flux entering from the river.

We recall that the initialization and beginning of the calculations are controlled by two variables defined in spreadsheet B_h; one is hini, which is the head value that will be assigned to all active cells during the initialization step, and the other one is Restart, which is used to initialize the heads to the hini value. When Restart is positive, all active cells are set to hini, when Restart is zero, calculations begin. There is a third variable, hundef, with an easily identifiable value to be given to the inactive cells within the rectangular range enclosing the aquifer; in this case, the value chosen is − 99.

The expression for the piezometric at cell I5 is more convoluted than for the confined case because it has to take into account the possibility that a cell becomes dry in one iteration, but it may have to rewet in the next one. With these considerations, the expression for computing the head value at generic cell I5 is

This expression reads as follows:

• If the current cell is inactive in the input data set, set the piezometric head to value hundef; else

  • If the cell corresponds to a prescribed head cell, set the piezometric head equal to the prescribed head value; else

    • If the Restart flag is set to a value larger than zero, initialize the cell with the value hini; else

      • If the cell head is undefined (because it was dry in the previous iteration) and the D_Wet flag is True, meaning that it should be wet (because the surrounding hell heads are above the cell bottom), then

      • Set the value to the rewetting value computed in D_ReWet; else

      • Compute the piezometric head using Eq. (1). For this calculation, function IFERROR is used: an error will occur if the denominator in Eq. (1) is zero, which would happen in the cells with zero saturated thickness. When such an error occurs the head is set to hundef, otherwise is set to the result of the calculation.

2.2 Workflow

Once the input data have been entered, the workflow starts by setting the Restart flag to 1, thus initializing the piezometric heads to the value hini. Then the Restart flag is set to 0, and the iterations will begin. It is convenient to track the convergence of the solution by checking both global and local balances.

The global balance check is set up in the same sheet where the heads are computed. For this purpose, the following variables are calculated: the flux extraction through wells, the recharge through infiltration, the recharge in and out of the river, and the flow entering the aquifer from prescribed head cells. The global error is computed as the sum of the four previous variables, and a relative error is computed as the ratio of the absolute error to the total flux entering from the prescribed head boundary.

The local balance check is performed by computing the fluxes entering each cell and verifying whether they sum up to zero. For this purpose, sheets F_QNorth, F_QSouth, F_QEast, and F_QWest are constructed. For instance, F_QNorth contains the fluxes entering each cell through its northern boundary, and it is computed using the expression

that reads as follows:

• If head is undefined leave QNorth as an empty string; else

  • If the current cell is active after the last iteration and the adjacent cell to the north is not, then set QNorth to zero; else

    • Compute QNorth using Darcy’s law.

Similar expressions are used to compute the rest of the flow components. Finally, sheet G_CellBalance is built with the sum of all the flow components entering a cell. This sheet should contain, after convergence, zero values, except for the prescribed head cells that will contain the flow exchange between the aquifer and the boundary. Iterations should continue (by hitting label “Calculate” at the bottom left corner of the sheet frame) until the error in the B_h sheet is zero (or very close), and all active cells show a zero value (or very close) in G_CellBal.

Figure 9 shows the resulting piezometric heads for the input parameters described before obtained after 8,000 iterations. A heat map is generated filtering out the hundef values (− 99) and spreading from the minimum computed head to the maximum one. This map is achieved using the conditional formatting feature of Excel. The procedure is as follows: select the rectangular area covering the aquifer I5:AO23, go to the Format-Conditional Formatting menu item and set two conditions; the first one is a “classic” style to format cells containing the value − 99 with a customized formatting style which consists in setting the font color to white, and the second one is a “three-color scale" using, for the minimum Formula = MIN(IF($I$5:$AO$23$<>$-99,$I$5:$AO$23)), which computes the minimum value of the head values excluding the cells with a value equal to − 99, for the midpoint the 50% percentile, and for the maximum the highest value (the colors chosen were dark blue, yellow, and red for the minimum, midpoint, and maximum, respectively).

Looking at Fig. 9, we can notice the three elliptical cones of depression formed around the three pumping wells. The elliptical shape is induced by the strong hydraulic conductivity anisotropy in the x direction. The solution has converged to an almost exact solution as indicated by an absolute error of 0 m$^3$ d$^{-1}$ after 8,000 iterations. The evolution of the piezometric heads during the iterations is shown in Fig. 10. After 1,000 iterations, the solution is very far from convergence, with a very large absolute error and a small total recharge but very large discharges through the wells, river, and boundary. After 2,000 iterations, the head pattern is very close to the final solution, in which anisotropy is noticed; the two cells with very high bottom values have already dried, but the absolute error is still high, with a river that seems to be getting water from the aquifer. After 3,000 iterations, the error is low, but it is not until after 8,000 iterations that the solution converges to that shown in Fig. 9 with a final balance of 23,000 m$^3$ d$^{-1}$ extracted by wells, 4,900 m$^3$ d$^{-1}$ recharged by infiltration, 2,510 m$^3$ d$^{-1}$ recharged from the river, 2,140 m$^3$ d$^{-1}$ discharged through the river, and 17,731 m$^3$ d$^{-1}$ entering the aquifer from the lake.

Fig. 9

Piezometric head solution of the groundwater flow problem for the unconfined aquifer

Fig. 10

Evolution of calculated heads as iterations progress

Figure 11 shows the aquifer’s saturated thickness, with the dry cells and the surrounding low-saturated-thickness cells, while Fig. 12 shows the flow entering the aquifer from the river where we can appreciate that the river is a losing river in the northern part, being disconnected for most of the upper stretch, and becomes a gaining river in the lower part. For the stretch along which the river is disconnected from the aquifer (i.e., the piezometric head in the aquifer is below the river bottom), the river recharge into the aquifer is only a function of the difference between the river stage and the river bottom. Since this difference has been set equal to 2 m for all river cells, the river recharge is constant and equal to 100 m$^3$ d$^{-1}$ at the disconnected cells. For the rest of the river cells, the recharge depends on the difference between the piezometric head at the aquifer and the river stage and shows values ranging from 94 to − 163 m$^3$ d$^{-1}$.

Fig. 11

Saturated thickness in the aquifer at the end of the iterations

Fig. 12

Flow entering the aquifer from the river

3 Two-Dimensional Vertical Cross-Section Unconfined Aquifer

With the model for the unconfined aquifer in plan view working, it is easy to modify it so that it can model a vertical cross-section of an unconfined aquifer, aiming at computing the phreatic surface, a boundary condition a priori unknown that renders the solution of the groundwater flow problem as highly nonlinear.

The model setup is for a vertical cross-section of a formation resting on an inclined impermeable bottom with a river crossing, a pumping well, and recharge by infiltration. Figure 13 shows a sketch of the aquifer, with three distinctive hydraulic conductivity zones, a well that is screened only in the lower part, and some relevant elevation data needed for the specification of the boundary conditions; the shaded area corresponds with the saturated zone; recharge coming from the surface will percolate and be added to the cells that intersect the phreatic surface. The numerical model will also contain 19 rows and 33 columns, but cells will be rectangular with a width of 10 m and a height of 5 m.

Fig. 13

Sketch of the vertical cross-section of the unconfined aquifer. The shaded area is saturated

3.1 Implementation

There is a need to introduce a number of modifications in comparison with the previous model.

Step 1. Define active cells and cells with prescribed heads The active cells represent the formation below the terrain surface and up to the impermeable bottom. Figure 14 shows the active cells in the cross-section. The phreatic surface at both vertical boundaries is at 72 m elevation, and the bottom tilts from 15 to 0 m as shown in Fig. 13. It is assumed that there is no vertical piezometric head gradient along the vertical boundaries, and therefore, the boundary conditions are prescribed heads along the first and last column equal to 72 m, the bottom is impermeable, and the a priori unknown phreatic boundary condition closes the boundary at the top.

The cells where the prescribed head intersects the phreatic surface must have a head value that coincides with their elevation, and there cannot be active cells above the phreatic surface. In this case, since the prescribed head is 72 m, the prescribed head cells intersecting the phreatic surface are in row 5, which has a bottom elevation of 70 m and a top elevation of 75 m.

Fig. 14

Sketch of the vertical cross-section of the unconfined aquifer. The shaded area is saturated

Step 2. Input data The cross-section is in the XZ plane. The cross-section is larger along the horizontal direction than along the vertical direction, implying that the discretization cells cannot be square. The cells not being square has two main implications, the first being that both $\varDelta x$ and $\varDelta z$ must be explicitly given (it is still considered that they are the same for all cells), and the second one is that the cell aspect ratio ($\varDelta x/\varDelta z$) plays a role in the calculation of the flows that enter each cell—more precisely, in the calculation of the interblock transmissivities—and consequently, in the expression to compute the piezometric head as a function of the surrounding heads, as will be explained next. (The cell aspect ratio is 1 for square cells, and for this reason, it does not appear in Eq. (1).)

The bottom and top of the cells are fully defined by the bottom elevation of the lower left corner of the rectangle enclosing the aquifer and $\varDelta z$. Therefore, only the lower left corner elevation must be input in spreadsheet A_Bot; the cell bottom and top values are automatically computed ($\varDelta z$ is given in B_h) for all cells.

Also, for the proper modeling of the recharge by infiltration and the river–aquifer flow exchange, it is necessary to specify which is the width of the aquifer perpendicular to the cross-section shown in the figures, that is, the value of $\varDelta y$. Neither recharge nor river inflow/outflow will be the same if the aquifer section is 1 m wide as if it is 100 m.

For these reasons, in spreadsheet B_h, values for deltax, deltay, and deltaz are provided, and in spreadsheet A_R, containing the river conductance, the values for the riverbed thickness and hydraulic conductivity are given. In this latter sheet, conductance is computed considering that the river covers the entire cell as $R = K' b^{-1}\varDelta x \varDelta y$, where $K'$ is the riverbed conductivity, and b is the riverbed thickness. There is a caveat in the calculation of the river infiltration: when the head is below the cell bottom containing the river, no river infiltration will be considered. It is possible to build a mechanism that would compute the infiltration from a disconnected river and carry it down to the phreatic surface, but it would complexify too much the workbook.

The recharge is computed in intermediary spreadsheet D_QN, where the active cells after the last iteration that intersect the phreatic surface are identified and assigned a flow entering by recharge equal to $Q_N = N \varDelta x \varDelta y$, where N is the infiltration rate. In practical terms, this is achieved using the following expression

that reads as follows:

• If the cell is inactive set the cell value to an empty string; else

  • If the cell has a prescribed head set the recharge to zero; else

    • If the cell above is inactive then assign a recharge equal to the infiltration recharge (named variable N1_) times the cell size in the x direction (named variable deltax) times the cell size in the y direction (named variable deltay); else set it to zero.

Step 3. Compute intermediate variables The saturated thickness intermediate variable refers to the saturated thickness on the plane of the cross-section, as opposed to the saturated thickness orthogonal to the plan view used in the previous case to compute the cell transmissivity. It will be equal to $\varDelta z$ for all saturated cells and equal to the difference between the calculated piezometric head and the cell bottom elevation for the cells intersecting the phreatic surface.

There is a need to account for the aspect ratio $\varDelta x/\varDelta z$ in the computation of the flows entering each cell and modify the calculation of the piezometric head accordingly. Equation (1) could still be used for the calculation of the piezometric head if the cell aspect ratio is incorporated to the calculation of the terms $T_N', T_S', T_W', and T_E'$ to transform them from interface transmissivities into what MODFLOW calls conductances and petroleum engineers transmissibility. The conductance is the value that multiplied by the head difference between two adjacent cells provides the flow crossing their interface; this value is equal to the interface transmissivity multiplied by the interface area and divided by the distance between cell centers. For this purpose, spreadsheets D_Tx, D_Tz, E_TN, E_TS, E_TW, and E_TE are modified as follows. D_Tx and D_Tz contain the transmissivity defined as the conductivity multiplied by the saturated thickness measured orthogonal to the cross-section; since the saturated thickness orthogonal to the cross-section is $\varDelta y$, D_Tx and D_Tz are equal to the values in A_Kx and A_Kz multiplied by $\varDelta y$. Then to obtain the conductances, first, the interface transmissivity is computed as the harmonic average of adjacent transmissivities (as before); then this value is multiplied by the cell saturated width (in the plane of the cross-section and orthogonal to the flow direction) and divided by the distance between cell centers (in the flow direction). This calculation will be different when computing conductances along the vertical axis, in which case, the cell width is constant and equal to $\varDelta x$, and the distance between cell centers is constant and equal to $\varDelta z$, from when computing conductances along the horizontal axis when the cell saturated width is variable and equal to the saturated thickness, and the distance between cell centers is constant and equal to $\varDelta x$. Figure 15 shows the saturated thickness at the end of the iterations.

Step 4. Iterative solution for piezometric heads Nothing has changed in the iterative calculation of the piezometric head. It must be noticed that when the aspect ratio $\varDelta x/\varDelta z$ is very large, the convergence slows down. Figure 16 shows the hydraulic head distribution for an aquifer’s vertical cross-section discretized in cells 10 m wide by 5 m tall, where there is a well pumping 35,000 m$^3$ d$^{-1}$. The resulting distribution shows a clear cone of depression around the cells where the well is open. Water flows from the prescribed head boundaries toward the river and toward the well. The total infiltration is very small since the exposed surface through which recharge would enter the aquifer is only $31\times 10 \times 10 \times 0.001$ (number of cells crossing the phreatic surface times $\varDelta x$ times $\varDelta y$ times the infiltration rate). The river drains 1,990 m$^3$ d$^{-1}$ from the aquifer.

3.2 Workflow

The workflow is exactly the same as before: enter the data in the corresponding spreadsheets, initialize the heads to start the iterations, and iterate until convergence.

Fig. 15

Cell saturated thickness (on the plane of the cross-section)

Fig. 16

Piezometric head solution for the vertical cross-section in an unconfined aquifer

3.3 Extension to the Simulation of Flow Through an Earth Dam

Seepage through an earth dam can be modeled as groundwater flow in an unconfined aquifer. Consider the dam sketched in Fig. 17 made up of a core with a conductivity of 0.01 m d$^{-1}$ and an embankment with a conductivity of 0.1 m d$^{-1}$. A vertical cross-section unconfined aquifer model is built using a discretization of 51 rows and 113 columns; each cell is 5 m wide and 2 m tall. The water level on the upstream side of the dam is at 80 m and on the downstream side at 16 m. The bottom is horizontal and impermeable at level 0 m. With these premises, the different spreadsheets are built. The final solution is shown in Fig. 18; where high gradients occur within the core, there is virtually no gradient on the upstream side of the embankment and some gradient on the downstream side. Considering that the dam length is 1,000 m orthogonal to the section displayed ($\varDelta y= 1,000$), the total seepage through the dam results in 320 m$^3$ d$^{-1}$.

Fig. 17

Sketch of an earth dam section with core and embankment

Fig. 18

Hydraulic head distribution within the earth dam

4 Conclusions

This paper extends the work presented by Gómez-Hernández (2022) to the solution of groundwater flow for unconfined aquifers both in the horizontal plane and in a vertical cross-section. The unconfined groundwater flow equation is highly nonlinear since the transmissivities depend on the piezometric heads, and also because cells may dry out, which requires the introduction of a mechanism that allows a cell to become dry (if the head goes below its bottom elevation) and to rewet (if, later during the iterations, the heads in the nearby cells suggest that the cell should be wet). For the case of the simulation of a cross-section, the phreatic surface is a boundary condition that is not known a priori, and that depends on the solution of the equation, which also introduces complications that have to be tackled during the implementation as explained.

Overall, the spreadsheet approach yields very good results that have been contrasted with the solution of the groundwater flow equation with MODFLOW (Harbaugh 2005), which are not shown here but that the reader can verify by making use of the Python Flopy scripts provided as companions.

Further work should focus on transient modeling. While the calculation of the piezometric time evolution would not be complicated, the management of the results aimed at producing time-dependent graphs is not trivial. Also, the possibility of having nonuniform discretizations with local grid refinements could be addressed in further works.

Finally, there is the possibility of building a three-dimensional model by having multiple horizontal layers placed in each spreadsheet. Again, this implementation should not be too difficult, although the handling of such a large spreadsheet may lose the attractiveness of the tool as a teaching/quick demonstration tool.