Program Fugco * ---------------------------------------------------------------------| * Calculation of fugacity coefficients for mixtures of H2O, CO2, * and NaCl. Warning—this program should not be used at temperatures * below 350 deg C, pressures below 500 bars, or pseudobinaries of * NaCl above 35 weight percent. The upper temperature limit is * roughly 600 deg C. * Bowers, T.S., Helgeson, H.C., 1985. Fortran programs for generating * fluid inclusion isochores and fugacity coefficients for the system * H2O-CO2-NaCl at high pressures and temperatures. Comput. Geosci. 11, * 203-213. * RELATED: Bowers, T.S., Helgeson, H.C., 1983. Calculation of the * thermodynamic and geochemical consequences of nonideal mixing in the * system H2O-CO2-NaCl on phase relations in geologic systems; equation * of state for H2O-CO2-NaCl fluids at high pressures. Geochim. * Cosmochim. Acta 47, 1247-1275. * Bowers, T.S., Helgeson, H.C., 1983. Calculation of the thermodynamic * and geochemical consequences of nonideal mixing in the system * H2O-CO2-NaCl on phase relations in geologic systems; metamorphic * equilibria at high pressures and temperatures. Am. Min. 68, 1059-1075. * X = Mole fraction, GAM = Fugacity coefficient, * R = Fugacity = X*GAM*PBAR. Species order is 1=CO2, 2=H2O, 3=NaCl * ---------------------------------------------------------------------| real*8 x(5), gam(5), r(3), tk, pat character*12 blank data blank/' '/ in1 = 10 iout1 = 11 open ( in1, file = 'furun.dat', status = 'old') open (iout1, file = 'fuout.dat', status = 'unknown') * Input minimum, maximum, and incremental values of temperature and * pressure. write (*,*) write (*,*) write (*,'(a15)') " Program Fugco" write (*,'(a54)') + " Computes fugacity coefficients at P > 500 bar and" write (*,'(a47)') + " T from 350-600 C in the sytem H2O-CO2-NaCl" write (*,*) write (*,'(a16)') " I/O file names" write (*,'(a31)') " -Input file = ""furun.dat""" write (*,'(a31)') " -Output file = ""fuout.dat""" write (*,*) write (*,'(a19)') " Input file format" write (*,'(a41)') " -Line 1: T(C), P(bars) [free format]" write (*,'(a39,a35)') + " -Additional lines: X(CO2), X(H2O),", + " and X(NaCl) triplets [free format]" write (*,*) write (*,*) read(in1,*) tc, pbar write(iout1,1000) tc, pbar write(*,1000) tc, pbar 1000 format (1X,F5.0,' DegC',2X,F6.0,' Bars',' ',/) write (iout1,1100) write (*,1100) 1100 format (4x, 'XCO2',12x, 'phiCO2',10x, 'fCO2', + 12x, 'XH2O',12x, 'phiH2O',10x, 'fH2O', + 12x,'XNaCl',11x,'phiNaCl', 9x,'fNaCl',/, + 4x,'------',10x,'------',10x,'------', + 10x,'------',10x,'------',10x,'------', + 10x,'------',10x,'------',10x,'------') tk = tc + 273.15d0 pat = pbar/1.01325d0 * X(1) = X(CO2), X(2) = X(H2O), X(3) = X(NaCl) 60 read(in1,*,end=999) x(1), x(2), x(3) call trkmix(tk,pat,x,gam) write( *,1110) (x(m),gam(m),x(m)*gam(m)*pbar, + m = 1,3) write(iout1,1110) (x(m),gam(m),x(m)*gam(m)*pbar, + m = 1,3) 1110 format(9(2x,e14.8)) goto 60 999 stop end SUBROUTINE TRKMIX(T,P,Y,GAM) *======================================================================| * Calculation of fluid mixture properties for the ternary system * using De Santis et al. (1974) * Modifications of the Redlich—Kwong equation together with * A and B parameters for H2O. CO2, and NaCl from Bowers and * Helgeson (1983). * P is in atmospheres, T is in deg kelvin, y(I) is the * mole fraction. Gam(I) is the fucacity coefficient for * the components in the mixture. * Written by T. Bowers (1982) * * To use subroutine trkmix with program denfind (1) change GAM to * VOL in subroutine card and (2) leave out lines 0740 through 1310 * ---------------------------------------------------------------------| real * 8 Y(5), GAM(5) real * 8 AI(3), BI(3) real * 8 AZCO2, BCO2, T, P, R, TCEL, R2T, RT, W, BH, BN, PHI real * 8 AH, AN1, AN2, AN3, BH2O, AZH2O, AH2O, ANACL, AH2OM real * 8 ACO2M, XK, CO2H2O, BSUM, ASUM, A2BMIX, BMIX, Z real * 8 CHI1, CHI2, CHI3, CHI5, CHI6, CHI7CW, DW *======================================================================| * CO2 parameters from De Santis et al. (1974) and Redlich and Kwong * (1949) AZCO2 = 46.d6 BCO2 = 29.7d0 * Define the gas constant in cm3*atm*deg-1*mole-1. R = 82.05d0 TCEL = T - 273.15d0 R2T = R*R*T**2.5d0 RT = R * T**1.5d0 IF(Y(3) .LE. 1.d-8) GO TO 11 * Calculate weight percent NaCl relative to H2O + NaCl W = 100.d0*Y(3)*58.448d0/(Y(3)*58.448d0 + Y(2)*18.016d0) GO TO 12 11 W = 0.d0 * Begin calculating water a and b parameters 12 BH = 14.6d0 BN = - .04420283d0 AH = 4.881243d0 + .1823047d-2*TCEL - .1712269d-4*TCEL*TCEL + + .6479419d-8*TCEL*TCEL*TCEL AN1 = .02636494d0 - .536994d-3*TCEL + .2687074d-5*TCEL*TCEL + - .4321741d-8*TCEL*TCEL*TCEL AN2 = .6802827d-2 - .948023d-4*TCEL + .3770339d-6*TCEL*TCEL + - .5075318d-9*TCEL*TCEL*TCEL AN3 = .5235827d-4 - .3505272d-7*TCEL BH2O = BH + W*BN AZH2O = EXP(AH + W*AN1 + W*W*AN2 + W*W*W*AN3)*10d5 AH2O = 111.3057d0 + 50.70033d0*EXP(-.982646d-2*TCEL) ANACL = -8.05658d0*EXP(-.982646d-2*TCEL) 20 AH2O = AH2O * 10.d5 ANACL = ANACL * 10.d5 AH2OM = AH2O + W*ANACL * CO2 parameter from Holloway (1977) * Interaction coefficient from De Santis et al. (1974) ACO2M = 73.03d0 - .0714*TCEL + 2.157d-5* TCEL*TCEL ACO2M = ACO2M * 10.d5 XK=EXP(-11.071d0+(5953.d0/T)-(2.746d6/(T*T))+(4.646d8/(T*T*T))) CO2H2O = XK * 0.5d0 * R2T CO2H2O = CO2H2O + SQRT(AZCO2 * AZH2O) Y(5) = Y(2) + Y(3) BSUM = BCO2*Y(1) + BH2O*Y(5) ASUM = Y(1)*Y(1)*ACO2M + Y(5)*Y(5)*AH2OM + 2.d0*Y(1)*Y(5)*CO2H2O A2BMIX = ASUM / (BSUM * RT) BMIX = BSUM / (R * T) CALL REDKW(T,P,BMIX,A2BMIX,Z) * Calculate V VOL = R*T*Z/P * For use of this subroutine with program denfind. Lines 0740 * through 1310 need not be included. * Begin calculation of fugacity coefficients of the components CHI1 = DLOG(VOL/(VOL-BSUM)) CHI2 = DLOG((VOL+BSUM)/VOL) CHI3 = BSUM/ (VOL + BSUM) DW = 58.448d0/18.016d0/100.d0 CW = 18.016d0/58.448d0/100.d0 IF(ABS(W-100.d0).LT.1.d-6) GO TO 16 GO TO 17 16 BI(1) = BCO2 BI(2) = BH2O - CW*BN*W*W*Y(5)/Y(3) BI(3) = BH2O AI(1) = 2.d0*Y(1)*ACO2M + 2.d0*Y(5)*CO2H2O AI(2) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O - + CW*W*W*Y(5)*(ANACL) + *y(5)/y(3) - Y(1)*Y(5)*SQRT(AZCO2*AZH2O)*(AN1 + 2.d0*W*AN2 + + 3.d0*W*W*AN3)*CW*W*W/y(3) AI(3) = 2.d0*Y(3)*AH2OM GO TO 86 17 BI(1) = BCO2 IF(W.LE.1.d-8) GO TO 85 BI(2) = BH2O - CW*BN*W*W*Y(5)/Y(3) BI(3) = BH2O + DW*BN*(100.d0 - W)*(100.d0 - W)*Y(5)/Y(2) AI(1) = 2.d0*Y(1)*ACO2M + 2.d0*Y(5)*CO2H2O AI(2) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O - + CW*W*W*Y(5)*(ANACL) + *Y(5)/Y(3) - y(1)*Y(5)*SQRT(AZCO2*AZH2O)*(AN1 + 2.d0*W*AN2 + + 3.d0*W*W*AN3)*CW*W*W/Y(3) AI(3) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O + + DW*(100.d0 - W)*(ANACL) + *(100.d0 - W)*Y(5)*Y(5)/Y(2) + Y(1)*Y(5)*SQRT(AZCO2*AZH2O) + *(AN1 + 2.d0*W*AN2 + 3.d0*W*W*AN3)*DW*(100.d0-W) + *(100.d0-W)/Y(2) GO TO 86 85 AI(2) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O AI(1) = 2.d0*Y(1)*ACO2M + 2.d0*Y(5)*CO2H2O IF(Y(2).LT.1.d-8) GO TO 91 AI(3) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O + + DW*(100.d0 - W)*(ANACL) + *(100.d0 - w)*y(5)*y(5)/y(2) + y(1)*Y(5)*SQRT(AZCO2*AZH2O) + *(AN1 + 2.d0*W*AN2 + 3.d0*W*W*AN3)*DW*(100.d0 - W) + *(100.d0 - W)/Y(2) BI(3) = BH2O + DW*BN*(100.d0 - W)*(100.d0 - W)*Y(5)/Y(2) GO TO 93 91 AI(3) = 2.d0*y(1)*CO2H2O BI(3) = BH2O 93 BI(2) = BH2O 86 DO 202 J = 1, 3 CHI5 = BI(J)/(VOL-BSUM) CHI6 = (ASUM*BI(J))/(RT*BSUM**2) CHI7 = AI(J)/(RT*BSUM) PHI = CHI1+CHI5-CHI7*CHI2+CHI6*(CHI2-CHI3)-DLOG(Z) PHI = EXP(PHI) 201 GAM(J) = PHI 202 CONTINUE 1000 CONTINUE RETURN END SUBROUTINE REDKW(T,P,B,A2B,Z) *======================================================================| * A subroutine to calculate compressibility factor * with the Redlich—Kwong equation following the method * of Edmister(1968). T is in kelvins, P is in atmospheres * written by John R. Holloway(1973) *======================================================================| * Hmmmm. Not a single variable declaration. Let's add those... Real * 8 T, B, P, BP, TH, A2B, RR, QQ, XN, XM, ARG, FP, Z Real * 8 TANPHI, COSPHI, THETA Real * 8 FAC, PHI, RH, R1, R2, R3, X, XN2, XMM, XNN BP = B * P TH = 1.d0 / 3.d0 IF(A2B.LE.0.0d0) A2B = 0.001d0 RR = -A2B * BP * BP QQ = BP * (A2B - BP - 1.0d0) XN = (QQ/3.d0) + RR - (2.d0/27.d0) XM = QQ - 1.d0 /3.d0 ARG = ((XN*XN) / 4.d0) + (XM*XM*XM/27.d0) IF(ARG.GT.0.0d0) GO TO 5 IF(ARG.LT.0.0d0) GO TO 15 FP = 1.0d0 Z = 1.0d0 RETURN 15 COSPHI = SQRT((XN*XN/4.d0) / (-XM*XM*XM/27.d0)) IF(XN.GT.0.0d0) COSPHI = -COSPHI TANPHI = SQRT(1.d0 - COSpHI*COSPHI) / COSPHI PHI =ATAN(TANPHI)/ 3.0d0 FAC = 2.d0 * SQRT(-XM/3.d0) R1 = FAC * COS(PHI) R2 = FAC * COS(PHI + 2.094395103d0) R3 = FAC * COS(PHI + 4.188790206d0) PHI = 3.d0 * PHI * Sort for largest root RH = R2 IF(R1.GT.R2) RH = R1 IF(R3.GT.RH) RH = R3 * Determine the value of theta IF (RH.EQ.R1) THETA = 0.0d0 IF (RH.EQ.R2) THETA = 2.094395103d0 IF (RH.EQ.R3) THETA = 4.188790206d0 Z = RH + TH GO TO 20 5 X = SQRT(ARG) F = 1.0d0 XN2 = -XN / 2.d0 XMM = XN2 + X IF(XMM.LT.0.0d0) F = -1.0 XMM = F * ((F * XMM) ** TH) F = 1.0d0 XNN = XN2 - X IF(XNN.LT.0.0d0) F = -1.0d0 XNN = F * ((F * XNN) ** TH) Z = XMM + XNN + TH 20 RETURN END