Modeling the effects of evaporation on acid rock drainage in mine rock samples

Acid rock drainage (ARD) is an environmental concern of great importance to the mining industry. Humidity cell testing is often used to evaluate the acid producing potential of rock piles. It was desired to develop a comprehensive model of the processes that occur in humidity cell testing. Specifically, it was desired to develop a model that considered both ARD reactions without neglecting temperature effects due to reaction or evaporation. To this end, a model of the humidity cell was developed using Comsol Multiphysics™. ID isothermal and adiabatic models of evaporation in a humidity cell were developed. Analytic solutions were used to validate these models. A ID nonad iabatic evaporative model was designed to simulate the actual conditions experienced during humidity cell testing. This model was verified using experimental data. A homogeneous 2D nonadiabatic model of evaporation in a humidity cell was developed and also compared to experimental data with good results. Heterogeneity was added to the 2D model in the form of an inclined layer with different permeabilities to investigate the effects of an irregular air-flow pattern on evaporation. The results showed some interesting flow fields and temperature distributions and a few numerical anomalies that needed to be addressed. Finally, a ID humidity cell model was developed which considered four reactions in addition to the leaching and evaporative processes involved in humidity cell testing. Comparison of this model with experimental data was also good. These models will be used to understand the process of mine rock piles weathering, particularly for a selected mining site. The results will be used to develop a more comprehensive field model that will help predict the stability and the geochemical changes that will occur in a rock pile into the future. Better environmental management of these rock piles will occur with this increase in understanding of the mechanisms involved.

MODELING THE EFFECTS OF EVAPORATION ON ACID ROCK DRAINAGE IN MINE ROCK SAMPLES by Paul Scott Evans A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science Department of Chemical Engineering The University of Utah December 2007 EV APORA TlON ,. Copyright Paul Scott Evans 2007 All Rights Reserved © ,. T H E U N I V E R S I T Y OF UTAH G R A D U A T E S C H O OL SUPERVISORY COMMITTEE APPROVAL Paul Scott Evans This thesis has been read by each member of the following supervisory committee and by u 1 Chair: Edward M. Trujillo /' 5 J I J r>~l y - - ^ - V 1_~ W , ( / ^ / t J o A n n Lighty Q~ ~(J Milind Deo THE UNIVERSITY GRADUATE SCHOOL of a thesis submitted by i majority vote has been found to be satisfactory. ( I ?Co/ I 5/r'/,0 --+ Dea • T H E U N I V E R S I T Y OF U T A H G R A D U A T E S C H O OL FINAL READING APPROVAL To the Graduate Council of the University of Utah: I have read the dissertation of Paul Scott Evans jn fm a[ form and fjave found that (1) its format, citations, and bibliographic style are consistent and acceptable; (2) its illustrative materials including figures, tables, and charts are in place; and (3) the final manuscript is satisfactory to the supervisory committee and is ready for submission to The Graduate School. Dat<# , <- ~ , ~ i / Edward M. Trujillo *r Chair: Supervisory Committee Approved for the Major Department T7~ Jo Arirfeighty'Y L" f} Chair/Dean ^ U Approved for the Graduate Council David S. Chapman^ Dean of The Graduate School THE UNIVERSITY UTAH GRADUATE SCHOOL APPROVAL Gr~duate in its final fonn ~ave I) spbmission ~<q, /O()7 Oat ! . <- Edw~rd M. Trujill~' ~ ,= --7 -= Chair: Supervisory Committee JoAnrk'igh;x" U ChairlDean U D'avid Cha;man \ (ARD) potential of the processes that occur in humidity cell testing. Specifically, it was desired to develop a model that considered both ARD reactions without neglecting temperature effects due to reaction or evaporation. To this end, a model of the humidity cell was developed using Comsol Multiphysics™. ID isothermal and adiabatic models of evaporation in a humidity cell were developed. Analytic solutions were used to validate these models. A ID non-ad iabatic evaporative model was designed to simulate the actual conditions experienced during humidity cell testing. This model was verified using experimental data. A homogeneous 2D nonadiabatic model of evaporation in a humidity cell was developed and also compared to experimental data with good results. Heterogeneity was added to the 2D model in the form of an inclined layer with different permeabilities to investigate the effects of an irregular air-flow pattern on evaporation. The results showed some interesting flow fields and temperature distributions and a few numerical anomalies that needed to be addressed. Finally, a ID humidity cell model was developed which considered four reactions in addition to ABSTRACT Acid rock drainage CARD) is an environmental concern of great importance to the mining industry. Humidity cell testing is often used to evaluate the acid producing !potential of rock piles. It was desired to develop a comprehensive model both or: 10 so~utions 10 non­adiabatic 20 20 different 10 I· the leaching and evaporative processes involved in humidity cell testing. Comparison of this model with experimental data was also good. These models will be used to understand the process of mine rock piles weathering, particularly for a selected mining site. The results will be used to develop a more comprehensive field model that will help predict the stability and the geochemical changes that will occur in a rock pile into the future. Better environmental management of these rock piles will occur with this increase in understanding of the mechanisms involved. leach ing deve lop fie ld stabi lity roek wi ll mechani sms v TABLE OF CONTENTS IV ix NOMENCLATURE xv ACKNOWLEDGEMENTS xvii Chapter 1. INTRODUCTION 1 1.1 ARD Chemistry 2 1.2 Humidity Cell Testing 2 1.3 Modeling 4 1.3.1 Surface drying 4 1.3.2 Drying from particle surfaces 5 1.3.3 In situ drying in porous media 5 2. THEORETICAL REVIEW 7 2.1 Properties of Porous Media 7 2.2 Flow in Porous Media 13 2.3 Evaporation in Porous Media 15 3. EXPERIMENTAL SETUP AND PROCEDURE 23 3.1 Apparatus 23 3.2 Procedure 28 3.2.1 Calibration of RTD probes 28 3.2.2 Sieve analysis 28 3.2.3 Preparation 29 3.2.4 Leaching phase 29 3.2.5 Dry air phase 31 TABLE ABSTRACT ......................................................................................................... iv LIST OF FIGURES ..................................................... .................................... .... IX LIST OF TABLES ................................................................................................ xiv ............................................................................................. ................................................. ............................... XV II ........................................................................ 1 1. 1 Chclni stry .................................. .............................. 1.2 Humidity Cell Testing .................................. .................... 2 1.3 Modeling ........................................................................... 4 ' 1.3.1 Surface dry ing .................................................... 4 1.3.2 Drying from particle surfaccs ............................ 5 1.3.3 In situ dry ing in porous rnedia .... ...................... .5 REViEW ........................................................... ................... ........................... ...................................................... ........................................... ......................... 3. 1 .... ....................................... ..... .......................... 23 .................... .................... .................. ................ 28 3.2. 1 ................................. 28 ..................................................... 28 ............................... .. ........................ 29 .................................................. .29 ,3.. 2.5 ...................................................... 31 3.2.6 Wet air phase 31 3.2.7 Standing phase 31 3.3 Future Experiments - Heated Cells 32 4. HUMIDITY CELL MODEL 33 4.1 Development of Model and Overall Assumptions 33 4.2 ID Isothermal Evaporative Model 34 4.3 ID Adiabatic Evaporative Model 37 4.4 1D Nonadiabatic Evaporative Model 38 4.5 2D Homogenous Evaporative Model. 40 4.6 2D Heterogeneous Evaporative Model 40 4.7 ID Reactive Model 44 s J 4.7.2 Weekly cycle 48 5. PRESENTATION OF RESULTS AND DISCUSSION 52 5.1 ID Isothermal Evaporative Model 52 5.2 ID Adiabatic Evaporative Model 55 5.3 1D Nonadiabatic Evaporative Model 62 5.4 2D Homogenous Evaporative Model 65 5.5 2D Heterogeneous Evaporative Model 73 5.5.1 High permeability layer 73 5.5.2 Low permeability layer 80 5.6 ID Reactive Model 80 5.7 Sensitivity Analysis 86 A s 5.7.2 Rate coefficient k s 2 87 kAo2 5.7.4 Solid heat capacity, CPs 88 5.7.5 Activation energy, Ecti 89 6. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 90 6.1 Conclusions 90 6.2 Future Work 91 A. DERIVATION OF EXPRESSIONS USED IN THE MODEL....93 B. ARD REACTIONS TO BE CONSIDERED IN FUTURE 97 vii o I J J-t 3.2.6 Wet air phase ..................................................... 3.2.7 Standing phase ................................................... Future Experiments .................................. HUMIDITY MODEL. ........................................................ ............ 1 D Model... ................................. 10 Model... ................................... 1 D Model... ............................. ModeL .............................. .40 Model... ......................... .40 1 D Model... ....................................................... .44 -4.7.1 Reactions ............................................................ 46 ..................................................... .48 .5. ............... 10 Model... ................................. 10 Model... .................................. ID:Model... ............................. 20 Model... ............................. 5.520 Model. ............................ 5.5.l ..................................... ...................................... 10 Model... ........................................................ .......................................................... '5.7.1 Specific surface area, As ..................................... 87 ks2 ............................................. 5.7.3 Oxygen adsorption rate kA02 ............................. 88 Cps ...................................... Ea] ....................................... ........................................................................... ....................................................................... ...................................................................... Appendices DERIV A nON MODEL. ... IN FUTURE WORK .......................................................................................... 97 1- Vll C. SCRIPT FILE USED TO RUN THE ID REACTIVE D. DETAILED RESULTS OF SENSITIVITY ANALYSIS 109 LEACHATE.... F. HEATING CONSIDERATIONS FOR FUTURE viii I D MODEL ....................................................... ... ............... ............... 101 ............ I09 E. PROCEDURE FOR THE ANALYSIS OF THE LEACHATE .... 120 EXPERIMENTS ....................................................... ............. ...... 127 REFERENCES ......... ............................................................................. .... ........... 133 , ,. VIII Figure Page 2.1 Conceptual diagram of the phases within a porous medium 8 'Conceptual diagram of capillary tube 11 3.1 Photograph of control cell 25 3.2 Photograph of cell F11 with RTD probes '. 26 3.3 Schematic diagram of air flow in experimental setup 27 3.4 Particle size distribution resulting from the sieve analysis of the two samples - crushed marble and a rock sample 30 4.1 Conceptual diagram of the ID evaporation model 35 4.2 The finite element mesh used for 2D homogeneous model of the dry-air phase 41 4.3 The finite element mesh used for the 2D heterogeneous model 42 close up of smoothed corner used for the 2D heterogeneous model 5.1 Saturation as a function of distance from inlet in meters after 50,000 seconds as predicted by the ID isothermal model 53 5.2 Isothermal model predicted saturation within the humidity cell as a function of time at 0.02, 0.05 and 0.08 m from the air jnlet for £4=s~l 5.3 Isothermal model predicted saturation within the humidity cell as a function of time at 0.02, 0.05 and 0.08 m from the air inlet for ^=100,000,000 s"1 56 LIST OF FIGURES 2. I : ............. g 2.2 ',Conceptual ................................. " ............ II ; cell... .............................................................. F II RTO .................................... : ..... ....................... the two samples - crushed marble and a rock sample ....................... 30 4.1 Conceptual diagram of the 10 evaporation model... ......................... 35 4.2 The finite element mesh used for 20 homogeneous model of the dry-air phase .. " ......................................................................... 41 model... .... .42 4.4 A 20 heterogeneous ................................................................................................. 45 after model... .............. inlet kA=100,000 s-I ............................. " ...................................... 54 kA=lOO,OOO,OOO s-' ...................................................................... I· 5.4 Isothermal model predicted saturation within the humidity cell as a function of time at 0.02, 0.05 and 0.08 m from the air inlet £4=100 s"1 5.5 Adiabatic model predicted saturation in the humidity cell as a function of time at 0.02, 0.05 and 0.08 m from inlet 58 5.6 Model predicted saturation in the humidity cell vs. distance from inlet at 15,000 seconds as predicted by the ID adiabatic model 59 5.7 Model predicted temperature drop after 50,000 s as a function • of distance from the inlet , 61 5.8 iExperimentally determined and model predicted /temperatures in humidity cell Fl 1, Week One 63 5.9 .Experimentally measured and model predicted temperatures 5.10 Experimentally determined and model predicted temperatures in humidity cell Fl 1, Week One 66 5.11 Experimentally determined and model predicted temperatures in Fl 1, 5.12 Experimentally measured and model predicted temperatures in the control cell, Week One 68 5.13 Experimentally determined and model predicted temperatures in the control cell, Week Two 69 5.14 Temperature and flow field in humidity cell after 200,000 seconds, as predicted by the 2D homogeneous model ; 70 5.15 Experimental and model predicted temperature drop in humidity cell Fl 1 as a function of time ...71 5.16 Saturation and flow field in humidity cell after 200,000 s 72 5.17 Temperature and flow field in the 2D humidity cell model with a layer of high permeability material after 1x10s s 74 5.18 Saturation and flow field in the humidity cell with a high permeability layer after l x l 0 5 s 75 I' I ,! for kA=1 00 s·' ..................................................................................... 57 inlet... ...................... from 10 model... ..... function ,of distance from the inlet... ......................... , ........................................ 61 ; Experimentally [F 11, ................................... Experimentally ~temperatures in the control cell, Week One ....................................... 64 Fll, ............................................................ humidity cell F 11, Week Two ............................................................ 67 ............................................................. ............................................................ 20 model... ........ : ........... F 11 ............................................................. 5.l6 ................. it 1 xl05 ............................ 1 1 05 ......................................................... x t 5.19 Temperature drop in the humidity cell with a high permeability layer 76 5.20 Temperature and flow field in the humidity cell with a layer of low permeability material after 1x10s s 77 5.21 Saturation and flow field in the humidity cell with a low permeability layer after 1x10s s 78 5.22 Temperature in the humidity cell with a low permeability layer as predicted by the 2D heterogeneous model 79 5.23 Temperature as predicted by the ID reactive model along t >with experimental data from Cell Fl 1, Week Two 81 5.24 Experimentally determined and model predicted pH of leachate 82 5.25 Total concentration of iron in leachate as predicted by the .ID '. 5.26 Experimentally determined and model predicted concentration of sulfate ion in leachate 84 5.27 Concentration of calcium ion in leachate as predicted by the ID reactive humidity cell model along with experimentally determined values from Farney (2004) 85 D.l The effect of a 50% increase in specific surface area on sulfate concentration in the leachate 110 D.2 The effect of a 50% increase in specific surface area on calcium 110 i D.3 The effect of a 50% increase in specific surface area on total iron concentration in the leachate 111 D.4 The effect of a 50% increase in specific surface area on the pH D.5 The effect of a 50% increase in the rate coefficient of reaction two on sulfate concentration in the leachate 112 D.6 The effect of a 50% increase in the rate coefficient of reaction two on the calcium concentration in the leachate 112 , xi layer. .................................................................................................. 1 xl 05 ............................................. lx105 ......................................................... modeL .................................. ';with Fit, ............................ ...... 'J D reactive humidity cell model along with experimental data from Farney ............... : ........................................................................ 83 .................................................................... 10 va~ues .............................................. .............................................................. concentration in the leachate .................................................. : ........... ll 0 ...................................................... of the leachate .................................................................................... 111 ...................................... I 12 50% .............................. {. Xl D.7 two on the total iron concentration in the leachate 113 D.8 two on the pH of the leachate 113 D.9 D. 10 . The effect of a 50% increase in the rate of adsorption of oxygen , ! leachate .• 114 D.l 1 The effect of a 50% increase in the rate of adsorption of oxygen into the aqueous phase on the concentration of iron in the leachate < 115 D. 12 The effect of a 50% increase in the rate of adsorption of oxygen into the aqueous phase on the pH of the leachate 115 D.13 The effect of a 50% increase in the heat capacity of the solid phase on the concentration of sulfate in the leachate 116 D.l 4 of a phase on the concentration of calcium in the leachate 116 D. 15 117 D. 16 The effect of a 50% increase in the heat capacity of the solid phase on the pH of the leachate 117 i D.l 7 reaction 1 on the concentration of sulfate in the leachate 118 1 D.l 8 reaction 1 on the concentration of calcium in the leachate 118 D.l 9 The effect of a 50% increase in the activation energy of reaction 1 on the concentration of iron in the leachate 119 D.20 The effect of a 50% increase in the activation energy of reaction 1 on the pH of the leachate 119 , xii 0.7 The effect of a 50% increase in the rate coefficient of reaction ............................ 0.8 The effect of a 50% increase in the rate coefficient of reaction ............................................................ 0.9 The effect of a 50% increase in the rate of adsorption of oxygen into the aqueous phase on the concentration of sulfate in the leachate ................................. ; ............................................................ 114 0.10 ' 50% oxygen ' into the aqueous phase on the concentration of calcium in the ;Ieachate ......... : .................................................................................... 114 0.11 ofa jnto ~leachate ......................................................................................... , ..... 0.12 .............................. 0.13 ......................... I 16 0.14 The effect ofa 50% increase in the heat capacity of the solid c?ncentration ....................... 1 16 0.15 The effect of a 50% increase in the heat capacity of the solid phase on the concentration of iron in the leachate ............................. 1 17 16 .............................................. ~ ........... 0.17 The effect of a 50% increase in the activation energy of .................. 1 18 0.18 The effect of a 50% increase in the activation energy of ................ 1 18 0.19 ...................... 1 19 .................................................. .119 XII E. 1 A schematic diagram of the ICP torch F. I Photograph of the insulated high temperature cells and air lines 129 F.2 Photograph of the filament heater used to pre-heat the air F.3 A schematic diagram of the airflow to the heated cells 131 u xiii I lep .. . ............. ....................... 126 ...... fi lament pre·before entering the humidity cells ........................... ... .. ..................... 131 ..................... ,. XIII LIST OF TABLES Table Page 4.1 .Parameters used in the ID isothermal evaporative model 36 'Parameters used in the ID adiabatic evaporative model 38 4.3 Parameters used in the ID nonadiabatic evaporative model for Cell Fl 1 and the Control Cell 40 4.4 "Darcy's Law parameters used in the 2D heterogeneous models 44 4.5 Parameters used in the ID reactive model 49 E.l Potential readings for pH solutions at given temperatures 121 E.2 Sulfate standards 123 E.3 A summary ofthe groups and standards used in ICP testing 126 i - 4. 1 Yarameters I D ................ 4.2 :,P arameters 10 adiaba tic model.. .................. I D evaporat ive ~Ce l l F II and the Control Cell. .............................................. ...... ....... .40 4.4 'Darcy's Law parameters used in the 2D heterogeneous models ........ 44 4.5 Parameters use,d in the I D reactive model ......................................... 49 E. I Potential readings for pH solut ions at given temperatures ................ 121 E.2 Sulfate standards .............................................................................. .. 12 3 E.3 A summary of.the groups and standards used in rcp testing ............ 126 ,. NOMENCLATURE Symbol Description Units cp • heat capacity J k g 1 K Deff effective diffusivity m2 s"1 h •height, or hydraulic head m hfg Jenthalpy of vaporization J kg'1 lumped mass transfer coefficient s'1 M .molar mass kg mol M 'Van Genuchten parameter N -Van Genuchten parameter m mass flux kg m'2: psat saturation vapor pressure Pa R universal gas constant mol"1 S saturation T temperature t time s u Darcy or superficial velocity of the vapor phase m s'1 V speed of propagation m s'1 X spatial coordinate m Greek Letters Symbol Description Units a Van Genuchten parameter Pa1 A change a property transformation variable m porosity P density kg m"3 mass concentration of water vapor in the pore kg m"3 matric suction or capillary pressure Pa DescriQtion Cp ,heat capacity DejJ : diffusivity : hfg {enthalpy kA coefficient molar Van .pal vap?r x DescriQtion in i· Jkg-I KI m 2 s- I J kg-I s -I morl kgm -2 -I s J morl KI K ms -I ms -1 Subscripts and Superscripts Symbol Description 0 initial condition a property of air d downstream condition e evaporation 1 denotes the il h component (usually water) i)f irreducible water content / liquid phase ref reference state s solid phase sat saturated conditions T ^thermal u ^upstream condition v vapor phase W .water nnd Superscripts o condition {I irr f ref T i1h usua lly irreduci ble 1/ V IV refercnce . sol id 'saturatcd ;thcrmal )upstream _water ,. XV I ACKNOWLEDGEMENTS I wish to thank MolyCorp for generously sponsoring this study. I would also like to take this opportufiity to express my gratitude to the faculty and staff of the Department of Chemical Engineering at the University of Utah, especially to my committee'members, Dr. JoAnn Lighty and Dr. Milind Deo. With special thanks to Dr. Edward Trujillo, my advisor, for all the guidance and support he has so graciously offered. Finally, I would like to thank my family, and in particular my wife, Shari Evans, for all her patience and encouragement. ",:ish tak!e opportuhity Departmel1t committee~ members, Dea. SpeCiffJ ~dvisor, Finally. family. One of the many environmental issues that must be dealt with at mining sites is acid rock drainage (ARD). When the sulfide bearing rock is disturbed, and comes into contact with air, it may be oxidized and form (among other things) sulfuric acid. The formation of sulfuric acid is of environmental concern. Additionally, the low pH levels associated with ARD promote leaching of metals from the rock which can further contaminate ground water. Therefore, the ability to understand and predict ARD is of great interest to the mining industry. The objective of this research is to further develop a model of the processes that occur in the humidity cell which can be used in conjunction with humidity cell testing as a tool to predict ARD. The model should be able to capture heating and cooling effects due to reaction, evaporation and heat loss to the surroundings. Additionally, the model should be able to predict concentrations of species of interest at any point time as well as in the leachate, and be verified with experimental data. I 1. I CHAPTER 1 INTRODUCTION , Onf I I~ contac,t in , . .. . 2 1.1 - ARD Chemistry Many reactions are involved in the ARD process. Appendix A enumerates many of them. The following reaction (oxidation of pyrite by dissolved oxygen) is the typical acid producing ARD reaction. 2FeS2 ( s )+ 7 0 2 ( a q ) + 2H20( 1 ) - 2 F e 2 + ( a q ) + 4S042 - ( a q ) + 4 H + ( a q ) [1.1] It is important to note that pyrite can also be oxidized through other mechanisms, most notably by ferric ion (Equation 1.2). FeS2(s)+UFe^)+m20(t) KSFPY )\5Fe^} + ISO^ + 1 6 / / ^ [1.2] Additionally, there are many possible "neutralizing" reactions which help to raise the pH of the drainage water. The reaction of calcium carbonate with hydrogen ion (reaction B.5 in Appendix B) is an excellent example of this mechanism. It is important, and not a trivial matter, to determine which reactions are dominant and take them into consideration for a particular mining site. 1.2 - Humidity Cell Testing Modeling of ARD requires knowledge of many rock properties including its acid producing potential, which is a function not only of composition, but also the proximity of certain minerals to the rock surface. Often, rock properties of interest are not known a priori. In order to determine the acid producing potential of a waste 1.1 ] notabjy D S 14 D 3+ 8H 0 KSFPY IS D 2+ 2S0-2 16H+ r e 2(s) + r e(aq) +. 2 (I) ) l' e(aq) 4(uq) (aq) [1.2] S rock pile, samples taken from the pile are sometimes put through an accelerated weathering process. This is referred to as humidity cell testing or kinetic testing. Rock samples are usually crushed, prior to testing, to a size (typically lA inch or smaller) that will fit into the humidity cell (an 8 inch tall, 4 inch diameter plastic cylinder). Typically, a one kilogram sample of dry rock is put through a weekly cycle designed to expose it to air and water. The cycle starts with a leaching phase, in which a predetermined quantity of water (usually 500 mL.) is allowed to percolate through the sample. The excess water, or leachate, is collected and concentrations of ions and metals in solution are measured. These data are used to evaluate which reactions are occurring and what the corresponding rates might be. The dry air phase begins on the day following the leaching. Some of the water is retained in the pores of the sample following the leaching. Dry air is blown up through the bottom of the sample for three days. Due to the small pore throat diameters and high surface area of the crushed rock, the pore water quickly vaporizes to its saturation vapor pressure. Thus the air exiting the humidity cell will typically have a relative humidity near 100% until evaporation has stopped. Temperature probes may be put in the cell in order to observe the magnitude of the evaporative cooling. The dry air phase lasts three days. The final three days of the humidity cell cycle may be a standing period, in which the sample is left undisturbed until the leaching phase begins again. Conversely, a wet-air phase may be used instead of standing. During the wet air phase air is first bubbled through water which is kept at a constant temperature ... 3 Yt tht phas~ surfac,e Conversely; 30 °et al. (2002). conditions! effects5 magnitude,of the temperature drop inside the cells is highly dependent on the flow rate of the inlet air. Flow rates of one liter per minute give temperature drops on the order of 5 degrees Celsius. 1.3-Modeling A humidity attempt to model evaporation. His model is isothermal and therefore does not account for changes in reaction rates due to temperature. Drying of surfaces, particles and porous media has been studied and modeled by many and is well documented in the literature. Surface 4 (typically 30°C) prior to being blown through the sample. This is done to keep the sample from completely drying out. Temperature effects can have a significant impact on the rates of reaction. Many of the reactions associated with ARD are highly exothermic. Temperatures in waste rock piles have been reported as high as 70°C by LeFevre Conversely, evaporation can have a significant cooling effect under certain conditions: Experimental data collected by Eldredge (2004) and Klinker (2005) show the effects of evaporation in humidity cell testing. Their results show that the magnitude.of .inlet I iter Celsiu~. 1.3 - Modeling reactive humid~ty cell model was developed by Lin (1996) who treats the humidity cell as 10 continuously stirred non-ideal reactors in series. Lin did not 1.3.1 Surfa·ce drying Jacobs and Verhoef (1997) used Sherwood numbers to estimate evaporation from soil when transpiration from plants plays a large role. This is a useful study, as Sherwood numbers are very important when estimating evaporation. However, this is not directly applicable because vegetation is not often found in large quantities on waste rock piles. Cairncross (1996) investigated the effect of competing rate processes, namely drying and chemical reaction. Their study dealt primarily with solvents and polymers. Analogies may be drawn to acid rock drainage however. Although polymerization does not occur in ARD, crystallization and precipitation do. Perhaps the most applicable study on evaporation from surfaces is the one done by Wilson et al. (1997). They were able to illustrate the effect of capillary pressure on the total drying rate using a soil water characteristic curve. They, however, did not actually model evaporation. The evaporative fluxes were simply compared to the evaporative flux from a sample in which there was no capillary pressure. 1.3.2 Drying from particle surfaces Viollaz and Suarez (1985) give a short review of mass transfer correlations used to predict evaporative fluxes under various conditions. They then use numerical algorithms to find solutions to the problem of drying shrinking particles. 1.3.3 In situ drying in porous media Silva al. (2000) studied porous media with very small pores (<5 nm.) in which Darcy's law was deemed invalid. They reported a method of estimating flow rates through the use of an effective diffusivity which was a function of many 5 mck et al. Petlhaps Wi lson then' Silv'a et af. «5 III 1- parameters such as porosity, tortuosity and pore size. Because waste rock piles have much larger pores ( » 5 nm.) it is not reasonable to use this approach. In 1986 Stanish et al. developed a numerical algorithm to model drying of wood. They assume local phase equilibrium. They also use an adaptive mesh technique to add resolution to the evaporation front. They have no forced convection in their model as is typically present in humidity cell testing. Stacey and Udell (1997) developed analytical solutions to the problem of drying unsaturated porous media. Their solutions were compared to numerical models and experimental data with reasonable agreement. They conclude that an isothermal -treatment will almost always over predict evaporation rates in porous media. The analytical solutions they developed are shown in Section 2.3. An evaporative humidity cell model was developed by Eldredge (2005). No reactions were considered in his model. His model illustrates the importance of properties of porous media such as porosity and initial saturation. This research will build upon his work. 6 .In sol:utions medi,a saturation, CHAPTER 2 'they matrix Perhaps the most basic property of a porous medium would be its porosity ($). Porosity can be defined as the volume fraction of a porous medium not occupied by the solid. Porosity, while very difficult to measure directly, can be easily determined with knowledge of the bulk (/?b) and solid (ps) densities of the medium. This relationship is shown in equation 2.1, Ps [2.1] where Vp is the pore volume and Vj is the total volume. The bulk density is defined as the mass of the solid (ms) divided by the entire volume of the porous medium (Vr). THEORETICAL REVIEW 2.1 - Properties of Porous Media Porous media, by definition, are heterogeneous. Under the simplest conditions they will be composed of two phases, a solid and a fluid phase. In more complex situations there may exist a gas phase as well as immiscible liquid phases in addition to the solid matnx (see Figure 2.1). To further complicate the issue, there exists a possibility of suspended solid particles in any or all of the fluid phases. (¢). C/Jb) (Ps) relationship is shown in equation 2.1, 2.1 ] V T defined ms) (VT). I Gas Phase t V g Liquid Phase 2 Solid Matrix 1 Figure 2.1 Conceptual diagram of the phases within a porous medium. Note that the pore volume (Vp) can be occupied by up to three fluid phases. 8 f .. . V g Vp I . Liquid Phase 1 VT , ! phases . [2.2] The solid density, also known as the grain density, is defined as the mass of the solid divided by its solid volume (Vs). m ."•'t • [2'31 When more than one fluid is present in the pores of the medium it is referred to as unsaturated. The saturation of a fluid phase (S) is defined as the volume of the fluid divided by the total volume of the pores. V S =-±- [2.4] Typically a liquid (Si) and a gas phase saturation (Sg) is considered. When two liquid phases are present, they are commonly referred to as water (Sw) or aqueous and oil (S0) phase saturations. The word oil can refer to any liquid phase that is immiscible in water. When S is written without a subscript, unless otherwise specified, it generally refers to the liquid phase saturation. Water in small unsaturated pores exists at a pressure below that of the surrounding atmosphere. This is due to the capillarity of the porous medium. The difference between the gas (P g) and liquid (Pi) phase pressures is known as the 9 (Vs). , m 'Ps =_s , Vs [2.3] Wh'en it'is (51 [2.4] SI) Sg) Sw) So) . g) PI) ,. 10 P c = - [2-6] Because the radius of curvature is rarely known a priori, it is useful to relate it rt). cosine of the wetting angle 0) which is the angle the interface makes with the solid surface at the gas-liquid interface (see Figure 2.2). cos 6 When the wetting angle has a value less than 90 degrees, the liquid is referred to as a wetting liquid. Conversely, wetting angles greater than 90 degrees correspond to nonwetting liquids. capillary pressure (Pc). Figure 2.2 illustrates this pressure difference. The capillary pressure at a distance h from the bottom of the meniscus is given by equation 2.5. Pe=Pg-Pt=pgh [2.5] Here, p is the density of the liquid, g is the acceleration of gravity and h is the height depicted in Figure 2.2. Capillary pressure is also related to the surface tension (cr) at the gas-liquid interface and the effective radius of curvature of the meniscus (refj). Pc). it]. Gapillary 0-) i (r efJ). [2.6] to the radius of the tube. (rt). Equation 2.7 shows that the two are related by the (B) r, r =-­elf cos () [2.7] liquid Figure 2.2 Conceptual diagram of a capillary 11 '- 1rel f l! l ft ; h • gas 1 ,. I tube. 12 pgh - cosO [2.8] r. Because porous media generally have irregular pore throats, this means that the smallest pores will have the strongest affinity for water. Capillary pressure, therefore, will vary with saturation. Many correlations have been developed to describe the relationship between capillary pressure and saturation. Two of the most basic and widely known relationships are the Brooks and Corey Model (Equation 2.9) and the van Genuchten Model (Equation 2.10). Seff ~ KPcJ [2.9] Seff4 + ( * P j \ M In equations 2.9 and 2.10 the parameters A, N, M, cib and are empirical constants and are dependent on the pore size distribution of the medium. In the Brooks and Corey Model ctb is often referred to as the air entry value. The air entry value can be thought of as the minimum pressure difference required for air to begin displacing the pore water. Both use an effective saturation (Sefj). Equation 2.11 Combining equations 2.5, 2.6, and 2.7 it can be see that water will rise highest in the pores of smallest diameter in a water wetting medium. 2a pgh = -cose r t , . P9res , , , S =(ab J~ eJJ p c [2.10] ab a apd ab wate~: (Self). 13 _ S S in Flo,w in porous media is most commonly described by Darcy's Law (equation 2.12). Developed by Henry Darcy in 1856, Darcy's Law says that the superficial velocity of a fluid in a porous media is directly proportional to its hydraulic head gradient. q=-KhVh [2.12] The proportionality constant, Kh, is known as the hydraulic conductivity and has units (h), therefore, is the sum of the pressure (y/) and elevation (z) heads. •h = i// + z [2.13] I-defines the effective saturation in terms of the upper saturation limit (Su) and the residual saturation (Sr). These two limits occur when phase saturation (Sg or Si) becomes so low, that the phase is no longer continuous. Thus below Sr, water can no longer flow. Su is often taken to be unity. defines Su) Sr). Sg Sl) Sr, Su S = S-Sr eJJ S-S u r [2.11 ] 2.2 - Flow Porous Media Flo}\' superficial of velocity. Because flow in porous media is typically very slow, the velocity component of the hydraulic head is typically assumed to be negligible. The total head If) [2.13] 14 ^ = - [2.14] PS Reynolds numbers. A Reynolds number of 10 is the typical cutoff value. In engineering applications (as opposed to earth sciences) Darcy's Law is often written ,q=~-VP [2.15] P =P + pgz < [2.16] following. k = ^ PS Darcy's Law in order to account for the smaller area available for flow. Relative permeability (kr ) is defined as the ratio of the effective permeability under the given conditions ke/f) to the saturated permeability. It is typically a number ranging from zero to one, although under certain conditions it can be greater than one. P If/ =­pg Because Darcy's Law neglects the velocity head it is only valid for low using permeability (k) instead of hydraulic conductivity (equation 2.15). , i _ k­, q =--VP Jl P+' The advantage of using equation 2.15 is that permeability is independent of the fluid, whereas hydraulic conductivity is not. The two are related by the following. k = KhJl pg [2.17] When the porous medium is unsaturated an additional term is needed in kr) (kejJ) from 15 K h.eff k K = K k [2.18] sat The subscript sat refers to the property under saturated conditions. 2.3 - Evaporation is typically a mass transfer problem. The porous media in humidity cells (crushed rock samples) is composed of small particles. This is, in essence, atpacked bed reactor. The evaporative flux (N) from the surface of the particles can be given by the following equation (Taylor and Krishna, 1993), where km is the mass transfer coefficient and c is the vapor phase concentration of the volatile component. km can be determined using the Sherwood Number (Sh). The Sherwood Number is usually a function of the Reynolds Number (Re) and the Schmidt Number (Sc). The analog for the Sherwood Number in heat transfer is the Nusselt Number (Nu) which is a function of the Reynolds Number and the Prandlt Number (Pr). [2.19] Nu / ( R e , P r) [2.20] Kh,ejf = Kh,sOI = ke/f ksol Evaporation in Porous Media Ev~poration IS c,ells a,km km = f(Re, Pr) 16 i Sh = • [2.22] . D Sc = - [2.23] D Pr = - [2.24] a Re = ^v [2.25] D represents the diffusivity. The product of the Reynolds Number and the Schmidt Number is known as the mass transfer Peclet Number (Pe,„). [2.26] Sh = f(Re,Sc) [2.21] The Sherwood Number can be estimated using empirical heat transfer correlations for the Nusselt Number by substituting the Schmidt Number for the Prandlt Number (Incropera and DeWitt, 2002). Nusse lt v · =- v ­a Re = lid v 2 .23] o diffus ivlty. Peclct (Pe rn ). lid Pem = D ,. 17 The Peclet Number is generally considered to relate the rate of advection to that of diffusion. Local phase equilibrium in a porous medium can be assumed when the mass transfer Peclet Number is much less than unity. Similarly, a heat transfer Peclet Number (Pe/,) can be defined as the product of the Reynolds Number and the Prandlt Number. Analogous to the mass transfer Peclet Number, low values of the heat transfer Peclet Number indicate local thermal equilibrium. When the porous' medium is heterogeneous it can be difficult to choose a characteristic dimension (d) because the medium is typically composed of particles of many different diameters. There are various approaches to the problem. Either an average or a median particle diameter can be used. Or one may choose to use a particle diameter taken from the sieve data such as d 10 (the sieve size through which only 10 wt. % of the particles can pass). Also, a porous media Reynolds Number (RepM) is often used in place of Equation 2.25 as a better representation of flow in the pores. ud [2.27] a Re ud [2.28] PM i - Peelet Peelet Peelet Peh) Pe h =- a Analogou~ Peelet Peelet partieles partiele partiele 10 partieles Re~nolds RepM = v(1 -¢]UIII 18 Nu = 0.00618 Re [2.32] Equation 2.28 includes porosity (</>) and a shape factor (0) in the determination of the Reynolds Number in porous media. Bird Stewart and Lightfoot (2002) report the following correlation for the Nusselt Number in a packed bed (for small Re). /VM = 2.19(RePr)^ [2.29] Thus, the Sherwood Number can be estimated using equation 2.30. •Sh = 2A9(KeSc)^ [2.30] Other empirical relationships for Nusselt and Sherwood Numbers have more recently been given in the literature. Equation 2.31, given by van der Ham and Browers (1998), gives the Sherwood Number as a function of the Peclet Number (Pe) and the characteristic particle diameter (d), and is said to be valid for a range of Peclet numbers between 5 and 60. \0'3mPe^MdLi2 2004, Ahn et al. reported a Nusselt number for hot nitrogen in a packed bed as the following. . I I' I I. I I I· (¢) B) Nu = Re Pr)13' ~herwood '= 2.19(Re Sc)13' relationsnips Sh = 10-303 PeQ88 dL82 [2.31] In , = 0.00618Re0 8 Pr 19 Sh' = ^ - [2.33] The most recent Sherwood number correlation that Chomsurin and Werth list is from Kennedy and Lennox (1997). This correlation is presented as valid for Reynolds numbers between 0.003 and 8.33, and is given by the following: 5A = 2.13Sc^Re" where n is a parameter given by Equation 2.34. 2 Z 3 n = 4.6 5 R e - o . 2 8 [2 3 5 ] 3 n - l Chomsurin and Werth illustrate the shortcomings of existing Sherwood Number correlations. They give experimentally determined Sherwood Numbers for NAPL but do not recommend a new correlation. Stacey and Udell (1997) studied evaporation in porous media. They developed an analytical solution for the speed of an evaporation wave (Ve) traveling through a partially saturated porous medium by transforming the one-dimensional Finally, Chomsurin and Werth have compiled numerous correlations for Sherwood numbers in columns packed with sand and/or glass beads in the context of non-aqueous phase liquid (NAPL) dissolution. Some of these correlations are for a modified Sherwood Number (Sh'). (2003) (Sh '). Sh'= kmd 2 :D , i (2003) K~nnedy Sh = 22.9 + 2. 13Scli Ren n 2 - 3n = 4.65 Re -0.28 3n-1 [2.34] [2.35] (2003) In Ve) I· V=-, dJ^ 1 [2.36] \SgdP,d+SldPld) In Equation 2.36 the subscript d represents the downstream condition, u is the superficial velocity of the gas phase, pi is the mass concentration of component i (the component being evaporated) and pi is the liquid density. Equation 2.36 can be used to give a .good estimate of the time it will take to completely dry an isothermal unsaturated sample. However, drying is rarely an isothermal process. Stacey and Udell (1997) also give analytical solutions for the temperature drop due to evaporation (AT) and the propagation speed (Vp) of the thermal wave. <t>V p,S,hf Ar = r e H l ' f g [2.37] VT - ua [2.38] pcP In the above equations the subscript a refers to a property of air and the subscript s indicates a property of the solid phase. h/g is the latent heat of vaporization, cP is the specific heat and m is the mass flux. Equations 2.36 through 2.38 can be solved simultaneously to give the solution for an adiabatic system mass balance into wave coordinates. Their solution (given below) assumes phase equilibrium and no migration of water due to capillary pressure gradients. 20 Equati0n condition. , superficia4velocity phase. Pi component being evaporated) and PI is the liquid density. Equation 2.36 can be used to give a .good estimate of the time it will take to completely dry an isothermal Ll1) VT) [2.37] [2.38] t~e hfg Cp m simyltaneously 21 assuming that conductive heat transfer is negligible. Equation 2.37 assumes a constant temperature boundary condition at the inlet. It also assumes that the heat capacity of the vapor phase is approximately the same as air's. This is a good approximation when the partial pressure of the volatile component is much less than the total pressure. Eldredge (2004) used three mass balances and an energy balance to model evaporation in a humidity cell numerically. Two of the mass balances used were for water vapdr and air and are given in equations 2.39 and 2.40. C w and Ca represent the concentrations of water vapor and air respectively. Da,w is the diffusivity of water vapor in air and s g is the volume fraction occupied by the gas phase. The third mass balance, equation 2.41, was for the volume fraction of liquid water ( s L ) . [2.39] dt [2.40] s g Be N-MW [2.41] dt PlVT In equation 2.41, represents the molar mass of water. Change in temperature (T) is calculated from the energy balance given in equation 2.39. assummg evaporatio:n vapor Cw CA DA,&g 2.41 ] MW ener~y 22 k ^ T - { P l C P J L + p g C F J ^ T - h - ^ - - = = V-J- [2.42] dt pCP Note that the volume averaged thermal conductivity (k) and heat capacity (pCp) are used. hfg'refers to the enthalpy of vaporization of water. In addition to the four preceding-mass and energy balances, Eldredge (2004) also used Darcy's law to predict pressure gradients across the humidity cells. In this chapter, the governing equations describing porous media, including flow and evaporation, have been presented. This theoretical background will serve as the basis of the model described in this work. conducti vity hJg .refers of, preceding 'mass Darcy 's pred ict pressure gradients across lhe humidity cell s. thi s ~chapter, now , .. CHAPTER 3 PROCEDURE Modeling of ARD requires knowledge of many rock properties including its acid producing potential, which is a function not only of composition, but the accessibility of air and water to the minerals in the rock. Often, rock properties of interest are not known and must be determined in the laboratory by conducting humidity cell tests. An ASTM (2000) standard procedure has been developed for humidity cell testing. We used a modified version of that procedure for our experiments that enabled us to measure temperatures around and within the rock sample continuously during the weekly cycle, particularly during the dry air phase. The humidity cell is an acrylic cylinder which has dimensions of roughly 4 inches inner diameter and a height of about 8 inches. A perforated plastic disc supports a polypropylene felt pad upon which the rock sample is loaded. Four Omega brand PR-17 RTD probes were fitted to each cell to in order to measure temperature at different locations in the humidity cell during the testing. According to the literature the probes give readings with an accuracy of ±0.3° C for a temperature range of -25° C to 85° C (Pace Scientific, 2004). The four RTD probes EXPERIMENTAL SETUP AND PROCEDURE AID 3.1 Apparatus bra':ld O.3° /. 24 were inserted into the humidity cells, at heights of 0.8, 3.3, 5.7 and 8.0 cm from the top of the polypropylene felt. It should be noted that the rock sample generally did not cover the top probe (see Figures 3.1 and 3.2). Cell F l l (Figure 3.2) contained waste rock from a mining site in the western United States, while crushed marble was used for the control cell (Figure 3.1). An Opto22 control system was used to constantly record the temperatures written toJa comma separated variable (csv) text file by Opto. This file is easily opened in most spreadsheet programs. Each cell has a hole in the bottom which has been tapped to accept a lA NPT fitting (depicted in Figure 3.3). During the leaching phase, an Erlenmeyer flask is connected to the bottom of the humidity cell to collect the leachate. This hole is plugged during all other phases of the humidity cell testing. During part of the humidity cell testing compressed air is fed through a flow meter and then through a Drierite column (Figure 3.3) before entering the humidity cell through a hole underneath the perforated disk (visible in Figures 3.1 and 3.2). If the wet air phase is to be conducted air is diverted from the Drierite column and into the humidifier. The water in the humidifier is maintained at 30°C by an electrical resistance heater. After passing through the humidifier, the air then enters the humidity cell as in the dry air phase. Fll Opt022 measured by the RTD probes at predetermined intervals (5 min.). Historical data are to)a t,4 Figu~e conwressed 25 I Figure 3. 1 Photograph of control ce ll. Note the location of the RTD probes, felt pad, perforated disk and the air in let port. ,. 26 Figure 3.2 Photograph of cell F II wi th RTD probes. 27 Figure 3.3 Schematic diagram of air flow experimental setup (Adapted from Eldredge, 2005). UU. _. Opto-22 System Tempe",,",, "obe Humidity Cell Regulator Bubbler .... ! Collec on E sk Onenla Column Ro,w Meier r---<~ -Flow Junction Air Source Spar €Irs Humidifier fl ow in ,. 28 3.2 Procedure RTD probes used (four per cell). The thermal bath was heated to a temperature of 35° C, which was measured by a mercury thermometer. The ice bath was prepared in a 1000-ml glass beaker. The beaker contained tap water about 3/4-full with several cubes of ice added to fill the beaker. Ten minutes of stirring brought the bath to an equilibrium temperature of 1°C according to the mercury thermometer. The probes were immersed in the hot water bath and allowed to come to equilibrium. Opto 22 was used to record the temperatures measured by the probes. The probes were then removed from the bath and allowed to come to equilibrium with the ambient air. Finally the probes were immersed in the ice bath. Again, sufficient time was given to allow the probes to come to equilibrium with the ice bath. The temperatures measured by the probes were all within the ±0.3°C uncertainty limits of the RTDs when compared to the mercury thermometer. 3.2.2 Sieve analysis Particle size distributions of each sample were measured using a sieve analysis. This was conducted with 11 sieves. The sieve openings ranged in size from 0.075 to 6.700 mm. The mass of each sieve is first determined and recorded. The sieves were loaded into the shaker and arranged in descending order with the most coarse on top down to the finest. A pan is placed underneath the finest sieve to catch any particles smaller than 75 microns. Next a one quarter kilogram portion of the sample is split off (using a splitter) and its mass is recorded. The sample is loaded 3.2.1 Calibration of RTD probes An ice bath and constant temperature water bath were used to calibrate the eight lOOO-4-full of the hot water bath and allowed to come to equilibrium. Opto 22 was used to record the temperatures measured by the probes. The probes were then removed from the bath and allowed to come to equil~brium with the ambient air. Finally the probes were immersed lin:its Sieve from ,- of particle sizes. The results are displayed in Figure 3.4. 3.2.3 Preparation Fl 1, ha.s F l l One kilogram is separated from the rock sample and loaded into the humidity cell whose weight has been recorded previously. The cell is then weighed again. 3.2.4 Leaching phase A separatory funnel is filled with 500 mL of water. The lid of the cell is removed and the valve on the funnel is opened slowly, allowing water to drip out at a rate of approximately one drop per second. Water is allowed to drain out the bottom of the cell into an Erlenmeyer flask. The leachate from the bottoms of the cell are collected and split into three samples. One is used for the measurement of the pH and the conductivity of the sample. Another is used to determine the sulfate concentration. The final sample is used for the inductively coupled plasma (ICP) analysis which was used to measure the concentrations of metals ions (iron, calcium, magnesium; aluminum, selenium, arsenic, and zinc). A more detailed description of the leachate analysis is available in Appendix E. The humidity cell is then allowed to 29 into the top sieve which is then covered. The shaker is allowed to run for 20 minutes. Then, each sieve is weighed again to determine the mass of the sample in each range Two rock samples were used. They were the control, which consisted of crushed marble chips, and F 11, an acid producing sample with high clay content. The control ha$ a much narrower particle size distribution than does FII (Figure 3.4). - lCP) magnesium~ leachat~ 30 0.010 0.100 1.000 10.000 sieve size, mm Particle distribution analysis of the control) l l ) . narrower of the control cell compared to cell F l l . 100.00 80.00 'c" '''Q""" . 6pOO - , C -u<II 4900 <II Q. 20.00 -- Control - marble ( , -- F11 - mste rock / . /) J , ~~ 0.010 0.100 1.000 10.000 Figure 3.4 Part icle size distribut ion (PSD) resulting from the sieve analys is two samples - crushed marble (contro l) and a rock sample ( F II ). Note the much PSD cel l cel l 11 , 31 drain for 24 hours. After the 24 hours the cells are weighed to determine the amount 3.2.5 Dry air phase 1pm 3.2.6 Wet air phase phase the air is no longer passed through the drying column. Instead, it is sparged through 30C water prior to entering the humidity cells. This was done in order to avoid complete dry out of the cells. The wet air phase continued for the final three days of the week, following which, the cell is weighed once again prior to the leaching phase. In the modified procedure this phase has been replaced by the standing phase. 3.2.7 Standing phase following which, the cell is weighed once again prior to repeating the leaching phase. of water retained by the sample. The air stream was passed through the drying column and into the humidity cells. The flow rate was held constant (usually at a value between 1-4 lpm per cell) for a period of three days. OPT022 is used to record the temperature measured by the RTD probes at fixed intervals. At the end of the dry air phase the cells are weighed again. In the past, humidity cell experiments have used a wet air phase. During this 30e ,Following the dry air phase, the humidity cells are left to sit for three days. The lids are left on and no air is permitted to flow in or out of the cells, following 32 Future experiments will involve higher temperature humidity cells, to evaluate the effect of higher temperatures which may be present in rock piles. In preparation for future humidity cell experiments nine heaters have been installed. One is a filament heater through which the air will flow prior to entering the humidity cell during the dry air phase. The other eight are band heaters which have been wrapped around eight humidity cells. Please refer to Appendix F for details concerning future experiments involving heated cells. '. I. 3.3 - Future Experiments - Heated Cells wi ll cel l refer, future experiment"s involv ing healed .. CHAPTER 4 HUMIDITY CELL MODEL Comsol Multiphysics™ (Comsol, 2007), a finite element based computer program, was used to create the humidity cell model. It was chosen because of its easy-to-use graphical user interface (GUI). As multiphysics software, Comsol has many modules which can be used to assist in modeling phenomena as varied as stress analysis, heat and mass transfer, and acoustics. It allows the user to easily include additional reactions and to change boundary and initial conditions. Additionally, Comsol is not a black box; it is a simple matter to modify the governing equations when necessary. Finally, the energy and mass balances used to model the humidity cell can easily be extended to modeling more complex geometries (such as a rock pile) in Comsol. 4.1 - Ambient conditions were assumed to be standard temperature and pressure. The thermal conductivities of air and water as functions of temperature were given by the expressions that Eldredge (2004) used. Other thermal conductivities were assumed to change insignificantly over the temperature range of interest. The enthalpy of vaporization of water was also assumed to change very little over the Corusol Gomsol Finally; 4.1- Development of Model and Overall Assumptions ... 34 temperature interval and was therefore treated as a constant. The Clausius Clapyeron equation was used to predict the saturation vapor pressure of water as a function of temperature. Water's normal boiling point was used for the reference condition. Additionally, the solid phase was assumed to be immobile. ID cell isothermal?case phase water. Concentration (c) of each species is determined by solving the continuity equation at each element. - + V-(uc-DVc)=-R 4.1] dt Evaporation was accounted for using the reaction term (R). The evaporation rate was determined from the flux of water vapor from the liquid surface (given in equation 2.19) which was then multiplied by the specific surface area (As) of the medium. Although the liquid surface area will change with saturation, it was assumed have a constant value equal to that of uniform cubic face center packed spheres. The sensitivity of the model was evaluated with respect to specific surface area and mass transfer coefficient. ,. . 4.2 - ID Isothermal Evaporative Model A '1 D model of the humidity c.ell was developed (see Figure 4.1). The isothennalfcase was run using two mass balances, one for liquid and one for vapor ac -+ V ·(uc- DVc)=R at [4.1 ] detennined coefficient. 4.1 of ID evaporation elements have been shown. 35 Air 6 water vapor out :1 Oiscretized Elements 'Scm 1 Air in Figure Conceptual diagram or the 10 model. Note, not all mesh clements have been shown. 36 because the equilibrium vapor pressure of water at the temperatures of interest is much less than one atmosphere. The liquid phase diffusivity was set to a very small boundaries. The concentration of the water vapor was set to zero at the inlet, and a convective flux boundary condition (no diffusive flux) was applied at the outlet. Table 4.1 provides a summary of the parameters used for the ID isothermal evaporative model, which are applicable to Cell Fl 1. F l l) Parameter Symbol Value Units Reference Porosity 0.36 - Eldredge, 2005 Initial saturation So 0.44 - Klinker, 2006 Initial & ambient temperature To 298 K Klinker, 2006 Air flow rate Qair 1 Ipm Klinker, 2006 Inner diameter of humidity cell D 10 cm Eldredge, 2005 of sample in cell h 2006 Mass transfer Peclet Number Pem 0.06 - calculated transfer kA s"1 2002 The liquid phase velocity was assumed to be zero, while the vapor phase velocity was set to the superficial velocity of the entering air. This neglects a small increase in vapor phase velocity due to evaporation. This assumption was made number, but not zero. This gave stability to the numerical solution, while simulating the wicking that arises due to capillary action. A zero flux boundary condition was given to the liquid mass balance at both 10 FII. Table 4.1 Parameters used in the ID isothermal evaporative model (Cell Fll) Parameter Symbol Value Units Reference Porosity ¢ 0.36 Eldredge, 2005 So & ambient temperature To rate qair of humidity cell 0 Height of sample in cell 8 cm Klinker, 2006 transfer Pedet Number Pem Lumped mass transfer coefficient 100,000 S-1 BSL, 2002 37 For the adiabatic case an energy balance was included. The convection and conduction in porous media module of Comsol was used to find temperature 7) at every point. CP>eq~-V-XeqVT = Q-CPfZ-VT [4.2] porous media, Cp,eq is the volume averaged heat capacity while Cpj is the heat capacity of the mobile fluid. Xeq is an equivalent (often volume averaged) thermal conductivity for the system and Q is a volumetric heat generation term. In this case Q is a heat sink resulting from the evaporative cooling (the product of the enthalpy of vaporization and the evaporation rate). For comparison with the analytical solution, the thermal conductivity was lowered four orders of magnitude in order to mimic the assumption made in equations 2.33 and 2.34 that there is no conductive heat transfer. The same boundary conditions were used as in the isothermal case. Additionally, the inlet was set to have a constant temperature, and a convective flux boundary condition (no conductive flux) was used at the outlet. The adiabatic model 4.3 - ID Adiabatic Evaporative Model (1) c, aT - V' . A V'T Q - Cu· V'T l,eq at eq PJ [4.2] In CP,eq CPJ Aeq conductivity for the system and Q is a volumetric heat generation term. In this case Q is a heat sink resulting from the evaporative cooling (the product of the enthalpy of vaporization and the evaporation rate). For comparison with the analytical solution, the thermal conductivity was lowered four orders of magnitude in order to mimic the flux uses the same parameters as the isothermal model (Table 4.1) with some additional parameters needed for the energy balance. These are given in Table 4.2. (. 38 Table 4.2 - Parameters used in the ID adiabatic evaporative model (Cell F l l ) hfg mot1 Ps kg nrf3 Cps mot1 K 1 Cpair J kg- 1 K_ 1 2007 water Cpwater J kg- 1 K"1 2007 2.5x10'4 W m 1 K 1 4.4 - ID Nonadiabatic Evaporative Model In order to have a numerical model to compare to experimental results, a non­adiabatic model was developed to account for heat transfer with the surroundings. Due to the one-dimensional nature of this model some simplifications had to be made with regard to heat transfer in the radial direction. Incropera and DeWitt (2002) give the steady state solution to the radial heat transfer rate (qr) as the following. T-T <lr = 77~r\ T^7\ [4-3] ln| r 2nLX, 2nLX... InLrh In Equation 4.3, the subscripts r and w refer to properties of the rock and the humidity cell wall, respectively, r, and r0 are the inner and outer radii of the cylinder and h is a convective heat transfer coefficient. The temperature of the surface of the cell was assumed to have a constant value equal to the ambient temperature (7^). This is the same as saying that the last resistance term in the denominator is much less than either of the other two terms. Additionally, the first resistance term vanishes iD Fll) Parameter Symbol Value Units Reference Enthalpy of vaporization h,g 44 kJ mor1 NIST, 2007 Density of solids 2700 kg m-3 Eldredge, 2005 Specific heat of solids Cps 858 J mor1 K1 Eldredge, 2005 Specific heat of air Cpa;r 1010 J kg-1 K-1 NIST, 2007 Specific heat of water CPwater 4180 J kg-1 K-1 NIST, 2007 Thermal conductivity of solids As 2.5x10-4 Wm-1 K 1 Eldredge, 2005 ~ iD non-adiabatic De Witt qr) [4.3] rand respectively. rj ro Ta,). 39 under the assumption that the rock has no radial temperature gradient (it is a ID problem). In order to be included in the volumetric heat generation term, qr must be divided by the total volume of the porous medium. Thus, as might be expected, the volumetric heat generation due to radial conduction (Qr) is independent of the length of the cell (L). The boundary conditions used were all the same as in the previous model, but with one exception. The temperature was not held constant at the inlet. Instead, a heat flux boundary condition was applied. The heat flux was determined from the enthalpy and flow rate of the incoming air. To model Cell F l l , the nonadiabatic model uses the same parameters as the adiabatic model, with the addition of the thermal conductivities of the humidity cell wall and the insulation used by Klinker (2005). These are given in Table 4.3. The control cell was also simulated using the ID non-adiabatic model. The differences in the thermal properties of the Control Cell as compared with Cell Fl 1 were insignificant. Therefore, when adjusting the model to represent the control cell, the only parameters that needed to be changed were the porosity and the initial saturation, which were 45% and 16%, respectively (Eldredge, 2005). [4.4] 10 qr Qr) (L). Qr = ~A)T-TJ k.2 In(;{. ) ,I r I Th~ FII, 10 F II 40 ID F l l W m_ 1 K"1 K W m"1 K"1 Ap W nrf1 K_1 A,- ' 4.5 - Model A two-dimensional axis-symmetric model of the evaporation phase was also developed.1 The mass balances used were the same as those already described for the ID problem. The energy balance is essentially the same as well, although Equation 4.4 was omitted. The temperature of the wall of the humidity cell was held constant. A zero flux boundary condition was applied to both mass balance equations along the cell wall. All the parameters used in this model were the same as those of the ID nonadiabatic model (See Table 4.3). Figure 4.2 shows the finite element mesh that was used in the 2D homogeneous evaporative model. 4.6 - 2D Heterogeneous Evaporative Model In the heterogeneous model an inclined layer was introduced to the domain in order to examine the effect of an irregular flow field (Figure 4.3). should be noted that this layer, in cylindrical coordinates, is actually shaped like a hollow cone. Two heterogeneous scenarios were run. In the first one, the inclined layer was given a permeability four orders of magnitude higher than the surrounding rock. In the Table 4.3 - Parameters used in the 1D nonadiabatic evaporative model for Cell Fll and the Control Cell Parameter Symbol Value Units Reference Thermal conductivity of solids As 2.5 Wm-1 K-1 Eldredge, 2005 Thermal conductivity of Wm-1 K-1 Incropera & DeWitt, plexiglass Ap 0.2 2002 Thermal conductivity of Incropera & DeWitt, Insulation Ai 0.08 Wm-1 K-1 2002 4.5"':' 2D Homogenous Evaporative Model developed:' I D energ):, It /. 41 0.08 0.07 0.06 0 05 3 £ 0.04 0.03 0.02 0.01 • • Rock Sample i Sim i V ? - . ' . i . v *• V ' Insulation Plastic Container : j-'-V"' '• 0.01 0 02 0.03 0.04 0.05 0.06 0.07 0.08 meters Figure 2D dry-air higher density elements scale in meters. ~ 2< ~ 0 .06 o.oy E 0,0,"f 0.0 1 o . , o . . . I , , ~. co. ' . I';. . " . ' .. . ' " 0 . . I , .. I I . , . , I' . .I· ' • 0 . I, Figu re 4.2 The finite element mesh used for the 20 homogeneous model of the dry­air phase. Note the densi ty of mesh clements along boundaries. The sca le is , 42 Figure 4.3 The finite element mesh used for the 2D heterogeneous model. Note the high density of mesh elements at the transition between zones of different permeability. The scale is in meters. 0.08 0.07 0.06 0.05 v~" ~ v 0.04 E 0.03 0.02 om 0 o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 meters lIsed 20 densi ty different penncability. ,. 43 Table 4.4 - Darcy's Law parameters used in the 2D heterogeneous models 1x1 cr5 2 1x1 tr9 2 Incropera DeWitt, Ma 1.8x10s Pas second, the permeabilities were switched such that the inclined layer had a lower permeability. This made it necessary to use Darcy's Law to find the flow field. For transient problems the following equation is solved. \p{\ - <f)+a^- + V--V(p + pgz) = Qs [4.5] dt jU P and represent the compressibility of the solid and fluid phases. The solid was assumed to be incompressible and air was given the compressibility of an ideal gas. The evaporation rate was used for the source term (Qs)- The parameters for the 2D heterogeneous evaporative model are the same as those for the ID nonadiabatic model (Table 4.3); also, this model requires additional parameters for Darcy's law which can be found in Table 4.4. The boundary conditions used for this model were the same as those used for the homogeneous model. Additionally a constant velocity boundary condition was used at the inlet and a constant pressure boundary condition, was used at the outlet for the momentum balance. A zero flux boundary condition was used on the cell wall. penneabilities penneability. field. oP k [P(l- rjJ) + arjJ]at + V . f.1 V(+ pgz) = Qs [4.5] f3 a r~resent , evapo~ation tenn Qs). 1 D condition. Parameter Symbol Value Units Reference High .permeability k 1 x1 0.5 m2 Schwartz & Zhang, 2003 1 x1 O·g m2 Schwartz & Zhang, Low permeability k 2003 Viscosity of Air J.la 1.8x10·5 I ncropera & DeWitt, Pa s 2002 44 corners that the energy balance predicted local temperature irregularities which did not appear to be physically realistic. This was likely due to a numerical effect of a single mesh element into which flow converges or diverges. No matter how fine the mesh, there would always be a single element in the vertex. This was remedied by "rounding the corners" as is shown in Figure 4.4. After smoothing the mesh in this manner the hot and cold spots disappeared. 4.7 - ID Reactive Model LD geochemical reactions. The ID model was used in order to save computation time and because the homogeneous model adequately described the behavior of the humidity cell during the experimental drying phase. The mass balances for liquid and vapor phase water were the same as those used in the nonadiabatic evaporation model. Additionally, component material balances were added for the ten chemical species considered in the four reactions. Thus, 12 mass balances were performed along with an energy balance. The reaction rates were coupled to the energy balance through and Arhenius rate expression, and an additional heat generation term was included for each reaction (the product of the rate and the heat of reaction). With this geometry, it was observed at the four comers of the inclined layer s~ots lD A lD model was used to model the complete humidity cell cycle including the 10 performed 45 0.06 0.06 0.06 B u E 0.0599 ! 0.0598 0.0598 0.0598 0.0498 0.0498 0.0498 0.0499 0.05 0.05 meters Figure 4.4 A close up of smoothed corner used for the 2D heterogeneous model. The scale is meters. J i 0.0599 / 0.0598 /" in ,. 46 •b,A, R(\-Y2^) [4.6] [OlM]D. ®;ksl[02M] Y [FeS2\ [4.7] bi is the stoichiometric coefficient of the mobile reactant (aqueous oxygen in this case) and As is the specific surface area of the porous medium (given in m2 per m3 or m"1). De ksj 02(aq)] the bulk liquid and 7/ is the ratio of pyrite concentration over its initial value. <t>\ is the fraction of surface area available to react. The second reaction is a neutralization reaction with calcite (CaCCh). CaCOm + 2Hlq) - ^ C a 2 ^ + C02(aq) + H20(l) [4.8] 4.7.1 Reactions The first reaction considered was the oxidation of pyrite by dissolved oxygen (Equation 1.1). The rate of reaction (r/) was modeled by the shrinking core rate expression (Levenspiel, 1972; Fogler, 1999) below (see Appendix A for derivation). r,) y. _ [FeS21 1- [FeSJo b, As m2 m3 mol). R is a representative particle radius and De is the diffusivity of the mobile reactant in the pores of the particle. ks' is the kinetic rate coefficient which has an Arhenius relationship to temperature. [02(aq)] is the aqueous oxygen concentration in Y, fhe 1>1 is CaC03). This reaction is also given a shrinking core rate expression. However, hydrogen ion is now the mobile reactant for this reaction. b2As [4.9] [CaC03] 2 \£aCO\ [4.10] The third reaction is the oxidation of ferrous ion. [4.11] r3 = ~ ¥ l 1 + [4.121 equation 4.12 d is an empirical constant (Lin, 1996). The final reaction is pyrite oxidation by ferric ion. Its rate is also modeled by a shrinking core rate expression; ferric ion is the mobile reactant. Water is assumed to be in excess. FeS2(s) +14Fel;q) + W20{l) - ^ - > 1 SFefo + 2SO^ +16H(aq) [4.13] 47 by - b2 As ~ ~-R-(1-_-y~_Y,T,)~~y~_Y,7,-- -"~~'=-~ + ~~'~~ [H' ]D, <!J ,k,, [H ' ] y _ [CaCO,] , - [CaCO,1 In etrlpirical Lin. 4.12J [4. 13J ,. 48 - ^ M ' „ - „ [4-14] [7V+]De 04 * i 4 [ / V +] Aqueous concentrations also change due to evaporation. The time rate of change in concentration of species i (C,) is given by equation 4.15 (see Appendix C). This was treated as an additional reaction rate. t dS r. . ,., r, = - - - S The overall reaction rate (rr) of a species was given by the following. vr, [4.16] v is the stoichiometric coefficient of the species (positive for products and negative for reactants) and n is the number of reactions in which the species is involved. The ID reactive model uses the same parameters as the ID nonadiabatic model (Table 4.3) with the addition of reaction parameters which are given in Table 4.5. 4.7.2 Weekly cycle In order simulate the weekly cycle, which is characteristic of humidity cell testing, three separate models were developed. One model, the drying phase, had [4.14 ] Ci) Ci r. =--- I P dt [ 4.15] reacti?n CrT) n rT=Lviri i~l v is the stoichiometric coefficient of the species (positive for products and negative for reactants) and n is the number of reactions in which the species is involved. The 10 reactive model uses the same parameters as the 10 nonadiabatic model (Table 4.3) with the addition of reaction parameters which are given in Table 4.5. ,. 49 kH mol kg- 1 bar"1 De 10'1 0 m2s"1 m o l 0 5 m"0 5s"1 ARi 2.6x10"1 0 ARZ mol m"1 s"1 2.3x101 0 AR4 mol m"1 s 1 Eai 5.7x104 mot1 Ea2 3 mol"1 Activation Ea4 9.7x104 mol"1 Reaction 6.3x10"9 mol m'3 s"1 d m"3 Heat AH1 2.8x106 J Heat of reaction 2 AH2 3.5x104 J AH3 9.9x104 mOl(OXygen) 1 AH4 1.8x104 J mOlfovrjte) 1 already been developed (other than the addition of the reacting species). This model simulates a time period of three days. In this model all of the liquid and aqueous components are stationary and have zero flux boundary conditions. The gas phase components, Cfyg) and H 2 0 ( V A P ) , both have the same convective velocity (u). Equation 4.17 accounts for the increase in velocity due to evaporation. " = J ° 1 4-17] uo is the velocity at the inlet and CT is the total concentration of the gas phase, which remains constant and is given by the ideal gas law. Both gas phase components have Table 4.5 - Parameters used in the ID reactive model Parameter Symbol Value Units Reference Henry's law coefficient for NIST, oxygen kH 0.0013 mol kg-1 ba(1 2007 Shrinking core diffusivity De 10-10 m2 S-1 Lin, 1996 Pre-exponential factor, reaction 1 AR1 7.5 molo.5 m-05 S-1 Lin, 1996 Pre-exponential factor, reaction 2 AR2 2.6x10-1O mol m-1 S-1 Lin, 1996 Pre-exponential factor, reaction 4 AR4 2.3x1010 mol m-1 S-1 Lin, 1996 Activation energy, reaction 1 Ea1 5.7x104 J mor1 Lin, 1996 Activation energy, reaction 2 Ea2 8.4x103 J mor1 Lin, 1996 Activati~n energy, reaction 4 Ea4 9.7x104 J mor1 Lin, 1996 Reacti~n 3 rate coefficient kf3 6.3x10-9 mol m-3 5-1 Lin, 1996 Reaction 3 rate parameter 1 mol m-3 Lin, 1996 Hept of reaction 1 ,1H1 -2. 8x106 J mol(pvrite) -1 Lin, 1996 He€lt ,1H2 -3.5x104 J mol(calcite) -1 Lin, 1996 Heat of reaction 3 ,1H3 -9.9x104 J mol(oxygen) -1 Lin, 1996 Heat of reaction 4 ,1H4 -1.8x104 J mol(pvriteL 1 Lin, 1996 mode I 9f 02(g) H20(vap), [ 4.17] Uo CT 50 concentration boundary conditions at the inlet and the convective flux boundary condition at the outlet. All other boundary conditions are the same as those used in the nonadiabatic model. The drying phase is followed by a three day standing phase. To model this, the gas phase velocity was set to zero and all the liquid components were given zero flux boundary conditions. The energy equation was also given zero flux boundary conditions'. should be noted that heat can still leave the system through radial conduction (Equation 4.4). The vapor phase concentrations were held constant at both boundary elements to allow for diffusion from the surrounding atmosphere. The leaching phase proved to be as difficult to model as the dry-air phase, if not more so. This is due to the exceptionally high concentration gradients that arise as the wetting front moves through the medium. This phase was modeled with a constant liquid phase velocity at the inlet for the period of one hour. It is worth noting that the inlet is now the top of the humidity cell. All species are given concentration boundary conditions at the inlet and convective flux boundary conditions at the outlet. The water concentration is set to its initial value (So). This assures the proper initial saturation at the beginning of the drying phase. The energy equation is likewise given a constant temperature boundary condition at the inlet and a convective flux boundary condition at the outlet. The leaching phase was followed by an additional 23 hr standing phase. These four models were developed using the graphical user interface and then saved as script files. The four models were combined into a single script and then set to run in a loop (Appendix C). Each iteration of the loop simulates 1 week. Following the I' i j 1 i I· I It gIven 'fil~s. " 51 completion of each phase, a binary file is created which can be opened in the GUI for visualization of flux effluent. species of interest in the leachate. This is to allow easy comparison with experimental data. 51 1 easy vi sualizat ion or the results. Following each leaching phase simulation, the script file performs a time integration on the nux of the species of interest at the bottom mesh element. This is "used to determine the total amount of each species in the humidity cell emuent. A text file is then created which contains the concentrations of ex perimental , I. ,,. CHAPTER 5 PRESENTATION OF RESULTS AND DISCUSSION Thjs chapter provides the results for the various models discussed in Chapter 4 and, where experimental data are available, compares the simulations with experimental data. The ID isothermal model will be discussed first, followed by the ID adiabatic model and the nonadiabatic model, and so on until the ID reactive model. Finally, the results of a sensitivity analysis will be presented to demonstrate the relative importance of various parameters. 5.1 - ID Isothermal Evaporative Model for a given set of parameters (see Table 4.1). The lumped mass transfer coefficient (kA) is the product of the. mass transfer coefficient and the specific surface area used mass transfer coefficient and surface area are those outlined in section 4.1. Figure 5.1 gives the water saturation profile within the cell at 5xl04 seconds and shows a relatively sharp front with the rock sample essentially drying out completely after the front moves through. For easy comparison with the analytical CHAPTERS Th\s - resul'ts :"-lD Isothermal Evaporative Model Figures 5.1 and 5.2 give the results of the one dimensional isothermal model coefficient to calculate total evaporative flux per unit volume using equation 2.19. The baseline 5xl04 53 Time=50000 Line: S 0.45 0.4 0.35 0.3 B 0.25 .2 - CO GO 0.15 0.1 0.05 0 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Distance from inlet, m axis) after 50,000 seconds as predicted by the ID isothermal model. The vertical line indicates the position of the evaporation front according to the analytic solution. Time-li ne: 0.45 OA 0.35 0.3 c 0 .~ -•"" 10.2 if> 0.15 0.1 0.05 0 0 om 0.02 0.03 0.04 0.05 0 .06 0.07 0.08 Distance from inlet, m Figure 5.1 Saturation (Y ax is) as a function of distance from inlet in meters (X axis) 10 ,. 54 predicted m for £4=100,000 s"1. essentially zero over a period of l.OxlO5 seconds at all three locations, indicating a relatively uniform drying front. The heavy vertical black line indicates the dry-out time according to the analytic solution. 5 0~4: F.:;;)"..'." •.' =.. . =... = .!.': : .:.;; ".~'=..!.= .=.. ~...• ~ ...: :1::.-..'• . -... '. ' ..- .. -... '. .- ...-. ..' . -... -...' ,,'. . -...'. . '==-----;~".';~::i:l 0.35 .~ .. . ....' . . ......• , .. .. 0.3 ... .... ; .. . ....... , ....... : ........ ~..... . ....... ...... " , ............... 0.25 0.2 0 .15 0.1 0.05 .... ........ ., ........... . ...... ,. . ... . . . ... ... : ..... . -....... ~ .... . . . . ... . • . .. . ... .. . .. .. ~ ... " .. : ...... .. ~ .. ... . ....... ) ....... ........ ,. .... .: .... .... : ..... . 1 1.5 2 2.5 3 3.5 4 4.5 5 Time xlO5 s Figure 5.2 Isothermal model pred icted saturation (S) within the humidity cell as a function of time (seconds) at 0.02, 0.05 and 0.08 rn from the air inlet f9r kA= IOO ,OOO 5. 1 , The initial saturation was 0.44 throughout. Note the drop in saturation to Ox 1 05 seco~ds relatively uniform drying front. The heavy vert ical black line indicates the dry-out time according to the analytic solution. ,. 55 solution, a vertical black line has been included in Figure 5.1 to indicate the position of the evaporation front at 5x104 seconds, as predicted by equation 2.36. can be shown, by comparing Figures 5.2 through 5.4, that the evaporation rate is insensitive to the mass transfer coefficient and the specific surface area above a threshold surface value. This indicates that, because the pore spaces are so small, the water vapor is able to come to equilibrium almost instantly. This is to be expected with such a small Peclet Number (0.06). This may not be the case for much larger rock particles or in situations with a much higher convective velocity. Using the parameters displayed in Table 4.1 (the properties of cell F l l ), equation 2.36 predicts that the time required for the evaporation front to move through the humidity cell under isothermal conditions would be 2.9x105 s. Figure 5.2 illustrates the movement of the evaporation front up the 8 cm length of the porous medium. Saturation is shown for three locations in the cell (2, 5, and 8 cm. from the inlet). Figure 5.2 shows model predicted dry-out time to agree perfectly with the analytical solution. 5.2 - Adiabatic Evaporative Model This section illustrates the results of the one dimensional adiabatic model described in Section 4.2. Inclusion of an energy balance greatly affects the dry-out time. Comparison between Figures 5.2 and 5.5 shows that neglecting heat effects will under predict the dry-out time by less than half. Figures 5.5 and 5.6 exhibit an accelerated evaporation front. This was discussed by Stacey and Udell (1997). They '" . 5x 104 It '0.06), 11), 236 cell' 2. 9x I 05 ID ,56 function from the kA= 100,000,000 s'1 . kA=\ 00,000 s"1. 0.45 0.4 0.35 0.3 0.25 V1 0.2 0.15 0.1 0.05 5 \ 1 · · .. ....... :··. ................. : .......i . . . . .. , ! , .. .. ......... : .. .... ~ ..• . •... ~ .......... .......................... . ··· . .... .. .. ... ·i ... .... : ........ ·, ........:.. . . ...... ~... . ...... ............................. . ........· .. .... ........ .···, .... .··v.... . . ·····,.....· · ·· .................. ............ . ........... .... ........·.. ... ..•........... .... .. ........ ........ ... ..... ...... . ....... : ....... ... .. ... : ...... ~ .... .... ~ ···· .... ... .................................... . . . . . . . . ~ . . .. . .. .. , ... ~ ...... ~ ...... .. ~ .......... .. ... .... ..... .... . . .; .... . ... .... , .; ... .... . ~ .. ....... . 1 1.5 2 2.5 3 3.5 4 4.5 5 Time xlO 5 s - 0.02 - 0.05 - 0.08 Figure 5.3 Isothermal model predicted saturation (S) within the humidity cell as a Function of time (seconds) at 0.02, 0.05 and 0.08 m From the. air inlet for IOO,OOO,OOO S·l. Note that the dry-out time is the same as with kA=IOO,OOO S·I . ,. Figure 5.4 Isothermal model predicted saturation (S) within the humidity cell as a function of time (seconds) at 0.02, 0.05 and 0.08 m from the air inlet for kA=[00 s". Note that the dry-out time is the same as with £4=100,000,000 s"1. 57 5 0.4 5 F"'<:=!==~;=!===:-r-'------'---'------'-----' ,=-----0'-"."0""20"" 1 ) ......... .. "'\ .. .... ... ....... ...... .. ... .......... ........... - 0.05 0.4 - 0.08 0.35 ....... : ....... ........ ~. 0.3 ... ). . . .. . ,.. .... . ...... {... .... . .... .. ........ ....... . 0 . 2~ ........·. . ..... ....... ". ............. ~ · . .................................... . V1 0.2 . 0.15 ..... ..·.. . .. .... ......... ...... 0.1 ...... . : ..... ..• , ..... . : .... . . ... ... ~ ... ... . 0.05 ....... ; .... '" .... . .. ;. ....... .; ...... . 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time pred icted cel l lime in let k.-t'" I 00 5. 1 , kA=\ 00,000,000 5. 1 , ,. 58 S Figure 5.5 Adiabatic model predicted saturation (S) in the humidity cell as a function of time at 0.02, 0.05 and 0.08 m from inlet. The time required for the thermal and evaporation waves to propagate through the humidity cell as given by the simultaneous solution of equations 2.36 and 2.37 are indicated by vertical lines (Gold for the thermal wave and black for the evaporation wave). Note the drop in saturation which accompanies the propagation of the thermal wave. s 0.45 I • - 0.02 0.4 - 0.05 - 0.08 0.35 0.3 ... ...... -,. .. .. ' .. 0.25 . , .: .......... , ....... .. ...... ; ..... . . . ,- . ~. . . . . ..... .... ,. .. _ .. .. . · . V1 ;' ···· ..... 0.2 .. .... _··· ....... . .. ...... ..... ... . .... . 0 .15 ..... ...... ......... .. ...................... ........ .. ................ . 0.1 ..... ... ~ .......... ~ ...... .. ~ ......... -: .......... ~ ... ....... .... .. .. 0 .05 . . ... . . . . ... .. .. ~ . .. . . ... . ..•.• . , .. . ; ... . .. ... ; .. ..... .. .; .. °0~~~1 --~2~~3----4~--5~--~6--~7~-U8 Time function which accompanies the propagation of the thermal wave. inlet ID accelerated , . 0.45 0.4 0.35 ;. 0.3 g ·E ; 0.25 .3 • • V> 0.2 0.15 0.1 0 .05 o 59 Time= 15000 Line: 5 o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Distance from inlet, m Figure 5.6 Model predicted saturation in the humidity cell vs. distance from in let at 15,000 seconds as predicted by the I D adiabatic model. Note the ae;ce lerated evaporation front. 60 illustrated how a thermal wave, when it propagates more quickly than the evaporation wave, will be accompanied by an accelerated evaporation front. This is due to phase equilibrium being shifted to a new value corresponding to the temperature difference across the thermal front. The ID adiabatic model gives reasonable agreement with, but fails to match perfectly the analytical solution. The temperature drop predicted by the numerical model (Figure 5.7) is about 2% less that.the 16.8°C predicted from the simultaneous solution of equations 2.36 and 2.37. This is almost certainly due to the fact that thermal conduction was not completely disabled. If the thermal conductivity were to be set to zero, the numerical solution of the energy balance would become unstable. This is why the thermal conductivity was lowered by four orders of magnitude instead of being set to zero. This lowered thermal conductivity will certainly serve to make the total heat transfer due to conduction much less than that of convection. However, wherever there are very high thermal gradients (i.e., on the thermal and evaporation fronts), heat transfer by conduction may become significant. Thus, the numerically predicted temperature drop is of slightly lower magnitude. Also, the simultaneous solution of equations 2.36 and 2.37 predict a dry-out time of 7.93x10s s while the numerical model's prediction is 10% sooner. This is probably caused by the same phenomena, although the magnitude of this discrepancy is significantly higher. This is likely due to the fact that the vapor pressure of water is exponentially related to temperature by either the Antoine or Clausius Clapeyron equations. The Clausius Clapeyron equation, for example, predicts a 20% increase in thennal difference front. . Fi$ure thaUhe , thennal thennal thennal ther~ thennal thennal 7.93 x 105 Time=50000 Line: T-TO 0 d -8 o -a 0 0.08 Figure 5.7 - Model predicted temperature drop after 50,000 s as a function of distance from the inlet. Notice the high temperature gradient at the evaporation front (left side of the curve). The thick horizontal line indicates the temperature drop given from the simultaneous solution of equations 2.36 and 2.37. o -2 -4 -6 "i ci -8 2 " ~10 i'! g E12 ~ -l4 T --16~~~~ o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Distance from inlet, m 61 di stance simultaneous sol ution of ~q uation s 2.36 and 2.37. ,. 62 downstream of the evaporation front). This increase in vapor pressure, will, of course, lead to lower dry-out times. Therefore, heat conduction should have a greater impact on dry-out time than on temperature drop. ID Nonadiabatic Evaporative Model F l l modeled nonadiabatic Fl 1 model. The experimental data from Klinker (2006) have also been shown for comparison. Klinker's data were chosen because he used insulated humidity cells in his experiments in order to minimize the effect of fluctuations in ambient temperature on temperature in the humidity cells. Klinker presented two sets of data for each humidity cell ( F l l and Control) which are referred to as Week One and Week Two, respectively. The following two figures display the results of the ID nonadiabatic model for each humidity cell along with the respective experimental data from Week One. The experimental data exhibits small fluctuations in temperature which are seen on both sets of data (Fll and Control). This is likely due to changes in ambient temperature. vapor pressure from a 1% increase in temperature at 8°C (the predicted temperature 5.3 - lD The humidity cell FII was mode1.ed using the parameters shown in Table 4.3. The ID nOhadiabatic model was also used to model the Control Cell. The parameters I· used to model the Control Cell were the same as used for Cell F 11 with the exception I I'- of the porosity and the initial saturation (Tables 4.1, 4.2, 4.3). Figures 5.8 and 5.9 show the results of the one dimensional non-adiabatic t? Fll F 11 ,- ... cell F l l , Week One. OA, 1A and 2A are the bottom, middle and top RTD probes respectively. Model points 0, 1 and 2 are the temperatures predicted by the ID non­adiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). Initial saturation was set at 44%. Note the spikes and dips in temperature (compare with Figure 5.9). i- ,. OA 4.00 +---1 1A ~ 2A ~ - model-O I: (I) -model-1 ~ ~ 0.00 +---1_ model-2 f-- ---------jj i5 CD 3 Of ~ ~400 I;: (I)' I- 63 -8.00 -1-~---_,-----,_----, o 100000 200000 300000 time, s Figure 5.8 Experimentally determined and model predicted temperatures in humidity FII , lA I 10 nOI1- adiabatic From sec 3. 1). I 64 -10.00 1 1 1 0 100000 200000 300000 t i m e , s Figure 5.9 Experimentally measured and model predicted temperatures in the control cell, Week One. OA, 1A and 2A are the bottom, middle and top RTD probes respectively. Model points 0, 1 and 2 are the temperatures predicted by the 1D non­adiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). Initial saturation was set at 16% in the model. Note that the spikes and dips in temperature appear to coincide with Cell Fl 1 (Figure 5.8). * - , . 2.00 lI:: Q) c" Q) -2.00 ~ ~ c DB Q) 1B ~ -:IIlI , 2B ~ -6.00 Q) -model-O Q. E - model-1 ~ - model-2 -10.00 +------,------,--------, o time, s 1 A bOltom, thc I D nonM ad iabatic 3. 1). Init ia l that' dips in temperature appear to coincide with Ce ll Fll (Figure 5.8). ,. 65 Because of this, it may be necessary to adjust the initial saturation and even sometimes the porosity, in order to match the experimental data. After making minor adjustments to the initial saturation, the nonadiabatic model was able to predict the temperature in the humidity cell with good agreement to the experimentally collected data. ID after adjusting the initial saturation, along with the experimental data collected by Klinker (2Q04) for both humidity cells on Weeks One and Two. The best fit was obtained by raising the initial saturation of the Control Cell from 16% to 19%. The initial saturation of Cell Fl 1 was lowered from 44% to 40%. It was not necessary to adjust the porosity for either of the humidity cells to obtain a good fit with the experimental data. the ID model. This suggests that a ID model is capable of adequately describing the system. FigUre inner edge of the humidity cell to be less than one degree Celsius at a time of 200,000 s after the start of the dry-air phase. Figure 5.15 shows that the time 200,000 s falls It has been observed that the humidity cell does not always retain the same amount of water (Eldredge, 2005). There can also be settling from week to week. Figures 5.10 through 5.13 show the results of the 1D nonadiabatic model, F 11 5.4 - 2D Homogeneous Evaporative Model Figures 5.14 through 5.16 display the results of the two-dimensional homogeneous evaporation model. The two-dimensional model gave results similar to 1D aiD Figure 5.14 shows that the difference in temperature from the center to the 9f ofthe falls ,. -8.00 1 1 1 0 100000 200000 300000 s Figure 5.10 Experimentally determined and model predicted temperatures in humidity cell F l l , Week One. OA, IA and 2A are the bottom, middle and top RTD probes respectively. Model points 0, 1 and 2 are the temperatures predicted by the 1D non­adiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). Initial saturation was set at 40% to match experimental data. , . ~ oj .uc, 4.00 +----1 ~ . OA 1A 2A - model-O -model-1 ~ 000 ;-t------11=~~~1-~-------J, c ~ -.". ,. iii-4.00 . Q. E ~ -8.00 +-----~-----~----___. o time, s 66 5. 10 FII , RTO respect ively. I I D ,. I CD o * O A • 1 A 2 A - modeled 0 - modeled 1 - 2 5 0 0 0 0 1 0 0 0 0 0 1 5 0 0 0 0 s 2 0 0 0 0 0 2 5 0 0 0 0 3 0 0 0 0 0 5 . 1 1 cell F l l , Week Two. OA, 1A and 2A arc the bottom, middle and top RTD probes respectively. Model points 0, I and 2 are the temperatures predicted by the 1 D non­adiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). Initial saturation was set at 40% to match the experimental data. • OA • 1A 4 2A ''"Ii . --momdoedeleledd 01 o 67 c , - modeled L m~ 0 ; -~==========~--------------------------------~I------~., ;;: o ~ ~ 1! 8. -4 E ,! o 50000 100000 150000 Time, s 200000 250000 300000 Figure 5. 11 Experimentally determined and model predicted temperatures in humidity I I, 1 A boltom, respecti ve ly. predi cted I D non­ad iabatic in let 3. 1). saturat ion 4.00 1 1 : 1 ' 0 100000 200000 300000 s Figure 5.12 Experimentally measured and model predicted temperatures in the OA, 1A and 2A are the bottom, middle and top RTD probes respectively. Model points 0, 1 and 2 are the temperatures predicted by the ID non­adiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). Initial saturation was set at 19% to match the experimental data. 68 4.00 -,------------------- OB 1B 2B -12.00 +--------.-----~--~-~ o time, s control cell, Week One. 1 A I 10 3. 1). sat ur~ti on ,. -12 / • 14 • 2B - 1 s Figure 5.13 Experimentally determined and model predicted temperatures in the control cell, Week Two. OA, 1A and 2A are the bottom, middle and top RTD probes respectively. Model points 0, 1 and 2 are the temperatures predicted by the ID non­adiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). Initial saturation was set at 19% to match the experimental data. '.",; u c ~ ~ ~ 'E c ~ ~ ~ 1ii ~ ~ Q. E .~.. 69 4 I , ,~ " . '. f -' o ' . ~ i OB -4 1B 2B - model 0 - ' model 1 - model 2 ·8 ,12 +-------_------_ ------_------_ _ ------_----~ o 50000 100000 150000 time, S 200000 250000 300000 determ ined cel l, I A arc I D non­ad iabatic correspo nding di stances in let 3. 1). 70 Time=2e5 Surface: Temperature [K] Max: 298 Arrow: Velocity field [m/s] «298 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08Min: 293.582 meters field cell homogeneous the radial temperature gradient occurs in the insulating layer. Both x and y axes are given in meters. . . b.08 0.07 0.06 0.05 ;""~; 0.04 E 0.03 0.02 0.01 0 Tlme=Tefl1)erature (298 297.5 297 296.5 296 295.5 295 294.5 294 o 0.Q7 O.08Min: . Figure 5.14 Temperature (color scale) and flow fie ld (arrows) in humidity cel l after 200,000 seconds, as predicted by the 2D homogeneolls model. Note that nearly all of insu lating Bolh , . 71 Figure 5.15 Experimental and model predicted temperature drop in humidity cell Fl I as a function of time. The temperature data comes from Klinker (2006) Week Two. • OA 5 • 1A 2A -modeled 0 - modeled 1 o ,, - modeled 2 , , -10 +-----~------~------~------~----~------~ o 50000 100000 150000 Time, s 200000 250000 300000 Fli , 72 Time=2e5 Surface: Saturation Arrow: Velocity field [m/s] Max: 0.184 I 0 . I 8 |o, 17 0.16 0.15 0.14 ao.12 o n I01 0.09 0 0.08^'^ 0 0 8 3 9 meters 5.16 color flow cell s. gradient saturation absence of any significant radial saturation gradient. Both the x and y axes are given in meters. I 0.08 0.07 0.06 ~ 0.05 8" 0 E 0.04 0.03 0.02 0.01 0 2eS m/s] 0.184 0.13 o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 O.08Min: 0.0839 Figure 5. 16 Saturation (co lor scale) and now field (arrows) in humidity ce ll after 200,000 Note the high grad ient in saturat ion at the evaporation front and the saturat ion in meters. , 73 difference between the temperature in the humidity cell and its surroundings is close to its peak. is ID uniform heterogeneous model described in section 4.5. Figures 5.17 through 5.19 were generated from the results of the "high permeability layer" model while 5.20 through 5.22 represent the "low permeability layer" model. cell than was exhibited by the homogeneous model. There is a difference in cell. This is roughly double that of the homogeneous 2D evaporative model. This is due to the fact that the flow field is no longer uniform in the radial direction (see Figures 5.17 and 5.18). Figure 5.18 shows how the irregular flow field gives rise to an irregular evaporation front, which explains the higher radial temperature gradient. roughly in the middle of the evaporation period. Thus, at this point in time, the Figure 5.16 shows that there is essentially no radial saturation gradient. This IS to be expected. There is, however, a very high saturation gradient at the evaporation front, just as was seen in the 10 model. The unifoml flow field is shown in both Figures 5.14 and 5.16. 5.5 - 2D Heterogeneous Evaporative Model Figures 5.17 through 5.22 illustrate the results of the two-dimensional 5.5.1 High permeability layer Figure 5.17 shows a higher radial temperature gradient within the humidity temperature of nearly two Kelvin from the center to the inner edge of the humidity 20 1.001e5 [K] Velocity [ m / s ] ™298 Figure 5.17 Temperature (color scale) and flow field (arrows) in the 2D humidity cell model with a layer of high permeability material after 1x10s s. Note that the air flow through the high permeability layer is nearly parallel to its orientation. Also note the modest radial temperature gradient in the medium. The x and y axes are given in meters. 0.08 0.07 0.06 0.05 ~ c ;E~; 0.p4 0:93 0.02 0.01 0 Time= l.OOleS Surface: Temperature (K) Arrow: Veloc ity field (m/s] o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 meters 74 Max: 298 Min: 289.379 co lor sca le) 20 permeabi lity I x I 05 The){ y meters. ,. 75 Time=1.001e5 Max: 0.196 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 . 0 8 M i n : 0.0841 meters Figure 5.18 Saturation (color scale) and flow field (arrows) in the humidity cell with a high permeability layer after 1x10s s. Note the irregular evaporation front. The x and y axes are given in meters. 0.08 0.07 0.06 0.05 -"0~ 0 0.04 E 0.03 0.02 0.01 0 Time=> L OOleS Surface: Saturation [moJ/m 3] Arrow: Velocity field [m/s) o O.OS O.OSMin: 0.OS41 meter!> Saturat ion sca le; celt permeabil ity a~er Ix 10 irregu lar evaporat ion ,. Figure 5.19 Temperature drop in the humidity cell with a high permeability layer as predicted by the 2D heterogeneous model. Note the irregularity introduced by the high permeability layer in the model. The time scale is given in seconds. 76 Temperature (K( 298 r---~--~~~~~--~---'r-~~~~ - (0,7.94e- 3) 297 ................. .. 296 295 , ~ 291 ~" 293 v . Co • E 292 ~ 291 290 ... ~ ......... , . ~ ..... . . ... : ... ... ~ ..... . .. ....... ~ ........... : ......... . · . ........ ~.... ........ . . . ... ... .. . ................................ . ..........······ _ . ....... .. -..., .. ..... .... . ..... ......•. ......... '......., ......... . ~ ........... ,.... . .... , ...... .. ... ~ ..... ... ...• ' · . ..... ~ ··· ........... ,...; . 289 2880~0-.5~ -~1- ~~15. -~2-~2.5~ ~3 ..... ; ... ......... ; ....... ..... ~ .. .. Time xlOSs -(0,0.03) - (0,0.06) 20 sca le second,s. ,. 77 1.001e5 K] Max; 0 0.01 0.06 M i n : 288.937 meters flow Field with a layer of low permeability material after 1x10s s. Note that the air flow in the low permeability layer is roughly perpendicular to its orientation. Also note the radial temperature gradient. The x and y axes are both given in meters. 0.08 0 .07 0.06 0.05 ~~" , ~E 0·9' 0.03 0·02 0.01 0 Time=l.DOleS Surface: Temperature [K) Arrow: Veloc ity field [m/s] o 0 .01 0.02 0.03 0.04 0.05 0 .06 0.07 Max ; 298 298 295 294 293 Min: 288.937 Figure 5.20 Temperature (color scale) and Ilow fie ld (arrows) in the humidity cell permeabi lity I x I OS tlow ,. 78 Time=1.001e5 Max: 0.190 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 . 0 8 M i n : 0.0841 meters Figure 5.21 Saturation (color scale) and flow field (arrows) in the humidity cell with a low permeability layer after 1x10s s. Note the irregular evaporation front. The x and y axes are given in meters. 0.08 0:07 0106 0.05 ~, g 0.1)4 ~ E 0.03 0.02 0.01 o l .OOl eS Surface: Saturation (mol/m3] Arrow: Velocity field [m/s] o O.OSMin: 0.OS41 meters now ce ll permeabil ity af!er Ix l05 S. irregu lar Figure 5.22 Temperature in the humidity cell with a low permeability layer as predicted the 2D heterogeneous model. Note the irregularity introduced the low permeability layer in the model. The time scale is given in seconds. 79 Temperature IKI 298'---~---'~-'--~~--~---'r~~~~~ --(O, 7.94e- 3) 297 .. ...... .............. ; ............ , ........ .. .. ;.. ..... -- (0,0.03) 296 ........... :. ·· ..........., .. ............; . ...... ... ....~ .. 295 ........•..........•. .. :;< ~ 29~ :J ~ 293 "a. , E 292 ~ 291 290 289 · . . · ···~···· · ······l···· .......... -.... .. .. .. ... ... ...... . . ....... _·· . .......... .... ...... . .··.. .... .. .... _..· . ...... .. ...... ......... . .. .. .... "· ............. ' .. . .. ... ·, . . ......... ". .............. ......... . · . .... .... ; ............ ~ · ........... ,.. . ........ . .... ) ... ...... ........ .... .. ; 2880~--0~.~5--~1----1~.-5--~2--~2~.5~~3 Time -- (0,0.06) humid ity by 20 by sca te scco nd ~ . 80 slightly different from the homogeneous case. 5.5.2 Low permeability layer would be expected. The air moves through the low permeability layer in a direction temperature gradient within the humidity cell with the low permeability layer is roughly the. same as that of the previous case (see Figure 5.20). Figure 5.21 shows that the evaporation front is also irregular. However, this irregularity is not as pronounced as it was in the previous case (with the high permeability layer). The transient temperature profile in the low permeability layer 2D heterogeneous evaporative model (Figure 5.22) appears to be roughly the same as with the high permeability layer. 5.6 - ID Reactive Model ID FU 5.11 shows that the temperature profile is basically the same as it was without inclusion of heats of reaction. Figure 5.23 shows that, for the case of cell Fl 1, the The time dependent temperature profile in the humidity cell (Figure 5.19) is also The flow fields predicted by the model in the heterogeneous cases are what perpendicular to the interface, while the air moves in a nearly parallel direction through the high permeability layer (Compare Figure 5.17 to 5.20). The radial th~ ID The 1 D reactive humidity cell model was used to better understand the reactions that occur in humidity cell Fil and their rates. Figures 5.23 through 5.27 give the results of the one-dimensional reactive humidity cell model. Experimental data have also been included for comparison. Comparison between Figures 5.23 and shows, F 11, 81 • O A 1 1 1 1 1 1 0 s ID experimental data from Cell F l l , Week Two (Klinker, 2006). OA, 1A and 2A are the bottom, middle and top RTD probes respectively. Model points 0, 1 and 2 are the temperatures predicted by the ID nonadiabatic model at the corresponding distances from the inlet of the humidity cell (see Section 3.1). • OA 5 • 1A 2A -modeled 0 - modeled 1 81 c: 0 i -f-.b- ..m.;o,=de=le=d ;2; L--------------I---- ::! ) c " -10 +-------~--------._------_.--------.---------r_------~ o 50000 100000 150000 Time, S 200000 250000 300000 Figure 5.23 Temperature as predicted by the 10 reactive model along with FII , IA RTO poi nts I 10 nonadiabalic cel l 3. 1). ,. • Measured - Model 10 ' 2 0 Figure 5.24 Experimentally determined and model predicted of leachate. The model appears to predict pi I well over the entire one-year testing period. 4 , ••• • • • •• • •••• • . • • ' a 2 • • i • • • • ••• •• • • • • • • ••• • • 82 .. • • t • o ~,------~--~----~------~------__ ------~-- o " 20 30 40 50 Time, weeks pred icted pH or pH , Figure 5.25 Total concentration of iron in leachate as predicted the ID reactive humidity cell model along with experimental data from Farney (2004). The model appears to adequately predict the concentration of iron in solution. e ~ ~ i' !O. 1000 ,-------- 100 • • • • 10 • • • • • • • • • • • • • • • •• •• • • 83 I • Measuled I-Model .. . .. • • 0 .1~----~-----~-----~----~----~~ o 10 20 30 40 50 Time, ..... eekS pred icted by 10 The sol ution. ,. 1500 T 1000 - E Q-O 500 - - - i - « - 10 20 30 40 50 Figure 5.26 Experimentally determined and model predicted concentration of sulfate concentration. "'" 1000 E ~ ~ 500 .. Weeks • . . 84 • Measured - Model sulrate ion in leachate (ppm). The model appears to slightly under predict sulfate ,. Figure 5.27 Concentration of calcium ion in leachate as predicted the ID reactive humidity cell model along with experimentally determined values from Farney (2004). Calcium ion is over predicted by the model. This is likely due the model's inability to model precipitation. 85 '~r--------------------------------------------r~~~~ I _ Measured l Model I ". '" f~ E • ~ ~ .' '" • !!. •• • • • • • • • • • " • • • ""'+------------r----------__ r-----------,-----------~-----------__ ~ " '" Time. weeks by I D experimenta lly dClcnnined va lues like ly model 's precip itation. three day dry-air phase provides nearly exactly the amount of time needed to completely dry the sample. The model's output is displayed with the Week Two experimental data collected by Klinker (2006) for comparison. The reaction rate coefficients have been adjusted to give a good fit with the experimental data. The values of these parameters are given in Table 4.5. For the purposes of this study, prediction of pH was deemed most important, followed by the concentrations of total iron, sulfate, and calcium, in that order. In Figures 5.24 through 5S7, the leachate concentrations predicted by the ID reactive model are shown along with experimentally determined values by Farney (2004). Figure 5.24 shows excellent agreement with experimentally determined pH values. Iron (Figure 5.25) also appears to be predicted well by the model. Figure 5.26 shows sulfate ion concentration to be slightly under predicted, but is still within the range of values measured. Calcium ion concentration (Figure 5.27) is over predicted. This may be accounted for by precipitation. The model does not take into account mineral precipitation or dissolution. It is likely that at high concentrations, calcium ion will precipitate out of solution. This would be especially prevalent when the sample dries out. It is likely calcium deposits are unable to dissolve quickly during the leaching phase, which would lead to lower than predicted values of calcium concentration in the leachate. 5.7 - Sensitivity Analysis A sensitivity analysis was conducted to determine the relative importance of various model parameters. This was done by changing the values of selected 86 an.d 5,27, 10 a~counted 87 results with those obtain in section 5.6. The results are briefly presented in this section. For more detail please refer to Appendix D. 5.7.1 Specific surface area, As An increase in specific surface area of 50% resulted in an overall increase in rapid 1 sulfate. pH dropped by 6% which led to an increase of calcium ion of 127%. By week 20, however, pH and concentrations of iron and sulfate were equal to that of the baseline case. At the end of the 52-week testing period concentrations of iron, sulfate and calcium were down about 95% from the baseline case, and pH was down 44%. 5.7.2 Rate coefficient kS2 increased by 50% the effect was relatively minor. The greatest effect was that of the leachate. Ca2 + only 16% greater than in the baseline case. The increase in pH was small. Initially it was 2% greater, but this value had dropped below 0.4% by week 52. The concentrations of iron and sulfate in this model were unchanged from the baseline case. This is to be expected. Sulfate ion (for the purposes of the model) is produced parameters by a fixed amount (50%) and running the model again and comparing the reaction rates as more surface area was available for reaction. However, this also led to more ra):>id depletion of pyrite which caused the reaction rates to drop off more rapidly. Week I showed an increase of 62% over the baseline for both iron and an,d ks2 When the kinetic rate coefficient for reaction 2 (Equations 4.8 and 4.9) was concentration of calcium ion in the leachate, Initially it was 37% greater than that of the baseline case. By week 52 of the simulation, however, Ca2 + concentration was 88 only by reaction 1 (Equation 1.1) and is not consumed by any reaction. Reaction 3 (Equation 4.11) might be affected by the depletion of hydrogen ion caused by an increase in kS2, but this will not affect the total iron in solution. 5.7.3 Oxygen adsorption rate kAo2 on aqueous oxygen is always present in its equilibrium concentration. 5.7.4 Solid heat capacity, Cps Raising the heat capacity of the rock should decrease the magnitude of the temperature drop during the dry-air phase of the humidity cell test cycle. This should, in turn, lead to higher rates of reaction and, therefore, higher ion concentrations in the leachate. A 50% increase in heat capacity of the solid phase, however, increased ion concentrations by less than 1% for all ions with the exception of ferric ion. The change in ferric ion concentration was actually negative, but still less than 1% different from the base case. This is because ferric ion concentrations, as predicted by the model, are always much less than one. The ferric ion appears to react as quickly as is formed. ks2, kA02 An increase of 50% in the rate of oxygen adsorption into solution had no impact on: the concentration of any of the species of interest. This suggests that, under the conditions that are present during humidity cell testing, the reactions are not limited by the amount of oxygen available for reaction. Therefore, it is likely that Cps t?e tum, 1 % 1 % it formed. 5.7.5 Activation energy Ea/ The rate of reaction is exponentially dependent on the activation energy. Thus an increase in activation energy of 50% should lead to half an order of magnitude decrease in the reaction rate coefficient. The model predicts that in the activation energy of reaction one of this magnitude would result in a 100% decrease in all ion concentrations except for hydrogen ion, which is present in solution even before the reaction. Thus, as would be expected, the reactions are extremely sensitive to the activation energy. 89 EUI ex ponentially acti vation Thlls magni tude coeffic ient. react ion onc cotlcentrat.ions {rhus, extreme ly sensi tive CONCLUSIONS FUTURE WORK weathering of mine rock piles could be obtained. Specifically, it was desired to develop a model that considered ARD (acid rock drainage) reactions without neglecting temperature effects due to reaction or evaporation. The Comsol Multiphysics™ software was used to create the humidity cell model. It was chosen because of its graphical user interface (GUI) which allows simple problem setup, and visualization of results. It allows the user to easily include additional reactions and to change boundary and initial conditions. Additionally, Comsol allows the user to modify the governing equations when necessary. Finally, the energy and mass balances used to model the humidity cell can easily be extended to modeling more complex geometries (such as a rock pile) in Comsol. analytical solution developed by Stacey and Udell (1997). The one-dimensional non­adiabatic evaporative model is capable of predicting the temperature in the humidity CHAPTER 6 CONCL USIONS AND RECOMMENDATIONS FOR WORK. 6.1 - Conclusions The purpose of this work was to develop a comprehensive model of the processes that occur in humidity cell testing so that a better understanding of the mode!. OUI) Comso!. The one dimensional isothermal model shows excellent agreement with the non-adiabatic mod,e. l 91 cell. has been shown that the evaporation rate is insensitive to the mass transfer coefficient and the specific surface area above a threshold surface area. This indicates that, because the pore spaces are so small, the water vapor is able to come to equilibrium almost instantly. Therefore, for the purposes of modeling evaporation, it is not important to determine the mass transfer coefficient or the surface area of the porous medium with a high degree of accuracy. This is not the case, however, where the kinetics of irreversible reactions is concerned. The composition of the leachate is dependent Jon the rates of the reactions that occur in the humidity cell. Many of these reactions occur on the exposed surface of the mineral and as such are highly dependent on the surface area. 6.2 - Future Work The groundwork has been laid to develop a humidity cell model that includes any and all reactions of interest (Appendix B). Additional reaction rates can easily be added. Additional mass balances can also be added as needed for consideration of additional species. Rates of dissolution and precipitation can also be included in like manner. Currently humidity cell experiments are underway to gather data about various samples taken from the Goat Hill North rock pile in New Mexico. Once collected, these data can give valuable information about the chemistry involved in each of the samples. These data will, therefore, provide a broader basis with which the humidity cell model can be verified. Thermal data taken from heated humidity cells can be used in conjunction with the 2D homogeneous evaporative model to help determine the thermal conductivities and heat capacities of the various rock samples. It kinetic:s dependentkm dependent.i~terest 20 92 Additionally, knowledge of the actual surface area of the samples will help to minimize unknowns. Finally, the insight gained from the humidity cell model can be used to generate a full scale model of the rock pile. Additiona ll y, wi ll he lp ins ight i ,. APPENDIX A DERIVATION OF EXPRESSIONS USED IN THE MODEL ,. 94 l The shrinking core rate expression is a common way to include mass transfer limitations in a kinetic rate law (Levenspiel, 1972). The following shrinking core rate expression is given by Levenspiel (p371), '1 dNE _ bCA dt 1 ' R(R-rc) R2 + - -+ - A . l ] k. rD r2k c s A(aq) + bB(S) -• Products [A.2] where S E X is the total external surface area and N b is the total number of moles of the solid species. R is the characteristic particle radius and rc is the radius of the shrinking core. Equation A.l can be rearranged to give the following. -bA, ra = R{R-rc) i R2 r.D.C, r2kCt [A.3] - 2 3 As is the specific area of the porous medium (given in m /m ). The term ksCA is the rate expression for a first order reaction. This will be replaced with the one half order rate expression which is generally used for pyrite oxidation. It is useful to define a variable Y as the ratio of the concentration of the solid species to its initial concentration. A.I - Shrinking Core Rate Expression P371), [ A.I] A(aq) + hE(s) ---* Products [A.2] Sex NB rc I As is the specific area of the porous medium (given in m2/m\ The term ksCA is the rate expression for a first order reaction. This will be replaced with the one half order rate expression which is generally used for pyrite oxidation. It is useful to define a variable Y as the ratio of the concentration of the solid species to its initial concentration. C B 0 K The concentration of the solid species is direction proportional to the cube of the particle radius. rc±Y^R [A.5] Thus a rate expression can be given in terms of concentrations and constants only. R /3 + A.2 - Concentration Change Due to Evaporation Concentration of aqueous species / is defined as the number of moles of i, divided by the total volume of the liquid (VL). In a porous medium the volume of the liquid is defined as the pore volume multiplied by the saturation. = vs [A.7] 1- 95 [A.4] concentrat ion rate- g iven [A.6] Eva poration i 'VL). mult ipl ied v,. = V, S C, =-Lj- [A.81 1 vps It is desired to find an expression for the change in concentration with respect to time which can be included in the rate term for each aqueous species. [A.„ dt dS dt - ^ = ^ A..0] as vps2 By plugging equation A. 10 into A.9 and substituting A.8 we get an expression for the to C, =V~, S 96 S] dC, = ac, as di as at 00, -N =--' as- v, s' [A.9] [A.IO] I 0 rate of change of concentration of an aqueous species due La evaporation. c, cis r =--- , S dt [A. I I] APPENDIX B ARD REACTIONS TO BE CONSIDERED IN FUTURE WORK • Il.EACTlONS ,. Oxidation of Pyrite by Dissolved Oxygen (Shrinking Core) 2FeS2 ( s ) + 7 0 2 ( a q ) + 2 H 2 0 ( I ) > 2Fe(^q ) + 4 S 0 4 ^ q ) + 4H+a q ) Oxidation of Ferrous Ion (both abiotic and biotic). Biotic f(bacteria 1), see Reaction 7 4Fe2+ +0 +4H+ KF(2) >4Fe3+ + 7H O Equilibrium Reaction of Ferric Hydroxide Oxidation of Pyrite by Ferric Ion (Shrinking Core) FeS2(s) + 1 4 r t f c } + UI20(l) KSFPY >\5Fe^} + 2SO;2 aq) + 1 6 / / ^ Neutralization Reaction of Calcium Carbonate (Shrinking Core) CaCOm +-2Hlq) KSHCA )Ca( 2 ; ) + C02(aq) + H20(l) Precipitation/Dissolution of Gypsum Cal;q) SO;lq) 2H20(I) +Jj&L+CaSO< • 2H20(s) Bacteria Xj growth rate as a function of ferrous, ferric and hydrogen ions d ^ =MU\*[X - * l]* [H+] KKM7A + [H+]j [Fe2+]*Mf Fe2+] K S l l K I l * [ Fe KDx,*[Xxf= [mg-rlh->] Neutralization Reaction of Potassium Feldspar (Shrinking Core) KAlSi3Os(s) + H+ (aq) + Ui20 KSHAW > i ^ 2 < 3 5 (Of7)4 + ^ + 10 Formation of Jarosite-K 3 ^ + 2 5 0 ; ^ + K+ + 6//2O( 0 *F ( 1 0 ) ) KFe3S04)2(OH)6(s)+6h 20 Equilibrium Reaction of Geothite FeOOH,., + 3HL, < ^ ( 2 0 ) > Fe,3l, + 2/7,0,, KF(20)/KE(20) 98 ID Reaction Description 1 2FeS2(s) + 702(aq) + 2H2°(l) KSOPY 2F 2+ 4S0-2 4H+ ~ e(aq) + 4(aq) + (aq) 2 = f(4Fe~a~) + 02(aq) + 4H(:q) KF(2) >4 F,e (3a+q ) + 2H2 0 (I) 3 , Fe(OH)3(S)'+ 3H(:q) ( KF(3) ) Fet~) + 3H20(l) t KF(3)/ KE(3) 4 , FeS2(s) + 14F~t:q) + 8H20(l) KSFFY ~ 15Fe~a~) + 2S0;(2aq) + 16H(:q) 5 , CaC03(s) +2H(:q) KSHCA ) C a(2a+q ) +CO2 (aq) + H 2 0 (I) 6 PrecipitationlDissolution Ca!;) + SO;(2aq) + 2H20(l) ( KF(6) ) CaS04 .2H2O(,) KF(6)/ KE(6) . 7 X 1 IOns d[X,] = MUI *[X 1*( [W1 ]- dt I KM7A+[H+] [ [Fe" 1 * M," [Fe2+] + KSI * ( 1 + Kll *[Fe - KDxI *[X 1]2 = [mg·rl·h-I ] ~ 8 KAISi30 8(s) + H(:q) + 2. H2O KSHAW ) l Al2Si20 s OH)4 K(:q) + 2 2 larosite-3Fe~:q) + 2S0~2aq) 6H2O(I) KF(lO) > KFe (SO) (OH) +6!J 3 4 2 6(s) FeOOH + 3H+ ( KF(20) ~ Fe3+ + 2H 0 (s) ,. (aq) KF(20)/Kr,(20) (aq) 2 (I) Hydration of Carbon dioxide C02{aq) + H20(l) < ^ ( 3 8 ) > H2CO, Equilibrium Distribution of Carbonate HCO~ < KF(39) > H+ +CO~2 Equilibrium Distribution of Bicarbonate H2C03 < K F ( 4 2 ) > HCO', + H+ KF(42)/A-£(42) K q> Dissolution/Removal of Carbon dioxide Dissolution/Precipitation of Fluorite Reaction of Epidote (Shrinking Core) CaFeAl, (Si04\ (OH)(s) H+ SH20 KSH > 2 + 2 + 3H4(Si04) Reaction of Anorthite ( Ca-plagioclase) (Shrinking Core) Ca Al2 Si2Om 7 Hlt) KSHAW ^ / + 3 Ca£, + tf^O" Reaction of Chlorite (Shrinking Core) M g ^ A l ^ O ^ O ^ + lSH^ lMgl+ aq) + 2Al+3 + 2 + Neutralization Reaction of Sodium Feldspar - Albite (Shrinking Core) NaAlSi3Om Hlaq) ^H20 KSHAW ) i A * ' 2 0 5 (OH), - Reaction of Illite (Shrinking Core) KAl2 (Si3Al)Ol0 (OH)2 10 Hlq) Klq) Al+2 # 4 S / 0 4 99 38 CO + H 0 ( KF(38) > H CO 2(aq) 2 (I) KF(38)/ KE(38) 2 3 39 HCO- ( KF(39) ) H+ + CO-2 3(aq) KF(39)/ KE(39) (aq) 3(aq) 42 H CO ( KF(42) > HCO- + H+ 2 3 KF(42)/KE(42) 3(aq) 45 DissolutionlRemoval , CO2(g) ( l KLAI >C O2(aq) KIA I , DissolutionlPrecipitation 49 ~ CaF2(s) § ~~li] Ca~) + 2 F-(aq) 50 Ca FeAl3 4 )3 OH\s) + 4 + 8H2O KSH > Fe+2 + Ca+ 2 + 3H4(Si04) Ca-plagioclase) 51 Alz SiZ0 8(s) + H(:q) KSHAW ) 2 Ar3 + Ca(;q) H7Si20~ 52 Mg3 Fe2AlzSi3 O\O(OH)8(S) 16H(:q) ~ 3Mg (az+q ) + 2At3 + 2Fe+ z + ". 53 NaAISi30 8(s) + H(aq) + 2. H2O KSHAW 1 . > - Al Na(aq) 2Sz20 s OH)4 + 2 2 54 KAt2 Si3At)OIO 2 + 1 0 H(:q) KS > K(:q) + 3 AI+3 + 3 H4Si04 100 of Smectite 55 Ca633[Mg0M A^KSOO^OHX^+llH^ + MTaO - S - 0.66 Mg2*q) -> 0.33Ca 3.34 At Reaction ofSmeclite (Shrinking Core) S5 Cao.ll [Mgo.66 All.34 ] (Si8)qO (OH )~(I)+ 1 2H(:tq)+ 8H20 KS) O.33Ca Mg~) + A/+ i APPENDIX C ID i APPENDIXC SCRIPT FILE USED TO RUN THE 1D REACTIVE MODEL ,. 102 % COMSOL Multiphysics Model M-file Generated by COMSOL 3.3 (COMSOL 3.3.0.405, $Date: 2006/08/31 18:03:47 $) % By Paul Evans tic-clear all clc flclear fern Start timer to track computation time version clear vrsn vrsn.name - 'COMSOL 3.3'; vrsn.ext = ' ' ; vrsn.maj or = vrsn.build = 405; vrsn.res ~ '$Name: $'; vrsn.date = '$Date: 2006/08/31 18:03:47 $', fern.version, = vrsn; Geometry\ solidl ['o, 0.07] ) ; Input height of rock in humidity cell % Analyzed geometry clear s s . ob j s= { gl),; s . name- { 111/ } ; s.tags=(1gl'}; fem.draw=struct( ' s ', s) ; fern.geom=geomcsg(fern) , weeks = 52; time = clock; mon = num2str(time(2)); day = num2str(time(3)); hr = num2str(time(4)); min = num2str(time(5)); filename = strcat(1concData out = fopen(filename, 'w+') fprintf(out,'Week H20,L Sulfate Ca2+\n'); Input simulation time, number of weeks Create file to record leachate concentrations ', day, hr, '_', min), pH Fe2 + Fe3 + total Fe % Constants fern, const = 1rhoRock' 1 CpWater1 1CpAir',' 'CpRock', 'kAirO',' 'kRock',' 'kWaterO' 'Hvap', fi','.36' ,'2700 , '4179 '1007 ' , ,'858', '26.3e-3' '2.5', ','.6131 (-241.826+285. Input model parameters 83)*1000' 'CO','cwlSat*S0 [mol/m3]', ... 'SO','.4', ... 'cwlSat','rhoWater/MWwater [mol/m3]', 'rhoWater', ' 998 1 , ... 'MWwater','18.02/1000', ... 'CvapO1,'PsatO/(R*T0)', ... 'TO', '298', ... 1Tboil1,'373', ... 'R', '8.314' , ... 'PsatO','P0*exp(Hvap/R*(1/Tboil-1/T0)) 'P0','101325', ... 'v_in', 1q/(Ac*effA)', ... 'q','1/60/1000', ... 'Ac', ' D"2/4*pi', ... 'D', '.1001' , ... q_in','v_rhoAir0*CpAir*T0', ... 'rhoAirO','P0/(R*T0)*MWair', ... file % COMSOL tic; clear fem % COMSOL version name ~ COMSOL major 0; 405; rcs = Name: ~ Date: $'; fem.versioq vrsni % Geometry: gl=solid1 ( [:0,0.07] ); 4.-1---'---------IL-_________________ ---' ,~ objs={gl).; name= ( 'I 1( ) ; s . tags= ( 'gl ' ) ; draw=struct('s',s) ; fem.geom=fem) ; ~---- 52; time(time(time(time(lconcData_', on, , I, , " I If min); w+'); out, 'n'); fern. canst = ('fi',' .36', 'rhoRock','2700', .. . Fe2+ Fe3+ 'CpWater','4179', .. . '1007', .. . CpRock','.. . kAirO', 26.3e-3', ... , kRock' , , 2 . 5 I, ... ....- ----11 kWaterO','.613', ... Hvap',' 241.826+285.83) *1000', CO', 'cwlSat*SO I SO I, I .4', cwlSat', 'rhoWater/MWwater mol/rhoWater','998', MWwater',' CvapO', 'R*TO)', TO ' ,'298', 'Tboil', '.. . R','S.314', .. . PsatO', 'PO*exp (R* (l/l/TO)) " PO', '101325', v 'q!(Ac*effA) " q', '1/60/Ac', 'D A pi', 0',' .1001', 'q in', 'v in*rhoAirO*CpAir*TO', rhoAirO'-;'PO/{R*TO) *MWair' , •MWair\ ' (.21*32+.79*28.01)/1000', . . . 'kGlass",'.2', ... •kins','.04*2', ... 'U', '2/D/(log(l/effD)/kRock+log(.056/.05)/kGlass+log(.076/.056)/kIns) ', ... ' alpha','.003747/1000', ... •n', '0.166699', ... •ra','4.845923', ... • Sr",'.05', ... 'As','3*3A.5/(4*dp)*400', ... 'dp','5e-4', ... 'Dvap','1.013e- 2*T0"1.7 5* ( (MWwater+MWair) / (MWwater *MWair) )"0.5/(P0*(19.7"(l/3)+13.1'-(l/3))"2)', ... 'Sc','mewAir/(Dvap'rhoAirO)', ... 'Re','v_in*dp*rhoAir0/(mewAir*(1-fi)'shape)', ... 'kMass','Dvap/dp*(1-fi)*shape*2.19*(Re*Sc)"(1/3) ', ... 'shape *,'.95', ... 'KLA','As/400*kMass', ... 'mewAir','1.789e-5', ... 'PsiO','1/alpha*(SeO"(-1/m)-1)*(1/n)', ... 'SeO',1(SO-Sr)/(1-Sr)', ... •effD','*', ... 'effA','e.ffD"2', ... 'bl','2/7'', ... ' kH','.0013*rhoWater/10"5', ... ' 02gO','.21*P0/(R*T0)', ... ' FeS20', '-.06*rhoRock/MWFeS2* (1-fi) ', ... 'MWFeS2'/'(26+2*16)/1000", ... 'Arl','10"(-0.6245)*1000A.5', ... 'Eal','lO"(7.7 559)/1000', ... •De','le-6/100"2', ... 'O2aq0','O2g0*R*T0*kH', ... 'DH1','-1.409e9/1000*2', ... 'Vwater','-0.5e-3/(Ac*fi-SO)/3600', ... •Dliq','le-9', ... •CaCO30','.05*rhoRock/MWCaCO3*(1-fi)', ... 'MWCaC03','100.09/1000', ... 'b5','l/2', ... 'Ar5','10A(-2.5914)/