HIBAL Hydrologic Isotopic Balance Model for Paleolake Systems ----------------------------------------------------------------------- World Data Center for Paleoclimatology, Boulder and NOAA Paleoclimatology Program ----------------------------------------------------------------------- NOTE: PLEASE CITE CONTRIBUTORS WHEN USING THIS DATA!!!!! NAME OF DATA SET: HIBAL Hydrologic Isotopic Balance Model for Paleolake Systems LAST UPDATE: 8/2002 (Original Receipt by WDC Paleo) CONTRIBUTOR: Larry Benson, United States Geological Survey IGBP PAGES/WDCA CONTRIBUTION SERIES NUMBER: 2002-054 SUGGESTED DATA CITATION: Benson, L. and F. Paillet, 2002, HIBAL Hydrologic Isotopic Balance Model for Paleolake Systems, IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series # 2002-054. NOAA/NGDC Paleoclimatology Program, Boulder CO, USA. ORIGINAL REFERENCE: Benson, L. and F. Paillet, 2002, HIBAL: a hydrologic-isotopic-balance model for application to paleolake systems. Quaternary Science Reviews, Vol. 21, Issue 12-13, July 2002, pp. 1521-1539. GEOGRAPHIC REGION: Pyramid Lake, Nevada with global applications PERIOD OF RECORD: Holocene DESCRIPTION: Benson and Paillet, 2002, HIBAL Hydrologic Isotopic Balance Model for Paleolake Systems ABSTRACT: A simple hydrologic-isotopic-balance (HIBAL) model for application to paleolake d18O records is presented. Inputs to the model include discharge, on-lake precipitation, evaporation, and the d18O values of these fluid fluxes. Monthly values of climatic parameters that govern the fractionation of d18O and d16O during evaporation have been extracted from historical data sets and held constant in the model. The ability of the model to simulate changes in the hydrologic balance and the d18O evolution of the mixed layer has been demonstrated using measured data from Pyramid Lake, Nevada. Simulations of the response in d18O to step- and periodic-function changes in fluid inputs indicate that the hydrologic balance and d18O values lag climate change. Input of reconstructed river discharges and their d18O values to Pyramid and Walker lakes indicates that minima and maxima in simulated d18O records correspond to minima and maxima in the reconstructed volume records and that the overall shape of the volume and d18O records is similar. The model was also used in a simulation of abrupt oscillations in the d18O values of paleo-Owens Lake, California. HIBAL fortran code: PROGRAM BENSON DIMENSION DD(1000),DFF(1000) DIMENSION ZZF(13),ZZE(13),ZZT(13),DUME(1000) DIMENSION DU(1000),DL(1000),ZZB(13),ZZR(13),ZZP(13) DIMENSION DARAG(1000),DOLAKE(1000) 491 FORMAT(I6) 492 FORMAT(F12.5) 301 FORMAT(F12.4,3X,'MODEL RUN ID NUMBER') 302 FORMAT(F12.4,3X,'EVAP IN METERS') 303 FORMAT(F12.4,3X,'INITIAL LAKE DO18') 304 FORMAT(F12.4,3X,'INITIAL LAKE DEPTH IN M') 305 FORMAT(F12.4,3X,'RIVER INFLOW DO18') 306 FORMAT(F12.4,3X,'PRECIP DO18') 307 FORMAT(F12.4,3X,'ADVECTED AIR DO18') 308 FORMAT(F12.4,3X,'FRACTION OF ADVECTED AIR') 309 FORMAT(F12.4,3X,'YEARLY PRECIP IN M') 310 FORMAT(F12.4,3X,'FRACTION OF INFLOW FROM TAHOE') 311 FORMAT(F12.4,3X,'TAHOE OUTFLOW DO18') 312 FORMAT(F12.4,3X,'DISCHARGE CUTOFF FOR TAHOE IN KM3') 313 FORMAT(I5,3X,'NUMBER OF MODEL YEARS') 314 FORMAT(I5,3X,'INITIAL MODEL YEAR') OPEN(UNIT=2, STATUS='OLD',FILE='FLOW.DAT') OPEN(UNIT=4, FILE='OUTPUT.DAT') OPEN(UNIT=3, STATUS='OLD',FILE='PARAM.DAT') ! PRECIP MONTHY FRACTION ZZP(1)=.146 ZZP(2)=.115 ZZP(3)=.089 ZZP(4)=.053 ZZP(5)=.070 ZZP(6)=.084 ZZP(7)=.030 ZZP(8)=.029 ZZP(9)=.030 ZZP(10)=.055 ZZP(11)=.051 ZZP(12)=.109 ZZP(13)=.129 ! MIXING DEPTH BY MONTH ZZB(1)=1000.0 ZZB(2)=1000.0 ZZB(3)=21.0 ZZB(4)=46.0 ZZB(5)=54.0 ZZB(6)=58.0 ZZB(7)=15.0 ZZB(8)=13.0 ZZB(9)=16.0 ZZB(10)=18.0 ZZB(11)=21.0 ZZB(12)=28.0 ZZB(13)=47.0 ! LAKE TEMPERATURE BY MONTH ZZT(1)=6.66 ZZT(2)=7.08 ZZT(3)=7.83 ZZT(4)=10.05 ZZT(5)=13.07 ZZT(6)=16.41 ZZT(7)=20.77 ZZT(8)=23.15 ZZT(9)=22.34 ZZT(10)=20.00 ZZT(11)=16.61 ZZT(12)=12.62 ZZT(13)=9.02 ! RELATIVE HUMIDITY FRACTION BY MONTH ZZR(1)=0.73 ZZR(2)=0.68 ZZR(3)=0.55 ZZR(4)=0.47 ZZR(5)=0.47 ZZR(6)=0.46 ZZR(7)=0.35 ZZR(8)=0.37 ZZR(9)=0.32 ZZR(10)=0.41 ZZR(11)=0.50 ZZR(12)=0.64 ZZR(13)=0.67 ! RIVER INFLOW MONTHY FRACTION ZZF(1)=0.0543 ZZF(2)=0.0625 ZZF(3)=0.0731 ZZF(4)=0.1127 ZZF(5)=0.1696 ZZF(6)=0.1543 ZZF(7)=0.0882 ZZF(8)=0.0559 ZZF(9)=0.0491 ZZF(10)=0.0443 ZZF(11)=0.0381 ZZF(12)=0.0436 ZZF(13)=0.0542 ! EVAPORATION MONTHLY FRACTION ZZE(1)=0.033 ZZE(2)=0.028 ZZE(3)=0.035 ZZE(4)=0.059 ZZE(5)=0.079 ZZE(6)=0.098 ZZE(7)=0.129 ZZE(8)=0.133 ZZE(9)=0.124 ZZE(10)=0.094 ZZE(11)=0.069 ZZE(12)=0.062 ZZE(13)=0.055 NTEST=9999 DFFDUM=0.5 106 FORMAT(' ENTER EVAPORATION IN METERS') ! READ INPUT PARAMETERS FROM FILE READ(2,491) NYR READ(2,491) IYR READ(3,492) CNUM READ(3,492) EVAP READ(3,492) DOINIT READ(3,492) DINIT READ(3,492) DOXIN READ(3,492) DOPRE READ(3,492) DOADV READ(3,492) FFAD READ(3,492) CPRE READ(3,492) TFRAC READ(3,492) TDO READ(3,492) FCUT XDUM1=DOXIN XDUM2=TFRAC*TDO+(1.0-TFRAC)*DOXIN IYR=IYR-1 ZPRE=CPRE ! INSURE AL INPPUT CONVERTED TO KM3 CPRE=CPRE/1000.00 RHUMO=1.0 RRAD=(DOADV/1000.0)+1.0 NDEX=10 105 FORMAT(F10.5) ZVAP=EVAP EVAP=EVAP/1000.00 ! INPUT RUNOFF INFLOW FROM FILE DO 28 I=1,NYR READ(2,492) DFF(I) 28 CONTINUE ! DEFINE AREA AND VOLUME AS FUNCTION OF LAKE DEPTH A0=9.522 A1=6.319 A2=-.01332 A3=-.0003194 A4=.0000025 V0=-.1405 V1=.0523 V2=.001974 DM=123.00 VMAX=V0+V1*DM+V2*DM**2 AMAX=A0+A1*DM+A2*DM**2+A3*DM**3+A4*DM**4 D=DINIT V=V0+V1*D+V2*D**2 A=A0+A1*D+A2*D**2+A3*D**3+A4*D**4 SUM=0.0 TEMPO=273.15 DUMOX=DOINIT DBOT=-10.00 VTOP=V VBOT=0.00 DOXUP=DUMOX DOXDN=DUMOX ! START YEARLY CYCLE DO 800 I=1,NYR IF(I.EQ.2)ALF1=ALFE IF(I.EQ.2)TEMP1=TEMP SUM1=0.0 SUM2=0.0 DOXIN=XDUM2 IF(DFF(I).LT.FCUT) DOXIN=XDUM1 ! START MONTHLY CYCLE DO 600 J=1,13 ! COMPUTE DO FRACTIONATION FACTORS BASED ON RH AND T RHUM=ZZR(J)*RHUMO TEMP=TEMPO+ZZT(J) TEMP2=TEMP*TEMP ALFE=1137.0/TEMP2-0.4156/TEMP-0.0020667 ALFE=EXP(ALFE) ALFK=0.994 ALF=ALFK/ALFE ALFARG=2.78*10.00**6 ALFARG=(ALFARG/TEMP2)-2.29 ! ADD RUNOFF AND PRECIP INPUTS TO LAKE CINA=ZZF(J)*DFF(I) CINB=ZZP(J)*CPRE*A CIN=CINA+CINB ! MIX RUNOFF AND PRECIP INPUTS DOXC=(DOXIN*CINA+DOPRE*CINB)/CIN SUM1=SUM1+CIN VBASE=V ABASE=A VOX=V V=VOX+CIN DDC=1000.0*CIN/ABASE DDC=DDC/10.00 DO 240 L=1,5 DO 200 K=1,50 D=D+DDC V=V0+V1*D+V2*D**2 DDV=V-VBASE IF(DDV.GT.CIN) GO TO 220 200 CONTINUE 220 D=D-DDC DDC=DDC/10.00 240 CONTINUE V=V0+V1*D+V2*D**2 VTEST=V A=A0+A1*D+A2*D**2+A3*D**3+A4*D**4 VBASE=V AINFL=1000.0*CIN/A VOLD=VTOP VTOP=V-VBOT ! MIX INPUTS WITH TOP LAKE LAYER DUMOX=(DOXUP*VOLD+DOXC*CIN)/VTOP DOXUP=DUMOX ! CHECK TO SEE IF LAKE SPILLS IF(D.LE.DM) GO TO 250 D=DM V=VMAX A=AMAX VTOP=V-VBOT 250 VOX=VBOT VDUM=V ! EVAPORATE FROM LAKE CVAP=EVAP*A CVAP=ZZE(J)*CVAP ! COMPUTE FRACTIONATION OF DO THROUGH EVAP RLK=(DUMOX/1000.0)+1.00 REV=ALF*RLK/(1.00-RHUM+0.994*RHUM) REV=(RLK/ALFE)-RHUM*RRAD*FFAD REV=REV/((1.0-RHUM)/ALFK+(1.0-FFAD)*RHUM) DOXE=1000.0*(REV-1.00) ! RECORD DO DIFFERENCE LAKE VS ATMOSPHERE FOR 10TH MONTH IF(J.EQ.10) DOLAKE(I)=DOXUP-DOXE SUM2=SUM2+CVAP DVAP=1000.00*CVAP/A V=VDUM-CVAP DDC=1000.0*CVAP/ABASE DDC=DDC/10.00 DO 340 L=1,6 DO 300 K=1,50 D=D-DDC V=V0+V1*D+V2*D**2 DDV=VBASE-V IF(DDV.GT.CVAP) GO TO 320 300 CONTINUE 320 D=D+DDC DDC=DDC/10.00 340 CONTINUE V=V0+V1*D+V2*D**2 VOLD=VTOP VTOP=V-VBOT ! MIX UPPER LAYER AFTER EVAPORATION ZDOX=(DOXUP*VOLD-CVAP*DOXE)/VTOP DOXUP=ZDOX A=A0+A1*D+A2*D**2+A3*D**3+A4*D**4 ADUM=A DUMD=D-ZZB(J) DBOT=DUMD VTOLD=VTOP VBOLD=VBOT ! MOVE BOTTOM OF TOP LAYER AND MIX ACCORDINGLY IF(DBOT.GT.0.00) GO TO 591 DBOT=0.0 VBOT=0.0 VTOP=V DUMOX=(DOXUP*VTOLD+DOXDN*VBOLD)/V DOXUP=DUMOX DOXDN=DUMOX GO TO 599 591 VBOT=V0+V1*DBOT+V2*DBOT**2 VTOP=V-VBOT VDU=VTOP-VTOLD VDL=VBOT-VBOLD IF(VDU.GT.0.0) DOXUP=(DOXUP*VTOLD+DOXDN*VDU)/VTOP IF(VDL.GT.0.0) DOXDN=(DOXDN*VBOLD+DOXUP*VDL)/VBOT ! RECORD ISOTOPE VALUES FOR 10TH MONTH 599 IF(J.NE.10) GO TO 600 DU(I)=DOXUP DL(I)=DOXDN DARAG(I)=DU(I)+ALFARG 600 CONTINUE ! END OF MONTHLY LOOP DD(I)=D DS=I*1.00 DUME(I)=SUM1-SUM2 800 CONTINUE ! END OF YEARLY LOOP IDEX=0 WRITE(6,107) 107 FORMAT(' YEAR WATER LEVEL DELOX') WRITE(4,301) CNUM WRITE(4,302) ZVAP WRITE(4,303) DOINIT WRITE(4,304) DINIT WRITE(4,305) DOXIN WRITE(4,306) DOPRE WRITE(4,307) DOADV WRITE(4,308) FFAD WRITE(4,309) ZPRE WRITE(4,310) TFRAC WRITE(4,311) TDO WRITE(4,312) FCUT WRITE(4,316) 316 FORMAT(7X,'YR',10X,'FLOW',7X,'DEPTH',8X,'VOL',9X,'DO',8X,'ARAG',8X,'DODIF') ! OUTPUT BY YEAR SENT TO FILE DO 900 I=1,NYR K=I VDUM=V0+V1*DD(I)+V2*DD(I)*DD(I) KYR=IYR+I CKYR=1.0*KYR IDEX=IDEX+1 IF(IDEX.LT.1) GO TO 900 IDEX=0 DS=I*1.00 IF(I.EQ.2) WRITE(6,104)K,DD(I),DU(I),DUME(I) WRITE(4,493) CKYR,DFF(I),DD(I),VDUM,DU(I),DARAG(I),DOLAKE(I) 104 FORMAT(I5,3X,3F10.3) 493 FORMAT(7F12.5) IF(I.EQ.2) READ(5,102) LDEX 900 CONTINUE CLOSE(UNIT=2) CLOSE(UNIT=3) CLOSE(UNIT=4) 101 FORMAT(2E12.4) 102 FORMAT(I3) 109 FORMAT(2I5) 999 STOP END