Posted by Rick Kossik
As part of the process of creating a very detailed online course for the GoldSim Contaminant Transport Module (for which a number of excerpts have been published in previous blog posts), I have been revisiting the most common questions and problems that we have been asked about by our users over the years in order to make sure to discuss these in the course (as well as to add more detailed descriptions, when warranted, in the User's Guide).In this blog post, I will first discuss the nature of the problem (when does it occur?) and how to address it in the current release of GoldSim (12.1). For those that are interested, I will then discuss the numerical reasons for this problem, as well as an alternative approach that will be available for addressing the problem in the next GoldSim release (expected later this year).
The Nature of the Problem
If the quantity of the medium (typically the volume of water) in a Cell is constant, 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., > 20% or so) 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. Moreover, 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).
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.
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 a particular algorithm implemented by GoldSim (the use of fractional timesteps). This algorithm 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 at the end of this article.
If the conditions that can lead to either of these two problems occur, warning messages 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) involved in an advective mass flux goes to zero.
Fortunately, when these problems do arise, their impact can be minimized as described below.
Addressing the Problem
The approach to addressing this problem is to prevent the volume from actually reaching zero, and reducing the errors produced as it approaches zero to an acceptable level by dynamically adjusting the timestep.
This can be done as follows:
Step 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). Such amounts should be modeled using a Pool (or Reservoir) element.
Step 2
Define a maximum timestep allowed for that Cell based on the Cell media amount, and the inflows and outflows of that media using the following logic:
Assuming that the volume of the Cell that we are concerned with is computed using a Pool element (named, for example, Volume1), we could express this logic in an Expression (e.g., MaxStep1) that defines the maximum allowable timestep as follows:
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 if the fractional change is too high. For most systems, the fractional change will likely need to be greater than 20% or so to produce a significant error.
Step 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. For example:
Step 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):
Dynamically adjusting the timestep in this manner is discussed in the Goldsim Help file here.
How Well Does This Work?
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 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).
Numerical Details: Why Do These Problems Occur?
For those of you who are interested in the numerical details of how GoldSim solves the transport equations, below I describe why the problems we discussed above occur. You do not need to understand these details to use the method discussed above. However, in this section, I will also discuss an alternative approach that will be available for addressing these problems in the next GoldSim release (expected later this year). While not suitable for all systems, this new approach will be suitable for some.
As stated above, there are actually two separate issues involved: 1) when the fractional change in the amount of media in a Cell is high (e.g., > 20%) over a timestep (this most commonly happens when the volume is very low and rapidly changing - either emptying or filling), computed concentrations in a Cell can show unrealistic spikes; and 2) if a Cell completely empties (the volume goes to zero), some mass will be left in the Cell after it has emptied.
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 a particular algorithm implemented by GoldSim (the use of fractional timesteps). This algorithm improves accuracy under most conditions, but can produce erroneous concentrations with rapidly changing media amounts.
Fractional Timesteps
Let's discuss the latter issue first. GoldSim employs a backwards difference (fully-implicit) algorithm to solve the Cell equations. This approach has the advantages of being unconditionally stable, conserving mass, and having bounded error levels. In order to enhance the accuracy of the backwards difference method, GoldSim automatically divides every full timestep (i.e., the timestep you specify in the Simulation Settings) into a number of fractional timesteps (up to 10, depending on the user’s selected precision level, as we will discuss below). A key point to understand is that these fractional timesteps are used only when carrying out the Cell computations, and are not applied to other GoldSim model elements outside of Cell nets (such as volumes or flow rates). These other variables stay constant over the full timestep and do not change during the fractional timesteps.
Using fractional timesteps is a very computationally efficient way to enhance the accuracy of the solution. This is because 1) only the Cell nets are updated during these fractional timesteps; and 2) numerically the fractional timesteps are computed via matrix calculations in such a way that it takes less computational effort than simply shortening the actual full timestep.
The mass transfer terms in the matrix equations solved by GoldSim are a function of the mass in the Cell, as well as a number of other parameters (e.g., volumes and flow rates). For example, GoldSim must compute a concentration for each solute being advected, and this requires the volume (and other parameters if sorption is occurring). For each fractional timestep, GoldSim updates the amount of mass in the Cell. To do so, at each fractional timestep GoldSim needs to compute concentrations (which in turn are used to compute mass transfer rates for that fractional timestep). However, as pointed out above, the values for the other parameters (most importantly, the media amounts such as the volume) are not updated every fractional timestep (we only have values for the current and previous full timestep). Moreover, the nature of the fractional timestepping algorithm is such that each fractional timestep within a full timestep must use the same value for the variables such as volume (we can’t, for example, interpolate between the two timesteps). However, using the value of the volume at the end of the full timestep to compute concentrations (used in turn to compute mass transfer rates) during a fractional timestep is not appropriate for all of the fractional timesteps (only for the last one). To address this, an approximation is used to handle systems in which the amount of media in the Cells is changing over a timestep. Under these circumstances, an effective value is used to represent the amount of each medium for the purpose of computing concentrations to compute mass transfer rates during the fractional timesteps. The effective value is a weighted average of the values at the current and previous timestep (it is weighted more heavily toward the current timestep).
The use of an effective value for the media amount reduces the error in the computed concentrations (and hence mass transfer rates) used in each fractional timestep, but does not eliminate it completely (i.e., it is better than using the value at the end of the timestep for all fractional timesteps, but it is still an approximation). If the volume is changing slowly, it provides a very accurate solution. However, when the media amounts do change significantly over a timestep, computed concentrations in a Cell can show jumps or drops (i.e., positive or negative “spikes”). For example, when the volume is dropping, not enough mass is moved out, so the concentration artificially increases, and when the volume is rising, too much mass is moved out, so the concentration artificially decreases. These can result in inaccurate mass transfer rates and can therefore propagate downstream producing additional inaccurate behavior. Such errors are a consequence of the fractional timestepping algorithm (and the requisite use of a weighted average for the media amounts).
The solution we discussed above minimizes any errors to an appropriate level by dynamically reducing the timestep to ensure that the media amounts do not change too rapidly over a timestep. Under these conditions, the fractional timestep algorithm produces accurate results.
An Alternative Solution to Rapidly Changing Volumes
If we choose to use fractional timesteps, we cannot avoid the error associated with rapidly changing media amounts (without dynamically shortening the timestep). However, if we do not use fractional timesteps, we can avoid the error, but in order to do so we also must no longer use a weighted average for the media amounts.
Within GoldSim there are three Solution Precision settings that you can select (from the Options dialog accessed via Model|Options): Low, Medium (the default) and High. These settings control a number of parameters associated with how GoldSim solves the contaminant transport equations. One of these is the number of fractional timesteps used. In particular, the value is 1 for Low (i.e., no fractional timesteps are used), 4 for Medium and 10 for High.
In the current version (12.1) as well as previous versions of GoldSim, a weighted average of media amounts is used for computing the mass transfer terms (as discussed above). While this is required for Medium and High precision, it is not, in fact necessary, for Low Precision (and the advantage gained by doing so is small). In the next release of GoldSim, a weighted average will no longer be used when running under Low Precision (instead, the current value for parameters like the volume is used, which is consistent with a fully implicit approach). Under these conditions (running with Low Precision), no spikes are produced due to rapidly changing volumes.
So based on this, you may think that the solution to avoiding concentration spikes (once the new version is available) is simply to run the model with Low Precision (as opposed to running with Medium or High Precision and using the dynamic timestepping approach outlined previously). For some models, that may indeed be an appropriate solution. However, just because solution does not produce spikes does not mean it is sufficiently accurate. As pointed out above, the fractional timestep algorithm was implemented in GoldSim because it provides an increase in accuracy in a very computational efficient way.
As a very simple illustration of this, consider a model in which we are computing the decay of 100 g of species with a half-life of 5 days in a Cell (with no inflows or outflows) for a period of 10 days, using a 1 day timestep. The solution is shown below:
As can be seen, there is a significant accuracy improvement as the Precision Setting increases. In this simple example, to achieve the same accuracy as the Medium Precision run, a Low Precision run would require 4 times more timesteps. But the key point is that the Medium Precision run is not as computationally intensive as a Low Precision run with 4 times the number timesteps (which are required to achieve the same accuracy). How much less computationally intensive? It depends on the nature of the model (e.g., how many other elements not associated with the Cell network exist in the model), but the advantage is likely to be at least 2 and could in some cases approach 4.
It is also important to understand that while in this simple decay example, the solution accuracy was only controlled by the number of fractional timesteps (1, 4 or 10), other variables controlled by the Solution Precision can play a role in more complex models (e.g., with different types of mass transport processes, more complex transient behavior, isotopes, solubility constraints, etc.).
Hence, this alternative approach using Low Precision (available in the next version of GoldSim ), while effective at avoiding spikes, can have a cost in terms of accuracy. For some models (particularly those with no solubility constraints, no rapidly decaying species and no isotopes) you may be able to simply run at Low Precision (in order to avoid spikes). However, you may then need to decrease the size of the timestep to accurately solve the mass transport equations (and this would likely be more computationally intensive than using Medium Precision and dynamic adjustment of the timestep discussed above). Certainly for complex mass transport models (e.g., with rapidly decaying species, solubility constraints and/or isotopes), however, using the Low Precision approach is not recommended.
Cells Going Empty and Leaving Mass Behind
The second problem we discussed was the fact that if a Cell completely empties (the media amount goes to zero), some mass will be left in the Cell after it has emptied. This issue of leaving mass behind in an empty Cell occurs whether or not fractional timesteps are used (i.e., it will also occur in the next version of GoldSim when Low Precision is specified).
The problem arises because when formulating the mass transfer terms for the matrix equations describing the Cell network, the implicit solution method by definition uses the mass (and the media amounts such as volume) at the end of the timestep (or, in the case of Medium and High Precision, a weighted amount). As a result, during the timestep at which the Cell goes empty, the mass transfer rate is either zero (for Low Precision, since we cannot compute a concentration if the media amount is equal to zero) or is too low (for Medium and High Precision, since we have used a weighted media amount that overestimates the media amount and therefore underestimates the concentration and hence the rate). Hence, during that final timestep before the Cell empties we either transfer no mass out at all (Low Precision) or too little mass (Medium or High Precision). That is, in all cases, some mass is left in the Cell (Low Precision leaving more mass than Medium Precision, and Medium more mass than High Precision).
The only solution to this problem is to set a non-zero Lower Bound. The amount of mass left behind can be minimized by setting the value appropriately low. It is also important to emphasize that this mass is not "lost". It still exists in the Cell, and in most models where a Cell goes dry, it very likely will subsequently fill again, so that the mass re-enters the system.
To keep this in perspective, 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 extremely complex, and hence it would very difficult (if not impossible) to model their concentrations accurately. This is because the processes controlling concentrations as a Cell empties are often difficult to accurately represent. In particular, in most cases, one of the key outflows for a Cell that is emptying (e.g., a pond) would be evaporation. In such a case, solutes would indeed precipitate out of solution before the Cell emptied (and therefore be left behind in the Cell). Moreover, during this period (and the period when the Cell was refilling), the geochemistry of this high concentration solution would be quite complex indeed.
No comments:
Post a Comment