Skip Links

USGS - science for a changing world

Scientific Investigations Report 2013–5016


Macrophyte and pH Buffering Updates to the Klamath River Water-Quality Model Upstream of Keno Dam, Oregon


Appendix B. CE-QUAL-W2 pH and Alkalinity Code Changes

The 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:
Director, Oregon Water Science Center
U.S. Geological Survey
2130 SW 5th Avenue
Portland, Oregon 97201
http://or.water.usgs.gov

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.

Accessibility FOIA Privacy Policies and Notices

Take Pride in America logo USA.gov logo 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