# User Guide

The in-depth guide on how to use Anaqsim

- Introduction
- Anaqsim Modeling Concepts
- User Interface
- General Modelling Sequence
- Menus
- File Menu
- Edit Menu
- Model Input Menu
- Plot Input Menu
- Analysis Input Menu
- Solve and Cancel Solve
- Make Plot Menu
- Analysis Menu
- License Menu
- Help
- About Anaqsim
- Plot View Options
- Computational details
- Checks of Anaqsim
- References

# Introduction

Anaqsim is an analytic element software for simulating groundwater flow. It uses subdomains as described in Fitts (2010), which gives it strong capabilities with respect to heterogeneity and anisotropy. It also employs high-order line elements, spatially-variable area sinks, and finite-difference time steps to allow multi-level aquifer systems and wide-ranging transient flow simulations. Anaqsim is a product of Practical Groundwater and was developed by Dr Charlie Fitts. It is coded in C#.

### Anaqsim Versions

There are two versions of Anaqsim: 'Student' and 'Licensed'. The installed software is the same for both, but the capabilities of the software vary depending on whether the software is Licensed or not.

By default the download provides you with an un-licensed version that will not solve any models.

The student version is not for commercial use, and incorporates a watermark. It can be used to run all of the tutorial models on the website.

The student version allows modeling of a wide array of situations, both steady and transient. When solving, its limits are as follows:

- Allows multi-level (3D) simulations with up to three aquifer levels.
- It allows steady state and limited transient simulations with up to two different time periods and up to five time steps in each period.
- The number of equations in the model's system of equations is limited to 2000. This allows fairly complex models, but users needing greater complexity, or for commercial use should purchase a license.

Activating a commercial license (either the free trial, monthly or annual subscriptions) allows you to solve larger, more complex models with these capabilities:

- Allows multi-level (3D) simulations with up to 15 aquifer levels.
- It allows steady state and fully transient simulations with any number of time periods and up to 20 time steps in each period.
- The number of equations in the model's system of equations is not limited, but somewhere around 40,000 to 50,000 equations is a practical upper limit based on computer memory and the time required to solve.

Anaqsim is licensed through either a monthly or annual recurring subscription that can be cancelled at any time. Please refer to licensing section for further details, and see www.anaqsim.com for license and pricing details.

It is possible to plot and analyse results from a model that has been solved on a licensed version of Anaqsim on an unlicensed version. See installation-manual for further details.

### System Requirements

- 64-bit Windows operating system.
- At least 500 MB of available disk space.
- At least 1 GB memory. More memory is needed for larger problems and most users have 8+ Gb.
- Microsoft .NET v 4.8 or higher.

### Installing Anaqsim

Please see the Installation documentation for details on how to install Anaqsim, and how to activate and de-activate licenses. Anaqsim licenses can be purchased as either a monthly or annually recurring subscription via the website.

### Release Version and History

The current Anaqsim release number and a history of Anaqsim releases, including lists of changes from one release to another, are posted on the Anaqsim website. The release numbering scheme gives the year followed by the release number in that year. New releases are posted to the website and may be downloaded, installed, and run by anyone with a valid license. See the next topic about how to update to a newer release.

You can see the release number of your installed Anaqsim by selecting About Anaqsim on the main menu.

### Updating to a Newer Release

As long as you are within the term of your Anaqsim license, you may upgrade Anaqsim to the current release. To update, first check that there is a newer release available from www.anaqsim.com/version-history. If there is and you want to update:

- Download the newer release installation file from the www.anaqsim.com/version-history to your computer,
- Uninstall the older release from your computer using the Windows Control Panel -> Add/Remove Programs dialog,
- Install the newer release as described in the Installing Anaqsim section.

Updating will not affect your license, as the Anaqsim license file is not removed by the uninstall operation.

### Help and Documentation

Help is available in many forms: this User Guide, tutorials, and example models available at the website.

The Anaqsim User Guide can be accessed by selecting Help from the main menu in Anaqsim.

Three tutorials that walk the user through construction of Anaqsim models of increasing complexity are available on www.anaqsim.com. A quick way to get familiar with Anaqsim is to work through these tutorials in sequence, building your own Anaqsim models like the ones shown in the tutorials. The tutorials lead you through creating these models step-by-step. The tutorials are pdf files with outlines at the first page that contain clickable links, so you can easily jump forward and backward to review, if needed. The completed tutorial models are included in the Documentation directory of the Anaqsim software installation:

- tutor1.anaq is the input file for a simple one-level steady model with irregular boundary conditions, a recharge basin, a river, and a well, and anisotropic K.
- tutor2.anaq is the input file for a steady model with heterogeneity and a 3D area with multiple levels.
- tutor3.anaq is the input file for a transient dewatering simulation model within a 3D area.

All three tutorials may be done with either the student or licensed Anaqsim.

### Contact and Support

Support is available for licensed Anaqsim users to make sure that Anaqsim is functioning properly on their computer. Please first check the Anaqsim User Guide (Help) and the tutorials on the website to see if your question may be answered there.

The contact information for support is support @anaqsim.com

# Anaqsim Modeling Concepts

This section does not repeat details that are published elsewhere in books and articles, but instead gives a quick outline of the techniques used in Anaqsim and cites the appropriate references for those interested in the details.

Anaqsim employs the analytic element method (AEM), which superposes analytic solutions to yield a composite solution consisting of equations for head and discharge as functions of location and time. The AEM is described in detail in books by Strack (1989) and Haitjema (1995). Shorter summaries of the method may be found in Fitts (2023) and Strack (2003). The AEM is fundamentally different than numerical methods like finite elements and finite differences, where the domain is broken into small blocks or elements and simple head distributions (e.g. linear) are assumed within these blocks or elements. In the AEM, boundaries of the domain are discretized, but the domain itself is not.

Anaqsim uses a variation of the AEM that divides the modeled region into subdomains, each with its own definition of aquifer parameters and its own separate AEM model (Fitts, 2010). The model for a given subdomain (called a domain in Anaqsim) includes contributions from elements inside and on the external boundary of the subdomain; elements beyond the subdomain do not contribute. Each subdomain model is written in terms of two-dimensional functions, but three-dimensional flow may be simulated using multiple levels in a model. In multi-level models, the resistance to vertical flow is accounted for in the vertical leakage between levels.

This subdomain approach allows for a high degree of flexibility with respect to a model's heterogeneity, anisotropy, and layering. For example, it is possible for a subdomain that is anisotropic to be adjacent to another subdomain that is isotropic or anisotropic with a different direction and ratio. The subdomain approach allows mixed layering schemes. For example, an area with multiple levels (subdomains stacked vertically and leaking vertically to each other) can abut an area with subdomains in just a single level. This allows the model to focus layering and computational effort in the area of interest, with a simpler single-level model for distant areas.

Another key aspect of Anaqsim is that it allows complete transient simulation capabilities by using finite difference time steps as suggested by Haitjema and Strack (1985). The transient term in the flow equations is handled in essentially the same manner as it is in finite difference programs like MODFLOW.

### Subdomains and Model Levels

Anaqsim uses one separate two-dimensional model for each subdomain. In Anaqsim, subdomain input is entered in the Model Input / Domains data tables (Domain in Anaqsim's menu system is short for subdomain). In each subdomain model, the resistance to vertical flow is neglected and the head is independent of elevation within the subdomain (Dupuit assumption). Resistance to vertical flow and three-dimensional flow are modeled by using multiple levels with vertical leakage between levels. To illustrate how subdomains and model levels are implemented in Anaqsim, we will go through several examples, starting from simple and working toward complex.

The simplest model would be one level (two-dimensional), and only one subdomain (homogeneous). A plan view of such a simple model is shown below.

The properties of the subdomain (hydraulic conductivity, porosity, elevations, etc) are defined in the Model Input / Domains menu. The spatial extent of the subdomain (blue) is defined by the line boundaries that are listed as being on the external perimeter of the subdomain. This is different from TWODAN or other analytic element programs where one domain is the infinite "background" domain, and a heterogeneity (domain) lies inside a polygon dedicated to defining the heterogeneity boundary. The scheme in Anaqsim allows different kinds of boundaries to define the limits of a domain, which is much more flexible. In the simple model shown above, there are two line boundaries that form the external boundary of the domain: one that is head-specified (red) and one that is normal-flux specified (black).

The algorithm that determines what points lie inside the domain works as follows: a point is inside the subdomain if going in the positive x direction (right) from the point, you cross a boundary of the domain an odd number of times. For most points inside the domain, there is only one boundary crossing in the positive x direction (like "a" above), but it is possible to have three or more crossings (like "b" above).

In general, subdomain boundaries should combine to form a perfectly closed polygon, with exactly matching points where the line boundary end/start points meet. If the line boundaries defining the boundary of the domain have a gap, it can lead to erroneous definitions of domain areas and odd model results. The boundary gap in the example below allows a strip to the left of the gap (white) that is technically not inside that domain (blue); points in this strip have zero boundary crossings to the right. A similar result happens if there is overlap in the location of two end/start points; where they overlap, the algorithm sees two crossings in the positive x direction.

Now consider a model that is still one level (two-dimensional), but is heterogeneous with three subdomains, blue, yellow, and green (below). This example has three types of line boundaries which are labeled: head-specified (hs), normal flux-specified (nfs), and inter-domain (id). Each different line boundary is shown in a different color and is labeled hs, nfs, or id. At all the common intersections where two or more line boundaries meet, the coordinates of the end points must match exactly, so that each subdomain region is defined without gaps or errors.

Along inter-domain boundaries, you specify which domain lies to the left and which domain lies to the right. For example, if the coordinates of the inter-domain boundary between blue and yellow are listed in order from bottom to top, blue is to the left and yellow is to the right. If the coordinates of the boundary between yellow and green is defined in clockwise order, yellow is to the left and green is to the right.

Now imagine that the model shown above is a single level in the blue and yellow areas (two-dimensional), but has four levels in the green area (three-dimensional). A vertical cross section along A-A' for such a model is shown below.

The convention for naming levels in Anaqsim is to start at level 1 at the top and increase the level numbers with depth in a multi-level stack. In the green area, level 1 is at the top and level 4 is at the bottom. Anaqsim assumes that there is vertical leakage between subdomains that exist in the same area but at different vertical levels. It is possible to skip a level number in a stack of subdomains; for example, the green area in the above section could have domains with levels 2,3,5 and 6 rather than levels 1,2,3 and 4. Anaqsim finds the next domains above and below, no matter what the level numbers are and even if there are gaps in level numbers. This is helpful in cases with complex layering schemes such as where a domain has limited extent compared to the domains above and/or below it. The resistance to vertical flow between levels is computed based on the vertical hydraulic conductivities and thicknesses specified for the domains involved.

The brown inter-domain boundary separating the blue and yellow subdomains would have one subdomain on the left (blue) and one on the right (yellow). The green inter-domain boundary separating the yellow and green subdomains would have one subdomain on the left (yellow) and four on the right (shades of green), assuming that boundary is defined with coordinates in clockwise order. Across inter-domain boundaries, there is approximate continuity of head and approximate continuity of the normal component of discharge. This approximation is discussed more in the Line Boundary Conditions section. Interdomain boundaries can connect subdomains in different levels (e.g. level 1 on one side and levels 2 and 3 on the other).

### Left / Right with Respect to Line Boundaries

With interdomain line boundaries and with normal flux-specified line boundaries, it is important to specify boundaries to the left and/or to the right of the line boundary. To explain what this means, refer to the following figure.

Assume that north is up and south is down in this map-view plot of an Anaqsim model. The different colored areas represent different subdomains, which are bounded by various line boundaries. Consider the red interdomain line boundary (id) that separates the blue and yellow subdomains. If the coordinates for this line boundary are in order from south to north, then the blue subdomain would be to the left and the yellow subdomain would be to the right. Think of it as though you were walking along the line boundary from the first vertex towards the last vertex in the coordinates. The "left" side is to your left as you walk along the boundary and the "right" side is to the right. Alternatively, if the coordinates of this interdomain boundary are listed in order from north to south, then the yellow domain is to the left and the blue domain is to the right.

Now consider the green interdomain boundary separating the green and yellow subdomains. If the coordinates of that boundary are specified in clockwise order, the yellow domain is to the left and the green domain is to the right. Alternatively, if the coordinates of that boundary are specified in counter-clockwise order, the green domain is to the left and the yellow domain is to the right.

When a normal flux-specified boundary is external or with a head-dependent normal flux boundary, the vertex coordinates must be specified in counter-clockwise order, with the domain to the left of the boundary as you proceed along it. Consider the purple normal flux-specified boundary (nfs) in the figure. It's coordinates must be specified in counter-clockwise order (west to east) with the yellow domain to the left. Likewise, for the blue normal flux-specified boundary, the coordinates must be specified counter-clockwise from southeast to northwest, keeping the yellow domain to the left. The normal flux is positive for flow across the boundary from left to right as you proceed from the start toward the end (flow out of subdomain). Normal flux is negative for flow from right to left as you proceed from the start toward the end (flow into subdomain).

### Recharge, Leakage and Transient Storage

Like in any flow model, the flow equation in Anaqsim is based on Darcy's Law and conservation of mass (and volume, with constant density). The conservation equation, in its simplest is:

where **∇Q** is the divergence of the two-dimensional aquifer discharge vector field and **γ** is the net extraction per area (sink term, units of L/T). The sink term **γ** may have contributions from leakage out the top of the subdomain (**L _{t}**), leakage out the bottom of the subdomain (

**L**), and transient discharge/area into storage(

_{b}**S ∂h/∂t**) . See equations 4-6 of Fitts (2010).

The vertical leakages **L _{t}** and

**L**are specific discharges proportional to the head difference between the domain and the head above or below (could be specified head or another domain at a different level), and proportional to the equivalent vertical hydraulic conductivity

_{b}**K**that is based on the vertical conductivities and the saturated thicknesses at the average heads specified for the domains involved. If there is leakage between two domains at different levels, the equation used for computing

_{e}**K**is

_{e}where **b1** and **b2** are the average saturated thicknesses of the two layers and **K1** and **K2** are the vertical hydraulic conductivities of the two layers. In Anaqsim, **K _{e}** is assumed to be constant and independent of head and actual saturated thickness.

Under restricted circumstances (uniform recharge, steady flow, and a single-level model),

**γ**is constant and independent of location. In such cases, the uniform

**γ**distribution may be modeled exactly using a uniform area sink.

In many practical cases, the model needs spatially-variable extraction (**γ** varies with **x**, **y**) due to spatially-variable vertical leakage and/or spatially-variable storage changes. When that is the case, the model needs spatially-variable area sinks to approximate the proper distribution of **γ**. The spatially-variable area sink functions in Anaqsim create a smooth, continuous, irregular **γ** surface within a subdomain. The model is using equation 13 of Fitts (2010) as a model of the distribution of **γ**, with **γ** equal to the perfect extraction/area (computed by equation 6 of Fitts (2010)) at each basis point, but approximating equation 6 of Fitts (2010) between basis points. The modeled distribution of **γ** satisfies the flow equation perfectly at basis points and approximates the flow equation between basis points. This approximation is more accurate if the basis point spacing is smaller.

To check the accuracy of this approximation, Anaqsim provides an analysis tool (Analysis Menu / Conditions Along a Line). If checks with this tool reveal a poor approximation, the basis point spacing is too large. If the approximation is fine, you may be able to decrease the basis point spacing and save some computation.

Spatially-variable area sinks can rack up a large number of equations and computational burden on the system, so use them sparingly. Where possible (often the far-field), use a single level in the model, which if it is steady means you can use the very efficient uniform area sinks instead of spatially-variable area sinks. Use the special well basis point spacings around wells to gain accuracy with minimal computation.

### Elements

The analytic elements used in Anaqsim include the elements described in the following references:

- The standard well element described by Strack [1989] is used in isotropic subdomains and the well element of Fitts [2006] is used in anisotropic subdomains.
- Line boundaries are represented by high-order line elements similar to those described by Jankovic and Barnes [1999]. Line elements are either linesinks, line dipoles, or line doublets. Line elements may have up to ten unknown parameters (degrees of freedom). For example, a linesink with one parameter has a constant discharge/length along its length, a linesink with two parameters has discharge/length that varies linearly along its length, and a linesink with three parameters has discharge/length that varies parabolically along its length. Strack [1989] gives a good overview of line elements.
- Uniform area source/sinks (constant over a subdomain) are modeled with the uniform recharge functions described by Strack [1989]. For these, the extraction/area is uniform over the entire subdomain. This is appropriate in cases with a single level, steady flow, and a uniform recharge rate or zero recharge.
- Spatially-variable area source/sinks are modeled with the multi-quadric radial basis function elements described by Strack and Jankovic [1999]. With these, the extraction/area varies with location. Spatially-variable extraction occurs when the extraction = recharge + leakage + transient storage flux is spatially variable. This is generally the case in multi-level models where spatially-variable leakage occurs and in transient models where spatially-variable storage flux occurs. See the previous section and Fitts (2010) for details of how these elements are used to represent recharge, leakage and storage flux.

### Wells

There are three possible kinds of pumping wells in Anaqsim:

- Discharge-specified wells in just one subdomain. With this type of well, you specify the known discharge of the well.
- Discharge-specified wells spanning multiple subdomains. This is for wells that are screened across multiple levels and subdomains in a multi-level part of a model. For example, a well could be screened across the two lowest levels in a 4-level part of a model. You specify the known discharge of the well, and Anaqsim imposes these conditions: a) the heads at the well radii in each screened subdomain are identical, and b) the total discharge of the well elements in the screened subdomains equals the specified discharge of the well.
- Head-specified wells in just one subdomain. With this type of well, you specify the known head at the well radius. Anaqsim determines the discharge needed to meet this condition.

### Interdomain Line Boundary Conditions

At the inter-domain boundaries, two different boundary conditions are enforced by Anaqsim. First, heads are matched across the boundary at specific points in all subdomains that intersect the boundary. For example with the green inter-domain boundary shown in the figures below, heads are matched in the yellow domain and the four green domains at control points along the boundary.

*Plan view (above) and section view (below) of a multi-domain model.*

Secondly the normal component of discharge is matched across line segments on the inter-domain boundary. For the green inter-domain boundary, this means that the sum of the normal components of discharge in the four green subdomains must equal the normal component of discharge in the yellow subdomain. This condition is enforced across intervals on a line segment. If each line segment on the inter-domain boundary has n parameters, then the condition is enforced across n intervals per line segment.

A good way to check that you have enough parameters specified on any line boundary condition is to use right-click/ Check Line Boundary Conditions, which creates a plot of the boundary condition accuracy along a line segment of the line boundary. If you don't have satisfactory accuracy, you can either (1) go under Line Boundaries and specify more parameters (control points) per boundary line segment, and/or (2) break a long boundary line segment into multiple, shorter segments.

One consequence of the head-matching condition is that there will be no vertical hydraulic gradient between different levels on the multi-level side of an inter-domain boundary. This is consistent with the fact that there is no vertical resistance to flow accounted for on the single-level side of an inter-domain boundary. Because of this, it is not appropriate to use one inter-domain boundary with multiple levels on both sides of the boundary; doing so would remove vertical head gradients at the boundary, which defeats the purpose of having multiple levels.

With inter-domain boundaries, there should generally be just one domain (level) on one of the two sides of the boundary. Anaqsim checks the input and gives a warning if there are two or more levels on both sides of an inter-domain boundary. To explain why, consider the case illustrated below, with two levels on the yellow side of the inter-domain boundary(s) and four on the green side. If there was just one inter-domain boundary there, it would be as though an infinitely conductive thin vertical boundary were inserted between the yellow side and the green side, and heads would match in all six domains that meet there, meaning there can be no vertical head gradient on either side of the boundary. The total normal component of discharge would match (total normal discharge in the two yellow domains would match the total normal discharge of the four green domains), but that normal discharge may be distributed oddly between different levels (e.g. a lot of flow from the lower yellow level going into the uppermost green level). A better solution would be to use two inter-domain boundaries:

- one with the lower yellow domain on one side and the two lowest green domains on the other, and
- one with the upper yellow domain on one side and the upper one or two green domains on the other.

This way, there can be vertical gradients between the yellow/green interface, and upper-level flows are matched with each other and lower-level flows are matched with each other.

### Pathline Tracing

Pathlines are traced in the horizontal plane using the aquifer discharge vector function in the domain and a numerical tracing procedure outlined in section 26 of (Strack, 1989). Three-dimensional pathline tracing is done using a finite-difference form of equation 24 in the paper by Strack (1984), which uses three-dimensional flow continuity to approximate the vertical component of pathlines. For transient models, pathlines are traced through a transient flow field that evolves through the duration of the transient simulation (this is new in release 2016-2; earlier releases traced pathlines through the flow field of the last time step in a transient model). From release 2016-2 on, pathlines in transient models start at the specified start time, and can be traced backward to the beginning of the simulation or forward until the end of the simulation. In steady models, pathlines are traced forward or backward until they either exit the model or until the specified total time is reached.

Where pathlines cross inter-domain boundaries, the elevation of the pathline is adjusted so that the fraction of the normal component (perpendicular to the boundary segment) of flow above and below the pathline on each side of the inter-domain boundary match. For example, say you have an inter-domain boundary with one domain on the left and three on the right. If the pathline comes to the boundary from the left with 45% of the normal component of discharge above the pathline and 55% below it, the domain and elevation of the pathline side will be determined so that this 45/55 ratio is preserved on the right side as well. There is typically a jump in the elevation of the pathline as it crosses an inter-domain boundary for the following reasons:

- The elevations of the domains on left and right do not necessarily match
- Where there are multiple levels on one side of the boundary, the distributions of the normal components of flow vary from one level to another due to K and flow field variations
- The total normal components of flow match only approximately from left to right side (see Fitts (2010)). Anaqsim is not matching these normal components of flow perfectly at every point on the boundary - it is matching the total normal discharge across segments of the boundary perfectly.

When the normal discharge components across an inter-domain boundary are in different directions in different levels (usually where flow is nearly stagnant) this algorithm breaks down, and under such conditions pathlines may terminate at the inter-domain boundary.

When pathlines intersect an internal line-sink boundary such as a river boundary, there are two possibilities: 1) the pathline could be consumed by the linesink and terminate, or 2) the pathline could continue on underneath the linesink. Whether option 1 or 2 occurs is determined based on the elevation of the pathline as it approaches the linesink, the discharge/length of the linesink at the intersection, and the normal component of the domain discharge at the intersection. It is assumed that the linesink is at the top of the domain and draws its water from the upper portion of the domain discharge. If the result is case 2 with the pathline continuing under the linesink, its elevation will generally jump as it crosses under the linesink, gaining elevation if the linesink is extracting water and losing elevation if the linesink is injecting water.

In multi-level areas of models, pathlines can cross from one level to another vertically.

Where pathlines exit the bottom or top of the model due to recharge or leakage (e.g. at the water table for pathlines traced upstream), Anaqsim draws a circle at the end of the pathline and if you move the cursor over the circle, it will display elapsed time, domain, elevation, etc.

# User Interface

# Overview

The user interface consists of one main menu plus three different tabs:

- Plot: for graphically displaying map-based inputs and outputs
- Data: for editing input data. The data view is automatically displayed when you choose to edit something under the Model Input, Plot Input, or Analysis Input menus.
- Log: for displaying text outputs. The run log displays information when files are opened or saved, when the Solve menu is executed, and in response to a number of choices under the Analysis menu.

More information about each of these tabs is given in the following pages. Change from one tab to another by clicking on the tab at the left. The current active tab is highlighted in blue. In the following image, the Plot tab is selected.

# Plot Tab

The Plot tab shows the model inputs and results plotted in map view. This view is automatically shown after opening an existing model, and after making a plot under the Make Plot menu. Most of the Plot view is a map view of the model that can display a basemap, model elements, simulated heads, flow vectors, pathlines, etc. The plot area has scroll bars that allow you to shift the view left/right and up/down while retaining the same scale. The scroll wheel on most mouse devices will cause the view to zoom in and out. Pressing down on the scroll wheel and moving the mouse will allow you to pan the view with most mouse devices.

On the upper left is a separate plot view menu (Plot File, View Manager,...Zoom Previous) that applies just to the plot. Choices in this menu allow you to zoom to a different view, save or print the plot, digitize and edit coordinates, or add annotations to the plot. See the tutorial videos at the website to see the various functions in the plot view menu explained and demonstrated. Right-clicking your mouse over the plot brings up a context menu with many choices for digitizing, editing line boundaries, and generating outputs related to the cursor location.

### Context Menu (right click over plot)

When the cursor is over the plot, you may bring up a context menu by clicking your right mouse button. This menu allows the following options, which are handy short cuts while building, editing, and analyzing a model:

- Edit Domain Properties - switches to the data view with a domain data table displayed and one highlighted row: the row of the domain where the cursor is located.
- Edit Nearest Well - switches to the data view with a well data table displayed and one highlighted row: the row of the well nearest the cursor location.
- Edit Nearest Line Boundary - switches to the data view with a line boundary data table displayed and one highlighted row: the row of the line boundary nearest the cursor location.
- Edit Area Sink that Applies Here - switches to the data view with an area sink data table displayed and one highlighted row: the row of the area sink that applies at the cursor location.
- Digitize - includes the most commonly used commands (Point, Polyline, Clear Digitizing Marks) from the Plot View/Digitize menu, but without the pop-up instructional windows. These are efficient for experienced users.
- Edit Line Boundary - allows the same operations (Insert Vertex, Delete Vertex) as are found in the Plot View/Digitize menu, but without the pop-up instructional windows. These are efficient for experienced users. One line boundary, spatially-variable area sink (SVAS) polygon boundary, or vertical leakage polygon must be selected before executing either of these menu items.
- Set Plot Window to Current View - this is the same as Plot Input/Set Plot Window to Current View, setting the window for subsequent plots to the current view.
- Set Plot Window to Entire Model - this is the same as Plot Input/Set Plot Window to Entire Model, setting the window for subsequent plots to the entire model.
- Check Nearest Well Head and Discharge - this writes the head and discharge of the nearest well to the Log tab.
- Check Nearest Line Boundary Condition - this allows you to check the accuracy of the approximation of boundary conditions along the particular segment of a line boundary that is closest to the cursor when this is selected. This causes a graph to be made of the conditions on that boundary segment. The graphs vary depending on the type of line boundary. You may export a bitmap graphic of the graph or the underlying data (see the Exporting X-Y Graphs topic). Prior to release 2016-2, this feature was under the Analysis menu.
- Graph Head Hydrograph(s) Here, All Levels - for a transient simulation, this causes hydrographs (graphs of head vs. time) to be made that show heads in all model levels at the location of the cursor. These graphs do not contain the initial heads (time=0) at these locations. If you want hydrographs that include initial heads, use Analysis Input/Hydrograph Points and Analysis/Head Hydrographs and Analysis/Drawdown Hydrographs.

### Data at Cursor Location

To the left of the plot and below the plot menu is an area that displays the coordinates of the cursor, model information and model results (X,Y Coordinates, Model Level, Domain Name, Head,...) at the location of the cursor. As you move the cursor, this data updates based on the cursor location.

Each item may be hidden or visible; to toggle this, click on the triangle to the left of the label. If the window is not tall enough to display all items, use the scrollbar immediately to the right of this area.

In transient models, the values reflect the time step of the selected period, which ends at the time listed. Head and interface values are at the end of the time step at the time listed. Discharge-related values apply over the duration of the time step listed.

The model results shown are described below.

**Head Above - Head**is the difference in head from the level above to the level of the plot.**Head - Head Below**is the difference in head from the level of the plot to the level below.**Interface Elevation**is the elevation of the fresh-salt interface (only applies in models with fresh-salt interface domains).**Domain Discharge**is the specific discharge times the saturated thickness, which equals the discharge in the domain per unit length normal to the discharge.**Flow Direction**is the direction of the specific discharge or average linear velocity vector, assuming the positive x axis is zero degrees.**Top of Model Condition**is the condition at the top of the topmost domain of the model at this x,y, defined by the area sink input.**Bottom of Model Condition**is the condition at the bottom of the bottommost domain of the model at this x,y, defined by the area sink input.**Modeled Extraction**is the extraction per area**γ**[L/T] that Anaqsim is modeling at the point, using the spatially-variable area sink functions. This quantity is defined by equation 13 of Fitts (2010).**Extraction from Heads**is defined as:

where**L**is leakage out the top of the subdomain computed using head differences and vertical hydraulic conductivities,_{t}**L**is the leakage out the bottom of the subdomain computed similarly, and_{b}**S ∂h / ∂t**is the transient discharge/area into storage computed using the head change over a time step and storativity. At spatially-variable basis points, this equals the modeled extraction, but between basis points, the two are unequal, but close if the basis point spacing is small enough. This quantity is defined by equation 6 of Fitts (2010). In a transient simulation, this cannot be computed for the first time step, since the initial head at time zero at the cursor location is not known and therefore**∂h**can't be computed.**Leakage out Top**is the leakage out the top of the domain,**L**as defined above [L/T]._{t}**Leakage out Bottom**is the leakage out the bottom of the domain,**L**as defined above [L/T]._{b}**Storage Flux**is the flux into storage in a transient model,**S ∂h / ∂t**, as defined above [L/T].**Leakage Factor Above**is the leakage factor (**LF**) computed for upward leakage from the subdomain. If the subdomain is at the top of the model with no overlying level, the leakage factor is computed as:

Where**T**is the transmissivity of the subdomain,**b**is the saturated thickness of the subdomain,**K**is the vertical hydraulic conductivity of the upper half of the subdomain. If the subdomain has another subdomain overlying it, the leakage factor is computed with the following equation:_{v}

Where

**T**is the transmissivity of the subdomain above,_{a}**b**is the saturated thickness of the subdomain above, and_{a}**K**is the vertical hydraulic conductivity of the lower half of the subdomain above. Leakage factors are used as guidance for determining appropriate basis point spacing for spatially-variable area sinks. See the discussion under spatially-variable area sinks._{va}**Leakage Factor Below**is the leakage factor (**LF**) computed for downward leakage from the subdomain. It is computed in a manner analogous to that described above for Leakage Factor Above.

# Data Tab

The data tab allows you to edit model inputs, plot inputs, and analysis inputs with a data grid component that displays data from the underlying data tables. To display and edit a table, make a selection under the Model Input, Plot Input, or Analysis Input menus.

There is a context menu that pops up when you right-click over the grid. Options in this menu include

- Paste New Rows - this pastes in new rows of data that are tab-delimited between columns. This format allows you to paste data copied from spreadsheets like Excel.
- Copy Selected Rows - this copies the selected rows to the system clipboard, which can then be pasted into spreadsheets like Excel.
- Copy All Rows - this copies all rows in the spreadsheet to the system clipboard, which can then be pasted into spreadsheets like Excel.

### Using the Data Grid

A data table is displayed in the grid when you select an item under the Model Input, Plot Input, or Analysis Input menus. The displayed data is linked to one of several database tables, and when you edit the displayed data, the underlying data table is updated.

The table of data is displayed with headers that define each column, such as Label, Domain, Parameters_per_line, ... as shown below.

You can move from cell to cell with the arrow keys or using the mouse. The current row is highlighted blue. You can enter new values by navigating to a cell and typing a new entry, or you can double-click on a cell to edit cell contents with a text editor, as shown below, where the contents of the cell with "101" is being edited.

When you enter a value in a grid cell, the underlying database is updated when you press enter or move to a different cell. At this point the cell value is checked to make sure it is compatible (e.g. a positive real number for hydraulic conductivity). If the value is incompatible, an error message is displayed and you must correct the cell entry. Be sure to remember to press enter or move to another cell after editing the value in a cell, otherwise the value will not be changed in the database.

New rows are created by editing the blank row at the bottom of the table. A new row of data is entered into the underlying data table only when you press enter after editing the row, at which point a new blank line appears below the line just entered. The following two screen shots shows a new 2nd row before it has been entered in the database (no blank row shows below it), and after (blank row below 2nd row).

In data tables that contain multiple rows, the leftmost field is often called Label, and it is always displayed even if you scroll far to the right. No entry is required in this field, and it accepts any text. It is wise to fill in a text label in this field (e.g. “PW-103” for a pumping well). The label will help you know which feature this row represents, and many analysis outputs make use of this label. Also you can sort the data based on entries in this column to easily find the row you want. The contents of the table can be sorted by clicking on the column header. Clicking a second time reverses the sort order. It is a good idea to choose labels that easily allow you to sort features. For example, if you want to easily find a group of wells on property A, you could give them labels such as "A_MW102", "A_MW105", "A_MW113"... so they would be grouped together after sorting by the label column.

Column widths are automatically adjusted to fit the contents. You can increase or shrink column widths by dragging the left or right the vertical line that separates columns in the header (top) row. Double-clicking on this vertical line automatically resizes the column width to fit the contents.

Some columns, like the Parameters_per_line column in the table shown above, are edited using a drop-down list of choices. To see the list, double-click the cell, then select the item you want.

Other columns, like the Coordinates column, contain buttons to edit or select data; these are edited by clicking on the button.

### Number Formats

All data grid cells that expect numerical input have common format constraints. You can input real numbers with formats such as the following:

- 1256.82
- 0.004
- 1.4e-2

The last one is scientific notation for 1.4x10-2.

You should not insert commas to mark thousands, millions, (e.g. 1,200,000) as the comma may be interpreted as a decimal mark. In North America, the convention is to use a period for the decimal marker. In Europe, the comma is often used as a decimal marker. There is a Windows operating system setting to switch between these modes. Often, European users need to adjust these settings to use Anaqsim.

### Editing Coordinates

In many of the data input tables, there are columns and cells that display an "Edit" button in the Coordinates column. When you click the button, a text box window pops up and you enter coordinate data there:

Often, you will digitize the coordinates in the plot tab and then paste the coordinates into this text box window. Alternately, you can just type coordinates in. The OK button records the edited coordinates and the Cancel button does not.

### Deleting Data Rows

Delete one or more rows of data in the data table by selecting rows and then pressing the Delete key. Row(s) are selected by clicking (and dragging for multiple rows) in the leftmost column of the grid. A dialog will ask you if you really want to delete those records from the data table.

### Importing and Exporting Data

To import data from Excel into a data table, highlight a block of data in an Excel sheet that corresponds to row(s) of data in a data table, copy that block in Excel, then right-click over the data grid and select Paste New Rows. This will add these copied rows to the data table. Make sure that the columns in the copied block match the columns in the data table. Data in Coordinates columns cannot be be pasted in due to their multi-line structure, but all other columns can be pasted in. In the case of a Coordinates column, a paste operation leaves that blank and you must enter the coordinates by clicking on the Edit button in that column.

To export rows of data to Excel or a text file, select rows of data (see section above) and then right-click over the data grid and select Copy Selected Rows. After doing this the rows of data are in the computer’s clipboard as tab-delimited data, which can then be pasted into Excel or into text files.

# Log Tab

The Log tab holds the run log, which is an area that displays text output from the program. The run log continues to accrue more text as you execute various tasks such as updating license information, opening a file, solving the system of equations, checking boundary conditions, checking calibration results, or closing a file. If the text in the run log gets long enough, a scroll bar will appear to let you scroll through the entire log, as shown below. You can select all or a portion of the text in the run log and cut, copy, and paste this text. This is an easy way to move text results to another document.

# Menu Keyboard Shortcuts

You can access common menu items with keyboard shortcuts by pressing the key sequences as listed below. Many are standard Windows shortcuts.

Ctrl-O | File/Open |

Ctrl-S | File/Save |

F12 | File/SaveAs |

Ctrl-W | File/Close |

Alt-FE | File/Exit |

Alt-S | Solve |

Alt-PA | Make Plot/All Selected Features |

Alt-PE | Make Plot/Elements Only |

Alt-V | Switch View |

Alt-H | Help |

# General Modelling Sequence

Creating a model follows this general sequence:

- You can start a new model either right after starting the program or after selecting File/Close which closes the current input and begins a new input data set. Once either of these steps is taken, you may edit data tables under the Model Input, Plot Input, or Analysis Input menus.
- Create model input using the Model Input menu. Make sure to define Domains before adding Well, Line Boundary, or Area Source/Sink elements. This sequence is necessary because the input for the elements includes specification of the domain(s) they are in. When adding elements, it helps a lot to use a basemap and digitize coordinates on top of the basemap.
- Define what you want displayed in plots with the Plot Input menu.
- Define what analysis features you want with the Analysis Input menu.
- Save your model frequently as you build your input.
- When the model input is complete, select Solve to solve the system of equations. This is required after making any model input changes and before making output plots or using the Analysis menu.
- View plots of the model results with Make Plot.
- Examine model results with the Analysis menu.
- Loop back through steps 2-8 to revise the model, re-solving the system after revision and before examining results.

# Menus

# File Menu

### Open

This selection opens a dialog that allows you to find and open existing input files (.anaq extension). These files store the data you edit under the Model Input, Plot Input, or Analysis Input menus in XML file format. XML is a common ASCII database file format. You could edit these directly with a text or XML editor, but that is not recommended since it risks corrupting input with improper values or format. When you open a file, the layout of elements in the model is drawn to the plot view.

If you want to be able to open .anaq files by double-clicking on them in Windows Explorer, in addition to opening them from the File/Open menu, the Windows operating system must associate .anaq files with Anaqsim. In case this association was not established during installation you can manually do it with Windows Explorer. To do this, locate a .anaq file in Windows Explorer. Right click on the file and then select Open with, then select Chose default program. In the dialog that pops up, check the box next to Always use the selected program to open this kind of file and then select the Browse button and browse to find the Anaqsim.exe file in the Program Files / Fitts Geosolutions / Anaqsim software directory. Now Windows will associate .anaq files with Anaqsim.exe, and you can open any .anaq file directly from Windows Explorer by double-clicking on it.

### Save, Save As

The Save As option brings up a dialog that allows you to save your input to a file with a new name. Using Save saves the input to the same file name. If you have yet to save input and have no filename, it will function like Save As. When you save, you save the input data tables to an XML format database file with the .anaq extension.

### Close

This closes the input you are working on and clears all the associated data tables in memory. After selecting Close, you may begin editing a new model.

### Save locations for Initial Transient Heads

A transient model needs initial heads so it can compute the head change that occurs during the first time step. These values are needed at the location of each basis point in each spatially-variable area sink, which account for storage fluxes. The initial head values come from a pre-existing model, which could be steady or transient. Initial heads are also retrieved for discharge-specified wells, hydrograph points, and transient line conditions (see Analysis Input Menu for the last two items).

The initial conditions model must have the same number and extent of layers as the transient model, at least in areas with basis points, wells, hydrograph points, or transient line condition lines. When the initial head locations are written, each location is identified by its x,y coordinates and its layer. This is new in release 2020-1; prior releases wrote the x,y coordinates and the internal domain number. The change made in release 2020-1 allows different domain configurations between the initial and the transient models, which can be helpful. Because of this change, you must not mix initial head location or initial head files created prior to release 2020-1 with a release 2020-1 or later model. To avoid incompatibility when you switch to release 2020-1, re-create the initial head location file and the initial heads file as outlined below using release 2020-1.

To create a transient model that has a proper set of initial heads, these steps are necessary:

- Make sure the transient model you begin to create has been saved to it's own unique file name, different from the file that contains the input for the model that will provide the initial condition heads.
- Set up the transient model (i.e. uncheck Steady under Model Input/General, establish the time step sequence under Model Input/Time Steps, adjust input for the transient case, set up all spatially-variable area sinks, make the appropriate settings under Analysis Input/Hydrograph Points and Analysis Input/Transient Line Conditions, etc.
- Select File/Save Locations for Initial Transient Heads, which saves the locations for transient starting heads from the transient model (these are the locations of all basis points associated with spatially-variable area sinks in the transient model and locations of hydrograph points, wells, and transient line conditions. This saves the level and the coordinates of each of these to a binary file with the .ihl extension.
- Close the transient model.
- Open the initial conditions model (the one that represents conditions at the start of the transient run). Solve it. Select File/Write Initial Transient Heads, click on the initial heads file you would like to use for this simulation, and click Open. This reads in the locations saved from the .ihl file saved in step 3.
- A dialog opens asking you to name the binary initial heads file. The default name is the same as the initial conditions input file but with the .hds ending. Anaqsim then writes the initial heads at the locations in the .ihl file to a binary file with the .hds file extension.
- Close the initial conditions model.
- Open the transient model. After checking that all model parameters are set correctly for the transient run, select Solve. At this point Anaqsim will ask you to select the .hds file containing the initial heads created in step 6.

When you solve the transient model, the heads are read in and used to determine the head change in the first time step at each basis point.

### Write Initial Transient Heads

See the discussion under Save Locations for Initial Transient Heads for an overview of setting up initial heads for transient models.

### Save Solution

This allows you to save the model solution after you have solved. Later, you can open the model input file, then load the saved solution, and avoid the "Solve" step. This is particularly handy for large models that have longer solve times, and allows you to save your solution and come back to it later for making plots or doing analysis of the solution. All model objects with their strengths are saved in a binary format, to a file that has the same name as the input file, but with the ".solu" filename extension instead of ".anaq".

### Load Saved Solution

This allows you to load in a previously saved model solution. To make use of this, first open the model input file for the model, then load the saved solution, which avoids the need for the "Solve" step. This is particularly handy for large models that have longer solve times, and allows you to save your solution and come back to it later for making plots or doing analysis of the solution. All model objects and their strengths are read in from a binary file that has the same name as the input file, but with the ".solu" filename extension instead of ".anaq".

### Export Input Data to Excel File

This causes the entire input database to be written to one excel file in the same directory as the input file (*.anaq), written to a file with the same name but the excel suffix (*.xlsx). The Excel file has multiple sheets, one for each data table. Each sheet contains the same headers as the data tables plus all rows of input. This is a handy way to document model inputs all in one readable file.

### Exit

This exits Anaqsim. A dialog asks if you want to save the current input before exiting. The same is achieved by clicking on the red "x" at the upper right corner of the Anaqsim window.

# Edit Menu

This is like edit menus in most other Windows applications with Cut, Copy, and Paste menu choices. These functions are also available with the usual keyboard shortcuts: control-x (cut) control -c (copy) control -v (paste).

# Model Input Menu

### General

This item has only one line of input. The first item is checked if the model is steady-state and not checked if the model is transient. To simulate storage fluxes in transient models, spatially-variable area sinks (SVAS) are required. If you try to solve a transient model without any SVAS, Anaqsim gives an error message.

The other three items are text values that document the length and time units used in the model, and provide comments to document the run. The model uses consistent length and time units throughout. For example, if you chose meters and days, then hydraulic conductivity, specific discharge, and average linear velocity are in m/day, well discharges are in m3/day, and time markers on pathlines are in days.

### Solution

The two data tables under this menu define settings involved in solving the system of equations in your model. The first defines parameters involved in solver iterations, and the second lays out the solution accuracy needed before iteration ceases.

#### Solve Settings

- Underrelaxation is a parameter that governs how unknown strength parameters are updated at each iteration. If underrelaxation is 1.0, the new strength parameter equals the parameter from the most recent iteration. If this value is 0.7, the new strength parameter is weighted 70% by the parameter at the most recent iteration and 30% by the parameter at the previous iteration. Lower values help damp out oscillations in parameter values from iteration to iteration and may improve convergence in some situations. Higher values close to 1.0 speed convergence when oscillation is not a problem. A good range for most cases is 0.9 to 1.0.
- Maximum Iterations. When solving, iteration continues until the tolerances specified in Check Settings are met at all boundary condition control points or basis points, or if those criteria are not met, iteration ceases when this maximum number of iterations is reached.
- Starting Heads. For transient runs, the first time step needs initial heads from some source. There are two possibilities here:
- Assign a constant value of head for each domain. The starting head at a point is set to the domain's average head, which is specified under Domains input. In the case of a Discharge-Specified (Multi-Domain) well, the starting head is set to the average head of the first domain listed. This option should only be used for very simple transient runs where the initial condition is a uniform, flat, constant head.
- Read initial heads from a file. This file is written by the initial (time zero) model, which is another Anaqsim model of the same area that may be steady or transient (see discussion of this file under the File menu).

- Almost_dry_fraction. This parameter affects the solution in cases with unconfined, confined/unconfined, unconfined interface, and confined interface domains, where the heads can drop to low enough levels that the freshwater saturated thickness of the domain approaches zero. Instead of letting the domain have zero or negative saturated thickness, Anaqsim has a minimum saturated thickness. In an unconfined, confined/unconfined, or unconfined interface domain, the minimum saturated thickness = (average head - base elevation) * Almost_dry_fraction. In a confined interface domain, the minimum saturated thickness = (top elevation - bottom elevation) * Almost_dry_fraction. When heads drop to or below a level that would create this minimum saturated thickness, the aquifer is treated like a confined aquifer with this fixed minimum saturated thickness. This prevents domains from actually going completely dry and helps models converge despite portions of some domains approaching "dry" conditions. Setting
- Almost_dry_fraction to a low value means that very little simulated horizontal discharge will occur in "'dry" portions of the domain. Heads in these "dry"' areas will drop to unrealistically low levels (below the base of an unconfined domain, for example), and have little real meaning. When contour plots are made of head, these unrealistically low values are neglected. These low heads do appear in other outputs such as the panel to the left of the plot and in profiles. Using higher values of Almost_dry_fraction may help convergence by limiting the magnitude of head gradients in "dry" areas, at the expense of allowing more actual discharge in "dry" areas. You can check the amount of discharge in a "dry" area by using the Analysis/Conditions Along a Line and making a profile of domain discharge along or across a line. Chose a low enough value of Almost_dry_fraction so the discharge in these "dry" areas is acceptably small.
- Interface_leakage_option. This specifies one of two options for computing vertical leakage at a spatially-variable area sink basis point where a fresh/salt interface is present. If unchecked (default), the vertical leakage is computed just like it is in all other cases: the vertical leakage rate is proportional to the difference in head from one level to the next. If this column is checked and an interface is present in the overlying layer, the head difference is computed by using the freshwater head that is at pressure equilibrium with static salt water at the base of the overlying layer; the head in the overlying layer is not used to compute the head difference. The head that is used to represent the overlying layer is the same as the head at the toe of the interface in the overlying layer. The unchecked mode is appropriate when the resistance to vertical flow between levels is due to the K3 of the domains themselves. The checked mode is appropriate when the resistance to vertical flow between levels is due to an aquitard between the levels that is not explicitly modeled (in this case the resistance of the aquitard must be incorporated into the K3 values specified in the upper and/or lower domain). See Fitts et al (2015) for comparisons of these methods and discussion.

#### Check Settings

These settings define the accuracy of boundary conditions required of the solution; iteration continues until these conditions are met. If when Solve is pressed the solution converged before reaching the maximum number of iterations, all boundary conditions were met within the tolerances specified here.

These settings also affect the function of Analysis/Check Boundary Conditions at Latest Iteration, which is used to check how well the solution meets specified boundary conditions. Such conditions include heads at head-specified wells and linesinks, extraction at spatially-variable area sink basis points, etc. When you select Analysis/Check Boundary Conditions at Control Points, each boundary condition is checked, and if the discrepancy between the specified condition and the model-simulated condition is greater than a threshold you specify here, the program prints the discrepancy to the run log. If the discrepancy is less than the threshold, nothing is printed. In cases where the solution did not converge to within these tolerances during the Solve process, the offending boundary conditions are listed along with their accuracy. This can be a help to home in on which boundary conditions are holding up the Solve process.

Four kinds of boundary condition tolerances are defined as follows.

- Head_check_tolerance is the threshold for the magnitude of discrepancy between specified and modeled heads, used at head-specified wells and line boundaries.
- Qn_check_tolerance is the threshold for discrepancies in the computed discharge per length along River line elements. It units are discharge/length [(L
^{3}/T)/L = L^{2}/T]. - Q_check_tolerance is the threshold for discrepancies in discharge, which are used along inter-domain line boundaries and normal-flux specified boundaries. In both cases, the condition being checked is the total discharge across a segment along the line [L
^{3}/T]. For inter-domain boundaries, it is the comparison of discharges on opposite sides of the boundary. For normal-flux specified boundaries it is the difference between modeled and specified discharges based on the specified discharge per length times the length of the segment. - Extraction_check_tolerance is the threshold for discrepancies in the modeled extraction (Equation 13 of Fitts (2010)) and the extraction computed from head values (Equation 6 of Fitts (2010)). The units of extraction and this threshold are discharge/area [(L
^{3}/T)/L^{2}= L/T].

### Time Steps

This input is used only for transient models. If you are doing a transient model, make sure you uncheck Steady under Model Input/General. Each row of input in the Time Steps table defines a time period, during which all boundary conditions are constant. For example, a model could have three time periods with different recharge rates, river stages, or well discharge rates in each of the three periods, but within each period the values remain constant.

For each time period, you specify the total length of the time period (Period_Length), the number of time steps the period is divided into (Steps_in_period), and the time step multiplier (Step_multiplier). The multiplier causes the length of successive time steps to grow by a factor equal to the time step multiplier. The following table illustrates the lengths of time four time steps for a period that is 100 time units long, using various time step multipliers.

Time Step | Multiplier = 1.0 | Multiplier = 1.5 | Multiplier = 2.0 |

1 | 25.00 | 12.31 | 6.67 |

2 | 25.00 | 18.46 | 13.33 |

3 | 25.00 | 27.69 | 26.67 |

4 | 25.00 | 41.54 | 53.33 |

Total Time: | 100.00 | 100.00 | 100.00 |

In all cases, the total time of the period is 100, but the lengths of the four steps change as the time step multiplier changes. This scheme is the same as employed in MODFLOW. Using a multiplier larger than 1.0 helps concentrate computing power early in the time period when there is more transient change occurring. Transient storage fluxes, which are part of spatially-variable extraction, are computed for each time step using a finite-difference approximation of the governing equation (Equation 6 of Fitts, 2010).

### Domains

The properties of each domain (called a subdomain in Fitts, 2010) are set with data tables under this menu. Different tables define the properties of different kinds of domains. A domain is a polygonal region of the model in a certain model level. Inside a domain, the aquifer properties (hydraulic conductivities, base elevation, storativity, porosity, etc) are homogeneous.

The boundary of a particular domain is defined by a combination of line boundaries that, in their input, are listed as external to the domain. Head-specified, normal flux-specified, and inter-domain boundaries can be external boundaries for domains. Other line boundaries like river and discharge-specified line boundaries are internal to domains and do not define domain boundaries. See the discussion under Subdomains and Model Levels for more detail and some examples.

All domain input data tables may be accessed through the main menu or by using a pop-up context menu when the cursor is over the plot.

#### Boundaries of Domains

The geometry of the boundary of each domain is not specified in the domain data tables, but is determined by the distribution line boundaries that define the external domain boundaries. Line boundaries that can be external domain boundaries are head-specified, normal flux-specified, and inter-domain. The data that is input for these types of line boundaries include information about which domain(s) that they bound. All domains should be completely bounded by such line boundaries, so that their geometry is unambiguous. See the discussion under Subdomains and Model Levels for more detail and some examples.

For the best accuracy, make sure the coordinates of the starting and ending points of adjacent external line boundaries match exactly (copy them). For example, if a domain boundary has two line boundaries defining it - one head-specified and the other inter-domain, make sure that the start/end points where these two line boundaries join have the exact same coordinates.

#### Input Common to All Domains

Several parameters are common and required input for any type of domain:

- The label_unique allows you and the model keep track of which domain is which. When adding wells, line boundaries, etc. to the model, this label defines which domain these features are in. These labels are required and must be unique (no two domains should have the same label). Once the labels have been declared and other features have been added to the model, do not change the labels because doing so would require changing the domain label of each well, line boundary, and area source/sink that is in that domain.
- The level of a domain refers to where this domain fits in the vertical sequence of model levels. In a multi-layered part of a model, the level begins at 1 at the uppermost level and increases in deeper levels. The level may be from 1 to a maximum of 15. If there are vertically stacked domains in a multi-level part of the model, vertical leakage is assumed to occur between domains that are at different levels but occur at the same x,y coordinate. For example, if three domains in a three-level part of the model were assigned levels 1, 2, and 4, there could be vertical leakage between levels 1 and 2 and between levels 2 and 4, if this area of the model has spatially-variable area source/sinks. When making plots of model results, plots show one level at a time.
- For Average_head, list your best estimate of the average head in this domain for your simulation. This then defines a constant that is added to the discharge potential for this domain. More details are given about this in the next topic.
- The value of Porosity is used in computing average linear velocity and advection travel times along pathlines. The average linear velocity = specific discharge/porosity. For a solute that adsorbs to the porous medium, the solute travel time will be be longer than the pure water advective travel time. If you want to plot travel times that factor in retardation of a solute, you may input Retardation Factor * Porosity instead of Porosity in this column and the travel times will reflect travel of a retarded solute plume.
- K1_horizontal is the horizontal hydraulic conductivity. In a domain that is anisotropic in the horizontal plane, K1 differs from K2, and these represent the principle hydraulic conductivities in the horizontal plane. K1 > K2, K1 < K2, and K1 = K2 are all possible.
- K2_horizontal is the horizontal hydraulic conductivity. In a domain that is anisotropic in the horizontal plane, K2 is in the direction normal to the K1_horizontal direction. To simplify inputs, you can enter "=K1" if K2=K1 and the domain is isotropic in the plane of the domain. If you want a fixed ratio of anisotropy, you can enter "=K1*D" where D is a real number. Using "=K1" or "=K1*D" is particularly handy for parameter estimation, to limit the number of parameters being estimated.
- Angle_K1_to_x is the angle, in degrees, between the x axis and the direction of K1_horizontal. Positive angles are measured counter-clockwise from the x axis.
- K3_vertical_top defines the vertical hydraulic conductivity of the upper half of the domain. This parameter is only used if there is vertical leakage with spatially-variable area source/sinks. To simplify inputs, you can enter "=K1" if K3=K1 and the domain is isotropic normal to the plane of the domain. If you want a fixed ratio of anisotropy, you can enter "=K1*D" where D is a real number less than or equal to 1.0. For example, make K1/K3 = 10 by entering "=K1*0.1". Using "=K1" or "=K1*D" is particularly handy for parameter estimation, to limit the number of parameters being estimated.
- K3_vertical_bottom defines the vertical hydraulic conductivity of the lower half of the domain. This parameter is only used if there is vertical leakage with spatially-variable area source/sinks. To simplify inputs, you can enter "=K1" if K3=K1 and the domain is isotropic normal to the plane of the domain. If you want a fixed ratio of anisotropy, you can enter "=K1*D" where D is a real number less than or equal to 1.0. For example, make K1/K3 = 10 by entering "=K1*0.1". Using "=K1" or "=K1*D" is particularly handy for parameter estimation, to limit the number of parameters being estimated.

#### Details about Average_head

In other two-dimensional AEM programs such as TWODAN, the flow region is open to infinity, and one unknown that needs to be solved for is the amount of flow that goes between the modeled area and infinity. In these programs, to generate an equation to solve for that additional unknown, you specify a head at one location ("reference head" in TWODAN).

In Anaqsim, each domain model is closed and finite, so there is not that extra unknown. You specify the average head in each domain. which in turn defines a constant that is added to the potential for that domain. Since there are linesinks that bound each subdomain, the flow field outside those linesinks does not matter (the flow to/from infinity doesn't affect the solution inside the domain boundary). You could specify a variety of different average head values, within a reasonable range (close to the actual average), and get essentially identical results.

Say you have a simple Anaqsim one-domain model that has head-specified boundaries all around the external boundary, with h=100. There is zero recharge, so h should be 100 everywhere inside the domain. If you specify the average domain h=100, the program adds a constant to the potential that is the potential corresponding to h=100. On solving, it will turn out that the boundary conditions are met perfectly everywhere on the boundary and the boundary linesinks all have zero discharge; the analytical model will boil down to the simple equation h(x,y)=100. With zero discharges in the boundary linesinks, there is no flow to or from infinity to the model boundary from the outside (even though you never see this part of a domain model, it exists).

Now imagine that instead you set the average domain head to 110, which adds a larger constant to the potential for this domain. Now, to achieve the boundary h=100, the boundary linesinks need to extract water to pull the head surface down. In this case the solution on and inside the boundary will still be approximately h=100, but there will be flow to the outside of the boundary linesinks from infinity. Likewise, if you set the average domain head to 90, the solution on and inside the boundary will be approximately h=100, but there will be flow from the outside of the boundary linesinks to infinity. When you change the average head for a domain, it changes the part of the domain solution that you never see - the part that lies outside the external boundary of the domain.

If you use long line elements with few parameters and the average head is not close to the actual average, the differences in the external, unseen part of the model may have some visible impact on the model within the domain boundary. The most likely manifestation will be some lumpiness in the head surface near those boundary elements. Correct this by choosing a more representative average head and/or shortening line boundary elements and increasing the number of parameters per line.

#### Confined and/or Unconfined

Starting with release 2015-1, confined, unconfined, and confined/unconfined domains are in the same data table, which allows the user to quickly switch between these domain types. For confined and/or unconfined domains, these parameters are needed in addition to those that are common to all domains:

- Domain_Type determines which type of domain is to be modeled. Confined, Unconfined, and Confined/Unconfined are the three options. For Confined, the domain is always a fixed saturated thickness equal to the top elevation minus the bottom elevation, and the domain's transmissivity is independent of head. Confined domains generate linear equations, while the other two options can generate nonlinear equations. It is often wise to begin your modeling with confined domains which tend to converge faster and be less prone to numerical issues associated with nonlinearity and drying up. Later, it is easy to switch to unconfined or confined/unconfined domain types. With the Unconfined domain type, the domain is always unconfined and the Top_elevation is not used (although some value must be input). The Confined/Unconfined domain type behaves like a confined aquifer when the head equals or exceeds the top elevation, but like an unconfined aquifer where the saturated thickness depends on head, when head drops below the top elevation (see Strack, 1989, section 8; Haitjema, 1995, section 3.1.3; Strack, 2003, equation 3).

Never put an unconfined domain beneath an overlying domain, because the unconfined domain saturated thickness is always computed as head minus base elevation. If you think an underlying domain may become unconfined, use confined/unconfined rather than unconfined

- Top_elevation defines the elevation of the top of the domain. Not used for unconfined domain type.
- Bottom_elevation defines the elevation of the bottom of the domain.
- Storativity (S) is the dimensionless elastic storage parameter. S normally is the saturated thickness times specific storage (Ss). Starting with release 2021-1, there is an option to input the value of specific storage instead of storativity. This can simplify input by avoiding the step of multiplying by saturated thickness. To input specific storage, enter "Ss:" in front of the value. For example, if you entered "Ss: 0.0037", the value of storativity in this confined domain would be set to 0.0037 * b, where b is the saturated thickness (top elevation - bottom elevation). See Storage Parameter Details for more on how different storage parameters apply for different domain types.
- Specific_yield (Sy) is the dimensionless storage parameter for the unconfined domain type. See Storage Parameter Details for more on how these parameters apply for different domain types.

For numerical stability where the saturated thickness of an unconfined or confined/unconfined domain approaches zero, Anaqsim imposes a minimum saturated thickness. When heads drop near or below the bottom, the domain reverts to a confined-type domain with a fixed minimum saturated thickness. This facet of Anaqsim is governed by a parameter called Almost_dry_fraction under Solution/Solve Settings.

#### Confined Interface

For confined interface domains, these parameters are needed in addition to those that are common to all domains:

- Top_elevation defines the elevation of the top of the domain.
- Bottom_elevation defines the elevation of the bottom of the domain. The saturated thickness is the difference between the top and bottom elevations.
- Storativity is the dimensionless storage parameter that applies in the confined portion of this type of aquifer, equal to saturated thickness times specific storage. As described under Confined and/or Unconfined domains, you may enter specific storage values here instead of storativity values.
- Salt_elevation defines the elevation of the surface of the salt water, which is assumed to be static. Typically this is about the elevation of sea level.
- DensityRatio is the ratio of the salt water density to fresh water density. This varies from place to place, but is often near 1.025.

Interface domains in Anaqsim are based on the Ghyben-Herzberg approximation:

- The salt water is hydrostatic - pressure is proportional to depth below Salt_elevation. This assumption is reasonable when the flow is roughly steady.
- The fresh/salt water interface is sharp, with no mixing.
- Fresh water within a domain is hydrostatic (Dupuit approximation); there is no vertical resistance to flow.

The confined interface domains are based on the techniques presented by Strack (1989) on pages 101-106 and in Fitts et al (2015). These domains are confined with fresh water from top to bottom when heads are high enough that there is no interface, or they are confined with an interface and some salt water if heads are low enough. Confined interface domains would go to zero fresh water saturated thickness when the fresh water head drops to a level where the fresh water pressure at the top of the domain equals the salt water pressure at that elevation. This occurs where the freshwater head = Top_elevation + (Salt_elevation - Top_elevation) * DensityRatio. For numerical stability where the fresh water saturated thickness approaches zero, Anaqsim imposes a minimum saturated thickness. When heads are low enough, the domain reverts to a confined-type domain with this minimum saturated thickness. This facet of Anaqsim is governed by a parameter called Almost_dry_fraction under Solution/Solve Settings.

Generally transient simulations with interface domains will be inaccurate because it is assumed that the salt water has a hydrostatic distribution of pressure on the interface. In most transient situations, the salt water is moving, and when that movement has a vertical component, the hydrostatic pressure assumption is violated. If you feel the hydrostatic salt water assumption is still reasonable, you may proceed with a transient simulation but Anaqsim will issue a warning. See Storage Parameter Details for more on how storage parameters apply to this domain type.

#### Unconfined Interface

For unconfined interface domains, these parameters are needed in addition to those that are common to all domains:

- Bottom_elevation defines the elevation of the bottom of the domain. The saturated thickness is the difference between the head and the bottom elevation.
- Specific_yield is the dimensionless storage parameter for an unconfined aquifer.
- Salt_elevation defines the elevation of the surface of the salt water, which is assumed to be static. Typically this is about the elevation of sea level.
- DensityRatio is the ratio of the salt water density to fresh water density. This varies from place to place, but is often near 1.025.

Interface domains in Anaqsim are based on the Ghyben-Herzberg approximation:

- The salt water is hydrostatic - pressure is proportional to depth below Salt_elevation.
- The fresh/salt water interface is sharp, with no mixing.
- Fresh water within a domain is hydrostatic (Dupuit approximation); there is no vertical resistance to flow.

The unconfined interface domains are based on the techniques presented by Strack (1989) on pages 108-111 and in Fitts et al (2015). These domains are unconfined with fresh water from the water table to the bottom when heads are high enough that there is no interface, or they are unconfined with an interface and some salt water if heads are low enough. Unconfined interface domains would go to zero fresh water saturated thickness when the fresh water head drops to Salt_elevation. To avoid dry conditions, keep all heads in the domain above Salt_elevation. For numerical stability where the fresh water saturated thickness approaches zero, Anaqsim imposes a minimum saturated thickness. When heads are low enough, the domain reverts to a confined-type domain with this minimum saturated thickness. This facet of Anaqsim is governed by a parameter called Almost_dry_fraction under Solution/Solve Settings.

Generally transient simulations with interface domains will be inaccurate because it is assumed that the salt water has a hydrostatic distribution of pressure on the interface. In most transient situations, the salt water is moving, and when that movement has a vertical component, the hydrostatic pressure assumption is violated. If you feel the hydrostatic salt water assumption is still reasonable, you may proceed with a transient simulation but Anaqsim will issue a warning. See Storage Parameter Details for more on how storage parameters apply to this domain type.

#### Storage Parameter Details

Storage parameters are defined differently for different domain types as explained below.

- In confined domains, the storage parameter is always S (storativity) regardless of head.
- In unconfined domains, the storage parameter Sy (specific yield) applies when h > the head at minimum saturated thickness (see Almost_dry_fraction under Solve Settings). When head is lower than the level at minimum saturated thickness, the storage parameter equals S * Almost_dry_fraction (usually a very small value). This scheme helps with convergence in cases where heads fall below the bottom of the domain, but makes storage changes small under these conditions.
- In confined/unconfined domains, storage is like a confined domain when h >= top elevation, and like an unconfined domain when h < top elevation.
- In confined interface domains, the storage parameter = S + (porosity / (DensityRatio - 1)) where an interface is present. Where the domain is confined with no interface present (higher heads), the storage parameter = S. Generally, storage contributed by interface shifts greatly exceed elastic storage; S << porosity / (DensityRatio - 1).
- In unconfined interface domains, the storage parameter = Sy + (porosity / (DensityRatio - 1)) where an interface is present. Where head is high enough that there is no interface (inland from the toe of the interface), the storage parameter = Sy. Generally the storage contributed by interface shifts is larger than water table storage; Sy < porosity / (DensityRatio - 1).

In all but confined domains, the storage parameter is a function of head. The head at the start of a time step is used to determine the storage parameter that applies for the time step, even though the head at the end of the time step may correspond to a different storage parameter. This approximation helps convergence and is minor if time steps are small enough.

### Pumping Wells

Pumping wells may be either discharge-specified or head-specified. The discharge-specified type may be screened in one domain or across multiple domains if the well screen spans multiple model levels. All well input data tables may be accessed through the main menu or by using a pop-up context menu when the cursor is over the plot.

#### Input Common to all Pumping Wells

With all types of pumping wells, the following parameters are required.

- Label is a text label that helps you keep track of multiple wells.
- X,Y defines the horizontal coordinates of the well. To graphically edit a well's coordinates, left click to select it. Once selected, the well will be enclosed in a purple square box as shown below

When highlighted in this way, the well can be moved by clicking on the purple box and dragging it. This will automatically alter the coordinates of the well in the model input table. To stop graphic editing, press Esc to de-select the well.

- Radius defines the radius of the well. If there is a high conductivity filter pack around the well screen, the radius should be the radius of the borehole, not the radius of the screen.

#### Discharge-Specified

With discharge-specified wells, negative rates are used for extraction from the aquifer and positive rates are used for injection into the aquifer. Additional parameters defined here are:

- Domain defines the domain that the well screen is in.
- Discharge is the discharge of the well in units of [L
^{3}/T]. In transient simulations, this parameter may vary from one time period to the next.

#### Discharge-Specified (Multi-Domain)

Use this type of well to simulate a well with a screen that spans multiple domains and levels in the vertical direction. Anaqsim computes the appropriate discharge from each domain spanned so that the total discharge equals the specified discharge, and the heads at the well radius in each domain match each other. With discharge-specified wells, negative rates are used for extraction from the aquifer and positive rates are used for injection into the aquifer.

- Domains defines the domains that the well screen spans.
- Discharge is the discharge of the well in units of [L
^{3}/T]. In transient simulations, this parameter may vary from one time period to the next.

#### Head-Specified

With head-specified wells, you specify a head that applies at the well radius and Anaqsim computes the discharge needed to achieve that head. The discharge of a head-specified well may be checked after solving from the Analysis menu.

- Domain defines the domain that the well screen is in.
- Head_at_well defines the head at the well radius. In transient simulations, this parameter may vary from one time period to the next.
- Off_Periods is used in transient simulations if you want the well discharge to be zero during certain periods. The periods you want the well off are delimited with comma(s). For example, the following input is for a transient model with five time periods. The well pumps at rates such that at the end of period 1 the head at the well is 80, the well is off during periods 2 and 3, and the well pumps so that the head at the well is 115 and 118, respectively, at the ends of periods 4 and 5.

### Line Boundaries

A variety of line boundary conditions are available in Anaqsim. Each line boundary is a multi-segmented line (polyline) and the user inputs a list of vertexes in sequence from one end of the polyline to the other. The line boundary condition is approximated using linesink elements similar to those described by Jankovic and Barnes (1999).

Most line boundaries have a parameter that varies from one value at the starting vertex to another value at the ending vertex. The interpolation scheme between the end points is described in the next topic.

Anaqsim approximates the specified boundary conditions along line boundaries, as discussed by Fitts (2010). You may check the accuracy of line boundary condition approximations under the Analysis menu.

For internal line boundaries, the coordinates of all polyline points must not be outside the subdomain boundary, otherwise numerical havoc will be wreaked! It is possible for the start or end point of an internal line boundary to coincide exactly with a corner point of an external line boundary.

All line boundary input data tables may be accessed through the main menu or by using a pop-up context menu when the cursor is over the plot.

#### Input Common to all Line Boundaries

The following input items are common to all of the line boundaries:

- Label is a text label that helps you keep track of multiple line boundaries. Please see River Line Boundary Discharges for special labeling of River line boundaries so that you may sum the discharges of groups of river boundaries.
- Each line segment has a line element with a certain number of unknown parameters that is defined in the Parameters_per_line column. If you choose 1 for this, the strength of the element is constant along the line, if you chose 2 for this, the strength of the element varies linearly along the line, if you choose 3 for this, the strength of the element varies parabolically along the line, etc. Using more parameters per line can increase the accuracy of the boundary condition approximation along the line, but at the cost of additional equations in the system of equations that is solved. You only need a high number of parameters per line if the heads or discharges along the element are expected to vary complexly. In many cases, 3 or fewer parameters per line is plenty. Experiment with this and reduce the number of parameters to the minimum that gives you reasonable boundary condition accuracy, which you can check using right-click / Check Line Boundary Conditions. In the case of a Discharge-Specified line boundary, the user does not define Parameters_per_line because it is always set to 1; the discharge/length of these line elements is constant and equal to the specified discharge of the line boundary divided by its total length.
- Coordinates defines the x,y coordinates of the polyline vertexes. When you click on a cell in this column, a small text box window appears. In this text box, enter/edit lines of coordinate data, one line per vertex. Each line should have x and y values separated by either a comma or a tab character. You may digitize the coordinates of a polyline in the plot view and then paste the coordinates into this text box. In the case of normal flux-specified external boundaries, the coordinates must be listed in counter-clockwise order with the domain to the left of the boundary as you proceed along it. Once input, coordinates can be edited graphically by selecting the line boundary and then moving the vertexes or inserting or deleting vertexes. To graphically edit a line boundary's coordinates, left click to select it. Once selected, the vertexes will be enclosed in purple square boxes as shown below

When highlighted in this way, the vertexes can be moved by clicking on the purple box and dragging it. This will automatically alter the coordinates of the line boundary in the model input table. To de-select a line boundary, press Esc. Procedures for graphic editing of line boundaries are covered in the tutorial videos on the website. To reverse the order of the vertexes (this can be handy if you digitized in the wrong direction), click the Reverse button.

Most line boundaries have additional parameters defined at the start and end points, such as heads for head-specified line boundaries. With all such parameters, the same algorithm is employed to interpolate the specified values along the line boundary:

- Apportion the intermediate values to the vertexes based on the number of line segments in the polyline. For example, if there are 4 segments and the end values are 100 and 110, the values at the vertexes would be 100, 102.5, 105, 107.5, and 110. If there are 5 segments, the values at the vertexes would be 100, 102, 104, 106, 108, and 110.
- Linearly interpolate values to the control points within a line segment, assuming a linear distribution from one end to the other.

The additional parameters specific to each type of line boundary are listed in the following topics.

#### Head-Specified

With head-specified line boundaries you specify these additional items:

- Domain defines the domain that the line boundary is in. Double-click on this cell to open a drop-down list and make a selection.
- Check the Domain_Boundary check box if this line boundary is on the external boundary of the domain.
- The h_start and h_end are the specified head values at the starting and ending vertexes. In transient simulations, these parameters may vary from one time period to the next.
- Off_Periods is used in transient simulations if you want the line boundary discharge to be zero during certain periods. This could be handy for a dewatering trench that is turned off and on. The periods you want off are delimited with comma(s). For example, the following input is for a transient model with three time periods. The line boundary discharges at rates such that at the end of period 1 the heads along it are 104, the line boundary is zero during periods 2, and the line boundary discharges at rates such that at the end of period 3 the heads along it are 106.

Boundary condition equations are written at each control point on each line segment. The equation specifies that the modeled head = specified head (interpolated between endpoints of the line boundary). The specified head condition is approximated between control points and may be checked graphically.

#### Normal Flux-Specified

With normal flux-specified line boundaries you specify these additional items:

- Domain defines the domain that the line boundary is in. Double-click on this cell to open a drop-down list and make a selection.
- Check the Domain_Boundary check box if this line boundary is on the external boundary of the domain.

When the boundary is external, it is necessary to list the coordinates of the polyline in counter-clockwise order with the domain to the left of the boundary as you proceed along it. If coordinates are specified in the wrong (clockwise) order, there is a check in the program that will detect and report this error. This error can occur if you mistakenly specified the coordinates in clockwise order around the outside of the domain. This message can also result if you have made errors in specifying additional external line boundaries either at these same coordinates (e.g. have two external line boundaries in the input that share the same vertex coordinates), or have erroneous external line boundaries to the right of this one (see discussion of left/right algorithm).

- The Normal_flux_start and Normal_flux _end are the specified normal flux values at the starting and ending vertexes. Normal flux is the component of domain discharge normal to the line segment, and has units of discharge per length [L
^{3}/T/L] = [L^{2}/T]. This may also be thought of as the normal component of specific discharge times saturated thickness. The normal flux is positive for flow across the boundary from left to right as you proceed from the start toward the end. Normal flux is negative for flow from right to left as you proceed from the start toward the end. In transient simulations, the normal flux parameters may vary from one time period to the next.

Boundary condition equations are written for sub-intervals in each line segment (e.g. three subintervals for 3 parameters/line). The equation specifies that the total discharge across the line over the subinterval equals the total to interpolated specified fluxes over the interval. The accuracy of this approximation may be checked graphically.

#### Head-Dependent Normal Flux (3rd type)

This allows head-dependent normal flux into or out of the model, depending on the modeled head at the boundary. Boundary conditions like this are sometimes called 3rd type, Robin, general head (GHB in MODFLOW), or mixed boundary conditions; they involve both head and flux. The boundary must be an external boundary of a domain. It is necessary to list the coordinates of the polyline in counter-clockwise order with the domain to the left of the boundary as you proceed along it. If coordinates are specified in the wrong order, there is a check in the program that will detect and report this error. This error can occur if you mistakenly specified the coordinates in clockwise order around the outside of the domain. This message can also result if you have made errors in specifying additional external line boundaries either at these same coordinates (e.g. have two external line boundaries in the input that share the same vertex coordinates), or have erroneous external line boundaries to the right of this one (see discussion of left/right algorithm).

This kind of boundary is illustrated conceptually in the vertical profile sketched below. It is as though there is a fictional domain beyond the boundary (orange) outside of the domain (blue) , and at a distance b outside the boundary, there is a fixed head h*. The head difference between h* and the head at the boundary (h) drives a component of discharge normal to the boundary.