Математическое моделирование миграции радионуклидов в поверхностном слое грунта Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

УДК 519.6

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ МИГРАЦИИ РАДИОНУКЛИДОВ В ПОВЕРХНОСТНОМ СЛОЕ ГРУНТА

А.Д. Истомин, С.А. Кораблева, М.Д. Носков

Северский государственный технологический институт E-mail: istomin@ssti.ru

Представлена математическая модель, описывающая основные физико-химические процессы, происходящие в системе газовая фаза - водная фаза - грунт. Модель описывает многокомпонентную фильтрацию водной фазы, движение газовой фазы, испарение (конденсацию) и кристаллизацию (плавление), диффузию, сорбцию и десорбцию радионуклидов, теплопроводность и конвективный теплообмен. На основе модели было создано проблемно-ориентированное программное обеспечение, которое позволяет прогнозировать миграцию радионуклидов в поверхностном слое почвы. Приведены результаты моделирования вертикальной миграции радиоактивных веществ в поверхностном слое почвы.

Изучение поведения радионуклидов в почвенном покрове является важной задачей с точки зрения экологической безопасности регионов, в которых расположены предприятия атомной промышленности. Радионуклиды могут попасть в почву в результате внештатных ситуаций, связанных с атмосферными выбросами и проливами радиоактивных веществ. Кроме этого в 50-70 годы XX в. для хранения жидких радиоактивных отходов (ЖРО), образующихся в процессе работы предприятий ядерного топливно-энергетического цикла, использовались поверхностные хранилища представляющие собой естественные (оз. Карачай) или специально созданные гидроизолированные водоемы [1, 2]. Миграция радионуклидов из естественных и поверхностных водоемов представляет собой потенциальную угрозу для экосистемы.

Миграция радионуклидов в почвенном покрове определяется значительным количеством взаимосвязанных физико-химических процессов, происходящих в системе: газовая фаза - водный раствор - вмещающая порода. Вследствие сложности этих процессов для описания поведения радионуклидов целесообразно использовать методы математического моделирования. С помощью математического моделирования возможно описать основные процессы, происходящие в системе, с учетом их нелинейного влияния друг на друга.

В настоящей работе представлена математическая модель, описывающая вертикальную миграцию радионуклидов в ненасыщенной зоне почвы с учетом фильтрации воды с поверхности почвы, фильтрации воздуха в почве, испарения (конденсации) и кристаллизации (плавления) воды, диффузии, сорбции и десорбции радионуклидов, радиоактивного распада, энерговыделения в результате радиоактивного распада, конвективного и кондук-тивного теплообмена. Кроме этого, учитывается внешнее воздействие со стороны окружающей среды (изменение атмосферного давления, температуры и влажности воздуха на поверхности почвы, а также атмосферные осадки). Приводятся результаты моделирования эволюции состояния поверхностного слоя почвы и вертикальной миграции радионуклидов.

В модели рассматривается одномерная вертикальная миграция радионуклидов в поверхностном слое почвы, который представляет собой многофазную систему жидкость - газ - твердая фаза. Рассматриваемая система состоит из семи частей: вмещающая порода (0), жидкая фаза (1), лед (2), газовая фаза (3), граница раздела фаз (ГРФ) порода - жидкая фаза (4), ГРФ порода - лед (5) и ГРФ порода - газовая фаза (6). Для количественного определения доли порового пространства, занимаемого частью системы Ф, используется понятие насыщенности £Ф:

V

г. _ е ф ^ф - —,

m V

(Ф- 1 - 3),

где т - пористость; VФ - объем занимаемый фазой Ф; V - объем среды. Части системы с 1 по 3 занимают весь объем порового пространства:

1 = £ ¿ф .

ф=1

Рассматривается поведение воды и радионуклидов (индекс воды равен 1, радионуклиды нумеруются с 2 по И). Полагается, что радионуклиды не переходят в газовую фазу при испарении воды. При замерзании воды радионуклиды остаются на поверхности вмещающей породы, контактирующей со льдом. Таким образом, радионуклиды могут находиться только в жидкой фазе и в адсорбированном состоянии на поверхности породы. В жидкой фазе количество /-ого радионуклида описывается объемной активностью а{:

А = А

a - — -

V пЯу

где А/Ф - активность /-ого радионуклида, находящегося в части системы Ф. При определении объемной активности /-ого радионуклида на поверхности породы а 'Ф+3, контактирующей с частью системы Ф= 1,2,3 полагается, что доля поверхности, с которой контактирует часть системы Ф, пропорциональна ее насыщенности:

A'

(1 - m)VS ф

(ф- 4 - 6).

Количество воды в газовой фазе описывается с помощью объемной концентрации С31:

«ф -

Сз _ V

М1

где Мф - масса /-ого компонента, находящегося в части системы Ф.

Вертикальная фильтрация жидкой фазы определяется капиллярными и гравитационными силами. Значение скорости фильтрации и1 жидкой фазы в ненасыщенной зоне рассчитывается с помощью уравнения Клюта-Ричардса [3]:

и'_-к£1г р -к )■

где к - абсолютная проницаемость; Д^) - функция фазовой проницаемости жидкой фазы; /л - вязкость воды; g - ускорение свободного падения; Т - потенциал капиллярного давления. В настоящей работе для определения зависимости фазовой проницаемости /И) от насыщенности используется степенная функция. Кроме этого полагается, что движение сплошной водной фазы может происходить при условии, что ее насыщенность превышает некоторое критическое значение которое называется остаточной насыщенностью [4]:

/ №) _

1 - 5*

0,

51 > 51

51 <

где у- параметр модели. В работе используется линейная зависимость потенциала капиллярного давления от свободной влажности [3, 5]:

/ 5 - 5* ^

Т_ Ь -1

т - Б*

из( 7) _-Ц +\\—Р-

•1,3 , •'зКТ

М 'Р

тБ3 дР

Т дг Р дг

тБ3 дТ

Л

dz.

Р1 _ ещ

М1 '

Плотность потока воды /¡,2 между жидкой фазой и льдом рассчитывается из баланса тепловой энергии.

Молекулярная диффузия радионуклидов в жидкой фазе и водяного пара в газовой фазе описывается на основе закона Фика. Диффузионный поток остается пропорциональным градиенту концентрации. Для учета замедления диффузии вследствие изогнутости пор и наличия в системе других фаз в законе Фика используется эффективный коэффициент диффузии /-ого компонента в части системы Ф Я ф Эф, связанный с коэффициентом молекулярной диффузии в фазе Ф Яф соотношением [6]:

Д

ф, Эф

Дф тБф X

где Ь, £ - параметры потенциала капиллярного давления.

Скорость вертикальной фильтрации газовой фазы из определяется из закона сохранения объема, записанного с учетом несжимаемости жидкой и твердой фаз:

где х - извилистость.

Плотность потока /,4 /-ого радионуклида между жидкой фазой и поверхностью вмещающей породы в результате адсорбции/десорбции рассчитывается в рамках приближения линейной кинетики, согласно которой /1/,4 пропорциональна отклонению от равновесного значения, определяемого с помощью изотермы Генри:

•,4 _У[,ААа - КА а,).

где К14 - коэффициент распределения для /-ого компонента между жидкой фазой и поверхностью вмещающей породы; 71,4 - постоянная скорости массообмена между жидкой фазой и поверхностью вмещающей породы для /-ого компонента.

Изменение насыщенностей частей системы с 1 по 3, будет сопровождаться переходом части ГРФ от одной фазы к другой и изменением объемных активностей радионуклидов в них. Для описания этого процесса вводятся эффективные массопото-ки /4/,5 и /4/,6. Значение потока /4/,5 пропорционально изменению насыщенности льда:

( ( яе Л ( ас Л А

, _д 52 •4'5 _ дг

в

д Б2 д г

а +в|

д г

где Т - температура; Р - давление; р/ - плотность /-ого компонента; /'ФТ - плотность потока /-ого компонента из части Ф в часть Т; / - время; Я - универсальная газовая постоянная; М/ - молярная масса /-ого компонента.

Плотность потока воды /1,3 между жидкой и газовой фазами в результате испарения (конденсации) пропорциональна отклонению парциального давления водяного пара в газовой фазе Р1 от равновесного значения Р ¿(Т):

•з _г!з ^1з(Р1(Т) - Р1),

где 5ФТ - удельная площадь контакта частей системы Ф и Т. Используя уравнения состояния Менде-леева-Клайперона для водяного пара, можно найти величину парциального давления Р1:

а /4,6 - изменению насыщенности газовой фазы:

, _д Бз

•4'6 _ дг

в

д Б3 д г

+ в|

д г

где в(х) - ступенчатая функция (в(х)=0 при х>0 и в(х)=0 при х<0).

Расчет динамики температурного поля основывается на предположении о локальном тепловом равновесии между всеми частями системы. Объемная теплоемкость среды определяется ее составом:

з

с _ Со(1 - т) + т^ Бфсф,

ф_1

где сФ - объемная теплоемкость части системы Ф. Теплопроводность в системы определяется уравнением Максвелла [7]:

3Д, - 2 т [в + £ ^фвф _V_ф=_

3 + т \Р0/X 5 фрф-1

Р=-

где вФ - теплопроводность части системы Ф.

Энерговыделение в системе происходит в результате радиоактивного распада и фазовых переходов. Мгновенное значение удельной мощности V энерговыделения в результате распада ''-ого радиоактивного радионуклида рассчитывается по формуле:

Г' =Х ™'аф,

ф

где V - энергетическая постоянная для ''-ого компонента. Удельная мощность энерговыделения V в результате испарения/конденсации и кристаллизации/плавления воды определяется формулой:

г1 = АХ - ьисЛ,

где Ьш и Ьисп - скрытые теплоты кристаллизации/плавления и испарения/конденсации воды.

Таким образом, изменение насыщенностей частей системы определяется следующей системой дифференциальных уравнений:

д/ "

д 52

т-

д Т - ¿1^ - ¿1 ' 1 и1 1 д 2 Р Р

т

д/

т

д 53 д/

= -А из + д 2

¿1,

р д/ д2 р

Изменение объемной активности ''-ого радионуклида рассчитывается из закона сохранения массы:

т

д5>; = д

(

д/

д 2

да4

-а и + д, эф

д ах

л

- ¿1,4 -х а1,

= ¿1,4 -Х а4 - ¿4,5 - ¿4,6.

да[ . 1 1 ,,■ т—1 = -Ха\ +¿4.

д/ 5 4,~

т

да6 д/

где X - постоянная распада ''-ого радионуклида. Динамика концентрации воды в газовой фазе определяется уравнением:

д53С3

г-3—3

д/

д д 2

-С3 Т + В1 ф

д 2

Л

+ ¿1,4-

Уравнение для расчета динамики температурного поля системы, записанное с учетом конвективного теплообмена, теплопроводности и энерговыделения в результате радиоактивного распада и фазовых переходов, имеет вид:

д (сТ)

= - —I (с1Т1 +с3Т3)Т д / д 2 { 3 д 2) £

Для прогнозирования миграции радионуклидов в поверхностном слое почвы необходимо учитывать внешнее воздействие со стороны окружающей среды. К основным факторам внешнего воздействия

, д Т

относят: изменение атмосферного давления, температуры и влажности воздуха, атмосферные осадки. Естественные факторы (давление, температура, влажность) подвержены периодическим колебаниям (суточные, годовые циклы), а также случайным изменениям в соответствии с меняющимися погодными условиями. В настоящей модели для определения атмосферного давления, температуры поверхности земли и влажности атмосферного воздуха используются периодические граничные условия, описывающие суточные и сезонные колебания этих величин и позволяющие учесть внешнее воздействие на рассматриваемую систему. Суммарное количество осадков, выпавшее на поверхность почвы за один месяц, соответствовало данным, представленным в работе [8]. При этом интенсивность, продолжительность и время начала осадков задавались произвольным образом и соответствовало средней полосе России. В модели учитывается увеличение инфильтрации воды с поверхности весной в результате таяния снега.

На основе численной реализации представленной математической модели было разработано проблемно-ориентированное программное обеспечение, позволяющее моделировать вертикальную миграцию радионуклидов в поверхностном слое почвы. С помощью программного обеспечения проведено моделирование миграции $г90 и С^37 с поверхности почвы. Исследовались динамика температурного поля с учетом радиогенного энерговыделения, влияние различных физико-химических процессов (инфильтрации воды с поверхности, взаимодействие с вмещающей породой) на вертикальную миграцию радионуклидов. В начальный момент времени радионуклиды находились в адсорбированном состоянии на поверхности почвы, в слое толщиной 5 мм. Суммарная активность радионуклидов на единицу площади составляла 0,5 мКи/м2. Коэффициенты распределения между жидкой фазой и поверхностью породы - ^4=1,35, ^4=40,5. Пористость почвы от=0,4. Остальные параметры моделирования определялись по данным [8, 9].

На динамику температурного поля основное влияние оказывают годовые и суточные изменения температуры поверхности почвы. Суточные и годовые колебания температуры поверхности приводят к формированию температурных волн, распространяющихся вглубь почвы. Прохождение температурной волны сопровождается колебаниями температуры грунта. По мере увеличения глубины температурные волны затухают. Интенсивность затухания тем больше (соответственно глубина проникновения волн тем меньше), чем меньше период колебаний температуры на поверхности. Суточные колебания температуры приводят к изменению температуры только в приповерхностном слое грунта. Амплитуда суточных колебаний температуры уменьшается с глубиной и становится пренебрежимо мала на глубине 30 см (рис. 1). Амплитуда годовых колебаний также уменьшается с глубиной. Поскольку период годовых колебаний больше периода суточных колебаний, то глубина про-

никновения температурных волн больше и составляет 1,5 м (рис. 2). Полученные результаты моделирования динамики температурного поля в поверхностном слое почвы согласуются с экспериментальными наблюдениями [8]. Моделирование с различной начальной активностью радионуклидов показало, что радиогенное тепловыделение приводит к заметному повышению температуры в поверхностном слое почвы лишь при значительных активностях (более 103 Ки/м2).

20 30

Время, час

Рис. 1. Изменение температуры на различной глубине в течение двух суток (июль): 1) на поверхности, 2) на глубине 5, 3) 10, 4) 15, 5) 20 см

30 П

Время, месяц

Рис. 2. Годовые колебания температуры на различной глубине (в 12 ч дня): 1) на поверхности, 2) 20, 3) 50, 4) 100 см

ч

Л- 0,6

и

Й 0,4

¡3 0,2

Рис. 3.

«на, м

N

2-

л н о о м я

Глубина, м

Рис. 4. Объемная активность стронция в жидкой фазе: 1) 1 год, 2) 3 года, 3) 5 лет, 4) 7лет, 5) 10 лет

я

о я я я

£

0,2

0,0

2

Миграция стронция с конвективным потоком жидкой фазы: 1 -4 - насыщенность водной фазыI, 5-8 - объемная активность стронция в водной фазе; 1, 5)до выпадения осадков; 2, 6) через 3 ч, 3, 7) через 40 ч, 4, 8) через 70 ч после начала выпадения осадков

Глубина, м

Рис. 5. Объемная активность цезия в жидкой фазе: 1) 1 год, 2) 10 лет, 3) 20 лет, 4) 30 лет

Вертикальная миграция радионуклидов происходит в результате молекулярной диффузии и конвективного движения жидкой фазы. Отсутствие жидкой фазы в пористой среде в результате ее замерзания или испарения приводит к прекращению миграции радионуклидов. Вследствие малой скорости диффузии радионуклидов в поверхностном слое почвы основным механизмом вертикальной миграции является массоперенос с нисходящим потоком жидкой фазы. На рис. 3 приведены результаты моделирования миграции стронция в результате инфильтрации воды с поверхности почвы после выпадения осадков. До выпадения осадков насыщенность воды в почве определяется капиллярными силами и составляет 10 % от объема пор (кривая 1). Объемная активность стронция в жидкой фазе имеет форму вала, сформировавшегося в результате миграции с водой в течение трех лет после попадания стронция на поверхность почвы (кривая 5). После выпадения осадков, в почве формируется область с повышенной насыщенностью водной фазы. Наибольшая насыщенность определяется интенсивностью осадков, а размеры области с повышенной насыщенностью - суммарным объемом выпавших осадков (кривая 2). Увеличение насыщенности водной фазы в области залегания радионуклидов вследствие десорбции приводит к увеличению их концентрации в жидкости (кривая 6). При движении водной фазы под действием гравитационных сил в результате сорбции происхо-

3

0

1.0-

0.8-

0,6-

0,4-

1,0

0,8

0,0

дит отставание фронта активности (кривые 7, 8) от фронта насыщенности (кривые 3, 4). Чем больше коэффициент распределения радионуклида между жидкой и твердой фазами, тем меньше его концентрация в жидкой фазе и соответственно тем медленнее происходит его вертикальная миграция. На рис. 4 и 5 представлены вертикальные распределения объемной активности стронция и цезия в различные моменты времени. Большой коэффициент распределения стронция между жидкой фазой и породой приводит к замедлению его скорости миграции. Таким образом, скорость вертикальной мигра-

ции радионуклидов определяется количеством ин-фильтровавшейся с поверхности воды и коэффициентом распределения с породой.

Разработанная математическая модель адекватно описывают вертикальную миграцию радионуклидов в поверхностном, ненасыщенном слое почвы. Созданное на основе численной реализации модели проблемно-ориентированное программное обеспечение может быть использовано для прогнозирования распространения радиоактивных примесей в поверхностном слое почвы, а также разработки рекомендаций по дезактивации почв.

СПИСОК ЛИТЕРАТУРЫ

1. Пронкин Н.С., Шарафутдинов Р.Б., Ковалевич О.М. и др. Классификация водоемов-хранилищ жидких радиоактивных отходов по опасности // Атомная энергия. - 2003. - Т. 94. -№ 6. - С. 449-457.

2. Вишневский Ю.Г., Ирюшкин В.М., Кислов А.И. и др. О регулировании безопасности при обращении с жидкими радиоактивными отходами, накопленными в водоемах-хранилищах ПО «Маяк», Сибирского химического комбината и Горно-химического комбината // Вестник Госатомнадзора России. - 2002. - № 3. - С. 3-12.

3. Вабищевич П.Н., Данияров А.О. Математическое моделирование промачивания зоны аэрации в условиях близкого залегания грунтовых вод // Математическое моделирование. - 1994. - Т. 6. - № 11. - С. 11-24.

4. Шейдеггер А.Э. Физика течения жидкостей через пористые среды. - М.: Гостоптехиздат, 1960. - 250 с.

5. Арутюнян Р.В., Большов Л.А., Зенич Т.С., Решетин В.П. Математическое моделирование вертикальной миграции в почве 137, 134С§ // Атомная энергия. - 1993. - Т. 74. - № 3. - С. 223-230.

6. Голубев В.С. Динамика геохимических процессов. - М.: Недра, 1981. - 208 с.

7. Кобранова В.Н. Петрофизика. - М.: Недра, 1986. - 392 с.

8. Хргиан А.Х. Физика атмосферы. - М.: Государственное издательство физико-математической литературы, 1958. - 476 с.

9. Таблицы физических величин. Справочник / Под ред. акад. И.К. Кикоина. - М.: Атомиздат, 1976. - 1008 с.

УДК 519.6

ГЕОТЕХНОЛОГИЧЕСКИЙ ИНФОРМАЦИОННО-МОДЕЛИРУЮЩИЙ КОМПЛЕКС ДЛЯ ОПТИМИЗАЦИИ ПРОЦЕССА ПОДЗЕМНОГО ВЫЩЕЛАЧИВАНИЯ УРАНА

А.А. Жиганов, М.Д. Носков, А.Д. Истомин, А.Г. Кеслер, Н.С. Невзорова

Северский государственный технологический институт E-mail: nmd@ssti.ru

Представлен геотехнологический информационно-моделирующий комплекс, предназначенный для оптимизации добычи урана методом подземного выщелачивания. Комплекс включает в себя геотехнологическую информационную и моделирующую системы. Геотехнологическая информационная система позволяет вводить, редактировать и отображать параметры, характеризующие состояние продуктивного горизонта, а также характеристики технологических скважин. Моделирующая система описывает физико-химические процессы, происходящие в продуктивном горизонте при подземном выщелачивании урана. Комплекс может быть использован для оптимизации разработки месторождений урана методом подземного выщелачивания с целью повышения рентабельности разработки месторождений, увеличения доли извлекаемого урана и снижения экологической нагрузки на окружающую среду.

Методом подземного выщелачивания (ПВ) разрабатывают экзогенные месторождения урана, которые находятся в хорошо проницаемых подземных водоносных горизонтах [1-3]. Извлечение урана из рудного тела происходит через систему технологических скважин. Через нагнетательные скважины в продуктивный горизонт нагнетается раствор веществ, способных растворять содержащие уран минералы. Образующийся в подземном водоносном горизонте продуктивный раствор извлекается через откачные скважины. Образующие-

ся после переработки продуктивных растворов маточные растворы доукрепляются выщелачивающими реагентами и снова подаются в нагнетательные скважины в качестве рабочих растворов.

Основными задачами управления геотехнологическим предприятием является повышение рентабельности разработки месторождения, увеличение доли урана, извлекаемого из продуктивного горизонта, и снижение загрязнения подземных вод. Для решения этих задач нужно уметь оценивать оставшиеся запасы урана, располагать информацией о