Assessments of irrigation water quality

[Charlie]

Has anyone used PHREEQC for assessing water quality in terms of irrigation purposes?  Are there any recommended methods they would suggest, or special factors to consider? Thanks.

[Charlie]

I have been reading a paper by Wanda et al, (2013). (Physics and Chemistry of the Earth 66 (2013) 51–59).  This gives a helpful and explanatory example of an assessment of irrigation water using PHREEQC.

The purpose of this study was to evaluate and ascertain the suitability of the groundwater for irrigation purposes. A number of researchers have proposed different methods of analyzing irrigation water quality data. These metrics are regarded as the most effective ways to communicate irrigation water quality.  They include:

#Total hardness (TH) calculated using mol/l.  Be aware that conversion of non-specific species data (lab data) to ppm is not compound specific (assumes only CaCO3 would ppt)

#TDS content of the water is considered satisfactory when it contains <1000 mg/l, fair if it contains between 1000 and 2000 mg/l, and inferior when >2000 mg/l.

#RSC. The excess sum of carbonate and bicarbonate in groundwater over the sum of calcium and magnesium also influences the unsuitability for irrigation. RSC value <1.25 meq/l is good for irrigation, a value between 1.25 and 2.5 meq/l is of doubtful quality and a value >2.5 meq/l is unsuitable for irrigation. Land irrigated with such water becomes infertile owing to deposition of sodium carbonate and will affect the crop yields. 

#Permeability index (PI).  The soil permeability is affected by the long-term use of irrigated water rich in Na+, Ca2+, Mg2+, and HCO3- and the soil type. Doneen (1964) classified irrigation waters in three PI classes. Class-I and Class-II water types are suitable for irrigation with 75% or more of maximum permeability, while Class-III types of water, with 25% of maximum permeability, are unsuitable for irrigation.

#Sodium concentration is important in classifying irrigation water because sodium reacts with soil to reduce its permeability. Excellent to good irrigation water (% Na+ < 60%) will not cause soil aggregates to disperse and subsequently cause reduction in permeability and are suitable for irrigation.

#The Kelly’s ratio of unity or less than one is indicative of good quality of water for irrigation.  Values greater than one suggests unsuitability of water for agricultural purpose due to alkali hazards.

#Electrical conductivity is a good measure of salinity hazard to crops, as it reflects the TDS in groundwater. High salt content in irrigation water causes an increase in soil solution osmotic pressure.  The osmotic pressure is proportional to the salt content or salinity hazard. The salts, affecting the growth of plants directly, also affect the soil structure, permeability and aeration, indirectly affecting plant growth. Electrical conductivity is considered the most influential water quality guideline on crop productivity. At higher EC, less water is available to plants.  Water are classified as excellent (C1, EC range = 100–250 μS/cm), good (C2, EC range = 250–750 μS/ cm) and doubtful (C3, EC = 750–2,250 μS/cm).

#Sodium adsorption ratio (SAR). SAR values <10 are excellent (S1 class) and the water is suitable for any crop and can be used for irrigation in almost all types of soils.

#SC Specific conductance (μS/cm).

###Although all these indices were evaluated in this study, the SAR is probably the only one in current use and is generally considered an effective evaluation index for most water used in irrigated agriculture.


Use SOLUTION_SPREAD for large water sample dataset

#####################################

USER_PUNCH   #calculations use meq/l, except TH.  Convert mmol/L to meq/l = mmol/L X charge (abs)

-headings TH_(mol/L)   SAR   RSC   EC(uS/cm)   PI   KR   %Na+ 
-start
10 TH = mol("Ca+2") + mol("Mg+2")        # Make sure L of water used in SOLUTION.

20 SAR = (mol("Na+")*1e3*1)/((((mol("Ca+2")*1e3*2) + (mol('Mg+2")*1e3*2))/2)^0.5)   
 
30 RSC = ((mol("HCO3-")*1e3*-1) + (mol("CO3-2")*1e3*-2))-((mol ("Ca+2")*1e3*2) + (mol("Mg+2")*1e30*2))
   
40 PI = ((mol("Na+")*1e3*1)+((mol("HCO3-")*1e3*-1)^0.5))/((mol("Ca+2")*1e3*2) + (mol("Mg+2")*1000*2) + (mol("Na+")*1e3*1))*100   

50 KR = (mol("Na+")*1e3*1)/((mol ("Ca+2")*1e3*2) + (mol("Mg+2")*1e3*2)) 

60 %Na+ = (((mol ("Na+")*1e3*1) + (mol("K+")*1e3*1))/((mol ("Ca+2")*1e3*2)+(mol("Mg+2")*1e3*2)+(mol ("Na+")*1e3*1)+(mol("K+")*1e3*1)))*100     

70 PUNCH TH SAR  RSC  SC  PI  KR  %Na+
-end

######

Graphical outputs:
United States Salinity Laboratory Staff (USSLS) plot used to classify groundwater in terms of degree of suitability for irrigation: a plot of EC vs. SAR

Wilcox diagrams used to classify groundwater in terms of degree of suitability for irrigation: a plot of EC vs. %Na+

[davidjoferreira]

Dear Charlie,

Thank you for this valuable topic!

I want to calculate the SAR index in PhreeqC, but I am experiencing some issues...

I only know the basics of PhreeqC, but as far as I understand, I can paste the code lines into the "input", after the solution description, and before "END".

But, where does the SAR value comes up? Is it the only way to achieve this calculation?

I would be tremendously grateful if you could help me regarding this.

Thank you!

David

[dlparkhurst]

Here is a revised version of Charlie's script. It demonstrates USER_PRINT to write values to the output file, and SELECTED_OUTPUT/USER_PUNCH to write to an external file, in this case, "sar.sel".

USER_PRINT and SELECTED_OUTPUT/USER_PUNCH must be defined in the same simulation as SOLUTION (or SOLUTION_SPREAD), or in a previous simulation (simulations end with the line "END").

You should probably check that the formulas are correct. TOT is the Basic function for total concentration of an element in mol/kg water; TOT("water") is mass of water in kg; SOLN_VOL is the solution volume in liters. You should use one of these databases--phreeqc.dat, Amm.dat, and pitzer.dat--which have the information to calculate specific conductance and solution volume. Further, you need analyses with all of the major ions to calculate these quantities. Of course, you can still calculate the SAR, RSC, etc without having an accurate value for SC or solution volume, so other databases could be used, but SC will not be calculated and solution volume will be assumed to be 1.0.

***Corrected according to next post***

Code: [Select]

SOLUTION 1  SEAWATER FROM NORDSTROM AND OTHERS (1979)
        units   ppm
        pH      8.22
        density 1.023
        temp    25.0
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Si              4.28
        Cl              19353.0
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
USER_PRINT
-start
10 Ca_meql = TOT("Ca")*2*1e3*TOT("water")/SOLN_VOL
20 Mg_meql = TOT("Mg")*2*1e3*TOT("water")/SOLN_VOL
30 Na_meql = TOT("Na")*1e3*TOT("water")/SOLN_VOL
40 K_meql =  TOT("K")*1e3*TOT("water")/SOLN_VOL
60 Alk_meql = ALK*1e3*TOT("water")/SOLN_VOL

100 REM ==Total hardness mol/L ======================================
110 TH = (Ca_meql + Mg_meql)/2/1e3
120 REM ==Total hardness mg/L as CaCO3 ==============================
130 TH_CaCO3 = TH * 1000 * GFW("CaCO3")     
140 REM ==Sodium absorption ratio ===================================
150 IF (Alk_meql > 0 AND (Ca_meql + Mg_meql) > 0) THEN \
    SAR = Na_meql/((Ca_meql + Mg_meql)/2)^0.5
160 REM Adjusted SAR is calculated after calcite equilibrium
170 REM ==Residual sodium carbonate index ===========================
180 IF (Alk_meql > 0 AND (Ca_meql + Mg_meql) > 0) THEN \
    RSC = Alk_meql - (Ca_meql + Mg_meql)
190 REM ==Permeability index ========================================
200 IF (Ca_meql + Mg_meql + Na_meql > 0) THEN \
    perm_index = (Na_meql + (Alk_meql)^0.5) / (Ca_meql + Mg_meql + Na_meql) * 100
210 REM ==Kelly ratio ===============================================
220 IF (Ca_meql + Mg_meql > 0) THEN \
    KR = Na_meql / (Ca_meql + Mg_meql)
230 REM ==Percent sodium ============================================
240 IF (Ca_meql + Mg_meql + Na_meql + K_meql > 0) THEN \
    pct_Na = (Na_meql + K_meql) / (Ca_meql + Mg_meql + Na_meql + K_meql)

250 PRINT "Na, meq/L:                            ", STR_F$(Na_meql,12,1)
260 PRINT "K, meq/L:                             ", STR_F$(K_meql,12,1)
270 PRINT "Ca, meq/L:                            ", STR_F$(Ca_meql,12,1)
280 PRINT "Mg, meq/L:                            ", STR_F$(Mg_meql,12,1)
290 PRINT "Alkalinity, meq/L:                    ", STR_F$(Alk_meql,12,1)

300 PRINT "Specific conductance, uS/cm:          ", STR_F$(SC,10,0)
310 PRINT "Total Hardness, mol/L:                ", STR_F$(TH,14,3)
320 PRINT "Total Hardness, mg/L CaCO3:           ", STR_F$(TH_CaCO3,10,0)
330 PRINT "SAR, sodium absorption ratio:         ", STR_F$(SAR,13,2)
340 PRINT "RSC, residual sodium carbonate index: ", STR_F$(RSC,10,0)
350 PRINT "Permiability index:                   ", STR_F$(Perm_index,10,0)
360 PRINT "Kelly ratio:                          ", STR_F$(KR,12,1)
370 PRINT "Percent sodium:                       ", STR_F$(pct_Na,12,1)
-end
SELECTED_OUTPUT 2
-file sar.sel
USER_PUNCH 2
-headings SC TH TH_CaCO3 SAR RSC PI KR Pct_Na Na_meql K_meql Ca_meql Mg_meql Alk_meql
-start
10 Ca_meql = TOT("Ca")*2*1e3*TOT("water")/SOLN_VOL
20 Mg_meql = TOT("Mg")*2*1e3*TOT("water")/SOLN_VOL
30 Na_meql = TOT("Na")*1e3*TOT("water")/SOLN_VOL
40 K_meql =  TOT("K")*1e3*TOT("water")/SOLN_VOL
60 Alk_meql = ALK*1e3*TOT("water")/SOLN_VOL

100 REM ==Total hardness mol/L ======================================
110 TH = (Ca_meql + Mg_meql)/2/1e3
120 REM ==Total hardness mg/L as CaCO3 ==============================
130 TH_CaCO3 = TH * 1000 * GFW("CaCO3")     
140 REM ==Sodium absorption ratio ===================================
150 IF (Alk_meql > 0 AND (Ca_meql + Mg_meql) > 0) THEN \
    SAR = Na_meql*1000/((Ca_meql + Mg_meql)*1000/2)^0.5
160 REM Adjusted SAR is calculated after calcite equilibrium
170 REM ==Residual sodium carbonate index ===========================
180 IF (Alk_meql > 0 AND (Ca_meql + Mg_meql) > 0) THEN \
    RSC = Alk_meql - (Ca_meql + Mg_meql)
190 REM ==Permeability index ========================================
200 IF (Ca_meql + Mg_meql + Na_meql > 0) THEN \
    perm_index = (Na_meql + (Alk_meql)^0.5) / (Ca_meql + Mg_meql + Na_meql) * 100
210 REM ==Kelly ratio ===============================================
220 IF (Ca_meql + Mg_meql > 0) THEN \
    KR = Na_meql / (Ca_meql + Mg_meql)
230 REM ==Percent sodium ============================================
240 IF (Ca_meql + Mg_meql + Na_meql + K_meql > 0) THEN \
    pct_Na = (Na_meql + K_meql) / (Ca_meql + Mg_meql + Na_meql + K_meql)

300 PUNCH SC
310 PUNCH TH
320 PUNCH TH_CaCO3
330 PUNCH SAR
340 PUNCH RSC
350 PUNCH Perm_index
360 PUNCH KR
370 PUNCH pct_Na
380 PUNCH Na_meql, K_meql, Ca_meql, Mg_meql, Alk_meql
-end


[davidjoferreira]

Dear dlparkhurst,
Thank you so much for this piece of script! It worked properly!

I've just made two slight changes:
- In the SAR expression, I've removed the *1000. The original expression doesn't have it and the results appear to be different than mine.

- In the line 330, I've changed for 2 decimal digits.

I've tried with the solution_spread as well and it worked fine!

Thank you so much for this valuable help!

[davidblevy]

Attached are two files I have used in the past...one calculates SAR and adjusted SAR (original method by Ayers and Westcot 1976) and the second file calculates adjusted SAR based on the subsequent 2009 paper by Lesch and Suarez. David