Domains¶

The properties of each domain (called a subdomain in Fitts, 2010) are set under Model Build → Domains. 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.

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 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.
K1 Angle to X Axis 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 (Top Half) 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 (Bottom Half) 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¶

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).

Warning

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.
Density Ratio 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) * Density Ratio. 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.
Density Ratio 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 / (Density Ratio - 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 / (Density Ratio - 1).
In unconfined interface domains, the storage parameter = Sy + (porosity / (Density Ratio - 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 / (Density Ratio - 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.