Jump to content
Geochemist's Workbench Support Forum

How to model temperature change and irrigation


Recommended Posts

My research work is related to deposition and precipitation of salts in soil under seawater irrigation.

Since I am considering soils in Abu Dhabi, the temperature becomes very important, as weather conditions here are extreme.

The models in GWB ( Virial ) can only model temperatures at 25 C.

I would like to compare results obtained using the Virial and the Bdot equations at different temperatures.

Is there anyway of comparing the models at different temperatures?

If not, how can I manipulate the temperature term.

Secondly, how can I use the platform to model precipitation or irrigation, so as to account for leaching of salts below the root zone.

Link to comment
Share on other sites

Hi Jubilee,

 

For applying the virial equations outside 25 C, try taking a look at thermo_phrqpitz.dat, which is included with the GWB. You should be aware that the data for temperature extrapolation is incomplete, however. You should see The PHRQPITZ/ PHREEQC documentation for more information on applying the virial methods outside 25 C. You might also take a look at the model of Yuan and Todd (1991), which is mentioned in Chapter 30 of the Geochemical and Biogeochemical Reaction Modeling textbook.

 

There are many chapters in the GBRM textbook which describe mineral precipitation. Once you familiarize yourself with the basics of chemical equilibrium and mass transfer, I would take a look at Chapter 24, Evaporation, and Chapter 27, Weathering.

 

Hope this helps,

 

Brian Farrell

Aqueous Solutions LLC

Link to comment
Share on other sites

Thanks Brian ,

I will work on this.

 

Please, why do i get this error message when i run my model in X1t?

 

Could not open X1t_plot.xtp file.

(Do you have write permission in this directory?)

 

I have been running the react model before now, but i just got access to the professional version of the GWB.

Thanks

Link to comment
Share on other sites

Dear Brian,

I am trying to model saturated groundwater flow in saline condition using X1t,

I keep getting the error message

"mineral volume in node zero is greater than node volume".

What am I doing wrong.

Also,is there any chapter of the GBRM text that treats examples of ground water flow using X1t

Finally, can Xplot be configured such that different minerals on the plot can be colour coded.

Regards,

Jubilee

Link to comment
Share on other sites

Hi Jubilee,

 

To figure the volume of a nodal block in X1t, just go to the Domain pane and multiply the node length (delta x) by the width (y) and the height (z). If you didn't set delta x, you can calculate it by dividing the total domain length by Nx. If you're specifying the volume of minerals in your system, make sure the sum of all minerals in a single nodal block does not exceed this value. If you still have a problem please post your .x1t file.

 

For transport in flowing groundwater, you should first look at Chapters 20 and 21 in the GBRM text.

 

Regards,

Brian

Link to comment
Share on other sites

Hi Brian,

I have a cylindrical laboratory soil column setup which i plan to operate under saturated conditions.The inlet fluid is concentrated salt water and the amendments in the soil are calcite and gypsum. The soil is mostly Quartz. I would like to use GWB to model the Precipitation of salt and also determine the change in concentration of ionic species at different levels in the soil column.

I have set up the process in X1t, but i keep getting errors. Please find below my sequence of commands. I must be doing something wrong.

Best regards,

Jubilee.

 

length = 60 cm

porosity = 40%

deltax = eqn "IF Inode < 4 THEN RETURN 5 ELSE return 10" cm

discharge = 100 pore_volume

time end = 1 yr

 

 

Scope = initial

Swap Gypsum for Ca++

Swap Calcite for HCO3-

swap Quartz for SiO2(aq)

react 20 volume% Gypsum

React 10 volume% Calcite

react 50 volume% Quartz

Ca++ = 46 meq

Mg++ = 144 meq

Na+ = 1129 meq

K+ = 20 meq

HCO3- = 3.0 meq

SO4-- = 1 meq

Cl- = 500 meq

pH = 8.0

 

Scope inlet

Ca++ = 9.4 mmol/kg

Mg++ = 52 mmol/kg

Na+ = 462 mmol/kg

K+ = 9.4 mmol/kg

SO4-- = 23 mmol/kg

Cl- = 521 mmol/kg

HCO3- = 3.0 mmol/kg

SiO2(aq) = 0.0001 mmol/kg

pH = 8.25

Link to comment
Share on other sites

Hi Jubilee,

 

There is a problem with the way you set your Initial system. You swap 3 minerals into the Basis (Gypsum, Calcite, and Quartz), but then you unswap the Gypsum and Calcite by specifying the concentration of Ca++ and HCO3-. You can specify either equilibrium with a mineral for a particular component, or the concentration of a dissolved species, but not both.

 

Later you use the commands "react 20 volume% Gypsum/ Calcite/ Quartz". If these 3 minerals are supposed to be in equilibrium with your system, you don't want to use the "react" command. Instead, you would say "Gypsum = 20 volume%". The "react" command is used to modify the Initial system by adding chemical mass to it. Commonly people use the react command to set up kinetic rate laws for non-equilibrium minerals.

 

I would recommend using units such as meq/kg or mol/l instead of just meq. That way the extent of the system doesn't matter when specifying a fluid's chemistry. If you are having trouble getting your model to even start, I would start by creating models in React for the Initial system and the Inlet fluid.

 

Hope this helps,

 

Brian

Link to comment
Share on other sites

  • 3 months later...

Hi,

 

The answer depends on whether you're accounting for reaction kinetics or not. Without kinetics, mass transfer is calculated in terms of reaction progress, which varies from 0 to 1. With a simple reactant, for example, you can set a value for delxi, the length of the reaction step. By default, this value is 0.01, meaning the program will take 100 steps to titrate in the entire reactant mass. You can set linear or log values, and use this command in conjunction with dx_init. The commands are explained in the GWB Reference Manual.

 

If you've set a kinetic rate law for mineral precipitation, however, the mass of mineral that precipitates in a single time step is calculated from integrating the rate law. The program chooses suitable time steps small enough to ensure a stable solution. You can optionally set some values, like delQ or theta, but in most cases the default works pretty well.

 

Hope this helps,

Brian

Link to comment
Share on other sites

  • 1 month later...

Hi Brian,

Thanks for the response.

Like I told you sometime ago, I am running a laboratory soil column with gypsum amendment of mass percentages varying from 1-5%.

I have set up my model to indicate the dissolution rate constants for Quartz and Gypsum ( the two major minerals making up the soil column).

The length of the column is 60 cm and I intend to monitor the composition of the pore fluid at intervals namely 5, 10, 15, 20, 30, 40, 50, and 60.

However, there seems to be a default setting in GWB that makes the initial reference x position 2.5 cm such that i cannot see what goes on at my preferred x positions. I have tried all methods but it seems i am not able to change this default setting.

 

Find below the code of commands.

 

time start = 0 day, end = 8.403 day
delta_x = table {5 5 5 5 10 10 10 10} cm
width = 8 cm
discharge = 25 pore_volumes
porosity = .38

scope = initial
H2O = 1 free kg
Ca++ = 1570 mg/kg
Mg++ = 21.2 mg/kg
Na+ = 4 mg/kg
K+ = 2 mg/kg
SO4-- = 2700 mg/kg
Cl- = 100 mg/kg
balance on Cl-
Br- = .12 mg/kg
pH = 7
SiO2(aq) = 1e-5 mg/l

scope = inlet
H2O = 1 free kg
Ca++ = 164.11 mg/l
Mg++ = 1639.24 mg/l
Na+ = 13133.2 mg/l
K+ = 503.26 mg/l
SO4-- = 3475.48 mg/l
Cl- = 20070.4 mg/l
balance on Cl-
Br- = 47.18 mg/l
pH = 8.25
SiO2(aq) = 1e-5 mg/l
kinetic Quartz 61 volume% rate_con = 4.2e-18 surface = 1000
kinetic Gypsum 1 volume% rate_con = 4.7e-9 surface = 1000
printout on
go

 

Am I missing out something?

 

Secondly, do i have to input dispersion values into the model of there is a range of default dispersion values in GWB that can be applied in this case.

 

Thanks,

Jubilee

Link to comment
Share on other sites

Hi Jubilee,

 

When we discretize a domain, we don't know anything about how chemical or physical properties vary within a single nodal block. In fact, the only place we know this information is at the nodal point, which is the center of a nodal block. By refining a grid (dividing the domain into more nodal blocks) we can solve our chemical and mass transport equations at more locations and potentially determine a more accurate solution, though at the expense of more computing effort.

 

To discretize a domain in X1t or X2t, you specify 1) the number of nodes, and 2) either the total length of the domain or the length of individual nodal blocks. Because you set delta x in your model to {5 5 5 5 10 10 10 10} cm, the boundaries between nodal blocks are located at 5, 10, 15, 20, 30, 40, etc., but the nodal points (the center of the nodal blocks where the equations are actually solved) are 2.5, 7.5, 12.5, 17.5, 25, 35, etc. These are the Xposition values you see. If you want to solve the equations exactly at the points you've mentioned above (presumably where you have sampling points in your column), then you should set Nx to 13 and deltaX to {2.5 5 5 5 5 5 5 5 5 5 5 5 5} cm. The initial offset will allow you to center your nodal blocks in your desired locations.

 

An alternative approach is to set a reasonably fine grid but not worry about the exact location of the nodal points. I suppose it depends on what you're trying to get out of your model. You can of course overlay scatter data onto your results in Xtplot using the exact positions of the sampling points and independent from the nodal points.

 

As for the dispersivity values, you should enter a value appropriate to your column, or whatever type of domain you are modeling. If you leave the dispersivity blank, GWB assumes a value equal to 1/100 of the domain's length. In X2t, the transverse dispersivity (if left blank) is assumed to be 1/10 times the value used for the longitudinal dispersivity. You should look into sections 2.5, 3.3, and 4.6 in the GWB Reactive Transport Modeling Guide for more on this topic.

 

Hope this helps,

Brian

Link to comment
Share on other sites

Thanks Brian,

This was helpful.

 

If i want to build a model for sweater infiltrating a porous medium that is composed of Quartz, Gypsum and Bio-char in the first 10-20 cm , and then only Quartz in the next 20 - 40 cm. Is there a way to carry out this kind of model in GWB?

I have reviewed the literature on modeling using X1t, but it seems this kind of scenario was not addressed.

Thanks,

Jubilee

Link to comment
Share on other sites

Hi Jubilee,

 

The answer is a little complicated. If the minerals are in equilibrium with the fluid in your domain (that is, they're swapped into the Initial pane) then they must exist at every node within the domain at the beginning of your simulation. You can click the + button next to entries on the Initial pane to set heterogeneous values throughout the domain, however. In this way, you could set large masses for Quartz, Gypsum, and Bio-char close to the inlet, but smaller (nonzero) values for Gypsum and Bio-char towards the outlet.

 

An alternative is to not specify equilibrium with these minerals directly (by swapping the mineral phase into the Initial system) and to instead choose appropriate values for concentration or activity of your aqueous species such that the fluid would be in equilibrium with each mineral near the inlet, but only Quartz towards the outlet. Again, you click the + button next to each species to set heterogeneous values. This won't allow you to set the masses of the minerals, however, only whether they end up being saturated or undersaturated in the initial system. SpecE8 or React will be useful for determining whether a mineral is saturated or undersaturated in solutions of a given composition.

 

A final approach is to set Gypsum and Bio-char as non-equilibrium reactant minerals from the Reactants pane. You might set them up as simple reactants, or as kinetic reactants. Again, click the + button next to the mineral's mass on the Reactants pane to set values that vary in space. In this case, you are allowed to set zero values for Gypsum and Bio-char towards the outlet.

 

You should look into the Heterogeneity Appendix of the GWB Reactive Transport Modeling Guide for more on setting values that vary in space.

 

Regards,

Brian

Link to comment
Share on other sites

  • 4 months later...

Hi Brian,

I have tried to run the script below several times using the thermo_Phrqpitz. However a couple of errors are displayer when i initiate the run.

Find my script below.

Regards,

Jubilee

 

# X1t script, saved Wed May 14 2014 by jadeoye
> data = "C:\Program Files\Gwb\Gtdata\thermo.dat" verify
> conductivity = "C:\Program Files\Gwb\Gtdata\conductivity.dat"
Reading conductivity data from: C:\Program Files\Gwb\Gtdata\conductivity.dat
> time start = 0 day, end = 275 hr
> delta_x = table {2.5 5 5 5 5 5 5 5 5 5 5 5 5 } cm
> width = 8 cm
> height = 8 cm
> Nx = 13
> discharge = .000265 cm3/cm2/s
> porosity = .38
> dispersivity = 30 cm
> heat_source = 0 J/cm3/yr, temp_min = 25, temp_max = 200
> scope = initial
> H2O = 1 free kg
> Ca++ = 1657.3 mg/l
> Mg++ = 1838.3 mg/l
> Na+ = 14776.9 mg/l
> K+ = 649.7 mg/l
> SO4-- = 5972.32 mg/l
> Cl- = 22800 mg/l
> balance on Cl-
> Br- = 101.343 mg/l
> pH = 8
> SiO2(aq) = 1e-5 mg/l
> TDS = 42000
> density = 1
> scope = inlet
> H2O = 1 free kg
> Ca++ = 579.3 mg/l
> Mg++ = 2036 mg/l
> Na+ = 14317.6 mg/l
> K+ = 603.05 mg/l
> SO4-- = 3834.1 mg/l
> Cl- = 22500 mg/l
> balance on Cl-
> Br- = 95.2 mg/l
> pH = 8.25
> SiO2(aq) = 1e-5 mg/l
> TDS = 42000
> kinetic Quartz 97 wt% rate_con = 4.2e-18 surface = 1000
> kinetic Gypsum 3 wt% rate_con = 4.7e-9 surface = 3000
> printout on
X1t >
Regards,
Jubilee
  • 0
Link to comment
Share on other sites

Hi Jubilee,

 

The thermo_phrqpitz.dat dataset doesn't include any silica species, so you won't be able to incorporate Quartz dissolution/ precipitation kinetics in your model as long as you're using that dataset. Thermo datasets utilizing the "Pitzer equations" can be very useful at high ionic strength, but there is generally less data available for them.

 

You might also try using one of the one-node programs (like React) to set up and verify the mineralogy and porosity of your system before starting the reactive transport model. I find it's easiest and most flexible to set mineral abundances in terms of volume fraction or volume %.

 

Hope this helps,

Brian

Link to comment
Share on other sites

Thanks Brian for the response.

I have gone through all the thermo data sets in the GWB. I tried to run my model making use of each one of them.

I was only able to run the model successfully with "thermo.dat and thermo.comV8.R6t.dat.

Are these the only two datasets with silica species.

May I know if there are alternatives for running the model using the Pitzer activity coefficient model.

Also, it is not clearly stated if the thermo.dat is a debye huckle method of Viral method.

Kindly clarify.

Regards,

Jubilee.

Link to comment
Share on other sites

Hi Brian,

My i also ask if there is a way to combine two activity coefficient models in GWB.

So, i know that at high ionic strength, the thermo.dat will not give me an accurate result. However, it contains the aqueous silicon oxide.

On the other hand, PHRQPITZ gives fairly accurate results at high ionic strength, but does not contain aqueous silicon.

I need a model that contains aqueous silicon oxide and also gives fairly accurate results at high ionic strength. My inlet water concentration is 42,000 ppm

 

Also, i would like to find out if any gwb user has complained about unreproducible results in the past. By this, i mean that, i ran the same codes a couple of months ago and now i am getting results that are totally different from what i got back then.

 

Regards,

Jubilee

Link to comment
Share on other sites

Hi Jubilee,

 

Your input to X1t needs to be tailored for the specific thermo dataset that you're using. For example, thermo_minteq.tdat and thermo_phreeqc.tdat also contain silica species, but the basis species H4SiO4 is used instead of SiO2. Since the basis species have different molecular weights and you've specified concentration in mg/l, you'll either need to convert to units involving moles or use elemental equivalents (i.e. mg/l H4SiO4 as Si). You might find GSS useful in this regard.

 

You can find the activity model on the third line of the thermo dataset if you're viewing it in an editor like Notepad or if you're using the TEdit thermo data editor it's on the header pane next to "Activity model". Section 2.3 (Thermodynamic datasets) of the GWB Essentials Guide describes the thermo datasets distributed with the software. The phreeqc, minteq, and wateq4f activity models are also variations on the basis debye-huckel method. Datasets labeled h-m-w or pitzer use the virial methods.

 

You cannot combine data from debye-huckel and virial equation based datasets. This is not a limitation of GWB or other similar modeling programs, but a result of the differences between the two theories. If you need silica in your model, I guess I would go with thermo.tdat or one of the other debye-huckel based datasets.

 

Hope this helps,

Brian

Link to comment
Share on other sites

Hi Brian,

Thanks for the response. This was helpful.

Can you also address my final question,

 

"Also, i would like to find out if any GWB user has complained about unreproducible results in the past. i have battled to reproduce the same result i obtained from a set of codes last month.By this, i mean that, i ran the same codes last month, and now, i am getting results that are totally different from what i got back then."

 

Regards,

Jubilee

Link to comment
Share on other sites

Hi Jubilee,

 

We haven't heard anything like this. Our last maintenance release for GWB9 was over a month ago, so I don't know why you'd see a difference in your results. Are you sure you ran the same exact input files? Remember, the thermo dataset is a part of your input as well. If you're loading different datasets or deleting species that don't exist in a particular dataset then you can't expect to get the same exact results. Do you have an example to illustrate this problem?

 

Regards,

Brian

Link to comment
Share on other sites

  • 2 weeks later...

Hi Brian,
I made the changes you suggested and was able to run the model using the different thermo data sets.
However, the experimental values i have obtained in the lab are not well predicted by the model.
Also, i have less accurate results using the phriqpitz and HMW datasets, despite the fact that my feed water has very high ionic strength (TDS = 42,000 ppm).
kindly let me know if there is a forum whereby i can get some hands on help, regarding the system i am trying to model.
Find below my model, in case one or two errors that may affect the accuracy of the model stick out.
Also, i have written a brief note to describe the system i am trying to model.

Code

X1t script, saved Sun May 25 2014 by jadeoye
data = "C:\Program Files\Gwb\Gtdata\thermo.com.V8.R6+.dat" verify
time start = 0 hr, end = 275 hr
delta_x = table {2.5 5 5 5 5 5 5 5 5 5 5 5 5 } cm
width = 6.3 cm
height = 8 cm
Nx = 13
discharge = 25 pore_volumes
porosity = .38
dispersivity = .409 cm
heat_source = 0 J/cm3/yr, temp_min = 25, temp_max = 25
scope = initial
H2O = 1 free kg
Ca++ = 1657.27 mg/l
Na+ = 14776.9 mg/l
Mg++ = 1838.17 mg/l
K+ = 649.7 mg/l
SO4-- = 5972.32 mg/l
Cl- = 22500 mg/l
balance on Cl-
pH = 7.5
SiO2(aq) = 1e-5 mg/l
scope = inlet
H2O = 1 free kg
Ca++ = 579.3 mg/l
Na+ = 14317.6 mg/l
Mg++ = 2036 mg/l
K+ = 603.05 mg/l
SO4-- = 3884.1 mg/l
Cl- = 22500 mg/l
balance on Cl-
pH = 7.75
SiO2(aq) = 1e-5 mg/l
TDS = 42000
density = 1.025
kinetic Gypsum 3 volume% rate_con = 4.7e-9 surface = 3000
kinetic Quartz 59 volume% rate_con = 4.2e-18 surface = 1000

 

Description

I have a cylindrical soil column (length = 60 cm, radius = 4 cm). The top 20 cm of the column is filled with sand (97% mass fraction) and gypsum (3% mass fraction).The rest of the column is filled with plain sand. The porosity of the soil column was calculated to be 38%. The column is orientated upside down, so that the 20 cm zone amended with gypsum is at the base.

I am flushing the soil column with seawater, from the base, at a rate of 2 mL/min using a peristaltic pump, and maintaining a hydraulic head at the top of the column to ensure that the system is operated under saturated flow. The pore volume will be displaced 25 times.

I wish to determine the change in concentration of calcium at different points in the soil column over time. This will help me monitor the dissolution of gypsum in the soil.

See the attached file, showing the model and experimental setup. I have also included the graphs in which i compared the model results and experimental results for the three gypsum mass fractions i considered (1%, 2%, 3%).

 

Regards,

Jubilee

Soil Colum study file.pdf

Link to comment
Share on other sites

Hi Jubilee,

 

Since Gypsum is a kinetic reactant (not an equilibrium mineral), you can set a value of zero for its mass in the downstream nodes to reflect the fact that it was only added to the first 20 cm. If Gypsum was initially in equilibrium with the system (swapped into the Basis on the Initial pane), on the other hand, you could set a very small mass in the downstream nodes. And if the mineral is initially in equilibrium but you set up a kinetic rate law for its dissolution, you could also set its rate constant to 0 in the downstream nodes. Please refer to the Heterogeneity Appendix to the Reactive Transport Modeling Guide for details on setting spatially varying field variables.

 

You'll want to actually calculate exactly how much Quartz and Gypsum are in your system. Right now, your input script does not reflect the solid phase mass fractions described in your attachment. Calculating the initial system ("go initial") and checking the mass and volume of the minerals, fluid, etc. should help you verify that your input matches the experimental conditions.

 

It looks like the flowrate of 2 mL/min is slightly different than the specific discharge resulting from the number of Pore Volume replacements you specify, but not by much. You might try adjusting the rate constant for Gypsum dissolution to see if you can get better agreement between your experimental and simulation results. If you measured other analytes in the fluid, you may be able to determine whether other reactions are occurring the in the column experiment.

 

Hope this helps,

Brian

Link to comment
Share on other sites

  • 2 weeks later...

Hi Brian,

 

Thanks for your respo

nse.

(1) can you expanciate on your previous comment (below).

 

You'll want to actually calculate exactly how much Quartz and Gypsum are in your system. Right now, your input script does not reflect the solid phase mass fractions described in your attachment. Calculating the initial system ("go initial") and checking the mass and volume of the minerals, fluid, etc. should help you verify that your input matches the experimental conditions.

 

How do i set the solid phase mass fractions. My current script is based on script from worked examples in the GWB transport modeling guide (rainwater infiltering a quartz aquifer). please share an example.

 

(2) How do i account for changes in mineral surface area (As) as the reaction proceeds. I did not consider this in my previous script and i believe it will be critical to the accuracy of the model. Please share an example if possible.

 

(3) is it possible to model evaporation in x1t. kindly share an example.

 

Regards,

Jubilee

Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...
×
×
  • Create New...