Tutorial 2

Building a multi-layer heterogeneous, anisotropic model with various boundary conditions

Introduction

This second tutorial builds on the model created in the first tutorial.  The first tutorial model was two-dimensional and steady, but this second tutorial will guide you through creating a model that has a multi-layer three-dimensional (3D) region in an area near the middle reach of the river.  The well will be changed from fully-penetrating to partially-penetrating.  This model may be built with the Demo version of Anaqsim, and will contain two levels in the multi-level region.  More levels are possible with the Deep Designer licensed version of Anaqsim and there is some instruction about that at the end of the tutorial.

Opening the First Tutorial Model

Either begin with the model you had after competing the first tutorial, or download and then open the following file tutor1.anaq from Anaqsim.  To open the file in Anaqsim, select File / Open from the main menu, then locate and open the tutor1.anaq file.  Next select File / Save As and save this under another name such as 'tutor2.anaq'.

Zoom in on the area where you will make the model multi-level (2D) so that the resulting view looks something like the following image.  In this case, the basemap (available in the C:\Program Files\Yellow Sub Hydro\Anaqsim\Documentation folder) is displayed in gray.  The model elements are displayed in green.

Adding New Domains for the Multi-Layer Area

In the area near the well and the middle river reach (h=144 to h=108), we will create a region with two levels instead of one.  The first step is to create new domains for this region.  Select Model Input / Domains / Confined and/or Unconfined to bring up the domains data table, which contains two domains 'all' for most of the model, and 'recharge basin' for the recharge basin area.  These two domains have the same properties, but differing recharge rates.

We will add 2 new domains for the multi-level region.  The quickest way to create the new rows is to duplicate the two existing rows.  Select both existing rows as is shown below.

Then right-click and choose Copy Selected Rows and then right-click and select Paste New Rows. Edit the first 5 columns of the new rows so that they look like the following table.  Columns 6+ will all be unchanged, since these domains all have the same properties, just different elevations.

Note that the bottom elevation of the new two-level area is 80, a bit lower than the bottom elevation of the 'all' domain at 82.  This is just to show that elevations may change laterally in a model if you want them to. The '3D upper' domain is in level 1 and is unconfined with a base at elevation 95.  The '3D lower' domain is in the 2nd level (level numbers increase with depth) and is confined with a constant saturated thickness of 15, going from elevation 80 to 95.

At this point, save your model with File / Save from the main menu or Ctrl-S.

Interdomain Line Boundary to Define the 3D Region

We have defined the two new domains that will apply in the 3D region, but still need to define the area where those domains exist.  Anaqsim does this with external line boundaries that refer to the domains.  In this case, we will create one new interdomain line boundary that will represent the lateral transition from the single-level 'all' domain to the multi-level 3D area.  At interdomain boundaries, the conditions imposed are continuity of flow and continuity of head across the boundary.

First digitize the interdomain boundary.  Select the Plot tab, and make sure that under the plot menu Snap Settings / Snap to Elements is checked.  This allows the digitized interdomain line boundary to 'snap to' the same vertex coordinates as the river line boundaries.  When the cursor moves over an existing vertex, an orange box will appear; if you click with the orange box showing, it will snap to the same coordinates. Digitize a closed polygon boundary  like the yellow one shown below, making sure to snap to the river vertexes labeled 'h=114' and 'h=108'. You can go around the boundary in either clockwise or counter-clockwise order, but it is important to remember which way you went.

Next, select Model Input / Line Boundaries / Inter-Domain and add a new row so it looks like this:

The 'ID' boundary defined the limits of the recharge basin, and the new '3D limits' one will define the boundary of the 3D area.  Click on the Edit button in the '3D limits' row and paste in the coordinates you just digitized, and then press OK.

Then, click on the Select button under Domains_left. If you digitized in clockwise order, select the 'all' domain and if you digitized in counter-clockwise order, select both '3D Upper' and '3D Lower'. Imagine walking on the ground surface along the line you just digitized; Domains_left would be to your left as you walk along the boundary and Domains_right would be to your right. This interdomain boundary makes a transition from two levels (3D upper, 3D lower) on one side to a single level (all) on the other. By doing this, interdomain boundaries can help you efficiently concentrate computing power and multiple levels in a limited area. Interdomain boundaries can also join a single domain on one side with a single domain on the other side.

Now click on the Select button under Domains_right. If you digitized in clockwise order, select both '3D Upper' and '3D Lower' and if you digitized in counter-clockwise order, select the 'all' domain.

To check your input, select Make Plot / Model Elements Only, then select level 1 and press OK, and you should see something like this in the area near the well.

The new interdomain line boundary is now green, since it is composed of elements in the model. Move the cursor over the vertexes of the interdomain line boundary to check that the left/right sense of the domains is correct. Move the mouse across the boundary, and observe the domain listed at the left. When you are outside the 3D area, you should see this:

When you are inside the 3D limits you should see this:

Now make a similar plot, but chose level 2 instead of level 1 and you will see something like this:

The river elements do not show up, since they are in level 1.  Move the mouse across the boundary, and observed the domain listed at the left. When you are outside the 3D area, there is no level 2 and you should see a blank under Domain Name:

When you are inside the 3D limits you should see this:

Putting the Well and River Reach into the New Domains

The two domains that occupy the 3D area are '3D upper' and '3D lower', but our input still lists the middle river reach and the well as being in domain 'all', which is no longer correct. We need to change the input for these elements and put them in the correct domains.

First, bring up the well data by moving the cursor near the well, right-clicking, and then selectin Edit Nearest Well. Then change the domain of the well from 'all' to '3D lower':

Now the well is screened only in the '3D lower' domain from elevation 80 to 95, not in the '3D upper' domain above. In this manner, Anaqsim simulates partially-penetrating wells.

Second, edit the middle river reach by moving the cursor near this line boundary, right-clicking, and then selecting Edit Nearest Line Boundary. This will open the river line boundary table and highlight the middle reach row. Change the domain of the well from 'all' to '3D lower':

The river now will be in the upper layer and there will be vertical resistance and 3D flow with the underlaying 2nd level.

Save your model now, selecting File / Save from the main menu or Ctrl-S.

Spatially-Variable Area Sinks in 3D Region

In the 3D area there will be vertical leakage between the two levels.  For each domain, this can be thought of as a distributed source/sink term like recharge, but it will vary with location.  For example, near the well that is in '3D lower' there will be vertical leakage from '3D upper' to '3D lower' with leakage increasing closer to the well.  To model this spatially-variable leakage, we need to use spatially-variable area sinks (SVAS). SVAS are required to simulate vertical leakage in multi-level areas of models, and also to simulate storage fluxes in transient models. In this steady model, SVAS are needed just in the 3D area where there is vertical leakage.  Outside the 3D area, the model is single layer, steady, and the uniform recharge can still be modeled efficiently with a uniform area sink so there is no need to modify that input.

SVAS by Domain

SVAS can be input over the area spanned by a domain or the area spanned by  polygon. For this model, we will use both. First, we will enter an SVAS that apples over the 3D area. Select Model Input / Area Source/Sink / Spatially-Variable, Domain, and complete a row like this:

This SVAS applies o the area spanned by the '3D upper' domain. SVAS are applied for all levels of a model, including those that underlie and/or overly the specified domain. Since '3D upper' and '3D lower' span the same area, we only need to create an SVAS for one of these domains. The result would be the same if we chose either '3D upper' or '3D lower' as the Domain. Condition_Top is set to Flux, meaning that the sink/source rate into the top of the the model is specified and equal to 0.0006, the same recharge rate that applies in the 'all' domain outside the 3D area.  Condition_Bottom is set to Flux = 0, so the bottom of the model (bottom of '3D lower') is impermeable.  The alternative top or bottom condition is head_dependent_flux, which is used for vertical leakage to a specified head at a surface water, for example.

Node_spacing is the approximate spacing between basis points in the SVAS. The leakage/storage fluxes are solved for perfectly at the basis points, and approximated with a smooth interpolated surface between basis points. Anaqsim will automatically populate the area of domain '3D upper' with basis points on a regular hexagonal grid of points, spaced about 80 units apart. To read more about SVAS, see papers by Fitts (2010) and Strack and Jankovic (1999).

To see how these basis points are distributed, select Make Plot / Model Elements Only, then select level 1 and press OK. You should see green '+' basis point symbols in the 3D area something like this:

A similar plot of level 2 would show basis points at these same locations. Move the cursor over some basis points to see the domain and top / bottom conditions at this point.

SVAS by Polygon

Where the rates of vertical leakage change over short horizontal distances, we need denser basis point spacing. We can use polygon SVAS and special well basis point spacing for this purpose. First, we will create a zone of denser basis points in the well / river vicinity using a polygon SVAS. Start by digitizing a polygon (either clockwise or counter-clockwise) that looks something like the yellow line in the following image. If you don't close the polygon, Anaqsim will close it for you from the first point to the last.

Then select Model Input / Area Source/Sink / Spatially-Variable, Polygon and edit the row so it looks like this:

Then click on the Edit button and paste in the coordinates of the polygon you just digitized. This will spread out basis points with a small spacing (25) inside the polygon. The domain basis points that were within this polygon will be automatically overwritten. The Nesting_level is used when you have multiple polygon SVAS that are nested within one another.  Polygon SVAS with higher nesting levels overwrite polygon SVAS with lower nesting levels.  We will input just one polygon SVAS, so nesting won't be an issue.

Save your input, and then select Plot Input / What to Plot and then check the SVAS_polygons box so that the SVAS polygons will be displayed in subsequent plots:

Select Make Plot / Model Elements Only and choose level 1 to see the distribution of basis points and the SVAS polygon (light blue):

Well SVAS Basis Point Spacing

Vertical leakage rates change especially rapidly near a partially-penetrating well, so Anaqsim has a built-in means of creating basis point spacing that gets denser and denser close to a well. See the User Guide to read more about how this works.  Select Model Input / Area Source/Sink / Spatially-Variable, Well Basis Points and create a row like this:

Then click the Select button and select the one well in the model:

This will create concentric circles of basis points, 6 per circle, with smaller and smaller circles until the innermost one which is just outside the well radius. The largest circle will be drawn such that the spacing between points is no more than Max_node_spacing. Generally set Max_node_spacing to a value close to the ambient domain or polygon basis point spacing.  Switch to the Plot tab, zoom in on the 3D area, and then right-click and select Set Plot Window to Current View, which will make subsequent plots in this smaller area. Now select Make Plot / Model Elements Only and choose level 1 to create a plot like this:

Zoom in on the well to see the special well basis point placement more clearly:

At this point, save your model.

Solving the System of Equations

Most of the elements in the model have unknown strengths that will be solved for by building a system of boundary condition equations based on conditions at line boundary control points, SVAS basis points, etc. All this happens when you select Solve on the main menu. This will initiate the solve process, switch the view to the Log tab and print out information as the system is solved. For my version of this model, here is what was written to the log:

The system has 439 equations and 439 unknowns. It took 5 iterations and about 1 second to solve within the convergence criteria listed under Model Input / Solution / Check Settings. Some of the equations in the model were nonlinear (vertical leakage equations at basis points, river equations) due to the nonlinear nature of the head-potential relationship in unconfined domains. Solving nonlinear equations requires iteration using updated linear approximations at each iteration.

What to do if your model has too many equations

The Demo and Quick Designer Licensed versions of Anaqsim are limited to no more than 2000 equations, and depending how you digitized your SVAS polygon, you may have exceeded this limit.  If this was the case, Anaqsim gave you an error message telling you to simplify your model and reduce the number of equations.  Here are some ways that you can simplify and reduce the number of equations:

• Graphically edit the SVAS polygon boudnary, making it cover a smaller area.
• Increase the basis point spacing for the domain SVAS and/or the polygon SVAS.
• Decrease the parameters per line for the line boundaries, particularly the interdomain boundaries, which generate more equations per segment than other types of line boundaries.

Use one or more of the above strategies, and then press Solve. Repeat until the model is within the 2000 equation limit and solves successfully.

Making a Plot of the Solution

Once the model is solved, Anaqsim offers many ways of examining the solution. The first to explore is a plan-view plot that shows head contours and velocity vectors. Before making the plot, we will adjust some of the plot input.

Plot Input Settings

Select Plot Input / What to Plot and check the boxes so they are like this:

Leaving the Window cell blank makes the plot cover the entire extent of the model. We will make a plot that shows the checked items.

Next select Plot Input / Contour Settings and adjust it like this:

The blank cells for Minimum and Maximum mean that the minimum and maximum contours will be determined such that they cover the full range of head (h) in the plot area. Select Plot Input / Vector Settings and adjust to look like this:

This will draw about 600 vectors proportional to average linear velocity. Read more about these and other plot input settings in the User Guide.

Make a Plot of Entire Model Area with Head Contours and Vectors

Make the plot by selecting Make Plot / All Selected Features and then select level 1 and press OK. Your plot will look like this:

Make Close-up Lots of 3D Area, Both Levels

To see more clearly the area of interest, zoom in so the 3D area covers most of your plot window, something like this:

This is just the same plot, but enlarged. You can easily make a more accurate plot of this close-up window. First, change the plot window to this view by right-clicking over the plot and then selecting Set Plot Window to Current View. Also, under Plot Input / Contour Settings, change the Increment to 0.2 instead of 0.5:

Now make a new plot by selecting Make Plot / All Selected Features selecting level 1 and pressing OK. The plot should look like this:

If you repeat this, but make the plot of level 2, you get something like this:

The results show different head patterns, with more pronounced connection the river in level 1 and deeper drawdown at the well in level 2 where it is screened. The anisotropy of the domains makes the drawdown ellipse-shaped, with the K₁ axis orientated 30 degrees north of east. The plot of the 2nd level is limited to the 3D area - there is no 2nd level outside the interdomain boundary that defines the limits of the 3D area. There is some drawdown near the well in level 1 due to strong downward vertical leakage between levels.

Examining Head Differences Between Levels

You can examine the difference in head between the two levels in a couple of ways. One is with the cursor information at the left. If you make a plot of heads in level 1, then move the cursor around the 2D area, you will see fairly large head differences listed at the left when you are near the well:

The positive number means the head in level 1 is greater than the head in level 2 and there is downward flow. If you move the cursor near some parts of the river, you may see a negative head difference indicating upward flow:

You can also make a contour plot of the head difference between levels. Select Plot Input / Contour Settings, and for Surface, select dh, which makes the contour plot a plot of the difference in head between this level and the one below:

If you select Make Plot / All Selected Features and choose level 1, the resulting plot of head differences looks like this:

Over most of the 3D area, the difference in head is less than 0.2, but there are larger positive differences near the well, and larger negative differences near the northern part of the middle river reach. The heavier contour is the zero contour. At the interdomain boundary that is the limit of the 2-level area, heads are forced to be equal by the boundary conditions, which is met at control points along the boundary and hence the zero contour weaves close to the boundary.

Making a Plot with Pathlines

Now we will make a plot with pathlines emanating from the well. Select Plot Input / What to Plot and edit the table so it looks like this:

The blank under Window makes subsequent plots cover the entire model, and we will display pathlines, but not vectors or SVAS polygons.

Select Plot Input / Contour Settings and edit the settings to show head contours with an increment of 0.5:

Select Plot Input / Pathline Settings and check the settings are as follows (similar to what we had in the first tutorial):

Select Plot Input / Pathlines / Well and edit the row so it is like the following, which will start 15 pathlines 50% of the way down in level 2 just outside the well radius, and trace pathlines upstream.

Now, select Plot Input / Pathlines / Line and if there is a row of data there from the first tutorial, uncheck the Show column since we only want to see well pathlines.

Now make a plot of level 1, which will result in something like this:

If you move the cursor over the pathlines near the well, you will see that they are in level 2, and as you move the cursor upstream, the pathlines generally rise because of recharge imposed on top of the model; some pathlines that enter the well from the northeast rise into level 1 within the 3D area. The pathlines east of the well migrates steeply downward in a stagnant zone east of the well (time markers are closely spaced here).

Analyzing the Solution in Detail

Anaqsim offers many ways to check and analyze the solution in detail. In this section we will explore some of them that apply to this model.

Examining Head at the Discharge-Specified Well

There are a couple of ways to check the head at the well. First, you could select Analysis / Discharge-Specified Well Heads / Write Heads to Run Log, which would display this in the Log tab:

Another way is to move the cursor near the well in Plot tab, right-click and select Check Nearest Well Head and Discharge, which will display a window like this:

Checking River Reach Discharge

Selecting Analysis / River Line Discharges / Write Discharges to Run Log writes the total discharge of all river reaches to the Log tab like this:

All three river reaches have negative discharges, meaning they are extracting water from the aquifer. You can see how discharge is distributed along a particular segment of a river reach by moving the cursor near a segment, right-clicking, and selecting Check Nearest Line Boundary Conditions. The following is the resulting plot for the segment just west of the well:

This segment loses water at one end and gains water at the other end.

Making Profiles along a Line

A useful tool in Anaqsim allows you to create a variety of profiles along a line.

Profiles Showing Head, Domain Boundaries, and Pathlines

Start by digitizing a line along which you will make profiles, similar to the one shown in yellow below, which goes approximately through the well.

Then select Analysis / Graph Conditions Along a Line and paste the digitized coordinates into the Coordinates box and edit the rest of the form as follows:

When you press the Make Plot button, you will see a profile of heads like this:

The domain boundaries are shown in gray, and you can clearly see the 2 layer 3D region roughly where 1200 < x < 1700. The red dot time markers are placed at the locations of each time marker arrow shown in the plan-view plot. The pathlines in this plot do several interesting things:

• Some are pulled up abruptly where they pass under a discharging river line (e.g. the lowest pathline at about x = 850)
• Some dive steeply as they pass under the recharge basin ( several at about x=1050)
• One dives steeply east of the well in a stagnant zone (about x=1650)
• Some jump vertically as they cross the interdomain boundary from the single-level 'all' domain to the 2-level 3D region (e.g. at about x = 1200). The jump occurs because pathline elevations are determined based on water balance. For example, if 40% of the discharge is above the pathline on one side of the interdomain boundary, 40% of the discharge is above the pathline on the other side of the boundary. Because the base elevation changes and because the flow field is different in the two layers, pathlines jump appropriately at these boundaries.
• Some pathlines appear to terminate above or below the water table (Head, level 1). This is because pathlines are projected onto the line, and most of these terminations are considerably to the side of the profile line and the water table elevation differs from the water table elevation on the profile line.

If you want to change the axes (for example, make the x axis start at a different value), you can do that by selecting Change Axis in the menu at the top of the plot. Read in the User Guide about saving and copying chart data and images.

Profiles Showing Domain Discharges

If you select Analysis / Graph Conditions Along a Line and edit the dialog so you check Domain discharge tangent to line and make the plot you will get a graph like this:

This shows the domain discharge, which is specific discharge times saturated thickness, tangent to the line. The discharge increases generally from x=0 to the well, as added recharge and impact of the well increase the discharge going east along the line. In level 2 where the well is screened, there is an abrupt shift in tangent discharge from positive to negative as you cross through the well. There is a similar, but muted shift in level 1 above the well. The steps you see in the level 1 profile are the result of river discharges subtracting (about x=850) and adding (about x=1400) to the discharge. At about x=1200 and x=1700, the profile crosses the interdomain boundary at the edge of the 3D area. Here the sum of the tangent discharges in levels 1 and 2 in the 3D area roughly equal the total discharge n the matching single-level area (the normal component of discharge is matched at interdomain boundaries, components in other direction can differ)

Profiles Showing Extraction, Evaluation SVAS Accuracy

If you select Analysis / Graph Conditions Along a Line and edit the dialog so you check Extraction and make a plot, you will see something like this:

Extraction is the net distributed flow out of the domain due to vertical leakage and storage fluxes, in units of discharge/area [L/T]. In the single-layer are of domain 'all', the only distributed flux is recharge coming into the domain a a uniform rate of 0.0006; since this is into the domain, it is a negative extraction (-0.0006 in the interval 0<x<1200). In the multi-level 3D area, in addition to recharge in level1, there is vertical leakage, which is mostly downward out of level1 and into level 2 where the well is screened. The vertical leakage is much larger magnitude than the recharge, so there is a large positive extraction in level1.  There is the opposite in level2, as the leakage from level 1 creates a significant distributed source.

That the modeled extraction closely  matches the extraction from heads indicates that the SVAS are working well. At this point, save your model, which is the final model for the second tutorial.

To see how basis point spacing can affect the accuracy of SVAS, we will tweak the model to make it have inadequate basis point spacing. Select Model Input / Area Source/Sink / Spatially-Variable, Well Basis Points, select he one row there and then delete it:

Then delete the SVAS polygon:

Now the model has sparse basis points in the 3D area, based on the domain SVAS settings.

Select Solve to solve this new, less accurate model. Then select Analysis / Graph Conditions Along a Line and make an extraction plot, which will look like this:

You can see with this profile that the comparison between modelled and heads-based extraction is poor in the 3D area. Use this tool to check if your basis point spacing is fine enough to give a good solution. See the user Guide for more details about extraction and SVAS. Don't save the model at this point, as it was only to point out how the extraction check can be used to evaluate SVAS accuracy.

Further Exploration for Anaqsim Deep Designer Users

This tutorial has been kept relatively simple so it could be done with the Demo and Quick Designer License. For Deep Designer licensed users, it is possible to expand and explore a bit more complexity. Here are some suggestions for modifying this model and adding some complexity:

Make the 3D area 4 levels

To do this, you will need to follow these steps:

• Edit the Domains table so you have two new domains that splice in between the top and bottom domains. Add two new domains (e.g. '3D 2' and '3D 3') and adjust the elevations and levels of these and the other domains in the 3D area. Don't change the names of the exiting domains because that will affect other elements that refer to domain names.
• Edit the '3D limits' interdomain line boundary so it refers to the new suite of domains on the multi-level side of this boundary.
• Edit the well and place it in the domain you want it in, also edit the well pathlines so they start at the desired level.
• Solve the system, make plots and examine the solution.

Add Head-Specified and/or Discharge-Specified Multi-Domain Wells

Anaqsim allows two other types of wells that you could experiment with. Read about and add head-specified and/or discharge specified multi-domain wells. The multi-domain well would have to be in the 3D area of the model, and could span 2 or more levels there. For multi-domain wells, the same head at the well radius is maintained in all domains spanned and the total discharge of the well elements in the spanned domains equals the specified discharge. Explore the tools in the Analysis menu for examining the discharges and heads at these types of wells.

Wrap up and Next Steps

This concludes the second tutorial. If you want to see the input file that was generated by this tutorial, then it can be found here: tutor2.anaq.

To explore transient and multi-level modelling with Anaqsim, please see the third tutorial, which edits this model to create transient simulations where the well discharge and recharge rates change with time.