«a: 1 Its... ”N7... z. .9... ya}... . 2%.: 3. Ian»... 5"! 5 an} .si...h..«th4?! Tie!» .«Nms... 5 £4»... £9 1 fiwfi flifila. 3LT X .t «ow Lnnuflvmh. «Hunwittf I... (Rikki . 11.... . .1918 .5..1:\(t«¢1 I. i... 3.. 7‘... .3 ; .. . . V . ‘ . P. inn}: 5 ”Ewfifiu ‘ .. . .m . ejflfifi 2 3 003 ii LIBRARY Michigan Stat University This is to certify that the thesis entitled Stream Transient Storage Modeling Based On Fractional-ln-Space Dispersion presented by Rammesh Padmanabhan Navaneethakrishnan has been accepted towards fulfillment of the requirements for the Master of degree in Civil and Environmental Science Engineering WW lMajor Professor’s Signature \zi \Li) 2007 Date MSU is an affirmative-action, equal-opportunity employer —-.-n-n--A-A-I-I-I-.-I-I-I-o-n-- -A- -.-.--I-I-I-l-I-I-I-0-O- 1-0-0-I-h-a-n-o-- —.-.-.-.--u-n-l-c-n-o-o-- PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/07 p:IClRC/DateDue.indd—p.1 STREAM TRANSIENT STORAGE MODELING BASED ON FRACTIONAL-IN-SPACE DISPERSION By Rammesh Padmanabhan N avaneethakrishnan A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Civil and Environmental Engineering 2007 ABSTRACT STREAM TRANSIENT STORAGE MODELING BASED ON FRACTION AL-IN-SPACE DISPERSION By Rammesh Padmanabhan Navaneethakrishnan Solute transport in natural streams is often described using the transient storage (TS) model in which two separate equations are used to describe the partitioning of the so- lute between the main channel and the storage zones. Several stream tracer studies have found heavy-tailed rising and falling limbs in the time—concentration breakthrough curves (BTCs) because of the long-range spatial correlation of the dispersion process as well as the exchange processes between the storage zones and the main channel respectively. Application of the TS model to field sites with multiple stream reaches is accomplished by estimating one set of parameters for each stream reach in order to describe the scale- dependent dispersion process. Dispersion and the storage zone processes are known to give rise to competing parameters in the TS model sometimes leading to singular convergence during parameter estimation. Therefore, the inability to describe dispersion accurately introduces uncertainty in the storage zone processes and parameters. In this work, we propose an alternative TS model based on the concept of anomalous diffusion described mathematically using fractional-in-space derivatives. The new Fractional-in-space Tran- sient Storage (FSTS) model has the ability to describe the observed early (non-Fickian) dispersion in some stream reaches and is able to describe the scale-dependence of the dis- persion process. The model is used to describe tracer transport in two Michigan streams using one dispersion coefficient but different fractional derivative exponents in different reaches. The new FSTS model was found to better constrain the TS model parameters by providing a superior description of dispersion. © 2007 Rammesh Padmanabhan Navaneethakrishnan All Rights Reserved Dedicated to my parents and my brothers Suresh and Srikantesh iv ACKNOWLEDGMENTS I would like to thank my research advisor and now a friend, Mantha Phanikumar, for his constant guidance, encouragement and support and for the numerous discussions we’ve had during the course of this research work. He had dedicated more time, effort and energy to this work than I could have asked for. I thank him for always being willing to meet me whenever I barged into his office and also answering to my never ending emails and phone calls, even during weekends. Thank you for always being there. I would like to thank my committee member, Shu-Guang Li, whose unconventional teaching style with emphasis on practical application of theory, encouraged me to learn and explore more about this field over a period of two years at MSU. His constructive feedback and comments have been significantly useful in shaping this work to completion. I would also like to thank my other committee members Roger Wallace and Irene Xago- raraki for their useful comments and suggestions, which helped to bring out clarity and focus to this work towards the end. Financial support from the following two funding agencies is gratefully acknowledged: (a) NSF Biocomplexity in the Environment: Coupled Biogeochemical Cycles (b) NOAA Center of Excellence for Great Lakes and Human Health I would also like to thank a number of other peOple in my everyday circle of colleagues who had enriched my professional life at MSU in various ways. First, my research group colleague and friend Chaopeng Shen for his inputs and companionship in RCE. Heather Stevens, for bugging her many a times with my never ending stories she used to listen to, when I used to get bored with programming. Finally, my roommate and friend, Sagar Rayepalli, for the same food he cooked and satisfied my unsatiable appetite during the course of this work! In the end, I would like to thank my parents for supporting my decision to go for graduate studies and not getting mad at me for not visiting them even once after I came to US. And obviously my two big bros, Suresh and Srikantesh for being a constant inspiration throughout my life. TABLE OF CONTENTS vi LIST OF TABLES ................................ viii LIST OF FIGURES ............................... ix 1 Introduction .................................. 1 1.1 Fractional Calculus and the Fractional Derivative .............. 6 1.2 Numerical Methods ............................... 8 1.3 Outline of thesis ................................. 11 Advection Schemes .............................. 13 2.1 The Weighted Essentially Non-Oscillatory (WENO) Scheme ................................ 15 2.1.1 One-Dimensional Reconstruction ................... 16 2.2 A Fourth-Order Compact Scheme with Spectral-Like Resolution ...... 20 2.3 A Test Case: Advection of a Gaussian Pulse ................. 22 2.4 FFT Analysis of Advection Schemes ...................... 24 2.5 Conclusion .................................... 26 Fractional Diffusion: Preliminaries and Numerical Approximations 28 3.1 Introduction ................................... 28 3.2 Fractional-in-space Diffusion .......................... 29 3.3 Fractional-in-time Derivative: Sub-Diffusion ................. 30 3.4 Fi-actional Derivative: Multiple Definitions .................. 31 3.4.1 The Riemann-Liouville Derivative ................... 33 3.4.2 The Caputo Derivative ......................... 34 3.5 Numerical Approximations: The GL and Caputo Derivatives ........ 35 3.5.1 The Caputo Derivative ......................... 36 3.6 Higher Order Approximations ......................... 37 3.6.1 Higher Order Caputo Fractional Derivative Using Integer Derivatives 39 3.6.2 A Test Example ............................. 40 3.7 An Exact Solution: The Method of Manufactured Solutions ........ 41 Fractional Advection Dispersion Equation ................ 48 4.1 Semi-Analytical Solution for the fADE .................... 48 4.1.1 Numerical methods ........................... 49 4.2 Operator Splitting Methods .......................... 51 4.3 Boundary Conditions .............................. 52 4.4 Solution of the fADE for Continuous Release ................. 54 4.5 Instantaneous slug release ........................... 55 4.5.1 High Peclet Number: Advection-Dominated Systems ........ 56 4.5.2 Low Peclet Number: Dispersion-Dominated Systems ........ 58 4.5.3 Effect of fractional derivative exponent(a) .............. 60 5 Stream Transient Storage Modeling .................... 62 5.1 Transient Storage Model ............................ 64 5.2 The Fractional-in-Space Transient Storage (FSTS) Model .......... 65 5.2.1 Comparison of the FSTS model with the Analytical Solution for a = 2 .................................. 68 5.3 Application of the FSTS Model to Describe the Field Data ......... 69 5.3.1 The Red Cedar River, Michigan .................... 70 5.3.2 The Grand River, Michigan ...................... 73 6 Conclusions .................................. 80 APPENDICES .................................. 82 A Caputo Derivative Operator ........................ 82 B Griinwald—Letnikov Derivative Operator ................. 85 C Griinwald-Letnikov and Caputo Weights ................. 87 Bibliography ................................... 96 vii 2.1 2.2 3.1 3.2 3.3 5.1 5.2 LIST OF TABLES Coefficients of a three-stage—third—order Runge-Kutta scheme, from Lowery and Reynolds (1986) .............................. Order of accuracy for the WENO—Roe and the compact schemes ...... Table showing the order of numerical approximation for Griinwald-Letnikov and h-Caputo derivative operators ....................... Table showing the weights for the Griinwald—Letnikov derivative for differ- ent a values ................................... Table showing the weights for the Caputo derivative for different a values . Storage zone parameters for the Red Cedar River using FSTS model . . . . Storage Zone Parameters for the Grand River for the three reaches with a constant D using the FSTS model ....................... viii 22 22 41 46 47 78 79 1.1 1.2 1.3 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 3.6 LIST OF FIGURES Illustration of transient storage mechanisms (a) when solutes enter small pockets of slow-moving water and (b) when solutes leave the main channel and enter the porous media that makes up the bed and banks of the channel. Arrows denote the solute movement between the main channel and the transient storage zone. [Runkel, 1998; Schmid, 2004] ............. Comparison of integer and fractional derivative of f (as) = 3:2 ........ Fractional derivative for f (:13) = x3 ...................... Comparison of fifth order WENO-Roe scheme with analytical solution for A2: = 0.1, CFL = 0.1 at t =1 ......................... Comparison of the compact, WENO and first-order upwind scheme for Ax=1,u=1,CFL=0.1att=5 ...................... Plot showing that comparison of RMSE error (log-scale) for compact, WENO and first-order scheme for different CFL and Art: values ...... Plot showing that comparison of the compact and the WEN O scheme for An: = 1, u = 1, and CFL = 0.1 with the analytical solution at t = 250. WENO scheme disperses while compact scheme oscillates .......... Plot showing that comparison of compact and WENO scheme for Ax = 1, u = 1, and CFL = 0.1 with the analytical solution at t = 250. Clearly the amplitude of the WENO scheme is lower compared to compact scheme for a crude grid size ................................ The random jump in space and time in a CTRW model ........... The left and right RL derivatives at a point x ................ Control volume approximation to the Caputo fractional derivative ..... Higher-order approximations: The generating function w as a function of 2 Comparison of RMSE error (log scale) for the fractional derivative with the analytical solution between GL and H-Caputo derivative (higher order Caputo) ..................................... Plot showing the variation of RMSE (log scale) with fractional derivative exponents (a) for two grid sizes A1: = 0.1 and Art: = 0.01 with D = 1 at time, t = 1 ................................... ix 4 19 23 25 29 33 37 39 41 3.7 4.1 4.2 4.3 4.4 4.5 4.6 4.7 5.1 5.2 Plot showing the variation of RMSE (log scale) with dispersion coefficients, D, for grid sizes of Ax = 0.1 and A2: = 0.01 with a = 1.8 and at time, t = 1 45 Breakthrough curves shown at a: = 10m from the point of release. u=ImS—1 and D = 1 ............................. 50 Concentration profiles for the continuous release of a tracer in a stream simulated using the fADE. The fADE was numerically solved using the WENO-GL and WENO-Caputo schemes and compared with the semi- analytical solution of [Benson et al., 2000]. The parameters in the solu- tion are Co=1ppb at t= 53, 103 and 15s with u = lms—l, D = 0.1m1°8/s, a = 1.8, CFL = 0.1 and Ax=0.1250m .................... 54 Plot showing the comparisons for all the four schemes for a P6 = 6.6 (Ax=0.25m, D = 0.05m1'8/s, u = nus—1, CFL =01, t=15s and a = 1.8). Numerical oscillations were observed in case of compact based approaches while numerical dispersion was observed in the case of WENO based ad- vection schemes ................................ 57 Plot showing the comparison for a Pa = 1.7 for all the four schemes (Az=0.25m, D = 0.05m1-8/s,u = inns-1 , CFL =01, t=15s, a=1.8). All the schemes show reasonably good fit, although the WENO based advec- tion scheme introduced small numerical dispersion ............. 58 Breakthrough curves for a low P6 = 0.057 (Ax = 0.5, CFL = 0.1, a = 1.8, fl = 1, D = 10m1'8/s, u = 1ms‘1) at a distance of 5m from the spill location. .................................... 59 Breakthrough curve (without OS) for a low Pe = 0.057 (A2: = 0.5m, CFL = 0.1, a = 1.8, ,6 = 1, D = 0.05m1'8/s,u = 1m/s) at a distance of 5m from the spill location. .......................... 60 Breakthrough curves comparison between GL and Caputo with Benson et al. [2000] for different values of a with Pa = 0.1 .............. 61 Numerical model showing the snapshot of concentration vs distance in the main stream at t = 6005 spatially for A = 40 m2, As = 4 m2, D = 1m1'8/s, u = 1ms‘1, e = 0.0 s—l, M = 1000g. The point of injection is at a: = 10 m 66 Breakthrough curve for the stream and storage zone concentration at a distance of x = 50m from the point of injection in (a) linear scale and (b) log log scale. In log scale, it can be observed that for the falling limb the major contribution of the concentration comes from the storage zone whereas, for the rising limb it is the main channel. (A = 40 m2, A3 = 4 m2, D =1m1'8/s,fi = 1, u = 0.5m,-1, c = 0.001 s—l, M =-. 1000 g, a = 2) 67 5.3 FSTS vs TS model breakthrough curve in the main channel for A = 40 m2, A3 = 4 m2, D = 1m1-8/s, e = 0.0001 ms-l, M = 1000 g, a = 2.0, A2: = 2.9950 In and a: = 1000 m ........................ 68 5.4 Comparison of compact-Caputo based fADE with the analytical solution given by De Smedt et 31. [2005] for Q = 16 m3s-1, A = 40 m2, AS = 4 m2, D =1 mas—1, M =1000 g, a = 2.0, A1: = 0.33 m and a: = 1230 m . 70 5.5 Comparison of compact-Caputo based fADE with the analytical solution given by De Smedt et al. [2005] for Q = 20 m3s-1, A = 40 m2, A; = 4 m2, 6 = 0.0001 mS_1, M = 100 g, a = 2.0, An: = 0.33 m and :1: =10 m . . 71 5.6 Comparison of compact-Caputo based fADE with the analytical solution given by De Smedt et al. [2005] for A = 40 m2, A3 = 4 m2, D = 1 mas—1, €=0.0001 ms‘1,M=100g,a=2.0, Az=0.33mand:r=10m .... 72 5.7 Probability plot showing the normalized concentration (or probability) as a function of the Z-scores of travel time for the Reach C of the Red Cedar River. Deviation from the straight line indicates a deviation from Fickian diffusion ..................................... 73 5.8 Red Cedar River Watershed showing the Red Cedar River and the sampling locations (adapted from Phanikumar et al. [2007]) ............. 74 5.9 Fluorescein breakthrough curves for the three reaches A, B and C for the Red Cedar River with a constant dispersion coefficient D = 4.14 using the FSTS model .................................. 75 5.10 Grand River Watershed and study region showing the sampling sites (adapted from Shen et al. [2007]) ....................... 76 5.11 Rhodamine WT breakthrough curves for the three reaches B, C and D of the Grand River with a constant dispersion coefficient D = 3.32 using the FSTS model (data taken from Shen et al. [2007 , in review]) ........ 77 xi CHAPTER 1 Introduction More than half of the world’s population lives in urban areas [Cohen, 2003] and stream ecosystems are undergoing significant changes due to rapid urbanization in many parts of the world. Changes in impervious cover and the concomitant changes in stream hydrol- ogy, geomorphology, water quality and physical habitat are a direct result of urbanization. Anthropogenic activities are the major cause of large increases in stream nutrient con- centrations measured worldwide, often resulting in declining water quality (often referred to as the “urban stream syndrome” [Meyer et al., 2005]). These changes have important implications for human health and the ecological integrity of ecosystems. Understanding how streams respond to natural and anthropogenic stressors is therefore critical to man- aging these resources effectively. Hydrological metrics are generally used as indicators of urban impacts on stream ecosystems [Roy et al., 2005; Booth, 2005]. These include, for example, a flashier hydrograph (i.e., increased frequency of high flows as defined by the stream flashiness index), higher nutrient concentrations but less efficient nutrient up- take rates and less taxonomic diversity or species richness when compared to unaltered streams. The dynamics of solute transport in steams plays a critical role because many solutes such as Nitrogen or Phosphorus are present in short supply and therefore regulate primary or secondary productivity. Solute transport modeling is important since it allows us to explore the relationship between physical characteristics of the stream (such as geo- morphology, nature of the substrate, flow etc.) and hydraulic variables such as dispersion, transient storage, and lateral inflow. A combination of tracer studies and solute transport modeling provides a framework to quantify important processes, separate physical and biological controls and allows comparisons to be made across different scales [Dangelo et al., 1993]. Transient storage (hereinafter TS) refers to multiple processes that con- tribute to solute retention in a stream resulting in a delay in the downstream movement of solute mass. A distinction is usually made in the literature between in~channel storage due to surface features such as vegetation, woody debris and eddies and pools (e.g., due to extensive meandering) within the stream (called surface storage) and the retention of solutes due to interaction with near-bed sediments (called hyporheic exchange). Unless otherwise specified TS usually refers to the combined effect of both types of storage. In stream solute transport models, the dispersion coeflicient D, the size of transient storage zones A S and the exchange rate a between the main channel and the storage zones are im- portant hydraulic variables that depend on stream physical characteristics. The variables D and a depend on pr0perties such as discharge, velocity distribution and the channel cross-sectional area, therefore they correlate with stream order. A 3 strongly depends on the presence or absence of in-stream complexity (e.g., the presence or absence of leaves, cobbles/ boulders). The lateral inflow q L which reflects the contributions from the sides as well as from groundwater depends on the geomorphology, streambank porosity and the location of the watertable [Dangelo et al., 1993]. Measuring stream health is one of the approaches to monitor and direct the sustainable use of streams and stream restoration efforts [Harbott et al., 2005]. Regulatory guidelines require that certain criteria be satisfied for water bodies to be considered healthy. In addition to the usual water quality parameters (such as suspended solids concentration, pH, temperature and dissolved oxygen), recent efforts have included stream organisms as indicators of stream health (e.g., fish and plant assemblages). Microbial communities respond rapidly to environmental change. Two important services are provided by mi— crobial communities - they are a food source for higher trophic levels in the food web and they regenerate nutrients through mineralization of organic detritus [Harbott et al., 2005; Allan, 1995]. Coliforms and nitrifying bacterial counts are typically higher in urban streams (due to sewage discharges) compared to native/unaltered streams. In contrast, urban streams are known to have decreased microbial metabolism due to increased sedi- ment metal concentrations associated with urbanization [Wei and Morrison, 1992]. One very important element of stream restoration, which had long been neglected and is still an area of active research is the hyporheic exchange zone. It is term derived from Greek roots -hypo, meaning under or beneath, and rheos, meaning a stream (rheo means to flow). Hyporheic exchange zone provides habitat for aquatic macro-invertebrates and provides conditions for heterotrophic primary production [Mulholland, 1997]. It is responsible for changes in surface water quality through physical and biogeochemical processes. There are number of variations to the definitions of hyporheic zone (Figure 1) depending upon the field they are used in as shown by the definitions given below: 1. The hyporheic zone is the zone below and adjacent to a streambed in which water from open channel exchanges with interstitial water in the bed sediments; 2. It is the zone around a stream in which the fauna of the hyporheic zone (the hy- porheous) are distributed and live; 3. It is the zone in which groundwater and surface water mix. A large number of studies have been reported in the last decade to understand the exchange between the hyporheic zone and main stream including physical mixing with groundwater [Constantz, 1998], chemical reactions [Bencala and Walters, 1983; Duff et al., 1998], microbially mediated chemical transformations [Duff and Triska, 1990] and trans- port of nutrients thorough groundwater-surface water interactions [Wondzell and Swan- son, 1996; Valett et al., 1996; Jonsson et al., 2003]. Stream solute transport modeling is an important activity in all of the efforts reported above as it provides a means to test hypotheses about processes based on observed data. The focus of this research is on the numerical solution of the TS equations which govern solute transport in streams, therefore we will study the hyporheic zone from a hydrologist’s angle, which often regards it as an extension of the stream channel, where water can spend longer duration than in the open channel [Jones and Holmes, 1996]. Conceptual models of the hyporheic zone and the processes that occur within it have been studied in many fields, but we will limit our study to the well-known TS model as presented, for example, by Bencala and Walters [1983]. Pool and Riffle Stream Meandering Stream . .3; — / \ I" «7’ Direction of \1’3 r n . l c :3’ Ground-we er Flow Influence of local and regional grou nd-water flow systems, hyporheic zone, and stream Figure 1.1. Illustration of transient storage mechanisms (a) when solutes enter small pockets of slow-moving water and (b) when solutes leave the main channel and enter the porous media that makes up the bed and banks of the channel. Arrows denote the solute movement between the main channel and the transient storage zone. [Runkel, 1998; Schmid, 2004] In order to understand the physical processes occurring in the TS systems, tracer studies are typically done in the streams. The tracer studies involves releasing a slug of tracer across the cross-section of the stream or river, and measuring concentrations at various downstream locations. Breakthrough curves (time vs concentration plot) at various down- stream locations helps in the understanding of transport of solute in the streams and their interaction with the storage zones. Typically, conceptual models are used to describe the transient storage processes using mathematical equations. Various parameters including dispersion coefficient, velocity and storage zone parameters are estimated by minimizing the error between the model estimates and the observed values. The estimated parame- ters from these models are useful in understanding the physical, chemical and biological function of streams. One of the most common mathematical models used in literature to describe the longitudinal solute transport in surface and sub-surface systems is the classical advection-diffusion equation (ADE) defined as [Chiba et al., 2006] ac 30 320 —_ = 1- at ”an; D8122 ( 1) where 0: concentration of the solute; u = average velocity of the solute; D = dispersion coefficient; a: 2 space coordinate in flow direction and t = time. The above equation as- sumes that the solute plume spreads at a rate consistent with the Fickian diffusion theory i.e., a particle’s motion has little or no long-range spatial correlation. But this is not al- ways true in transport of solutes in rivers and streams as has been observed in many field dispersion studies. The observed concentration breakthrough curves (BTCs) (concentra- tion vs time) at a given spatial location have been found to have a flat long tails stretching upstream, demonstrating a greater variance than predicted by the Fickian solution. It has also been observed that that BTCs exhibit heavy leading edge, which is indicative of faster-than—fickian processes, which usually occurs because of preferential flow paths in a stream. These behaviors are referred as anomalous diffusion in the literature because it does not follow the classical bell shaped breakthrough curve usually associated with Fick- ian diffusion. One of the possible reasons for the heavy falling edges in the BTC is that the solute particle may be sorbed to solids or get diffused into regions where the advective flux in negligible (e.g. dead zones) resulting in heavy-tailed falling limb BTCs. Therefore, in order to seek an explanation for the physics underlying the observed non-Fickian disper- sion a wide range of models have been proposed in the past. Some studies have focussed on the use of TS models based on the solute exchange processes between the main channel (e.g. stream, river) and the storage zone (eg. dead zone, hyporheic zone) [Bencala and Walters, 1983; Davis et al., 2000; Fischer, 1979; Lees et al., 1999], while some models have used scale-dependent dispersion (i.e., increase in dispersion coefficient indefinitely with downstream distance) models to describe the faster-than-Fickian processes, also called super-diffusion [Berkowitz and Scher, 1995; Worman, 1998]. One of the main reasons for using scale dependent parameters were in order to describe the long-range spatial and/ or temporal correlation of particle movement observed in various field studies. However, the vast majority of these models, assume, either explicitly or implicitly, an underlying Fick- ian transport at some scales [Sposito and Jury, 1986], which assumes that the particle’s motion has little or no spatial correlation. One of the approaches used in the application of the TS models based on Fickian dispersion involves the use of a set of constant param- eters for each breakthrough curve within a selected stream reach. The difliculty with this approach is that since each breakthrough curve is treated separately by fitting the model to the observed data, the estimated dispersion coefficients contain errors associated with the non-Fickian behavior of the solute plume. This situation is further complicated due to the fact that other parameters in the TS model produce an effect that is essentially similar to the dispersion coefficient (e.g., changing the peak in the breakthrough curve). The result is a set of parameters that do not describe the observed data. Examples of sit- uations in which the TS model produced physically unrealistic parameters are described in Fernald et al. [2001] and Phanikumar et al. [2007]. One of the objectives of this work is to propose an alternative method of describing dispersion using the concept of fractional- in-space diffusion in-order to describe the non-Gaussian rising-limb of the BTCs observed in field studies as well as better constraining the stream solute transport models. When the second order diffusion term is replaced with a fractional order term, the classical ADE becomes the fractional ADE (fADE). We use the term “fADE”with the understanding that only the dispersion term is described using a fractional-in—space derivative. In a similar way, when the Fickian diffusion term is replaced with a fractional diffusion term in the well-known TS model [Bencala and Walters, 1983], we get the fractional-in-space TS model which we refer to as the FSTS model in this work. More discussion about the fADE, FSTS and their mathematical properties as well as their numerical solution will follow in Chapters 3, 4 and 5. FYactional diffusion relies on the concept of fractional derivative Operators which we discuss in the next section briefly. 1.1 Fractional Calculus and the Fractional Derivative Fractional calculus has a long history, having been mentioned in a letter from Leibniz to L’Hospital in 1695 [Weilbeer, 2006]. In his letter Leibniz writes: “You can see here sir, that one can express a term like (1235 or dl‘zx—y by an infinite series, even though it seems to be far from the geometry, which usually only considers the differences of positive integer exponents or the negatives with respect to sums, but not yet those, whose exponents are fractional. It is true that it is still to show that it is this series for dl‘zx ; but not only this can be explained in a way. Because the ordinates a: are expresses in a geometric series, such that by choosing a constant dfi it follows that da: =xdfl: a ,or (if one chose a as unit) dz = xdfl , meaning ddx would be 23.332, and c132: would be = 2:333 etc. and den: = 233-58 . And thus the differential exponent has been changed by the exponents and by replacing d6 with da: : z, yielding den: = mean” Thus it follows that (12:1: will be equal to :1: m. It seems like one day very useful consequences will be drawn from this paradox, since there are little paradoxes without usefulness.” Other leading mathematicians like Fourier, Bernoulli, Euler and Laplace were among the many who had dabbled with fractional calculus and the mathematical consequences of it [Nishimoto, 1991]. The fractional integral and derivative operators are an extension of definitions from integral values to real values, similar to fractional exponents from integral exponents. There exist numerous definitions of fractional derivatives including the definitions of Riemann-Liouville, Griinwald-Letkinov, Weyl, Caputo, Marchaud, Riesz, and Miller Ross [Kumar and Agrawal, 2006]. The most popular definitions in the world of fractional calculus are the Riemann-Liouville, Griinwald-Letnikov (GL) and the Caputo definitions. A number of textbooks have also been published over the past few decades in the field of fractional calculus which treats the subject in-depth [Miller and Ross, 1993; Oldham and Spanier, 1974; Podlubny, 1999]. One of the main differences between the integer-order derivatives and fractional-order derivatives is that while integer-order derivatives are local, fractional-order derivatives have a non-local behavior. Due to this fundamental difference, the fractional derivative at a point depends upon the characteristics of the entire function and not just the values in the vicinity of a point; hence a differintegral operator (function consisting of both differential as well as integral operators) is used to define them [Ochoa-Tapia et al., 2007]. As the function is defined using an differintegral operator a lower limit and an upper limit is required to define the fractional derivative. The plot shows the comparison between fractional derivative and integer order derivative of a function f (x) = 2:2 and f (:c) = 33. We present the formal definitions, various properties and numerical approximations to the fractional derivatives in Chapter 3 4 I I I , ’I ’I 305- ’ - I ’I a: I 3L Derivative N n a: I I I I I . : 7. N I I I I I I I I I I I \I |\ I ‘ \ I s l \ I \ | \ I I I I I I I I l L a (II I 9 II .3 UI \ I 0.5 - , ’ - I I I 0 ’ l a 1 0 0.5 1 1.5 2 x Figure 1.2. Comparison of integer and fractional derivative of f (1:) = 3:2 1.2 Numerical Methods Closed-form solutions of the governing equations exist only for limited cases (i.e., bound- ary and initial conditions). For the TS model, such solutions are even rare and the few solutions that exist need to be evaluated numerically by evaluating infinite series and com- plex integrals. The computational time and effort required to evaluate complex analytical solutions is comparable to the time it takes to numerically solve the governing equations, therefore numerical solution of the governing partial differential equations (PDEs) is the preferred approach in recent years, especially with the advent of high-speed processors and compilers that exploit the architecture of these processors. The main advantages of using numerical methods are that they are not limited to constant parameters or simple initial and boundary conditions common to analytical solutions. Additionally, they are com- putationally effective and easier to program and more effective when running parameter estimation models. One can use any of the three approaches to solve the ADE numerically - Eulerian, La- grangian, or mixed Eulerian-Lagrangian [Neuman, 1984] approaches. Similarly the fADE can be solved using the above three approaches but most of the work reported in the liter- ature used an Eulerian approach [Schumer et al., 2001; Meerschaert and Tadjeran, 2004; Deng et al., 2004] or a mixed Eulerian-Lagrangian method. In the Eulerian approach, the advection and fractional diffusion terms are solved simultaneously by numerically approx- imating them individually. Since, we find many situations in surface as well as sub-surface flow are advection dominated; a finer finite difference mesh is required to minimize the errors and capture the advection front accurately. Since, hydrologic simulations are run for large time intervals, any reduction in grid size will translate into a significantly lesser computational time. By using a higher order scheme for solving the advection term better accuracy can be achieved using a larger grid size. In the literature, there are many higher order numerical schemes (e.g. compact, WENO, spectral methods) to solve the hyper- bolic system of partial differential equation (e.g. advection ) as well as parabolic system of partial differential equation (e.g. diffusion) but these schemes have not been used for solving the fADE. Operator splitting methods in which the advection and dispersion terms are solved using separate numerical schemes (best suited for their class) are an attractive alternative but such methods have not been solved in the context of the fADE or the FSTS models. One of the objectives of this work is to fill this gap. We solve the advec- tion and the fractional diffusion equations independently using Eulerian techniques and combine them using Operator-splitting methods. We restrict our study to the usage of two higher order accurate advection schemes: a weighted essentially non-oscillatory (WENO) scheme [Liu et al., 1994] and a fourth-order accurate compact scheme [Demuren et al., 2001] for solving the advection terms. The fractional diffusion term is solved using the Griinwald-Letnikov [Oldham and Spanier, 1974] and Caputo [Caputo and Mainardi, 1971] definitions. Although fractional calculus has a history that is as old as the calculus itself, numerical solution of fractional differential equations is a relatively new area of research. Some of the traditional approaches to handling boundary conditions and approximating derivatives are not directly applicable while dealing with fractional PDEs. For example, depending on the definition of the fractional derivative being used, the (fractional) deriva- tive of a constant may not be equal to zero. This result has immediate consequences from the point of specifying boundary conditions to transport problems based on conservation laws. This means that certain definitions of the fractional derivative, although mathemat- ically well defined, can not be used in describing solute transport in streams. Therefore a systematic evaluation of different definitions of the fractional derivative is important to understand the relative merits of various approximations. The main objectives of this study are as follows: 1. To use fractional-in-space dispersion in the TS model (i.e., the FSTS model) to describe the non-Gaussian rising-limb of the BTCs, in order to better constrain the stream solute transport model; 2. To evaluate operator-splitting approaches applied to the fADE and compare the performance of two higher order accurate advection schemes (WENO and compact); 3. To compare the numerical approximation of a fractional order dispersion term using the Caputo and the Griinwald -Letnikov definitions to model fractional diffusion; 4. To demonstrate that the new FSTS model can be used to describe field data by estimating the TS model parameters in two Michigan streams. The next section gives the brief outline of the thesis. 10 1.3 Outline of thesis In this thesis, we systematically analyze all the model equations (or sub-components) of the FSTS model: the hyperbolic advection equation, the parabolic fractional dispersion equation, the fADE and the fractional transient storage (FSTS) model . This approach allows us to examine the numerical approximations and the errors involved by comparing numerical solutions with analytical solutions for the respective sub-models. In Chapter 2, a brief overview of the finite difference techniques to solve one- dimensional hyperbolic equations is discussed. A comparison of the WENO and fourth- order compact scheme is reported using a test case whose solution resembles the move- ment of a slug of tracer in a river. Our main objective here is to discuss the two advection schemes and identify the better method for solving the advection part of the fADE equa- tion. In Chapter 3, we discuss the two well known definitions of approximating the fractional derivative numerically: GL and Caputo. A higher order numerical approm'mation to the Caputo derivative will be discussed. Additionally, different forms of fractional diffusion equations will be discussed. A comparison of the CL and Caputo for the space-fractional diffusion equation will made with the analytical solution using method-of-manufactured solutions (MMS) approach. In Chapter 4, we propose the new fADE model based on the Operator-splitting (OS) approach to solve the fADE. The higher order numerical scheme for hyperbolic PDE (WENO and compact) discussed in Chapter 2 will be used for solving the advection term while the fractional derivative operators (Caputo and GL ) discussed in Chapter 3 will be used to solve fractional dispersion using OS method. A comparison of the the new numerical model will be made for describing the movement of a slug of tracer in a river. In Chapter 5, we use the new numerical fADE model proposed in Chapter 4 for solving the solute transport using the fractional transient storage model (FSTS) and compare it with the analytical solution by De Smedt et al. [2005] for second-order diffusion case. We test our FSTS model by estimating various field parameters for the Red Cedar River and the Grand River. 11 To as H AR: .8“ m>3e>flop Enosoflm .mé 9.5mm“.— To ma x¢mow. x F o mN Wm m6 mN; xmomm. x 1.5.. .30 p o m; m.~ m.m BAIIEAIJBG |EUOI139.I:| 12 CHAPTER 2 Advection Schemes Advection and diffusion processes fall into hyperbolic and parabolic systems of partial dif- ferential equations respectively. The combination of these classes of equations is among the most widespread in all of science, engineering, finance and other fields where mathe- matical modeling is involved. Very often one of the processes is dominant, as is the case in surface and sub-surface contaminant transport where advection is dominant. Advection- dominated transport processes are characterized by sharp fronts which are often difficult to capture without using a higher order numerical method or fine grid resolution. Since numerical methods that work best for diffusion (parabolic class of PDEs) are not always well-suited for describing advection (hyperbolic class of PDEs), we take advantage of Operator-splitting approaches to solve the fractional advection dispersion (fADE) equa- tion by using different numerical schemes to describe the two processes. In what follows, we limit ourselves to one spatial dimension due to our interest in stream solute transport in the longitudinal direction. Extensions to multiple dimensions are generally straight- forward and can be based on the concept of dimensional splitting [Khan and Liu, 1998; van der Houwen and Sommeijer, 1997]. The unsteady one-dimensional advection equation written in conservative form is given by: (9c 80 —— — = 0 0 < t > 0 2.1 6 t + ”82: , :1: < oo, ( ) where c denotes the concentration, 21 the mean velocity, :1: the distance and t denotes time. The analytical solution to this problem is simply the advection of the profile given 13 by the initial condition without any attenuation in the peak. For a finite-difference nu- merical scheme to be useful for solving the above equation, it must satisfy certain re- quirements such as accuracy, stability, consistency, convergence and efficiency [Roache, 1998]. Lower-order schemes introduce significant smearing while describing sharp fronts due to excessive numerical diffusion. Numerical errors generally build up with time and contaminate the solution, particularly in situations where transport simulations must be run for large time durations. Since, numerical (or false) diffusion results from an error in the finite-difference approximation of its continuum counterpart (i.e., discarding the higher-order terms in a Taylor series expansion), grid refinement can be used to obtain accurate solutions using lower-order schemes. However, this is not an attractive Option for describing field data since parameter estimation (which requires thousands of model runs if global search algorithms are used) is often an important component and lower- order schemes are unattractive as the computational effort involved can be prohibitive. Spectral, compact and essentially non-oscillatory (ENO) schemes are attractive classes of high resolution numerical methods that are known to have higher resolving power (points per wavelength or PPW) compared to lower—order schemes. Upwind and high-resolution schemes belong to an active area of research and a comprehensive review of all these methods is beyond the scope of this work. Hussaini et al. [1997] provides the background and a review of some of the earlier deve10pments in this field. In this chapter we consider and evaluate two types of numerical methods in their ability to describe the hyperbolic advection process - a fifth-order accurate Weighted Essentially Non-Oscillatory (WENO) scheme [Shu, 1997] and a fourth-order accurate compact scheme [Demuren et al., 2001]. Although spectral methods are attractive in situations where the imposition of periodic boundary conditions does not pose a problem (e.g., in atmospheric science), they are gen- erally less flexible in specifying boundary conditions for stream solute transport modeling; therefore these methods are not considered in this work. An important class of numerical methods for solving the advection equation includes the Lagrangian and semi-Lagrangian schemes [Wallis, 2007]. These methods are particularly attractive as they do not exhibit false numerical diffusion and their time-step is not limited by the Courant condition. The 14 Courant condition can be restrictive when spatial variations in velocity are encountered (e.g., groundwater). Solute transport equations for streams are often applied on a reach basis using a constant average velocity for each reach. Therefore, the Courant condition may not be very restrictive depending on the scheme used. Lagrangian advection schemes are therefore not considered in this work. 2.1 The Weighted Essentially Non-Oscillatory (WENO) Scheme The weighted essentially non-oscillatory (WENO) finite difference schemes have become one of the most popular methods in solving equations based on hyperbolic conservation laws. These methods are often used in computational fluid dynamics (CFD) to simulate incompressible and subsonic compressible flows. The primary motivation for using these schemes has to do with their ability to reduce or eliminate spurious oscillations near discontinuities (Gibbs phenomenon). Traditional finite difference schemes (e.g., upwind and central schemes) are based on fixed stencil interpolation which works for globally smooth problems. For example, in a simple finite difference scheme information for cell i is based on information at neighboring cells (e.g., ....2—2, i— 1, 2', 2+1...) It is well established in literature that the oscillations encountered in using the fixed stencil approximations suffer from numerical instabilities in nonlinear problems containing discontinuities. The oscillations do not decay even when the finite difference mesh is refined. In order to remove these oscillations the essentially non-oscillatory scheme (ENO) was first pr0posed by Harten et al. [1987]. In the END scheme instead of using a fixed stencil for interpolating the concentration at cell i, a functional criterion based on the size of the Newton divided differences between the candidate stencils on the left ( e.g., i—1,z'—2, ...) and right (e.g., i+1, 7+2, ....) is used to determine the local stencil for interpolation. This technique is more robust compared to other functional criteria such as highest degree divided differences among all candidate stencils and picking the one with the least absolute values. Even though this approach is robust for a wide range of grid sizes, Act, the scheme is usually of lower order accuracy, because all candidate stencils are not used in the final interpolations. In order to address 15 the problem of achieving higher order acuracy without any oscillations Liu et al. [1994] pr0posed a new scheme which uses all the candidate stencils by assigning specific weights to each of them. This class of schemes are called weighted essentially non-oscillating schemes (WENO) in which a convex combination of all candidate stencils is used instead of just one as in the original ENO. The WEN O schemes thus achieve uniformly higher- order accuracy throughout the domain which is preserved for piecewise-smooth functions as well. The WEN O schemes have been extensively used, in particular to simulate shock turbu- lence interactions [Pirozzoli, 2002], relativistic hydrodynamics [Dolezal and Wong, 1995], gas dynamic problems [Serna and Marquina, 2004] as well as in image processing [Burgel et al., 2002; Fedkiw et al., 2003]. Most of the problems solved using WENO are of the type in which solutions contain both shocks and rich smooth-region structures [Cadiou and Tenaud, 2004; Sebastian and Shu, 2003; Tai et al., 2002]. In environmental engineer- ing and the earth sciences, shocks can occur while describing the transport of conservative and reactive solutes in streams and ground water, sediment transport [Tsai et al., 2004], dam break situations, chemotaxis and bioremediation [Gallo and Manzini, 1998]. 2. 1 . 1 One-Dimensional Reconstruction Consider the model equation for unsteady, one-dimensional advection in conservative form: 80 8(uc)_ 87+ 8:1: —0 (2'2) In the above equation uc represents the advective flux. Replacing the flux (uc ) with the function f (c) we get 80 + (9 f (c) at 82: = 0 (2.3) which is the one-dimensional conservation law [Costa and Don, 2007]. We discretize the space into uniform intervals of size A2: and denote z,- 2 film. Variables evaluated at the spatial location 2:,- will be identified by the subscript 2'. The conservative finite difference 16 form of equation 2.2 can be written as: g: _Aix (13% — 12%) (2'4) where the numerical flux fi+§ is the average value of the flux over cell 2'. This is de- termined using WENO reconstruction. For a (2k - 1)th order WENO scheme, we first compute k numerical fluxes given by:1 =Zlf,_ _,.+,- r=0,. .,-—k 1 (2.5) j==0 f2+2 Since we are considering a fifth order scheme, 1: = 3 and the concentration in cell 2' is given by: 124% = f(C2—2,---2C2'+2) (2-6) where f () is determined using the WEN O algorithm, which takes into account any shocks observed in neighboring cells using nonlinear weights. For the fifth order WENO scheme the numerical flux is estimated by using a convex combination of all the three numerical fluxes (k = 3) as shown below * __ + 2+ 7+ 11+ + + ++ f,,,—wo (aft—2'sf'—1+sf')+w1(‘ifu—‘fiiif +5.1) ++ (2.7) +1”2 (lift +6fi+1“ 6fi+2) where the nonlinear weights w; (k = 0,1,2) are defined in the following way + + + 20+ 2 00 20+ = 01 211+ = a2 (2.8) 0 abl+aif+cxéF 1 0:6" +¢3£1I+Oz§F 2 ag+af+a§ where 2 2 2 + _ I 1 + _ 6 1 + __ 3 I 2 9) a —- a — a — - 0 m (6+ISOZF) 1 To (EMS?) 2 It) (6+152ZF) ( Here the I 5",, (k = 0, 1, 2) are the so called “smoothness indicators” which measures the smoothness of the function f (2:) in the stencil and e is a small constant used to avoid the denominator to become zero and is typically taken as 10-6. The smoothness indicators (15),) are given by: 2 [Sah==lfg(fi_2—2fi_.1+fi) +%(f2-2 4f11+3fi2) [5+ = g (fi_1 — 2fz. + fi+1)2 + % (fH -fz.+1)2 (2.10) [8+ = i301“ 2fi+1+ fi+2)2 +I(4f2+1+ f2’+2)2 17 Similarly, f;_% can be found by replacing 2 with 2' — 1. Substituting 124% and f;_% in equation 2.3, we get an approximation to the derivative (ac/8t). For stability, it is important that upwinding is used in constructing the flux. Therefore, in order to determine the direction of flow, the easiest and most inexpensive way is by determining the Roe speed [Shu, 1997] which is given by: . 1 = Eli—ff (2.11) 2+2 C2+1 - C; If a.+1 Z 0 , the flow direction is from left to the right and we would use f __ 1 while 1 2 2+2 for a. 1 < 0, the flow direction is from right to left and we would be using f .— . In 2+2 2+1/2 the above discussion we estimated f .+ In a similar manner one can estimate f2_+1/2’ 2+1 2' which is similar to the above method with different weight coefficients [Shu, 1997]. So far, we have discussed only spatial discretization. In order to achieve higher order ac- curacy in time, a third order Total Variational Diminishing (TVD) Runge-Kutta method [Shu, 1997] was used. TVD-based time discretization is used since non-TVD based Runge- Kutta time discretizations can generate oscillations even for WENO—based spatial dis- cretization [Gottlieb et al., 2006]. The following TVD Runge—Kutta algorithm was found to produce excellent results [Shu, 1997]. c1 = c” + uAt (__af(C;(n)) 02 = gen +21[c1 + uAt (_6f 2(1)) (2.12) c”+1 = §cn + §c2 + guAt (__g)_8f(cx(2)) The above algorithm is stable for a CFL < 1/4 [Gottlieb and Shu, 1997]. The scheme was found to be uniformly fifth order accurate including at smooth extrema [Jiang and Shu, 1996]. We considered a test case to examine the order of the above scheme with periodic boundary conditions by considering the following example in which a sine wave is advected in the positive x—direction with a constant speed. 80 BC a+ua—x_0, —1<:1:0 = max ( Comm — Cnuml) (2.16) where N is the number of grid points and Candy and Cnum are the analytical and numerical solutions, respectively. The fifth order accuracy was achieved at about 40 grid points (Table 2.1) which is in agreement with the results of [Shu, 1997]. We therefore conclude that in order to achieve fifth-order accuracy for a WENO scheme, a minimum grid resolution is required. Based on our test case this requirement can be expressed as: Ax S 0.0250 for u = 1. 19 2.2 A Fourth-Order Compact Scheme with Spectral-Like Reso- lution To produce nth-order accuracy, most numerical schemes require a stencil of (22 + 1) grid points. compact schemes, on the other hand, require fewer than (n + 1) points to achieve the same level of accuracy (hence the name compact). A well-known example involves the central difference approximation of the second derivative of a function. To achieve second- order accuracy (in space or time), the scheme uses information from three grid cells 2', 2+1 and 2'- 1. Compact schemes are attractive as they can achieve higher-order accuracy using just two points. At least two other reasons make the compact schemes attractive - the two-point stencil is ideally suited for making computations on highly irregular or stretched grids (e.g., to resolve the structure of a turbulent boundary layer) and the schemes have exponential error convergence similar to the classical spectral methods. By contrast, the error convergence rate for the WENO schemes discussed above is linear, which limits their ability to achieve higher-order accuracy in large computational domains. Compact schemes achieve higher-order accuracy by treating the derivative terms in the governing equation as unknown functions to be solved (in addition to the function itself). This means that for the advection equation, the concentration c and its first deriva- tive c’ will be solved at the end of every time step. A fourth-order compact scheme was introduced by Gupta, Manohar and Stephenson in the early 80’s for solving the advection- diffusion equation (ADE). It has also been used to solve the Navier-Strokes equation [Li; Phanikumar, 1994] and Euler equations [Abarbanel and Kumar, 1988]. Lele [1992] found that implicit compact schemes have better fine-scale resolution, and yield better global accuracy than standard higher-order finite difference schemes. Here we briefly describe the algorithm used to solve equation 2.1 using a fourth-order implicit compact scheme. Additional details and examples applications to other areas (e.g., aero-acoustics) are avail- able in [Demuren et al., 2001]. This scheme was also used to model solute transport in soil columns [Phanikumar and Hyndman, 2003]. The first derivative term appearing in equation 2.1 was obtained by solving the following tri-diagonal matrix system [Lele, 1992]: I I I “I OICi_1+Cz+O!ICi+1 = é—A—E(Cz+1—Ci_l) (2.17) 20 where c; = 5%, a] = 11;, and a1 = 154—. The LHS of equation 2.17 contains the unknown derivatives at grid points 2' and 2' i 1, while the RH S contains the known function values (e.g. concentration) at at the grid points 2 :l: 1 . In matrix form the system of equation 2.17 can be written as AxC’ = Bx (2.18) where Ax is a tridiagonal N3 x Na; matrix (Nx is the number of grid points), and Bx is a vector. The resulting tridiagonal system of equations was solved using Thomas algorithm [Press, 2002]. For nonperiodic boundaries, one sided finite difference approximations are required to close the system of equations at the boundary points: 2’ = 1 and 2' = N. A third-order approximation was used at 2' and 2' = N given by [Demuren et al., 2001] C] + 2c’2 = A—lx. (IE—5C1 + 202 + $03) (2.19) A similar approximation was used at 2' = N. Storage of variables could potentially become an issue while using the compact schemes. Therefore, we used a three-stage—third order Runge-Kutta scheme proposed by [Lowery and Reynolds, 1986] with low-storage requirements for temporal differencing. By low-storage we mean only two storage locations (one for time derivative and one for the variable itself) are required for time advancement, which is achieved by continuously overwriting the storage location for the time derivatives and unknown variables at each substage for M = 1, 2 and 3 as * (M) _ M : M—l Hz. _a H2 (2.20) cE-MILI) = c?! + bM+1uAtHz-(M) (2.21) M 3,94) = 71,9”) _ ufl< ) 2.22 8,, < > where coefficients aM and 0M are given in Table 2.1 and M -1) ..(M—1) _ 1c} H,- — 08% (2.23) 21 Table 2.1. Coefficients of a three-stage—third-order Runge-Kutta scheme, from Lowery and Reynolds (1986) M b Table 2.2. Order of accuracy for the WENO-Roe and the compact schemes or or 2.3 A Test Case: Advection of a Gaussian Pulse To evaluate the performance of the two advection schemes and to assess the nature and magnitude of errors introduced by the schemes, we used a test case for which the solution is similar to the instantaneous release of a slug of tracer in an infinite domain. The slug will maintain the same initial profile (in our case Gaussian) for all times since dispersion is zero. The initial condition (which is also equal to the exact solution for all times) is given by: C(x, 0) = 05 exp [— (32 111(2)] , -20 g a: g 450; u = 1 (2.24) This solution will test the time-advancement of the initial profile using these schemes with the equation 2.1. Table 2.2 below summarizes the comparison of the results obtained using the WENO and the compact schemes. All the comparison were done at time t = 5. We used two criteria to determine which scheme will be preferred for solving the fADE equation. Since our interest was in running the models for large times (while describing field data), schemes that give stable results for the largest times step size At with the least numerical dispersion will be preferred. Additionally, the domain sizes for simulations used in surface water modeling are typically large, therefore, schemes that give reasonably accurate solutions with a larger step size will be preferred because of the 22 lower computational time. The stability of the numerical schemes is determined using CFL (Courant-Friedrichs-Lewy) number defined by: uAt CFL -—- E (2.25) For a scheme to be stable the general criterion is that CFL S 1. For an explicit scheme, if CFL 2 1, the scheme will become unstable and will not converge. 0-5 I r I I I I O Analytlcal 0-5 ' —-*-- Weno ' - -0- - Compact 0.4 - —>— First-Order Upwlnd - 0.3 - _ o 0.2 0.1 -0.1 . -15 -10 -5 0 5 10 15 20 x Figure 2.2. Comparison of the compact, WEN O and first-order upwind scheme for A2: = 1, u=1,CFL=0.1att=5 We initially used a crude grid size and found that for CFL = 0.1, both WENO and compact schemes perform reasonably well in describing the peak with the WENO scheme giving a lower Root Mean Squared Error (RMSE) (Figure 2.3). The first order upwind scheme introduced significant numerical dispersion and was not able to describe the peak accurately (Figure 2.3). Comparisons for a larger time using the same crude grid showed that while the compact scheme tends to produce oscillations, the WENO scheme intro- duces numerical dispersion. It is therefore, clear that a more refined grid size should be used for accurately describing the peak. In a second simulation using fine grids, we found 23 Variation with A x 10 W I l I F 1 I 4 —l—WENO 10 - g \ - ... compac‘ 5 ’~.. + +FIrsI—ordor -. ~ .- - - - 10 " ---- -o- ______ 104 I I I I I I I I 01 0.2 03 04 0.5 06 07 08 09 1 Ax Variation with CFL number 1 I I I I I I t u- u- ------------ - - .... ---------- ‘- '— I I I I I I I 03 04 05 0.6 07 08 09 1 CFL Figure 2.3. Plot showing that comparison of RMSE error (log—scale) for compact, WENO and first-order scheme for different CFL and A2: values that both compact and WENO schemes were able to describe the peak without any os- cillations. As the grid size (Ar) decreased the compact scheme was found to give a lower RMSE error (see Figure 2.3) while as CFL number decreased both of them approached similar orders of RMSE error. 2.4 FFT Analysis of Advection Schemes Numerical solutions of differential equations generally contain both dissipative (i.e., errors in the peak or amplitude) and dispersive (i.e., phase errors) errors even if the model equation is non-dispersive as is the case for the pure advection test case considered here. Numerical solutions can be compared and analyzed based on a number of metrics such as the truncation errors and rates of convergence of the schemes, as well as dispersive and dissipative behavior. It is often a challenge to select criteria that fairly evaluate different 24 0 _ I I I I _ .5 . 0 Analytical .’ ‘ a. +Weno 0-45 ‘ '0 ‘, - e - Compact ' I \ 0.4 - 'I ‘. H b 0.35 - e “ - I I 0.3 P ‘ - o 0.25 - ._ 0.2 - I .. l 0.15 F g -‘ 1’ ‘ 0.1 - 0 - I \ a a. 0.05 - ' a _ c .’ 3 o 0-e' I I J I J 240 245 250 255 200 distance (x) Figure 2.4. Plot showing that comparison of the compact and the WENO scheme for A2: = 1, u = 1, and CFL = 0.1 with the analytical solution at t = 250. WENO scheme disperses while compact scheme oscillates classes of numerical methods. A general methodology which can be used to compare numerical solutions obtained from different schemes is a mathematical framework based on Fourier analysis since the method provides a great deal of information about the errors involved. Fourier analysis is most commonly used to find the frequency components buried in a noisy signal in the time domain (e.g., time series of discharge data at a gaging station). To analyze the resolution properties of the schemes, the concentration profiles as a function of distance at a given time can be transformed from the physical space to the wave number space and analyzed using FFT analysis. Results of the FFT analysis are shown in Figure 2.5. We noticed that the WENO scheme introduced significant numerical dissipation in addition to producing considerable phase errors at high wave numbers. 25 2.5 Conclusion Based on the above results, we conclude that the compact scheme provides a superior description of advection and is therefore the preferred choice for solving the fADE since WEN O scheme for large time simulations seems to disperse higher compared to the com- pact scheme. It is important to capture the peak accurately since it allows us to under- stand various effects of fractional dispersion. If numerical dispersion is introduced with an advection scheme, it will be diflicult to understand whether the dispersion observed is due to numerical dispersion or due to fractional dispersion. Our results indicated that the compact scheme produced superior solutions without any oscillations as long as low CFL numbers were used to maintain stability. 26 one 3% 33.8 a How 2538 S8988 3 @8388 $32 3 080:8 02mg 2.: .8 meaning 2: 3330 .omm H A as .8538 Rompbag on» 55 Hd H QED was A H a .H M ed .8.“ 2558 02m? 28 3388 Ho aoflneaaoo 3% @539? BE .n.n 0.53..” DOQEDCO>N>> thEDCO>m>> P wd 0.0 To Nd o P md ed To N6 0 . W 4 . miop . A q - .vl r . r . r ml m, .An 2 . .1 ._ NI . . .u v __ T w dl d u. _ n o w. w. m _ _ _ - -_ —- l —. -_ __ - __ I u —- —- 1 N __ _. H " «umanu “ I _ _ I m 0531 I I u . - . . 2 OF . . . . w 27 CHAPTER 3 Fractional Difl'usion: Preliminaries and Numerical Approximations 3. 1 Introduction Transport of solute particles in surface and subsurface water can be viewed as a purely probabilistic problem. The classical Fickian diffusion model is based on the assumption of Brownian motion. Brownian motion assumes that the particle’s motion has no long- range spatial correlation i.e., long walks in the same direction are rare. But what happens to particles with long-range spatial correlation? One can answer the above question using fractional calculus and a class of probability density functions (PDFs) described using continuous time random walks (CTRW). Using CTRW, the diffusion process can be viewed as a result of random walk in space and time. In other words, a particle can undergo jumps of random size at random times. The random walk in time and space can be represented by the joint probability density function P(:r, t), which describes each particle “transition” over a distance, x — 23’, in time t—t’, using the master equation given by [Montroll and Weiss, 1965]. Pa, t) = 6(x) / w(t')dt' t Contribution from particles that have not moved during (0, t) (3.1) t 00 + ft/J(t — t’) [ f Ma: — :r')P(:r’ — t’)d$’] dt’ 0 -oo Contribution from particles located at x and jumping to x’ during (0, t) This equation may be used to describe various kinds of diffusion equations, depending 28 In 4W Tn = waiting time 01(1) = waiting time PDF Cm =Iump Aft) = jump size PDF Figure 3.1. The random jump in space and time in a CTRW model upon how the waiting time (7n) and jump distributions (A(()) are chosen. Fickian diffu- sion is the limiting case of CTRW. The assumptions being that the waiting times between the jumps (Tn) are exponentially distributed and the jump distances ( Cn) are normally (Gaussian) distributed. But in case of anomalous diffusion, which is a more general case, the jump distributions are non-Gaussian with divergent moments to account for anoma- lously large particle movements also known as “Lévy flights”. Similarly, the waiting times can have an arbitrary distribution instead of an exponential distribution for a second-order diffusion. An increasing number of processes are being described by anomalous diffusion where the random waiting times and jumps do not follow the Gaussian distribution. In these cases, the overall transport is better described by steps that are not independent (long range) and can have arbitrarily large” waiting times. 3.2 Fractional-in-space Diffusion Fractional-in-space diffusion is used to describe faster-than—Fickian growth rates, skewness and heavy leading limbs of BTCs in streams and heterogenous aquifers. The diffusive flux for fractional-in-space diffusion can either be based on a fractional Fickian flux [Zhang et al., 2007] or on a regular integer order flux (the well-known Ficks law) [Schumer et al., 2001] In the first case, the diffusive flux is fractional order and non-local, and is given by [Yong et al., 2006]. —1 ac_e[()aa C (3.2) 5? .- 51.: :17 8130—1 While for the the second case, the integer order divergence in the above formulation was 29 replaced by a fractional analogue resulting in fractional-in—space diffusion given by: ac 30-1 ac B—t ‘- a—mafi [Mtg] (33) Both equations 3.2 and 3.3 are equivalent for the case, when D(:r) is constant and inde- pendent of the spatial variable a: [Yong et al., 2006], reducing it to g)” 300 at _ 82:“ (3.4) 3.3 Fractional-in-time Derivative: Sub-Diffusion Fractional-in-time derivatives are used to describe long waiting times between the ran- dom jumps, usually accompanied with travel times of particles longer than the expected classical Fickian diffusion, called sub-diffusion. For example, Haggerty et al. [2000] used fractal, or power law to describe a distribution of rates instead of using a single rate model. This can be thought to be as scaling in time, where the coefficients on the time operators need to be adjusted with time. A fractional-in-time derivative can be used to describe rate coefficients which are scale invariant in time. One classic example of a sub- diffusive process involves solute transport in a stream in which the exchange process is between the main channel and the storage zones. The storage zones trap the contaminant particles, hence the longer time, characterized by heavy falling limb in BTCs. In a stream transient storage (TS) models, this exchange process is usually described as a first-order mass transfer process, which can be viewed as a special case of the processes described by a fractional-in—time derivative. The governing equation for the sub-diffusion is given by Gorenflo et al. [2002]. 050 320 W =3 DEE-f Where 0 < )6 $1 (3.5) Other variations of the fractional-in-time diffusion are as follows 300 a 30 Tia _ 5:; (195) (3.6) One form of the diffusion equation that is fractional in both space and time can be represented as aa'c aa-l ( 30) ad = W a; (3”) 30 where a and a’ are fractional exponents in space and time, respectively. Fractional diffusion relies on the concept of fractional derivative operators which we discuss in the next section. We limit our discussion to the definitions and theory relevant to this study. We introduce formal definitions of fractional derivatives and discuss the numerical methods to approximate them. Additionally, a comparison between the two well known definitions of fractional derivatives (Caputo and Griinwald-Letnikov) will be presented. 3.4 Fractional Derivative: Multiple Definitions Although the integer derivatives of a function are mathematically well-defined with well- known physical and geometrical interpretations, the fractional derivatives are much less explored and understood. For example, it is well-known that there are no geometrical interpretations of a fractional derivative as of today [Podlubny, 1999]. While there are a number of areas in which physical interpretations of the fractional derivative are meaning- ful and useful, one has to come to terms with the fact that there are multiple definitions of the same fractional derivative. While some definitions are well-suited for certain types of applications, other definitions are not. Therefore, as mentioned in chapter 1, we consider (and evaluate) three main approaches to defining the fractional derivative: the Riemann- Liouville definition, the Griinwald-Letnikov (GL) definition and the Caputo definition. Here we introduce their formal definitions and in the next section we introduce numerical methods to approximate these derivatives. In order to understand the fractional deriva- tive it is important to understand the concept of a fractional integral, both of which are closely related. This section summarizes some of the prOperties of the fractional order integral and differential operators relevant to our work. Additional details can be found in Podlubny [1999] and Oldham and Spanier [1974]. The traditional expression of repeated integration of a function for an integer order is given by: a: $1 $2 $n-1 aDgnflx) =/dx1/da:2/dx3... / dccnf(1:n) (3.8) 31 The order of integration can be interchanged and equation 3.8 can be written as aD;"f(:r) = (...—:17? [x (x — y)"‘1f(y)dy (3.9) a To extend the above equation to non-integers, the integer is replaced with a real number and the factorial in the denominator is replaced with the more general gamma function I‘(), which is an extension of the factorial function to complex and real number arguments. The gamma function is defined as 00 I‘(a:) = /e_tt$—1dt (3.10) 0 and satisfies the following prOperty I‘(n) = (n —1)! (3.11) Replacing the factorial with the more general gamma function, we get the Riemann- Liouville (RL) fractional integral defined as [Oldham and Spanier, 1974] 1 a: —u _ _ u—l Fractional derivative operators of order or act as an inverse of the fractional integral of order a, i.e., aDgaDEO‘Nl‘) = f (as) (3-13) dN aDg‘flx) = ——W [aDEVf(m)] where l/ = N — a (3.14) d2: d” —(N—a) d0 Guam) = M [an f(x)] = 25:3 (x) (3.15) In a similar way, the RL derivative can be defined as [Samko et al., 1993] I‘ZnI—al (£3) (1; (x—ujlguln+1 du n — 1 S a < 12. Mrs) = (3.16) gm) a = n 32 3.4.1 The Riemann-Liouville Derivative Fractional derivative of a function f at a point a: can be defined based entirely on infor- mation either on the left-hand side or the right-hand side of a spatial location in question (Figure 3.2). If the independent variable is time instead of space, then the left fractional derivative can be interpreted as an operation on the past states of the function (f) while the right derivative can be interpreted as an operation performed on the futures states of f [Podlubny, 1999]. The left RL derivative at a: is defined as a: a _ 1 an f(u) anf(x) _ l"(n __ a) 5312. (:1: _ u)a—n+1 a. du (3.17) Left Right anaf(X) xDbaf(X) r-A-v—H 4. : : a x b Figure 3.2. The left and right RL derivatives at a point a: while the right RL derivative at 2: is defined as a _ (_1)n an f(u) mob f(x) - r / (u_ du (3.18) (n — a) 32:" m)a—n+1 where n - 1 g a S n and I‘() is the gamma function. Although the RL derivative is a well-defined mathematical function, one can easily encounter problems when applying it to physical problems. In particular, the fact that the RL derivative of a constant is not zero can cause problems while specifying boundary conditions (e.g., continuous slug release of a tracer of a known constant concentration in a stream). For example, consider the fractional derivative of a constant C', C anC = r(1— a)(:r — a)“ (3.19) where 1 S a S 2. As observed, the RL derivative is singular at the lower limit. This presents a problem when using the RL expression for estimating the derivative at the initial boundary since . a _ mama) — oo (3.20) 33 f"(a) = o (3.21) where k = 0, 1, 2...n—1 and 1 g a g 2 . Similarly, the Laplace transform of the RL deriva- tive depends on the fractional derivative at zero, which introduces initial value problems (refer Podlubny [1999] for details). These drawbacks of the RL derivative Operator can be resolved by using the Caputo definition of the fractional derivative. 3.4.2 The Caputo Derivative As the RL derivative, Caputo derivative is also a dz'flerintegral expression but the opera- tion of integral and derivative is interchanged in case of Caputo derivative. The Caputo derivative can be defined as Podlubny [1999] Car 1 “’ f(")(U) u anf( )—P(n_a) a (z_u)a-n+1d (3.22) where n —— 1 g a g n and the superscript C is used to distinguish the Caputo derivative from the regular RL derivative. One of the main differences between Caputo and RL definitions is that the initial value of the Caputo derivative is in terms of integer-order derivative in Caputo unlike in RL. For example, consider the case of 1 g a S 2 commonly encountered in describing the transport of solutes, for which the Caputo derivative is defined as aCDé'flx) = 172—133 / (x fufflmdu (3.23) The Riemann-Liouville fractional derivative requires initial conditions expressed in terms of initial values of fractional derivatives of the unknown function [Podlubny, 1999; Samko et al., 1993] while the Caputo derivative, in comparison, requires that the initial condition are expressed in terms of initial values of integer order derivatives. Therefore, Caputo derivative can be used in physical problem where one can interpret the integer-order initial values whereas in case of RL derivative the fractional initial value is difficult to interpret and not known easily. This can be avoided by treating only the cases of zero initial conditions. However, for zero initial conditions the Riemann-Liouville and Caputo fractional derivatives coincide [Podlubny, 1999]. 34 Horn the above definitions of non-integer fractional derivatives, it is apparent that these derivative are non-local operators and these derivatives find important applications in various fields. A fractional derivative at a certain point in space or time contains information about the function at earlier points in space or time, respectively. Thus, fractional derivatives possess a memory effect, which they share with several materials like viscoelastic materials and polymers and this property is also important to applications in hydrology such as anomalous diffusion. 3.5 Numerical Approximations: The GL and Caputo Deriva- tives In this section we discuss the finite difference approximations to a fractional derivative using the Griinwald-Letnikov and Caputo derivative operators. The GL based fractional derivative operator is a fractional formulation using backward difference for any arbitrary function f (as), although the convergence of the infinite sum cannot be ensured for all functions. There exists a link between the Riemann-Liouville and GL approaches to fractional derivative operators and the exact equivalence between the two approaches has been shown in Podlubny [1999]. The GL approximation to the fractional derivative is given by . "" a) . aDgflx) = N11m———— (1)}: CUP—— —j-——F(+ f(a: — 3h) (3.24) where N = (.2: — a) / Am and is a positive integer] . The above equation can be simplified to the following sum of series given by: D“ f :2 w] f (x — jh) (3.25) j=0 where the weights wj are defined as a —1 wj = (1 — j ) wj_1 (3.26) The right sided GL fractional derivative can similarly be defined as 1 NPU-a) $013705): “in mf($+jh) 35 3.5.1 The Caputo Derivative The Caputo derivative can be approximated for the fractional diffusion equation using finite volume method. The scheme described below was adapted from Zhang et al. [2006] who used this method for solving the fractional advection dispersion equation (fADE). We use essentially the same approach but we only consider the fractional diffusion equation. The transport domain [0, L] is divided into N equal elements of size Ax. Using the mass- balance equation based on the diffusive flux during time period At the fractional diffusion equation can be written as Uzi-Hm — C: _ Qt+At/2 _ Qt+At/2 At — i—1/2 2+1/2 (3.27) where the subscripts 2' :l: 1 and i :1: 1/2 represent locations xiil and mad/2’ respectively as show in Figure 3.3. The two fractional diffusive fluxes are estimated using the Caputo derivative given by: xi—l/Z “2+1/2 QHI/ZZFTZ—{D—a) o/ ($¢+1/21—y)a’1%dy (329) which can be simplified and written in finite difference form as follows. Further details are provided in Appendix A. F', - ijf—j - Cf—j-r) (3-30) co l H D Qi-l/Z = F(2 — a)Aa:a ll H 1 D Qi+1/2=p(2_ a)A$a :7ij(0 2+1—j Cf—j) (3-31) The above value of diffusive flux can be substituted back in equation 3. 27 to get the final form of the discretized fractional diffusion equation as (Ct+At _ 2Ct+At+ +Ct+At) z'+1 62”“ -0f 0 + >: w (C.+1_—C:_-> ‘ __ = i=1 .7 J (3.32) At F(2 — a)Axa 1w where w)- = (j + I)“: — (3)2"a (3.33) The above equation is solved using a semi-implicit approach described by Lynch et al. [2003], where the Gaussian terms (:7 — 1, i and, i + 1) are solved implicitly while the “tail” terms (2' — 2,i — 3....1, 0) are solved explicitly. The above equation reduces to the classical Fickian diffusion for a = 2. x'-l xi xi+1 b j .g - ‘ ...— ? V xi-l/Z xi+1/2 Figure 3.3. Control volume approximation to the Caputo fractional derivative 3.6 Higher Order Approximations Most numerical solutions of differential equations involving fractional derivatives are first order accurate (either in time or space or both depending on which of the derivatives is being approximated). This is true for both the CL as well as the Caputo definitions. Higher-order accurate approximations are attractive when solutions need to be obtained with minimum computational effort, especially when parameters are being estimated. In order to implement higher order finite difference approximations we need to approximate the weights in GL or Caputo to a higher order. Podlubny [1999] showed that the weights 10):!) can be obtained as the coefficients of the power series expansion of a generating function w. The generating function w1(z) generates the coefficients for the first order derivative to first order approximation while its a-th power, given by the function w‘f generates the coefficients for the first order approximation of the a—th fractional derivative. Here 2 denotes a complex number. Once the generating function is selected, the weights can be obtained by evaluating the Fourier integral of the generating function. For example, the weights for the first order formulas presented in the previous sections can be obtained using the generating function for the a—th derivative (to first order approximating) using 37 the function wl (2) = (1 — 2)a: a °° P(a +1) k °° (a) k (1 — 2) = (—1)k 2 .-= w 2 (3.34) [:20 F(a — k +1) 1:230 1‘ Substituting z = e—i‘f’ we have n m o (1 — e-up)°‘ = Z mime—21w (3.35) k=0 (a) and the coefficients wk can be obtained from the relation: (a) 1 27r wk =-,—,—,;O/fa( 0 for a given function using the trapezoidal numerical approximation is as follows. Consider a function f (1:) between interval [0, a] divided into k subintervals [xj,:cj+1] of equal width h = f: by using the nodes .2]- : jh for j = 0,1,2,3, ....,k. The Caputo fractional derivative is given by [Odibat, 2006] [(k - 1)n—a+1 - (k - n + a — 1)kn-a] f"(0) hn—a +fn(a) I‘(n + 2 — a) k—l (k — j + 1)n—a+1 (3.42) + .2 -2(k - 22““ 1mm) 3:1 +(k __ j _ 1)n—a+1 S’Dgfe) = where n — 1 < a < n. The integer order derivatives (fn) in the above expression can be approximated using any regular finite difference numerical approximation. We used the fourth order compact scheme discussed in Chapter 2 to approximate the first derivative ( for 0 < a g 1) and second derivative (for 1 < a S 2). 3.6.2 A Test Example We tested the new numerical method by taking the a-th fractional derivative (a = 0.5) for the test function f (x) = sin(:r). Using the new numerical method discussed above, we estimated the first derivative of the function f (x) = sin(a:) using the fourth-order compact scheme at all the finite difference mesh points between 0 S :1: S 1. The analytical solution to the fractional derivative of f (2:) = sin(a:) is given by [Odibat, 2006] 00 . . CD” sina: 2 331—0: E (_1)2$22 for 0 < a < 1 (3.43) 0 3’ i_ I‘(2i+2—a) 0° i+1 2i+1 C a - 2—a (“1) ‘13 D = E . 0 Isrnx :1: {—0 I‘(22’+4—a) for 1? vEBQmHU 83.9 o “someway .8“ ozfizuoo >ofi53éwmkgu0 23 H8 Saws? 2: 9:32? 33a. .Né oEaH. 46 ooood Hoooo moooo moooo noooo woooo oHooo oHoo.o ooooo ooooo oH-H ooooo moooo woooo ooooo ooooo HHooo mHooo mHooo HHooo noooo o-“ oooo.o moooo moooo woooo HHooo vHooo mHooo oHooo vHooo ooooo w-“ oooo.o voooo uoooo HHooo mHooo wHooo omooo omooo wHooo HHooo n-H oooo.o moooo oHoo.o oHooo Hmooo mmooo hmooo nmooo mmooo mHooo o-H ooooo woooo oHooo mmooo omooo omooo omooo wmooo Nmooo omooo m-“ ooooo mHooo mmooo smooo pvooo mnooo wmooo onooo Stood wmooo v-H oooo.o mmooo ovooo ooooo Nwooo moooo ooooo ooooo whooo Yvooo m-H ooooo Hmooo noooo mmHoo voHoo owHoo meoo ooHoo mmHoo Rood m-“ oooo.o omHoo vwmoo owmoo vaoo mogo omvoo momoo Homoo onoo H-H oooo.m- mwmo. H- 23. H- mwon. H- moon- wmwn. H- mewwH- mmum. H- owmmH- mmmH . H- _ ooooH ooooH ooooH ooooH ooooH ooooH ooooH ooooH ooooH ooooH Hi o.N o.H oz" 5H ozH mH TH ”H NH H.H nomumooq A3 ado—893 o>5w>maoc EnoEowfim - 3AM??? 3:960 8:13 o ”.58on no“ w>sw>ioo 3250 23 n8 Saws? 2: mnEofi Bash. .né £2.59 47 CHAPTER 4 Fractional Advection Dispersion Equation In this chapter we present the numerical methods for solving the one-dimensional fractional-in-space advection dispersion equation (hereinafter fADE) based on operator splitting. The fADE is given by: a 30 aC_D(flaaC 30) ‘57 “5:" 526+(1‘5)_3(_3)a (4.1) where C is the concentration, t is the time, :1: is the spatial coordinate, u is the mean velocity, (1 is the order of fractional derivative, D is the dispersion coefficient and B (0 3 fl 3 1) is the skewness parameters that controls the bias of the dispersion. We will discuss the operating-splitting (OS) technique to solve equation 4.1. 4.1 Semi-Analytical Solution for the fADE Closed-form solutions to the fADE do not exist in the literature but semi-analytical solu- tions to boundary value problems can be found using the Laplace or the Fourier transform in a manner similar to what was described in Ogata and Banks [1961]. Benson et a1. [2000] gave a semi-analytical solution for the fADE equation for boundary conditions similar to an instantaneous release of a slug of tracer in a stream (the instantaneous release can be mathematically represented using the Dirac-delta function). The solution for the fADE for an infinite domain given by [Benson et al., 2000] is C(k, t) = exp -;—(1 — fi)(-—z‘k)0‘t + —;-(1 — fi)(—z’k)aDt — ikut] (4.2) where i = \/—1 and —1 _<_ 3 S 1 indicates the relative weight of the forward versus backward transition probability and 1 < a < 2 is the scaling exponent in the one- dimensional space, 11 is the velocity, D is the dispersion coefficient with units LOT—1. 48 With the notational simplification B = [cos(7ra/2)D|, and with the identities i = cur/2 and em = cos 0 + sin 0, equation 4.2 can be rewritten as 6%. t) = expi-Bt Ikl" [1 + w_ 0 (4.12) while the RHS boundary (at an infinite distance) was simulated using zero prescribed flux for all times and is given for the classical ADE by BC III—'00 Prescribed Flux boundary: The LHS and RHS boundaries were simulated using the prescribed flux boundary conditions for an infinite domain :5 6 (—oo, 00) with instanta- neous release of slug of tracer at a: = 0. For the classical ADE with an infinite domain, the prescribed flux at the two boundaries is equal to the dispersive flux at the inlet and outlet boundaries and can be written as BC —D = 0 35 ave-00 (4.14) The above is not the correct boundary condition in case of fractional dispersion, since the dispersive flux at the two outlets is a fractional dispersive flux unlike the integer-order flux in the ADE. The fractional-order flux at the two boundaries is given by: -D(fi%Z—2i¥+<1-fi>%i%) x-r-oo (4.15) _ aa-lc _ aa-lc _ The two dispersive fluxes at the two boundaries are approximated using the control volume approach described by Zhang et al. [2006] and are given by (details in the Appendix A) flw0(C-N+1 — C—N) Q N = _ D N—l _ +1/2 I‘ZZ—aiAxa + 20 (1 - fi)(C-N+1+j — C—N) J: N —1 Q _ _ D . flwj(CN—j “— CN—j—l) N—l/2 — I‘[2—a[Aza 1:0 +(1- 5)(CN—1 — CN) We equated both the above fluxes to zero, since we were using infinite boundary conditions. (4.16) 53 4.4 Solution of the fADE for Continuous Release For simulating the fADE using a prescribed-concentration boundary for the case of con- tinuous release of a tracer, we chose a domain of size L = 20m with zero concentration everywhere in the domain initially. At t = 0 a conservative tracer with a constant con- centration of CO was imposed at the LHS boundary at a: = 0 and the boundary condition remains constant throughout the period of simulation. The RHS boundary at a: = L was simulated by using zero-flux prescribed flux given by . We used fl = 1 , a case also studied by Benson et al. [2004] and Liu et al. [2004]. The fADE was solved using strang OS approach, with the advection term solved using WENO numerical approximation and fractional dispersion term was solved using GL and Caputo derivative approximations (see Figure 4.4). 1.4 1 I l l + Benson et al. (2000) 1.2 - - - - Grunwaid - t=5s t=10s —-Caputo Concentration (ppb) 9 P o on d P a P N o 5 1o 15 20 25 Distance , x (m) Figure 4.2. Concentration profiles for the continuous release of a tracer in a stream simulated using the fADE. The fADE was numerically solved using the WENO-CL and WENO-Caputo schemes and compared with the semi-analytical solution of [Benson et al., 2000]. The parameters in the solution are Co=1ppb at t: 55, 103 and 15s with u = lms‘l, D = 0.1m1'8/s, a = 1.8, CFL = 0.1 and Ax:0.1250m We find that the comparison based on the CL approach does not match with the semi- 54 analytical solution by Benson et al. [2000]. This has also been observed by Zhang et al. [2006], where the fractional dispersion was solved using the Riemann-Liouville definition which has been been shown to have exact equivalence to the CL approach used here by Podlubny [1999]. The inability of the GL—based numerical solution to correctly describe solute transport as given by the analytical solution was attributed to the fact that the fractional derivative of a constant (in our case the prescribed concentration C0) does not vanish (that is not equal to zero) in the case of the CL approach. This was not an issue while using the Caputo derivative based numerical scheme which agreed well with the analytical solution. We therefore conclude that the Caputo derivative is the preferred approach for simulating continuous release in a stream. 4.5 Instantaneous slug release We simulated the instantaneous release using a prescribed-total-flux-boundary at a: = 0 and a free drainage outlet at a: = L, commonly encountered in hydrology by using the following boundary condition, described earlier in Chapter 3 30-10 60-10 x——L 30-10 601-10 where an initial slug of conservative mass of tracer equal to m = M0 was released at t = 0 at the spatial location a: = 0. The domain of simulation is a: = [—L, L]. We simulated this boundary using an “infinite” domain with the tracer slug released at :1: = 0. We wanted to understand which combination of numerical schemes works best in describing the instantaneous slug release, which will be used later to describe field data. Therefore, we tested both WENO and compact schemes for solving the advection process while the CL and the Caputo derivative operators were used for approximating fractional dispersion. Hence, we tested all the four possible combinations of schemes, i.e., WENO- GL, WEN O-Caputo, compact-CL and compact-Caputo for simulating this case using the OS technique. 55 During our initial runs, we found that the RMSE error was not a true measure of whether or not the scheme was reasonable enough to capture the peak or the tail ac- curately. Therefore, we ran many simulations and visually checked to see within what reasonable range the schemes perform well. We wanted to know if the numerical solution can describe both the peak and the tail and if it is in reasonable agreement with the semi- analytical solution by [Benson et al., 2000]. In order to accurately capture the effects of various flow conditions occurring in the field including high advection or high dispersion, we used the cell Peclet (Pe) number as a metric to represent various conditions. For example, using the cell Peclet number we can capture the effects of fractional derivative exponent (a), velocity (u), dispersion coefficient D and the mesh size Aa: using a single parameter. The cell Peclet number for the fADE can be defined as: a—1 Pa = ."—A—:)—— (4.19) Note that for a = 2, the fADE cell Pe number is equivalent to the ADE cell Pe number. All the breakthrough curves were plotted at a distance of a: 2 10m from the point of injection of the tracer. We compared all the four schemes by plotting the BTCs and comparing then with the semi-analytical solution by Benson et al. [2000]. We used forward skewed fractional dispersion (fl = 1), since in an advection-dominated flow (e.g., streams) dispersion will be biased forward which was also noted by Benson [1998] and Zhang et al. [2006] 4.5.1 High Peclet Number: Advection-Dominated Systems We simulated a high Peclet number condition to compare the four schemes listed above with the analytical solution. A comparison is shown in Figure 4.3. We found the that the compact-based advection scheme gave oscillations in the tail region of the breakthrough curve using both the CL as well as the Caputo based fractional derivative operators. This supports the well-known fact (see Demuren et al. [2001]) that the compact scheme introduces oscillations at high Peclet number. The WENO based advection scheme was not able to simulate the peak or the tail correctly using any of the fractional dispersion schemes but did not suffer from oscillations because of its non- 56 + Benson et al. (2000) f. -a - o -WENO-GL - v WENO-Caputo ->- Compact-Caputo -1o -0- Compact-CL Concentration (ppb) I I l I I :I: I I 3.5 4 4.5 5 5.5 6 6.5 7 Time (s) Figure 4.3. Plot showing the comparisons for all the four schemes for a Pa = 6.6 (A2:=0.25m, D = 0.05m1'8/s, u = 1ms‘“1, CFL =0.1, t=153 and a = 1.8). Numerical oscillations were observed in case of compact based approaches while numerical dispersion was observed in the case of WEN 0 based advection schemes oscillatory nature. The WEN O scheme is a little dispersed as it requires a minimum grid size to achieve the fifth-order accuracy as illustrated in Chapter 2. Since a crude grid size of Ax=0.25 is used in the above simulation, dispersion error is noticed. Our runs show that for P6 < 2 all the four schemes were able to simulate the peaks well but the GL-based fractional dispersion did not simulate the tail well. This is in agreement with the higher RMSE errors for the CL derivative operator compared to the Caputo-based approximation which we had discussed in Chapter 3. For Pe < 2 both the compact and the WENO—based fADE solutions were not able to simulate the analytical solution well. We find that as the Pe number increases, there exists a significant difference between outputs from the two models. We decreased the Pe number and found that all the four schemes give a reasonably good fit without any significant difference. The advantage of using the compact based advection for P6 < 2 is that it performs the best in explaining the peak compared to the WENO for various crude grid sizes. This is important since the 57 10 l l I l l 0 Benson et al. (2000) - 0 - Weno-GL - 4 - Weno-Caputo —-¢— Compact-Caputo —°— Compact-6L .L O o l .5 .3 O I 0 Concentration (ppb) 3 s, 10“ , , I I I I I 5 10 15 20 25 30 35 Time (s) Figure 4.4. Plot showing the comparison for a Pa = 1.7 for all the four schemes (Ax=0.25m, D = 0.05m1'8/s,u = Ims’1 , CFL =0.1, t=15S, 01218). All the schemes show reasonably good fit, although the WENO based advection scheme introduced small numerical dispersion tracer studies are performed in rivers over long reaches. A crude grid size will result in a faster simulation which in turn will result in a lesser computational time during parameter estimation. On the other hand, the WENO scheme performs satisfactorily even for high Pe number without oscillations although significant dispersion errors are introduced for crude grid sizes. We conclude that for simulating advection-dominated systems, the grids should be refined such that the condition Pe < 2 is satisfied. 4.5.2 Low Peclet Number: Dispersion-Dominated Systems In order to understand whether Caputo or GL based fractional dispersion scheme explains the dispersion better for a spill case, we ran simulations for another limiting case, i.e. high dispersion. For a high dispersion situation , there will be no significant difference between the advection schemes because the physical dispersion far exceeds any numerical dispersion in WENO or compact schemes. A test case for D = 10m1'8/s, u = 1m/s with Fe = 0.0574 was used. We found that using the GL fractional diffusion, the peak 58 was higher than the analytical solution, while the Caputo based fractional diffusion was able to simulate the analytical solution (both tail and peak) accurately. Both WENO and compact based advection schemes were able to explain the analytical solution well with almost no apparent difference as expected. Therefore, the moving front was highly smeared (due to physical dispersion), so the effect of any numerical smearing was negligible in comparison. It was also observed that as one approaches lower Pe number,the GL based fractional dispersion was unable to explain the peak. This is consistent with the results in chapter 3 (Figures 3.6 and 3.7) where we demonstrated that the Caputo derivative provides a superior description of fractional dispersion. We found that the GL dispersion peaks and tail concentration were higher compared to the analytical solution as shown in Figure 4.5 below 0.45 - 0.4 '- 0.35 - - + Benson et al. (2000) 0.2 - - -Compact-Caputo 4 — Compact-6L 0.15 .. Concentration (ppb) 3 (II I 0.1 '5 0.05 .. 0 I I I I I I I I I 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (s) Figure 4.5. Breakthrough curves for a low Pa = 0.057 (A1: = 0.5, CFL = 0.1, a = 1.8, fl = 1, D = 10m1'8/s, u = lms‘l) at a distance of 5m from the spill location. In order to confirm whether significant errors are introduced due to OS, a first order up- wind based advection scheme using GL and Caputo was used to solve the fADE without any OS using the same parameters (see Figure 4.6). As observed from figure 4.6 both OS 59 and the no—OS approaches gave similar results, which implies that the OS does not intro- duce any significant errors. This implies that Caputo based fractional-in-space provides a superior description compared to CL, supporting the earlier results in Chapter 3. 0-05 I I I I I I 1 I I + Analytical 0 045 _ No Operator-6L ' - - - No Operator-Caputo 0.04 - - I 0.035 - - l E . ' 3 0.03 - - 5 : E 0.025 - - E I 8 I 5 0.02 - - o l l I 0.015 b - l I 0.01 b .. I 0.005 - - 0 I I I I I I I I I 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (s) Figure 4.6. Breakthrough curve (without OS) for a low Pa = 0.057 (Ax = 0.5m, CF L = 0.1, a = 1.8, fl = 1, D = 0.05m1'8/s,u = 1m/s) at a distance of 5m from the spill location. 4.5.3 Effect of fractional derivative exponent(a) To understand the effect of the fractional derivative exponent a with P3 number less than 1, we plotted the breakthrough curves for three different values of a. It was found that as 0: tends towards 1, the fractional dispersion term becomes similar to the advection term. It was found that both caputo and GL based fractional dispersion schemes are not able to simulate the tail. Also, there seems to be a phase error introduced as a approaches 1. The observed difference between the analytical and numerical solutions as a -—> 1.0 is 60 attributed to a change in the nature of the partial differential equation from parabolic (for any value of a > 1.0) to hyperbolic (for a = 1.0). Therefore the point a = 1.0 may be treated as a singularity as far as the validity of the numerical solution is considered. 0-14 I I I I + Benson et al. (2000a) °‘ = 1'9 - - -Compact-Caputo 0.12 - / — Compact-6L - 0.1 " q or = 1.5 Concentration (ppb) 0L=1.1 + - +++*+ "’1- . I I _ h..- ..... ( 0 10 20 30 40 50 Time (s) Figure 4.7. Breakthrough curves comparison between GL and Caputo with Benson et al. [2000] for different values of a with P6 = 0.1 In conclusion, the operator splitting approach was found to produce accurate numerical solutions to the fractional advection dispersion equation as long as the cell Peclet number was sufficiently low (typically less than 2.0). The combination of a compact scheme for advection and the Caputo-based approximation of the fractional dispersion was found to produce the best results (that is, describing the peak as well as the tail in the breakthrough curves). The WENO scheme did not introduce oscillations (not a surprising result con- sidering the nature of this scheme), but significant dispersion errors were introduced at high Peclet numbers. 61 CHAPTER 5 Stream Transient Storage Modeling It has been observed in various studies that the dispersion process traditionally described using Fickian diffusion fails to explain the field measurements. Data collected from nearly 100 streams and rivers shows that the unit-peak concentration tends to attenuate in proportion to the travel time with the 0.89 power, not the 0.5 power predicted by second- order dispersion [Jobson, 2001]. This is indicative of a non-Fickian diflusion process sometimes also called as anomalous diffusion. This can result in a heavy leading and falling limbs in the time-concentration breakthrough curves (BTCs) [Schumer et al., 2003]. The heavy leading limb is indicative of faster-than—Fickian growth rate, usually because of preferential flow paths. On the other hand, heavy falling limb is indicative of the delay of solute in the main channel because of the exchange processes between the main channel and storage zone. In order to understand both these processes, we examine solute transport processes in a stream. Solute dispersion in rivers or streams primarily consists of contributions from four mechanisms: molecular diffusion, turbulent diffusion, shear—flow dispersion and scaling dispersion [Fischer, 1979]. The most important among them being the shear-flow dis- persion, which typically involves mixing due to spatial variability in the flow field. This process produces progressive spreading of the dissolved and suspended substances that causes the variance of the concentration distribution with time. Taylor [1954] showed that after sufficient mixing distance, the variance grows linearly with time and then the concentration distribution follows Fick’s law of diffusion. However, this is not always 62 true because of the substantial water influx across channel boundaries, which results in a far less uniform velocity field than what has been previously assumed while using a second-order dispersion in these systems [Bencala and Walters, 1983; Harvey and Gore- lick, 2000]. Additionally, the exchange of solute in the stream with the hyporheic zone has been shown to produce a consistent delay in the solute transport relative to the main stream flow, leading to a heavy falling limb in the BTCs. On the other hand, the heavy leading or rising limb in the BTC is observed because of the long-range spatial correlation of the dispersion process. This leads to an enhanced diffusion (super-diffusion) than what the second-order diffusion solution estimates [Sahimi, 1993]. It is a result of the infinite- variance particle jump distributions that arise during transport of solute in heterogenous media [Schumer et al., 2001]. There have been many approaches in the past to describe the non-Fickian processes. For example, to model the long-range spatial correlation of the dispersion process, some of these models have used variable dispersion coeflicient, by varying the dispersion coefficient linearly with downstream distance [Berkowitz and Scher, 1995; Wheatcraft and Tyler, 1988]. However, it is difficult to justify this approach, unless there is a some definite trend in the heterogeneity of the stream. There have also been many cases where these adjusted models gave physically unreasonable parameters for the observed field data [Constantz, 1998]. One of the other approaches used in the past is to use the fractional-in-space advection-dispersion equation (fADE) to model the long-range spatial correlation [Deng et al., 2004]. The drawback of using fADE is that it does not account for the transient storage of solutes in the river banks or the hyporheic zones. These exchange processes can be of an arbitrary order depending upon the type of processes occurring. For example, the hyporheic exchange process typically follows a power law residence time distribution (RTD) because of the large waiting times of solute along the subsurface flow paths and due to adsorption. On the other hand, surface storage exchange processes typically follow an exponential law RTD because of comparatively shorter waiting times. The fADE model does not account for these processes which includes adsorption and desorption, and the effects of vegetation and reactions [Deng et al., 2006]. To account for the transient storage 63 of the solute in these storage zones, we use the transient storage (TS) model in this study. 5.1 Transient Storage Model The transient storage model describes solute transport in streams taking the physical, chemical and biological characteristics into account. Most of the previous TS models used the standard advection-dispersion equation (ADE) for solute transport with addi- tional terms accounting for transient storage, lateral inflow and other processes (e.g., decay, sorption) depending upon the scope and complexity of the problem. Generally, the TS model is used in conjunction with the observed field studies to estimate the hydro- logic parameters affecting solute transport. The standard one-dimensional TS model is described using equations 5.1 and 5.2 for the main channel and the storage zones respec- tively [Runkel and Broshears, 1991] ac ac _ 1 a ac qL a7 + “a; — 237.; (ADE) + 7(01 0) + 5‘03 C) ‘5'” 803 A where C = cross-sectional average of concentration in the main stream (ML-3); C; = con- centration in the storage zone (ML—3); C L = lateral inflow solute concentration (ML-3); u = cross-sectional average of water velocity (LT‘l); D = coefficient of longitudinal dis- persion (L2T); a: = space coordinate in the flow direction (L); t = time (T); e = a first- order storage exchange coeflicient (T-l); A = main channel cross-sectional area (L2); A3: storage zone cross-sectional area (L2) and (IL 2 lateral inflow rate (L3T"1L-1). These equations are applicable to conservative (non-reactive) solutes such as tracers, but nonconservative (reactive) solutes may be accounted for by adding chemical reaction terms (e.g., kinetic sorption and decay) to equations 5.1 and 5.2. This conceptual model assumes a first-order mass exchange between the main channel and the storage zones. However, field studies have reported breakthrough curves with heavier power law tails as opposed to the exponential tails modeled using a first-order exchange process [Becker and Shapiro, 2000; Haggerty et al., 2000]. This deviation has been attributed to the longer or deeper hyporheic flow paths [Worman, 2002; Marion 64 et al., 2003] which are generally not described properly using first-order exchange kinetics. Additionally, equations 5.1 and 5.2 assume that the subsurface is a well-mixed reservoir, while in reality hyporheic exchange is generally characterized by a strong exchange due to spatial variations of properties such as permeability in the subsurface [Marion and Zaramella, 2005]. Some researchers have used a fractional-in-time derivative to describe more complex exchange patterns between the main channel and the storage zones in aquifers [Schumer et al., 2003], however, time fractional derivatives is beyond the scope of the current work. In this study, we use the existing TS model assumptions, of describing the heavy-tailed falling limb, using the first-order exchange process between the main channel and the storage zones. Additionally, we replace the second-order dispersion using fractional-in-space dispersion to model the heavy rising limb observed due to faster-than- fickian processes. We refer to this new model as the fractional-in-space transient storage model (FSTS). 5.2 The Fractional-in-Space Transient Storage (FSTS) Model The FSTS model can describe the non-Gaussian rising and falling limb in a BTC by using fractional-in-space dispersion and a first-order exchange process between the main channel and the storage zones respectively. The FSTS model is based on fADE instead of ADE and can be rewritten as 3C BC 60C q L -3—t+u—3—1_:-:D5x—a+71—(CL—C)+E(Cs—C) (5-3) 3C5 A where a is the fractional derivative exponent and fl is the skewness parameter that controls the bias of the dispersion. We ran a number of simulations in order to understand the differences between the FSTS and TS models. All the simulations were run subject to D = 1, for a concentration profile skewed forward, a case also studied by [Benson et al., 2004] to represent preferential transport along pathways such as fractures and macr0pores. We simulated the injection of a tracer in a stream, subject to the following initial and 65 boundary conditions C(O, 0) = co and C(z,0) .—. 0,v x(x 75 0) and t = o (5.5) -D (0.2% + (1 — mW—SB—C; ) _+_ = 0 30-10 804:, °° for o g t < 00 (5.6) Figure 5.1 (the figure is a spatial snapshot of concentration at a given time) shows that higher concentration is observed for a < 2 compared to a = 2 at larger downstream distances in the main channel at t = 6003. This is due to the faster-than-fickian process discussed earlier, which results in solutes moving faster than the standard second-order dispersion based TS model. 10° - - _1 cl=1.4 ;\ 10 / (131.8 10" Concentration (ppb) 10‘ 104 I I I I I I I I 10 100 200 200 300 400 500 600 700 800 Distance , x(m) Figure 5.1. Numerical model showing the snapshot of concentration vs distance in the main stream at t = 600s spatially for A = 40 m2, A, = 4 m2, .0 = 1m1-8/s, u = 1ms’1, e = 0.0 s—l, M 2: 1000g. The point of injection is at a: = 10 m We also ran the FSTS model with a = 2.0, to understand how the concentrations in the stream and the storage zones change with time. In Figure 5.2 we plot the distribution of concentration in the main channel and the storage zone for a = 2 and storage exchange coefficient, 6 = 0.001. We can observe that the major contribution to the rising limb of 66 the BTC comes from the main channel (stream). However, the opposite is true for the falling limb of the BTC (Figure 5.2 log scale). This is because of the large waiting times of the solute in the storage zones before releasing back into the main channel. 0.9 . . 10° 03 I. Stream .. Stora e Zone _ 0 7 _ 9 _ 10 2 . Total Cons A . (Stream + A '8, 0-6 ' Storage '8, 3 Zone) 3 10" - g 0.5 - - g .2. j 5' C 0.4 - C 0 o -0 5 51° ' o 0.3 o 0.2 10" . 0.1 -10 l 0 10 0 200 400 600 1 o“ 1 0° 1 05 Time since release (s) Time since release (s) Figure 5.2. Breakthrough curve for the stream and storage zone concentration at a distance of a: = 50m from the point of injection in (a) linear scale and (b) log log scale. In log scale, it can be observed that for the falling limb the major contribution of the concentration comes from the storage zone whereas, for the rising limb it is the main channel. (A = 40 m2, A, = 4 m2, D = 1m1-8/s, a = 1, u = 0.5ms-1, e = 0.001 s-1,M = 1000 g, a = 2) We also compared the FSTS model with the TS model in order to compare the behavior of the BTCs in the main channel for these models. We ran the model for the parameters mentioned in Figure 5.2 and observed that the heavy falling limb for both of these models gave almost similar concentration while significantly higher concentration for the rising limb in case of FSTS model. This reinforces our earlier discussion that both these models uses first-order exchange processes between the main channel and the storage zone, which controls the falling limb, and therefore give similar concentration magnitude. On the other hand, the rising limb is controlled by the fractional-in-space dispersion, therefore 67 the FSTS model gives a heavier rising limb compared to the TS model. 1° I I I I i I —FSTS(01=1.5) o TS 104 .....'...o o -2 o .' A 10 0 " E. o E o - o E 10 ° 0 é o 8 ..- 10 o o 10" o . o 10" I I l I l I 10“ 1031 10” 1o” 10” 10” Figure 5.3. FSTS vs TS model breakthrough curve in the main channel for A = 40 m2, A, = 4 m2, D = lml's/s, e = 0.0001 ms-l, M = 1000 g, a = 2.0, A2: = 2.9950 111 and a: =1000 In 5.2.1 Comparison of the FSTS model with the Analytical Solution for a = 2 We used the analytical solution by De Smedt et al. [2005] to evaluate the accuracy of our FSTS model for the special case of Fickian diffusion, i.e., a = 2. There exists no analytical solution for the FSTS model, hence we compared our numerical model for a = 2, which effectively reduces it to a TS model. Analytical solution for the TS model exists and is given by De Smedt et al. [2005] for an instantaneous injection of a tracer in a river with uniform flow in an “infinite” domain. The analytical solution to equations 5.1 and 5.2 is given by [De Smedt et al., 2005] t 2.2—11272 1 e t—T. s + ———2— — — E J 8’)", C(12, t) = f 40" 2? ) ( ) C1(:I:,T)d1' (5.7) 0 —EJ (SEED, 87') where B = A3/ A and C1(a:,t) is the classical solution to the ADE for an instantaneous 68 slug release and is given by: M (x — 'ut)2 C1(:l:,t) = m exp— (W) (5.8) The J function in the above equation can be evaluated using the following relation 12. n-l a 00 m J(a, b) = 1 — e—b/e—AIO (2M) d/\ = 1 — e_a-b Z a. (rln-l (5.9) n=1 0 ' F o where 10 is the modified Bessel function of zero order. One of the drawbacks of the above analytical solution is the problem of convergence for higher values of 5 due to the large computational time required to evaluate the J functions. Therefore, if this equation is used for parameter estimation it will increase the computational time significantly which makes the numerical model more attractive. The above analytical solution was also found to agree well with other numerical models like OTIS [Runkel, 1998]. The FSTS model was found to be in excellent agreement with this analytical solution for a = 2, for different values of e, D, and u as shown in Figures 5.4, 5.5, and 5.6 respectively. 5.3 Application of the FSTS Model to Describe the Field Data We applied the FSTS model to describe tracer transport in two Michigan streams. One of the streams had a large reach length (The Grand River, Reach Length - 42 km) while the other had a smaller reach length (The Red Cedar River, Reach Length - 5 km). In order to apply the FSTS model to the streams, we first confirmed whether the BTCs of the observed tracer concentration for these sites exhibited non-Gaussian behavior during early or the late time. This was done by plotting the concentration on a probability scale against the Z—scores of the travel time. A deviation from the straight line is indicative of the non-Gaussian behavior. One can observe from the Figure 5.3, that the rising limb for the Reach C of the Red Cedar River shows a deviation from the straight line, which indicates non-Fickian early breakthrough. This justifies the use of fractional-in-space dispersion based TS model to describe the field data for the site. 69 e=0 -———' Numerical ' 0 Analytical 8 = 0.001 0.15 1- 0.1 0.05 Concentration (ppm) Concentration (ppm) 10- ' . I I I ‘ . 2000 2500 3000 3500 4000 4500 Time since release (s) Figure 5.4. Comparison of compact-Caputo based fADE with the analytical solution given by De Smedt et al. [2005] for Q = 16 mas-1, A = 40 m2, A, = 4 m2, D = 1 mas’l, M = 1000 g, a = 2.0, Ax = 0.33 m and x = 1230 m 5.3.1 The Red Cedar River, Michigan The Red Cedar River (RCR) is a fourth-order stream in south central Michigan. It originates as an outflow from Cedar lake, Michigan, and flows into East Lansing and Michigan State University (MSU). The RCR meanders through the MSU campus over a stretch of 5km (our study reach) and has an average slope of 0.413 m/km. The study reach was bounded by Hagadorn Bride on the East and the Kalamazoo Street Bridge on the West (Figure 5.8). Out of the many tracer studies done on this river in 2002 we used the tracer test conducted on 19th March, 2002 [Phanikumar et al., 2007] with a discharge rate of Q = 19.89 m3/s. Fluorescein dye was injected at the Hagadorn Bridge and samples were collected at downstream locations at the Farm Lane (1400m), Kellog (3100m) and Kalamazoo Bridge (5079m), respectively (Figure 5.8). The total mass injected was equal to 10.70 kg Phanikumar et al. [2007]. Since it was a small reach, the influence of the lateral inflow (contribution from tributaries and baseflow) was negligible [Phanikumar 70 — Numericalupresent work) 0 Analytical(De Smedt et al.,2005) _ 12 - <—o=0.25mzls A 10 - q E D. 5 _§_ 8 " D=0.5mzls ' é . 8 D=1m Is 5 O 100 150 Time since release (s) Figure 5.5. Comparison of compact-Caputo based fADE with the analytical solution given by De Smedt et al. [2005] for Q = 20 m3s‘1, A = 40 m2, A3 = 4 m2, 6 = 0.0001 ms’l, M = 100 g, a=2.0,Ax=0.33mandx=10m et al., 2007] and therefore, we used q L = 0 for all the simulations. Phanikumar et al. [2007] applied the TS model on a reach-by-reach basis for the three reaches A, B and C as shown in Figure 5.8. They used a global optimization algorithm to estimate parameters for each of the three reaches (total of 15). They were able to estimate parameters for the reaches A and B , but for the reach C parameters were found to suffer from singular convergence [Runkel, 1998; Fernald et al., 2001], likely due to competition between processes. Our FSTS model uses a constant D and variable a values for each of the reaches, to capture the long-range correlation of the dispersion process. Using this approach we were able to describe the scale dependency of the dispersion process and better constrain the model parameters. We estimated the a value on a reach~by-reach basis to account for the intra-reach variability that can not be captured using a constant a. In the FSTS model we estimated the six parameters (i.e.,A, A3, 9L: 5, D and (1) instead of five for the TS model for each of the reaches because of the introduction of an extra parameter, a. For the three reaches A, B, and C where we tested our FSTS model, 71 10 I I — Numerical(Compact-Caputo) 0 Analytical(De Smedt et al.,2005) Concentration(ppm) h -2 1 l 0 50 100 150 Time since release (s) Figure 5.6. Comparison of compact-Caputo based fADE with the analytical solution given by De Smedt et al. [2005] for A = 40 m2, A, = 4 m2, D = 1 mas“1, e = 0.0001 mS—l, M = 100 g, a=2.0,Ax=O.33mandx=10rn the total number of parameters were equal to 16 (5 for each reach and a constant D). All the parameters in the FSTS model were estimated using global optimization procedure based on genetic and pattern search algorithms as implemented in MATLAB [Goldberg, 1989]. Parameter estimation was done by minimizing the root-mean-square error (RMSE) between the observed data and the simulated values. We found that the FSTS model (see Figure 5.9) was able to describe the observed tracer concentrations satisfactorily for all the three reaches. A constant D = 4.14 was used to describe all the three reaches A, B and C. The a values decreased with downstream distance from 1.98 to 1.75. This was because of the increase in the magnitude of deviation from the Fickian diffusion as the plume travels a larger downstream distance. The higher values of the (As/A) with downstream distance are attributed to the alluvium storage and the sediment characteristics (gravel and coarse sand) in reach C. Table 5.1 gives the estimated values for all the three reaches with the normalized RMSE error between observed and simulated values. 72 I F I I I 1 ' Red Cedar River: Reach c J” " 0.8 - . 0.5 - - O Q U 0.4 - - 0.2 e . Fickian o - l I I l I - -0.75 -0.5 -0.25 0 0.25 0.5 0.75 (t—tlh/t—f Figure 5.7. Probability plot showing the normalized concentration (or probability) as a function of the Z—scores of travel time for the Reach C of the Red Cedar River. Deviation from the straight line indicates a deviation from Fickian diffusion 5.3.2 The Grand River, Michigan We also tested the FSTS model on the Grand River in order to see whether it was able to simulate the scale-dependent dispersion process on a significantly larger reach of 42 km compared to the 5 km reach for the Red Cedar River. This study reach extended from the city of Grand Rapids to the town of Coopersville. The data and analysis based on the TS model (i.e., with a second-order dispersion term) are described in Shen et al. [2007]. The surficial geology of the Grand River Basin is dominated by rivers crisscrossing the moraines and outwash plains formed by extensive glaciation during the Pleistocene. Sampling was carried out at four downstream sites / bridges from the point of injection at Ann Street Bridge near downtown Grand Rapids (see Figure ??). The downstream sampling locations were Wealthy Street (Site 1), 28th Street( Site 2), Lake Michigan Drive (Site 3) and 68th Street (Site 4). The study reach was sufficiently long to make watershed influences important Shen et al. [2007], hence we accounted for the lateral inflow in the FSTS model. The average discharge for the river during the tracer study measured by 73 / . «535359134 Legend 7‘ 1 , —Streams N ‘- If, -.::c-b I“ 2 H ----- ' 13:] Lakes in — In ‘ 3 ¢ 11.. S Michigan 2 I P . I— O 5 ‘ S" o 5 10 20 9 E Kilometers I I 34'35'0'W 84’1 O'O'W I— c : B i A 1 Mb” {—— Flow 2 I Kalamazoo Ubrary _ g ue Hagadorn . Farm Lane 309 5 0.3 0.6 Kilometers — I I M'30'0'W 84'20'40'W “'ZT'SG'W Figure 5.8. Red Cedar River Watershed showing the Red Cedar River and the sampling locations (adapted from Phanikumar et al. [2007]) the USGS streamflow gaging station was 91.43 cubic meter per second. Further details can be found in Shen et al. [2007]. In order to correctly estimate the mass of the tracer injected in the main stream channel, we iterated the BTC for reach A and used it for estimating FSTS model parameters for the reaches 8,0, and D. The parameters for this site were estimated similar to the Red Cedar river. The RMSE error and the estimated parameters for the reaches is shown in Table 5.2. A constant D = 3.32 was used to describe all the three reaches A, B and C. The a values decreased with the downstream distance from 1.99 to 1.62 similar to the Red Cedar river indicating an increase in the magnitude of deviation from the Fickian diffusion with downstream distance. 74 70 I I I I I 0 Observed 60 _ : Rae" A FSTS-Simulated _ , / a =1.98 50 - l - E 5 40 - . Reach B a =1.83 " C .2 *5 30 - - *5 Reach C on =1.75 0 § 20 - - U 10 - '- V \ 0 ‘ - _10 I I I J I 1 2 3 4 5 6 Time (hrs) Figure 5.9. Fluorescein breakthrough curves for the three reaches A, B and C for the Red Cedar River with a. constant dispersion coefficient D = 4.14 using the FSTS model 75 Grand River Watershed Michigan . r t N r + i """ e R Z . i :c .‘-= . :9 68th St. Bridge *3 8 Ann St. Bridge i; H‘s/”x "...\ “~- USGS ”19000 Lake Michigan DrNe k, in}, z { ySt. - . ‘ =6 “ , r . I? .3 ‘r: I“. ..i “‘"u E:- - .0“ Legend ‘-"":j__':» muse '- .. 1- 3,. v I injection Site "52‘ _ 7 I i' ‘g A USGS Gage ' if - 0 Sampling Sites ‘w‘ 34.. ........... z — Grand River 5...; =9 ff? Catchments -.,, ' ' 2' N V 86°0'0"W 85°48'0"W 85°36'0"W Figure 5.10. Grand River Watershed and study region showing the sampling sites (adapted from Shen et al. [2007]) 76 cl Q .1 -| .i l 0 Observed 16 - Rm" 3 FSTS-Simulated - 14 - - .... 12 - .. .D E v 10 - Reach c '- 5 a=1£9 3 g 3 ' 0 ~ Reach D ' E a = 1.62 g 6- r o = 1.69 o 4 l- O a -i ‘e \ 2 - - o 3 ' l. _2 I I I I 0 5 10 15 20 25 Time (hrs) Figure 5.11. Rhodamine WT breakthrough curves for the three reaches B, C and D of the Grand River with a constant dispersion coefficient D = 3.32 using the FSTS model (data taken from Shen et al. [2007, in review]) 77 E3 mg as 3m 3 m: is 2.: 5% a $3. woe was as so a: 3s NE. mm? m Se mod 9: new as a: Es Se $3 < $505 Alev Aaxwwmgv Nadav HNEV ANSV $sz a 3% ...L: xt as a Q as w flag E58 mBmm mama: 63d 380 com e5 8% Emuogfida ween omfiopm Ad 03¢? 78 a: a: as? E SE a: New as 3.2: o as; one so 3 .83 a: mom 3. an: 0 EB. m3 a; so 83 a: N3 3.. ENE m 3.505 AHlmv A8x\m\mmv Am\dav dwav ANEV $382 e is. are xm as a Q as w ease _ocoa mfimm m5 wfims Q £53980 a HEB 8:32 $53 2: 8.“ 33m seesaw 2: .8.“ mampmfieamm BEN owgopm .Né Each. 79 CHAPTER 6 Conclusions The major conclusions of this work are summarized below. 1. The fractional derivative based on the Caputo definition was found to provide a superior description of the fractional-in-space dispersion for the case of instantaneous slug release in a stream. In the past, space-fractional Caputo derivatives were shown to provide a better description (compared to the Griinwald-Letnikov definition) for continuous slug release in a stream; however, this issue was not examined in detail for the case of instantaneous release. The superior performance of the Caputo derivative is likely due to the distribution of the weights as shown in chapter 3. 2. Of the two advection schemes examined in this work, the fourth order accurate compact scheme was found to provide better solutions. The WENO scheme was found to introduce phase and amplitude errors for large wave numbers. Although these schemes have a Courant number restriction, it was not too stringent for the applications considered in this work since constant reach-averaged velocities are used in the models. For the more general case in which the velocities are spatially non-uniform, schemes that do not have a Courant number restriction may be more attractive (e.g., Lagrangian or semi-Lagrangian schemes). 3. We have successfully implemented an operator splitting approach for solving the governing equations. The approach was based on the Strang OS method and the advection and dispersion equations were solved using separate numerical methods. 80 Comparisons shown in chapter 4 indicate that the OS scheme produced excellent results. . One of the objectives of this work was to test whether improved dispersion modeling will enable us to better constrain the TS model parameters such as the storage zone exchange rate. The motivation was provided by the fact that in some stream reaches dispersion was not needed to describe the observed breakthrough data while in other reaches the classical ADE was found adequate. Examination of tracer breakthrough in the Red Cedar River showed that in one of the reaches (Reach C), the early breakthrough was clearly non-Fickian which indicated that a fractional order derivative was needed to describe the transport. Application of our fractional transient storage (FSTS) model showed that one dispersion coefficient but different a values described the observed data well without resulting in any false convergence problems noted earlier for this reach using a standard TS model. . Higher order approximations of the fractional derivative are attractive for several reasons (e.g., more accurate solutions on a relatively crude grid means less computa- tional effort while estimating parameters). In this work we described one framework that can be used to obtain higher order accurate solutions based on the Caputo def- inition of the fractional derivative by evaluating the integer order derivatives to higher order accuracy using a compact scheme. 81 APPENDIX A Caputo Derivative Operator The one-dimensional fADE equation to simulate the solute movement in surface and sub- surface system without the source term is given by: at +u$=D 3—0 8" (5:S+<1—B>—affga) (A1) The finite difference method described below is based on the method preposed by Zhang et al. [2006] and modified for an infinite domain [—L, L]. The finite difference form of equation A.1 can be written as _z—Et—i'l' A A (A2) t 05+ t—Cifl t _ Qt+At/2 _ Qt+At/2 “ A1: — i—1/2 i+1/2 where the fractional dispersive flux can be approximated using the Caputo derivative given by: a: i—1/2 Q. = _:D__ f 1 diy (A 3) 2—1/2 m2 _ a) (xi—U2 _ war—1 3y ' “5241/2 —D 1 BC Qi-l—l/2 : F(2 _ (I) / ($i+1/2 .. y)a_1‘a—;dy (A.4) Assuming that C in [:c — zi+1/2,x — xi_ 1/2] is linearly distributed, the concentration gradient aC/Br can be approximated by [C(x — $i+1/2) — C(x — :rz-_1 /2)] /Az. Sub— 82 stituting it in equation A.3 to get (22-4/2 as i-l 3:0 22 I‘(2—a) + .53 which is equal to D Z (Ci—j -Cz'—j—1) f N—z' J=0 xi_1/2 fl0,($i—1/2—y) a—l dy y —L L-z- + 2—1/2 __ 0' , f (1 a) (ac,_1/2+y)dy xj+1/2 31-1/2 $j+1/2 (Cm — Ci+j-1) f 31—1/2 Qi—l/2 = _I‘(2 — a)A:c°‘ where w = 0+1)“ — (212—0 Similarly the flux Qi+1/2 can be approximated as F D Qi+1/2 = —F(2 — a)Aa:a i 1.230 'Bwj(0i+1—j ‘ Ci-j) N—i—l + 1.120 (1_’6)wj(Cz'+j—l-1— 1— ya" fidy C. 2+j) (A.5) (A.8) (A.9) In equation Al the advection term can be approximated using central differencing given by [Zhang et al., 2006] 3_C = Ci+1 - Ci—l 8:1: 2A$ 83 (A.10) Substituting equation A10 in equation A.1 the finite-difference form of fADE is given by: t+At t t+At_ t+At C. —C. C. C._1 2+1 2 ‘J—A—H + “ 2A2: ' t+At t+At t+At ‘ (Ci—_l ‘26} +Ci+1 ) 1. t t +j§13wj(cz'+1-j — Ci—j) i—l . t t A.11 D _J=1flw](0z J —Ci_j_1) ( ) = I‘l2—crlAzZ2 l > N—z—l t t + 3'21 (1 — fi>wj(0i+j+1 ‘ Ci+j) N—z' t t ‘ 3.21 (1 ‘ mijifi ‘ Ci+j—1) t = J The above equation is a semi-implicit form, since the Gaussian dispersion terms (RHS -first three terms) and the advection terms are solved implicitly while the non-Gaussian terms (remaining RHS terms) are solved explicitly. One can also approximate the advec- tion term to a higher order accuracy using either a WENO or compact scheme discussed in Chapter 2. A prescribed concentration boundary or a prescribed-flux boundary with inlet boundary at x = —L and outlet boundary at a: = L can be used to solve the above equa- tion. The prescribed-flux boundary commonly encountered in hydrology can be modeled as a prescribed flux at a: = —L and a free drainage at the outlet :1: = L. The boundaries can be modeled as 30—10 30-10 _ uC' — D (3W + (1 - B)6(—$)Q"I)$=_L — UC—L (A 12) -1 —l _D (fl—TZa—C + (1 -— B)———I8f’:)a€ x = 0 C'_ L represents the concentration at the inlet boundary. This approach uses an integer- order flux boundary instead of fractional-dispersive boundary given by [Zhang et al., 2006] D 311:3(C;-N+1 " C—N) Q—N+1/2 = ‘W + £0 (1 — fl)(C_N+1+j - C—N) N_13— (A.13) D Z fiwflCN—j _ CN-j”1) QN—1/2 = “W i=0 +(1- [3)(CN—1 — CN) 84 APPENDIX B Griinwald-Letnikov Derivative Operator The value of a fractional derivative Operator (based on Griinwald) acting on the function C(x, t) is an infinite series given by [Oldham and Spanier, 1974] (900 1 N—l a) . 3335 Msmm 2:0“ p0- ——+—j1)0($-Jh.t) (3.1) where I‘()=gamma function, Ax = :r/N; N = positive integer. Oldham and Spanier [1974] presented another definition to the above definition which has a better convergence property and is given by: N— 1F BO‘C_ 1 (j — ___a) , In analogy, the backward-finite difference form of the above equation is given by Griiwald- Letnikov by induction N—l 800 1 a _= ' _— c _ A , 313a N133» Axum—a) ‘72—: (j) (CE J :r,t) (B 3) where the coefficients a over j = weighting factors, reflect the length of the memory of the fractional derivative and is given by [Podlubny, 1999]: —1—a wq_ .7 ma ,7 ——J—,——w j“ _1 , 1118 =1, j=1,2,3... (B4) The coefficients of the equation BI, 82, and 8.3 are equivalent, and can be expressed as aac~ axa Naia :0 Mom] (13.5) jzo 85 which is the Griiwald-Letnikov definition while N+1 8‘10 1 t 330: z Am“ 2 w?CN+1—j (B.6) j=o is the Deng-Singh-Bengtsson definition [Deng et al., 2004]. 86 APPENDIX C Griinwald—Letnikov and Caputo Weights We estimate the normalized weights of the CL and Caputo derivative operators to compare the relative influence of the functional values at various spatial locations for the derivative operators. While the CL derivative is approximated using the shifted GL derivative operator, the Caputo weights are approximated using a finite volume approach. The shifted GL derivative is given by: N+1 as: R. A? Z “’J'CNH—a' (0'1) i=0 where a —1 . = 1__ -_ , :1, . 21,2,3". C-2 ui7 ( 2' )w'7 1 W0 J ( ) In case of the CL, since a finite-difference approach is used to approximate the derivative, the coefficients of the series given by equation 0.1 is equal to the weight itself and is given by equation 0.2. In case of Caputo fractional derivative, the derivative is defined using a finite volume approach given by: OCDgC = Qi—l/Z ‘ Qi+1/2 ((13) where the flux Q2._1 /2 is given by: z—l D Qi-1/2 : ‘r(2 — ammo: 3&2..ch ‘ Ci—j—I) (Q4) and the flux (2241/2 is given by: D ' C Qi+1/2 = —[‘(2 _ 0,an Z wj(Cz'+1—j — Ci—j) (G5) 87 The weights to? are given by: w; = (j + I)“ — (2')“ (0.6) for 1 S a < 2. Substituting the flux values in equation 0.3, we get D i z-l F(2 —— a)Axa Z w§(0i+1-j - Ci—j) - Z w§(0i_j - Ci—j—l) ((3.7) .=0 j==0 The above equation can be rewritten and simplified as F i - CH1 _ 202' + Ci—l + .2 w§(Cz'+1—j — Ci—j) 01300 = 0 3:1 0 “3 I‘i2—alAaza i—l " .2 w§(Cz'—j “ Ci—j—l) L 3:1 , 1 q 1— C CH1 ‘ 202' + 01—1 + 3.5211” '(Ci—H—j — gci-j + Ci—j—l) __ D — I‘l2—alea J +wf(Cl — CO) Collecting the common terms from the above equation and rearranging it, we get 11180241 + (-2w8 + w‘f)Ci + (1118 -— 2112? + w§)Cz-_1+ D (“If — ng + w§)Cz-_2 + (112% - 210% + wfi)Cz-_3+ C C C C C 630%: = (0.9) We replaced constant 1 by 1118 in order to have a consistent subscript notation. 88 BIBLIOGRAPHY S. Abarbanel and A. Kumar. Compact high-order schemes for the euler equations. Journal of Scientific Computing, 3(3):275—288, 1988. J. David Allan. Stream ecology : structure and function of running waters. Chapman & Hall, London ; New York, 1st edition, 1995. J. David Allan.ill. ; 26 cm.Includes bibliographical references (p. 343-37 7) and indexes. M. W. Becker and A. M. Shapiro. Tracer transport in fractured crystalline rock: Evidence of nondiffusive breakthrough tailing. Water Resources Research, 36(7):1677—1686, 2000. K. E. Bencala and R. A. Walters. Simulation of solute transport in a mountain pool-and- riflie stream: A transient storage model. 19(3), 1983. D. A. Benson. The fractional advection-dispersion equation: Development and application. Dissertation of doctorial degree, University of Nevada Reno, 1998. D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert. Application of a fractional advection-dispersion equation. Water Resources Research, 36(6):1403—1412, 2000. D. A. Benson, C. Tadjeran, M. M. Meerschaert, I. Famham, and G. Pohll. Radial fractional-order dispersion through fractured rock. Water Resources Research, 40(12): -,2004. B. Berkowitz and H. Scher. On characterization of anomalous-dispersion in porous and fractured media. Water Resources Research, 31(6):1461—1466, 1995. D. B. Booth. Challenges and prospects for restoring urban streams: a perspective from the pacific northwest of north america. Journal of the North American Benthological Society, 24(3):724—737, 2005. A. Burgel, T. Grahs, and T. Sonar. From continuous recovery to discrete filtering in numerical approximations of conservation laws. Applied Numerical Mathematics, 42 (1-3):47—60, 2002. A. Cadiou and C. Tenaud. Implicit weno shock capturing scheme for unsteady flows. application to one-dimensional euler equations. International Journal for Numerical Methods in Fluids, 45(2):197—229, 2004. M. Caputo and F. Mainardi. New dissipation model based on memory mechanism. Pure and Applied Geophysics, 91(8):134, 1971. 89 Ryuichi Chiba, Sergei Fomin, Vladimir Chugunov, Toru Takahashi, Yuichi Niibori, Toshiyuki Hashida, Kazuyuki Tohji, Noriyoshi Tsuchiya, and Balachandran Jeyade- van. Numerical simulation for non-fickian diffusion into fractured porous rock, 2006. J. E. Cohen. Human population: The next half century. Science, 302(5648):1172—1175, 2003. J. Constantz. Interaction between stream temperature, streamflow, and groundwater exchanges in alpine streams. Water Resources Research, 34(7):1609——1615, 1998. B. Costa and W. S. Don. High order hybrid central - weno finite difference scheme for conservation laws. Journal of Computational and Applied Mathematics, 204(2):209—218, 2007. D. J. Dangelo, J. R. Webster, S. V. Gregory, and J. L. Meyer. Transient storage in ap- palachian and cascade mountain streams as related to hydraulic characteristics. Journal of the North American Benthological Society, 12(3):223—235, 1993. P. M. Davis, T. C. Atkinson, and T. M. L. Wigley. Longitudinal dispersion in natural channels: 2. the roles of shear flow dispersion and dead zones in the river severn, uk. 4 (3);355—371, 2000. F. De Smedt, W. Brevis, and P. Debels. Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage. Journal of Hydrology, 315:25—39, 2005. A. O. Demuren, R. V. Wilson, and M. Carpenter. Higher-order compact schemes for nu- merical simulation of incompressible flows, part i: Theoretical development. Numerical Heat Transfer Part B-Fundamentals, 39(3):207-230, 2001. Z. Q. Deng, V. P. Singh, and L. Bengtsson. Numerical solution of fractional advection- dispersion equation. Journal of Hydraulic Engineering-Asce, 130(5):422—431, 2004. Z. Q. Deng, L. Bengtsson, and V. P. Singh. Parameter estimation for fractional dispersion model for rivers. Environmental Fluid Mechanics, 6(5):451—475, 2006. A. Dolezal and S. S. M. Wong. Relativistic hydrodynamics and essentially nonoscillatory shock capturing schemes. Journal of Computational Physics, 120(2):266—277, 1995. J. H. Duff and F. J. Triska. Denitrification in sediments from the hyporheic zone adjacent to a small forested stream. Canadian Journal of Fisheries and Aquatic Sciences, 47(6): 1140—1147,1990. J. H. Duff, F. Murphy, C. C. Fuller, F. J. Triska, J. W. Harvey, and A. P. Jackman. A mini drivepoint sampler for measuring pore water solute concentrations in the hyporheic zone of sand-bottom streams. Limnology and Oceanography, 43(6):1378—1383, 1998. R. P. Fedkiw, G. Sapiro, and C. W. Shu. Shock capturing, level sets, and pde based methods in computer vision and image processing: a review of osher’s contributions. Journal of Computational Physics, 185(2):309—341, 2003. 90 A. G. Fernald, P. J. Wigington, and D. H. Landers. Transient storage and hyporheic flow along the Willamette river, oregon: Field measurements and model estimates. Water Resources Research, 37(6):1681—1694, 2001. Hugo B. Fischer. Mixing in inland and coastal waters. Academic Press, New York, 1979. C. Gallo and G. Manzini. A mixed finite element finite volume approach for solving biodegradation transport in groundwater. International Journal for Numerical Methods in Fluids, 26(5):533—556, 1998. David E. Goldberg. Genetic algorithms in search, optimization, and machine learning. Addison-Wesley Pub. Co., Reading, Mass, 1989. R. Gorenflo, F. Mainardi, D. Moretti, and P. Paradisi. Time fractional diffusion: A discrete random walk approach. Nonlinear Dynamics, 29(1-4):129-143, 2002. D. Gottlieb and C. W. Shu. On the gibbs phenomenon and its resolution. Siam Review, 39(4):644—668, 1997. S. Gottlieb, D. Gottlieb, and C. W. Shu. Recovering high-order accuracy in weno com- putations of steady-state hyperbolic systems. Journal of Scientific Computing, 28(2-3): 307—318,2006. R. Haggerty, S. A. McKenna, and L. C. Meigs. On the late-time behavior of tracer test breakthrough curves. Water Resources Research, 36(12):3467-3479, 2000. E. L. Harbott, M. R. Grace, J. A. Webb, and B. T. Hart. Small-scale temporal variation and the effect of urbanisation on extracellular enzyme activity in streams. Journal of Environmental Monitoring, 7(9):861—868, 2005. A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy. Uniformly high-order accurate essentially nonoscillatory schemes .3. Journal of Computational Physics, 71(2):231-303, 1987. C. Harvey and S. M. Gorelick. Rate-limited mass transfer or macrodispersion: Which dominates plume evolution at the macrodispersion experiment(made) site? 36(3):637— 650,2000. M. Yousuff Hussaini, B. Van Leer, and J. Van Rosendale. Upwind and High-Resolution Schemes. Springer—Verlag Telos, 1997 . G. S. Jiang and C. W. Shu. Efficient implementation of weighted eno schemes. Journal of Computational Physics, 126(1):202—228, 1996. H. E. Jobson. Predicting river travel time from hydraulic characteristics. Journal of Hydraulic Engineering-Asce, 127(11):911—918, 2001. J. B. Jones and R. M. Holmes. Surface-subsurface interactions in stream ecosystems. Trends in Ecology 63 Evolution, 11(6):239—242, 1996. 91 K. Jonsson, H. Johansson, and A. Woman. Hyporheic exchange of reactive and conser- vative solutes in streams - tracer methodology and model interpretation. Journal of Hydrology, 278(1-4):153—171, 2003. L. A. Khan and P. L. F. Liu. An operator splitting algorithm for the three-dimensional advection-diffusion equation. International Journal for Numerical Methods in Fluids, 28(3):461—476, 1998. L. A. Khan and P. L. F. Liu. An operator splitting algorithm for coupled one-dimensional advection-diffusion-reaction equations. Computer Methods in Applied Mechanics and Engineering, 127(1~4):181-201, 1995. P. Kumar and O. P. Agrawal. An approximate method for numerical solution of fractional differential equations. Signal Processing, 86(10):2602—2610, 2006. J. M. Lees, L. A. Camacho, and S. Chapra. On the relationship of transient-storage and aggregated dead zone models to solute transport in streams. 36(1):213, 1999. S. K. Lele. Compact finite-difference schemes with spectral-like resolution. Journal of Computational Physics, 103(1):16—42, 1992. Ming Li. Numerical solutions for the incompressible Navier-Stokes equations. PhD thesis. F. Liu, V. Anh, and I. Turner. Numerical solution of the space fractional fokker-planck equation. Journal of Computational and Applied Mathematics, 166(1):209—219, 2004. X. D. Liu, S. Osher, and T. Chan. Weighted essentially nonoscillatory schemes. Journal of Computational Physics, 115(1):200—212, 1994. S. L. Lowery and W. C. Reynolds. Numerical simulation of a spatially-deveIOping mixing layer. Technical Report Rep. TF-26, 1986. C. Lubich. Discretized fractional calculus. Siam Journal on Mathematical Analysis, 17 (3):704—719, 1986. V. E. Lynch, B. A. Carreras, D. Del-Castillo—Negrete, K. M. Ferreira-Mejias, and H. R. Hicks. Numerical methods for the solution of partial differential equations of fractional order. Journal of Computational Physics, 192(2):406—421, 2003. 747ZX Times Citedz21 Cited References Count:18. G. I. Marchuk and N. N. Yanenko. Solution of multidimensional kinetic equation by splitting method. Doklady Akademii Nauk Sssr, 157(6):1291—1310, 1964. A. Marion and M. Zaramella. Diffusive behavior of bedform-induced hyporheic-exchange in rivers. Journal of Environmental Engineering-Awe, 131(9):1260—1266, 2005. A. Marion, M. Zaramella, and A. I. Packman. Parameter estimation of the transient storage model for stream-subsurface exchange. Journal of Environmental Engineering— Asce, 129(5):456—463, 2003. 92 M. M. Meerschaert and C. Tadjeran. Finite difference approximations for fractional advection-dispersion flow equations. Journal of Computational and Applied Mathe- matics, 172(1):65—77, 2004. J. L. Meyer, M. J. Paul, and W. K. Taulbee. Stream ecosystem function in urbanizing landscapes. Journal of the North American Benthological Society, 24(3):602—612, 2005. Kenneth S. Miller and Bertram Ross. An introduction to the fractional calculus and fractional diflerential equations. Wiley, New York, 1993. E. W. Montroll and G. H. Weiss. Random walks on lattices .2. Journal of Mathematical Physics, 6(2):167—175, 1965. P. J. Mulholland. Dissolved organic matter concentration and flux in streams. Journal of the North American Benthological Society, 16(1):131-141, 1997. S. P. Neuman. Adaptive eulerian-lagrangian finite element method for advection- dispersion. 20:321—37, 1984. K. Nishimoto. An Essence of Nishimoto’s Fractional Calculus (calculus in the 213t Cen- tury): Integrations and Difi'erentiations of Arbitrary Order. Descartes Press 00., 1991. J. A. Ochoa-Tapia, F. J. Valdes-Parada, and J. Alvarez-Ramirez. A fractional-order darcy’s law. Physica a-Statistical Mechanics and Its Applications, 374(1):1—14, 2007. Z. Odibat. Approximations of fractional integrals and caputo fractional derivatives. Ap- plied Mathematics and Computation, 178(2):527—533, 2006. A Ogata and R. B. Banks. A solution of the differential equation of longitudinaldispersion in porous media. Geological Survey Professional Paper,Movement in EarthMaterials., 1961. Keith B. Oldham and Jerome Spanier. The fractional calculus; theory and applications of differentiation and integration to arbitrary order. Mathematics in science and engi- neering ; v. 111. Academic Press, New York,, 1974. M. S. Phanikumar. Thermosolutal convection in a rectangular enclosure with strong side-walls and bottom heating. International Journal of Heat and Fluid Flow, 15(4): 325—336, 1994. M. S. Phanikumar and D. W. Hyndman. Interactions between sorption and biodegrada- tion: Exploring bioavailability and pulsed nutrient injection efficiency. Water Resources Research, 39(5):—, 2003. M. S. Phanikumar, I. Aslam, C. P. Shen, D. T. Long, and T. C. Voice. Separating surface storage from hyporheic retention in natural streams using wavelet decomposition of acoustic doppler current profiles. Water Resources Research, 43(5), 2007. S. Pirozzoli. Conservative hybrid compact-weno schemes for shock-turbulence interaction. Journal of Computational Physics, 178(1):81-117, 2002. 93 Igor Podlubny. Fractional differential equations : an introduction to fractional deriva- tives, fractional differential equations, to methods of their solution and some of their applications. Mathematics in science and engineering ; v. 198. Academic Press, San Diego, 1999. William H. Press. Numerical recipes in C++ : the art of scientific computing. Cambridge University Press, Cambridge ; New York, 2nd edition, 2002. Patrick J. Roache. Fundamentals of computational fluid dynamics. Hermosa Publishers, Albuquerque, N.M., 1998. A. H. Roy, M. C. Freeman, B. J. Freeman, S. J. Wenger, W. E. Ensign, and J. L. Meyer. Investigating hydrologic alteration as a mechanism of fish assemblage shifts in urban- izing streams. Journal of the North American Benthological Society, 24(3):656—678, 2005. R. L. Runkel. One-dimensional Transport with Inflow and Storage (OTIS): A Solute Transport Model for Streams and Rivers. US Dept. of the Interior, US Geological Survey; Information Services distributor, 1998. R. L. Runkel and R. E. Broshears. One-dimensional transport with inflow and storage (otis)a solute transport model for small streams: Boulder, colo., university of colorado, 1991. Muhammad Sahimi. Fractal and superdiffusive transport and hydrodynamic disper- sion in heterogeneous porous media. Transport in Porous Media, 13(1):3—40, 1993. 10.1007/BF00613269. S. G. Samko, A. A. Kilbas, and O. I. Marichev. Fractional integrals and derivatives : theory and applications. Gordon and Breach Science Publishers, Switzerland ; Philadelphia, Pa, USA, 1993. B. H. Schmid. Simplification in longitudinal transport modeling: case of instantaneous slug releases. Journal of Hydrologic Engineering, 9(4):319*324, 2004. R. Schumer, D. A. Benson, M. M. Meerschaert, and S. W. Wheatcraft. Eulerian derivation of the fractional advection-dispersion equation. Journal of Contaminant Hydrology, 48 (1-2):69—88, 2001. R. Schumer, D. A. Benson, M. M. Meerschaert, and B. Baeumer. Multiscaling fractional advection-dispersion equations and their solutions. Water Resources Research, 39(1):—, 2003. K. Sebastian and C. W. Shu. Multidomain weno finite difference method with inter- polation at subdomain interfaces. Journal of Scientific Computing, 19(1-3):405—438, 2003. S. Sema and A. Marquina. Power eno methods: a fifth-order accurate weighted power eno method. Journal of Computational Physics, 194(2):632—658, 2004. 94 C. Shen, M. S. Phanikumar, T.T. Fong, I. Aslam, S.L. Molloy, and J.B. Rose. The transport of bacteriophage p22 relative to rhodamine wt in suraface waters: The grand river, michigan. Environmental Science 63 Technology, 2007. C. W. Shu. Essentially Non-oscillatory and Weighted Essentially Non-oscillatory Schemes for Hyperbolic Conservation Laws. Institute for Computer Applications in Science and Engineering, NASA Langley Research Center; National Technical Information Service, distributor, 1997. G. Sposito and W. Jury. Fundamental problems in the stochastic convection-dispersion model of solute transport in aquifers and field soils. (1)277—88, 1986. Y. C. Tai, S. Noelle, J. M. N. T. Gray, and K. Hutter. Shock-capturing and front-tracking methods for granular avalanches. Journal of Computational Physics, 175(1):269—301, 2002. G. Taylor. Conditions under which dispersion of a solute in a stream of solvent can be used to measure molecular diffusion. Proceedings of the Royal Society of London Series a—Mathematical and Physical Sciences, 225(1163):473—477, 1954. C. W. Tsai, K. Z. Kuai, M. Bursik, and D. Hess. Nonlinear modeling of dam-break flood wave induced sediment transport and morphological evolution in rivers. AGU Fall Meeting Abstracts, 54:02—02, 2004. H. M. Valett, J. A. Morrice, C. N. Dahrn, and M. E. Campana. Parent lithology, surface groundwater exchange, and nitrate retention in headwater streams. Limnology and Oceanography, 41(2):333—345, 1996. P. J. van der Houwen and B. P. Sommeijer. Splitting methods for three-dimensional transport models with interaction terms. Journal of Scientific Computing, 12(3):215— 231, 1997. 10.1023/A:1025645326705. S. Wallis. The numerical solution of the advection-dispersion equation: A review of some basic principles. Acta Geophysica, 55(1):85—94, 2007. C. Wei and G. Morrison. Bacterial enzyme-activity and metal speciation in urban river sediments. Hydrobiologia, 235:597-603, 1992. Jm308 Times Cited:10 Cited References Count:15. M. Weilbeer. Efiicient Numerical Methods for Fractional Differential Equations and their Analytical Background. PhD thesis, Institut computational mathematics, 2006. S. W. Wheatcraft and S. W. Tyler. An explanation of scale-dependent dispersivity in heterogeneous aquifers using concepts of fractal geometry. Water Resources Research, 24(4):566—578, 1988. S. M. Wondzell and F. J. Swanson. Seasonal and storm dynamics of the hyporheic zone of a 4th-order mountain stream .2. nitrogen cycling. Journal of the North American Benthological Society, 15(1):20—34, 1996. 95 A. Worman. Analytical solution and timescale for transport of reacting solutes in rivers and streams. Water Resources Research, 34(10):2703—2716, 1998. A. Worman. Effect of flow-induced exchange in hyporheic zones on longitudinal transport of solutes in streams and rivers (vol 38, pg 1001, 2002). Water Resources Research, 38 (6) :—, 2002. Z. Yong, D. A. Benson, M. M. Meerschaert, and H. P. Schefiler. On using random walks to solve the space-fractional advection-dispersion equations. Journal of Statistical Physics, 123(1):89—110, 2006. Y. Zhang, D. A. Benson, M. M. Meerschaert, E. M. LaBolle, and H. P. Scheffler. Ran- dom walk approximation of fractional-order multiscaling anomalous diffusion. Physical Review E, 74(2), 2006. Y. Zhang, D. A. Benson, M. M. Meerschaert, and E. M. LaBolle. Space-fractional advection-dispersion equations with variable parameters: Diverse formulas, numeri- cal solutions, and application to the macrodispersion experiment site data. Water Resources Research, 43(5):—, 2007 . 96 IIIIIIIIIIIIIIIIIIIII uultllllljuljjlujln