Frequently Asked Questions

1. In the Clog.file the time step is repeated and a truncation error is stated. What is possibly causing the truncation error?

   The message you encountered is related to the automatic time-step algorithm. See page 13 of the HST3D 2.0 documentation and page 115 of the HST3D 1.0 documentation. The program calculates the length of the next time step based on how much the pressure, temperature and/or concentration changed in the previous time step. It compares these with the "Max change in pres. (DPTAS)", "Max change in temp. (DTTAS)" and "Max change in mass frac. (DCTAS)" to calculate the length of the next time step. After it has advanced to the end of the next time step it checks how much pressure, temperature and concentration changed in that time step. If they changed too much, it will go back to the beginning of the previous time step and try again using a shorter step length. It will do this up to 5 times. The "Number of repeats of time step to achieve truncation error" is the number of times it had to repeat a time step to achieve a time step length for a suitable dependent variable change. When the program repeats a time step, it prints a message to the Clog file. Thus, the message you saw does _not_ mean there is anything wrong with your model. It is for your information.

2. Can HST3D simulate pure heat flux at an internal cell boundary?

   A flux boundary condition can be located on any cell face in principle. However, this feature has not been extensively tested.

3. How can I get a restart (Rst.xxx) file to be created?

   Use Data.form line 3.9.3

4. Is there a GUI for HST input processing and output display?

   Yes, there is an HST-PIE (plug in extension) for use with ARGUS-ONE software. It is available at: HST3D GUI.

5. Why do I get strange null values for concentration when using the direct solver with an excluded region?

   There is a bug in the reordering algorithm. Use the iterative solver, SLMETH=3 or 5, as a workaround until a fix is made.

6. In the HST3D documentation it is stated that the default boundary condition is no advective flux through the boundary faces of a cell. In a one layer model there exist vertical velocity components. Does there exist flux through the upper and lower horizontal boundary of the model?

   If you have a specified flux boundary on the upper or lower surface of a one layer model, there will be a flux through the upper or lower boundary of the model. You could also have an evapotranspiration flux through the upper or lower surface of a one layer model. If you don't have any boundary conditions on the upper or lower surface of the model, there is no flux through the upper or lower surface of the model.

7. What is the difference between specified temperature b.c. and associated temperature at specified pressure b.c.?

   At specified temperature boundaries, the fluid will be set to the specified temperature at the boundary regardless of the amount or direction of fluid flux. For example, If there were a buried sealed container of radioactive material, the heat of the radioactive material might cause the fluid to be heated near the container. However, there wouldn't be any fluid flow associated with that source of heat. At a specified pressure boundary, fluid may either leave or enter the model region. If the fluid enters the model region, the temperature of the fluid is determined by the "associated temperature at specified pressure". If the fluid leaves the model region at the specified pressure boundary, the "associated temperature at specified pressure" is ignored.

8. I want to simulate pumping of water from one interval of cased or uncased borehole (say the bottom) and recharge the used water back to the same borehole at another interval (say near the top). The pumped water is to be used as heating/cooling water, so it will be returned to the borehole at a different temperature. In short, I want to be able to simulate mixing of water discharged and simultaneously recharged to the same borehole. I need to account for convective heat exchange within the fluid due to buoyancy and mixing in the borehole plus conductive heat exchange to the surrounding rock. Is HST3D capable of doing this?

   This situation is beyond the current design of HST3D. You might be able to trick the simulator into doing your well configuration by defining two wells, one for injection and one for production and writing code to do the re-injection of the produced water.

9. When I compare the output from the examples to the provided output, I see differences. Are they significant?

   You should get agreement, except for the representation of zeros, for the simple examples with fixed time steps. The more complex examples with automatic time steps may not take exactly the same sequence of steps and may give results that only agree to two or three significant figures for some values. This is a consequence of different implementations of floating-point arithmetic on different operating systems and CPU chips. If significant differences occur, they will have to be investigated on a case-by-case basis, so please contact the author.

10. I noticed something missing from version 2.0.5 of HST3D. The beta release had a line 2.21.4A, which allowed one to define a geothermal-gradient initial condition. Now there seems to be no gradient initial condition option. Was that a mistake or was it left out deliberately?

    The ability to specify a geothermal-gradient initial condition is now included under the general options for data input by x,y,z range. There is a linear interpolation option (number 5) which allows for all sorts of piecewise linear initial conditions and boundary conditions.

11. I'm having a problem with restart. My data file includes lines 1.3 and 3.93. My Out.xxx file indicates that I have generated restart planes, but when I try to restart the file I get the following message: FORTRAN I/O ERROR 920: OPEN OF FILE CONNECTED TO DIFFERENT UNIT ATTEMPTED FILE: RST.casec.7 UNIT:9 LAST FORMAT: (A) My data file is casec.8.dat My output files are identified by casec.8 At the prompt for restart, I enter casec.8.

    The input restart file for the continuation run needs to be the name generated from the ancestor run, Rst.casec.7, in your file naming scheme. The input file to initiate a restart simulation has to be a new file that includes only Read 1 and Read 3 data. Using your original input file by simply changing lines 1.3 and 3.93 will not work.

12. I have set up a steady state model, that creates a restart file (Rst.out). This restart file is then used as input together with a new input file to initiate pumping. The following occurs: " run-time error F6508: READ(Rst.out) -too many bytes read from unformatted record". I have attached the help for the error F6508 and for the Keyword RECL.

    Ensure that the dimensions used to compile HST3D for the ancestor run that is being restarted are the same as those used to compile the executable that you are currently using.

13. How is the direction of gravity defined with respect to the coordinate system, which can be chosen freely? Is it implied that gravity is always positive with respect to the z-axis?

    Although the z-axis need not be parallel to the direction of gravity, the layers in HST3D should be arranged so that the index of node layers increases from the bottom to the top of the simulation mesh. Thus the component of the gravity vector along the z axis must point in the negative z direction.

14. Why can't the HST simulator re-wet a column of cells after the entire column has gone dry?

    The present algorithm only allows for a dry cell to start refilling through its bottom boundary. Lateral flow into a cell only occurs in subsequent time steps after re-wetting. Thus dry columns of cells can not re-wet. A workaround is to lower the bottom boundary of the region so that the column never becomes completely dry.

15. Is it possible to simulate the effect of tides on seawater intrusion in a coastal aquifer?

    The missing feature is the ability to use a function of time to specify the sinusoidal pressures at a sea boundary. Presently, HST can only represent transient boundary conditions as piecewise constant values. One would have to approximate the tidal variation as a set of steps.

16. I have read that "A change in volume occurs when fluids containing different amounts of salt are mixed. For example, mixing 1L of freshwater (density 1.0 g/cm^3) with 1L of brine (density 1.2 g/cm^3) yields a total mass of 2.2 kg of salt water, but the volume is 2.075 L, not 2.000 L, which is about a 4-percent difference." How is the mixed volume of 2.075 L calculated?

    Basically, when you mix fluids of different density and solute concentration (expressed as a mass fraction), mass is conserved but volume of solution is not. The steps are: 1) Take fluid 1 with density 1 and solute mass fraction 1 and fluid 2 with density 2 and solute mass fraction 2; 2) Mix a volume of fluid 1 with a volume of fluid 2 by first converting each volume to the equivalent mass of solution; 3) Then calculate the mass fraction of solute of the mixture using linear interpolation between the solute mass fractions of the two solutions interpolating by the relative mass percentage of each solution used; 4) Then use a formula or table to find the density of the mixture from the solute mass fraction. This will show that, for saline solutions, you can not just volumetrically add the two volumes of fluid that you mixed and arrive at the correct density. This is due to non-ideal behavior of mixing on a volumetric basis. For density differences of less than about 0.03 g/cm^3 or so things are pretty close to being additive and the error is small. For a wide density difference the more rigorous calculation is needed. That is why we use conservation of mass equations not conservation of fluid volume for variable-density flow.

17. On the seaward boundary of a coastal aquifer, where both specified pressures and specified mass-fractions exist, what is the difference between (1) specifying a mass-fraction value associated with a specified pressure boundary and (2) specifying the mass-fraction value by itself? Should both mass-fraction values be defined at the same boundary location?

    No. You must choose which type of boundary condition you want for a given part of a boundary. See the HST User's Manual Version 1, p.50. For case (1), the associated mass fraction is that which specifies the concentration of any incoming fluid. Outgoing fluid will have the concentration of the aquifer fluid at the boundary cell concerned. For case (2), the specified mass-fraction sets the concentration of the boundary cells concerned. This is not as physically realistic as case (1), especially where there is significant outflow of groundwater. You can not specify both an associated mass fraction and a specified mass fraction at the same boundary location.

18. What are the criteria for defining the HST3D cell layer thicknesses: 1) when the model represents a single hydro-stratigraphic unit; 2) when the model represents different hydro-stratigraphic units; 3) when the model represents a two-aquifer system separated by a semi-confining layer?

    It is important to discretize in the vertical direction sufficiently fine so as to allow the model to capture essential vertical variations in concentration or temperature. Buoyancy effects will be lost if too coarse a vertical discretization is used. Models that simulate ground-water flow only can often be discretized much more coarsely in the vertical direction. A sufficiently fine discretization of each aquifer unit must be done in order to represent density profiles for proper simulation of buoyant flow. Several nodes will be required in a semi-confining layer, especially near the boundaries where large contrasts in permeability are found.

19. I have a vertical-slice model region with boundary conditions that are uniform in the y-direction. However my computed velocity field is asymmetric in the y-direction. Why is this happening?

    A common cause of this problem is having the convergence tolerance for the linear equation iterative solver set too large.

20. I have a model where I've specified SI units and time units of seconds. I want to change the model to use time units of days. If I leave the units as SI and change the time units to days, I have to change the time values and flow rates. Should I also change the viscosity, reference pressure, etc. to use days instead of seconds?

    For SI units, viscosity stays in Pa-s always, pressure is Pa always. Output of velocity, mass and energy fluxes and flow rates will follow the simulation time units. Molecular diffusivity follows the units of time marching in versions since 2.09.

21. I get the error message:"142 - WBBAL - Wellbore excess residual flow for one or more wells". Could you briefly describe why this error might occur?

    This error message means that the sum of the layer flow rates does not equal the total rate for an injection or production well. This could mean a poorly designed problem if it occurs at the first time step. Or it could mean there is trouble converging to a solution between the flow equation and the well flow equation if allocation by mobility and pressure difference is being used. One possible cause is too large a time step.

22. I need to consider the effect of the geothermal gradient on the transport of heat. Would the heat conduction boundary condition be an appropriate way of doing this? Could the specified flux boundary be a more appropriate way of including the geothermal gradient? What happens when specified heat flux is used at an open flow boundary?

    The heat conduction b.c. does not provide a specified thermal gradient or specified heat flux. Instead, it responds to the temperature gradient at the boundary of the aquifer. So it is not appropriate for a geothermal flux boundary. The heat flux b.c. can vary in space. The heat flux b.c. is intended for boundary faces with no flow. It could be used with open boundaries but in most cases the heat transported advectively with ground water flow would overshadow the conductive heat flux. Along open boundaries, use the associated temperature to advect fluid with a given heat content into or out from the simulation region.

23. When I use centered (Crank-Nicolson) disretization in time, I get a fractional imbalance of 0.5 on the first time step. Why does this happen?

    The mass balance error of 0.5 on the first time step using centered-in-time discretization is a limitation of the simulator as written. The application of any change in the b.c. flow occurs at the end of the time step. This includes the change from no flow to some flow after time zero. With centered-in-time differencing, you have b.c. flow with the value from the beginning of the time step acting over the first half of the step and the value from the end of the time step acting over the second half of the step. Thus the amount of fluid from a flux b.c. that starts at zero is only that flow rate times half the first time step. Usually the first time step is small relative to the duration of the simulation and the error becomes insignificant.

24. When using the BCFLOW program, I get an error saying--Input statement requires too much data. What is wrong here?

    In your data-input file for BCFLOW, you must flag each type of boundary condition used in the simulation for which you are obtaining flow rates and cumulative flows. That is, at line 1.3, you must enter a 1 for each type type of boundary condition that was present in the simulation. However you are NOT required to define a patch for totalization for each type of boundary condition. This non-intuitive data requirement is a design flaw in the handling of data from HST3D to BCFLOW.

25. When using the linear interpolation method for specifying the values for a pressure boundary condition, the nodes at the upper end of the interpolation range do not get set to the desired specified pressure value. What is wrong here?

    The linear interpolation method does not extrapolate beyond the range of the two coordinate points used to define the straight line. Sometimes roundoff error will cause a coordinate point to be slightly less (or greater) than the end node that is to have a value assigned for the boundary condition. The remedy is to increase (or decrease) the coordinate point location slightly so that the range of the two coordinate points spans the range of nodes desired.

26. My model has three layers of different porous media. Although I defined them to be contiguous, I get two layers of zero conductance values in the z-direction. Thus, no flow takes place between the media layers. What is happening?

    You defined your grid spacing in the z-direction with a uniform distribution of nodal planes. Then you defined your zones such that the z-plane zone boundaries do not coincide with a z-plane of nodes (at least to 5 or 6 digits). Your zone boundaries lie between two planes of nodes. Because HST3D does not have a "snap-to-grid" feature (like many drawing programs have), your simulation region has two layers of zero z-conductance that prevents fluid and heat from moving into or out from the upper and lower zones. You can see these layers if you print out the O.kd file using 'prikd' print control variable. You need to change your data file to have the porous media zones start and stop at planes of nodes. You can either redefine the zones or adjust the mesh.

27. I have set up an HST3D model for groundwater heat storage using two wells. One well pumps warm water in and out of the aquifer and the other pumps cool water in and out of the same aquifer. During the first simulation period, one well injects warm water for 41 days while the other well is idle. The uniform initial condition temperature is 11.5°C and the temperature of the injected water is 17°C. After 41 days the temperature over much of the region has dropped below 11.5°C by as much as 0.30°C! What could be the cause of the observed energy-loss?

    This unexpected decrease in temperature is caused by the common assumption that the pressure-volume work term in the enthalpy transport equation is negligible for compressed water systems. However, the table of enthalpy as a function of pressure and temperature does contain the pressure-volume energy term. If we replace the enthalpy table with enthalpy as a linear function of only temperature, the simulation yields a slight temperature increase of up to 0.22°C above the initial condition temperature over much of the region.

email: klkipp@usgs.gov

[Home]
[Previous]