Home

User's manual for the particle tracking model ZOOPT

1. Figure 20 Example dxf files for visualisation of particle path lines in a x y plane b x z plane and c y z plane 36 References BLANDFORD T N AND HUYAKORN P S 1991 WHPA A modular semi analytical model for the delineation of wellhead protection areas version 2 0 Report to the U S Environmental Protection Agency Washington DC FARAHMAND RAZAVI A 1995 Theory and practice of contaminant transport modelling using the random walk Unpublished Ph D Thesis University of Birmingham FRANZ T AND GUIGER N 1990 FLOWPATH two dimensional horizontal aquifer simulation model Waterloo Hydrogeologic Software Waterloo Ontario GOODE D J 1990 Particle velocity interpolation in block centred finite difference groundwater flow models Water Resources Research Vol 26 No 5 pp 925 940 JACKSON C R 2001 The development and validation the object oriented quasi three dimensional regional groundwater model ZOOMQ3D British Geological Survey Internal Report IR 01 144 JACKSON C R 2002a The representation of the variation of hydraulic conductivity with depth in the object oriented regional grou
2. Figure 19 Order of data written in heads txt 32 Table 5 Flow component Order of variables listed on line of flowbal txt Active inactive integer flag Inflow from left x direction Outflow to right x direction Inflow from below y direction Outflow above y direction Inflow from base z direction Outflow through top z direction Specified inflow to node if boundary 1 2 3 4 5 6 7 8 9 Increase in storage Recharge River leakage Leakage Pumping Spring flow 0 inactive dry 1 active wet Positive into node Positive out of node Positive into node Positive out of node Positive into node Positive out of node Positive into node Positive values are equivalent to an outflow in terms of the flow balance calculation Positive into node Positive out of node Positive out of node Positive into node Positive out of node Considering that the flow components 2 to 14 in Table 5 are represented by the terms Q Osee Q3 Q then the following equation represents the flow balance for a node Q Q Q Q Q Q7 Qs Qe Qio Qi Qn Q4 Q O 5 6 CONSIDERATIONS WHEN UNDERTAKING STEADY STATE PARTICLE TRACKING RUNS When the user specifies that the program is to be used to perform a steady state particle tracking simulation by entering s on line 2 of the input file
3. Q and Q are the flow rates entering and exiting the cell in the x direction mda D xl x2 g g y K z is the hydraulic conductivity m day in the x direction at elevation z Zg K K z dZ is the mean hydraulic conductivity m day in the i z zi Zi x direction 0 is the porosity of the node Ay is the width of the node in the y direction m and z and z are the elevations of the bottom and top of the node m Similar equations are written for the component of velocity in the y direction However because the hydraulic conductivity in the z direction is considered uniform throughout the node no modification is made to way in which the z component of velocity is calculated ZOOPT recognises when a particle enters a node in which hydraulic conductivity varies with depth and then invokes the use of the Runge Kutta technique and the application of Equations 3 4 and 3 5 Q xl DI gt Q x2 Z X X5 gt K z Figure 8 Illustration of the interpolation of velocity in VKD nodes 19 4 Running ZOOPT To install ZOOPT on a Windows PC copy the executable zoopt exe into suitable directory such as c Program FilesZOOM Then add this directory to the Windows system PATH variable Control Panel System 3 Advanced Tab Environment Variables No installation procedure is run in which ZOOPT program files are added to the system registry All the input files required by ZOOPT must be located in
4. lt Vo a Vo lt 0 strong sink Figure 3 Cell wall velocity conditions needing consideration during particle tracking In Figure 3a the interfacial velocities are the same In this case Equation 2 14 is undefined and the time of travel to exit the cell must be estimated by At x x Va if v gt 0 2 16a At x x v if v lt 0 2 16b In Figure 3b a local groundwater divide exists In this case if the particle is to the left of the divide it exits to the left otherwise it exits the cell to the right Finally Figure 3c shows the situation in which flow is toward the cell centre from both x directions In this case the particle cannot leave the cell in the positive or negative x direction If this is simultaneously true of the flows in the y and z directions then the node is termed a strong sink and the particle is terminated at the cell 2 3 0 Numerical integration techniques Equation 2 2 can also be solved numerically In fact if a higher order interpolation scheme is used then it can only be solved using numerical methods Numerical integration methods involve the movement of the particles in discrete tracking steps along the path line Of the numerical integration techniques Euler s method is the simplest In this method the velocity at the current point is extrapolated over the tracking step The particle is moved along the path line using the equations x 2 x t V At 2 172 y 2 y VA 2 17b Z Z
5. cod OR vkdky dat per vkd scheme 2 vkdzp map amp vkdzp cod OR vkdzp dat per vkd scheme 22 vkdgrad map amp vkdgrad cod OR vkdgrad dat per vkd scheme 23 zbase map amp zbase cod OR zbase dat per layer 24 zoomg3d dat 25 ztop map amp ztop cod OR ztop dat per layer Note the names of these files are fixed The names of the remaining files can be specified by the user in the input file zoomq3d dat Table2 List of the ZOOMQ3D output files required by ZOOPT 1 heads txt 2 flowbal txt Note the names of these files are fixed Table3 List of the additional input files required by ZOOPT zoopt dat particles dat porosity map amp porosity cod porosity dat per layer porosity a map amp porosity a cod OR _ porosity a dat per layer Note the names of these files are fixed 2 Particle tracking theory Particle tracking is commonly used to define the path lines of solute particles under purely advective transport The technique is often applied for the definition of borehole catchments and associated source protection zones the identification of recharge and discharge areas and the visualisation of groundwater flow patterns However the method also forms the basis of a number of solute transport models which simulate the effects of hydrodynamic dispersion Random walk methods Prickett et al 1981 Farahmand Razavi 1995 and the method of characteristics Konikow and Bredehoeft 1
6. 1 j m jl e e e e Finite difference node Cell wall J e e vn y A j1 9 e v gt X i 1 i i Figure 1 Inter nodal velocity components in the x direction Given the inter nodal velocities obtained from Equation 2 3 the x component of the velocity at any arbitrary location within the cell i j can be calculated using linear interpolation between the two opposite cell walls The component of the velocity in the x direction is calculated at an arbitrary x co ordinate between x and x using the following linear interpolation equation V x 2 A s i V 2 4 2 where v run 2 5 i i M mS Similar equations to Equations 2 3 2 5 are defined in the y and z directions by replacing the terms x and 1 by the terms y and j or z and k respectively These are used to calculate the remaining two directional components of the velocity vector at any arbitrary position in the model domain 2 5 PATHLINE DEFINITION 2 3 1 Semi analytical technique If linear interpolation is used to calculate the x y and z directional components of the velocity field at a particular position in the finite difference model domain as described above then the integral in Equation 2 2 can be solved analytically Pollock 1988 called this method of particle tracking the semi analytical technique because of the combination of a numerical velocity interpolation routine and an analytical path line definition procedure Considering the x compone
7. Mansour of the British Geological Survey and P J Hulme of the Environment Agency for their help in reviewing the ZOOPT software Additionally the author would like to acknowledge the assistance of A T Williams of the British Geological Survey for reviewing this document Preface to the second edition The production of the second edition of the ZOOPT manual coincides with the release of version 1 03 of the code This version of the code incorporates one only change to version 1 02 DIFFERENCES BETWEEN VERSION 1 03 AND 1 02 OF ZOOPT e All executables should now be placed in a suitable directory e g c Program Files ZOOM and this folder should be added to the Windows system PATH variable ZOOPT can then be run from any working directory by typing the name of the executable followed by the path of the working directory e g ZOOPT c myDirectory Alternatively this string could be placed in a batch file and the batch file run from the command line ii Contents i EorewWord ieeiiiossiseccsrcresceocesosscstnesteosaosoos esteos cooo rs sse ioo coas ro esise poose ae insere sieer p Acknowledgements eecscevexenn ve essa sev ae a aka ev Von ta rea veo paren vx E Eee Ya char neu rinse e ei Pe ci er oce ve reverso E MIHI IS AA mR oa HARP BG AR AONE TSE EYEE t Een Nu OH E 1 i Introduction iudei Ion ea debo e Coe doner c e Yee eo odor eva AA AAA MENA WAA HANA eode gosse da ANA KUA KAN KE AE lb Perm npology ene etae ep
8. The quasi layer Ola is below layer 01 of finite difference nodes Consequently there are never porosity input files for a quasi layer beneath the bottom finite difference layer The file numbering scheme is illustrated in Figure 18 30 Finite difference node Layer 01 e e e e e Layer 01 e e e e e Quasi layer 01a Layer 02 Layer 02 e e e e e e e e 9 e Zh a b X Figure 18 Layer numbering scheme used for specification of porosity in a a model without quasi layers and b a model with quasi layers 5 5 NODE BY NODE HEADS AND FLOWS Two space delimited ASCII text output files are produced by ZOOMQ3D if the user specifies that particle tracking simulations are to be performed after it has been run These are named heads txt and flowbal txt The first contains simulated groundwater head values for each model node for each time step of the simulation The second contains the flow rates m day into and out of each node for each time step during the simulation These ZOOMQ3D output files are required by ZOOPT as input The formats of each of the files are described next 5 5 1 Output file heads txt Heads are contained in this file for each model node at the end of each time step of the ZOOMQ3D simulation The heads are written in the following order during the flow simulation 1 For each time step of the simulation 2 Foreach layer in the model Starting with th
9. be specified Therefore for example if twenty six letters were used in the map file for a specific grid i e the letters a to z twenty six parameters values must defined within the code file for the corresponding grid In the above description it is stated that there must be an equivalent number of parameter values specified in the code file to the number of letters used in the corresponding character array in the map file However this is not strictly the case In fact there must be at least the same numbers of parameters specified in the code file as are used in the map file for the corresponding grid Therefore it is allowable for example to specify parameter values for all fifty two characters in the range a b y z A B Y Z but only use a subset of these letters in the map file The final rule to follow when editing the map and code files is that the minimum number of parameters to enter in the code file is determined by the highest letter in the range a b y z A B Y Z For example if C is the highest letter used then twenty nine parameters values must be defined in the code file Similarly if m is the highest letter used then thirteen parameters values are required This means that if only the two letters a and z are used to define a grid character array in a map file at least twenty six parameters values must still be defined in the corresponding code file for the corresponding
10. files except each of the character arrays is replaced by a set of data defining the values of the letters contained within the array The code file which corresponds to this example model and map file is shown in Figure 16 Similarly to the map file the code file contains four sets of data separated by comment lines The comment lines contain the information required to identify to which grid each data set relates Again ZETUP the set up program for ZOOMQ3D produces this file and writes the co ordinates of the bottom left and top right hand corners of the mesh within the comment lines Consequently the user need only adjust the parameter values for each grid In Figure 15 eight letters are used in the base grid so eight parameter values must be defined in the code file for this grid On the second grid listed in the map file only two zones are defined by the letters a and b and consequently only two parameter values must be assigned in the code file for this grid This relationship is the same for the final two grids in 27 the map file Re iterating because eight letters are defined on the base grid in the map file eight corresponding parameter values must be defined in the code file for this grid The first integer number in each block of data between the comment lines in the code file is the number of zones or letters specified in the character array After this integer value an equivalent number of parameter values must
11. input file and file format 24 i or time variant t Comment line Two character flags Comment line Integer Comment line Two decimal numbers Comment line Charac Commen Decimal Comment 1 Character Comment line Decimal Comment line Three integers 5 3 FILE FORMATS FOR INPUTTING SPATIAL DATA POROSITY Spatial information can be entered into the model using two different methods or sets of input files These methods employ either 1 map and code input files or 2 numeric input data files In method 1 a pair of ASCII text files is required for each model layer with the file extensions map and cod A map file contains a single array or multiple arrays if grid refinement is incorporated of characters in the range a b y Zz A B Y Z to define zones of different parameter values That is the range is composed of the lower case alphabet followed by the upper case alphabet Each array represents one of the grids contained in the model These characters define zones across the model mesh at which a particular value of a parameter is specified The values of the parameters are specified in the corresponding file with the cod extension A maximum of fifty two zones can be used to define the spatial distribution of model parameters when using this method The name of each of these pairs of files is suffixed with a two digit number that represents the model layer that the data applies to For example t
12. mesh composed of four grids and b representation of the grid hierarchy 26 The structure of the map file for this example model is shown in Figure 15 As described above this contains four character arrays separated by comment lines The comment lines contain information required to identify to which grid each array relates ZETUP the set up program for ZOOMQ3D produces this file and writes the co ordinates of the bottom left and top right hand corners of the mesh within the comment lines Consequently the user need only adjust the letters in the arrays in order to modify the model s parameter values Map for grid Bottom left rights 700 700 aaaaaaaa bbbbbbbb ccccccoc dddddddd eeeeeeee ffffffff Map for grid Bottom left 100 100 Top right 300 300 bbbbb bbbbb bbbbb bbbbbbbb bbbbbbbbb aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa Map for grid Bottom left 400 400 Top right 600 600 abcde abcde abcde abcde abcde Map for grid Bottom left 150 150 Top right 250 250 aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa aaaaaaaaa Figure 15 Example map file for the entry of spatial data within a layer Each letter in the array represents a value of a parameter that is specified in the corresponding code file i e the file with the same name but with the cod extension Code files have the same structure as map
13. procedure is illustrated in Figure 6c A full derivation of this correction term is presented by Zheng 1994 3 3 2 Unconfined aquifer layers In phreatic aquifer layers the elevation of the water table determines the vertical thickness of the finite difference node The elevation of the water table is defined as the simulated groundwater head Consequently in unconfined model layers the nodes are distorted vertically again The position of particles that are tracked through unconfined nodes are corrected in the same way as described for fixed but vertically distorted grid nodes 3 4 QUASI THREE DIMENSIONAL LAYERING Not all hydrogeological layers are always included explicitly in a groundwater model Consider that a groundwater system is composed of a sequence of high permeability horizontal layers that are separated by a series of low permeability aquitards Often these aquitards are not represented by a series of finite difference nodes in a groundwater model This is because the assumption can be made that the flow in the low permeability layers will be essentially vertical Consequently low permeability layers are commonly modelled by adjusting the vertical conductance between the two adjacent aquifer layers When this approach is adopted the model is stated to be quasi three dimensional Whilst this method of representing aquitards in groundwater models is computationally efficient the particle tracking routine has to be modified because t
14. the particle on release local or global depending on the layer number Layer is the number of the layer in which the particle is released and t is the time of release of the particle days The z co ordinate can be specified as either a local co ordinate within a layer or as a global co ordinate If the layer number specified in the file is 0 the z co ordinate is global and is related to the elevation of the datum of the model as defined by the user If the layer integer number is greater than zero the z co ordinate is local and must be between 0 0 and 1 0 i e it is defined as a fraction of the layer thickness If it is 0 0 then the particle is placed on the bottom of the layer If it is 1 0 then it is placed at the top of the layer If the program is tracking particles backwards the time of release t is automatically re set to Zero days 23 zoopt dat File format Forward f or back tracking b steady state s time instant bs Number of intermediate tracking points within cell gt 0 10 Runge Kutta safety factor lt 1 0 and error criterion m 0 9 0 000001 Enforce the use of Runge Kutta y or n n 1 2 3 4 5 6 7 8 9 DXF drawing z scale factor gt 1 0 10 0 Draw grid in DXF file y or n y Length of time for which to track particles days 100000 0 Time step stress period amp block for time instant tracking 3 12 4 Figure 12 Example zoopt dat
15. 05 0 05 0 05 0 05 0 05 0 05 U5 0 05 0 05 0 U 5 0 4Uo 0 05 0 05 05 0 05 0 05 0 05 0 05 0 05 0 05 s05 0 05 0 05 0 05 0 405 0505 0 05 05 0 05 0 05 0 05 0 05 0 05 0 05 0 0 0 0 15 15 15 15 15 05 0 05 0 05 0 05 05 0 05 0 05 05 0 05 0 05 0 05 05 0 05 0 05 05 0 05 0 05 0 05 05 0 05 0 05 05 0 05 0 05 0 05 05 0 05 0 05 OOOOOoo0000 zo00000zo0o000o000o00zxoooooooo Qo 0 0tf10 0o QUO C O0 GO OO 0 O 1 0 OOo OO OO 0 Oooo QUO uo D oO Q0 0 Oo Figure 17 Example numeric data file for the entry of spatial data within a layer 29 5 3 4 Definition of data at coincident points on multiple grids It is apparent from the ezample model and input files described above that parameters values at particular finite difference nodes are specified more than once in an input file when the model mesh is locally refined This occurs at those nodes on a refined grid that are coincident with the nodes on a coarser grid Where this occurs the data defined on the refined grid will always override that defined on the coarse mesh 5 4 SPATIAL DEFINITION OF POROSITY Layers of finite difference nodes in ZOOMQ3D and ZOOPT may abut or may be separated by a quasi layer This quasi layer is not represented by a layer of finite difference nodes Rather flow in this layer which for example could represent a low permeability aquitard between two aquifers is assumed to be vertical and depends on the vertical conducta
16. 978 Zheng 1990 use particle tracking to describe the advective component of solute transport 2 1 MATHEMATICAL BACKGROUND Assuming that fluid density is uniform the path lines of contaminants under advection alone are governed by the equation dp vlp t 2 1 q Yt 2 1 where p xi yj zk is the position vector and Vxi vyj vzk isthe seepage velocity vector Koh V 0 x K is the hydraulic conductivity in the x direction 0 is the porosity and h x y Z t is the groundwater head which is a function of space and time The solution of Equation 2 1 for the position of a particle at time t is p t plte fv p t at 2 2 where p t is the initial position of the particle at time t The solution of Eguation 2 2 reguires the evaluation of the velocity field at any given time and position in the model domain If an exact solution for the velocity field exists then Equation 2 2 can be solved analytically for p t However this is generally not the case and Equation 2 2 must then be solved numerically In a finite difference model velocity components are only known at specific locations in the aquifer that is at the position of the cell wall between two adjacent nodes Consequently an interpolation scheme must be used to evaluate the velocity field at arbitrary positions and times This means that an analytical solution to Equation 2 2 cannot be calculated Furthermore the selection of the interpolat
17. AA AA AA AA AAA NON eU ENTER NU KA AAA d GJ Output file tracks OUU os osito CDM EOS DU i pae i DR E tu duc ic tdt tiun 34 6 2 Output file piend Out ases ca d ti ua 34 6 3 Visualisation of particle paths oce eicere ts enne Ue genes ise eda pu tese pei Qe sas ie edu 35 6 4 Error file AP ii wawi 35 References ossec oossoo n dosse asese Seose eose EEDE E A eo e SDEsse eE SAT Oi is i a OL iii FIGURES Figure 1 Inter nodal velocity components in the x direction ssesmmmnmemamanmnemmamawa 5 Figure 2 Particle track through a two dimensional cell eene 7 Figure3 Cell wall velocity conditions needing consideration during particle tracking 8 Figure 4 Intermediate steps taken during fourth order Runge Kutta method 10 Figure5 Illustration of a well as a distributed sink at the centre of a finite difference node showing surrounding cell walls in a three dimensions and b plan view 14 Figure 6 Correction of particle elevation in vertical distorted model layers 16 Figure 7 Parameters used to define VKD profiles in ZOOMQ3D ee 18 Figure8 Illustration of the interpolation of velocity in VKD nodes sess 19 Figure9 Starting a command line window from the Windows start menu 20 Figure 10 Example of changing the working director
18. Centre Harrier Way Sowton Exeter Devon EX2 7HU T 01392 445271 Fax 01392 445371 Geological Survey of Northern Ireland 20 College Gardens Belfast BT9 6BS 4 028 9066 6595 Fax 028 9066 2835 Maclean Building Crowmarsh Gifford Wallingford Oxfordshire OX10 8BB 01491 838800 Fax 01491 692345 Sophia House 28 Cathedral Road Cardiff CF11 9LJ T 029 2066 0147 Fax 029 2066 0159 Parent Body Natural Environment Research Council Polaris House North Star Avenue Swindon Wiltshire SN2 1EU 01793 411500 Fax 01793 411501 www nerc ac uk Foreword The development of the modelling software within the ZOOM family of which ZOOPT is a part has been undertaken through a continuing tripartite collaboration between the University of Birmingham the Environment Agency and the British Geological Survey The development of the flow model ZOOMQ3D was initially undertaken at the University of Birmingham between 1998 and 2001 but continued after this time as a collaborative project between the three partner organisations Since the inception of the collaborative project the development of the software has been directed by the ZOOM steering committee the members of which are University of Birmingham Dr Andrew Spink Environment Agency Steve Fletcher Paul Hulme British Geological Survey Dr Denis Peach Dr Andrew Hughes Dr Chris Jackson Acknowledgements The author would like to acknowledge the assistance of A G Hughes and M M
19. a single directory All the output files produced by ZOOPT will be created in the same directory ZOOPT should be run from the command line in a MS DOS console box and not started from Windows Explorer To start a command window select Run from the Windows start menu and type cmd in the drop down list box Figure 9 The user should then change directory to that of the working directory where the input files are located Figure 10 For help on the commands used to change directory type help cd within the console box Figure 10 To run the model type zoopt followed by the path to the working directory on the command line e g zoopt c myDirectory Alternatively this string can be placed in a batch file a text file with a bat extension e g runzoopt bat and the name of this batch file can be typed on the command line omit the extension when doing this e g type runzoopt AZ 15 Type the name of a program Folder document or Internet resource and Windows will open it For you Open cmd Cancel Browse Figure9 Starting a command line window from the Windows start menu In the event that an error occurs messages are written to the screen If ZOOPT is run from Explorer it may terminate before the user is able to read the error messages The program reads data from ASCII text input files which must be located in a single working directory The formats of all the input data files
20. ameters used to define VKD profiles in ZOOMQ3D The value of the VKDGrad parameter may be either negative zero or positive Consequently in addition to an increase in hydraulic conductivity with depth above Z hydraulic conductivity can be specified to decrease or remain constant To calculate transmissivity the following equations are used T Ki h Zorrom 0 5 VKDGrad h Z T K h Zyorrom 0 5VKDGrad h Z Y for h gt Z and T K h Zsorrom T Ky h Zgorrom for h lt Zp where h is the water table elevation At those nodes of the finite difference grid where hydraulic conductivity varies with elevation the Runge Kutta particle tracking technique must be used to define path lines The integral in Equation 2 2 cannot be evaluated analytically because the horizontal velocity varies in the z direction that is towards the top of the VKD profile the horizontal velocity is greater than towards its base Consequently the semi analytical method cannot be used Within ZOOPT the assumption is made that the horizontal velocity of the particle is proportional to the horizontal hydraulic conductivity at its location With reference to Figure 8 the component of the velocity in the x direction at the cell walls at x and x2 are given by Valz Ke B 3 4 Ay z z Ke Qs V 2 z K N 2 3 5 where V z and V z are the x components of velocity m day on the cell walls at elevation z
21. and z are the x y and z co ordinates of the cell walls a K i K is the ratio of hydraulic conductivities in the x and y directions and 0 is the porosity of the node This interpolation scheme forces particles either to converge towards the well or to leave the finite difference node through one of its interfaces Particles are terminated if they enter within the radius of the well The scheme is implemented in ZOOPT and is an elegant solution to the problem of weak sink wells At wells that generate strong sinks this interpolation scheme is not applied and all particles are terminated as they enter the node Vz Q Vv WELL Va Aa e i uides xp Yp Xp yp Zp A Vx gt Vy Vaa gt Vy x ji v A Pi Vy a Vn b Vz Figure 5 Illustration of a well as a distributed sink at the centre of a finite difference node showing surrounding cell walls in a three dimensions and b plan view 33 GRID CONSIDERATIONS 3 3 4 Distorted vertical discretisation Because of efficiency considerations finite difference models are often constructed in which the elevation of the top and bottom of model layers varies over the model domain This variation generally results from the need to approximate the changing shape of hydrogeological units However vertical distortion of the mesh causes problems because the particle tracking solution is based on a fixed orthogonal grid For example consider that it is calculated that a particle l
22. ath analysis Bulletin 69 Illinois State Water Survey Champaign IL ZHENG C 1989 PATH3D a groundwater path and travel time simulator version 3 0 user s manual S S Papadopolous amp Associates Inc Bethesda MD ZHENG C 1990 MT3D A modular three dimensional transport model for simulation of advection dispersion and chemical reactions of contaminants in groundwater systems Report to the U S E P A Ada OK ZHENG C 1994 Analysis of particle tracking errors associated with spatial discretisation Ground Water Vol 32 No 5 pp 821 828 ZHENG C AND BENNETT G D 1995 Applied contaminant transport modelling Van Nostrand Reinhold New York 37
23. ation Shafer 1987 Bi cubic Runge Kutta Pollock 1989 Linear Semi analytical Zheng 1989 Linear Fourth order Runge Kutta Franz and Guiger 1990 Reverse distance Euler integration Blandford amp Huyakorn 1991 Linear Semi analytical or Euler Linear velocity interpolation is implemented in ZOOPT As stated above this enables the analytical solution of the integral in Equation 2 2 The approach also maintains consistency with the finite difference mass balance equations is generally more accurate than other methods in heterogeneous media and is computationally efficient These issues are discussed by Zheng and Bennett 1995 who state that simple linear interpolations schemes are generally preferable to multi linear interpolation schemes Furthermore the approach can easily form the basis of a semi analytical time variant particle tracking technique presented by Lu 1994 The calculation of the velocity at the cell wall between two nodes is based on the inter nodal volumetric flow rate calculated by the flow model ZOOMQ3D Considering the component in the x direction as illustrated in Figure 1 the velocity at the cell wall at position i j denoted by V is 2 Q 1 i V 2 m da za 2 3 si T gayag nd 2 3 where Q is the inter nodal flow rate calculated by the flow model m day i gt j 2 J O is the aquifer porosity of cell i j Ay y wan Ya y 2 m and Az is the aquifer thickness of cell
24. cles path lines are plotted as dxf files The plots can be stretched in the z direction by entering a number greater than 1 0 This number is used to multiply the z values of all the particle locations when writing the dxf file Line 12 A single character flag is entered on this line This determines if the model grid is written to the dxf files for visualisation in addition to the particle path lines Line 14 On this line the length of time is specified days for which to track the particles The distance that particles move during this period is dependent on the porosity of the model layers Line 16 On this the last line of the input file three integers are defined These are the time step stress period and block used to define time instant particle tracks see description of line 2 The block number is a counter and not the actual value of the block i e if the blocks represented years which had values starting from 1970 a value of five would be entered to define the fifth year not the number 1974 52 INPUT FILE PARTICLES DAT The number of particles to be tracked and their starting positions are defined in the input file particles dat This file has the following format n x y z layer t one line of data for each of n particles where n is the number particles to be tracked x is the x co ordinate m of the particle on release y is the y co ordinate m of the particle on release z is the z co ordinate m of
25. cular model features for example weak sink nodes or nodes which exhibit a vertical variation of hydraulic conductivity particles have to be tracked using the Runge Kutta method The switch between the use of the semi analytical and Runge Kutta methods is made automatically within the model code Though the Runge Kutta technique is only implemented occasionally the user can enforce its continuous use This option provides an alternative to the semi analytical technique though the semi analytical method is computationally more efficient and should be used in preference to Runge Kutta tracking ZOOPT is designed to track particles in models where the horizontal hydraulic conductivity varies with depth VKD and this is one feature that requires the application of the Runge Kutta method In this case the horizontal velocity depends on the hydraulic conductivity at the elevation of the particle within the node Consequently the integral in Equation 2 2 cannot be determined analytically and a fully numerical tracking method must be employed In addition to the application of the particle tracking code to VKD nodes the code is compatible with the local grid refinement technique incorporated within ZOOMQ3D Local grid refinement enables the zooming of the mesh within discrete areas of a model grid to increase accuracy or model detail ZOOPT tracks particles through these locally refined grids The occurrence of weak sinks which is a problem assoc
26. direction last Within each row nodes are scanned from left to right in the positive x direction Data for a node within a subgrid is only output to the file if it does not also exist on a grid at a coarser grid level i e if it is not located on a parent grid line Because the data is not structured in grid arrays dummy flow values relating to nodes outside the model boundary are not written to the file Table 5 shows the flow components that are listed within each line of nodal data The first parameter output is an integer flag which specifies if the node was active during the time step i e if it was wet 1 or dry 0 If the node was dry a zero is written on the line for the node but its flow components are not Flow data is only output if the node was active wet during the time step Base grid array Layer 1 Subgrid 1 array Subgrid 2 array Base grid array Time step 1 Layer 2 Subgrid 1 array Subgrid 2 array Base grid array Layer 3 Subgrid 1 array Subgrid 2 array Base grid array Layer 1 Subgrid 1 array Subgrid 2 array Base grid array Time step 1 Layer 2 Subgrid 1 array Subgrid 2 array Base grid array Layer 3 Subgrid 1 array Subgrid 2 array Base grid array Layer 1 Subgrid 1 array Subgrid 2 array Base grid array Time step 3 Layer 2 Subgrid 1 array Subgrid 2 array Base grid array Layer 3 Subgrid 1 array Subgrid 2 array NB Example file structure for three time steps only for a model containing three layers and three grids
27. e etai Ret actress 1 1 2 UNI CONVENTION 3x1 eee des oce e Certe eot reve e des voee eoe ree te tia eve C bared Seats 1 Particle tracking he Ory ssssssssssissssssssssesossissssssssssesvesssisssssss osossssesosassssesoscs ssssso osos ssrsssosssssst J 2 1 Mathematical background 5 rt ikala UI tei iiia 3 2 2 Meloems mterpolatioli 252 5 tote dtes hes eode ti uid t ii eat d aaa o eas tuca fep abeant 4 2 9 Patbline definition uie rr WS RET au ENS RUNE GER T ux 6 Capabilities of the particle tracking code ZOOPT eeeeee eene eee ee eese eeeeess L2 ex Michi n ber waid dn Satats a aamaney gurl aes waenk dete E EENE 12 3 2 Velocity calculation senna n css wade sia ass Oates 13 3 9 Ord eonsideratlons ss Aa ka ia 14 3 4 Quasi three dimensional layering eese esent mwema 15 3 5 Vertical variation of hydraulic conductivity with depth 17 ZOOPT input tiles ia Set InpubHle ZoopLUdat esse r ra ties fue abe a as o Eni 22 5 2 Input file particles dat aa 23 5 3 File formats for inputting spatial data porosity essere 25 54 Spatial definition of porosity ec ee Tete tie S eI dac a USERS e Ye RENT ae End do Re seas 30 2 3 Node by node heads and flows sicisiccicsicsssscastecscesegicdsancndenseedesasdues ene aoo e enar etse e Ee noon 31 5 6 Considerations when undertaking steady state particle tracking runs 33 ZOOPE ontput files disco
28. e mesh within the model layer that the file relates to An example numeric data file is shown in Figure 17 This specifies exactly the same model data as that specified in the example above in which a pair of map and code files are used for defining a hydraulic conductivity distribution Map for grid on level 1 Bottom left 0 0 Top right 700 700 10 0 10 0 10 0 10 0 10 0 10 0 10 212505 012 0X12 0 12 0212 0 12 012 14 0 14 0 14 14 0 14 14 14 20 0 20 0 20 0 20 0 20 0 20 0 20 01 0 01 0 01 0 01 0 01 0 01 0 01 02 0 02 0 02 0 02 0 02 0 02 0 02 03 0 03 0 03 0 03 0 03 0 03 0 03 07 0 07 0 07 0 07 0 07 0 07 ap for grid on level 2 Bott 20 0 20 0 20 0 20 0 20 0 20 0 0 0 0 0 0 0 0 1 gt O O GOG aO OOO O OGO ao e m le 20 20 20 20 20 10 T0 t 200 100 Top tight 300 300 7 20 0 20 0 20 0 20 0 20 0 20 20 0 20 0 20 0 20 0 20 20 0 20 0 20 0 20 0 20 20 0 20 0 20 0 20 0 20 10 0 10 0 10 0 10 0 10 10 0 10 0 10 0 10 0 0 20 0 0 0 0 0 10 0 10 0 10 0 10 0 o 0 0 0 0 20 20 10 T0 10 0 10 10 04 10 0 10 2 Bottom le OOooOooocooococoocooococoocoocooco 0 0 0 0 0 0 0 0 0 CO 000 000 Oo LO LO 10 0 10 0 10 0 10 ap for grid on level e Hj 0 00 0 Oo 6S Fh t 400 400 Top right 600 600 209 05208 0 07 09 0 08 0 07 09 0 08 0 07 09 0 08 0 07 0 09 0 08 0 07 ap for grid on level Bottom left 150 150 Top right 250 250 05 0
29. e of the particle y co ordinate of the particle z co ordinate of the particle Total travel time days DI on oP CES YA Total distance moved m 6 2 OUTPUT FILE PTEND OUT This output file contains information on the starting and finishing positions of each particle and the total distance and time of travel The file contains nine space delimited columns that contain the following data listed below Each line of the file relates to a single particle 1 Particle number x co ordinate of the end position of the particle m y co ordinate of the end position of the particle m z co ordinate of the end position of the particle m Total travel time days Total distance moved m x co ordinate of the start position of the particle m y co ordinate of the start position of the particle m NA 09 m SOY A a QS z co ordinate of the start position of the particle m 34 6 3 VISUALISATION OF PARTICLE PATHS ZOOPT produces dxf files for visualisation of the particle path lines This is a standard file format used by CAD software which can be viewed using for example AutoCAD ArcView ArcMap or Surfer If these commercial products are not available to the user dxf viewers can be downloaded from the internet Three dxf files are produced by ZOOPT which are the same except for the orientation of the Cartesian axes tracks_xy dxf tracks_xz dxf and tracks yz dxf Only one of these files is required if the software pac
30. e top layer layer 01 3 For each model grid in the layer Listed in the same order as they are in the grids dat input file The base grid is the first grid listed Grids are then listed in grid level order Heads are output as an array representing each grid This format of this file is illustrated in Figure 19 Some points on a grid may be located outside the model boundary however a groundwater head value is written for all elements of the rectangular array or matrix of grid points For the nodes outside the model boundary a dummy head value of 999 0 is output This head value is also output at those nodes that have de watered during the time step 5 5 2 Output file flowbal txt All the components of flow into and out of a node are contained in this file for each model node at the end of each time step of the simulation The flows are written in the following order during the ZOOMQ3D simulation 31 1 For each time step of the simulation 2 For each layer in the model Starting with the top layer layer 01 3 For each model grid in the layer Listed in the same order as they are in grids dat The base grid is the first grid listed Grids are then listed in grid level order 4 Data for each node within the model are written on separate lines The grids are scanned seguentially by row Data are written for the nodes in the top row furthest in the positive y direction first and bottom row furthest in the negative y
31. eaves a model cell in the horizontal and positive x direction If the particle leaves towards the top of the cell and the grid is distorted then it is possible for the particle to enter a node contained 14 in an upper layer This situation is shown in Figure 6a If this occurs the particle position must be corrected before the next particle move is made When the semi analytical solution is adopted the particle s elevation is adjusted on the interface between nodes that is when its passes from one node to the next At this point the local z co ordinate of the particle prior to the correction with respect to the top and bottom elevations of the first node must be equal to its local z co ordinate in the node it is entering after the correction This calculation is shown in Figure 6b Because of this correction the particle path can appear unsmooth when plotted This problem is inherent in the representation of three dimensional models as vertically varying layers When the Runge Kutta method is used the particles position must be modified after the tracking step within the node the particle has entered The required vertical correction is calculated by considering that if the vertical velocity component was zero the particle s local z co ordinate within each cell would have remained the same The correction factor Az is calculated using this assumption and is given by _ Az Az AZ z z z z 3 3 The correction
32. ed on a step wise process a set of objects are defined that exchange messages The user of ZOOPT can think of an object in abstract terms as any distinct entity that stores data and performs tasks In ZOOMQ3D and ZOOPT objects are defined to represent real world features For example a pumped well is represented by an object Pumped wells are described by data such as a depth and radius and have the capability to pump water out of an aquifer References are made in this manual to particles and these are also represented by objects Each particle is characterised by a position and can move by advection through the aquifer 1 2 UNIT CONVENTION All lengths in ZOOPT are specified in metres The unit of time is specified as days Table1 List of the ZOOMO3D input files required by ZOOPT K anisotropy map amp anisotropy cod OR anisotropy Ht dat per layer K aquifer map K boundary dat K clock dat K fixedheads dat grids dat K K hydcond map amp hydcond cod OR hydcond dat per layer leakage dat K 1 2 3 4 5 entry_method dat 6 7 8 9 K 10 noflow map per layer 1 pumping dat 12 recharge dat 13 recharge cod amp recharge map amp OR recharge rates dat K 14 rivers dat K 15 springs dat 16 vcond map amp vcond cod OR vcond dat per layer 17 vkd cod amp vkd map 18 vkd dat 19 vkdkx map amp vkdkx cod OR vkdkx dat per vkd scheme 20 vkdky map amp vkdky
33. el features the Runge Kutta technique is implemented in order to solve some specific problems associated with particle tracking The problem of particle termination at weak sink nodes is solved by the application of the special velocity interpolation scheme presented by Zheng 1994 This approach enables the definition of borehole catchments around wells that induce weak sinks which is not possible with many other widely used particle tracking codes ZOOMQ3D incorporates the representation of the vertical variation of hydraulic conductivity with depth VKD within finite difference nodes This has been implemented in the flow model to enable the more accurate description of the variation of hydraulic conductivity in limestone and particularly Chalk aquifers ZOOPT is fully compatible with VKD models ZOOMQ3D also enables the local refinement of the finite difference grid for example around pumping wells Again ZOOPT is fully compatible with this model feature and can be used to track particles through such refined meshes ZOOPT has been rigorously tested through its comparison with an analytical solution and another particle tracking code and through the inspection of path lines generated using numerous test models Jackson 2002b 1 Introduction ZOOPT is the particle tracking code associated with the groundwater flow model ZOOMO3D The program can track the advective movement of particles through steady state flow fields in both the f
34. grid This is because the parameter values defined in code files must be listed in the order a to z then A to Z i e in alphabetical order with the lower case alphabet preceding the upper case alphabet in the fifty two character scheme Codes for grid Bottom right 700 700 8 0 10 0 12 0 14 0 20 0 01 0 02 0 03 0 07 Codes for grid Bottom 100 100 Top right 300 300 2 0 10 0 20 Codes for grid E Bottom 400 400 Top right 600 600 5 0 13 0 15 0 09 0 08 0 07 Codes for grid Bottom 150 150 Top rights 2507 250 1 0 05 Figure 16 Example code file for the entry of spatial data within a layer 28 5 3 3 Specifying spatial data using numeric data files Instead of using a pair of map and code files to specify the spatial variation of a parameter within a model layer a single data file can be used containing raw numeric data This has the same name as the map and code files but has the dat extension For example the file porosity04 dat could be used to specify porosity values in layer 04 of the model i e three layers down from the top layer instead of the files porosity04 map and porosity04 cod The numeric data file has the same structure as the map file except that each character array is replaced by an array of numbers representing the parameter values that will be applied at the corresponding nodes of the finite differenc
35. haa Survey NATURAL ENVIRONMENT RESEARCH COUNCIL User s manual for the particle tracking model ZOOPT Groundwater Systems amp Water Quality Programme Internal Report IR 04 141 BRITISH GEOLOGICAL SURVEY GROUNDWATER SYSTEMS amp WATER QUALITY PROGRAMME INTERNAL REPORT IR 04 141 User s manual for the particle tracking model ZOOPT C R Jackson The National Grid and other Ordnance Survey data are used with the permission of the Controller of Her Majesty s Stationery Office Ordnance Survey licence number Licence No 100017897 2004 Keywords Particle tracking ZOOPT ZOOMQ3D Bibliographical reference JACKSON C R 2004 User s manual for the particle tracking model ZOOPT British Geological Survey Internal Report IR 04 141 46pp Copyright in materials derived from the British Geological Survey s work is owned by the Natural Environment Research Council NERC and or the authority that commissioned the work You may not copy or adapt this publication without first obtaining permission Contact the BGS Intellectual Property Rights Section British Geological Survey Keyworth e mail ipr bgs ac uk You may quote extracts of a reasonable length without prior permission provided a full acknowledgement is given of the source of the extract NERC 2004 All rights reserved Keyworth Nottingham British Geological Survey 2004 BRITISH GEOLOGICAL SURVEY The full range of Survey publications is ava
36. he pair of files porosityOl map and porosity01 cod would be used to specify the spatial distribution of porosity in layer 01 the upper model layer In method 2 a single ASCII text file is required for each model layer with the dat file extension Each of these files contains a single array or multiple arrays of numeric data with each array representing one of the grids contained in the model These numbers specify the value of a specific parameter at each finite difference node within the model directly The name of each of these files is suffixed with a two digit number that represents the model layer that the data applies to For example the file porosity03 dat would be used to specify the spatial distribution of porosity in layer 03 i e two layers below the top layer of finite difference nodes 5 3 1 Selection of spatial data entry method Either of the methods described above can be used to assign values for each of the parameters that can vary spatially across the model domain To specify which data entry method is to be used the user must adjust the input file entry method dat The format of this file is shown in Figure 13 It consists of a series of integers each of which is preceded by a comment line There is also an additional comment line at the top of the file Each comment line except for the first notifies the user of the related parameter e g specific yield The entry mode of the specified parameter is entered on
37. here are no grid nodes associated with these low permeability layers ZOOPT recognises the presence of quasi three dimensional layers and moves the particles vertically through them The time of travel through the low permeability layer is calculated from the leakage between the two adjacent simulated aquifers however the user must specify the porosity of each aquitard 15 a Uncorrected particle path b Semi analytical solution c Runge Kutta solution Correction at nodal interface Correction after tracking step A A i Correction i i Pone A A i A A A Z 1 v A Z 1 z fA AA Z a AZ AZ Y v 52 G zi v v i v Z Z Z n l YA 6z Az Z aH Z at sp z z zi Z5 3 Za Az AZ Az Figure 6 Correction of particle elevation in vertical distorted model layers 16 3 5 VERTICAL VARIATION OF HYDRAULIC CONDUCTIVITY WITH DEPTH The flow model ZOOMQ3D incorporates the variation of hydraulic conductivity with depth within layers of the finite difference grid This representation of the vertical variation of hydraulic conductivity provides an alternative to the development of multi layer models in which individual layers are characterised by uniform horizontal hydraulic conductivity in the vertical direction The approach has been developed to enable the more accurate description of the variation of hydraulic conductivity in limestone and particula
38. iated with particle tracking is circumvented by ZOOPT Weak sinks are commonly generated when for example an abstraction well which distributes its effect over the volume of the cell is not sufficiently strong to cause groundwater to flow into the associated finite difference node through all of its faces In this case it is not possible to determine whether a particle should leave the cell through one of its walls or whether the particle should terminate at the well An approach is adopted in ZOOPT that eliminates the problems associated with weak sinks This is discussed in Section 3 2 12 3 2 VELOCITY CALCULATION 3 2 1 Weak sinks A problem that can be associated with particle tracking codes is that of weak sinks In finite difference models groundwater can flow out of a node through a sink that is distributed throughout the whole volume of the cell For example abstraction wells rivers leakage nodes and springs are all simulated in ZOOMQ3D using such sinks If the discharge rate of the sink is sufficient water will flow into the node across all of its six faces This is termed a strong sink However if the discharge rate of the sink is insufficient to cause inflow across all sides of the node that is water flows out of the cell across one or more of its six faces then a weak sink is generated Weak sinks present a problem because it is not possible to determine whether a particle leaves the node through a cell wall or through the
39. ilable from the BGS Sales Desks at Nottingham Edinburgh and London see contact details below or shop online at www geologyshop com The London Information Office also maintains a reference collection of BGS publications including maps for consultation The Survey publishes an annual catalogue of its maps and other publications this catalogue is available from any of the BGS Sales Desks The British Geological Survey carries out the geological survey of Great Britain and Northern Ireland the latter as an agency service for the government of Northern Ireland and of the surrounding continental shelf as well as its basic research projects It also undertakes programmes of British technical aid in geology in developing countries as arranged by the Department for International Development and other agencies The British Geological Survey is a component body of the Natural Environment Research Council British Geological Survey offices Keyworth Nottingham NG12 5GG T 0115 936 3241 Fax 0115 936 3488 e mail sales bgs ac uk www bgs ac uk Shop online at www geologyshop com Murchison House West Mains Road Edinburgh EH9 3LA 0131 667 1000 Fax 0131 668 2683 e mail scotsales bgs ac uk London Information Office at the Natural History Museum Earth Galleries Exhibition Road South Kensington London SW7 2DE 020 7589 4090 020 7942 5344 45 Fax 020 7584 8270 email bgslondon bgs ac uk Forde House Park Five Business
40. ion method determines which numerical integration techniques can be used to define the path line These considerations are discussed next 2 2 VELOCITY INTERPOLATION Different velocity interpolation methods have been used in particle tracking codes of which a number are listed in Table 4 Each velocity interpolation scheme has its advantages and disadvantages though the selection of a method is often based on the comparison between linear and multi linear interpolation techniques The benefit of using simple linear velocity interpolation in each co ordinate direction is that the technique satisfies finite difference cell by cell mass balances Goode 1990 and preserves velocity discontinuities at cell boundaries in heterogeneous systems A disadvantage of the method is that it can produce less realistic path lines in homogeneous aquifers when compared to higher order interpolation methods such as bi linear interpolation However a significant benefit of the use of linear velocity interpolation is that it allows Equation 2 2 to be solved using a semi analytical method which is computationally efficient The efficiency of the method is discussed in Section 2 3 1 Table 4 Velocity interpolation and solution methods used in particle tracking codes Particle Tracking Code Author Interpolation scheme Particle movement technique MOC Konikow amp Bredehoeft 1978 Bi linear Euler integration RANDOM WALK Prickett et al 1981 Bi linear Euler integr
41. kage used to open them enables visualisation in three dimensions such as AutoCAD If such software is not available and the files can only be viewed in two dimensions tracks xy dxf shows the path lines in the x y plane tracks xz dxf shows the path lines in the x z plane and tracks yz dxf shows the path lines in the y z plane An example of the three different views in two dimensions is shown in Figure 20 This shows the advective movement of particles towards two abstraction boreholes When time variant particle tracking simulations are performed each particle s path line is drawn using four colours Each coloured section of the path line represents the distance the particle moved during one time step of the simulation As described in Section 5 1 the user can specify if the model mesh is drawn within the dxf files by adjusting one of the parameters in the zoopt dat input file In this case the model mesh is placed in a different dxf layer to the particle paths The dxf layers can be turned on and off when using certain software packages to view the file for example AutoCAD 6 4 ERROR FILE ZOOPT ERR The file zoopt err contains a log of errors that are encountered during the particle tracking This might for example include messages informing the user that the co ordinates have been incorrectly defined for certain particles and that consequently these are not created 35
42. l in Jackson 2001 2002a and 2002b and Jackson and Spink 2004 In addition to representing hydrogeological features such as rivers and pumped wells that are commonly simulated using groundwater flow models ZOOMQ3D incorporates the vertical variation of hydraulic conductivity with depth and local grid refinement Particle tracking within ZOOPT is compatible with all of these mechanisms At this stage of development ZOOPT enables steady state particle tracking under advective transport in both the forward and reverse directions The code can be applied to the definition of borehole catchments and associated source protection zones the identification of recharge and discharge areas and the visualisation of groundwater flow patterns for example Back tracking is easy to implement within a particle tracking code Steady state particle tracking can be applied easily at any time step of a time variant simulation The flows calculated at the end of a time step are used to define an approximate borehole catchment by implementing steady state particle tracking routine This approach can be useful to analyse for example the approximate change in shape of a borehole catchment during a seasonal recharge or abstraction cycle In addition to tracking particles under steady state conditions the code can be used to forward track particles in unsteady flow fields Particles are tracked using the semi analytical technique described above However around parti
43. l through which it exits the node For example consider Figure 2 which shows the path of a particle in a two dimensional model grid from its initial position Xp yp at time tp to the point at which it exits the cell Xe ye at time te In this example the cell wall velocities are denoted by Vx and V in the x direction and Vy and Vy in the y direction for simplicity e e Finite difference node Vy2 ya t Cell wall AUS gt Ye ste e Particle track o LJ nd MI Vx2 Xp yp tp e ud y y 1 f A KI X25 e o X Figure2 Particle track through a two dimensional cell To illustrate Pollock s algorithm the assumption is made that all the cell wall velocities are greater than zero Then if we also assume in a first instance that the particle leaves the cell through the wall at x2 that is in the positive x direction then Equation 2 10a can be used to calculate the length of time At that the particle takes to travel from x to x2 Equation 2 10a gives X X Lv expla at v 2 11 since x t x and V 1s the x component of the velocity at the point Xp yp Rearranging Equation 2 11 gives A x x Vua 7 Vp eXp A At 2 12 From Eguation 2 4 A Gu x V v 2 13 x2 and therefore by substituting this in Equation 2 12 we obtain 7 At Hefta 2 14 V x xp If it is assumed that the particle leaves the cell through the wall at y then by a similar p
44. nce that is assigned between the numerical layers In ZOOPT values of porosity must be assigned to all of the finite difference layers and to the quasi layers This is achieved though the use of two sets of data files one set for the finite difference layers and one for the quasi layers Input file names and format As described above the porosity of nodes in the finite difference layers must be entered into the model using either 1 pairs of map and code files named porosity map and porosity cod or 2 numeric data files named porosity dat The porosity of nodes in the quasi layers must be entered into the model using either 1 pairs of map and code files named porosityftta map and porosity a cod or 2 numeric data files named porosity a dat The symbols must be replaced by a two digit 01 to 99 number representing the layer to which the data files refer The upper layer is layer 01 and layer numbers are incremented from the top to the bottom of the model Either a pair of map and code files is required for each of the model layers or a single numeric data file is required for each of the layers The format of these files is described in Section 5 3 The user specifies whether map and code files or numeric data files are used in the input file entry method dat Note that file names relating to the input of porosity in the quasi layers are appended with the letter a after the layer number
45. ndwater model ZOOMQ3D British Geological Survey Commissioned Report CR 02 152N Environment Agency Technical Report NC 01 38 1 JACKSON C R 2002b Steady state particle tracking in the object oriented regional groundwater model ZOOMQ3D British Geological Survey Commissioned Report CR 02 210N Environment Agency Technical Report NC 01 38 2 JACKSON C R AND SPINK A E F 2004 User s manual for the groundwater flow model ZOOMQ3D British Geological Survey Internal Report IR 04 140 KONIKOW L F AND BREDEHOEFT J D 1978 Computer model of two dimensional solute transport and dispersion in ground water USGS Water Resources Investigations Book 7 Chapter C2 LU N 1994 A semi analytical method of path line computation for transient finite difference groundwater flow models Water Resources Research Vol 30 No 8 2449 2459 POLLOCK D W 1988 Semianalytical computation of path lines for finite difference models Ground Water Vol 26 No 6 POLLOCK D W 1989 Documentation of computer programs to compute and display path lines using results from the U S Geological Survey modular three dimensional finite difference ground water model U S Geological Survey Open File report 89 381 PRICKETT T A NAYMIK T G AND LONNQUIST C G 1981 A random walk solute transport model for selected groundwater quality evaluations Bulletin 65 Illinois State Water Survey SHAFER J M 1987 GWPATH interactive ground water flow p
46. nt of the velocity field only then the equation of the particle track Equation 2 1 is written Sx V or Lk dt 2 6 dt V Substituting Equation 2 4 into 2 6 and integrating between two arbitrary times t and t2 gives x t 1 t dx dt 2 7 Malan J l SS is x t and x t are the particle co ordinates at arbitrary times t and tp where Equation 2 7 is integrated to give sp S d TV A At 2 8 X J Vy 1 In 2 where kit and x t are the particle x co ordinates at time t and t and At t t Noting that from Equation 2 4 Nue A O a V 1 2 9 then Equation 2 8 can be re arranged to give x t x 1 excess Ar v 2 10a 2 x hy Equivalent equations to Equation 2 10a can be derived in the y and z directions These are of 4 y t Y a v t expla At vy 2 10b 2 yL 3 Zt ez gt v t exp A At v Q 10c 2 zL F3 d Eguations 2 10a to 2 10c are used to delineate the path line of the particle as it moves through the model domain However these eguations are only applicable when the linear interpolation coefficients A A and A are constant Consequently a particle cannot be allowed to cross a cell wall between the time t and time tz Pollock 1988 presents an efficient algorithm which eliminates the possibility of this event by calculating the length of the tracking step that is required for a particle to travel from its current location to the cell wal
47. of the finite difference node The selection of the appropriate wall to which the flow is assigned is based on physically justifiable assumptions for example groundwater recharge arrives at the water table from above However with regard to abstraction wells such an assumption is not justifiable because groundwater is drawn towards wells from all directions Consequently abstraction wells present a more significant problem To identify if a particle terminates at a weak sink well Zheng 1994 uses a special velocity interpolation scheme within the corresponding finite difference node This is based on the superposition of an analytical solution for radial flow to a well and a solution for unidirectional regional groundwater flow The interpolation scheme alters the velocity components in the x and y directions but continues to use linear interpolation in the z direction With reference to Figure 5 the horizontal components of the particles velocity are given by V Q va x a Va Va 3 1 2n z z e xz a y 2 V Q Ja Ye T v m V 3 2 m 2n z z 8 x a y 2 l where x and y are the particle s x and y co ordinates with respect to the centre of the node 13 V and V are the x and y components of the particle s velocity Va Vi V4 Vj V4 and V are the x y and z components of the cell wall velocities OF Quas 1a Va Guys Qwez is the abstraction rate of the well Ki X5 Yi Yo Z
48. omponents at the points x pi gt Y n fori 1 to 4 At is the length of the tracking step and x satt You is the final position of the particle at the end of the tracking step 9 Equations 2 18a and 2 18b cannot be applied directly because the velocities v and Vin depend on the co ordinates of x and y Hence the co ordinates of the particle at the end of the tracking step are calculated iteratively Pe Qu Y P A yA ji S A NO X4 P SD A e x y p 1 Figure4 Intermediate steps taken during fourth order Runge Kutta method The velocities at each of the points pi p2 ps and p4 are calculated using linear interpolation of the cell wall velocities after the co ordinate of each previous trial point has been calculated The co ordinates of these points are calculated by performing the following steps The co ordinates of p are calculated using Xj Xy Vy At 2 2 192 Yp2 Ypi Vy At 2 2 19b from which the co ordinates of the point p3 are subsequently calculated using Xy Xy Vigo At 2 2 20a Yp Yp Vyoo At 2 2 20b and then finally the co ordinates of the point p4 are determined using Xy Xp Vyp3At 2 21a Ypa Yor VypsAt 2 21b By repeating this procedure the particle is moved through the model domain until it reaches a discharge point The above equations are easily extended for application of the technique to three dimensional problems as in ZOOPT the particle tracking code develo
49. orward and backward directions It can also be used to forward track particles through non steady groundwater flow fields The particle tracking model is straightforward to run and only reguires a few input files in addition to those reguired by the flow model ZOOMO3D All input to ZOOPT is in the form of ASCII text files ZOOPT produces ASCII text and dxf files as output To run ZOOPT a subset of the ZOOMQ3D input files is required Table 1 in addition to two of the output files produced by the flow model Table 2 A further small number of input files are required which are specific to ZOOPT Table 3 Detailed descriptions of the files listed in Table 1 which form input to both ZOOMQ3D and ZOOPT are presented in the ZOOMQ3D manual Jackson and Spink 2004 Consequently they are not described here This document describes the development of the particle tracking code ZOOPT which is based on the object oriented groundwater flow model ZOOMQ3D Jackson and Spink 2004 Particle tracking methods are described prior to the description of how to run the program 11 TERMINOLOGY ZOOPT is written using an object oriented programming language Whilst the users do not need to concern themselves with what this means the term object oriented is used within this manual and consequently a brief explanation is required The object oriented method is an approach to structuring and developing software applications Instead of an application being bas
50. ped here Whilst the Runge Kutta technique incorporates a higher order of accuracy than the simpler Euler s method it is computationally less efficient Consequently it is important to optimise the length of the tracking step during the procedure to both maintain accuracy and minimise computational effort This is performed in ZOOPT using the step doubling procedure presented by Zheng and Bennett 1995 In this procedure the tracking step is performed twice First the tracking step is made using a time interval of At and then it is repeated by taking two steps of half the length i e At 2 The distance As between the two points calculated using a full step and two half steps is used to adjust the full length of the tracking step As Zheng and 10 Bennett 1995 state if the fourth order Runge Kutta technigue is used the tracking solution is 1 accurate to the fourth order hence At can be scaled as As Equation 2 22 is used to calculate the required tracking step size At that will yield an error less than As given an initial calculation of As using an initial tracking step of At The term f is a safety factor and is given a value slightly smaller than one e g 0 9 1 At f At RI 2 22 11 3 Capabilities of the particle tracking code ZOOPT 31 INTRODUCTION The particle tracking code ZOOPT has been developed for use in conjunction with the object oriented groundwater model ZOOMQ3D ZOOMQ3D is described in detai
51. prof iles username programs start menu hich is what you would have to type if extensions were disabled zoopt_project gt _ Figure 10 Example of changing the working directory within a console window The size of the console box can be adjusted by clicking on the icon in the top left hand corner of its window and selecting Properties from the menu list Suitable values for the width and height of the window and its associated screen buffer are shown in Figure 11 C WINNT system32 cmd exe Properties Options Font Layout Colors Window Preview r Screen Buffer Size Width f 20 Height 300 gt r Window Size Width I 20 H Height 50 m Window Position Left 8 Top 46 IV Let system position window Figure 11 Changing the properties of the console window 21 5 ZOOPT input files Most of the input files for ZOOPT are also required to run the flow model ZOOMQ3D Only a small number of additional input files need to be created to run the particle tracking model The ZOOMQ3D input files that also form input to ZOOPT are listed in Table 1 These relate to the structure of the numerical model to the aquifer s hydraulic parameters and to the stresses applied to the system The format of these files is discussed in detail in the ZOOMQ3D user s manual Jackson and Spink 2004 and consequently they are not described within this document In addition to the ZOOMQ3D inpu
52. rly Chalk aquifers in which higher hydraulic conductivity values are often associated with the zone of fluctuation of the water table The method circumvents numerical difficulties that are related to the de watering of layers in multi layer models The variation of the horizontal hydraulic conductivity with depth is defined by profiles such as that shown in Figure 7 A VKD profile describes the change in hydraulic conductivity with depth at a particular point in the aquifer Profiles are defined by two sections In the lower section between Zporrom and Z in Figure 7 hydraulic conductivity is constant In the upper section between Z and Z4op hydraulic conductivity increases linearly with elevation Because different values of hydraulic conductivity can be specified in the two orthogonal horizontal directions x and y six values are used to parameterise the profile 1 The elevation of the base of the profile Zorro The elevation of the top of the profile Zo The elevation of the point of inflection Z The hydraulic conductivity in the x direction K below Z oh rum gt The hydraulic conductivity in the y direction Ky below Z 6 The gradient of the profile above Z VKDGrad This is equal to the increase in hydraulic conductivity per metre rise in elevation Ztop groundwater head h Zp Tx Kx h Zpgorrom r 0 5 VKDGrad h Zp y ZgoTTOM P LA Hydraulic conductivity K x gt Figure 7 Par
53. rocess the length of time At that the particle takes to travel from y to y can be derived V Ab xxx 2 15 y yp The comparison of At and At defines through which cell wall the particle exits The time for the particle to exit the cell At is taken as the smaller of At and At If At is smaller than At the particle will exit the cell through the interface at x x and vice versa If At At then the particle will exit at the corner of the cell through the point x2 y This method of moving the particle from cell wall to cell wall is easy to implement and computationally efficient If the particle track needs to be defined in greater detail Equations 2 10a to 2 10c can be used to define intermediate points along the path line within the cell This is achieved by dividing At by the number of intermediate steps within the cell and then using multiples of this smaller time step to calculate the intermediate points using Equations 2 10a to 2 10c In the above discussion it has been assumed that the interfacial velocities are all greater than zero for simplicity However there are three other possible flow conditions that must be identified before the semi analytical solution algorithm can be applied These are shown in Figure 3 a Equivalent interfacial velocities Va gt gt V Va Vo Al b Flow divide in cell z S Va lt 0 a X gt Vx l E Vo20 HI c Interfacial flows in Vaz 0 opposite directions Va gt
54. sink This is because the calculation of the velocity within the cell is based on the discharge rates across the cell interfaces only and does not take account of the effect of a distributed sink on the velocity field Weak sinks caused by abstraction wells are dealt with separately in the next section With regard to the other model features listed above that can cause weak sinks the problem can be circumvented by assuming that the discharge to the sink actually occurs through one of the cell walls In effect therefore the sink is removed from the node and one of the cell wall velocities is re calculated In ZOOPT the flows between the aquifer and rivers head dependent leakage nodes and springs which are represented as distributed sinks are all assigned to the upper face of the corresponding finite difference node The velocity is then recalculated across the node s upper face Recharge is dealt with in a similar manner Recharge is assumed to fall vertically onto the aquifer Consequently the velocity across the upper face of the node is adjusted to account for this inflow Abstraction wells are the other features of ZOOMQ3D that can generate weak sinks These are less straightforward to deal with and are thus discussed separately in the next section 3 2 Weak sinks caused by abstraction wells As stated in the last section the creation of weak sinks by mechanisms other than wells is dealt with by assigning the discharge to one of the walls
55. solution methods used in particle tracking codes 4 Table 5 Order of variables listed on line of flowbal txt esee 33 Table 6 ZOOPT autput Bes 955 Sets e eO Pete en tuos i Tot leves te pub e Rut Lacs o ua Repeated 34 Summary This report describes the development of a steady state particle tracking code for use in conjunction with the object oriented groundwater flow model ZOOMQ3D Jackson and Spink 2004 Like the flow model the particle tracking software ZOOPT is written using an object oriented approach to promote its extensibility and flexibility ZOOPT enables the definition of steady state and time variant path lines in three dimensions Particles can be tracked in both the forward and reverse directions in steady state flow fields enabling the rapid definition of borehole catchments recharge and discharge areas and the visualisation of groundwater flow fields for example The program also enables the visualisation of steady state particle tracks that are based on the node by node flows at a specific instant of a time variant simulation For example this capability allows the examination of the changing shape of an approximate borehole catchment over an annual recharge or abstraction cycle Particles can currently only be tracked in the forward direction in dynamic or time variant flow fields Path lines are defined using the semi analytical method Pollock 1988 however around particular mod
56. t files listed in Table 1 two output files produced by the flow model are required as input for ZOOPT heads txt and flowbal txt These contain the groundwater heads and components of flow simulated at each model node for each time step of the simulation The formats of these two files are discussed in Section 5 5 Finally the files listed in Table 3 are also required as input to ZOOPT These are specific to the particle tracking model The formats of these input files are discussed in Sections 5 1 to 5 4 5 1 INPUT FILE ZOOPT DAT This input file is the main control file for the particle tracking simulation It is used to specify the direction of tracking and whether the flow field is steady or dynamic in addition to parameters relating to the particle tracking technique An example file and its format is shown in Figure 12 It is composed of seven pairs of lines the first of each being a comment line and the second a data line A description of each of the data lines is presented next Line 2 Two character flags are specified on this line separated by a space The first specifies whether the tracking direction is forward f or backward b The second specifies whether the run is a steady state s time instant 1 or time variant t simulation If it is a steady state run the particles are tracked using the flow field simulated during the first time step of the flow model simulation Steady state particles tracks can also be prod
57. the next line as an integer value If this integer is equal to 1 a pair of map and code files is required for each model layer If the integer is equal to 2 one numeric data file with the dat extension is required for each layer 25 Entry method 1 Map and codes files 2 Raw data Base elevations Top elevations Specific storage Specific yield Hydraulic conductivity Anisotropy VKD Vertical conductance Wetting thresholds Post wetting heads Porosity Figure 13 entry method dat file format 5 3 2 Specifying spatial data using map and code files In this method the spatial distribution of parameter values is defined using maps These maps are contained in files with the map extension For example porosity02 map would be used to define porosity in layer 02 of the model i e one layer down from the top layer The map files contain one character array for each grid in the model For example in Figure 14 the model mesh is composed of four grids the coarsest base grid two child grids on grid level 2 and one grandchild grid on grid level 3 Consequently four character arrays are required in each of the map files for this example model These are listed in grid level order and are separated by a comment line Grid level 1 Base grid 8 columns by 8 rows Grid level 2 5 columns by 5 rows Grid level 2 9x9 Grid level 3 9x9 a b Figure 14 a Example
58. uced for any other time step of the flow model simulation This is referred to here as time instant tracking The steady state flow field can be based on the nodal flows and heads for any time step of the flow model simulation see line 16 Finally particles can be tracked through unsteady flow fields across a number of model time steps Time variant backward tracking has not been implemented yet Line 4 single integer is entered on this line This specifies the number of points at which the particles position will be calculated between the cells walls If a zero is entered the particles path lines will only be defined at those points on the interface between two finite difference nodes In this case the path lines may appear non smooth when visualised Line 6 Two decimal numbers are entered on this line They are the Runge Kutta safety factor and error criterion which are discussed in Section 2 3 2 of this manual The safety factor should be a assigned a value slightly smaller than one e g 0 9 The error criterion is a distance in metres Line 8 single character flag is entered on this line which can be used to enforce the use of the Runge Kutta technique to track the particles If this is not enforced the semi analytical method is used to track the particles except within weak sink nodes where Runge Kutta is implemented automatically 22 Line 10 A single decimal value is entered on this line Three dimensional plots of the parti
59. v 2 17c where x y and Z are the co ordinates of the particle s new location X y and z are the co ordinates of the particle s current location V4 V and v are the components of the particles velocity at its current location and xp At is the length of the tracking step Whilst Euler s method is straightforward to implement the length of the tracking step At must generally be small to maintain accuracy This is because the velocity is extrapolated over the tracking interval A numerical method with a higher order of accuracy is that of the Runge Kutta technique This is implemented in ZOOPT in addition to the semi analytical technique The Runge Kutta method moves a particle over a tracking step At by combining information from a number of Euler type steps It is generally more accurate than Euler s method but computationally less efficient however the Runge Kutta method need only be implemented in ZOOPT in a few specific situations These are discussed in detail in Section 3 of this report In the Runge Kutta method the velocity is calculated four times for each tracking step once at the current particle location twice at two trial midpoints and once at a trial end point With reference to Figure 4 a two dimensional case the process is defined using the following equations Vick 2V at 2V3 H V Xi X Cart Baga Dt Yat ye 2 18a You Yn Von M st M 2 18b 6 where v vs are the velocity c
60. which are specific to ZOOPT and not also required by ZOOMQ3D are described in detail in the relevant section of this manual Running the particle tracking ZOOPT model is straightforward The following procedure is undertaken e Run the flow model ZOOMQ3D e Copy the ZOOMQ3D input files listed in Table and the ZOOMQ3D output files heads txt and flowbal txt into the ZOOPT directory e Create and edit the ZOOPT input files listed in Table 3 e Run ZOOPT 20 SA KAA AAA Ta KATA gt cd zoopt project Nzoopt project gt help cd isplays the name of or changes the current directory HDIR D drive lL pathl HDIR 1 cD drive path ae Specifies that you want to change to the parent directory ype CD drive to display the current directory in the specified drive ype CD without parameters to display the current drive and directory se the D switch to change current drive in addition to changing current irectory for a drive If Command Extensions are enabled CHDIR changes as follows he current directory string is converted to use the same case as he on disk names o CD C TEMP would actually set the current irectory to C Temp if that is the case on disk HDIR command does not treat spaces as delimiters so it is possible to D into a subdirectory name that contains a space without surrounding he name with quotes For example cd winnt prof iles username programs start menu is the same as cd winnt
61. y within a console window 21 Figure 11 Changing the properties of the console window eee 21 Figure 12 Example zoopt dat input file and file format sese 24 Figure 13 entry method dat file format iai eh etna oen icd evtl 26 Figure 14 a Example mesh composed of four grids and b representation of the grid hierarchy em taet 26 Figure 15 Example map file for the entry of spatial data within a layer 27 Figure 16 Example code file for the entry of spatial data within a layer 28 Figure 17 Example numeric data file for the entry of spatial data within a layer 20 Figure 18 Layer numbering scheme used for specification of porosity in a a model without quasi layers and b a model with quasi layers eene 31 Figure 19 Order of data written in heads txt essere 32 Figure 20 Example dxf files for visualisation of particle path lines in a x y plane b x z plane anc ain ah do De Msn did loda AA LE herren tr 36 TABLES Table 1 List of the ZOOMQ3D input files required by ZOOPT eee 2 Table2 List of the ZOOMQ3D output files required by ZOOPT eee 2 Table3 List of the additional input files required by ZOOPT sere 2 Table4 Velocity interpolation and
62. zoopt dat it reads head and flow data from the files heads txt and flowbal txt for a single time step only However if ZOOMQ3D was used to model the steady state flow field by running time variantly to steady conditions using multiple time steps the data relating to the steady state conditions i e to the last time step of the time variant flow model simulation will be located at the end of the heads txt and flowbal txt files and not at their beginning In this case time instant tracking must implemented in order to track the particles through the steady state flow field This is achieved by entering an i on line 2 of zoopt dat and by specifying the use of the last time step s heads and flow data on line 16 of the input file 33 6 ZOOPT output files ZOOPT produces six output files two giving information on particle tracks three DXF files used to plot the path lines for visual examination and an error message file These are listed in Table 6 and described below Table6 ZOOPT output files tracks out ptend out tracks_xy dxf tracks_xz dxf tracks_yz dxf zoopt err 6 1 OUTPUT FILE TRACKS OUT This output file contains information describing the movement of the particles during the tracking process The file contains six space delimited columns that contain the data listed below Each line of the file relates to a single position of a single particle Particle number x co ordinat

User's manual for the particle tracking model ZOOPT

Contents

Download Pdf Manuals

Related Search

Related Contents