شبیه­سازی عددی جریان و انتقال آلودگی در آبهای زیرزمینی مطالعه موردی:

 آبخوان دشت نهاوند

حسین بانژاد*¹، حمید محب زاده²، محمد حسین قبادی³ و مجید حیدری⁴

تاریخ دریافت: 16/02/91  تاریخ پذیرش: 15/08/91

^1- دانشیار گروه مهندسی آب، دانشکده کشاورزی،  دانشگاه بوعلی سینا همدان.

^2- کارشناسی ارشد آبیاری و زهکشی، دانشکده کشاورزی، دانشگاه بوعلی سینا همدان.

^3- دانشیار گروه زمین شناسی، دانشکده علوم پایه،  دانشگاه بوعلی سینا همدان.

^4- استادیار گروه مهندسی آب، دانشکده کشاورزی، دانشگاه بوعلی سینا همدان.

^*مسئول مکاتبه: Email: hossein_banejad@yahoo.com

چکیده

شبیه­سازی عددی جریان آب­های زیرزمینی به دلیل تخمین پارامترهای هیدرولیکی و هیدرولوژیکی ابزاری مهم برای مدیریت منابع آب آبخوان­ها می­باشد. این تحقیق نتایج حاصل از یک مدل ریاضی شبیه­سازی جریان آب­های زیرزمینی در آبخوان دشت نهاوند واقع در غرب ایران را نشان می­دهد. برای انجام این کار از کد عددی MODFLOW-2000 که در قالب نرم افزار مدل­سازی آب­های زیرزمینی (GMS) قرار گرفته است، استفاده شد. پس از جمع­آوری اطلاعات مورد نیاز از قبیل اطلاعات زمین شناسی، هیدرولوژیکی، هیدروژئولوژیکی و نقشه­های توپوگرافی ابتدا مدل سه بعدی هیدروژئولوژیکی دشت با استفاده از لاگ چاه­های منطقه و اطلاعات ارتفاعی لایه سطحی دشت تهیه شد و سپس MODFLOW برای شبیه­سازی جریان مورد استفاده قرار گرفت. پس از شبیه­سازی اولیه جریان، مدل با استفاده از سعی و خطا و روش تخمین پارامتر و اطلاعات تراز آب چاه­های منطقه برای حالت پایدار کالیبره گردید. نتایج حاصل از کالیبراسیون نشان می­دهد که خطای میان تراز مشاهده شده و محاسبه شده با توجه به حد مجاز 20± متر در دامنه مطلوب می­باشد همچنین تراز محاسبه شده توسط مدل نشان دهنده حرکت آب­های زیرزمینی در جهت شیب غالب منطقه یعنی از جنوب شرقی به شمال غربی است. درنهایت از کد عددی MODPATH که انتقال توده­ای ذرات را شبیه­سازی می­کند برای تخمین مسیر حرکت آلودگی و همچنین مبدأ آلودگی استفاده شد.

واژه­‎های کلیدی: انتقال آلودگی، جریان آب­های زیرزمینی، GMS، MODFLOW

Numerical Simulation of the Flow and Contaminant Transport in Groundwater, Case Study: Nahavand Plain Aquifer

H Banejad^1*, H Mohebzadeh², AG Ghobadi³ and M Heydari⁴

Received: 5 May 2012   Accepted: 5 November 2012

^1-Assoc. Prof., Dept. of Water Engin., Faculty of Agric, Bu-Ali Sina Univ, Hamedan, Iran.

^2-M.Sc., Dept. of Water Engin., Faculty of Agric, Bu-Ali Sina Univ, Hamedan, Iran.

^3-Assoc. Prof., Dept. of Water Engin., Faculty of Agric, Bu-Ali Sina Univ, Hamedan, Iran.

^4-Assist. Prof., Dept. of Water Engin., Faculty of Agric, Bu-Ali Sina Univ, Hamedan, Iran.

^*Corresponding Author E-mail: hossein_banejad@yahoo.com

Abstract

Numerical simulation of the groundwater flow in order to hydraulic and hydrologic parameters estimation is an important tool for management of aquifers. This study presents results of a mathematical model developed for simulation of groundwater flow in Nahavand plain aquifer in southwest of Hamedan province. For this purpose Groundwater Modeling Software (GMS) was used which supports the MODFLOW-2000 code. After gathering required data same as hydrological, hydrogeological data and topography maps, first 3D hydrogeological model of plain was constructed using boreholes and surface elevation data. Then MODFLOW was used for simulation of flow. After initial simulation of flow, model was calibrated in steady state with trial-and-error and parameter estimation methods and groundwater table monitoring data. Results of calibration show that error between observed head and computed head is in allowable range (±20m). Also results of computed head with model show that groundwater flow is in direct of dominant slope (southeast to northwest). Finally MODPATH code which simulates advective transport of particles was used for estimation of flow path and source of contaminants.

Keyword : Contaminant transport, Groundwater flow, GMS, MODFLOW

مقدمه

مدل­سازی عددی آب­های زیرزمینی ابزاری مهم برای مدیریت منابع آب در آبخوان­ها می­باشد. این مدل­ها می­توانند برای تخمین پارامترهای هیدرولیکی و همچنین مدیریت منابع آب و پیش بینی چگونگی تغییر یک آبخوان در مقابل تغییرات آب و هوایی و پمپاژ استفاده شوند (رجلی و همکاران 2003). از آنجایی که آب­های زیرزمینی در معرض آلودگی­های شدید قرار دارند، استراتژی­های موثر برای مدیریت و حفاظت منابع آب­های زیرزمینی برای اجتناب از اثرات محیطی برگشت ناپذیر مانند کاهش شدید کیفیت این آب­ها و نابودی آنها لازم است (هاردن 2000). از این رو مدل­سازی انتقال و سرانجام آلودگی­های مختلف در آب­های زیرزمینی در سال­های اخیر مورد توجه زیادی قرار گرفته است. در اثر همین افزایش توجه، مدل­های عددی مختلفی برای تحلیل چگونگی حرکت آلودگی در آب­های زیرزمینی توسعه یافته است.

امروزه استفاده از مدل­ها برای به تصویر کشیدن واقعیات و فهم بهتر آن­ها و همچنین اتخاذ تصمیمات درست در مورد هر پدیده­ای بسیار متداول شده است. مدل­های شبیه­سازی جریان آب­های زیرزمینی یکی از این مدل­ها می­باشند که توسط متخصصان در خیلی از کشورها مورد استفاده قرار می­گیرند و دقت و صحت نتایج آن­ها نیز به اثبات رسیده است (هارباگ 2005).

یائوتی و همکاران (2008) از مدل MODFLOW  برای بررسی تغییرات شرایط هیدروژئولوژیکی و شبیه­سازی رفتار جریان، تحت فشارهای مختلف در آبخوان غیر محصور بو،آرگ[1] استفاده کردند و با استفاده از روش سعی و خطا به شبیه­سازی هدایت هیدرولیکی پرداختند و به این نتیجه رسیدند که نوسانات بار هیدرولیکی به تغییرات فصلی تغذیه از بارش و آبیاری بستگی دارد.

لارون و همکاران (2005) تأثیر ساختارهای زمین شناسی را بر روی جریان آب­های زیرزمینی آبخوان صحرای جودین[2] واقع در فلسطین اشغالی مورد ارزیابی قرار دادند و با استفاده از مدل MODFLOW نحوه جریان آب­های زیرزمینی و همچنین نفوذ آن از یک آبخوان به آبخوان دیگر و از یک زیرحوضه به زیرحوضه دیگر را مورد ارزیابی قرار دادند.

بررسی کمیت و کیفیت آب­های زیرزمینی دشت نهاوند به دلیل تأمین بخش عمده­ای از آب کشاورزی و آشامیدنی از اهمیت خاصی برخوردار است . با دانستن کیفیت آب­های زیرزمینی این منطقه و شبیه­سازی انتقال آلودگی­های احتمالی موجود در این آب­ها می­توان به جهت و سرعت انتقال آلودگی پی برده مناطقی را که در سال­های آینده با خطر آلودگی آب­های زیرزمینی مواجه می­شوند را تعیین کرد و برای جلوگیری از انتشار بیشتر آلودگی اقداماتی را انجام داد. همچنین با شبیه­سازی جریان آب­های زیرزمینی می­توان به خصوصیات هیدرولیکی منطقه از قبیل هدایت هیدرولیکی ساختارهای زمین شناسی، تراز سطح آب­های زیرزمینی، مقدار تغذیه آب­های زیرزمینی از طریق بارش و میزان تخلیه از این منابع پی­برد.

مواد و روش ها

اطلاعات مورد استفاده

موقعیت منطقه مورد مطالعه

منطقه مورد مطالعه قسمتی از دشت نهاوند در جنوب استان همدان می­باشد که مساحت km² 285 از کل دشت را به خود اختصاص می­دهد و در موقعیت °48 تا "15 °48 درجه شرقی و "10°34 تا "25°34 درجه شمالی  قرار گرفته است. این دشت یکی از دشت­های حوضه رودخانه کرخه است. شیب عمومی دشت از جنوب شرقی به طرف شمال غربی بوده، مرتفع­ترین و پست­ترین نقاط آن به ترتیب 1650 و 1450 متر از سطح دریا ارتفاع دارند. این دشت از سمت شمال به خط الرأس ارتفاعات چال قیله، آردوشان، شادمانه، کمرزرد، گوین و از سمت جنوب و غرب به کوه گرین محدود می­گردد. شکل 1 مکان دشت مورد مطالعه در استان همدان و خصوصیات ارتفاعی آن را نشان می­دهد.

شکل 1- نقشه دشت نهاوند در استان همدان و موقعیت چاه­های لاگ برداری

 اطلاعات هیدرولوژیکی

در این مطالعه به دلیل کمبود ایستگاه­های باران­سنجی در منطقه از اطلاعات بارش ثبت شده ایستگاه وسج در فواصل سال های 1376 تا 1387 برای انجام محاسبات مربوط به تغذیه آب­های زیرزمینی از طریق بارش استفاده شد. نمودار 1 به ترتیب بارش تجمعی ماهانه برای سال 1387 و میانگین ماهانه آن را برای دوره آماری 1376 تا 1387 نشان می­دهد. مقدار کل بارش در سال 1387 و همچنین مقدار آن از سال 1376 تا 1387، 7/407 و 1/412 میلی­متر در سال می­باشد. نمودار رسم شده برای دوره آماری 1376 تا 1387  نشان می­دهد که بیشترین بارش در ماه­های آبان تا اردیبهشت و کم­ترین آن در ماه­های خرداد تا مهر اتفاق می­افتد.

رودخانه گاماسیاب که بزرگ­ترین رودخانه استان همدان است در این محدوده جریان دارد. این رودخانه پس از دریافت سایر جریانات سطحی دیگر حوضه­های مجاور نظیر ملایر، تویسرکان، اسد آباد و کنگاور در نهایت در محلی به نام دوآب واقع در جنوب غربی استان همدان از منطقه خارج می­گردد.

شکل2- بارش تجمعی ماهانه سال 1387، میانگین بارش تجمعی ماهانه دوره آماری 1376 تا 1387

اطلاعات زمین­شناسی

دشت نهاوند از نظر زمین­شناسی در منطقه زاگرس رورانده یا مرتفع قرار گرفته است. ارتفاعات منطقه عمدتاً متشکل از سنگ­های رسوبی و یا دگرگونی هستند. سنگ­های رسوبی منطقه از جنس سنگ آهک­های کارستی شده وکنگلومرا می­باشند. تشکیلات دگرگونی شامل آهک­های متبلور و شیل­های کم متامورف است که بیشترین گسترش را از خود نشان می­دهند . به طور کلی می­توان گفت که سازندهای شمال شرق نهاوند بیشتر مربوط به قبل از ژوراسیک می­باشند. این سازندها در شمال، سنگ­های آهکی کریستالیزه و طبقات آهکی و شیلی دگرگون شده می­باشند که آثاری از توده ­های آندزیتی کرتاسه نیز در آن­ها به چشم می­خورد. در جنوب غرب نهاوند سازندها عمدتاً آهکی تا آهک کریستالیزه بوده و در سمت جنوب دارای آثار فرعی آتشفشانی می­باشند. رنگ این آهک سفید خاکستری تا خاکستری تیره بوده و در مواردی لایه­های آهکی مارنی تا مارنی در بین آن­ها مشاهده می­گردد. به سبب فعالیت زمین ساختی شدید حاوی درز و شکاف فراوان در جهات مختلف هستند که بعضاً توسط مواد رسی و کلسیت پر شده­اند. بر روی این آهک­ها کنگلومرای پلیستوسن قرار می­گیرد. مقطع زمین­شناسی منطقه مورد مطالعه نشان می­دهد که لایه­هایی از سنگ آهک و سنگ آهک دولومیتی و سنگ­های پیروکلاستیک همراه با سنگ آهک بعنوان سنگ بستر در زیر رسوبات شنی قرار گرفته اند ( شکل 3).

شکل 3- مقطع زمین­شناسی منطقه مورد مطالعه(بی نام، 1991)

با توجه به مقطع زمین­شناسی دشت و موقعیت سازندهای قرار گرفته در زیر آن، ضخامت آبخوان حدود 100 متر می­باشد که فرض می­شود این ضخامت در تمام مساحت دشت همین مقدار است.

اطلاعات هیدروژئولوژیکی

             در این مطالعه از اطلاعات مربوط به تراز آب و همچنین لاگ 19 چاه مشاهده­ای که در سطح منطقه پراکنده­اند استفاده می­شود. متأسفانه اطلاعاتی در مورد دبی برداشتی از این چاه­ها در دست نمی­باشد. عمق این چاه­ها از 15 تا 70 متر متغیر بوده و اطلاعات بدست آمده از لاگ­ها نشان دهنده این امر می­باشد که دشت نهاوند شامل لایه­هایی از رس و سیلت، شن و ماسه، ماسه ریز، شن دانه درشت، آهک و شیست بوده و با توجه به لاگ چاه­ها جنس غالب دشت از لایه­های رس و سیلت ،شن و ماسه می­باشد. به دلیل عدم آزمایش پمپاژ اطلاعاتی در مورد هدایت هیدرولیکی مواد تشکیل دهنده دشت در دسترس نمی­باشد. موقعیت چاه­ها در شکل 1 نشان داده شده است.

ساخت مدل عددی و ورود اطلاعات مورد نیاز به آن

              در این مطالعه از کد عددی MODFLOW که در قالب نرم افزار  [3]GMS 7.1قرار گرفته است برای مطالعه فرایندهای هیدروژئولوژیکی آب­های زیرزمینی دشت نهاوند استفاده شد. MODFLOW  مدلی از نوع تفاصل محدود[4] می­باشد که قادر است به شبیه­سازی سه بعدی جریان در حالت پایدار[5]  و گذرا[6] در آبخوان­های با لایه­های ناهمگن[7]  بپردازد (یائوتی و همکاران 2008).

           در این مطالعه ابتدا اطلاعات هیدرولوژیکی شامل میزان بارش وارده به منطقه، اطلاعات هیدروژئولوژیکی و زمین­شناسی شامل نقشه زمین­شناسی منطقه همراه با مقطع زمین­شناسی آن و موقعیت چاه­های منطقه به همراه لاگ این چاه­ها و در نهایت نقشه­های توپوگرافی، کاربری اراضی و گروه هیدرولوژیکی خاک­های آن برای اجرای مدل از سازمان­های مربوطه (سازمان آب منقطه­ای، سازمان هواشناسی همدان و سازمان زمین­شناسی ایران) جمع آوری شدند. سپس مرز منطقه مورد مطالعه مشخص و مدل مفهومی[8] آن که شامل موقعیت رودخانه گاماسیاب، میزان تغذیه آب­های زیرزمینی و موقعیت چاه­های مشاهده­ای بود، ساخته شد. برای وارد کردن ویژگی­های هیدروژئولوژیکی مانند هدایت هیدرولیکی افقی مواد تشکیل دهنده دشت اقدام به ساخت مدل سه بعدی هیدروژئولوژیکی[9] شد. برای کالیبره کردن مدل مقادیر هدایت هیدرولیکی و تغذیه با استفاده از روش سعی و خطا  و سپس روش  PEST [10] تغییر داده شدند تا نتایج حاصل از مدل با مشاهدات بدست آمده از تراز چاه­های مشاهده­ای بهترین انطباق را داشته باشند.مدل برای سال 1376 کالیبره و میزان آب ورودی و خروجی به آبخوان دشت در این سال محاسبه شد.

تخمین تغذیه حاصل از بارش به آب­های زیرزمینی

            نتایج حاصل از چاه­های مشاهده­ای دشت نهاوند حاکی از این است که آبخوان این دشت آبخوانی آزاد یا نامحصور[11] بوده و به همین دلیل قادر است به طور مستقیم توسط بارش تغذیه شود. در این مطالعه برای محاسبه تغذیه از روش بیلان آب استفاده شد اما بدلیل در دسترس نبودن اطلاعاتی از تبخیر در منطقه از این فاکتور صرفنظر و تغذیه از طریق محاسبه رواناب حاصل از بارش و تفاضل آن از مقدار بارش محاسبه شد.

           برای محاسبه میزان رواناب از روش CNІ که اصلاحی از روش NRCS-CN [12] می­باشد و  توسط پاتیل و همکاران (2008) ارائه شده، استفاده گردیده است.

که در این معادله P: بارش سالانه بر حسب میلی­متر و S: ضریب نگهداشت خاک می­باشد. برای محاسبه رواناب به روش CNІ ابتدا با استفاده از نقشه­های کاربری اراضی منطقه و گروه هیدرولوژیکی خاک های آن نقشه [13]CN با استفاده از نرم افزارHec-GeoHms به صورت رستر[14] تهیه شد. سپس با استفاده از معادله 2 که رابطه­ای را میان ضریب نگهداشت خاک و مقدار CN بر قرار می­کند مقدار ضریب نگهداشت خاک نیز به صورت نقشه رستر تهیه شد.

پس از محاسبه ضریب نگهداشت خاک، با استفاده از معادله 1 که رابطه­ای میان این پارامتر و باران تجمعی برقرار می­کند مقدار رواناب و به تبع آن مقدار تغذیه حاصله از بارش قابل محاسبه خواهد بود.

محاسبه تغذیه برای دوره آماری

برای محاسبه تغذیه در این حالت از اطلاعات بارش کل دوره آماری (1387-1376) استفاده شد. طبق محاسبات میزان بارش تجمعی سالانه برای ایستگاه وسج 1/412 میلی­متر می­باشد. ابتدا با استفاده از این مقدار و با توجه به روش CNІ سعی در محاسبه میزان رواناب حاصل از این بارش شد. پس از انجام محاسبات مقدار رواناب 377-226 میلی­متر در سال و با کم کردن این مقدار از بارش سالانه میزان تغذیه 185-34 میلی­متر در سال یا 0005068/0-0000952/0 متر در روز بدست آمد.

توسعه مدل عددی برای شبیه­سازی جریان آب­های زیرزمینی

مدل جریان در آب­های زیرزمینی با چگالی ثابت و مواد متخلخل با شرایط متعادل توسط معادله دیفرانسیل جزئی (معادله 3) نمایش داده می­شود (مک دونالد و هارباگ 1998).

که در این معادله  ،   و   به ترتیب هدایت هیدرولیکی در جهات x، y و z بوده و دارای بعد ( ) می­باشند. h تراز پیزومتریک(L)  می­باشد.

جریان آب­های زیرزمینی در دشت توسط کد عددی MODFLOW-2000  شبیه­سازی شد. این کد جریان را با استفاده از روش تفاضل محدود حل می­کند. لایه­ها در آن می­توانند به صورت محصور[15] و غیر محصور در نظر گرفته شوند و همچنین قادر است استرس­های وارد بر جریان از طریق چاه­ها، تبخیر، زهکشی و رودخانه­ها را شبیه­سازی کند (هارباگ 2005).

در این مطالعه از MODFLOW که در قالب نرم افزار GMS 7.1 قرار گرفته است استفاده شد. نرم افزار GMS نرم افزاری با واسط گرافیکی[16] است که با داشتن ابزار متفاوت کاربر را قادر می­سازد تا به ساخت مدل مفهومی، شبکه­بندی و حل معادله جریان و در نهایت نمایش نتایج حاصل بپردازد (لای 2002).  شبیه­سازی در این نرم افزار در سه مرحله پیش پردازش[17] ، پردازش[18] و پس پردازش[19] صورت می­گیرد.

          مدل سه­بعدی هیدروژئولوژیکی

برای ساخت مدل سه­بعدی هیدروژئولوژیکی دشت از اطلاعات لاگ 19 چاه موجود استفاده شد. بدلیل عدم انطباق اطلاعات چاه­ها با یکدیگر، سعی در ساده­سازی اطلاعات شد و لایه­هایی که در اکثر چاه­ها وجود داشتند به عنوان لایه­های تشکیل دهنده آبخوان در نظر گرفته شدند. همانطوری که قبلا مطرح گردید این لایه­ها رس و سیلت، شن و ماسه می­باشند.

در اطلاعات لاگ چاه­ها محل سنگ بستر مشخص نشده بود به همین دلیل برای تعیین محل سنگ بستر از مقطع زمین­شناسی دشت استفاده شد. طبق این مقطع لایه­های زیرین دشت به طور عمده از جنس آهک کریستالیزه می­باشند که تقریباً در عمق 100 متری قرار گرفته و دارای نفوذ پذیری پایینی هستند به همین دلیل این لایه­ها به عنوان سنگ بستر در نظر گرفته شدند و مدل هیدروژئولوژیکی با ضخامت تقریبی 100 متر بر روی سنگ بستر ایجاد شد.

ابتدا اطلاعات لاگ چاه­ها وارد GMS شد و پس از تعریف لایه­های مختلف دشت مدل سه بعدی زمین شناسی آن تهیه گردید. پس از آن شبکه سه­بعدی [20]اولیه با توجه به اطلاعات موجود و شرایط هیدروژئولوژیکی در ابعاد 150 و 100 سلول در جهت x و y ساخته شد این شبکه دارای 91506 گره[21]  و 75000 سلول می­باشد. مدل در جهت عمودی به 5 لایه تقسیم گردید تا نمایش گرافیکی نتایج بهتر صورت گیرد. شکل 4 (الف، ب وج) به ترتیب مدل سه­بعدی زمین­شناسی، شبکه سه­بعدی اولیه و محدوده منطقه مورد مطالعه را نشان می­دهد. در ادامه شبکه سه­بعدی اولیه محدود به مرز منطقه شد (شکل4-د) و مدل سه­بعدی زمین­شناسی به آن انتقال پیدا کرد و در نهایت لایه سطحی آن توسط اطلاعات ارتفاعی منطقه و از طریق روش وزن­دهی معکوس فاصله[22] تصحیح شد (شکل 4- ه).

مدل مفهومی

دقت مدل هیدروژئولوژیکی وابسته به صحت تعریف شرایط فیزیکی منطقه و تخمین صحیح آب ورودی و خروجی به آن است (لارون و همکاران 2005). مدل مفهومی شامل شرایط هیدرولیکی، تغذیه و تخلیه و محدوده منطقه مورد مطالعه می­باشد همچنین می­توان با استفاده از آن خصوصیات هیدرولوژیکی منطقه مانند محل چاه­ها و رودخانه­ها را نیز تعریف کرد.

مدل مفهومی در این مطالعه جهت ورود اطلاعات به مدل هیدروژئولوژیکی در چند پوشش[23] مجزا ساخته شد. پوشش اول مربوط به میزان تغذیه به آب­های زیرزمینی می­باشد که این مقدار از محاسبات روش NRCS-CN که پیشتر ذکر شد به عنوان ورودی برای تغذیه وارد شد. مقدار ورودی میانگین حد پایین و بالای محاسبات یعنی 000301/0 متر در روز در نظر گرفته شد و فرض شد این مقدار به طور یکنواخت در کل منطقه اتفاق می­افتد.

طبق مشاهدات انجام شده از منطقه، رودخانه گاماسیاب در اکثر نقاط توسط آب­های زیرزمینی تغذیه می­گردد. (شکل 5) در نتیجه این رودخانه در مدل مفهومی به صورت زهکش[24] با قابلیت انتقال[25] مشخص در نظر گرفته شد.  پوشش سوم مربوط به تراز آب 18 چاه مشاهده­ای موجود در منطقه می­باشد که لازم است اطلاعات تراز آب این چاه­ها نیز برای استفاده از آن­ها در کالیبراسیون مدل و مقایسه تراز مشاهده­ای و محاسبه­ای وارد مدل مفهومی شود. محدوده هر سه پوشش در لایه اول تعریف و سپس مدل مفهومی وارد مدل سه­بعدی هیدروژئولوژیکی شد.

شرایط مرزی و شرایط اولیه

خصوصیات هیدرولوژیکی مجاور و درون محدوده مدل باید توسط شرایط مرزی نمایش داده شوند. خطوط مرزی دشت به صورت مرزهای بدون جریان[26] در نظر گرفته شدند در حالی که محل ورود رودخانه به مدل به صورت مرز با تراز ثابت[27] در نظر گرفته شد. تراز اولیه آب­های زیرزمینی برای اجرای اولیه مدل با استفاده از توپوگرافی منطقه تعریف شد زیرا طبق یافته های هابرت (1940) تراز آب از سطح زمین تبعیت می­کند شکل 6 نحوه تعریف شرایط مرزی را برای مدل نشان می­دهد.

بعد از تعریف شرایط مرزی با استفاده از بسته [28]LPF که در MODFLOW برای حل جریان و تراز آب­های زیرزمینی قرار داده شده است، استفاده شد. این بسته به صورت مجزا هر ساختار هیدرولوژیکی را شبیه­سازی کرده و به محاسبه جریان در هر سلول می­پردازد (یانگ و همکاران 2009). مقادیر اولیه هدایت هیدرولیکی افقی طبق اعداد ارائه شده توسط تاد (2005) برای لایه­های رس و سیلت، شن و ماسه به ترتیب 1 و50 متر در روز در نظر گرفته شد.

شکل4- مراحل ساخت مدل سه­بعدی هیدروژئولوژیکی- الف. مدل سه­بعدی زمین­شناسی، ب. شبکه سه­بعدی اولیه، ج. مرز منطقه مورد مطالعه، د. محدود شدن شبکه سه­بعدی به مرز منطقه، ه. وارد شدن نقشه زمین­شناسی سه­بعدی به شبکه سه­بعدی و تصحیح لایه سطحی با اطلاعات ارتفاعی

شکل 5- رودخانه گاماسیاب، نمایی از زهکشی آب­های زیرزمینی توسط رودخانه

نتایج و بحث

 کالیبراسیون

به طور متداول مدل­ها توسط فرآیند سعی و خطا که در آن پارامترهای مدل در محدوده معقولی از یک شبیه­سازی به شبیه­سازی دیگر تعدیل می­شوند تا بهترین برازش برای مدل بدست آید، کالیبره می­شوند (یائوتی و همکاران  2008).

مدل­ها توسط روش­های معکوس[29] نیز می­توانند کالیبره شوند. در این روش­ها مقادیر بهینه یک پارامتر برای ساختار پارامتر داده شده با استفاده از تکنیک­های ریاضی مانند رگرسیون غیر خطی[30] تعیین می­شوند(کولی و نف 1990). این تکنیک اغلب به عنوان تخمین پارامتر نامیده می­شود. در این مطالعه از برنامه تخمین پارامتر به دلیل توانایی آن در محدود کردن دامنه مقادیر استفاده شد و عمل کالیبراسیون برای سال 1376 انجام شد.

کالیبراسیون برای سال 1376

مدل برای حالت پایدار با فرض بارش ثابت، تخلیه و تراز ثابت آب­های زیرزمینی کالیبره شد. در این مطالعه فرض شد که هدایت هیدرولیکی افقی و تغذیه از پارامترهای ناشناخته جریان بوده و روش PEST برای تخمین آن­ها به کار گرفته شد. عمل کالیبراسیون برای سال 1376( شروع دوره آماری) انجام شد و مشاهدات استفاده شده برای انجام عمل کالیبراسیون تراز آب چاه­های منطقه در سال 1376 بودند.

 بعد از انجام کالیبراسیون مقدار تغذیه 0008211/0متر در روز بدست آمد که این مقدار 72/2 برابر مقدار تغذیه ورودی برای مدل می­باشد. مقدار هدایت هیدرولیکی افقی لایه­های رس، سیلت و شن، ماسه نیز به ترتیب 5 و 52/26 متر در روز بدست آمد. پس از محاسبه هدایت هیدرولیکی و تغذیه با استفاده از روش سعی و خطا، سعی در محاسبه میزان هدایت رودخانه گاماسیاب برای رسیدن به بهترین برازش میان محاسبه و مشاهده شد. مقدار این پارامتر بعد از سعی و خطا 5 بدست آمد. شکل 7 نموداررسم شده برای مقایسه میان تراز مشاهده­ای و محاسبه­ای را نشان می­دهد. در این شکل شاخص کالیبراسیون بیانگر چگونگی انجام عمل کالیبراسیون می­باشد. این شاخص دارای یک حد  پایین و بالا است که میزان حد مجاز برای تخمین بیش از حد[31] و تخمین کم تر از حد[32] را نشان می­دهد که در اینجا 20±  متر در نظر گرفته شده است. قسمت میانی این شاخص میزان خطای موجود بین محاسبه و مشاهده را نشان می­دهد. خطای قرار گرفته در دامنه مجاز به رنگ سبز می­باشد (یانگ و همکاران 2009).

نحوه کالیبراسیون با پارامترهای آماری ریشه مربع میانگین خطا[33] ، میانگین خطای مطلق[34] و میانگین خطا[35] بررسی شد که این پارامترها به ترتیب 44/10، 1/9 و 86/2 می­باشند. با توجه به مقدار میانگین خطا مدل تمایل به تخمین بیش از حد تراز آب دارد.

میزان آب ورودی و خروجی به آبخوان برای سال 1376 در جدول 2 آمده است. میزان جریان بر حسب متر مکعب در روز بیان می­شود و مقادیر مثبت و منفی به ترتیب بیانگر آب ورودی و خروجی به آبخوان هستند. تعادل میان آب ورودی و خروجی با فرض حالت پایدار مطابقت دارد.

شکل8 (الف وب) به ترتیب تراز آب محاسبه شده توسط مدل برای لایه 1و2 را نشان می­دهد. مناطق غیر اشباع[36] در لایه 1 به صورت سلول­های خشک[37] تبدیل شده اند. همانطور که در شکل دیده می­شود جریان آب از شیب عمومی دشت تبعیت می­کند یعنی از جنوب شرقی به شمال غربی در حال حرکت است.

شکل 6- نحوه تعریف شرایط مرزی برای مدل

شکل 7- مقایسه میان ترازهای مشاهده شده و محاسبه شده

جدول 1- مقادیر آب ورودی و خروجی به آبخوان در سال 1376(m³/day)

شکل 8- تراز آب محاسبه شده توسط مدل برای لایه 1و 2 در سال 1376 (سلول­های خشک با رنگ سفید

نمایش داده شده­اند)

مقادیر تراز آب در جنوب شرقی و شمال غربی به ترتیب 1564 و 1475 متر می­باشد که با توجه به فاصله تقریبی 28 کیلومتری باعث ایجاد گرادیان هیدرولیکی 3/0% در دشت می­شود که این امر خود نشانه سرعت کم جریان آب­های زیرزمینی می­باشد.

شبیه­سازی انتقال آلودگی توسط ردیابی ذرات

در روش ردیابی ذرات گروهی از ذرات فرضی در مکانی خاص قرار داده شده و جهت و زمان حرکت آن­ها بر اساس سرعت جریان آب مشخص می­شود. با مدل­سازی حرکت ذرات می­توان حرکت رو به جلو و همچنین رو به عقب آن­ها را جهت مشخص کردن مبدأ و مقصد آلودگی شبیه­سازی کرد (بری و همکاران 2009). در این مطالعه برای انجام این کار از کد عددی MODPATH که یک برنامه پس­پردازش ردیابی ذرات^[38] است و توسط پولاک (1994) ارائه شده است جهت تعیین مسیر و زمان انتقال ذرات استفاده شد.

ذرات برای بررسی حرکت رو به جلوی آلودگی در دو طرف دشت قرار داده شده و انتقال آلودگی توسط ردیابی حرکت آن­ها شبیه­سازی شد (شکل 9). به دلیل حرکت کند ذرات در آب­های زیرزمینی و عدم تأثیرتناوب ماه­های تر و خشک در حرکت ذرات این شبیه­سازی برای تراز محاسبه شده توسط مدل در سال 1376 انجام شد. زمان انتقال بین 30 تا 1000 سال بود و شبیه­سازی در لایه دوم انجام شد. همانطوری که در شکل 9 دیده می­شود نحوه حرکت ذرات منطبق بر گرادیان آب­های زیرزمینی بوده و آلودگی از هر دو طرف دشت به سمت مرکز دشت حرکت می­کند. دلیل این امر تفاوت ارتفاعی در مرزها و مرکز می­باشد که باعث ایجاد گرادیان و هدایت ذرات به مرکز دشت می­شود. طولانی­ترین مسیر طی شده ذرات در این بازه زمانی 12232 متر می­باشد.

برای بررسی مبدأ آلودگی از شبیه­سازی رو به عقب حرکت ذرات استفاده شد. برای این کار 24 ذره در لایه دوم و در مرکز 5 چاه انتخابی در سطح دشت قرار گرفت و حرکت افقی و رو به عقب ذرات برای 10، 20 و 100 سال شبیه­سازی شد (شکل 10). این شبیه­سازی نیز برای حالت پایدار سال 1376 انجام شد زیرا تغییرات ماهانه تراز و گرادیان آب تأثیری بر روی حرکت افقی ذرات که در طی ده ها سال انجام می­شود، ندارد. در طول مدت شبیه­سازی بیشترین مسیر طی شده توسط ذرات از مرکز هرچاه به سمت مبدأ آلودگی تقریباً 8366 متر می­باشد.

طبق یافته­های یائوتی و همکاران )2008( کلرید (Cl^-) بهترین یون برای انجام عمل ردیابی مسیر ذرات می­باشد زیرا این یون در طول مسیر حرکت هیچ گونه واکنشی با عناصر اطراف خود ایجاد نمی­کند. نتایج حاصل از ردیابی ذرات به صورت رو به عقب نیز نشان دهنده این امر می­باشد که در صورت بروز آلودگی در چاه­های انتخابی منشأ این آلودگی از نقاط اطراف دشت می­باشد که پس از طی یک مسیر طولانی  وارد چاه­های قرار گرفته در مرکز آن می­شود. به این ترتیب با دانستن مسیر حرکت آلودگی و شناسایی مکان­های آلوده  می­توان با اقدامات مدیریتی از انتقال آن به نقاط دیگر دشت جلوگیری کرد.

شکل9- شبیه­سازی مسیر جریان ذرات انتخابی قرار گرفته در دو طرف دشت به صورت حرکت رو به جلو برای شرایط پایدار (1376)-زمان حرکت 30 سال (الف)، 100 سال (ب) و 1000 سال (ج)

شکل10- شبیه­سازی مسیر جریان ذرات  قرار گرفته در مرکز 5 چاه انتخابی به صورت حرکت رو به عقب برای شرایط پایدار (1376)-زمان حرکت 10 سال (الف)، 20 سال (ب) و 100 سال (ج)

نتیجه­گیری کلی

آبخوان غیر محصور دشت نهاوند با استفاده از کد عددی MODFLOW برای شبیه­سازی سه­بعدی جریان آب­های زیرزمینی در حالت شرایط پایدار مورد ارزیابی قرار گرفت. برای ساخت مدل ریاضی پارامترهای مختلفی جمع آوری شدند. استفاده از مدل مفهومی کمک زیادی در تعریف شرایط مرزی و ویژگی های فیزیکی آبخوان کرد. نتایج حاصل از کالیبراسیون مدل نتایج معقولی میان تراز مشاهده شده و محاسبه شده در سال 1376 نشان داد. مقادیر ریشه مربع میانگین خطا، میانگین خطای مطلق و میانگین خطا به ترتیب 44/10 ، 1/9 و 86/2 بود. مقدار محاسبه شده تغذیه آب­های زیرزمینی توسط بارش 0008211/0 متر در روز بود که این مقدار 72/2 برابر مقدار اولیه یعنی 000301/0 متر در روز می­باشد. مقادیر هدایت هیدرولیکی لایه­های رس،سیلت و شن، ماسه نیز به ترتیب 5 و 52/26 متر در روز محاسبه شد. طبق محاسبه تراز آب برای مدل، جریان آب­های زیرزمینی از جنوب شرقی به شمال غربی است که این جریان با شیب غالب منطقه مطابقت دارد. مقادیر تراز آب در جنوب شرقی و شمال غربی به ترتیب 1564 و 1475 متر می­باشد که با توجه به فاصله تقریبی 28 کیلومتری باعث ایجاد گرادیان هیدرولیکی 3/0% در دشت می­شود.

نتایج حاصل از MODPATH نشان داد که حرکت آلودگی از دو طرف دشت به سمت مرکز دشت می­باشد اما این حرکت، حرکتی کند است (1000 سال). شبیه­سازی رو به عقب ذرات از مرکز چاه­های انتخابی نیز نشان دهنده این موضوع است که با استفاده ازاین مدل می­توان نحوه ورود آلودگی به چاه­های منطقه را از نظر مکانی و زمانی مورد ارزیابی قرار داد.

سپاسگذاری

در اینجا لازم است تا از سازمان آب منطقه­ای همدان، سازمان زمین­شناسی ایران و سازمان هواشناسی همدان به خاطر در اختیار قرار دادن اطلاعات تراز آب­های زیرزمینی، بارش، نقشه­های زمین­شناسی منطقه و... تشکر نمائیم.