Scientific Investigations Report 2013–5016
Appendix B. CE-QUAL-W2 pH and Alkalinity Code ChangesThe most important code changes made to the CE-QUAL-W2 pH and alkalinity routines are included in this section. For a definitive list of all changes, those who are interested should compare the new and old source files electronically. The enhanced pH buffering routines and the nonconservative alkalinity algorithms can be turned on or off by the user through two new global input parameters in the control file on the MISCELL card; that modified card is shown below:
MISCELL NDAY PHBUFC NCALKC
100 ON ON
The PHBUFC variable (ON/OFF) controls the use of the enhanced pH buffering routines; the NCALKC variable (ON/OFF) controls the use of the nonconservative alkalinity algorithms. If the enhanced pH buffering routines are turned on, then the model will read a new input file named “ph_buffering.npt” that takes the following form:
Enhanced pH Buffering Input File for CE-QUAL-W2
BUFTYPE NH4BUFC PO4BUFC OMBUFC
ON ON ON
OM TYPE OMTYPE NAG POMBUFC
DIST 2 OFF
DENSITY SDEN SDEN SDEN SDEN SDEN SDEN SDEN SDEN SDEN
0.14 0.10
pK VALS PK PK PK PK PK PK PK PK PK
4.5 9.6
STD DEV PKSD PDSD PDSD PDSD PDSD PDSD PDSD PDSD PDSD
1.2 1.0
--------------------------------------------------------------------------------
Input Variables:
NH4BUFC ON/OFF, specifies whether ammonia/ammonium is included in pH buffering
PO4BUFC ON/OFF, specifies whether phosphoric acid is included in pH buffering
OMBUFC ON/OFF, specifies whether organic matter is included in pH buffering
OMTYPE DIST or MONO
where DIST specifies one or more Gaussian distributions of pKa values,
or MONO specifies a collection of discrete pKa values
NAG the number of acid/base groups to model, either as the means of
Gaussian distributions of pKa values or as discrete monoprotic acids
POMBUFC ON/OFF, specifies whether POM is included in OM buffering
where ON indicates that OM buffering includes both DOM and POM
OFF indicates that OM buffering includes only DOM
SDEN site density, in moles of acid/base sites per mole of carbon in OM
PK the pKa values (negative log10 of the acid dissociation constant),
specified either as the mean of a distribution or a discrete value
PKSD the standard deviation for a Gaussian distribution of pKa values
(ignored when specifying an OMTYPE of MONO)
When the model is run with the enhanced pH buffering capability turned on, the code will generate a new output file named “ph_buffering.opt” based on the new inputs from the ph_buffering.npt file. The object is to echo some of the inputs and provide the user with information on the distribution of site densities and pKa values for the organic matter acids. The following is the output that corresponds to the input file shown above: Enhanced pH buffering output file Ammonia buffering: ON Phosphate buffering: ON OM buffering: ON POM buffering: OFF OM buffer type: DIST Inputs: Group Site density pKa std.dev. 1 0.1400 4.500 1.200 2 0.1000 9.600 1.000 Modeled: Group Site density pKa 1 0.0001 0.500 2 0.0003 1.000 3 0.0010 1.500 4 0.0027 2.000 5 0.0058 2.500 6 0.0107 3.000 7 0.0164 3.500 8 0.0213 4.000 9 0.0233 4.500 10 0.0213 5.000 11 0.0165 5.500 12 0.0107 6.000 13 0.0060 6.500 14 0.0033 7.000 15 0.0032 7.500 16 0.0059 8.000 17 0.0110 8.500 18 0.0167 9.000 19 0.0199 9.500 20 0.0184 10.000 21 0.0133 10.500 22 0.0075 11.000 23 0.0033 11.500 24 0.0011 12.000 25 0.0003 12.500 26 0.0001 13.000 27 0.0000 13.500 The “ph_buffering.npt” input file is read by using new code in the input.f90 source file, including the following:
! Initialize variables for enhanced pH buffering ! entire section ! SR 01/01/12
IF (CONSTITUENTS .AND. PHBUFC == ’ ON’) THEN
OPEN (298,FILE=’ph_buffering.npt’,STATUS=’OLD’)
READ (298,’(///8X,3A8)’) NH4BUFC, PO4BUFC, OMBUFC
READ (298,’(//8X,A8,I8,A8)’) OMTYPE, NAGI, POMBUFC
ALLOCATE (SDENI(NAGI),PKI(NAGI),PKSD(NAGI))
READ (298,’(//(:8X,9F8.0))’) (SDENI(J), J=1,NAGI)
READ (298,’(//(:8X,9F8.0))’) (PKI(J), J=1,NAGI)
READ (298,’(//(:8X,9F8.0))’) (PKSD(J), J=1,NAGI)
CLOSE (298)
AMMONIA_BUFFERING = NH4BUFC == ’ ON’
PHOSPHATE_BUFFERING = PO4BUFC == ’ ON’
OM_BUFFERING = OMBUFC == ’ ON’
POM_BUFFERING = POMBUFC == ’ ON’ .AND. OM_BUFFERING
IF (OM_BUFFERING) THEN
SDENI = ABS(SDENI)
IF (OMTYPE == ’ DIST’) THEN
IF (ANY(PKSD <= 0)) THEN
WARNING_OPEN = .TRUE.
WRITE (WRN,’(A)’) ’WARNING -- PKSD inputs in the ph_buffering.npt file must be greater than zero.’
WRITE (WRN,’(A/)’) ’Please fix your inputs. For now, PKSD values of zero will be set to 1.’
END IF
DO JA=1,NAGI
IF (PKSD(JA) <= 0) PKSD(JA) = 1.0
END DO
NAG = 27
ALLOCATE (SDEN(NAG),PK(NAG),FRACT(NAG))
SDEN = 0.0
DO J=1,NAG
PK(J) = 0.5*J
END DO
DO JA=1,NAGI
SUM = 0.0
DO J=1,NAG
FRACT(J) = EXP(-0.5*(((PK(J)-PKI(JA))/PKSD(JA))**2))
SUM = SUM+FRACT(J)
END DO
DO J=1,NAG
SDEN(J) = SDEN(J)+SDENI(JA)*FRACT(J)/SUM
END DO
END DO
ELSE
ALLOCATE (SDEN(NAGI),PK(NAGI))
NAG = NAGI
SDEN = SDENI
PK = PKI
OMTYPE = ’ MONO’
END IF
END IF
OPEN (299,FILE=’ph_buffering.opt’,STATUS=’UNKNOWN’)
WRITE (299,’(A/)’) ’Enhanced pH buffering output file’
WRITE (299,’(2A)’) ’Ammonia buffering: ’, ADJUSTL(TRIM(NH4BUFC))
WRITE (299,’(2A)’) ’Phosphate buffering: ’, ADJUSTL(TRIM(PO4BUFC))
WRITE (299,’(2A)’) ’OM buffering: ’, ADJUSTL(TRIM(OMBUFC))
IF (OM_BUFFERING) THEN
WRITE (299,’(2A)’) ’POM buffering: ’, ADJUSTL(TRIM(POMBUFC))
WRITE (299,’(2A)’) ’OM buffer type: ’, ADJUSTL(TRIM(OMTYPE))
WRITE (299,’(/A/A)’) ’Inputs:’,’Group Site density pKa std.dev.’
DO JA=1,NAGI
IF (OMTYPE == ’ DIST’) THEN
WRITE (299,’(1X,I3,5X,F8.4,4X,F6.3,3X,F6.3)’) JA, SDENI(JA), PKI(JA), PKSD(JA)
ELSE
WRITE (299,’(1X,I3,5X,F8.4,4X,F6.3,3X,A)’) JA, SDENI(JA), PKI(JA), ’ N/A’
END IF
END DO
WRITE (299,’(/A/A)’) ’Modeled:’,’Group Site density pKa’
DO JA=1,NAG
WRITE (299,’(1X,I3,5X,F8.4,4X,F6.3)’) JA, SDEN(JA), PK(JA)
END DO
END IF
CLOSE (299)
ELSE
AMMONIA_BUFFERING = .FALSE.
PHOSPHATE_BUFFERING = .FALSE.
OM_BUFFERING = .FALSE.
POM_BUFFERING = .FALSE.
END IF
In the water-quality.f90 source file, the PH_CO2 subroutine was modified to include the optional buffering by ammonia, phosphoric acid, and organic matter. The modified code is:
ENTRY PH_CO2 ! Enhancements added for buffering by ammonia, phosphate, and OM ! SR 01/01/12
! pH and carbonate species
DO I=IU,ID
DO K=KT,KB(I)
T1K = T1(K,I)+273.15
CART = TIC(K,I)/12011. ! SR 01/01/12
ALKT = ALK(K,I)/50044. ! SR 01/01/12
AMMT = NH4(K,I)/14006.74 ! SR 01/01/12
PHOST = PO4(K,I)/30973.762 ! SR 01/01/12
OMCT = (LDOM(K,I)+RDOM(K,I))*ORGC(JW)/12011. ! moles carbon per liter from DOM ! SR 01/01/12
IF (POM_BUFFERING) OMCT = OMCT + (LPOM(K,I)+RPOM(K,I))*ORGC(JW)/12011. ! SR 01/01/12
!**** Ionic strength
IF (FRESH_WATER(JW)) S2 = 2.5E-05*TDS(K,I)
IF (SALT_WATER(JW)) S2 = 1.47E-3+1.9885E-2*TDS(K,I)+3.8E-5*TDS(K,I)*TDS(K,I)
!**** Debye-Huckel terms and activity coefficients
SQRS2 = SQRT(S2)
DH1 = -0.5085*SQRS2/(1.0+1.3124*SQRS2)+4.745694E-03+4.160762E-02*S2-9.284843E-03*S2*S2
DH2 = -2.0340*SQRS2/(1.0+1.4765*SQRS2)+1.205665E-02+9.715745E-02*S2-2.067746E-02*S2*S2
DH3 = -4.5765*SQRS2/(1.0+1.3124*SQRS2) ! extended Debye-Huckel for PO4 ! SR 01/01/12
DHH = -0.5085*SQRS2/(1.0+2.9529*SQRS2) ! extended Debye-Huckel for H+ ion ! SR 01/01/12
H2CO3T = 10.0**(0.0755*S2)
HCO3T = 10.0**DH1
CO3T = 10.0**DH2
PO4T = 10.0**DH3 ! SR 01/01/12
HT = 10.0**DHH ! activity coefficient for H+ ! SR 01/01/12
HPO4T = CO3T ! tabled values similar to those for carbonate ! SR 01/01/12
OHT = HCO3T ! tabled values similar to those for bicarbonate ! SR 01/01/12
H2PO4T = HCO3T ! tabled values similar to those for bicarbonate ! SR 01/01/12
NH4T = HCO3T ! tabled values similar to those for bicarbonate ! SR 01/01/12
NH3T = H2CO3T ! neutral species, set coefficient to same as that for carbonic acid ! SR 01/01/12
H3PO4T = H2CO3T ! neutral species, set coefficient to same as that for carbonic acid ! SR 01/01/12
!**** Temperature adjustment
KW = 10.0**(-283.971 -0.05069842*T1K +13323.0/T1K +102.24447*LOG10(T1K) -1119669.0/(T1K*T1K))/OHT
K1 = 10.0**(-356.3094 -0.06091964*T1K +21834.37/T1K +126.8339 *LOG10(T1K) -1684915 /(T1K*T1K))*H2CO3T/HCO3T
K2 = 10.0**(-107.8871 -0.03252849*T1K + 5151.79/T1K + 38.92561*LOG10(T1K) - 563713.9/(T1K*T1K))*HCO3T/CO3T
KAMM = 10.0**(-0.09018 -2729.92/T1K)*NH4T/NH3T ! SR 01/01/12
KP1 = 10.0**(4.5535 -0.013486*T1K -799.31/T1K)*H3PO4T/H2PO4T ! Bates (1951) ! SR 01/21/12
KP2 = 10.0**(5.3541 -0.019840*T1K -1979.5/T1K)*H2PO4T/HPO4T ! Bates and Acree (1943) ! SR 01/21/12
KP3 = 10.0**(-12.38) *HPO4T/PO4T ! Dean (1985) ! SR 01/01/12
!**** pH evaluation
PHT = -PH(K,I)-2.1
IF (PH(K,I) <= 0.0) PHT = -14.0
INCR = 10.0
DO N=1,3
F = 1.0
INCR = INCR/10.0
ITER = 0
DO WHILE (F > 0.0 .AND. ITER < 12)
PHT = PHT+INCR
HION = 10.0**PHT
F = CART*K1*(HION+2.0*K2)/(HION*HION+K1*HION+K1*K2)+KW/HION-ALKT-HION/HT ! SR 01/01/12
IF (AMMONIA_BUFFERING) THEN ! SR 01/01/12
F = F + AMMT*KAMM/(HION+KAMM) ! SR 01/01/12
END IF ! SR 01/01/12
IF (PHOSPHATE_BUFFERING) THEN ! SR 01/01/12
F = F + PHOST*( KP1*KP2*HION + 2*KP1*KP2*KP3 - HION*HION*HION ) &
/( HION*HION*HION + KP1*HION*HION + KP1*KP2*HION + KP1*KP2*KP3) ! SR 01/01/12
END IF ! SR 01/01/12
IF (OM_BUFFERING) THEN ! SR 01/01/12
DO JA=1,NAG ! SR 01/01/12
F = F + OMCT*SDEN(JA)*( 1.0/(1.0+HION*(10.0**PK(JA))) - 1.0/(1.0+(10.0**(PK(JA)-4.5))) ) ! SR 01/01/12
END DO ! SR 01/01/12
END IF ! SR 01/01/12
ITER = ITER+1
END DO
PHT = PHT-INCR
END DO
!**** pH, carbon dioxide, bicarbonate, and carbonate concentrations
HION = 10.0**PHT
PH(K,I) = -PHT
CO2(K,I) = TIC(K,I)/(1.0+K1/HION+K1*K2/(HION*HION))
HCO3(K,I) = TIC(K,I)/(1.0+HION/K1+K2/HION)
CO3(K,I) = TIC(K,I)/((HION*HION)/(K1*K2)+HION/K2+1.0)
END DO
END DO
RETURN
In the water-quality.f90 source file, a new ALKALINITY subroutine was added to realize the nonconservative components of alkalinity. The new code is:
!********************************************************************************************************************
!** A L K A L I N I T Y **
!********************************************************************************************************************
ENTRY ALKALINITY ! entire subroutine added ! SR 01/01/12
! According to Stumm and Morgan (1996), table 4.5 on page 173:
! Utilization of ammonium during photosynthesis results in an alkalinity decrease: 14 eq. alk per 16 moles ammonium
! Utilization of nitrate during photosynthesis results in an alkalinity increase: 18 eq. alk per 16 moles nitrate
! Production of ammonium during respiration results in an alkalinity increase: 14 eq. alk per 16 moles ammonium
! Nitrification of ammonium results in an alkalinity decrease: 2 eq. alk per 1 mole ammonium
! Denitrification of nitrate (to nitrogen gas) results in an alkalinity increase: 1 eq. alk per 1 mole nitrate
! Alkalinity is represented as mg/L CaCO3 (MW=100.088). CaCO3 has 2 equivalents of alk per mole.
! Nitrogen has an atomic mass of 14.00674. These numbers account for the factor of 50.044/14.00674 used below.
DO I=IU,ID
DO K=KT,KB(I)
ALKSS(K,I) = (50.044/14.00674) * ( 14./16.*(NH4AP(K,I)+NH4EP(K,I)+NH4ZR(K,I)+NH4MR(K,I)-NH4MG(K,I)) &
+ 18./16.*(NO3AG(K,I)+NO3EG(K,I)) &
- 2.*NH4D(K,I) + NO3D(K,I) + NO3SED(K,I)*(1-FNO3SED(JW)) )
END DO
END DO
RETURN
|
First posted March 1, 2013 For additional information contact: Part or all of this report is presented in Portable Document Format (PDF); the latest version of Adobe Reader or similar software is required to view it. Download the latest version of Adobe Reader, free of charge. |
U.S. Department of the Interior |
U.S. Geological Survey
URL: http://pubsdata.usgs.gov/pubs/sir/2013/5016/section11.html
Page Contact Information: GS Pubs Web Contact
Page Last Modified: Wednesday, 27-Feb-2013 13:32:17 EST
