Dealing with Changing Media Amounts in a Cell

If the quantity of the medium involved in advecting mass (typically the volume of water) in a Cell is constant or changing slowly, the solution algorithm is quite accurate. However, rapidly changing Cell media amounts can result in errors in computed concentrations (and mass transfer rates).

These errors are most common and significant when the amount of media in a Cell is low and rapidly changing (either emptying or filling), in which case the fractional volume change is by definition high. In this case, computed concentrations in a Cell can show unrealistic jumps or drops (i.e., positive or negative “spikes”). These can, in turn, result in inaccurate mass transfer rates and can therefore propagate downstream producing additional inaccurate behavior. The errors will be small as long as the amount of media does not change significantly (e.g., > 10 or 20%) over a full timestep, but can become significant if it does.

If a Cell completely empties (the media amount goes to zero), an additional problem occurs. In particular, when a Cell empties, some mass will be left in the Cell after it has emptied. More importantly, if the Cell is empty but still has an inflow (with a matching outflow), mass will accumulate in the Cell. If the Cell then refills, any mass left behind (or accumulated) in the Cell would be released suddenly (and this itself could result in unrealistic spikes).

Note: Of course, if the Cell emptied at least partially due to evaporation (which would not be unusual, for example, for a surface water body that went dry), some mass realistically would in fact be left behind in the Cell, since evaporation does not transport mass (and mass would eventually precipitate out of solution as water evaporated). The discussion above, however, refers to the fact that even if the Cell empties due to outflows that transport mass (and there is no evaporation), mass will be (unrealistically) left behind for numerical reasons.

To keep this problem in perspective, however, it should be noted that in the real-world, the behavior of solutes in a system in which a Cell was emptying would likely be very complex, and hence very difficult to model accurately. This is because the processes controlling concentrations as a Cell empties are often difficult to accurately represent. In particular, in many cases, as pointed out above, one of the key outflows for a Cell that is emptying (e.g., a pond) would be evaporation. In such a case, solutes would precipitate out of solution before the Cell emptied (and the geochemistry of such a high concentration solution would be quite complex indeed).

Note: Numerically, these two problems arise for two different reasons. The error associated with mass being left behind in an empty Cell is a direct outcome of using an implicit approach to solve the equations. The error associated with rapidly changing Cell volumes is a result of the fractional timestep algorithm implemented by GoldSim (which improves accuracy under most conditions, but can produce erroneous concentrations with rapidly changing media amounts). A detailed explanation of why these problems arise is provided in Appendix B of the Contaminant Transport Module User’s Guide (in the sections entitled Solving the Matrix Equations Numerically and Minimizing Errors for Systems with Rapidly Changing Media Amounts).These sections also provide a more detailed explanation of the two methods presented below for addressing these problems.

If the conditions that can lead to either of these two problems occur, a warning message will be written to the Run Log. In particular, GoldSim monitors the volume of the Reference Fluid and writes a warning message if this is changing in such a way that it could cause erroneous results. It also writes a (separate) warning if the amount of any media (typically the Reference Fluid) goes to zero.

Fortunately, when these problems do arise, their impact can be minimized. There are two different approaches for doing so.

One approach is to simply use Low Precision when running the model. As pointed out above, the error associated with rapidly changing Cell media amounts is a result of the fractional timestep algorithm implemented in GoldSim. This algorithm, however, is not used for Low Precision runs. As a result, the “spiky” behavior will not manifest itself. (In fact, if you run at Low Precision, GoldSim will not produce the warning message mentioned above resulting from rapidly changing Reference Fluid volumes).

However, the problem associated with leaving mass behind in an empty Cell still applies for Low Precision runs. This can be addressed by defining a (small) non-zero Lower Bound. Some mass will still be left behind in the Cell, but the error can be minimized by making the Lower bound suitably small.

Although this approach is certainly easy to implement, it comes with a potential cost. The fractional timesteps (used in Medium and High Precision runs) provide significant accuracy improvements (with minimal computational expense) during periods when the media amounts are not changing rapidly. Moreover, other variables controlled by the Solution Precision can play an important role in more complex models (e.g., with different types of mass transport processes, more complex transient behavior, isotopes, solubility constraints, etc.).

So for some models (particularly those with no solubility constraints, no rapidly decaying species and no isotopes) you may indeed be able to simply add a Lower Bound and run at Low Precision to avoid spikes. However, you may then need to decrease the size of the timestep to accurately solve the mass transport equations. For complex mass transport models (e.g., with rapidly decaying species, solubility constraints and/or isotopes), however, using Low Precision is not recommended, and the approach discussed below should be used.

In this approach, we take advantage of the accuracy improvements of the Medium and High Precision setting (when Cell media amounts are not changing rapidly) while also avoiding the inaccurate behavior (when they are changing rapidly).

The approach to doing so is to reduce the error to an acceptable level by dynamically adjusting the timestep only when these problems occur.

1. Create a non-zero Lower Bound for the media amounts that are rapidly changing in the Cell (typically the Reference Fluid water, but could be another Fluid, or a Solid if you were simulating advection in these media).

2. Define a maximum timestep allowed for that Cell based on the Cell media amount, and the inflows and outflows of that media.

Assuming that the volume of the Cell that we are concerned with is computed using a Pool element (named, for example, Volume1), we would define an Expression (e.g., MaxStep1) that defines the maximum allowable timestep as follows:

if(Abs(Net_Out) > 0 m3/day and Volume1.Total_Outflow* Default_DT / Volume1>Allowed_Fraction, Allowed_Fraction*Volume1/abs(Net_Out), Default_DT)

In this equation, Default_DT is the default scheduled timestep and Net_Out is the total outflow rate minus the total inflow rate:

The Allowed Fraction is the maximum allowable fractional change that we select to appropriately minimize the error. GoldSim provides a warning message (when running with Medium or High Precision) if the fractional change is greater than 10%. For many systems, the fractional change will likely need to be higher than this to produce a significant error, but the warning is conservatively displayed at 10%. Experimentation will be required to determine the appropriate value. You want to use the largest value that minimizes the error to an acceptable level (the smaller the value, the higher the computational expense).

3. If there are multiple Cells that may potentially be changing rapidly in this way, repeat step # 2 for each Cell and then compute the minimum of these values.

4. Enter this value into the field defining the Maximum time between updates in the Advanced Time Settings (accessed via the Time tab of the Simulation Settings dialog):

This algorithm works very well and is generally very efficient (the timestep is only minimized rarely and when it needs to be).

The two parameters that control the accuracy are the Allowed Fraction and the Lower Bound. The smaller these values, the smaller the error. Of course, as these values are reduced, the timestep becomes smaller (and the computational effort increases). As a result, some experimentation may be required to determine the most appropriate values for a particular system. As pointed out above, for most systems, the fractional change will likely need to be greater than 10% to 20% or so to produce a significant error. Of course, a larger Lower Bound also results in more mass left behind in an “empty” Cell, but in many cases, this may not pose a problem (e.g., it quickly re-enters the system when the Cell refills).

The higher the net outflow (Net_Out above) as the volume approaches zero, the greater number of timesteps that will be required (e.g., the smaller Allowed Fraction may need to be). It should be noted, however, that in the real-world, outflow rates generally decrease as the volume decreases (due to natural feedback loops). For example, as the volume of a pond decreases, the surface area of the pond decreases, and the evaporation rate decreases. Similarly, as a pond volume decreases, the pumping rate will often be lowered (e.g., either by turning off one of several pumps or adjusting the pumping rate of variable speed pumps). These kinds of feedback loops result in a more gradual approach to the Lower Bound (and it is strongly recommended that you represent them in your model if they are present).