- AI..- ... A .(.A undiw..'. .AA A . 0A.-..Z .l‘oA .0...- ’2“... . A A . A 3%. .A ...:..:‘..1.'.. ... ... ....u.‘..~' .xln.. . '9 ‘1 I .A.... 0.0 ‘_.A....IA “duazruflflfzfrio ... . . . . A .0 ‘0 311.10 1.. .X. ......‘J‘AOUAIA. I . .. _.£—u.fi.’ur...ni’ .A....A... .A. ......chc .IA. ..c. ".../30.90%“. . . II . I ......‘1 Q3 ..¥.-§._OAA~A...$A. 3..-... ..- . ... ... .33....Aommnsyfl. . . .. . . .. . 0 _ , . ..A . . .. . ....A 3.33. ...»...A“..x...n....; . O .‘ .'I A O .9331 A. . o . . Atop. . o 4: ’ 1A.".Jb. 'gAAzu. ..4 9. .u... s I‘ ’8“ 7 ...KA. 42v... . VGA-.4... ..s .9". . A ...... ..A .A. .9 ’0 J o J. . .. ...—3‘ AI .- OVIA... .000“... r (a a. ._ a. #2335. . p ulgi. : A A A 2 I .IA) 9.50 4 IA... AAA-.A.. A AA... .0 3.13593 . ....I ILA-Olfloo‘ ..Acfu’. . A. ‘01.! A ..I.—0'I.9....o..A.A-.... ... A. 040...... 90.0.4... . z . A . II .0... A .84? .l ‘0' .A Attic.-. .A.... c A... 2.9.7... .0... . at.) A .A.. o). .... 0.. A .. ... ... l ... .AA 30.... l..’Af.A..I..... . to..A . .A. . 03$. 0 3 .. .00. 0:91 AIAI'. I .u A 9'. v . “b.3§i.c . I.- otJh. .8. I... :02. 1...... .- A 0’..\A’ [are 0’7... At.- 0...o,.’i ..0 . . I. . I . I A. ......'I. a A 70.0. A 0A .- . .. 8.. ... A?!..AAIVI. A 9.03.3 AVA}! 9.0 .A. 15.13.“..2. . . ’22.... . .2......:.... 001...... v. A A.I9Awt,;3 i A 0A0 ...- ’0: A... 0 .020? . AA. ’{AA.. .A. z I 0.. .500: ...... 48 I . I 9. v I A. . ..I. A. 9.13.... ..AAAIO..A .1: .A. -..r. ......A ... . .I0!.,A.A.. ...! . ail. 4.. ‘A' Dr AVAI. O ’ .A. A‘rO-n A... 003'.» .A V. “A A d .I ...-0.! 0.....- .0". ‘0 . ... Ail-‘95.. .l'l..AA . .’ 1118.... . a .A. . A ......A..AoA—A...JIA . .A.:oo...‘ AAA . «GOA-gv ioo‘OI'... .. o. {‘30. . ...; .A.-A... . A“ ..I. "...-n... .. §.L O... . .. ..04 .oA .. .....QAA ..I u A 10.10. A ...-Tynan...” on”... .3. 6.0.. . . . .A. .... £39.73 5.... .‘.t AAA. .3, .7 Ag. .A.-$3.07. . . in “IA .1u4o00 ,J . Uo’ .0 A . . ,A .08». A? .3703 . CA AI. ... ..t":.. 7... ...‘Q. ’ I“ d.....' ... 4.. J03...‘ - '3‘ . A.- A ...-.A.." A. .0440 ..4 . A....AA1 9. .AA. IA_.«.A..A.I. .- fl A A. F. .. 4 — 2A.. c.9aA‘a.-4 C .A. .Pafl; . vdwt . A .¢A Ii . ..o’t. .... .A.. . nfi ....A R .t .8: I . A .A. .17... ... , ’1- .‘I.. .a .A...A.AA..... ,. .. . . n .30.: I. . 5. 363:! 0'... . ‘é.§f." our-i... 0 I... ... .. ASL; "’8. .A....A A; A. .0 A. A ..A: $3. ...-.... AA . .. . t . lOAHAN‘OI .u;1..o ‘1 EAA. 9111 n? .- ...I: . A . o... . v . A t... not.) .03.;31...‘ .7.“ . A .. '0 .A. n . . _ . .A I .‘o. .. .I a A .94 A A . A on": 11.0 . A 10.2... V9.51. cJoAnAAAOA. . Au .2. .... .1... an ..0A {I A no... .. .A" 00:; 10“. .c .. . AA...'I8 AAA- on.s’u..«o...x.43 .AAA‘AP‘ 0' . .A.—,AA'I: ... .A‘ a. t .1 ...... 0.. .IA '3. . . AAA). ..‘2009’ r I ...!!..Am. . 1.1.3 .I A. .. AAA...&A\...I.O.0I *A...¥IJ~...0J..A0¢¢A. A A ‘ A. .. ......Q.o .I": A . 0... :I A .. ..0."...1A!. A. .r on ”’9': ”A“... ...... "hu.. .0. ‘0 . . , A ..A . . l p o :A ... «XI...A. '01..) .. .A.} ..I .0. .A. . . . . ... 0-1... , I... . ,AAA.WO.LA‘J.”A.!0 :0"...A HAN-0.00.2”. .03.”. 3.0.0... Hog... . ...? 1.. '41.... .A.-AAI’Ah—O: 01-150 1'59... .1“. ..3 II! . A .IA ...-A. o.'0‘“-1Aoc0.’0o.. .A. 00.5. E.‘ I .5 9... "...-......A.... 30.0.. 0.. 10A... A.A..‘.. .A-.02.O'.Q.a.d ..‘I . - A ..A ....w. 2. ...-(dun... \Ig-u A.W.A.A?.. 3.4.1.. .A. Aft: a? A-.. .A.? . 4 .Au. :3... .0 “.....Ir A... AA.A.. A64... .AAI JAIL A. A ".0. $021... .....II....‘.. VALVAAHC‘A .4-.. :A I9A~Lh 0 A . ”BEA - . I . . A Col . ... .OA . .0. .A. 0.. 6"; .l . .. so... A. I. ..l . . . IA . .. U :0; 9’ .IhAVAv. Olindcuwh“ . .I‘JH‘ 0)..” . .. AahhfluA'uo .0... .lfltVJo ArlleAAQ AA... .O.l .- ..A . at... A. .. A >1 A . 7. . ,1 2.0.010. A. _ ..A C . ’w-L .. .o .QACIIPIDAI A ...er 'A ~ .A '0. 0.. .”.AA 0v 1' ‘6 ‘0 .A .u. {.A.}. ... ...... .. ...»...rif .1. . .A : -..: .ur..f..§........... . . .. ...}... .... ... .... ...... 933...3......A...>....1.a.fi..u..”award...“ “3.... .... .A.... . ..- .. ....,......nu.:... .A.... . A . . AI .. i. . A . A .. . .. . . . L . .0 I A n :A; '0' A . I. 0.. .I . . . A. . I. I. ..A -.0 . o . c. .. . .‘wr..XAA.Vu.«.A.wo I :19. A AAAA.auH..cuSAnAAh‘ ......IA Vuhhmtn”. w.23...."AuquuuAALO......-AI 0A9. . .uA “If... ..A‘. MuAfA “3.40! as. 3.1.3. AA... .A.. 1.]. I'........... ...... {$2.10, .... ...... IO, AH...A ......AAFNAAWAs .A.-v. W332». A7... .0. IICZq. 1...: A fih . A ... .1."¢C.!. c. YO.AAQI.‘. .IA. . 9. II. A. ..‘ISQAAA’OAA. 4....IA’ A2030 '30....2‘ .00 .A.A..0A .(AAOVA AAA.<.$.A ISIAl‘ I...A...A. ..o .A. A"8 A '1‘}... JA‘A...A .IA...” < .AQ‘AA .. .A ...0-7! ......I...~0.A.A.0.AAA. 0.]... AI“ AAA \0 .93.}3u. ion-...: at ’37.. AAA... ...“..‘1'... ... ’ .... anti.’ 3.0.3.3... A32... .4 I 4"! ....AAA 10., Catt-.A. - no 60...:- A ..m? ... A AAA. .0 $33.. ... .A. ...A, AA. 8.. A- . . . .A. A... . A AA 00A.A . A0: ”A. . 0. .IA 3 . . A6. .A.. ... $.01. - . Iii . .... 0...... ...-321.005.. 0—3.. 3811.... , (3:.A2..\..AA .A.A .. oAI‘. Av . ‘3 '3. L LA 17"“... It... 9 A I. ‘37.: 0013.30130 OJ. 2- .....AIAAAAIAAAA‘...I 3'30, 3 .. 3.1.0 .A . .....P... J‘..I.sc..:.A'... ‘ ”I... 00 .. .A..”a.’ AI 1‘. A]... )0 I .0029 Ir. $01. .A. . n..- ..A! A. .36: .. a 94L. d?~.ni£i0¢. ‘A’ . 0.0.3.1- AAfidflé...v ,IA......¢.... o.17.0h3.1¢h.?_.o .7 11'... ......AA... .0 A... sci-9.9. \iln... (I!!! ‘30.. 2.93.3 .‘uo:.0hm.A\uh gnu“. u... 9. WNW-2AA.- nWA.?AIOM.A.« ...!I,.uAA..HA'tAu.AA..-na LL... ..Yn.0.<.”ou.ro A...» At... ...3 . .31..“ A.” 90 0.0... A. .A.}.LAOINflIAfl-HA’A AA A A A - A 0‘: .J‘ Ail: I . ...... A. 07 .30....AQA...’ . ‘ .A.; 0A 3A.}. .o l A ‘A . II b. .‘0... 3:0 0 A 100 A .-vo k _ A [CIR I A .A . . . _ . . A .. .3 . I!.-.. ‘. .. 1.1 . . A030 A3. . 0., II . 1 .o .0“ gr. . A01. 3 . . . 4 . I . . A . . .. . . 0kA .A . . . v . .f. V . AA ...... 3.10 5;! A. In. A .A c A 3. . in! I v I . .. Azibv'o A .. .A . . LA... .... r n»! . ....x... . I : ”......r. in” . ”uric hard ....un.uk.«r.AI.!A.n.Ae.H.I..rr..me 53..."... A an A An- ...nrr ..A .....A: . .731... .A.». “.A.“..uvn...’ u. um... n... 39.2. . .3....:. . .31.... ... . .... .. ... .. as... .53.}?! .. .... .54...AI:..n....AA. .. Azrtio. . . ... AuanW A. .... food“....b..umc.v.« «A A 9...”? . .90.. .72.. .A.»..Aol 9.0.010! A .. AM: 3.2.2.... DAJAOIO. ... 3A3!J.§'.3.ofi. I5 {I}. 30.. ...»! ..A. AA..~AA A " “‘7 An. I... Ayn”? m.’,, 33.0...“ ...-“o"...r: A0 Agi..n A ..AA ..HHHAAioAAAMmr. -hwlh .0...An.QA-JAAHIAIA.AJOAAIO¢A..A. IIZLoA.0AlI.A.2.’ .9'00. .10 “3401-“? ..u.s..lw1..0_01!.ou. .. .. , .. . . . v . . , . . . «A. .A...vA_P. :07 .iIoI. SIC-I .AA _ ? AI ’v .. .A.. . I V .A. .. A . 61‘” A .72....) . - . l 3'. . . ... All! ..A I... L a J ..0W .7 . .. 35. A: 7A. .v "A: v.d.\.. It $0.3 . AA .A.. ”A 31:. Aoqqouk IA. .‘AAAA . A 3.3.233. .1 AAA? 0 A I. A; 9:. ... Alt .9. ......s A. . fAIA ... 1 -I A» A. .A .15. I.AA ....A .....Iof . 1.... .IA It... 0‘ A 2A».' .0 4.0 . i .0... a“ 1 An. ....A‘ at. . .A. o I ' .. 3.... I o d .A‘ZAJMAZE. 3.. .. .IAAA ..AlfOAAo ... . . A01. .30. . . A 5. . .20... . . 4AA . A A. A.. . A . AA 2); AA. . 0AA. . .A.."QA‘A ...,A. bf} '08.. ’A . P IAA . I . . . Qty .20! A. .4 4A AA. A .A... 3...“ o. .A 9' A II :0 It... 00...... o .‘v. I: ‘3 I A" o 0 ’ .A.. .A . A o o x .4 .A. .A A A , ...AP AOAAQ‘ A .1}. .A. A . A . ...: . LI .9870..1n~o!.?: 0’. 0. Arflzl* 3A.}. .0. . . A . . . . . . ”A . Lit-10 .AAAAD ’o‘r..\c A3 I. A P: IIAAOAa‘A-A.’ .6 3‘ CI . ~33! A‘A... . I . 0.051....3 1. 6 -.. tr. 9:. ECA..’AA .. . A A0,... , . . A I. 15.0.1.0; ..AuAi! AA... 017.03.31.0C‘A‘Awll ... ’ Arm}. ..fl 0.0.- ...? A A. o.‘ >.! A. 1-21 .A .05 An 2 A9... ._ . . . - .. . . . ... A .... AA. AAA. 3.3.0.1... A. .A. a. .A.... ..AA. 012‘. .A \. AaI .11; I . . . .A.. .. o.Ic.II..A.H ... JAIAA _ .. . .A a: As _ .c .... V, II? .I. ova. A... A .... . . . . .. .3 0. A A .A «9'- QA 503...! '18:... .IAA’IAIA .A.... 9.3.3....1 3.1440AA .A....A..In .A. . sh.’¢u..\’ ... Cur-3.. P0 . C. ~ .‘gL.V.Wo.WAIAINHAA ~_ . .... 9A.. AAAAAAAI. 0 .!‘A.. 0.0.- .A... 3' . . AA in”. ... J0 .4 A .4 on .. .. - Aer .3"; ... AA. A. .A., .1... A.. ...! .0120 lot. .AI. 1.3!}... . A. I. . I. I 0. . . A... fll , .. . DIAAAI: A 40.00.. A... A .A 0000.. 6 1A.. . . . .. . , . . . . , . . .. :13 ..A. A A. D. . .II~ A A‘A Io - A A ,. A‘ AA. . . . I I . A . Al I... 0.. 32! A53 . I." . .0€.3N . A. 51.. ...A.. ._ ..A JAAA . . 3.3 . ... 6. .I . . . . I. As AAAA {'4 . 0.8 .A..I- A.) .... 03.1“.5. A-.. ... . 0A.)...8. ....A‘ A... A. .A.. .A.. . so... I‘ . 0.2.0 1329.. ..A .3 o . . .. . . o .93. 1:33.51 3.54.2.0. ’97.).03.’ . . . .. . . 0’. .... a . .A. 04. . . . A ... A...I....IA...AA.‘.Q.A.‘,..A.AA¢AAY, .ArAl. . . . .. A. u I .A V1.00.“ ‘hhai ”ob.“ A .'v ”0 vhf-0,00 Jr A. . v. 1M".AV.AAu.)AA.HI.. ‘IA.-.}0.0.fi .1” :0. ‘40.".I}‘A‘.D‘~.o.’ A. . . . . . . .u 6. Av ...A-m“ ”1.. ”Yr.&.tm.fi.oh(".fl.“"og..n A...I “£0 . AAAI‘M'QPLAKIIIO Ina”...0.’4' 50.". 0A. 0...... .049... 9A.} A ..I’.I A” Isa-.... Y‘hVi- 0.“...4‘”... 0t c'l A. A. ... Abs-1 DAMN.” ... n r104... AA .....A' . J A: . I .A. AA. y . :3“... {Heft 49.01.}. . . . . . . . . c. A A... 1‘0... Av .A A. A .A .203... 31...: ...-A52. . .0. A.. . A... .01.... A . AA... A A . .A.... A. 5A.... ..A .A. .... .I .A.: ... 0‘9...“ “A. .35. ...A A; '91....04 U ...! . 1 . ‘16-. n...‘.. 2.0:. A 1. .f. A‘. .A ... A." 17’ A." .V. 0 I) a {Ail-Amt. o... I .A. . . . . 4&1. . ...-0A .1033...“ A .....anA AA . A I‘d". .\I..9AAA0.A AA. VI QAAA. .10: 0.1.? ..CA...‘ Lin}! .90.. .7... Ax AAA ’.‘.V I . A. IVAA 2" 0.90" AA ..‘AIHAAADHHA. fr. .0: “...! £135.... . .. A «A (.70 ......A {on A... . . slut é\z.6.A0AA.0A 3 .. '1 A. .. ...—0A a. 7 A ‘. 0.1 I .50. . . . . .A.-..OAQAAA. A1.$Io \AA... of. 01.. .0 Av. AAAA .- 0!. A'? cl.‘ .... ..I (.A.-....AA .’ 1.. 3..-. A ... .0 a... a. A... AA. —.v.-. .A. . . ... O AQAAaO A 6A.. ft I. “.0 .. ¢ III. .A),’.HAI,I'.’ . 0‘ .- A.. Jun». ‘_Al.8 LWAW .IA-...JV.‘ .....udAuIAoquA. 13¢...f. AMI. ..CflTa. .A . . .. 0...”...kf. $104 35.0% s 3.....AAAw‘ . . 3‘ 1(..AAAA-AAH.NA.Q::.. ”’39... ..Sif‘lviu07 ...“..I0u... AA)... .23.. int .A.. .00) .A...."fl.... A! 6...... 3.011 .....Lf ..AA.. Nani?” .33“. ”...“...lAr... 00$th5.“. . AluRAAAAu. .-....1, ...Aarn . ... 01!. A...) A. 0'19”».1FAWIIHA. . 3A A; . I . . C . .1 I IA ... I. .0) .- A. A15- 9. . . [1.0 an. AA; .A.... . J AVA’A‘I- I :vbozz’AocA-A‘. U a . .A.-05'3" 7.3.0 ‘33:; A. . A .. A... . (AA. .A.. A.. .. I? A. . A. . . . . .P .A.. A... A .0 0! u I A ......0 o 4.”.6' .1. 0 .A A 52.0.1... 0.0.. AA. 3(- QA A.A. It: (.0090 ... . .0. on ‘?’..3I.AA 0.. AA. .A.. .A . I. L. .. a .0 s ‘A . 53:»... ..A! OCAAA VI A .V .I .A.|A’. .A.»..AAA AA-I. A..I. .' .A.. .Alo' \uOAIA! . At '0-‘ A .‘. .3.0. .A... .. .. . . AHA. . .....7. . . 0.... A .A. .8433; A...‘ A‘AA‘AAI . N70“ A£A.o A. 0.0... .A.. . ...-AA. .. ..J 0..- A . An , . .A. o . A 2 34 9 .7"; o -.4.‘ I SA. .A. 0. .1 A .93! .. .A. .A.. Au? . .. .15.. A62... . A .A I... ..."...Aib . A 0 .IA. .AAOA. ..A!.I!fl.t..o..ut .. A. .AA . . . . I . . . .. . :QAi, AA ‘. .IA}..- N :5.- ..A... A 9’1A 5’0... 1.. AUX-.A.: :3 A I! 1...‘ 04‘] .A!‘ A. v IAAA. a . . A . 3.3;? . . . . I .A.... Alli. 0 I... on .A. . .. oul!.....oA.0A...r.9... . b .. .1 . Ali 0. X! .T’n .AJIAA... IAAI..XAA.I'A.uI AchuAI... 0.. $52 .A. . . . . .... A. . .A .. , . AA ..AA A. Amiv . . .19..-.A'... 3.....I8'....!2.A-. 051.....Q... v.1 1...: ...IcAA. A . lilr..I‘I.an I. .....k. >p.. ...A. ...A A.) Afr-fllulr ... .4 A. .A . ALA. .w1AI ..I. 3.9.5. A... 34% AoA AA. AA 94! AMAAAOAJa’Am’A . . .A. . ...... .A. .. p.379... AA «3' ’3 O! .A. I. ... I A ..I A .A. .. .A 9A. ... . 01A...A A’ . 0...! A. A .r A I. .. .AA .A AAA. 0‘. I.. ll I23. ... A. 1.7 IA. . O .01).... Aiv... A... A .)J.>-’.IAA..P . ..AAH‘H.. .nA. A3... .«I. o. .38... . .. ..ASA...’A .. A .-I. . . A .A. .A.... .A .. . QAOAOAAJwDIIA. AJ-..¢JAAP. ..L AJ.)...II I... AvaAhAAA A.” A“ .A.”.AAH..0nl..AA.0 I’.0A.2‘..Aw... A.I9.AA0..0.I0II A,.AL0..II..A'O.A9A.M A.. ......2 art... A .J 2’23... .AAA ... 1;”, ’0'} J _ 0 ...r .04 on... .0. 14.. IN...“ Q." o’. A A no. “I. IA‘ 0... 4". A“ . .P'. A‘AJ HA0... OC.£I ! mad-000w“...- . . . . . a 5‘“ ic3~cvl90A ‘9. a 0 |.4A:l A. ”3.... .IA 1!,“ .4 x .0. ...DA .1 A 0.. .....‘ra'l [10. A3 0". DAOA‘. .A. A A. A . 5... Ar 00A...910’v|1 ..l‘! , A. .‘o .A..JI .‘30 A! '00. . n.8'431 .IOQN‘” I. .‘l:IA.. Al.- _ ..j. “hf n0.‘......4 Au‘ J’s‘. .I 40.." ..A.Ao\.0, oc"5.§.0’ 4. A . . . .. . .A 0.- Q... . I .A. A...’O.A.A ‘W ‘0'. A... '1 A A 0.. .. Afl . A . .9. n.|.'..- .3.» . .AA. .o‘...o-.O. .- . cf...“ ..-. 1 l . b 2... . ....Aa....eu.n..._m...n ..-... .... ....u.......L.r........... ... a . ......hfntuwawfifi..?n1m_ . {...}... Amunrfivfifiukd- ....n..§i..:.u..w£.a...”v3.9.3............A.W.A.....u..£.......-..........a.... “nix“... mafia}. 3......9...»$A.J3........_u .Af......w....... .... . . o I AA. 110-?- I. ... I A .. . . .. A A... 1.. .. . A! ... I it .IA . 9.3.. -... , ...-At 0.2 .III I... ms‘AAO.A AA .0? . A A A ‘I ...! . .A.. 9A AJI.A»G.. o.- ..A.A . A. .A..A .. . . . a. . .. . . .. s. .o .A. .A. 3.0235 .A.... AM! A. . I. AA A. . .A . A . . ..AI__. A.. . . . . . . .. . '0A A... . A . A... u . , ...-A . I ....AAAvA .... .. Aft-4.131,. ..AA...\00 .?I.Ao 04.01 . . ._ .A ...... ... 2... . .. .. x... uni"... .3... .u........... cram...u....“..ui... .... .. . _ . if”... ....A... ..«fl...3.............u-n..r ......k a an... I. "an... ..u...“ .....1 ...- . .u 4.33.“! hr... ...... z ... ......5... I.......n.. .r .....1. ...3. ...... L. . .... .. 2. ...A to}. .. A3. 1.... A... 0'8 .. A ..A {...-1.3.1.0.3U .3112 o. A: A... . . . . .A.... A.A A. . pout“. {.A. )4. A..A3.A,fl-AA .A. .0 .1 5...... 0A. A. .A.. . A 7 f. A A... *no ...-0...... . I . I. A ... I A 7. Li. .....S . . .... ..- .U AL. .P, _ 33.1.}... a”. .A.... A. . . L... A ...‘r A... 2.047.... .A. ..2, a .....le H.- A: ...“:6‘!‘ SJ] ... 1A ...-I93}, A910...A..AA..’. "r’cfi (I: .-....i’fio. A 6A ... .....A .3123! .... .A. .A.; to... Aha}... .A .... .. A. "3": . ....AAIQA 0 . .0 o 33...: 3...... d 6237 ... .A.. ...-Am: J 3. ‘. .AI AA .3 I4. .324 73...... ..A. A AIIQQIAIIFYB. u. A 9.... 3.553.319. AAAAAI: .1 I A. II...L.AJ.A.A ...: ..I A. I... ... .A 0. .A. ... . A. A. A .I. . .IA: ...! . 33...... .l. 9 ... . .253 o Q ".A....JA 22:294....36. ”Atelirnviv «2.032.03AAAAAI3A‘AACA 11". A. .A A i A. .‘00 9AAA .C Q. . {...-......P’QAA’A. ‘3.’oi!di¢ A...Al.rho'¢!md.0o'1 .U‘LKI.A“"A. 0).. t IzA.AA..Nn.Ner~A9.A .JHAIOAH \HA..POIFJ.".JA.AIA...PAA.AA"JUIIAH’A'.. ”.30.!“ I, A . ... .IA {.A. ..- ....AA. 2 ..A" AA 0 -..v 00A. A'AAAAIA.A.A VOA. - on? .5 ‘A- ...... . . if. . ...-o 3“ 11.15.09 0 9.1:! .koAA . .a on .0. . . .v I . . . . .A, . AIV'A. 0. A . . ..._. . . . .A . A o , ... . .. .. . . . . A . . . A, A. . $0 :32... ...OIA'..,.~......0. '. ..’t P a". . .... .... ... ,n um...» ......4 H.915...) ...... ...»..- .....u... ... Jfiran. .. ....-.ut....vn.m..:...9.nw,mvv . . ”...... u.............n..é.w . 2......» ..J. u”. away”... .... . . .....-u...%...., “.A.... ..."...nkr... .uquKJI... .... henna? . . .. A . ,. .. . A AAA. _ . . ...-7.72.7...‘5Aooz... A I A... s 7. ...-A. ?. 2... .c .A.... . .333. 3. 10.... .. . 1 AAAQAAAI. ...... .‘A.P.. 11.11%... A." ... 0...-.. .. . AA. A ..I - . , . A. AA... A 3.1.0»... .....I.A.ou.ro.u,u...a-I. .. 2*? .A.... 40 .A It... . l. k A. ... .....AA. 0 :0- r. A A . . 2!-..‘0..§l .93... A. .AJAAo... J .... ...llllrzf vol“: .A.-AA..." .0. A.‘ o c Afi ' .I- . :DQAA. 0.4.01 ’O’AL - 0’. ”AV“. '0 o... .0. A.” . 1‘ A Act‘s .01 1.00.0- ... .0 A CQfi'O’OAI ”0 c. _ ‘p ......I [A QO’CFA A...AN.’ Co. a o O. .I;.I‘0’}0 0”. I: ‘ Agni-......AAP'fJ: A.. A A .‘I‘.’ .30. A0".. -... ..P. .A ...-0? A... .A. . ... (A A .A.. d .03. . Ah: . . l.!. 31%.... J Ital-0117.". . .I .37.. . . . . . .5 I a . .0 .9 .033 ...... .AA. . . ASA "I. I... Q 1...... I27. .A.. ?A 192.31.. 2 .I. A.. .I ....-. A ... A. t....... 3.1%.... v.9!- “ 0 r5.- N-AA. . .r. A . A. ... . A“. ..7‘ .I- O . In... I . . . is ....o o. “no! . ..I’A‘...’ .0. .A....AA ...... ... 1......” .. .Au A4. ‘ Y, . A AA 5...! . .8 AA. -A a .A. .I A 9 .Acirv.A.1 .A 4.0.”... ..I.. nL...‘ a... ...... 20" «...-o In. A A .. . “ILA...“uufmr ....A LCMF A. fldrc‘oroza. I a..g"A.0vu_u0A’AAfAn.nI!IA-. A o . _: «...-.A.), . . .IA! ".AA‘ANA ......Al .Anfifwwcun. ”Honour“ VA”. VA .....A...A.”AAU A Loy: .I ....An A! .AAA... . A ... l...“ .A.... h": ..Am Li A .‘Acua . A: ..M... AA 3‘. .A.... . 3.3.1.0. £11.31... III A33 "A... .3777... .IA! «.3 A0... ...-AAA...A.. .Al. 0.... . AAA?!- VArA‘A .A.. . A . .IA. .A...A.Al..fio.f .10.: . .A.. 0.3- A... .IA _ A . . ... . in ...... .A 2.12... .A....z.’ J... A. . ...... .A. .A...A. .A.: 7:. 02.2190. .A ... 3.5.. AAA... {...-T . .. . . .33.? AA A of... .0. ID. ..AO..'!..!A..\ $ . .. .A. .A.. A ..-... . . ... 5.....0FI: A .IA .9.‘ A ...A- to... .aI .o . .AA 0.1.83.3... 00...}... ...A.AAA!..II .. "....IAAAA AA: 1:. ...... ,.A A . .30. A .. ‘0’... 9.... AAA .AA.A.I\.IAA 2! 2 .A. A...“ 0.. ...AAAA... . A; .. .95....2- ... A... A. ....A: .A ..A. ..A .03.}. ...!!.$ ....A. ...I.......M...A.AAI :3 .A....Tflt.’ ..- .I . ...Au...!....A.1.. . 1.2.... .A. AA :5"? -.....z. A. A... .. s. "VIA”. a}. Java." no: It“... . A 3.0%! . .... .... A ....A. . woof... ...: “A A.AI ...... It}? is 1.3.7.091... :1on .A.: ......t. , . ...: . .... A. .. Itch“... . A .A. D .. . (.....- .0». «SI 3 I 10‘s.... .. (.20... :0 A... 5.52.. 0’ a I o. x I: . .A. 3AI~A.!.8.AI:..0A\A..3. .A.. AI 19.... . . . , . .033. 119...!!! AAA . l . I. 0AA . . .4 .....Q.... 0.35". 2 ..-.1. it . 235 .....o. A‘A..AA..A.’ . I.» 4.. .. .A.- 2.4 ....o. ...-4,. ... AA... ...-(A... . . .. I . L cl. I A]. u ...A . 1 I. . A01. It ..A 6...... {A ‘2 ......cl. ..A...o¢..A A ...... A...IA3. 101.. n.3,. AAA-a"... IA..- A .....I. .01.: .A .... . .A.» . . a. A. . . .31...- IA... A... 9‘ ....A. .00.. Kl...IIlA ».(A :5... . .....AAA 113.. ...-.A.: AAA AA: (.6071. .3.- “f... .A\ A6 £5751”. . lift: 0:. .A ....ap. .A. ..l: ‘4; ...J... ,IA......AI.I. f. 1.11.0312... c. . . 13 A? -I “h!” 3... AA“.W..UAA. .90.. ..Al. .. Om... In... ..fl.’.AH . ...... . fail—A c.0020V ”IO! . .An.: A 3:: 07.. ‘9‘ AA. A . .I ,. . . I. A A. A’l; .l .. .U o . A. . A . ...)..{AIL3413 DOA .m!. I... 3.0.5! 3.0.". A. A . .3. m ufi§A.DOA..WAI .”.3-"05’.A'.IA‘" .01w DA... ‘06 . . MAT}! .A....O ..J. .r. I... . 0‘ .0..\9w0¢o....|..,. ’5'.‘ ..A It. 113‘. A... a .A. .9. AA}... (...-AAA 0 0:9 0... It .03... A3». *I. .1. IA." ”fluid... . .A. . ...»... A -. I)?! ... )v . )1 AI: . - .....A I . Holt! A, A I. A. . I. O . . ‘6 ...-0.0.17. .A.‘vo..o‘u .A...‘ A‘ ...-"u". 0 0...... ... A. 3.33.001... . . ' .87 A... 3"“? .0}? o. .. . A A 3-5.1.170. the." «A 007i! 3... A A... . A .A. .A- .V\ . 0"! A. 0-... . A A .0... 0&0:- ouv .' cl. A. ..3‘1 I‘iAA. . A . A 01!" u.’..".1 ‘ A A o *30.’. .. 0- to 8AAI¢ In I AA. .0 0. .....u. “an... .. . ..AA ...,‘7 at... . . o 0%-, AA ... A AAA; .49.... I... "M3,: ”.0... .AQ .0A A. ‘0. uu.... . A. .§A A A ‘0 .‘X.Ao.i ... ..n A u. (a a A A. . A... .. AAA... .30- Ad: ...... Hm. A 33.... . .0. A.” J. ....AA‘AQIIJIISA ARIA}. . . . I «33 ’ . . .— ....z ‘u” .A. . a... . A1. 04:... 3’... .I .00 .....I . .. “a A A AI '50....4 0. A- )2} o . .0730. .0 A A A I. .405 0AAAA,...A.A.. A. 239-. A A1;- 0.1A.’f0AH.I. 1......10. A. AAA. A r ‘01... . .0 I‘ .A. 17:... . 0A.... A. A. Q I'.‘ A... A .A. - . I ...... Auras-01.... A. ..A. ‘3' A”..PAA.A .002?! D I A A. O. I A. A]... AA -4 O . .7 I ...! A’AO...0.OA . .A. 0 .9303..AGAA.AAAAAAAA. A 0.- A... A A... A 3 . .A.... AA...A.AA« .n..:c\li ....A If... A... ..f A... A ..flm...0I:A.-.I.A.AC. A333... ....I I. I I . .A. . AAA .0 A .0. I. .39 A... A A. A A .A. 0. I .3 0...“: lo‘AA'Izld 03Vv.9‘ .A 1... A A. . .A. AltrA...1A If: 90 . A A .OIoAl_ . . Am .A .A. .I A A. 09.0- ...A 3,333.; ..A 3.0.-0A £1... .- Auam A AIA‘ f3... .33.... A I. .‘oAZ A. .A.... AQIW' .20. u ”A. ... . 4... AAA. ..II A A AAA. -IA. 0”: A u . . AA. 0....AA Iv $03.2... 13.040? IAJoAg I 0.. ..AO‘;50A.A A l IVS. A A. ...... IA. ..0 AAA, . ”.07.. o *4 o Aofp..9.,.4nh.‘¥o, '3.» AA. A. 0’2 . A . AA A . ... ...A...10A.A..IA.A..:A.7A. AA. 1., . 0...... C .A.: .01.: A... .A.. .n A. A..." .I L... A . 3 1.5.. .A.... {Ali .16)... 0AA. .v’flnr. A, 0;: A’.’ .I'AI A0310. .... A. 0.! ..Ayvoxo‘ AA 230.. (.0. A. I . .... 0.. . ...-310.91.30 II A . A . . o. . ._ A .A.. ..0I000 1.0! I‘m. A ...-0A..."- _ . T A V.- .o A... .71.... . A A9... «- ...... DA .A A. 10".... A .. AQQII . ......9a A ...: "SA. A .. .. . . 0.5000!- 0... A AA 0.09.0... .A. "A. .3; .890. A... no... .~.AI.A.AA..1. A 19-..?!A4 .A....ai ...‘AA.A..' .... A .. .. .0 .0. A .A 3A.. . .IA... .. .6123"... . '1!!. 9.. .A..- .d A. A A3: ...... .A I}. .I '9? AA ... . ._ . I ...... ’I...... A. . .. .... O A .A.».AAIA... .3... . A As AAA}... . 9...... 1.... .. one. ... . . 0.60. A . 1.9.21 7.139.131; A. l. Aowv A... v A. A .A.Q.I.9AA. AAA: A4... . A? .‘i‘quit.’ .A. A. 99.1,... :50: ...... 9A: A . A J.) A. 3-9.... 30:03.... A92... AAA. r10 .A.: (3.10,. .Aaooc .A.. .’3.8.. 9 .A. .A. A AA... «A: .A.—09.9.... A. .A. A.... .7... 012-6.. .A.. .A. . - .... .....v..v.|A . .I .A 2 .A A- I700 10.. 5.01.0...{4 ..AA h 42.. ’37:. 4A. 1: . ..PI.AA.A.X.I.’OA c .0203 . .003 0... .A .. A .A .A. info... AFAAnoivoi ...: . .. 0.,AAA ..9 ....-. 0.. 3"! .11 .11 .0. 8 A A. AAA). “AliAf. 20A,: . . :3 A... 1 .A..) .9. In... 3.3.. 39...- ... . . A 0.... 0.... ..Av‘... 010 “A. h. (A. ...? A. AAA... 0.. o .O‘aS.II .QIOC A to, ..A Al. '1' AA ADJ . .A. .A ... ......QIA ... .AAA. .... ...... 9 . .74... Av ....I A? AA 0. .l on: A . .A ...A 23.90.... IA... .09. A; I. . .. . A . A . .. 0.. .....IA. LAAAA A. .A .A.! . II... cn".P.OA A . 9 ’..AJ ".0...vm AAA. Ad... 0.11:3"3. A v’ .3..1A ‘1 L00. .1... A.<‘.A§ 10.... Alice. 0-. ..z.. 00 A J 0 0A 9. A... 00A .2 . AAA. ,0! . .o) .2 .IA .30., 0‘. .a-olc A; I o. A.» .r’AIL’O‘... '0 A... 00.0.”: . 90A b . AOAIOA'A.0.?~A..I¢D A. .9. 0:, t A _ o «J .. A . A H can“. 0;. .0! . .V . ,I .A . ...A .A 3. - OI. At... A ‘2'.)- o .0 . A A .OAt gAu‘A 7.7.196 .A.... ...;0 AAA... ...... ., .. Io«0.r.'.... AA A .A.. . A...\n..I..‘.o. g’3. 30¢... 50 . a .0... P .....A0 12A . A A: 24’01... I. .- ..IQA 9'. I. .3..0.A.A.Ao. ... \ .A. 001 . II. . . A. .rA.I Itch.” .A.” tumor... .0. .A. 5A. 3.0.1.3... A. ’10.”.- ... . and. .o ’ AR .0. .0; 2.0.1.6 32.03030? .-AA.A~ r. 8.3.1330. ...}...r. . .A..-. I . .Q‘ .A.. AAAs‘ffoA'! '01. A..... l :15... .3." ..3. ._ .33.: A . . . il‘. . .A.. AOAA’- :1 ......1 ...! .A AAA...I I? I ...; .IIS . l .3... 03 _ 5 .’.AI 3:5“- "u~AM¢.VVA A . »u .r. 1’3 3‘1 0A..“ A» A..“.A A0 A .A. AI....AI’70.0’A A130: Avou A! AAI... I“... no t. .n . .....A’A’AOf... {A 3A.. ,. ...AIOO I32; .Av 3...”: 38.3...‘0‘ . .A. .0A. 0.. .. I: .. ..AKAAAII . .. .A.:‘A “IA. ..0 o A .. v .20”... AA...“ Rug.“ rout” A.” h A . . . - ._ .. .0. A A .9 . .00. .‘ o. .A.! . A . .-.. III. . . A ’.. l A... A A . . 0". o _.v . . . Ad I . 0 0A,. . AI AA . .A. 9 . . o . . . . A . A A . A ......AAVA .A.. Uv A. AA}. ”it“: ......Ififlo. u“. 1.3.... Atom. .3...- .IA-A...” 8.3.4.“ 7.. . .RHN.“ 16.0) ‘23 AAA...J ..z N. . .A....iIoA .«I . ....AA.M.... (IA. NANAA‘JAA. I ....n: I.u.'Au..AA..!JAA no. A" ..uA.... ”.1 ..AMAPu -...hiAQAH. ”...AAAA . A3... .A. u“ Aunt...“ no. A? .. .... z AA. AAA 0 .. Ar." ... on. . I.I NJ»! ... . . A. . A0. A0! A... o .100 0 If! {4- -.A V I. .. A) 9.: . . on ui.‘..’.0 . . 20"... Q. 0 a A, .0: . 0 I O. .0 ..ACIA 'J ST) 94 0.. ‘00 .A.. I. 00 .09 _ I 0 . A 0...... . . . euAu....c. Alcrl; ‘ .0.‘ Or: A A. I AA...’ J u q H K AA J ,o 4 A '0 A .A‘ .A. A A . . AA. A ..I. y o . 1A.. A .. . ...-A . . . AA. v A !.,.0 .A .A. ' ... .0. 1A.. \ ....A . .‘AII . . A. A A... .A. A n.- A .. . . . . .. 3.1.. 90. 4.. A 9A A 0 AAA . Z . .A....nu. .44.... M3... ..0.AA¢.o. :18 3:10.! ....SA.. a. 0.... ....A. .0 2.3.1:? A. 9.7.9.1!“ A4... A... ...”..-A .A. a“: u. . . . ALIA-4.9. "Aw-u tn... $.13. IaAIU.... n" 3'1: AA. .32.... .A.. .0...“ A EzA ....A 2..- ... ..A . .. .. .AAA..I.~Q....W .. .l. A A" ’A a: A. p r! AI 0A AAA’I. .130 A . .6539: o. .A.. .5A ... AAA .. .IA. a .A A? A .rflo H. AAIV. ......A. 3.1.. . . .Aa‘QAA..A\—r.y . . ... . . . _ . . .III'AA I . it .9557 . . . . LIBRARY Michigan State University ). ) S) K 3 II This is to certify that the g dissertation entitled 5 HIERARCHICAL PATCH DYNAMIC MODELING OF COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE SCALES presented by SOHEIL AFSHARI has been accepted towards fulfillment of the requirements for the Ph. D. degree In CIVIL ENGINEERING :1) Major Professor’ Signature ' _,,/ 2-6- 06 Date MSU is an Affirmative Action/Equal Opportunity Institution HIERARCHICAL PATCH DYNAMIC MODELING OF COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE SCALES By Soheil Afshari A DISSERTATION Submitted to Michigan State University in partial firlfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Civil and Environmental Engineering 2006 ABSTRACT HIERARCHICAL PATCH DYNAMIC MODELING OF COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE SCALES By Soheil Afshari In this research, I apply a hierarchical patch dynamics paradigm (HPDP) for solving complex groundwater problems that are difficult to solve using traditional modeling methods. The HPDP takes advantage of hierarchy theory, mimics how biological systems “divide and conquer” complexities, and decouples subsurface dynamics hierarchically. The new paradigm enables modeling complex groundwater systems in high resolution without having to solve large, ill-posed matrix systems, significantly alleviating the infamous “curse of dimensionality” and the associated modeling bottlenecks in large-scale groundwater modeling. My first application is a synthetic example patterned after a field situation concerning the regional management of water resources in the presence of site—specific constraints associated with groundwater contamination, wellhead protection, sustainable yield evaluation, and local groundwater surface water interactions. The problem is difficult to solve using traditional approaches because of the large domain size, the multiple scales of complex stresses, sources/sinks, and natural variability, and the need to characterize both regional and local dynamics. The application demonstrates that, by hierarchically decoupling the system and then dynamically integrating the sub-solutions, the HPDP enables integrated watershed-scale analysis that honors the local data in a way that is physically-based. The example also shows that the HPDP not only makes it possible to solve this large problem, but also enables solving it on the fly and interactively. It allows “zooming” in real-time incrementally from the regional scale into local hotspots of interest. The accuracy of the hierarchical solution is verified through a comparison with a single fine grid solution. Empirical guidelines for selecting hierarchical boundaries and resolutions are provided. My second application involves predicting well drawdown dynamics in large complex systems. Detailed drawdown is needed for many water management problems. Typical applications of numerical models to large problems generally require large grids that can seldom accommodate cells as small as the actual well diameter. Locally applying analytical solution within the well grid block in a numerical model represents a useful alternative, but it only applies for situations when the flow is 2D, steady, isotropic, homogeneous and free of other sources and sinks within the regional well cell. This example shows the HPDP enables accurate prediction of drawdown dynamics under general conditions in a large region. Furthermore, this accurate prediction can be achieved even on a small computer. The accuracy of the hierarchical solution is verified for a 2D well field problem in which an analytical solution exists and for a 3D unsteady problem through a single fine grid numerical solution. My final example is for a real site. The HPDP is applied to evaluate hydraulically a bioremediation delivery system for a contamination site in Schoolcrafi, MI. In particular, the HPDP is used to model the interplay of the flow dynamics at a number of scales present at the site. The example demonstrates that the HPDP makes it possible to characterize the detailed 3D delivery dynamics and evaluate its robustness in the presence of larger scale heterogeneity across the biocurtain and irrigation pumping in the region. DEDICATION I wish to thank god for giving me strength, patience, and courage for concluding my Ph.D., and delivering a thriving chapter in my life. I wish to dedicate this dissertation to: My father, for believing in me and making a dream become true My mother, for her love My family, especially my sister and brother for their support and understanding iv ACKNOWLEDGEMENTS I gratefully acknowledge the contribution of those who have made this dissertation possible. My deepest gratitude goes to Dr. Shu-Guang Li, Associate Professor, Department of Civil and Environmental Engineering, Michigan State University, my dissertation director and my advisor during completing my graduate study. I appreciate his fatherly advice, patience, and continual encouragement throughout my entire graduate program. I wish to extend my appreciation to Dr. Michael J. Dybas, Center for Microbial Ecology, Michigan State University, for serving in my dissertation committee and also giving me the opportunity to work on Schoolcrafi project, and teaching me numerous scientific, technical and practical lessons in real world remediation and treatment applications. I also wish to extend my appreciation to Dr. Manta S. Phanikumar, Assistant Professor, Department of Civil and Environmental Engineering, Michigan State University, and Dr. Roger B. Wallace, Associate Professor, Department of Civil and ' Environmental Engineering, and, for serving in my dissertation committee and for valuable comments in this study. Also, I wish to thank Dr. Qun Liu for his patience and restless help during my simulation study. My research was funded by National Science Foundation (NSF) and Department of Environmental Quality (DEQ). TABLE OF CONTENTS PAPER NO.1: AN OBJECT-ORIENTED HIERARCHICAL PATCH DYNAMICS PARADIGM (HPDP) FOR MODELING COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE-SCALES ........................................................................................ 1 ABSTRACT ........................................................................................................................ 1 INTRODUCTION .............................................................................................................. 2 A HIERARCHICAL PATCH DYNAMICS PARADIGM ................................................ 3 Hierarchical Equations ................................................................................................... 4 Patch Boundary Identification ........................................................................................ 6 Patch Grid Resolution ..................................................................................................... 7 Significance of the HPDP ............................................................................................... 8 The Achilles ’ Heel ........................................................................................................... 9 THE HPDP ENVIRONMENT ........................................................................................... 9 A “PARALLEL COMPUTING” PARADIGM AND REAL-TIME STEERING ........... 10 ILLUSTRATIVE APPLICATION ................................................................................... 13 CONCLUSIONS ............................................................................................................... 26 ACKNOWLEDGMENTS ................................................................................................ 26 REFERENCES ................................................................................................................. 27 PAPER NO. 2: A METHOD FOR PREDICTING DRAWDOWN AT THE RADIUS OF A PUMPING WELL FOR LARGE COMPLEX SYSTEMS .......................................... 28 ABSTRACT ...................................................................................................................... 28 INTRODUCTION ............................................................................................................ 29 Local Grid Refinement .................................................................................................. 31 Local Analytical Correction ......................................................................................... 3 I Local Nested Numerical Correction ............................................................................. 33 vi A DYNAMICALLY INTEGRATED HIERARCHICAL PATCH DYNAMICS APPROACH ..................................................................................................................... 34 ILLUSTRATIVE EXAMPLES ........................................................................................ 36 Single Well Example ..................................................................................................... 3 7 Well Field Example: ..................................................................................................... 44 Three-dimensional Example for Single Well: ............................................................... 51 SUMMARY ...................................................................................................................... 59 REFERENCES: ................................................................................................................ 61 PAPER NO. 3: ASSESSING THE HYDRAULIC PERFORMANCE OF A BIOREMEDIATION CURTAIN FOR PLUME G AT THE SCHOOLCRAFT SITE, MI USING AN INTEGRATED HIERARCHICAL GROUNDWATER MODEL ............... 63 ABSTRACT: ..................................................................................................................... 63 INTRODUCTION ............................................................................................................ 65 Scale Dependent Problem ............................................................................................. 68 Research Questions and Challenges ............................................................................. 69 Scope of Work ............................................................................................................... 72 OVERVIEW OF THE SCHOOLCRAFT FIELD SITE ................................................... 73 Site Characterization and Pilot Scale Tracer Test ....................................................... 73 Site characterization ................................................................................................. 73 Field tracer test ......................................................................................................... 77 CONCEPTUAL AND NUMERICAL MODELS ............................................................ 80 Tracer Scale Numerical Model ..................................................................................... 8] Calibration Process and Discussion ............................................................................. 84 BIOCURTAIN SCALE MODEL ..................................................................................... 87 Sensitivity Analysis and Discussions ............................................................................ 96 REGIONAL MULTI-SCALE MODEL ......................................................................... 127 HIERARCHICAL PATCH DYNAMIC PARADIGM (HPDP) .................................... 128 Regional Multi-scale Model Construction .................................................................. 129 Results and Discussions .............................................................................................. 133 vii SUMMARY AND CONCLUSIONS ............................................................................. 139 ACKNOWLEDGEMENT .............................................................................................. 141 REFERENCE .................................................................................................................. 142 SUMMARY OF ORIGINAL CONTRIBUTIONS ......................................... 145 viii LIST OF TABLES PAPER NO.1: AN OBJECT-ORIENTED HIERARCHICAL PATCH DYNAMICS PARADIGM (HPDP) FOR MODELING COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE-SCALES Table 1. Governing patch-dynamic equations .................................................................... 5 Table 2. Aquifer properties and boundary conditions for the example problem .............. 17 Table 3. Surface water bodies modeled in the example ................................................... 17 Table 4. Pumping wells and rates used in the example .................................................... 18 PAPER NO. 2: A METHOD FOR PREDICTIN G DRAWDOWN AT THE RADIUS OF A PUMPING WELL FOR LARGE COMPLEX SYSTEMS Table 1- Presents the numerical parameter set up for all the models in three dimensional single well example ................................................................................................... 54 Table 2- Numerical parameter setup for different submodels .......................................... 5 7 PAPER NO. 3: ASSESSING THE HYDRAULIC PERFORMANCE OF A BIOREMEDIATION CURTAIN FOR PLUME G AT THE SCHOOLCRAFT SITE, MI USING AN INTEGRATED HIERARCHICAL GROUNDWATER MODEL Table 1. Details of tracer experiments (source Graulau 2003). ...................................... 78 Table 2. Input physical and numerical parameters for tracer scale model ...................... 84 Table 3. Input physical and numerical parameters for biocurtain scale model. ............... 91 Table 4. Spatially averaged relative concentration with minimum, and maximum value obtained in the typical control volume in biocurtain scale model for 5 different realizations with 3 different hydraulic conductivity variances for bottom layer. Second column represents the case which effective porosity for bottom layer was varied spatially correlated to natural logarithmic of hydraulic conductivity. ......... 125 Table 5. Distance and pumping rate of existing high capacity wells in the vicinity of biocurtain. ............................................................................................................... 133 Table 6. Regional multi-scale input parameters for different submodels. Physical parameters for this model are same as biocurtain scale model. .............................. 133 ix LIST OF FIGURES Images in this dissertation are submitted in color PAPER NO.1: AN OBJECT-ORIENTED HIERARCHICAL PATCH DYNAMICS PARADIGM (HPDP) FOR MODELING COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE-SCALES Figure 1. Hierarchical decomposition of a groundwater system: lefi) recursive subdivision, right) hierarchical patch network ............................................................ 4 Figure 2. A typical model-patch and its interaction with the parent daughters. ................. 6 Figure 3. Modeling region for the example problem. Solid circle represents a well or a group of wells. The shaded areas (irregular polygons) denote surface water features. Note the locations of many wells are too close to be resolved and the sizes of many hydrological features are far too small to be visible in this regional map. ............... 16 Figure 4. a) Hierarchical patch network, b) Hierarchical patch network solutions. ....... 19 Figure 5. Comparison of the hierarchical and single model solutions along the y-sect A- A in patch M16, : a) regional scale, b) local scale. The location of the y-section A-A is defined in Figure 4. Note the hierarchical solution M :4. (open diamond) reproduces the single model solution (the solid line). ................................................................. 24 Figure 6. Comparison of the hierarchical model solutions along the x-sect B—B in patch area M 38, . The location of the x-sect BB is defined in Figure 4. ............................ 25 PAPER NO. 2: A METHOD FOR PREDICTING DRAWDOWN AT THE RADIUS OF A PUMPING WELL FOR LARGE COMPLEX SYSTEMS Figure 1. Presents the hierarchical-network and simulated heads after 50 days pumping for the single well example problem. Each model has a grid that is 21 rows by 21 columns in size .......................................................................................................... 39 Figure 2. Illustrates the drawdown comparison at the end of 50 days between the numerical model and analytical solution. Each model has a grid that is 21 rows by 21 columns in size ..................................................................................................... 40 Figure 3. Presents the hierarchical-network and simulated heads after 50 days pumping for the single well example problem. Each model has a grid that is 51 rows by 51 columns in size .......................................................................................................... 42 X Figure 4. Illustrates the drawdown comparison at the end of 50 days between the numerical model and analytical solution. Each model has a grid that is 51 rows by 51 columns in size ..................................................................................................... 43 Figure 5. The distribution of the pumping wells and their pumping rates for the well-field example. .................................................................................................................... 46 Figure 6. Hierarchical modeling layout which illustrates the relationship between submodels and their parent models. [In label of Submodel M 1:- a represents the generation level (grandparent), b is the parent index, and c is the kid index which is related to parent index]. ............................................................................................ 47 Figure 7 - The modeling results are compared to the analytical solution on a cross-section profile (A-A) for different submodel levels and grid resolutions (a)- large scale (b)- local scale. ................................................................................................................. 50 Figure 8- Definition sketch for three-dimensional example for single well. .................... 52 Figure 9- Illustrates the drawdown comparison at the well screen between the numerical model and Base Model solution ................................................................................ 56 Figure 10- Represents the drawdown comparison between the numerical model and Base Model solutions ......................................................................................................... 58 xi PAPER NO. 3: ASSESSING THE HYDRAULIC PERFORMANCE OF A BIOREMEDIATION CURTAIN FOR PLUME G AT THE SCHOOLCRAFT SITE, MI USING AN INTEGRATED HIERARCHICAL GROUNDWATER MODEL Figure 1. Location of study area (source Graulau-Santiago 2003) ........................... 66 Figure 2.Groundwater contour map in the VOC contaminated region[source Graulau 2003] ........................................................................................................................ 74 Figure 3. Cross section parallel to groundwater flow at vicinity of tracer test [10m distance] [source Graulau 2003]. .......................................................... 75 Figure 4. Measured hydraulic conductivity for specific depth [source Graulau 2003]. 76 Figure 5. Layout of the tracer injection system. source Graulau-Santiago 2003. .......... 79 Figure 6. Conceptual model set up for tracer scale model ................................................ 82 Figure 7. Comparison of model results and data for breakthrough curve at extraction and injection wells. .......................................................................................................... 86 Figure 8- Well spacing layout for extraction injection well system for potential full scale biocurtain design (Dybas et a1. 2003). ...................................................................... 89 Figure 9. A typical control volume [black box-solid line] for comparison purposes. Light blue zone is treatment area ........................................................................................ 90 Figure 10. Conceptual model set up for a 4 set of 10 m pumping and injection wells ..... 92 Figure 11. Steady state head contour map for 4 sets of 10 meter extraction injection well at biocurtain scale model. ......................................................................................... 93 Figure 12. Concentration map for 4 sets of 10 meter extraction injection well at the end of one day pumping. ...................................................................................................... 94 Figure 13. Breakthrough prediction at extraction and injection wells for one and four sets of 10 meter well spacing model. Notice the difference between breakthrough results due to squeeze effect which occurs in the 4 set of 10 meter model results. ............. 95 Figure 14. Comparison breakthrough curves at the extraction well [point E2 Figure 19] for different hydraulic conductivity variances and multiple realizations with homogenous [constant k] case. Notice that homogenous case result agrees with ensemble mean for all the realization. Next 4 figures presents each breakthrough curve separately. ....................................................................................................... 98 xii Figure 15. Comparison breakthrough curves at the extraction well for different realizations of hydraulic conductivity with variance of0. 15 with homogenous [constant k] case. Relative concentration value varies from %49 to %60. ............... 99 Figure 16. Comparison breakthrough curves at the extraction well for different realizations of variable effective porosity positively correlated with hydraulic conductivity variance of 0. 1 5 with homogenous [constant k] case. Relative concentration value ranges from %45 to %65. ....................................................... 100 Figure 17. Comparison breakthrough curves at the extraction well for different realizations of hydraulic conductivity with variance of 0.5 with homogenous [constant k] case. Relative concentration value ranges from %46 to %70. Notice that the breakthrough range for variance of 0.5 is bigger than the range with variance of 0.15 .......................................................................................................................... 101 Figure 18. Comparison breakthrough curves at the extraction well for different realizations of hydraulic conductivity with variance of 1 with homogenous [constant k] case. Relative concentration value ranges from %41 to %61. Notice that the breakthrough range for variance of 1 is bigger than the range with variance of 0.15 and 0.5. .................................................................................................................... 102 Figure 19. Conceptual sketch for a 4sets of 10 meter extraction injection model representing location of profile cuts and points E2 (extraction well) and W (potential weakest concentration coverage point). Point W is located along the biocurtain length in the middle of extraction and injection well and horizontally it is located in the middle of two different sets of extraction injection. It is located in the intersection of profile cut D-D and C-C. ................................................................ 104 Figure 20. Comparison breakthrough curves at point W [Figure 19] for different hydraulic conductivity variances and multiple realizations with homogenous [constant k] case. Notice that the results for homogenous case doesn’t agree with ensemble mean of all the realizations. Due to high uncertainty at point W [typical weakest point] the breakthrough range is significantly bigger than the range at extraction well [point E2 Figure 19]. ...................................................................... 105 Figure 21. Comparison breakthrough curves at point W [Figure 19] for hydraulic conductivity with variances of 0.15 with homogenous [constant k] case. Homogenous case predicts the lowest breakthrough value %9, however the range for other realizations varies from %13 to %69. ............................................................ 106 Figure 22. Comparison breakthrough curves at point W [Figure 19] for different realizations of variable effective porosity positively correlated with hydraulic conductivity variance of0. 15 with homogenous [constant k] case. Homogenous case predicts the breakthrough value %9, however the range for other realizations varies from %6 to %7 9. Notice that one of the realizations drops below the homogenous case’s result. ............................................................................................................ 107 xiii Figure 23. Comparison breakthrough curves at point W [Figure 19] for hydraulic conductivity with variances of 0.5 with homogenous [constant k] case. Homogenous case predicts the lowest breakthrough value %9, however the range for other realizations varies from %9 to %88. ....................................................................... 108 Figure 24. Comparison breakthrough curves at point W [Figure 19] for hydraulic conductivity with variances of 1 with homogenous [constant k] case. Homogenous case predicts the lowest breakthrough value %9, however the range for other realizations varies from %9 to %100. Notice that by increasing the variance the range of breakthrough results increases. ................................................................. 109 Figure 25. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable effective porosity positively correlated with variable hydraulic conductivity with variance of 0.15 and also the homogenous case [constant k] [4 sets of 10 meter extraction injection well model]. Notice that pick concentration values drops slightly comparing to the case of only k variable with variance of0.15 ....................................................................................................... 111 Figure 26. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable hydraulic conductivity with variance of 0.15 and also the homogenous case [constant k]. Notice that minimum concentration value obtains at the middle of two extraction injection set. Range of uncertainty in the results drops when we move toward the injection well ...................................................... 112 Figure 27. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable hydraulic conductivity with variance of 0.5 and also the homogenous case [constant k] [4 sets of 10 meter extraction injection well model]. Notice that range of predicted concentration value increases by increasing the variance from 0.15 to 0.5. ................................................................................. 113 Figure 28. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable hydraulic conductivity with variance of 0.5 and also the homogenous case [constant k] [4 sets of 10 meter extraction injection well model]. Notice that range of predicted values increases by increasing the variance. ................................................................................................................................. 114 Figure 29. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous case with variance of 0. 1 5 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that results for homogenous case predicts the lowest. The range of predicted realizations at pick vary from %20 to %105. ............................................................................. 117 Figure 30. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous case with variance of 0.5 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that xiv due to high uncertainty the range of predicted values are very large. The range of predicted realizations at pick varies from %20 to %122. ....................................... 118 Figure 31. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous case with variance of l and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that homogenous case predicts the lowest concentration across the profile and range of results for other realizations vary fi'om %20 to %128. ........................................... 1 19 Figure 32. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous porosity and conductivity with variance of 0.15 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that one of the realization predicts value lower than homogenous case and the range of prediction varied from %9 to %118. .................................... 120 Figure 33. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous case with variance of 0.15 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. Notice that there is a good agreement with homogenous case prediction and other realizations. ............................................................................................................. 121 Figure 34. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous porosity positively correlated with hydraulic conductivity with variance of 0.15 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. .......................... 122 Figure 35. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous case with variance of 0.5 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. Notice that by increasing the variance the range of prediction at the pick area increases. ................................................................................................................. 123 Figure 36. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous case with variance of 1 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. Notice that by increasing the variance the range of prediction at the pick area increases. Results for case of variance 1 has higher predictions and higher ranges comparing to cases of 0.15 and 0.5. ........................................................................ 124 Figure 37. Concentration layout and profile cut for a) 24 hour period during the recirculation system, b) 10, and c) 14 days after shutting down the wells. ........... 126 Figure 38. Regional multiscale modeling domain and existing high capacity well [circled] at vicinity of biocurtain. .......................................................................................... 130 XV Figure 39. Hierarchical modeling layout for regional multi scale model which illustrates the relationship between submodels. ...................................................................... 132 Figure 40. Breakthrough comparison for an extraction well from regional multi scale model during the one day recirculation for cases of with and without regional inigation. Notice that results for two cases are very close. ................................... 135 Figure 41. Profile cut D-D from regional multi scale model presents the concentration distribution along the biocurtain for cases of with and with our regional pumping. ................................................................................................................................. 136 Figure 42. Head comparison across the profile A-A Figure 39 for four different scenarios a) natural gradient, no pumping b) presence of irrigation and no recirculation 0) presence of irrigation and recirculation d) presence of recirculation and no inigation. ................................................................................................................................. 138 xvi PAPER NO.1: AN OBJECT-ORIENTED HIERARCHICAL PATCH DYNAMICS PARADIGM (HPDP) FOR MODELING COMPLEX GROUNDWATER SYSTEMS ACROSS MULTIPLE-SCALESit Shu-Guang Li, Qun Liu, Soheil Afshari Department of Civil and Environmental Engineering Michigan State University, East Lansing, MI 48824 Accepted for publication in Environmental Modeling and Software ABSTRACT In this paper, we present an obj ect-oriented hierarchical patch dynamics paradigm (HPDP) for modeling complex groundwater systems across multiple scales--one that has the potential to significantly alleviate the infamous “curse-of-dimensionality” and the associated computational bottlenecks. The HPDP takes advantage of hierarchy theory, divides and conquers complexities, and decouples scale-dependent dynamics hierarchically. The object-oriented HPDP represents a generalization of the “telescopic- mesh-refinement” (TMR) techniques, providing dynamic model coupling, visual interactive steering, and freeing the modelers fi'om the impractical task of having to interact offline and iteratively with potentially large numbers of modeling patches. The HPDP provides a valuable tool for understanding scale-dependent processes and for ‘ Latest version of paper is available at www.cgrmgreduhlishug 1 practical groundwater investigations. We illustrate the effectiveness of the HPDP and its potential for multi-scale groundwater modeling using an example. Software Availability Software name: Interactive Groundwater (IGW) Developers: Shu-Guang Li, Qun Liu Hardware/software required: IBM-compatible PC/Windows- NT4/95/98ME/2000/)(P/Server 2003 Programming language: Visual-Basic 6; Visual-FORTRAN 5; Borland C++ Builder 6; Borland Delphi 5 Program size: 150MB Availability: downloadable with documentation at http://www.egr.msu.edu/igw. License: free trial version 3.x INTRODUCTION Despite the exponential growth of computational capability over the last two decades— one that has allowed computational science to become a powerful tool for scientific discovery - our ability to model large groundwater systems and accurately represent complex processes is still limited because of the following computational challenges: 0 The machine bottleneck. Modeling computation increases exponentially with the domain size and the level of details simulated and quickly becomes prohibitively expensive for large problems. This is especially the case for multi-scale modeling, coupled processes modeling, inverse modeling and parameter estimation, uncertainty analysis, hypothesis testing, and data sufficiency evaluation. 0 The algorithmic bottleneck. Modeling based on a single numerical representation of a complex system faces an algorithmic bottleneck. The high-dimensionality, especially when combined with distorted-grids representing heterogeneity, a multitude-of-scales, and complex structure, translates into ill-posed matrices and causes a host of numerical problems, including solution failure. A HIERARCHICAL PATCH DYNAMICS PARADIGM In this paper, we address systematically these difficulties or the so-called “curse-of- dimensionality” in numerical modelling. In particular, we present a hierarchical patch dynamics paradigm (HPDP) for modelling complex systems across multiple-scales. The basic idea is to move away from developing a single monolithic system that would solve a large and ill-conditioned problem, and to employ a hierarchy of simpler sub-systems that cooperatively solve the problem. In such hierarchical systems, each sub-system is simpler to solve than the monolithic system. The integrated system based on the whole population can generalize better than any single subsystem in the hierarchy. Specifically, we propose to represent a complex system hierarchically, breaking it “vertically” into many levels and “horizontally” into many patches (Figured-left). The dynamics in different patches are related efficiently to each other through a hierarchical patch network (Figurel-right). The patches at higher-levels are used to simulate larger- scale dynamics whereas nested-patches at lower-levels are used to simulate smaller-scale processes. The upper-level exerts constraints (e. g., as boundary conditions) to the lower- level. The conceptual basis for the HPDP emerges from hierarchy theory (Simon, 1973) and a diversity of studies in various disciplines, including management science, biology, ecology, and system science. The hierarchy theory has been significantly expanded in the context of evolutionary biology and ecology (Wu and David, 2002). level 1-1 p. ‘\. Coarse ‘\ . \_ r" fl level I / \\ ., A Wjilhhmc Finest scale patches Figure 1. Hierarchical decomposition of a groundwater system: left) recursive subdivision, right) hierarchical patch network Hierarchical Equations Mathematically, the HPDP involves solving a set of hierarchically-formulated groundwater equations (Table 1). These equations describe the space-time distribution of head h, velocity u,-, and concentration c for a typical patch at level I defined on Q’ , and its interaction with its “parent” patch on Q“ at level [-1, and various “daughter” patches (Qf', k=1,2,.., ) at level 1+] (Figure 2). Equations (1-2), (1-3) and (1-5) and (1-6) provide patch boundary and initial conditions, reflecting the coupling with the parent at level 1-1. The symbols ne, S,, R, K,-,- and ng are effective porosity, specific storage, retardation factor, conductivity-tensor, and dispersion- tensor; 8 and r are the water and solute source/sink term. We assume that groundwater is 4 incompressible, homogeneous, and isothermal and the medium saturated and non- deformable. The equations in Table 1 can be solved sequentially for l=0, 1,2, , based on a scale- dependent discretization through a systematic process called “downscaling”. This process begins with modeling at the highest-level on a coarsest grid/timestep and proceeds downward to finer-scales with increasing space-time resolutions. Successive solution of Eq.1 down the hierarchy provides an approximation of the system dynamics and a mechanism to pass field information from high to low-levels. I 1 I ah ah ah I a 1 1 1 1 I 1-1 - —=— K..— ; .=-—,1<..— Flow Eqs. SS at 6x.[ 1} 6x . ]+s u! "I I] 6x. I j e .1 1 5”] 14‘3"]-1 1 1—1 Flow-BC K..—=K.. ——or h =h 1-2 U 6x . U 6x . J J Flow-1C h1(x,0) = h’ ‘1(x,0) 1-3 16c] aufc’ a D, ad 1 .. R -— —-—=— .. TWW“ Eq’ 5: + 6x. ax. y ax. H M J I J l l—l ac dc ufclni —D;. a ni =uf_1cl_1nl. —Di:._1—a—ni or _ x . x . - Transport BC 1 J 15 61:61—1 Transport-1C c;(x,0)=c;—l(x,0) 1-6 Table 1. Governing patch-dynamic equations Coarser [-1 1+1 Finer k=1, 2, 3,... Figure 2. A typical model-patch and its interaction with the parent daughters. Patch Boundary Identification The HPDP represents a generalization of the frequently-used telescopic mesh refinement (TMR) techniques (e. g, Mehl and Hill, 2002) and the marriage of hierarchical modeling and patch-dynamics (Wu and David, 2002) to specifically address the infamous “curse- of-dimensionality” and the associated computational bottlenecks in multi-scale modeling. The way to specify patch boundaries are essentially the same as that in a basic TMR. In other words, patch boundaries are identified empirically based on the parent solution behavior, system characteristics, and data densities. Modeling patches are created in localized areas exhibiting much more rapidly-varying dynamics. The patch boundaries are selected such that they are located outside the “hot-spots” or in areas where the parent head varies relatively slowly and is deemed accurate. In many cases, the patch definition may be significantly facilitated by the fact that groundwater usage, data distributions, impacts from industrial activities, and the areas of concerns are often naturally-clustered spatially at different scales. For example, a groundwater basin may include multiple regions of much higher groundwater uses. Each of these regions may consist of multiple wellfields, which may in turn include many well- clusters made of individual wells. This naturally-nested clusteredness is ideally-suited for hierarchical modeling. In general, identifying patches for local grid refinement may have to be iterative. A patch solution is deemed accurate if it becomes insensitive to changes in the patch boundary location. One may also potentially avoid the need to iterate or minimize the number of trial runs needed if modeling patches are defined “conservatively” or made larger than needed. Patch Grid Resolution To ensure that any model in the patch hierarchy is well-posed and can be easily solved we always limit the patch grids to relatively small sizes that can be easily handled using the computer available (e. g. 50 by 50 by 3 for a geological layer). Highly variable local dynamics that cannot be resolved are left to additional nested-submodels. This nested modeling process is applied recursively for as many times as desired until one obtains needed details. Given the patch grid resolution, the number of nested patches needed generally increases with the problem size and complexity. Given the problem size and complexity, the nrunber of nested patches decreases with the increasing patch grid resolution. Significance of the HPDP The generalized HPDP makes it possible to: 0 model approximately large systems in high resolution without having to solve large matrix systems and thus significantly alleviate the machine bottlenecks. 0 model heterogeneous dynamics incrementally, thus obviating the need to use a highly-distorted grid representing multiple-scales of variability, and significantly alleviating the algorithmic bottleneck. As the nested grid is refined, the geometries of the sources and sinks and parameter heterogeneity can be represented more and more accurately. Discretized lakes and rivers approach their actual shape and discretized wells become asymptotically “point” sources/sinks at their actual locations. Detailed local measurement clusters can be also resolved more faithfully. We feel that the HPDP is essential to understanding large, complex systems. It is hardly justifiable theoretically and overwhelmingly difficult technically to translate information directly between two distant-levels, when ignoring relevant intervening-levels. Although it is possible to scale down from regional processes at basin-scale to detailed dynamics at site-scale, successful approaches most likely have to be hierarchical. The Achilles’ Heel Before one can reap the potential benefits of the HPDP we have to cross another critical hurdle. The following are inevitable questions facing practical implementation: 0 How can one manage the large numbers of modeling patches? 0 How can one tolerate the potential enormous numbers of files and book-keeping necessary? 0 How can one know a priori where to define the patch boundaries? 0 How can one effectively route the information across the hierarchical network? Indeed, a general and flexible implementation of the HPDP under the traditional, batch- based modeling flarnework can be difficult, if not impractical. To transform the powerful concept into a tangible impact, we need a new way in modeling — one that allows accommodating flexibly-and-seamlessly the complex interactions of large numbers of disjointed-models, fleeing the modeler flom having to interact offline and perhaps iteratively with individual patches and large numbers of data-streams. THE HPDP ENVIRONMENT We have recently developed such an approach and a software to support the HPDP. The new environment enables flexible and interactive patch-creations and simulation of evolutionary patch-dynamics. It provides dynamic integration of patch-models at each time-step during simulation and makes the modeler feel as if they were using a single model that can provide efficiently hi gh-resolution dynamics. Results flom the parent- models (e.g., boundary-conditions), and any changes to the parent-models propagate automatically to all children, without the need for offline post-processing. The HPDP becomes a process of successively and visually zooming into more details. A modeler can insert interactively and recursively a hierarchy of patches into a parent-model in order to provide greater detail where it is required. Modeling under the HPDP continually displays results that have been intelligently processed and integrated. A “PARALLEL COMPUTING” PARADIGM AND REAL-TIME STEERING We achieve these interactive hierarchical modeling capabilities by adopting a “parallel computing” paradigm and implement it in a single, seamless and obj ect-oriented program. The term “parallel computing” here does not mean modeling on multiple-processors but, rather, a new way of structuring computation. The basic concept is simple. Instead of treating flow and transport separately, we model them concurrently. Instead of treating the various scales of patch-dynamic modeling (regional, subregional, local, site, hot spots) as different phases in a long sequential batch-process, we couple the multi-scaled processes and model them simultaneously. Instead of relegating the graphical result presentation and analysis for each patch of interest to the “post-processing” phase, at the end of a time-consuming sequence of many disjointed-steps, we incorporate them into a single program along with the simulation, to 10 permit the interpretation of results as soon as they become available, at the end of each time-step. To achieve this, the software utilizes a number of dynamic-linking-libraries (DLLs) and dynamically-embedded software components, including a single-step flow modeling DLL, a single-step transport modeling DLL, a sparse-matrix solver DLL, a Geostatistics DLL, a 2D-graphics component, a GIS mapping component, and a 3D-visualization component. At a particular time-level tn and given the head h,.-1 and concentration c,.-1 at time-level t,,-1, the software performs the following tasks: 1. Call the flow DLL to compute head h,, and velocity u,, given h,,-, . If particles are released, compute particle coordinates at tn given those at t,,.,. If a concentration plume(s) is introduced, call the transport DLL to compute 0,, given ca. 1, If a patch(es) is introduced, perform single-step subscale flow modeling with the boundary/initial conditions extracted flom the parent flow-models. If particles/contarninant(s) are also introduced, perform single-step subscale particle- tracking/transport-modeling, with boundary/initial conditions flom the parent transport-models. Call the 2D graphics component to visualize overlays of model output/input distributions (head, velocity, particles, concentration, aquifer-properties) for a selected layer or cross—sections. If requested, perform and visualize water/solute budget over interactively-defined zones/volumes in selected patch-models 11 7. If requested, call the 3D visualization component to visualize aquifer flamework, head, velocity, plumes, particles, scattered-observations, and monitoring-network in selected patch-models. The “parallel computing” strategy makes the HPDP practically implementable since it eliminates the strangling I/O, offline-processing/ file organizations and the disconnect between potentially large numbers of patch-models, analyses, and visualizations. The visual interactive implementation enables computational steering since it provides flexibility and allows tracking user-actions. Therefore, the modeler is always "in-the- loop". At anytime, the modeler can interrupt/pause the program to assess the solution, and based on this assessment, to insert patches in selected areas and layers of interest. Since the modeling-hierarchies are dynamically-coupled, the modeler can easily perform sensitivity-analysis, making changes to patch boundaries, resolutions, conceptual representations, solvers, and assessing iteratively the influence of these changes on the solutions. The obj ect-oriented “parallel” implementation allows the investigator to interactively steer the hierarchical computation, to control flexibly the program execution-sequence, to guide the evolution of the aquifer- dynamics, to control the visual data representation during processing, and to dynamically modify the hierarchical computational process during its execution. 12 ILLUSTRATIVE APPLICATION We illustrate in this section how the generalized, dynamically-integrated HPDP significantly enhances our ability to model a large complex region across multiple scales with improved accuracy, efficiency, and flexibility. Our example is patterned after a real-field situation concerning integrated, regional management of water resources in the presence of scattered, localized, and site-specific constraints associated with groundwater contamination, wellhead protection, sustainable yield evaluation, local groundwater surface water interactions, and monitoring network design. The problem is difficult to solve using a traditional approach because of the large domain size, the multiple scales of variability, and the need to characterize both regional and scattered local dynamics. The issues involved are characteristic of many current water resources and environmental problems in the nation and around the globe. Our example considers a region that encompasses an area of 10,000 square kilometers (figure 3), including many wells of very different pumping capacities (table 4) and a network of sloughs, wetlands, wildlife habitats, and lakes of vastly varying sizes (Table 3). Table 2 presents key parameters describing the regional aquifer system features. The wellfields which provide drinking water for major cities nearby, lies in the proximity of several industrial contamination sites, two of which have caused significant contamination of the underlying groundwater. The local communities are concerned how the various contamination sites may impact the drinking water supply and how the city 13 pumping may complicate the on-going characterization and remediation activities at the various contamination sites and potentially draw more contamination into the aquifers flom the overlying pollution sources. The local communities are also concerned about potential water rights and conflict issues and how city pumping may in general stress the ecosystem, potentially disrupting the wetland and groundwater links and drying up wildlife habitats. To enable addressing these issues in a systematic and integrated fashion, we apply here the new, dynamically integrated HPDP that allows modeling the aerially expansive area while at the same time easily and fleely zooming interactively into local subareas or patches across the entire regional aquifer system . To demonstrate the advantage of the approach and verify its accuracy, we also utilize the single model approach to simulate the entire region in highest resolution possible on our 2GHz/2GB Pentium 4 computer and compare the single model solution with the HPDP-based predictions. Figure 4 presents the final, complete hierarchical patch network solutions to the example problem. The HPDP-based solution allows characterizing the regional groundwater system across multiple scales in response to both regional and local stresses and heterogeneity. The hydraulic conductivity is regionally defined by 130 scattered values and a number of isolated data clusters, low permeability lenses and preferential channels of different sizes. The aquifer stresses are characterized by large numbers of sources/sinks of very different scales and strengths. The dynamically coupled patch solutions and especially their sensitivity to different integrated management options provide valuable information for integrated analysis and allow addressing both coupled 14 regional issues and local problems and site specific constraints related to, e.g., contamination, remediation, water conflicts, pollution control, groundwater and surface water interaction, and wellhead protection. The interactive, flexible HPDP allows capturing regional and subregional hydrogeological controls and at the same time resolving detailed local data, aquifer features, local heterogeneity, preferential channels, plume structures, and even drawdown dynamics around individual low capacity wells without having to solve large, ill-conditioned matrix systems or bottlenecked by the limited computer system capacity. Under the HPDP, we model, visualize, and analyze the region incrementally and hierarchically. Figure 4b provides a graphical representation of the hierarchical patch dynamic network for the overall groundwater system. We begin with modeling the heterogeneous system using a uniform coarse grid and then refine the dynamics in areas where we believe the predictions are inaccurate and where we are interested in knowing more. We create modeling patches, and patches within a patch recursively until we are able to obtain sufficient details in the areas of critical interest (e. g., the local contamination areas, wells, wellfields, and other local “hot spots”). The modeling patches are interactively created based on a visual examination of the parent solution and the patch boundaries selected in areas where the predicted parent dynamics are deemed accurate or properly resolved. 15 . (A . 1A . 0 Lin. 60000 'r‘np 0 L34; Wellficld B i L7, .‘p. . W3 : d“ 1: .,L1 *‘= ° . L3 -3 40000(m) ~ ’ o . Wellfield A L5 0 0 60000011) 80000 1 l Figure 3. Modeling region for the example problem. Solid circle represents a well or a group of wells. The shaded areas (irregular polygons) denote surface water features. Note the locations of many wells are too close to be resolved and the sizes of many hydrological features are far too small to be visible in this regional map. 16 Aquifer type Unconfmed Domain size 100000 m by 100000 m Head on north boundary 10 m Head on south boundary -90 m Aquifer bottom elevation Variable —50 m to —150m Recharge 1.392 x 10'5m/day Pumpingwells (30 wells) 1500 to 10 gpm Number of surface water bodies simulated (lakes, streams, wetlands) 20 (see table 3 for more detailed description) Variable, more data in plume areas, regionally Hydraulic conductivity described by 130 scattered points, locally by zone- based heterognrw. Average K =15m/day Effective porosity 0.3 Dispersivity 0 concentration Contamination source area / Plume B: 560 m2 / 100 ppm Plume A: 300 m2/ 100 ppm Transport time step 50 days Maximal plume simulation time 3 years Table 2. Aquifer properties and boundary conditions for the example problem Surface water bodies Elevation (m) Leakance Area (m2) Comments see Figures 3 and 4) (1/day) L1 Lake -26 l 2.221e7 Head dependent flux L2 Lake -48 1 9.541e7 Head dependent flux L3 Lake -32 0.1 1.878e7 Head dependent flux L4 Lake -17 0.05 1.3 39e6 Head dependent flux L5 Lake -31 0.1 1.893e7 Head dependent flux L6 Lake -90 NA 6.249e8 Const head L7 Lake -10 to 20 0.1 7.597e7 Head dependent flux L8 Lake -21.5 0.05 3.23e5 Head dependent flux L9 Lake -22 0.01 4.9e4 Head dependent flux L10 Lake -22.6 0.01 5.265e4 Head dependent flux L11 Lake -24 0.05 5.688e5 Head dependent flux L12 Lake -17.5 0.05 2.297e6 Head dependent flux L13 Lake -18 0.05 4.296e6 Head dependent flux 81 Stream -90 to ~30 l 8.071e7 Head dependent flux 82 Stream -23.4 to 26 0.01 1.765e5 Head dependent flux S3 Stream -23.2 to —22.6 0.01 5.835e4 Head dependent flux S4 Stream -24 to -18 0.05 2.678e6 Head dependent flux W1 Wetland -18 0.01 5.501e5 Drain W2 Wetland -l9.5 0.01 3.51e3 Drain W3 Wetland -27 0.01 1.049e6 Drain Table 3. Surface water bodies modeled in the example 17 Well X(m) x 10‘ Y(m) x 10‘ Rate ID (gpm) 1 4.8000 3.4000 1000 2 2.4000 5.4000 1000 3 5.428125 5.46875 10 4 5.4250 5.46875 10 5 5.6000 5.6000 100 6 5.6375 5.6625 100 7 6.0500 6.4000 10 8 6.3500 6.0000 100 9 6.0000 6.4000 100 10 5.3500 6.3000 100 1 1 5.2000 6.6000 1500 12 5.4625 5.4625 100 13 5.4750 5.4375 200 14 5.0000 3.6000 1000 15 6.0000 6.4500 460 16 5.4000 6.3000 800 17 6.4000 6.0000 100 18 5.6125 5.6000 10 19 5.465625 5.465625 10 20 7.6000 2.8000 1000 21 2.8000 5.2000 1000 22 4.4000 2.6000 1000 23 7.4000 3.0000 100 24 3.6000 3.2000 40 25 3.8000 5.8000 100 26 5.4000 5.4375 20 27 5.48125 5.43125 20 28 5.528125 5.39375 10 29 5.5375 5.4000 40 30 5.54375 5.390625 10 Table 4. Pumping wells and rates used in the example 18 b) patch network solution 2 a) Hierarchical patch M network Grid and x-Domain size M;; Ax = 2000m;Lx =100000m M2. Ax =SOOM'L =50000m ’1' ’ X 3 M..; Ax =250m',L =30000m U X M..; Ax =125m;Lx =15000m M..; Ax =62.5m;Lx =6000m M..; Ax = l5.625m',LX = 3000m 7 M..; Ax =3.9m;L = 200m ll X M8'Ax 09 'L 100 0., — .m, X- m Figure4. a) Hierarchical patch network, b) R Hierarchical patch network solutions. Plume A WHPA 3 yrs. Plume B 19 For this illustration, we have created a total of 8 different hierarchical levels and 11 patch models zooming into 3 different areas of critical interest (including 2 local contamination sites, 1 wellhead protection site). We have purposely limited the maximal grid size of each modeling patch to about 100 by 100 in order to drive home the point that we can solve a large problem even on a relatively “small” computer (one that is only powerful enough to solve small 0(100 by 100) matrix systems), if we approach the large problem hierarchically and take advantage of how the patch dynamics are correlated. In general, the number of levels and patches needed to achieve a given resolution in hierarchical modeling depends on how large and complex the problem is, how powerful the computer is, and how much data is available. Obviously, the number of hierarchical levels required decreases as the problem sizes decrease and the computing power available increases. We can “zoom” into the details needed in one “shot” if we were to have a computer of infinite speed and memory and if numerical solvability were not an issue however heterogeneous, anisotropic, singular, or complex the aquifer system in question is. Figure 4-M111 shows the coarsest solution at the top of the hierarchical model tree. The symbol M ,1 represents a patch model at level a whose parent is b and is the cth children of b — a represents the generation level (grandparent), b is the parent index, and c the child index. This highest level solution provides the big picture of the general dynamics and overall water budget balances. Figure 4-M211 presents a patch model that provides more details in describing the subregional dynamics and the interactions between the 20 wellfields. Figures 4-M311, M411, and M5 11 present additional hierarchically nested patch models zooming incrementally into an environmentally sensitive region. These modeling patches allow capturing more detailed surface and groundwater interactions and the interactions among the wells and wellfields. Figures 4-M611 , M7ij , and M8”- present successions of patch models further zooming incrementally into two contamination hot spots and a wellhead protection delineation project that play critical roles in the overall scheme of integrated regional groundwater management. The initial coarse model M111 at the highest level covers a large area with a grid resolution of 2000 m. The modeling patches at the finest level (level 8) cover much smaller areas with a resolution being three orders of magnitude higher at approximately 1 m. The grid resolution used at level 4 corresponds to the highest resolution we can achieve using the single model with a uniform grid of 125m on our Pentium 4 computer. At this resolution, it takes 36 hours for the single model approach to solve just the groundwater flow equation. By contrast, the hierarchical patch approach allows providing flow and transport solutions at multiple resolutions and integrated visualizations, mapping on the fly. Of course, the most important advantage of the object-oriented hierarchical modeling approach lies in that it allows at any time zooming into new patch areas of interest or further into more details within existing patches. Figures 4-M51 1, M6“, M7,], and M8,]- show four additional levels of local patch solutions at a grid resolution of respectively 62.5m, 15.7m, 3.9m, and 0.9 m. These successive patch models provide scale dependent 21 representation of the stresses, heterogeneity, site features, locations of the wells, and increasingly more detailed descriptions of local dynamics of interest within the regional system. Figures 5a and 5b present a detailed comparison of the predicted water level along the A- A cut (see Figure 4 for the cut location) using the single model (solid line, grid resolution = 125m) with the corresponding hierarchical patch solution (flom M411 of same resolution, represented as open diamonds). The results clearly show that, despite the hierarchical decoupling, the local patch model W11 provides virtually identical solution to that flom the single fine model. The finer patch models M511 and M6“ provide more detailed, deeper drawdown cone of depression than that can be achieved using the single fine model. The results also show that the coarser patch model M311 (dotted line) is unable to provide the necessary details in the area around the wells but is sufficient to provide the boundary condition for M4” away flom the wells or wellfields. The coarsest patch model M111 (dash line) totally misses the local dynamics although, as illustrated in figure Sa, it is adequate at the regional or subregional scale. Figure 6 presents a comparison of the solute concentration profile along the B-B cut (see Figure 4 for the cut location) obtained flom the three fine patch models M611, M713, and M831. The solute plume is conservative and originates flom a continuous source with a steady concentration of 100 ppm. At the finest level, the patch model M831 is able to characterize the detailed plume structure and reproduce the maximal source concentration even with a relatively simplistic transport solver (the modified method of characteristics with linear backward interpolation). The coarser patch models M713 and especially M6” 22 suffer flom significant numerical dispersion and underestimate the maximum concentration. The coarsest patch models at levels 1 through 5 cannot even pick up the plume because of the large grid size. 23 a) 20- Single fine solution; I L A l l A A J A Ax =125m E40 j " ‘ ‘ ' M111; Ax =2000m '3; l """"""" M121; Ax =500m :40 ............................ M131; Ax=250m <> M141; Ax =125m '60 i x M5 Ax=62.5m 11; - ..... M161; Ax =15.625m 100000 58000 Y(m) Figure 5. Comparison of the hierarchical and single model solutions along the y-sect A- A in patch M f; : a) regional scale, b) local scale. The location of the y-section A-A is defined in Figure 4. Note the hierarchical solution M 1‘, (open diamond) reproduces the single model solution (the solid line). 24 100 90 80 . 70 9 p ‘3 ... . ,1 E 60 - / l " 1 1‘ 6 . = E50 1 ‘ M,,,Ax 0.9m ] I l — -—o— — M13;Ax=3.9m 5 4° i l \ u : I \ 1:1 M g, ;Ax=15.625m 30 “l I ‘ 20 l l l : I 1) ,\ 4 \ 10 1 I \ I, \\ ~ [:1 2 Cl .- -. 0 — / 1 1 1 1 0’ 55100 55125 55150 55175 X(m) Figure 6. Comparison of the hierarchical model solutions along the x-sect B-B in patch area M f, . The location of the x-sect B-B is defined in Figure 4. Finally, we want to stress that our numerical experience shows that the HPDP is consistent with how human brain works and processes information. Although the total number of submodeling patches can be quite large daunting for large problems, especially on a small computer, they are created in an incremental, flexible, and intuitive fashion and in a way consistent with how we naturally assimilate information. In fact, even if the computer we have were infinitely fast, none of our five senses has sufficient bandwidth to pennit conveyance and interpretation of the immense amount of data 25 obtained flom a single large-scale simulation. This simple physiological constraint still requires that we process, visualize, analyze, and assimilate the results hierarchically and in patches — beginning flom large-scale, general dynamics to ultimately finest-scale, critical details. CONCLUSIONS We have presented an interactive, dynamically-integrated HPDP for modeling complex groundwater systems and illustrated it with an example. The preliminary results show the effectiveness of the HPDP and its potential to significantly alleviate the computational bottlenecks in large-scale groundwater modeling. We are currently conducting additional testing on the HPDP, especially with respect to its ability to model more complex and heterogeneous systems. The results will be reported in future publications. ACKNOWLEDGMENTS This research is funded by the NSF under grant CISE-0430987. 26 REFERENCES Li, 8.6 and Q. Liu. (2006). “Interactive Ground Water (IGW)”, short communication, Environmental Modelling and Software. 21(3), 417-418. Mehl, S., and M. C. Hill. (2002). “Development and evaluation of a local grid refinement method for block-centered finite—difference groundwater models using shared nodes”. Advances in Water Resources 25 (5):497-511. Simon, HA. (1973). “The organization of complex systems. In: Pattee, Om H.H. (Ed.), Hierarchy Theory: The Challenge of Complex Systems”. G. Braziller, New York, pp.1e27. Wu, J ., David, J. (2002). A spatially-explicit hierarchical approach to modeling complex ecological systems: theory and applications. Ecological Modeling 153, 7e26. 27 PAPER NO. 2: A METHOD FOR PREDICTING DRAWDOWN AT THE RADIUS OF A PUMPING WELL FOR LARGE COMPLEX SYSTEMS. Soheil Afshari] , Richard Mandlez, Qun Liu', and Shu-Guang Lil 1Department of Civil and Environmental Engineering Michigan State University, East Lansing, MI 48824 2Michigan Department of Environmental Quality, Lansing, MI 48909 ABSTRACT One of the challenges in groundwater modeling is the prediction of hydraulic head in close proximity to a pumping well using a regional-scale model. Typical applications of numerical models to field-scale problems generally require large grids that can seldom accommodate cells as small as the actual well diameter. Several methods have been used to simulate a more accurate head at the well scale. The three primary methods are: 1) local grid refinement, 2) local analytical correction, and 3) local nested numerical correction. However, all these methods have limitations, particularly for applications that involve the development of numerical models for large-scale hydrogeologic systems with many pumping wells. In this paper, we present a method utilizing the hierarchical patch dynamics paradigm (HPDP) (Li et al. in press) to predict head at the well scale. The HPDP enables converting a large, complex problem into a network of hierarchically nested and dynamically coupled patch models that can be easily solved. The performance of the method is verified against the analytical solution for a single well, against a ° To be submitted to Ground Water, latest version of paper is available at www.cgr.rnsu/~lishug 28 superposition of analytical solutions for a well field, and against a numerical solution in a 3D heterogeneous system involving a well. The results show that the HPDP-based method is capable of matching an exact solution of well drawdown and providing an accurate and efficient representation of head and groundwater velocities in a well field in large-scale hydrogeologic systems. INTRODUCTION In groundwater modeling, it is usually impractical to employ grids that are comparable in size and dimension to a pumping well. However, predicting head, or drawdown, in close proximity to a pumping well (well scale) may be important for many groundwater flow and solute transport modeling applications, especially in complex environments. In groundwater flow modeling, long term drawdown prediction at the well is critical for proper well design and for evaluation of groundwater management, pollution control, remediation, and sustainability strategies. In solute transport modeling, the estimated velocity field is derived flom head values at the model nodes; this is needed to predict the advection component of the transport of solutes. In finite-difference models, many discrete cells, that typically have relatively large spatial dimensions, are used to represent the aquifer system. A point source or sink of water is injected, or extracted, over the volume of aquifer represented by the cell that contains the point source or sink. However, the diameter of the well is typically much smaller than the dimensions of the cell. F ield-scale problems generally cover large geographic areas 29 that require grids having cells with large spatial dimensions and can seldom accommodate cells as small as the actual well diameter. The resulting simulated heads are generally not a good approximation of heads or hydraulic gradients in close proximity to the pumping well; numerous small cells are needed to simulate the relatively steep hydraulic gradients near a point source or sink accurately. However, regional model- derived heads may be correct at nodes located away flom the point source or sink (Anderson and Woessner. 1992). Several methods have been adopted in an attempt to improve the accuracy of simulated heads at the well scale. Three commonly used methods in finite-difference models are: 1) local grid refinement (LGR), 2) local analytical correction (Prickett, 1967; Peaceman, 1978; Pritchett and Garg, 1980), and 3) local nested numerical correction (e.g., Mehl and Hill, 2002; Ward et al., 1987). However, all these methods may have limitations, particularly for applications that involve the development of numerical models for large- scale, complex hydrogeologic systems with many pumping wells. The objective of this paper is to address the limitations in the traditional methods and present an accurate and efficient method for modeling complex well dynamics utilizing the hierarchical patch dynamics paradigm (HPDP) (Li et al. in press). We begin with a very brief review of the traditional methods for calculating hydraulic heads at the well scale and point out their practical implementation difficulties. 30 Local Grid Refinement Local grid refinement represents the most commonly used method to predict detailed well dynamics in a numerical model. The relatively large finite-difference grid cell flom a model is subdivided or refined with progressively smaller spatial dimensions. This results in a more accurate estimation of hydraulic head or drawdown at the well scale. In other words, there will be less approximation error involved with averaging the predicted head for that particular well node. In using local refinement for simple or small-scale problems, the solution is obtained quickly, and consistency between the regional and local area around the wells is maintained. However, for large-scale, regional groundwater models where there may be a significant increase in the number of nodes (e. g., millions rather than thousands), the cost of computation increases exponentially and the process can become prohibitively expensive. Local Analytical Correction An alternative approach to model detailed well dynamics in a regional model is to use local analytical correction within the cell containing the well. The model calculated head at the well node can be thought to represent the head at an effective distance (re ) flom the well node. An estimation of the head in the well can be obtained flom formulas based on the following steady-state Thiem equation (Thiem 1906; Anderson and Woessner, 1992» 31 Q r h =h. .__WL1n _e_ (1) W la] 271T ’IV where hw is the head in the well; hi . is the head computed by the finite-difference 3 model for the well node (ij); QW T is the total pumping or injection rate flom the well; T is the transmissivity in the well cell; re is the effective radial distance measured flom the node at which head is equal to hi . ; and ’iv is the radius of the actual well. 9 Prickett (1967), Peaceman (1978, 1983) and Trescott et al. (1976) provided an equation to approximate re based on different model grid sizes. Trescott et al. (1976) examined the performance of the corrected drawdown based on Eq. (1) against the analytical solution. Results revealed a reasonable improvement for predicting the drawdown in the pumping well. Trescott et al. (1976) also stressed that the local analytical correction based on (1) is applicable only if the following assumptions are satisfied: (1) flow to the well is within a square finite-difference cell (well block) and can be described by a steady-state equation with no source term except for the well discharge, (2) the aquifer is isotropic and homogenous in the well block, (3) only one well, located at the cell center, is in the well block, and (4) the well firlly penetrates the aquifer. 32 Local Nested Numerical Correction A more general approach for modeling detailed local well dynamics is through “local numerical correction” or nested grid modeling. In this approach, grid-dependent information flom the regional model is used to construct a separate model with finer grid spacing to obtain more information around the area of interest (Ward et al. 1987). The model with the finer grid is called a “submodel” or “local model.” In order to obtain detailed information for the local model, it is necessary to interpolate grid-dependent information flom the regional model at the finer grid spacing of the local model. The transformation of information flom regional model to local model and the selection of boundary conditions and starting conditions for the local model are the most important issues using this numerical methodology (Townley and Wilson 1980; Ward et al. 1987; Buxton and Reilly 1986). A nested grid approach avoids the potential difficulties encountered with local grid refinement by separating the regional model flom the local model, and solving each model individually. Therefore, instead of solving very large, complex matrices, a problem is solved using multiple, smaller-scale local models. The local model will derive its boundary and starting conditions flom the parent model. Once the local model is created, it performs as an independent model. However, a major drawback in implementing the nested models is that the interaction between the parent and local models (of which there can be many) depends on the offline 33 analysis and processing of model modifications or simulation results flom the parent model to obtain the boundary and starting conditions for the local model. For example, once a new simulation is completed using the regional model, new boundary conditions and starting flow and solute transport conditions, if applicable, are determined for each local model at the local model grid spacing. Making modifications to models or processing simulation results for use in different scales of models can be very time consuming; especially when the problem is a transient-flow or transport condition, or there is uncertainty in selecting the boundary conditions. The effort involved may become impractical when the offline conceptual changes must be made iteratively or in more than one model. Because of this, applications of the nest grid approach are limited, in most cases, to very small number of submodels (e.g., 1 or 2), and are implemented with little flexibility. A DYNAMICALLY INTEGRATED HIERARCHICAL PATCH DYNAMICS APPROACH Li et al. (in press) recently presented an obj ect-oriented, dynamically integrated hierarchical patch dynamics paradigm (HPDP) for modeling flow and transport in complex groundwater systems. The new paradigm is particularly suited for modeling complex, clustered dynamics in a regional model. The HPDP represents a generalization of the nested grid approach and enables modeling of a large complex system as a network of hierarchically-nested patch models. For a 34 particular patch model at a particular scale, the HPDP is just a nested grid approach (Li et al. in press). The way to breakdown the problem and to relate a particular submodel to its immediate parent model is the same as that in a basic nested grid approach. Prior methodological research or empirical modeling experience related to nested grid can be directly used in the generalized HPDP. But under the HPDP, the patch dynamics models are dynamically coupled. The interaction between the regional model and all submodels are seamless. In other words, regional-model-simulation results (e. g., boundary conditions), and any changes to the regional model propagate automatically to all submodels, without the need for offline post-processing of data or simulation results. This transfer of information between the regional model and submodels is accomplished for each time step. It gives the modelers the perception of using a single model that provides high resolution dynamics at the speed that is very similar to that of a low-resolution coarse-grid model. The dynamic- integration eliminates the disconnection in a traditional disjointed nested grid approach and makes generalized hierarchical patch dynamics modeling practical. The obj ect-oriented implementation enables flexible, interactive creation of hierarchical patch models. The hierarchical modeling process becomes a recursive process by which a modeler can create incrementally and interactively a hierarchy of nested submodels in order to “zoom in” more and more details where it is required. 35 Our experience with hierarchical modeling and submodel development has shown that to maintain computational efficiency and robustness in modeling detailed flow dynamics around wells in a very large complex system it is best to keep the number of grid cells relatively small and grid spacing relatively uniform. This may result in the development of many submodels with successively-smaller grid spacing until it is possible to estimate head accurately at the well scale. However, with the HPDP, this process of creating an appropriate number of submodels with successively-smaller grid spacing is accomplished easily and in a way that is naturally intuitive (Li et al., in press). ILLUSTRATIVE EXAMPLES In this section we illustrate this HPDP-based method for modeling near well dynamics in large regional models through three concrete examples. Our first example considers a simple situation involving unsteady flow toward a single well. Our second example extends our examination to a more complex well-field system and the simulation of detailed flow dynamics around wells in a large regional system. Finally, we apply the HPDP to a 3D example involving cyclic pumping flom a well in a heterogeneous unconfined aquifer. 36 Single Well Example The simulation results obtained using the HPDP are compared to the analytical solution given by Theis (1935) to demonstrate the ability of the hierarchical modeling approach to predict drawdown at the radius of a pumping well. The analytical solution assumes that the confined aquifer is homogenous and isotropic, with impermeable boundaries bounding the top and bottom of the aquifer. It firrther assumes that the aquifer is two- dimensional and has infinite lateral extent. The objectives of this example are to 1) illustrate the process and concept of hierarchical modeling and verify its accuracy in predicting the drawdown at the radius of a pumping well, and 2) examine the relationship between the number of nested modeling levels and the submodel grid resolution in hierarchical modeling. In the finite-difference model developed to solve this problem, the aquifer is represented by a rectangular layer whose extent is 10,000 meters by 10,000 meters. No-flow boundary conditions are imposed along each face of the model. The aquifer parameters are as follows: the transmissivity and storage coefficient of the aquifer are 17.3 mz/day and 0.001, respectively; the discharge rate of the pumping well is 518.4 m3/day; the well radius is 0.1 meters (0.1m) and the duration of pumping is 50 days. Figure 1 presents the hierarchical-network and solutions to the single well example problem. The regional model and the six submodels have grid resolutions of 500m, 100m, 30m, 9m, 2.7m, and 0.81m respectively. Each model has a grid that is 21 rows by 37 21 columns in size. Figure 2 illustrates the drawdown comparison between the numerical model and analytical solution. Drawdown data flom every other node were used to prepare this graph. It is apparent, in examining the results obtained flom hierarchical- modeling-approach, that the numerical model closely-approximates the analytical solution. Also, it is important to point out that the results flom each submodel were obtained and visualized instantaneously. And, with hierarchical modeling, different drawdown values may be recomputed, displayed, and analyzed very quickly whenever the model stresses or other model parameters are changed. 38 M :1; Ax=500m 10000 8000 : x. Ylm) e. 000* zooox 0 11 114L411 0 ‘2000 4995 5000 5005 Figure 1. Presents the hierarchical-network and simulated heads after 50 days pumping for the single well example problem. Each model has a grid that is 21 rows by 21 columns in size. 39 10 - """"""""""" S A _ ”,—””" E 20 ra”"' l: 1 g — xxr” Exact solution All" , g d ,e”/’ “-0-" M1,;Ax=500m g 30; ’ —---El---- M3,;Ax=100m 4 ----A---- M?,;Ax=30m ----e---- M1,;Ax=9m 40 ----x---- M5;Ax=2.7m 11 —---o---- M?,;Ax=0.81m 50 411111111 1 11111111 4 11111111 1 11111111 14111111] 10'1 10° 101 102 103 10‘ Distance from well node (m) Figure 2. Illustrates the drawdown comparison at the end of 50 days between the numerical model and analytical solution. Each model has a grid that is 21 rows by 21 columns in size. 40 An alternative hierarchical model applied to the same single well problem is shown in Figure 3. In this case, each model has a larger grid that is 51 rows by 51 columns in size, requiring only three nested submodels to accurately predict the drawdown at the well scale. Figure 4 illustrates the drawdown comparison between the analytical solution and numerical hierarchical models. On this graph, only the simulated drawdown results flom every other model node are plotted. The alternative hierarchical solution show that there is more than one way to subdivide a larger regional model into several smaller submodels. The decision on how to subdivide a regional model depends primarily on the computational ability of the computer available to perform the analysis. One may choose to divide the regional model into larger numbers of smaller submodels when using a less powerful computer. Likewise, an approach using a smaller number of larger submodels would be taken when using a more powerful computer. This flexibility eliminates the longstanding, infamous “curse of dimensionality” in large—scale groundwater modeling and allows the modeler to simulate detailed flow dynamics, even around low capacity wells, in a large, regional groundwater system. 41 M :1; Ax=200m 10000 8000 I 6000 - Y(m) 4000 - 2000 - M :1; Ax=0.46m I I I I I I.- ‘n I' I I I I I I 111111111 4990 5000 501 0 Figure 3. Presents the hierarchical-network and simulated heads after 50 days pumping for the single well example problem. Each model has a grid that is 51 rows by 51 columns in size. 42 ’1-Illll T ..fll 1 J J 10‘ A _- Ea- """ A 'l E 20-4 c 7 x- g _ ”I, o ’g- '0 2 " ES 7 Exact solution a 30— 1 _ ----El---- M,,;Ax=200m : ----A---- Mf,;Ax=24m I 40% ---—X---- M3,:Ax=2.88m 4 ------.. M:,;Ax=0.46m .4 50_ 111111111 1 11L11111 1 11111111 1 11111111 1 11111111 10'1 10° 101 102 103 10“ Distance from well node (m) Figure 4. Illustrates the drawdown comparison at the end of 50 days between the numerical model and analytical solution. Each model has a grid that is 51 rows by 51 columns in size. 43 Well Field Example: The objective of this example is to illustrate the capabilities of the hierarchical-modeling- approach in calculating hydraulic heads for a large grid containing a well field with multiple wells each pumping at different rates. This example applies for groundwater- management problems, delineating the wellhead protection area (WHPA) for a well field, or prediction of groundwater-flow directions and velocities for solute-transport modeling. Computational difficulties encountered in modeling well fields center around calculating heads and groundwater flow rates where there are clusters of wells that pump at different rates and with different pumping schedules. In this case, the system may be very dynamic with spatially- and temporally-variable heads and groundwater flow. In areas far removed flom the well field, predictions using a regional model are usually accurate; however, closer to the well node, a finer grid is needed to model the system. Model grids of different extent and cell size are needed to model the system accurately so that the flow dynamics in a well field with multiple well clusters and variable pumping rates are resolved. Another difficulty in modeling a well field is representing the exact location of the individual pumping wells, especially those in close proximity to each other. With a regional model, or coarse grid submodel, clusters of closely spaced wells are grouped together in a single well node. A locally refined grid is needed to simulate the impact of individual wells on heads and flow rates. 44 The model simulation results are compared to an analytical solution for drawdown in a well field obtained by superimposing Theis solutions for each pumping well to verify the hierarchical-modeling—approach for a well field in a regional model. The distribution of the pumping wells and their pumping rates for the well-field example are shown in Figure 5. The regional-model domain is 100 km by 100 km with a uniform grid size of 4000 meters. As with the single-well example, the outer boundaries are represented by no- flow boundaries,. The no-flow boundaries in the regional model have been placed at this distance so they do not affect the simulated drawdown in the well field, satisfying the assumption of an infinite aquifer. The aquifer parameters are as follows: the transmissivity and storage coefficient of the aquifer are 200 mz/day and 0.001 , respectively; the well radius is 0.1 meter (0.1m); and the duration of pumping is 64 days. The discharge rates of the pumping wells are variable and are shown in Figure 5. 45 1 00000 75000 - Q E o E: Q 9 o Q12 00 1 50000 — Q 5 .0 Q11 Q 7..Q 2 0 Q4. Q 10 00 3 L 25000 Well ID Q1 Q2 Q, Q, Q5 Q, Q7 Q, Q9 Q1(, Q,,Q12 Pumping rate 0.25 0.25 0.15 1 0.25 2 0.2 0.5 0.05 0.15 0.5 0.5 Unit = x1000m 3/day 0 1 1 4 1 1 l 1 0 25000 50000 75000 100000 X (m) Figure 5. The distribution of the pumping wells and their pumping rates for the well-field example. 46 M f, ; Ax=15.85m M11 ;Ax-3.91m Figure 6. Hierarchical modeling layout which illustrates the relationship between submodels and their parent models. [In label of Submodel M ,1 a represents the generation level (grandparent), b is the parent index, and c is the kid index which is related to parent index]. 47 One possible example of the application of hierarchical modeling is shown in Figure 6. In this example, different levels of submodels are developed for each pumping well. Starting flom a coarse model we can successively approach the area of interest and obtain detailed information as close as the effective well radius. The diagram shown in the upper right quadrant of this figure illustrates the hierarchical relationship between parent models and subsequent submodels. The remainder of the figure shows the different model domains, pumping wells and simulated drawdowns. Lines showing the relationship between parent model and dependent submodel domains were not included for figure clarity reasons. However, an examination of the hierarchical tree and the different model domains should be sufficient to understand the relationship between models. The grid resolution flom the regional model to finest well-scale model varies between 4000 m to 0.49 m. The number of rows or columns in these models is typically 30 to 100, depending on the size of the model area and the uniform grid spacing for that model. Figures 7a and 7b show the comparison between modeling results for different submodel levels and grid resolutions and the analytical solution obtained by superimposing Theis solutions for the different pumping wells. The results are plotted along a cross-section profile (A-A) drawn through the regional model whose location is shown in Figure 6. In the drawdown comparison we show the results flom the regional model (M 11) through the finest-grid submodel (M 181). The area of interest in this problem is the pumping well shown in the center of model M18I (pumping well of interest). Figure 7a shows the 48 simulated heads along a profile through the entire well field flom model M12, including the pumping well of interest. Figure 7b shows simulated heads flom model M ,8, along the same profile through the pumping well of interest. However, in this figure, the lateral extent of the profile is limited to the extent of the domain of model M?1 in the X- direction. As with the single well example, these results show that, as the grid spacing is reduced, the numerical solution approaches the analytical solution. In addition, the results flom each submodel were obtained and visualized instantaneously. With hierarchical modeling, different drawdown values may be recomputed, displayed, and analyzed very quickly whenever the model stresses or other model parameters are changed. 49 a) Exact solution M111; Ax =40000m M121; Ax =1000m 3. - M12, Ax -250m 4 . _ M22, Ax —62.5m 5 . _ M42, Ax -15.65m 8 \ / 6 : A M41, Ax =3.9lm 7.5 J x M271; Ax =0.98m 3 8. _ 1 0 M11, Ax —0.49m 7 1 + l 4 A l . . . 1 1 1A 1 g. l a n A . l . TQM—I 20000 30000 40000 50000 00000 70000 80000 X(m) b) 10 - 9.5 E g 9 I 8.5 8 . 1 1 . I . . . . I 1 . 1 1 1 1 . 1 I J 57500 57750 58000 58250 58500 X(m) Figure 7 - The modeling results are compared to the analytical solution on a cross-section profile (A-A) for different submodel levels and grid resolutions (a)- large scale (b)- local scale. 50 Three-dimensional Example for Single Well: In this example we apply the HPDP to illustrate the capability of the hierarchical modeling approach in calculating hydraulic heads for a pumping well with a variable pumping rate in a large fully three-dimensional grid having heterogeneous media. A sketch of the problem is shown in Figure 8. The example problem consists of two conceptual layers that each represents an aquifer system. Conceptual layer 1 is an unconfined heterogeneous aquifer with a uniform thickness of 120 m and a mean hydraulic conductivity of 0.1 m/day. This upper layer has a lower overall hydraulic conductivity than the underlying layer. Conceptual layer 2 is a semi-confined aquifer having a uniform thickness of 80 m and a mean hydraulic conductivity of 1 m/day. Heterogeneity for each layer is defined by using 100 scattered points within the model domain having a hydraulic conductivity within one standard deviation of the mean and a correlation scale of 400 meters in the x-and y—directions. The steady state solution for the initial and boundary conditions shown in Figure 8 were obtained prior to conducting the transient simulation. The pumping well is located in Conceptual layer 2 and has a screen length of 2 meters. The well pumping rate is 545 m3/day (100 gpm) with a recycle period of 50 days. That is, the well pumps continuously for 50 days and is idle for 50 days. The transient solution is simulated for three recycle periods, a total of 300 days. A uniform time step of 4 days was used. 51 Qw A f 10000 m p a No flow boundary 0 '9 <9 0” 9° 00 I; ‘0 or? Z A 4,090 1 I v I i x / NoJ'flow boundary N I Mean k ,=0.1 mlday | I I I B,=120m I I s,=0.0011/m ' 1 . ‘ 89 ' 1» Mean k2=1mlday 1 | Bz=80m : l s,=0.0011/m T Figure 8- Definition sketch for three-dimensional example for single well. A finite-difference model with a uniform grid spacing of 100 meters and four computational layers for each conceptual layer was developed for the example problem against which the results of the hierarchical modeling at the well-scale were compared. 52 The final grid size for this uniform-grid model (Base Model) was 100 rows by 100 1 columns by 8 layers. The parent hierarchical model (M 11 , see Table 1) used a coarse grid having dimensions in the x- and y-direction of 500 meters and two computational layers for each conceptual layer, resulting in four computational layers in the parent model. As with the two-dimensional hierarchical approach presented earlier in this paper, additional submodels, consisting of finer girds are developed. The process of selecting each submodel area for the three-dimensional example problem is same as with the two-dimensional hierarchical modeling approach. In this problem, because of the three dimensions the grids of submodel are refined in the vertical direction (2 direction) by adding computational layers. The boundary condition for each new submodel in horizontal (x-y) direction is represented by a constant head boundary, whose values are interpolated hydraulic heads flom the coarse-model nodes flom the preceding parent model. The head values are also interpolated in the vertical dimension for new computational layers. A linear interpolation scheme is used for nodes where the aquifer is confined and a zero-order extrapolation method is used for nodes where the aquifer is unconfined. At every time step, the interpolated head value in each submodel node is updated flom each preceding parent model. A harmonic averaging scheme is used to represent hydraulic conductivity in three dimensions. The numerical grid configuration for each model is shown in Table 1. 53 Model name Lecrlrgth in x and y Grid'sizein x and y Number of computational layers in 1rectlon (m) direction (m) each conceptual layer Base Model 10000 100 4 M 1‘, 10000 500 2 M ,2, 6000 100 4 M ,3, 2000 50 8 M ,4, 1000 12.5 12 M f, 600 6.25 16 (only in aquifer 2) Table 1- Presents the numerical parameter set up for all the models in three dimensional single well example. Each submodel is defined within the preceding parent model. We can continue creating new submodels until the desired grid resolution is achieved. In this example, for comparison and illustration purposes, the grid spacing in horizontal direction has been reduced to 6.25 meters, which is about 100 times finer than the coarse model grid resolution. The grid spacing in the vertical direction has also been reduced to 5 meters, compared to 40 meters for the coarse-grid parent model and 10 meters for the Base Model. The simulation results for each hierarchical model are shown in Figure 9. For comparison purposes the solution for each hierarchical model are plotted with the solution flom the Base Model (shown in Figure 9 by the solid line). The simulation results shown in Figure 9 represent the transient drawdown and recovery at the cell in which the pumping well screen is located within conceptual layer 2. Results are plotted 54 flom the regional coarse model (M l',) to finest submodel level (M151). As a result of computational difficulties, we were not able to discretize and obtain a solution for the Base Model was at grid resolution that was less than 100 meters in x-y direction with 4 computational layers for each conceptual layer. The results show that coarse grid solution (M 1',) does not have the resolution necessary to resolve the drawdown at the well-scale during each pumping event. The solution using submodel (M121): which has same resolution as the Base Model, is better able to capture the impact of the pumping stress, these results also match the Base Model results fairly well. By introducing additional submodels (M ,3, and M 141) we improve the prediction of drawdown and recovery in the vicinity of the well screen. The A2 listed in the legend in Figure 9 represents thickness of each computational layer for aquifer 2 (second conceptual layer) in each model. In the finest submodel (M151), the addition of computational layers is performed only for the second conceptual layer of parent submodel. This allows us the ability to simulate only the area of interest without refining and modeling any area not necessary for this solution. 55 95 ._.-- l . |1 I I u 'I ‘ I I I u I ' 1 I 1 ., :1 . \tL l. “x I. I‘\ ' 1 85 B‘D—GQ It DBB‘GB'E :Gfl-Gaflu I 1 : : 1 80 I :I : ' l' . I II I I II .1 u '1 1' I ll l‘ 11 '1 1' 11 11 I' A 75 “ :1 I: u I. :1 E l .1 l t: :: 1 E “ '1 1l 11 .1 '1 w " '1 1' 11 .1 :1 :ii 70 1'8 'l' .8 'i ” '1 1 1 : 609% .‘060900 1090090 1 I I I 1 I ' 1 1 l 1 I 55 i I I I ' I I l l 1 1 l ' l l l I I : 1 1 1 1 : 60 1 i l l : . I 1 I I I ' ' 1 1 1 1 I ' I I 1 1 I l 1 l I I ' 55 ‘ 1 x l \ l x‘x-i-x-lr-il \"-X-x-~x-x—' x‘X-x-x—xgi 50 7 1 1 l 1 l I l 1 1 l l l 0 50 100 150 200 250 300 Time(day) — Base Model ----- M :, ;Ax=500m; Az=40m o M i, ;Ax=100m; Az=20m --1:i-- M2, ;Ax=50m;Az=10m --0-- M :, ;Ax=12.5m ; Az=6.7m —-*-- M ?, ;Ax=6.25m ; Az=5m Figure 9. Illustrates the drawdown comparison at the well screen between the numerical model and Base Model solution. 56 An analysis to investigate the impact of selecting the extent of submodel boundaries was performed. For comparison purposes, the results were compared to those obtained using the Base Model. In this analysis, three submodels each having equal grid spacing in three dimensions (x, y, and z) and but with different submodel dimensions in the x- and y- direction (Table 2) were created. The Base Model and each submodel have a uniform grid spacing of 100 meters. The parent model grid spacing is the same as that used in the previous analysis. Model name ”352;: 0:? y Grid spacing (m) Base model 10000 100 M f, 10000 500 M ,2, 6000 100 M ,2, 8000 100 M ,2, 4000 100 Table 2- Numerical parameter setup for different submodels The simulated drawdowns and recoveries for the different submodels (M 121 , M 122, and M 123), and the Base Model are shown in Figure 10. 57 92 - I ’ ‘ 1 '1 g 1 91.6 — i E l 1 i i 91 2 H e‘ . Base MOdOI I I I“ ‘ l- ...... Mi. : LX=B°°°m I E 1), .................. M :2 ; Lx=6000m 1’ 'o “I“ ‘ 8 908 4 \{x __ ..... — M I: i m l I L"! ‘ 90.4 4 \‘:§.:~>..: ~ - 1 \;\\:“~ 3‘ e 5 ~\ ..... *----~““....‘._.‘..7‘_.:~--3- ' i 90 _ --~--..-.\,;';j_'”41:77:75 7i l L L l 1 l 1 I 1 I 89.6 o 10 20 30 40 50 Time(daY) Figure 10. Represents the drawdown comparison between the numerical model and Base Model solutions. These results show that the smaller the submodel extent the greater the difference in drawdown and recovery flom the Base Model solution. The difference in simulated drawdown and recovery values between the submodel solutions and the Base Model solution is partially due to loss of information at the boundary. But even for the 58 submodel with the smallest extent, the accuracy gained through nested grid refinement is much more significant than that that lost due to errors introduced through inaccurate submodel boundary conditions. SUMMARY Nested grids are frequently applied in groundwater modeling but most applications are limited to a very small number of modeling patches (e. g., 1-2 patches) and nested levels (e. g., 1-2 levels). Even for a simple nested system with just one submodeling patch, practical implementation proves to be difficult when iterative conceptualizations and evaluations are needed. This difficulty is caused by the fact that the submodels are disconnected and to iterate back and forth among different disconnected models (which may require in themselves iterative implementation, analysis, and postprocessing) is difficult, if not impractical. In a recent paper, Li et al. (in press) generalized the nested grid approach to a systematic, dynamically integrated hierarchical patch-dynamics paradigm (HPDP). The HPDP solves a complex system as a network of hierarchically-organized, small patch models. Since solving many small, uniformly gridded models are much easier than solving one single large model on potentially strongly nonuniform gn'ds, the HPDP is very efficient, substantially relaxing the infamous “curse of dimensionalit)?’ and the associated computational bottlenecks in modeling complex systems across multiple scales. 59 In this paper we apply the dynamically-integrated HPDP (Li et al. in press) to address a difficult groundwater modeling problem — modeling well-scale drawdown and detailed near well dynamics in a large, regional groundwater model. We employ three examples to demonstrate the HPDP’s ability in modeling well-scale drawdown and the detailed near well dynamics in large regional models. The results show that the HPDP has the following advantages for modeling complex well dynamics: 0 The HPDP allows simulating detailed well dynamics in a large regional model without solving large matrix systems. 0 The HPDP eliminates the restrictive assumptions inherent in the analytical correction approaches. 0 The HPDP eliminates the strangling I/O and intermediate offline operations in the traditional, nested grid-based nmnerical correction approaches. 60 REFERENCES: Anderson, M. P. and W. W. Woessner. 1992. Applied groundwater modeling: simulation of flow and advection transport. San Diego: Academic Press. Buxton, H. and T. E. Reilly. 1986. A technique for analysis of ground-water systems of regional and subregional scales applied to Long Island: US. Geological Survey,U.S. Geological Survey. Li, S.G., Q. Liu, and S. Afshari. In press. An Object-Oriented Hierarchical Patch Dynamics Paradigm (HPDP) for Groundwater Modeling. Environmental Modeling and Software. Mehl, S., and M. C. Hill. 2002. Development and evaluation of a local grid refinement method for block-centered finite-difference groundwater models using shared nodes. Advances in Water Resources 25 (5):497-511. Peaceman, D. W. 1978. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation. Society of Petroleum Engineers Journal 18 (3): 1 83-194. Peaceman, D. W. 1983. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Non-Square Grid Blocks and Anisotropic Permeability. Society of Petroleum Engineers Journal 23 (3):531-543. Prickett, TA. 1967. Designing pumped well characterstics into electrical analog models. Ground Water 5 (4):38-46. Pritchett, J. W., and S. K. Garg. 1980. Determination of Effective Well Block Radii for Numerical Reservoir Simulations. Water Resources Research 16 (4):665-674. Thiem, G. 1906. Hydrologische Methoden: Leipzig, Gebhardt. Thies, C. V. 1935. The relation between the lowering of the pizometric surface and the rate and duration of discharge of a well using groundwater storag. Trans. Amer. Geophys. Union. (2):519-524. Townley, L. R. and J. L. Wilson. 1980. Description of and user's manual for a finite- element aquifer flow model AQUIFEM-l. Cambridge Massachusetts: Ralph M. Parson Laboratory Technology Adaption Program. Trescott, P.C., G. F. Finder, and SP. Larson. 1976. Finite-difference model for aquifer simulation in two dimensional with results of umerical experiment. US. Geological Survey, Chap. C1, Bk.7. 61 Ward, D. S., D. R. Buss, J. W. Mercer, and S. S. Hughes. 1987. Evaluation of a Groundwater Corrective Action at the Chem-Dyne Hazardous-Waste Site Using a Telescopic Mesh Refinement Modeling Approach. Water Resources Research 23 (4):603-617. 62 PAPER NO. 3: ASSESSING THE HYDRAULIC PERFORMANCE OF A BIOREMEDIATION CURTAIN FOR PLUME G AT THE SCHOOLCRAFT SITE, Ml USING AN INTEGRATED HIERARCHICAL GROUNDWATER MODEL. Soheil Afshari', Shu-Guang Li', Qun Liu', Phanikumar Mantha', and Michael J. Dybasm 1Department of Civil and Environmental Engineering Michigan State University, East Lansing, MI 48824 2Center for Microbial Ecology Michigan State University, East Lansing, MI 48824 ABSTRACT: This paper describes three-dimensional multi-scale flow and transport modeling effort for a contaminated aquifer at the Schoolcrafl site in western Michigan. The ultimate objective is to design a broadly applicable “biocurtain” which be applied to remediate this site using a recirculation well gallery installed normal to groundwater flow. Characterization of flow and transport phenomena during the delivery and treatment process is crucial for designing the biocurtain. But this process is significantly complicated by the potential complex interaction of the flow dynamics across a number of relevant spatial scales present at the site —the “delivery unit-scale” 0(10m), the “biocurtain scale” —0( 100m), the “site scale” —O(1000m), and the “regional scale” — ' To be submitted for publication, latest version is available at www.cgr.msu.edu/~lishug 63 0(10,000m). At the unit scale, the closely-spaced injection and extraction wells create complex and rapidly varying head variation and a grid resolution on the order of 10 centimeter is required to resolve the local spatial dynamics. For the accurate modeling of transport and delivery process an even smaller grid size may be required to drive the numerical dispersion to an acceptable level. Creating such high resolution, 3D unsteady model is only feasible for a relatively small area (e.g., on a unit scale). On the other hand the presence of heterogeneity can potentially induce significant unit-to-unit interaction and evaluation of this interaction requires modeling simultaneously multiple delivery units or the entire biocurtain on a fine grid. Additionally, potential impact from seasonal irrigation further complicates the site characterization since it may require expanding substantially the modeling domain from the site-scale to on the order 10,000 meter. In this paper, we apply the hierarchical patch dynamic paradigm (HPDP) to model efficiently the complex interplay of the flow dynamics across multiple scales in the vicinity of the Schoolcrafi site. The resulting integrated hierarchical modeling system is used to characterize, evaluate, and optimize the performance of the biocurtain and the three dimensional delivery system with and without the presence of large-scale heterogeneity and seasonal inigation. 64 INTRODUCTION Over the last decade, bioremediation has become increasingly common and accepted in applied groundwater remediation efforts (Azadpour-Keeley et al. 2001). In situ remediation is one of the alternative strategies that have been in a great deal of interest in groundwater remediation technologies. The in situ remediation process carries on by injecting and delivering nutrients and microbes to contaminated zone and degrade the contaminants to harmless end products (Phanikumar et al. 2005). Our research case study in this paper is, a volatile organic compound (V CC) contaminated Plume G site at Schoolcraft Michigan (Figure 1). In situ remediation strategy for clean up the plume G was chosen by series of extraction injection well gallery installed across the width of the plume and normal to groundwater flow. The objective of remediation design is to deliver nutrients and microbes to a treatment zone, which is called as “biocurtain”, and maintain a uniform microbial concentration across the vertical and horizontal extent of biocurtain during the treatment period. Numerical models have been used to design subsurface delivery strategies (Lang et al., 1997; Hyndman et al., 2000; Scheibe et al., 2001 , Phanikumar et al., 2005). Due to limited access to field data and characterization of contaminated sites, we apply the numerical models to maximize the use of valuable data and provide design guidance and predict the delivery efficiency. Utilizing numerical models will allow us to perform sensitivity analysis on variable physical parameters and evaluate design efficiency for a 65 .Amoom owaucmméfisfic 00803 33 beam .«o cosmos ._ 0.5me L 38 335m , mmmhmz 7: mq iili YI 66 long run, under different circumstances. Due to high grid resolution that is required to accurately model the flow and transport phenomena at delivery zone, numerical models are dimensionally limited to small modeling area. Such models can be very useful for modeling labrotary scale or small field scale. Ifthe regional flow is uniform we can take advantage of symmetry and only model one unit of extraction injection system, which has smaller dimensions and will allow us to use smaller grid size for modeling. Creating a high-resolution model for a small area is affordable. Then we can apply the single unit design to the full biocurtain scale. This is a traditional single scale modeling approach. However, in reality delivery can be a significant problem because the transport of contaminant and agents needed for remediation are profoundly affected by heterogeneous properties of aquifer sediments (Hyndman et al. 2000). Due to heterogeneity of media and complexity of flow dynamic of each system, full-scale biocurtain flow dynamic could be nonuniform. In this case a single unit design will not be adequate and we need to expand the scale of problem to entire width of the plume or full biocurtain scale. We need to model full scale biocurtain [300 m] and the influence area of biocurtain to predict interaction between individual sets of pumping and injection systems [unit set includes one pumping well and two injection well]. Additional complicating factor is possible seasonal irrigation events that could cause a temporal and spatial nonuniform flow dynamic and generate a unit set to set flow and transport variation. Impact of irrigation complicates the issues and adds an additional level of complexity with a different scale of nonuniformity across the full biocurtain 67 scale. To address the impact of irrigation on fiill biocurtain scale we need to expand the problem boundaries to a larger scale [10000m] and model beyond the biocurtain scale. Scale Dependent Problem In general this problem is a complex interplay between multiple scales of flow and transport dynamic systems. Multiple scales in this problem includes 1) Unit scale that includes a single pump and two delivery well system which covers a small area (Order of 10m) with high small scale variation. 2) Biocurtain scale which includes the entire pump and injection well gallery and it extents to the full width of plume (Order of 300m) 3) Local scale which includes the area beyond full biocurtain scale and captures the influence area of biocurtain scale (Order of 1000m) 4) Regional scale which includes possible irrigation pumping region that could impact the biocurtain flow dynamic (order of 10000m). Each of these problem scales is very complex and difficult to solve, they also interact with each other at different scales. It requires a detail flow and transport model to capture the multiple scales of variability and complex flow dynamic. A higher level of difficulty is to model the interaction between these scales and ability to address the interrelated issues between multiple scales simultaneously. Therefore we have to deal with different levels of difficulties in this problem. Another fact is that in most of the cases treatment zones are designed to remediate for span of large time (0 decades). Therefore our design stability and delivery for a long period of time is highly crucial for treatment efficiency; specially if some large scale or seasonal flow dynamics occur during the treatment period. Having the capability to model multi-scale flow and 68 transport variability is essential to investigate the delivery process across spatial and temporal scales. Research Questions and Challenges A major question arises that how can we simultaneously predict the impact of large flow dynamic variation on small scale flow and transport phenomena? How can we address the delivery process in small scale and mean while model the large scale flow dynamics? Is it possible to decouple the flow and transport phenomena at different scales and be able to address them at individual scales, and still keep them connected across the scales? Is it possible to fill the gap between large scale and small scale flow and transport dynamic and provide a ladder to transit the information across multiple spatial and temporal scales? These research questions are very challenging using the traditional approach (Li et al. in press). Traditional approach includes 1) Local grid refinement, which employs a set of non uniform girds to solve a multi-scale large problem. Large number of non uniform grids are needed to model a large multi-scale problem, which cause an ill-posed matrix systems that leads to convergence and solution failure. 2) Telescopic Mesh Refinement method (TMR) (Ward ed al. 1987; Mehl and Hill 2002) which breaks down a large multi- scale problem to individual single scale smaller problems. Advantage of TMR is that it makes modeling a large multi-scale problem possible by solving multi-single scale problems. However, disadvantages of TMR are that it only applies to few numbers of levels. In addition TMR requires post processing the data for flow and transport model, 69 and works offline for multi-scale problems. In general TMR could only apply to limited numbers of local models and data post processing and visualizing the data in a real time is very time consuming. For problems that consist of huge scale disparity it is impossible to apply the traditional approach to solve the problem. Modeling complex interaction between flow and transport phenomena at different scales is still one of the advanced challenges in numerical modeling research. In this research paper we present an innovative hierarchical patch dynamic paradigm (HPDP) (Li et al. in press) modeling tool to address the impact of multi-scale flow and transport variation on the delivery process. The combination of delivery design process for this specific research site and applying the innovative HPDP tool to address the impact of multi-scale variability across spatial and temporal scales on delivery process makes this research work very interesting and unique. Remarkably in this research project we had the opportunity to utilize site characterization data and field tracer data from previous research at Schoolcraft, MI site (Hyndman et al. 2000; Dybas et al. 2002; 2003; Graulau 2003). Our modeling strategy in this paper is organized in three different categories. Basically due to insufficient regional detailed information and high uncertainty about large scale data, we begin our modeling practice with simplified assumptions. Then gradually we build on top of the base case model and investigate different levels of complexity in the 70 model. In the ultimate regional multi-scale model we assess the most complex and nonuniform flow situation which involves in multiples scales of variability across different models. We test several hypotheses in the multi-regional scale model. In the first base case model, we applied uniform flow and took advantage of symmetry to model the local scale flow and transport dynamic. In the first numerical model, a tracer scale model [0 10m] created to calibrate against the tracer data. The objective of this model was to better characterize the hydrogeological aspect of contaminated site and obtain calibrated parameters to employ for larger scale model so called biocurtain scale. In the biocurtain scale model [0100m] we added additional level of complexity by predicting the delivery process for multiple sets of extraction and injection wells. Purpose of this model is to predict the delivery at full scale spacing design, and investigate the impact of multiple sets of extraction injection wells on each other. Additional level of complexity in this model, canied on by imposing small scale heterogeneity. Number of realizations for different heterogeneous scenarios simulated in the biocurtain scale model to better understand the complex hydraulic system between extraction and injection wells during the recirculation event, and also provide guidelines for design purposes. Finally the ultimate regional multi-scale model was created to assess impact of large scale flow dynamic on local scale delivery process. In this model we applied HPDP capability and exceeded the traditional single scale modeling approach for design purposes. Traditional single scale modeling approach evaluates the simplified local scale 71 model and uses a factor of safety to design the recirculation system and model the delivery process. However it fails to provide a big picture from regional scale dynamic and to connect its impact on local scale delivery process and also to test design stability. The purpose of regional multi—scale model was to address more complex and nonuniform large scale regional flow dynamic and interactively predict the impact of large scale flow dynamic on biocurtain scale delivery process. Scope of Work Specific research modeling objectives and efforts in this paper are: 0 Create a three dimensional flow and transport model and calibrate the model against field tracer data fi'om previous research study [tracer test scale]. 0 Expand the modeling domain to full biocurtain scale and investigate efficiency of delivery at heterogeneous media for a larger scale [biocurtain scale]. 0 Utilize an innovative HPDP modeling tool to predict the impact of multi-scale variability across different spatial and temporal scales on delivery process and design stability [regional multi-scale]. 72 OVERVIEW OF THE SCHOOLCRAFT FIELD SITE The Village of Schoolcrafi is a small rural community located approximately 16 km south of Kalamazoo, MI, USA (Figure 1). The unconfined glaciofluvial aquifer underneath the village has been contaminated with organic and metal compounds as result of previous industrial and commercial activities in the village. Schoolcraft Plume G is 1.6 km long and 360 m wide and it extends from about 9 to 27 m below ground surface (bgs). The water table at this site is roughly 5 m bgs and the base of the aquifer is defined by the top of a thick clay unit at 27.3 m bgs. For more information regarding regional and local hydrogeologic conditions refer to Lipinski 2002; Hyndman et al. 2000; Dybas et al. 1998; Mayotte et al. 1996. Site Characterization and Pilot Scale Tracer Test Site characterization Previous research presents a fairly smooth and uniform head distribution with gradient of 0.0011 on Schoolcrafi site (Figure 2) (Graulau 2003). Figure 3 shows the main features found during the visual classification of hydraulic conductivity analysis (Graulau 2003). There are two distinct layers with in the 27m (82 ft) vertical extent of aquifer. A top layer which consists of (19m) 62ft is mostly fine sand changing to medium sand as one moves deeper into the formation (Graulau 2003). Deeper layer with average thickness of (6m) 20 ft consists of mixture of coarse and medium sands with gravel and cobbles at specific depth (Figure 3, and 4) (Graulau2003). 73 A N LEGEND —262.0- Potentiometric head in meters Figure 2.Groundwater contour map in the VOC contaminated region[source Graulau 2003] 74 MW-8 MW-5 MW-2 0m _ WELL GRADED GRAVEL (GW) m CLAY (CL) [:1 WELL-GRADED SAND (GS) _ SORTED SANDS WITH GRAVEL (SP) LEGEND 75 Figure 3. Cross section parallel to groundwater flow at vicinity of tracer test [10m distance] [source Graulau 2003]. Depth (m) L 141 :1 l l l l l 1 i l I l 4 L L L ¥ l l I o ‘ 0.1 0.2 1 0.3 J 0.4 0.5 Conductivity (cm/sec) Figure 4. Measured hydraulic conductivity for specific depth [source Graulau 2003]. There are partially some high conductivity layers with in the (6m) 20ft thickness layer but due to uncertainty in extent of this layers we decided to use the entire 20 ft thickness layer as one conceptual layer. Later in sensitivity analysis section we will investigate impact of heterogeneity on delivery process. The variogram analysis from previous research studies shows a horizontal correlation scales of 18m and vertical correlation 76 scale of 4.4 m with variance of 0. l 3 (Graulau 2003). Average total porosity from laboratory repacked samples was 0.331 with a standard deviation around the mean of 0.050. Maximum and minimum values were 0.504 and 0.208, respectively (Graulau 2003) Field tracer test Field tracer test was conducted to estimate flow and transport parameters that influence distribution of solutes within this aquifer section. The main objective of these experiments was to study the impacts that the stratigraphy exert upon contaminant migration and distribution in this region. An “inj ection-extraction” strategy was used to deliver tracer solution to the saturated zone of the aquifer. To achieve uniform tracer delivery in the geologic formation, a 4hour injection time was chosen since preliminary numerical models predicted an approximate 80% tracer breakthrough in the extraction well at the end of the pumping period. The screened interval for the extraction and injection wells in the network is between 18.3 to 25.0m bgs [62-82ft]. Details of the field tracer experiments are provided in Table 1 (Graulau 2003). Tracer breakthrough at delivery wells reached a maximum relative concentration (C/Co) of 1.8 at 4hr of pumping. Relative concentration in delivery wells are greater than 1 as soon as tracer breaks through in the flux control well. When this condition occurs, the injected concentration will increase by the relative concentration coming from the flux control well (Figure 5) (Graulau 2003). 77 F luorescein Tracer injected concentration (t=0) 100 ppm injection into DWI-D W2 pumping rate 1.3 x 10’ m3/s (each well) 4 hrs pumping time Tracer injected into DW's by adding the solute to the pumped water fi'om F C W ’s, Tracer test procedure .see Figure 5 for well notation Table 1. Details of tracer experiments (source Graulau 2003). 78 meow $333-52:ch 8.58 .8393 838 w: Bog 05 .«o 3984 .m onE G N _ (PA—E 6 o _ - . o 3 Injection E .1 0 o 5 _ .( No-tlowboundnry L A l . . . l L O O 5 10 15 X(m) Figure 6. Conceptual model set up for tracer scale model 82 Selection of numerical parameters for delivery modeling requires a good understanding of physics of the problem. Delivery process is a transient process therefore modeling result is extremely sensitive to grid resolution and time step selection. These two parameters should be selected such that to reflect the physics involved in the tracer delivery and mixing process. As we mentioned earlier in order to accurately capture the flow dynamic we need a grid resolution in order of cm. Because the well diameter installed in the test was 15 cm (6 inches) we set the model grid resolution at 15 cm in horizontal direction to create a model which numerically be a close representative of physics of the problem. In the tracer test, injection process happens almost instantly therefore we need to use a very small time step to model the delivery process. Transport modeling is very sensitive to numerical parameter selection. We select time step such that to minimize the numerical dispersion introduced into the modeling result. We examined that by utilizing a time step of 1.2 min, transport solution converges. For flow model we compared the results between transient flow and steady state, a good agreement between steady and transient flow model results was obtained. Therefore we modeled the flow at steady state condition. Numerical and physical parameter set up for this model is presented in Table 2. 83 Regional gradient 0.001 Aquifer-type Unconfined Aquifer-top layer thickness 62 ft (18.9m) Aquifer-bot layer —thickness 20 fl (6.1m) Hydraulic conductivity top layer 0.027 cm/sec Hydraulic conductivity bottom layer 0.17 cm/sec Effective porosity top layer 0.3 Effective porosity bottom layer 0.175 Physical digaersion in longitudinal direction 0.05 m Physical dispersion in transverse direction 0.005m Domain-size 15m by 15m Grid resolution in horizontal direction 0.15 m Time step, total simulation time 1.2 min, 240 min Table 2. Input physical and numerical parameters for tracer scale model. Calibration Process and Discussion After creating the numerical model now the task is to calibrate the model against tracer data. Calibration process is always challenging; it requires a good understanding of conceptual physics and numeric associated with the problem. In the breakthrough calibration process one can always argue that there is a unique set of physical parameters to fit the data. In the calibration process we alter hydraulic conductivity, effective porosity and physical dispersion to produce the breakthrough curve at extraction and injection well. Our major target was to calibrate against the data from extraction well, because the data from extraction well were obtained from underground but data from injection well were measured before injecting to ground, and also numerical model predicts the results from underground. 84 During the calibration process we observed that effective porosity and hydraulic conductivity control the pick magnitude and can shift the entire curve up or down. Specially breakthrough curve was very sensitive to effective porosity, because effective porosity can vary the pore velocity and in this problem pore velocity during the pump and inject even dictates the delivery performance. We also had to take into account the impact of physical dispersion. Although during pumping event advection portion of transport has major role in the delivery but due to averaging the hydraulic conductivity over the entire conceptual layer we had to add the physical dispersion to compensate the hydraulic conductivity averaging. Altering the physical dispersion can impact the early arrival in the breakthrough curve, it also smooth out the curve and yields to a less pick concentration at the end of the curve. Given our general physical understanding of the problem, after testing multiple scenarios and adjusting the physical parameters to fit the breakthrough curve at the extraction well we came up with the physical parameters which is illustrated in Table 2. Figure 7 presents the breakthrough comparison between model and data, which illustrates a relatively good match between model results and data. 85 200 - 180 1001 4O 80 120 160 200 240 80 Lill.r.l.r#l.r J ..l..LL 0 4O 80 120 160 200 240 Tin. (mil) Figure 7. Comparison of model results and data for breakthrough curve at extraction and injection wells. 86 BIOCURTAIN SCALE MODEL Due to complex interplay of different scales of variability, modeling the delivery process could be very challenging. In scaling up the model, we investigate the problem for a larger pumping rate and longer recirculation period and wider spatial spacing between the wells. Therefore we need to know the breakthrough value after certain recirculation period to come up with the design parameters. In the design process of biocurtain scale in order to minimize the treatment cost we need to expand the distance between extraction and injection wells as far as possible with in the range of pumping rate of 100gpm and recirculation period of 24-48 hours. In the meantime we have to maintain a uniform delivery concentration for the entire extent of biocurtain in vertical and horizontal direction. The tentative desire concentration delivery is 50 percent relative concentration achievement in 24 hours of recirculation period. Dr. Dybas in his report to department of environmental quality (DEQ) (Dybas. et al, 2003) proposed using a triangle layout with two injection wells and one extraction well. Figure 8 illustrates the spacing and the geometry between the wells. We will evaluate performance of suggested design with given spatial distance between extraction and injection wells with 100gpm pumping and 50 gpm injection. The ultimate objective of delivery is to create a uniform microbial concentration treatment area across the width of plume (Figure 9) (for more information regarding treatment area see Kim et al. in review). Plume G is the largest plume in the series of 87 plumes in Schoolcraft, MI. It is spread 1200fi (366m) perpendicular to groundwater flow. In order to treat about 300 m width of the plume and employing proposed spatial distance in extraction injection wells set up (Figure 8), we potentially need to install l7sets of 10 m pumping and injection wells to cover core portion of width of the plume G. The major purpose of this model is to investigate the impact of different sets of extraction injection wells on each other. An interesting issue in the biocurtain scale model is the interaction between different units of extraction injection systems. The “squeeze effect” between different units may change the flow dynamic along the biocurtain and also it may impact the biocurtain scale delivery process. To better understand the squeeze effect we need to model multiple units of extraction injection wells. In the biocurtain scale model we take advantage of flow uniformity and only model 4 sets of extraction injection wells. We chose 4 sets because two sets at the edges will represent the edge condition of the biocurtain and the other two sets in the middle will represent the condition in the center of the biocurtain. 88 _._ Injection Well Extraction ‘6— Well Injection __ Well Flow E [s Injection _..L_ Well ' E Extraction O Well F Injection __ Well Figure 8- Well spacing layout for extraction injection well system for potential filll scale biocurtain design (Dybas et al. 2003). 89 04 N()\ 2005 J-D Model Results 100 80 60 Y(m) 20 20 40 Concemratimmpm) 160 140 120 100 80 60 4O 20 60 80 100 X(m) Figure 9. A typical control volume [black box-solid line] for comparison purposes. Light blue zone is treatment area. 90 A conceptual sketch of 4 sets is represented in Figure 10. We applied the same physical parameters from tracer scale model in the biocurtain scale model; the numerical model set up parameters is shown in Table 3. The head contour map and delivery concentration plots are illustrated in Figure 11, 12. Regionaljradient 0.001 Aquifer-type Unconfined Aquifer-top layer thickness 62 ft (18.9m) Aquifer-bot layer ~tlrickness 20 ft (6.1m) Hydraulic conductivity top one 0.027 cm/sec Hydraulic conductivity bottom layer 0.17 cm/sec Effective porosity top layer 0.3 Effective porosity bottom layer 0.175 Physical dispersion in longitudinal direction 0.05 m Physical dispersion in transverse direction 0.005m Domain-size 120m by 120m Grid resolution in horizontal direction 0.75 m Time step, total simulation time 0.2, 24 hour Table 3. Input physical and numerical parameters for biocurtain scale model. To address the “squeeze effect” issue we compared the breakthrough results from the middle extraction well in the biocurtain scale with a single (one set ) 10 m case. The breakthrough curve for a 24 hour extraction and injection is shown in Figure 13. Comparison reveals a significant difference in breakthrough magnitude. The relative concentration for the breakthrough curve at the extraction well in the single 10 meter model set reaches to % 30 (c/co), however for the same physical parameters the 4 sets of 10 meter model reaches to %58 relative breakthrough concentration at the extraction well. 91 100fl 80- OE b 2' 5° ‘5 e " E g .3 ' 0E A ‘5 i h g 3 E .g 5 >40 __3 _ I. .E E .0 I. g . E 0 a l. E 20 - a V4: i 0 -< No-flowboundnry L .19 r L. r . . .1 r r 20 40 BO 80 100 Figure 10. Conceptual model set up for a 4 set of 10 m pumping and injection wells. 92 Head(m) 260.1 5 260.1 260.067 260.05 260.036 260 259.95 259.9 259.86 259.85 259.829 259.8 259.75 259.7 259.65 259.6 259.55 Y(m) O 20 40 60 80 100 X(m) Figure 11. Steady state head contour map for 4 sets of 10 meter extraction injection well at biocurtain scale model. 93 100 +- 80— 60— Y(m) 20- .4..41..A1...1...19#.1-n 0 20 4O 60 80 100 X(m) Figure 12. Concentration map for 4 sets of 10 meter extraction injection well at the end of one day pumping. 94 60- 40? / Concentration (ppm) 0: o l 20L / r / 10. ———-4setsof10m I Onesetof10m 0 b A_ A_ A__ A A A A A J A A A l A A A J 0.2 0.4 0.6 0.8 1 Time(day) 160 f F / : / 150 r / ’ t / . / 140 L / I E / .. / A 130 L , ’ v , / .5 120 j / 3E ‘ / 5110 I ’ - / c . _ / 8 : , 100 f ’ 90 t —— Onesetof10m C ———-4setsof10m 80 4 A A l A A A l L A J l A A A l A A A l 0 0.2 0.4 0.6 0.8 1 T‘Irne(day) Figure 13. Breakthrough prediction at extraction and injection wells for one and four sets of 10 meter well spacing model. Notice the difference between breakthrough results due to squeeze effect which occurs in the 4 set of 10 meter model results. 95 One alternative way to explain the results is that squeeze effect or interplay complex flow dynamic between different sets of 10 meters create higher gradient. And higher gradient dictates the tracer transport. For the single case of 10 m, local gradient imposed by pump and inject is less than 4 set case, therefore it takes longer to pickup the concentration and it reaches to a smaller pick value. Therefore for homogenous design purposes we can use this breakthrough range for the design. Sensitivity Analysis and Discussions Heterogeneity could enhance complexity of the interaction between different units. Impact of heterogeneity on delivery process has been illustrated by various researchers (Freyberg 1986; Hess et al. 1992; LeBlanc et al. 1991; Garabedian et al. 1991 ; Adams and Gelhar 1992; Boggs et al. 1992; Rehfeldt et a1. 1992; Hyndman et al. 2000). To better understand and investigate the impact of heterogeneity on delivery process we simulated tracer model for heterogeneous media. From previous research we used the correlation scale of 20m in horizontal direction (Hyndman et al. 2000; Graulau 2003). Comparison results for different breakthrough curves for multiple variances and realizations provides a flavor of the range of uncertainty in the breakthrough curve (Figures 14 through18). The range is between %40 to %60 relative concentration at the extraction well. Results also show that the bigger the variance the larger the range of breakthrough value. Different breakthrough curves for heterogeneous case shows that the homogenous (constant k) breakthrough curve is always close to ensemble mean curve. 96 This may suggest that for design purpose the homogenous breakthrough curve potentially could be adequate for design purposes. We may investigate more stochastic analysis in feature research. In addition we investigated the impact of effective porosity on transport phenomena. We altered the effective porosity positively correlated with the natural logarithmic of hydraulic conductivity (In k). In general when effective porosity varies positively correlated with In k we expect a drop in breakthrough concentration. However, after comparing the results to the case where In k was the only variable, results didn’t agree with our expectation. A higher breakthrough range, for variable effective porosity which was positively correlated with In k was obtained. 97 —+— mum-«10m VII'1 —-°—- Wk4'“°“°"' “to.“ 7o . ---------- vnr1 Runs ’ 7o ---------- woman's f _____ W1 m " ————— No.15 RIM :____ 1m , ———W0.15RLm r:- _ ........... VI” [“31 74'” 50 - ----------- W015” . ./-- ~—---»--~-~- «Amman , .-.-.-._.--.._........ ”1 Rm I: I. m E '4 1; I E 50 ———-o—— EW'I“. 'u 40 , ’ 40 so i so 20 20 10 1 o 1 L L 4 1 0.2 0.4 0.0 0.8 1 TIME-Y) —~— Comhntk4utlot10m arms —°-— Wk4o¢bd10m l'lVlfOJS 7° .......... "053L233 70 _ ---------- wOJGnVIRLZlG _____ wObRLZIG I —-—-— No.15anLZIO ' _____ “053m 801—__ No.15nvnrRLM // _ .......... “05m“ : - ----------- v.0.15nvnrRLZI1 ---_._._._.._.__ ”0.5 Rm " " ’ ----—-~————~ VIOJGI‘IVI’RLM ’ —— Minoan / I (mu) 4° . so I ' 2° . 10 0' "' - l . 1 A . 1 J 0.2 0.4 0.6 0.8 1 Figure 14. Comparison breakthrough curves at the extraction well [point E2 Figure 19] for different hydraulic conductivity variances and multiple realizations with homogenous [constant k] case. Notice that homogenous case result agrees with ensemble mean for all the realization. Next 4 figures presents each breakthrough curve separately. 98 70- 60 '- ’1 . —o— Conttantlduhof‘wm XI ------ monsatzss ,I’ ; —— - ——- varo.1s thts /’ . “,v 50 — — — var0.15RI.Z#2 l,’ . , / ',/ " — ----- - var0.15RLZ#1 ,” , ", A I ’, 0 g varO.15Rl.Z#4 ,’ ,/ . / a. ' —0—- Ensemble moanvaro.15 ,’ //,_” 7: 4° - ,.//// s ' ,7, g '- ’1’ I 04’”, / c ' I / 04 8 30 '— I, I ll,” / .. I ,/ g " III X" / 0 I ' i I: / 2o_— [I 7/ 1o; o ‘ . 1 1 l I n A n l l 1 1 1 4 1 1 l 0.2 0.4 0.6 0.8 1 Time (day) Figure 15. Comparison breakthrough curves at the extraction well for different realizations of hydraulic conductivity with variance of 0.15 with homogenous [constant k] case. Relative concentration value varies from %49 to %60. 99 70 : I 60 ~ —e—Constantk4utsof10m /" _ —————— var0.15nvarRLZ#3 - » —-—var0.15nvarRLZ#5 ° 50 L — — var0.15nvarRLZ#2 t — ----- - var 0.15 n var RLZM A [ var 0.15 n var Run // a _ ——+— EnumblomannvarwiflrkvarOJS .. 40 ~ ”I c . ” s + g . 3 30 — c . ° .- o . 20 — r L. 10 — L. i. t 04 1 1 ‘ L %‘ I 0.2 0.4 0.6 0.8 1 Timo(day) Figure 16. Comparison breakthrough curves at the extraction well for different realizations of variable effective porosity positively correlated with hydraulic conductivity variance of O. 15 with homogenous [constant k] case. Relative concentration value ranges from %45 to %65. 100 60 _ —e— Const-MkAeeuoHOm ’ . ------ var 0.5 RLzes ,’ _ — - — var0.5 RLZ#5 / f— . — — var0.5 RLzsz ’ , ¢ ’ 50 j _ ----- - var0.5RLZ#1 ” /./' _ var0.5 RLzu ,’ .~ ’ ’ - —.___ Ensemble meanvar0.5 I” / , i / Concentration (ppm) 0.2 0.4 0.6 0.8 1 Time (day) Figure 17. Comparison breakthrough curves at the extraction well for different realizations of hydraulic conductivity with variance of 0.5 with homogenous [constant k] case. Relative concentration value ranges from %46 to %70. Notice that the breakthrough range for variance of 0.5 is bigger than the range with variance of 0.15. 101 70 T r I r I \ —e— Constantk4setsof10m I . 60 X 17 F ------ ver1 RLZ#3 ,r x - —-—— var1RLZ#5 ,I’ ’K/ ’ : — — ver1 RLG r’ x 50 _ —- ----- - var1RLZt‘l /’ ’ var1RLZ#4 I” , E _ —-—— Ensemble mean var1 ,’ a I a . ... 4o - / s - / E r / c ' ’l 8 3° ~ I’ , / C > ’I O r ’ / o . ,’ ' . I / 20 — I / 10 r- 0 k ‘ A 1 l k ' l ‘ 1 l o 0.2 0.4 0.6 0.8 1 Time (day) Figure 18. Comparison breakthrough curves at the extraction well for different realizations of hydraulic conductivity with variance of 1 with homogenous [constant k] case. Relative concentration value ranges from %41 to %61. Notice that the breakthrough range for variance of 1 is bigger than the range with variance of 0.15 and 0.5. 102 The breakthrough curve at the extraction well predicts the delivered concentration at a specific point. Our objective in the biocurtain design is to provide a desire relative concentration along the biocurtain horizontal and vertical extend. To better investigate the delivered concentration during the recirculation process we investigated the breakthrough results at a possible weakest points in terms of concentration coverage. This point is located between two edges of each unit system (point W Figure 19). As we expect the breakthrough results for homogenous case at point W shows a smaller concentration value (Figure 20 through 24). The interesting part is that the homogenous breakthrough value shows the lowest breakthrough value comparing to most of the realization’s predictions. This is different from the homogenous breakthrough curve at extraction well where the homogenous k case was in a good agreement with ensemble mean value. Additional interesting point is that in point W due to high level of uncertainty we obtain a higher range in breakthrough magnitude. We need to highlight that for the given correlation scale and physical parameters selection, homogenous case for design purposes, is a conservative design approach and leads to potential worst case scenario performance. 103 ,_— __—,—~__._______..__ ._ A 100 — 0 ow Flow 80- —————-——l ................................................................. Y(m) 1 bob 0 a Mini wr- ° ; 11 2°.” W o u Cori-hut Medicaid-Iv '5 I - I ”TE .E - I 0 Is . E I . I0I D 2°‘ " .4”: Is“ I I 0“ 0| No-llaw ry 1 I .All. ..1. 1 1 0 20 40 80 80 100 Figure 19. Conceptual sketch for a 4sets of 10 meter extraction injection model representing location of profile cuts and points E2 (extraction well) and W (potential weakest concentration coverage point). Point W is located along the biocurtain length in the middle of extraction and injection well and horizontally it is located in the middle of two different sets of extraction injection. It is located in the intersection of profile cut D- D and CC. 104 —.— zooomumofmn‘ nVlr0.16 —._ MMRAseteollom VII'1 100 ---------- moanwmsRLzss 10° ---------- 2009m1RI-Z“ ,- 90 —— EnsemblemeenZOOOnwkverOds 90 —_ MMMRN1 ,I/ - ------- zooen ver0.15 RLZM - ------- 2009 "'1 Rl-Zfl ,,-’ I,» 80 - ----------- zoos nver0.15Rl.Zl1 /. 80 ------------ 2009 human ‘,-’ f/ A 70 ——— zoos "NO.15RLZ'2 l/ . 70 ———— zooowmm ,-’ _,-- / g ————— zoos nvuO.15Rl.ZlG ,o’ g 60 ' 60 50 50 40 4O 8 so 30 20 20 10 1O 0 o ' " 41 - 0.2 0.4 0.6 0.8 1 Tin-(div) ——+— zooecommuuuorwm “'05 _._ MIMsehoHOm are.“ 100 ---------- 2000w0.5RLZi3 100 ---------- zooewo.1sRLzss 90 EneenblemeenZOOOkvero-G 90 —~— EneemblermenZOOOka-ts - ------- soon was run ,,-' — ------- zoos varo.15RLzu so _ ........... zooowosRLzm // 80 - ........... 2000ver0.15RLZl1 A 70 ——— 2000 VII’O.5RLM /./' I; A 70 ——— 2000 No.15“ g ————— 2009 was Rm x/ g ————— zoos vero.15RLZl§ ,,-’ . . ' so . 50 4O 30 20 10 O Figure 20. Comparison breakthrough curves at point W [Figure 19] for different hydraulic conductivity variances and multiple realizations with homogenous [constant k] case. Notice that the results for homogenous case doesn’t agree with ensemble mean of all the realizations. Due to high uncertainty at point W [typical weakest point] the breakthrough range is significantly bigger than the range at extraction well [point E2 Figure 19]. 105 100 :- 90 _— 80 — 70 l A I —o— oonstantk4setsof10m ./ E . ------ zooovaro.15RLzss // 3 6° : __._ Ensemble mean 2009kvar0-15 /' I: I —-—-- 2009 var0.15RLZ#4 /’ 3 50 _ — ----- - zoos varo.15RLzs1 /' E — — 2009 varO.15RLZ#2 ,/ ../ / 3 _ — - — zoos var0.15RLZ#5 // /.o / g 40 L— /, './' / . , ,/ .4 , o - '/ './‘ a / 30 ;- ’/ ,./ ’ / . / 'o/ / ./ . / ./ ’ / 2o — / ,/‘ ,/ ./ / .- / ‘/' '/ I : o/ a/‘/’ ‘—"‘ _ / ,/' , ....... 1o _ /’ /, W.-::2:::'.'.:::::::::'.22‘-"""1", ,/"»” ‘_”,,','-"""'.'.';;zczttizzzztz ....... 0 s ''''''''' :7"..',‘.".".‘.1'.'.'.'.‘. ....... l %l A A I 02 04 0.6 08 1 Time (day) Figure 21. Comparison breakthrough curves at point W [Figure 19] for hydraulic conductivity with variances ofO. 15 with homogenous [constant k] case. Homogenous case predicts the lowest breakthrough value %9, however the range for other realizations varies from %13 to %69. 106 100 - 90 '— 80 l , I '/ : / 7O ;- ’/ E i —e— 2009constmtk4setsof10m // P ------ 2009 nvar0.15 RLZ I 3 ’ Q 60 — / ,/ v : —.— Ensemble meen2000nvarkver0.15 /' é ; —-—-- zoosn varo.15RLzu /' / g 50 _— — ----- - zoos nver0.15Rl1#1 / . c ; —— — zoos nver0.15RLZ#2 /' g 40 L — - — zoos nver0.15RLZ#5 ’/ o " ./ , ’ o : ./ / ’ / 30 L ,/ . / , / t '/ , l/ / I / 20 — / ' , . /. ’ / . /' ’ / A. ’ / 10 r / l / .....'12:::::::::1:313:33”, ” ’./ ’ / “...31131133231133“: ____________ / / ..---'2133:”;;:::::::.U ----------- 0 "113-:1-t‘ze'l‘ue‘ ‘1 ————— A 1 A A I A A I 02 04 0.6 08 1 Time (day) Figure 22. Comparison breakthrough curves at point W [Figure 19] for different realizations of variable effective porosity positively correlated with hydraulic conductivity variance of 0.15 with homogenous [constant k] case. Homogenous case predicts the breakthrough value %9, however the range for other realizations varies from %6 to %79. Notice that one of the realizations drops below the homogenous case’s result. 107 100 - C so :— - ’/' .- ’/ so _— ’/ ' —o— 2009mstantk4setsof10m // 7o - ------ zoos varo.5 RLzss /’ /' A —.— Ensemble mean 2009 kvaro-s /’ ,./ i so . _-—-- zoos varo.5 RLzu /’ // / 3: ~ -— ----- — zoos ver0.5 RLZ#1 /' / / g . —— —- zoos varo.s Run '/ /" / 50 f — - — 2009 var0.5 RLZ” ./ x/ / ’ g . / ,/ / ’ s : ./ / c 40 — / a'/ I/ o ' . ,/ o _ / / / ./ l- / /" ’/ 30 F . ,. . / : / / - , : / '/ I/ 20 - /”,/ _ / : '/ ' ’ : '9' '/ 10 ~ ..// /,/ : "//a I/W‘_.;2113331112122112212 ........ . fl' / ,....'-.'::::::"’—'.:211311;” """" O b _-:.':,:_‘.j,,'-,'.'-?"" """"" A l A A l A I 02 04 0.6 08 1 Time (day) Figure 23. Comparison breakthrough curves at point W [Figure 19] for hydraulic conductivity with variances of 0.5 with homogenous [constant k] case. Homogenous case predicts the lowest breakthrough value %9, however the range for other realizations varies from %9 to %88. 108 100 — /, l/‘ / 90 ;- ’/ E ./ . / / 80 - ———0—— 2009constentk4seuof10m /. / r ------ 2009var1RLZ#3 /' / / t —._ Ensemblemean 2009 kvar‘l /' /' 7° 7 —-—-- zoos varmtzn /' ./ / E ' - ----- - zoos var1RLZt1 /’ ,/ / n. * _ ' / e .. - _-_::::::::: / ,/ g i m /' / / / g 50 :— /' /'a/ / ’/ . I .' / g r ’/ '/ / /, c 40 _ / ’/ / l 8 /I /o / : /' ,/ / / 30 '- i /' .‘ /I : //x , C //o' , / 20 - I ’ C // ,/ I ..//'/ / 10 "" //l ’ -::::> - r /I/ .....-:::::;:‘.:: ............. b l/ / ...—-—-'"'"..—:;::':;::::::'.‘-15'-""’ """" P /" ’M.-".‘.'.1'.'.'.C ----- O ...... :..’:' ””” .:::_:1.;££. '''' l l L l A l 02 04 0.6 08 1 Time(day) Figure 24. Comparison breakthrough curves at point W [Figure 19] for hydraulic conductivity with variances of 1 with homogenous [constant k] case. Homogenous case predicts the lowest breakthrough value %9, however the range for other realizations varies from %9 to %100. Notice that by increasing the variance the range of breakthrough results increases. 109 Basically in this problem variability in hydraulic conductivity and effective porosity will be in favor of design, because it will help the particles to randomly migrate around the boarders of weak point and that will lead to a higher concentration breakthrough curve at weak area (point W). This may not be true for all the cases but for most of the cases could be applicable. To better investigate discontinuity of the biocurtain along width of the plume and highlight the inadequate concentration coverage area we compared the concentration results at the end of day one delivery, along the biocurtain (profile cut D-D in figure 19). Different realization comparison shows that for most of the predictions homogenous case predicts the lowest concentration (Figures 25 through 28). If tentative desire standard concentration coverage along the biocurtain is %40-50 of relative concentration, a gap of 1 meter long with relative concentration of %10-35 occurs at the edges of unit sets (weak area Figures 25 through 28). Profile cut D-D shows that, due to high uncertainty at weak coverage area, variability concentration values for the weak area is much higher than the concentration range at the pick area. 110 .‘\ . nver0.15 —— Constant i'\7-\‘.r.::,;-\ — — RLZ#nver0.15 \.,-\-~:. r.\\,_‘ —-—---—-- RLZ s n var 0.15 30 _ “\ — — RLZ# n ver 0.15 \ ................ RngnwoJ‘r, , .1; — ----- RLZ#nver0.15 ...----""‘""'....4’ '. «('N" """ 51' 1|: ' . . 1 \)_;= 1] I‘.' w.- Inm — _én— 60* 1 '5 1 .i s \x-(‘T’Rx 1 7,-_r\ 3 3 .. ’ \ —————' E l .... ...:— .. ‘\ I,’“' ‘.\ > -— \ .1 .' 40- ...——»—-"‘"/‘" <._ “" _————' ___/ . .3 . ...“ )...:’I; _V M. 20— 3"" ” ' )‘2 1 /./l b /. :7!”- -1"‘"=‘3"/ 3W o 40 so 120 160 ConcentratloMppm) Figure 25. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable effective porosity positively correlated with variable hydraulic conductivity with variance of O. 1 5 and also the homogenous case [constant k] [4 sets of 10 meter extraction injection well model]. Notice that pick concentration values drops slightly comparing to the case of only k variable with variance of 0.15. 111 _ . _____ . Constant - fl _ — -— RLZ#var0.15 so _ \‘ “-.--_-....-.... RLZ # var 0.15 _—— RLZ#var0.15 ...................... RLZ # var 0.15 ~\U _ -------- RLZ#var 0.15 ------ ,_ Y(m) o— _ f 201 '3' 0 40 80 120 160 ConcentratloMppm) Figure 26. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable hydraulic conductivity with variance of 0.15 and also the homogenous case [constant k]. Notice that minimum concentration value obtains at the middle of two extraction injection set. Range of uncertainty in the results drops when we move toward the injection well. 112 , “I,” ——.—_ Constant \\ _ — — RLZ# var 0.5 o . K 3".“~~\ -.-.-.-..-.-.-._..... RLZ# var 0.5 8 ._ \\‘ \\ \-\‘ — -— — RLZ# var 0.5 ,,.....-.—-- -.. ”a“; ....................... RLZ# var 0.5 ‘ .A. --~~-.... _ ........ RLZ# var 0.5 zukxx‘ _ \ 1.1.x r/ .N- 60— ,5) >- " - 40— ——-—-.‘.'.'" ') \2 i2 J 20— 1’; ‘ 7; /-‘}1 r .-' / , ’3'; 0' A J l l l I 1 41 o 40 80 120 160 Concentratlonmpm) Figure 27. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable hydraulic conductivity with variance of 0.5 and also the homogenous case [constant k] [4 sets of 10 meter extraction injection well model]. Notice that range of predicted concentration value increases by increasing the variance from 0.15 to 0.5. 113 \ \ _.___ Constant """"""" - —— — — RLZ#var1 80 _ \“\\ -~----—~---- RLZ# var1 1 , — — — RLZ#var1 «M:} 33 "J ................... RLZ#var“ .L““ ' ;-\\ — ----- - RLZ#var1 7 *\ .. \ s‘ .2 I. v “F J; ) .— ..hflf/ . %‘k 2 :5: i \\ >. \ 4o — ’- .......... .31 . - .plr‘fl‘ . . zo— - _ 1' /—"';; é ...» .y’ / r ‘ / / / c . ‘ I 0 40 80 120 160 ConcentratioMppm) Figure 28. Profile D-D presents the concentration distribution along the biocurtain for multiple realization of variable hydraulic conductivity with variance of 0.5 and also the homogenous case [constant k] [4 sets of 10 meter extraction injection well model]. Notice that range of predicted values increases by increasing the variance. 114 Additional point is that for the case of variable effective porosity, we obtain about %10 less relative concentration comparing to the same variance and realization case where only hydraulic conductivity was variable. The prediction agrees with our general understanding from impact of variable effective porosity with hydraulic conductivity on delivery performance. By comparing other profile cuts (Figures 29 through 36), we can conclude that uncertainty by velocity variability dictates the performance of the design and delivery process. At the extraction point there is less uncertainty variability in velocity and that explains why we obtain acceptable concentration coverage from all the cases. However, for design purposes the concentration at the weak spot may drop below the minimum standard design concentration coverage. Results reveal that deterministic homogenous case could be utilized for a conservative design approach. To firrther investigate effectiveness of delivery, we compared the average delivery concentration at the end of one day for a typical control volume (see Figure 9). Typical control volume is the zone which covers both highest and weakest area of coverage at delivery area. We ran 5 different realization for three different variances. The average concentration with minimum and maximum concentration obtained in the typical control volume are illustrated in Table 4. Tentative desire delivery objective is to reach an average relative concentration of 40-50% at the end of 24hours at a typical control volume; however the average relative concentration for constant k case and average of all realizations are %116 and %121. Since effective porosity is also one of the major key 115 parameters in the delivery process; we investigate impact of variability in effective porosity. Results show that by varying the effective porosity which is correlated with In k [variance of O. l 5], we can cover an average of % 119.5 relative concentration. 116 100 r C-C var 0.1 5 Constant — — — RLZ#1var0.15 -—-——-—-— RLZ#Zvar0.15 __— RLZ#3var0.15 ....................... RLZ#4varO,15 — -------- RLZ#5var0.15 80 l Concenhatlon(ppm) Figure 29. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous case with variance of 0.15 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that results for homogenous case predicts the lowest. The range of predicted realizations at pick vary from %20 to %105. 117 C- C ver 0.5 120 P r —-e—— Constant ; 2 _ _ — RLZ # var 0.5 100 __ ————--—- RLZ#var0.5 _ __ — RLZ#var0.5 ....................... RLZ # var 0,5 80 _ — -------- RLZ # var 0.5 Concentrationmpm) Figure 30. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous case with variance of 0.5 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that due to high uncertainty the range of predicted values are very large. The range of predicted realizations at pick varies from %20 to %122. 118 C-C var 1 —<>— Constant 120 ~ — - — RLZ#var1 ' ‘ --~—-- RLZ#var1 L ——— RLZ#var1 ....................... RLZ#var1 100 F — -------- RLZ#var1 Concentratiomppm) Figure 31. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous case with variance of 1 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that homogenous case predicts the lowest concentration across the profile and range Of results for other realizations vary from %20 to %128. 119 c-c nvar0.15 100 60 Concentration(ppm) 40 —s—— Constant — — — RLZ#nvar0.15 -——-———--— RLZ#nvar0.15 —-——— RLZ#nvar0.15 ....................... RLZ#nvarO.15 — -------- RLZ#nvar0.15 Figure 32. Profile cut C-C presents the concentration distribution for different realizations of heterogeneous porosity and conductivity with variance of O. 15 and also homogenous case. This cross passes the possible weakest zone of concentration coverage. Notice that one of the realization predicts value lower than homogenous case and the range of prediction varied from %9 to %118. 120 B-Bvar0.16 160 - ——e— Constant _ _ — RLZ#1var0.15 140 _ ——--—-- RLZ#Zvar0.15 ———— RLZ#GvarOAS ....................... RLZ#4varO,15 120 — -------- RLZ#GvarOJS ConcenbatloMppm) on o so 4o 20 ° 20 " 40 so so 100 X(m) Figure 33. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous case with variance of O. l 5 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. Notice that there is a good agreement with homogenous case prediction and other realizations. 121 B-B n ver 0.15 —o— Constant 160 - — — —- RLZ#nvar0.15 -—---———— RLZ#nvar0.15 140 - ——— RLZ#nvar0.15 ' ' ....................... RLZ#nvaron 120 ~ '1 1 \\ — -------- RLZ#nvar 0.15 .. 1 [.1 g 100 — ill \ ‘=' . ,1! 3 --1 ‘ ‘E 80 — .11, ‘ r 3:1 1;, 1 40 h 1 20 - H1 1 1 I ’ ‘, H 0 0 20 I 40 60 80 100 Figure 34. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous porosity positively correlated with hydraulic conductivity with variance of 0.15 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. 122 160 140 120 Concentrafion(ppm) on o 60 40 20 T l l 1 1 l B-B var 0.5 —s— Constant _ _ — RLZ # var 0.5 ......z—w RLZ # var 0.5 _ _ — RLZ#var0.5 ....................... RLZ # var 0,5 _ ........ RLZ#var 0.5 Figure 35. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous case with variance of 0.5 and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. Notice that by increasing the variance the range of prediction at the pick area increases. 123 B-Bver1 1 + Constant — — — RLZ#var1 160 _ -~—-———-—- RLZ#var1 ——— RLZ#var1 140 - ..................... RLZ#varl 1 — -------- RLZ#var1 120 - %— 100 ~ 5 so — i 1 0 60 _ 40 - 20 - o . . . 0 20 80 100 Figure 36. Profile cut B-B presents the concentration distribution for different realizations of heterogeneous case with variance of l and also homogenous case. This cross passes a possible strong zone of concentration coverage [injection zone]. Notice that by increasing the variance the range of prediction at the pick area increases. Results for case of variance 1 has higher predictions and higher ranges comparing to cases of 0.15 and 0.5. 124 We further investigate the delivery in treatment zone after recirculation. Figure 37 presents the concentration contour for three different time periods, one day during recirculation, and 10, 14 days after shutting down the wells and modeling transport under natural gradient [no reaction]. Comparison results show that center of the concentration moves with speed of 0.8m/day. It also shows that biocurtain readjusts itself after shutting down the recirculation system. Transport under natural gradient will help the delivered concentration to spread more towards the weak points and that will create a more uniform biocurtain in the treatment zone. Results show that after 10, and 14 days we have an average concentration of % 104 and %108 at the typical boundary volume in the treatment zone (Figure 9). Conductivity n correlated with Variance 0.15 0.15 variance 0.5 1 Average Average Average Average min, max min, max min, max min, max Realization #1 127 128 126 127 46,165 58,161 55,163 46,165 - Realization #2 121 133 122 120 48,155 63,162 61,155 52,154 Realization #3 1 12 91 1 12 109 84,156 2.5,143 13.4, 162 7,165.20 Realization #4 127 126 130 135 58,162 63,161 46,167 74,170 Realization #5 118 123 l 18 1 16 38,152 41,157 37,149 36,147 Table 4. Spatially averaged relative concentration with minimum, and maximum value obtained in the typical control volume in biocurtain scale model for 5 different realizations with 3 different hydraulic conductivity variances for bottom layer. Second column represents the case which effective porosity for bottom layer was varied spatially correlated to natural logarithmic of hydraulic conductivity. 125 8) 1]] Concentration(ppm): 50 100150 I'll-don (ll Constant head boundary ‘em' ‘ntneadboundau Arr” Constant head boundary N O A_AA4 A4 0 . 1 Combat h‘ead‘ boundary 1W (I) .- A- A .— ‘ ‘160 Constant head boundary A A A A A E b A AA160A 3-3 100— l Mad w 0 50100150 0 50100150 Concentration (ppm) 1— A A A D A A A A A r— A .A.A.A..LA..A.A.LA.AA.AJ 0 50100150 0 50100150 Concentration (ppm) M 0 50100150 0 50100150 Concentration (ppm) v—I.v V'VW 1 3......“ 0 0 CH: M 50100150 0 D-D W 50100150 0 50100150 M 50100150 50100150 0 50100150 Figure 37. Concentration layout and profile cut for a) 24 hour period during the recirculation system, b) 10, and c) 14 days after shutting down the wells. 126 We will utilize the information obtained in sensitivity analysis to provide guidance for next phase of field full scale tracer test. In addition we can also provide detail information to set the location for piezometers to measure and identify the weak coverage points from extraction and injection well. REGIONAL MULTI-SCALE MODEL So far most of our modeling effort is obtained in a relatively local scale and closed flow system with simplified assumptions. Basically we have applied regional gradient and included local sources and sinks posed by recirculation system. We have ignored any possible large scale flow nonuniformity (i.e: irrigation event) and its impact on delivery in biocurtain scale. Specially delivery design system will be operated for a long term (i.e 25 years), therefore we need to have the capability to address a large scale flow dynamic and its impact on biocurtain scale delivery and evaluate efficiency of the design during and after recirculation event. One may raise a question that how can we model large scale dynamic and mean while address the local scale flow and transport phenomena? In the introduction of this paper we mentioned to the challenges associated with multi-scale modeling and the status of practice in multi-spatial and temporal scale modeling. In the next section we will introduce the HPDP (Li et al. in press) tool which is incorporated in IGW (Li et al. 2003; 2006). Then we will apply the HPDP to create a regional multi-scale model [OlOOOOm] which is extended beyond the biocurtain scale. We use regional multi- scale model to assess the impact of large scale flow dynamic on biocurtain scale delivery during and after recirculation. 127 HIERARCHICAL PATCH DYNAMIC PARADIGM (HPDP) In this project, we apply a hierarchical patch dynamic paradigm (HPDP) (Li et al. in press) to evaluate the biocurtain recirculation system and model the complex interplay of different dynamics in multiple scales. HPDP represents a powerful generalization of TMR, which could arbitrary apply to number of levels and recursively solve multiple local models or submodels. HPDP breaks down a large multi-scale problem to number of small problems and simultaneously model flow and transport dynamics at different scales. HPDP fi'ees us from any post data processing and allows us to refine the gird to the extent of required resolution and enables to visualize the coarse model and local models results in multiple scales in a real time. HPDP allows us to explore large number of scenarios, which is important for design purposes. Utilizing HPDP we can perform a sensitivity analysis for different design parameters at different scales and address the effectiveness of bioremediation and clean up process at different scales. HPDP satisfies required highlighted points for modeling the interplay of multi-scale problems for this specific site. However difficult portion of the work is how to implement the HPDP in a modeling fi'amework and apply it to a large multi-scale problem? HPDP is adopted in a new software system so called Interactive Ground Water (IGW) an innovative approach to groundwater modeling and model-based design at the plume site. Parallel computing paradigm which is implemented in IGW software system make the HPDP approach possible. The new modeling approach makes it possible to quickly 128 explore a large number of alternative complex modeling scenarios, enabling rapid groundwater modeling, immediate display of results to allow visualization, a process that greatly facilitates group discussion and analysis. This approach has tremendous advantage in that the new paradigm has the potential to significantly reduce lead times for model implementation and visualization during the design phase. In addition we can utilize the high performance 3D visualization tool in IGW to better visualize and understand 3D delivery process in the real time. This will allow us to test multiple scenarios and provide intuitive guidelines for design purposes. Regional Multi-scale Model Construction To create the hierarchical regional multi-scale model, and estimate the regional modeling domain, initially we need to identify regional major sources and sinks in vicinity of biocurtain area. Figure 38 shows three high capacity wells at vicinity of biocurtain. Detailed information about the wells and their distance from center of the plume is presented in Table 5. We considered a possible scenario that all the wells are at a full pumping rate and reached steady state. This was a reasonable assumption to predict the long term regional influences on local delivery. Modeling area and input parameters are tabulated in Table 6. We applied the same physical parameters fi'om biocurtain scale model and applied the regional gradient to the model. We also created the full well gallery across the width of plume. Well gallery consists of potential 17 extraction well and 34 injection well at the rate of 100gpm pumping and 50 gpm injection. After creating 129 ~ , -+-"“‘"*"""3|‘ ‘ e I ‘ 1‘ V ‘ ‘ 3‘ " MM” ‘2‘. . Production \ Plume F&G '1 1:> 200 m scale Regional groundwater flow Figure 38. Regional multiscale modeling domain and existing high capacity well [circled] at vicinity of biocurtain. 130 the regional large scale model, we need to create and add submodels to get to the vicinity of biocurtain. We need to create multiple levels of submodel to transfer the information from regional scale model to local biocurtain scale model. For more information about HPDP method and model construction see Li et al. in press. Figure 39 shows the level of submodels and their patch network. It shows that how different submodel at different spatial scales and different grid resolution can related to each other and transit information fi'om large scale model to a local scale model. Input parameters for submodels are illustrated in Table 6. By applying HPDP method we can break down a regional scale model with grid resolution of 120m and construct smaller area model with finer resolution. The finest scale model in this problem has resolution of 60 cm which is comparable to a 4 set of 10 m scale model. Time step in this model is 0.2 hour. The uniqueness of the regional multi-scale model is that it provides regional scale flow information and mean while predicts the local flow and transport phenomena at smaller scale. It allows us to assess the large scale impact on local scale model and also monitor the delivery process across the biocurtain scale for different modeling scenarios. Regional multi-scale model allows us to test the complicating factors of regional flow in local scale delivery. I-IPDP capability provides the opportunity to assess several complex scenarios in regional scale modeling and predicts their impact on local scale delivery. 131 Y(m) VA“, .. .\‘ I \\ I ‘\B M - 25 , WM: . m f 1 1 :1 1 L” J 1‘13; I; ‘ 1 ‘ {.-. \ ‘._ 1 r." \ I l ‘ {'1 A .. J.- \ " t 1 ~. g \ A‘ J) \ . A \' ‘. J. \ “ ‘4’ A _.. ._L-tu_»~ _.‘ .. ,- . .-" ' .-" . 9360 1 9380 w . l I o A l _ r . ‘ I ~~‘ ~v . 1‘, I 7‘;'~. | . ~ < . . ‘ ._ _ \A__“ E .1 5“ AA .A " ”TV “2" YIV‘FI' ‘1 ’ .- I" A . _‘ Figure 39. Hierarchical modeling layout for regional multi scale model which illustrates the relationship between submodels. 132 Table 5. Distance and pumping rate of existing high capacity wells in the vicinity of biocurtain. . Distance from Well # Pumpmg rate gpm biocurtain (m) [RI 900 1426 IRR2 900 563 Production 300 927 Model name Model length in x Grid size in x and y and y direction (m) direction (m) Mr‘r 11910, 9765 120,122 M121 4876,4267 60, 61.5 M131 2175,1700 25.2,25.7 M141 495,712 6.2,8 M151 200,430 2.5,3.9 M161 50,38 06,06 Table 6. Regional multi-scale input parameters for different submodels. Physical parameters for this model are same as biocurtain scale model. Results and Discussions To investigate the impact of large scale flow dynamic on local scale delivery process, we compared the breakthrough curve at a typical extraction well for the cases of with and without irrigation event. Figure 40 shows the breakthrough curve results from the 133 hierarchical regional multi-scale model. For the case of presence irrigation pumping, we obtain a relatively close curve to the case of no irrigation scenario. However the impact of large scale dynamic raised the curve by %2-3. We also compared the concentration distribution along the profile cut DD for the cases of with and without irrigation pumping during the recirculation period. We need to mention that for the case of regional multi-scale model we set the maximum injection well to %100 of relative concentration and we didn’t add the additional concentration from extraction well. Results are not comparable with 4 sets of 10 meter scale model, however they are comparable with cases of with and without irrigation event for the hierarchical model. Figure 41 shows similar results to the results of biocurtain scale model, at the edges of each extraction injection system, concentration may vary to the range of %10-40 relative concentration. However this gap is in the range of one meter. In general, results reveal that for this test case [homogenous] during the recirculation system local flow dynamic is very dominant and large scale flow dynamic doesn’t have much impact on delivery process. Although this statement is NOT an absolute conclusion, we know that in reality the homogeneity assumption is not valid, and real aquifer is heterogeneous. Results for different physical parameters may vary. To further investigate the impact of regional scales dynamic on local scale delivery we may need to test flow nonuniformity and heterogeneity scenarios to set an uncertainty bandwidth for our generalized conclusion. 134 607 so — A 40 — E r 5 _ C 0° '- 5 so — 1: . 8 C o _ ° 20 — 10 L 1. 1 0 . A . A l A A A l A L A J A A A l 0.2 0.4 0.6 0.8 1 Time(day) Figure 40. Breakthrough comparison for an extraction well from regional multi scale model during the one day recirculation for cases of with and without regional irrigation. Notice that results for two cases are very close. 135 2995 2990 2985 2980 2975 Y(m) 2970 2965 2960 2955 0 —o— with pumping : E — — - without pumping E /‘j i =37 : — : j— ’ E F E _. f “AAAIAAAIAAAI..Al...|..A| 20 40 60 80 100 120 Concentration(ppm) Figure 41. Profile cut D-D from regional multi scale model presents the concentration distribution along the biocurtain for cases of with and with our regional pumping. 136 We performed additional comparisons to investigate the impact of irrigation on local gradient at vicinity of biocurtain after shutting down the extraction injection systems. Figure 42 presents the change of head across the profile cut A—A (Figure 39) for multiple pumping scenarios. The cut passes between two sets of extraction injection systems, and shows the change of head for cases of a) natural gradient b) presence of irrigation and no recirculation c) presence of irrigation and recirculation (1) presence of recirculation and no inigation. The results show that the natural gradient is 0.001 (m/m) and with presence of irrigation gradient changes to 0.0016 (tn/m). Therefore, local gradient increases by %60 due to background irrigation. The increase in local gradient due to background irrigation has a significant potential in biocurtain transport. It may impact the design of recirculation periodic time. In overall big picture design we need to consider the impact of large scale pumping event on local delivery. 137 260.8 —- I 1' 260.6 260.4 260.2 — 260 - E I — — — ProuncooflRR nomch'cuiaflon v . — - — - - Natural Gradient '3 259.8 - —o-— Only recirculation No IRR 0 Z ——t—— Presence of IRR and recirculation I r- 2596 - 259.4 _ 259.2 ‘ 259 L. 258.8 i- 1 k l 1 14L 1 l 1 1 1 1 L 1 1 1 1 l 1 1 1 1 l 1 9350 9360 9370 9380 9390 X (distance m) Figure 42. Head comparison across the profile A-A Figure 39 for four different scenarios a) natural gradient, no pumping b) presence of irrigation and no recirculation c) presence of irrigation and recirculation (1) presence of recirculation and no irrigation. 138 SUMMARY AND CONCLUSIONS In this research paper we investigated the impact of large scale and local scale flow dynamic on delivery process for an in situ remediation project at Schoolcrafi, MI. Our modeling strategy was to create a simplified base case model and investigate more complex cases by scaling up the numerical model. We constructed three levels of numerical model in this project to address different aspects of delivery process. In the simplified tracer scale model we utilized hydraulic conductivity and tracer test data from previous research to characterize the site. We applied IGW (Li et al. 2003;2006) to calibrate a local 2 meter extraction and injection well system. We further employed the calibrated physical parameters to a larger scale [4 sets of extraction injection ] biocurtain scale model. We found that “squeeze effect” or interplay interaction between extraction and injection well system can create a significant difference in the breakthrough curve at the extraction well. We further investigated the efficiency of delivery during recirculation for homogenous and heterogeneous cases. Results reveal that the uncertainty range due to heterogeneity based on different points at biocurtain varies between %40-80 relative concentrations. However results for homogenous case predicts the lowest breakthrough value in the weakest point. Therefore, for design purposes choosing the homophonous value is a conservative approach. In this problem impact of heterogeneity on delivery process is in the favor of design. In the heterogeneous case for a typical control volume we obtained an average value of 121% relative concentration. Results satisfy the goal of 139 design which was delivering more than %50 relative average concentration to the delivery zone and maintains this concentration for 10 to 14 days to the treatment zone [no reaction, pure forward transport]. We further investigated impact of large scale flow dynamic on biocurtain scale delivery. We applied an innovative HPDP tool implemented in IGW to construct and model a multi-scale regional model. HPDP capability allowed us to monitor the delivery process for cases of with and without regional pumping during recirculation and after shutting down the wells. Results reveal that during the recirculation system, given the input physical constant parameters impact of large scale dynamic doesn’t have much impact on biocurtain scale delivery. We need to emphasize that this is not a general conclusion. To better evaluate the impact of large scale dynamic on local scale delivery we may need to run additional variable scenarios to better predict the impact of large scale dynamic on local scale delivery process. We further assessed the impact of large scale pumping on local scale delivery after shutting down the recirculation system. Results reveal that due to irrigation the local gradient increase about %60, which directly impacts the biocurtain retention time which may yield to a different recirculation timing design during the irrigation season. HPDP methodology enhanced our research capability in assessing the impact of nonuniformity of flow and variability of physical parameters on local scales delivery phenomena. 140 ACKNOWLEDGEMENT Jaime A. Graulau-Santiago This research was funded by MDEQ, and NSF. 141 REFERENCE Adams, BE, and L. Gelhar. (1992). Field study of dispersion in a heterogeneous aquifer, 2. Spatial moments analysis. Water Resour. Res. 28, no. 12: 3293-3307. Azadpour-Keeley, A., Keeley, J. W., Russell, H. H., and Sewell, G. W. (2001). "Monitored natural attenuation of contaminants in the subsurface: Processes." Ground Water Monitoring and Remediation, 21(2), 97-107. Boggs, J. M., S. C. Young, L. M. Beard, L. W. Gelhar, K. R. Rehfeldt, and E. E. Adams. (1992). Field-Study of Dispersion in a Heterogeneous Aquifer .1. Overview and Site Description. Water Resources Research 28 (12):3281-329l. Dybas, M. J ., Barcelona, M., Bezborodnikov, S., Davies, S., Fomey, L., Heuer, H., Kawka, 0., Mayotte, T., Sepulveda-Torres, L., Smalla, K., Sneathen, M., Tiedje, J ., Voice, T., Wiggert, D. C., Witt, M. E., and Criddle, C. S. (1998). "Pilot-scale evaluation of bioaugmentation for in-situ remediation of a carbon tetrachloride contaminated aquifer." Environmental Science & T echnologv, 32(22), 3598-3611. Dybas, M. J ., Hyndman, D., Heine, R., Tiedje, J ., Linning, K., Wiggert, D., Voice, T., Zhao, X., Dybas, L., and Criddle, C. (2002). "Development, operation, and long- term performance of a full-scale biocurtain utilizing bioaugmentation." Environmental Science & Technology, 36(16), 3635-3644. Dybas, M. J., Voice, T., Zhao, X., Hashsham Syed, Phanikumar, M.S. , Li, S.G., Long, D., and Tjedje, J. (2003). “Evaluation of Remediation in Schoolcrafi Plumes G and F Contaminated with Chlorinated Solvents and Metals.” Report submitted to Michigan Department of Environmental Quality (MDEQ). Freyberg, D.L. (1986). “A natural gradient experiment on solute transport in a sand aquifer .2. spatial moments and the advection and dispersion of nonreactive tracers”. Advanced Water Resources Research. 22 (13): 2031-2046 Garabedian, S. P., D. R. Leblanc, L. W. Gelhar, and M. A. Celia. (1991). Large-Scale Natural Gradient Tracer Test in Sand and Gravel, Cape-Cod, Massachusetts .2. Analysis of Spatial Moments for a Nonreactive Tracer. Water Resources Research 27 (5):911-924. Graulau, Santiago J. (2003). " Development and application of a methodology to evaluate natural attenuation of chlorainated solvents using conceptual and numerical models" Ph.D. Desertion, Michigan State University, East Lansing, MI. Hess, K. M., S. H. Wolf, and M. A. Celia. (1992). Large-Scale Natural Gradient Tracer Test in Sand and Gravel, Cape-Cod, Massachusetts .3. Hydraulic Conductivity Variability and Calculated Macrodispersivities. Water Resources Research 28 142 (8):2011-2027. Hyndman, D., Dybas, M. J ., Fomey, L., Heine, R., Mayotte, T., Phanikumar, M. S., Tatara, G., Tiedje, J ., Voice, T., Wallace, R., Wiggert, D., Zhao, X., and Criddle, C. S. (2000). "Hydraulic characterization and design of a full-scale biocurtain." Ground Water, 38(3), 462-474. Kim, H., and Dybas, M.J. (in review) “Targeted Biostimulation: Stimulation of High-rate Complete Dechlorination Activity by Pulse Feedings.” Submitted to Applied and Environmental Microbiology, Jan 2006. Leblanc, D. R., S. P. Garabedian, K. M. Hess, L. W. Gelhar, R. D. Quadri, K. G. Stollenwerk, and W. W. Wood. (1991). Large-Scale Natural Gradient Tracer Test in Sand and Gravel, Cape-Cod, Massachusetts .1. Experimental-Design and Observed Tracer Movement. Water Resources Research 27 (5):895-910. Lang, M. M., P. V. Roberts, and L. Semprini (1997). “Model simulations in support of field scale design and operation of oremediation based on cometabolic degradation”, Ground Water, 35(4), 565— 573. Li, S.G and Q. Liu. (2003). “Interactive Ground Water (IGW): An Innovative Digital Laboratory For Groundwater Education and Resear”. Computer Applications in Engineering Education, 1 1 (4): 179-202. Li, S.G and Q. Liu. (2006). “Interactive Ground Water (IGW)”, short communication, Environmental Modelling and Software. 21(3), 41 7-418. Li, S.G., Q. Liu, and S. Afshari.(in press), “An object-oriented hierarchical patch dynamics paradigm (HPDP) for modeling complex groundwater systems across multiple-scales”, Environmental Modelling and Software. Lipinski, B. A. (2002). "Estimating natural attenuation rates for a chlorinated hydrocarbon plume in a glacio-fluvial aquifer, Schoolcrafi, Michigan," M.S. thesis, Michigan State University, East Lansing, MI. Mayotte, T. J ., Dybas, M. J ., and Criddle, C. S. (1996). "Bench-scale evaluation of bioaugmentation to remediate carbon tetrachloride-contaminated aquifer materials." Ground Water, 34(2), 358-367. Mehl, S., and M. C. Hill. (2002). Development and evaluation of a local grid refinement method for block-centered finite-difference groundwater models using shared nodes. Advances in Water Resources 25 (5):497—51 1. Phanikumar, M. S., Hyndman, D. W., Zhao, X., Dybas, M. J. (2005) “A Three- Dimensional Model of Microbial Transport and Bioremediation at the Schoolcrafi, Michigan Site Water Resources Research Vol. 41, W05 01 1 . 143 Rehfeldt, K.R., M.J. Boggs, and L.Gelhar. (1992). Field study of dispersion in a heterogeneous aquifer, 3. Geostatistical analysis of hydraulic conductivity. Water Resour. Res. 28, no. 12: 3281-3291. Scheibe, T. D., Y.-J. Chien, and J. S. Radtke (2001), “Use of quantitative models to design microbial transport experiments in a sandy aquifer”, Ground Water, 39(2), 210— 222. Ward, D. S., D. R. Buss, J. W. Mercer, and S. S. Hughes. (1987). “Evaluation of a Groundwater Corrective Action at the Chem-Dyne Hazardous-Waste Site Using a Telescopic Mesh Refinement Modeling Approach”. Water Resources Research 23 (4):603-617. 144 SUMMARY OF ORIGINAL CONTRIBUTIONS Taking advantage of the innovative hierarchical patch dynamics modeling paradigm [Li et al., in press; Li and Liu, accepted], I have made in this dissertation research a number of original contributions. Specifically, o I have developed a hierarchical modeling approach that enables integrated, model- based groundwater management analysis that provides both a regional solution and detailed local characterizations at critical hotspots; o I have developed a hierarchical modeling approach that enables predicting the clustered drawdown dynamics around wells and the drawdown at the radius of a well in a complex groundwater model under general conditions; and o I have developed a multi-scale groundwater model that can be used to evaluate the complex 3D delivery dynamics and the hydraulic performance of the bioremediation curtain for the plume G at the Schoolcrafi site, Michigan in the presence of large scale spatial heterogeneity and seasonal irrigation pumping in the region. 145 IllilllllllHillllllllllllllllllllllllllllllllllllllllllllll 31293 02845 1148