CC BY
32
Вычислительные технологии
Том 22, Специальный выпуск 1, 2017
Численное исследование распространения примеси в пресном озере на основе распределения мутности воды
Б. О. Цыденов
Национальный исследовательский Томский государственный университет, Россия Контактный e-mail: [email protected]
Описана математическая 2.5Б-модель для воспроизведения гидродинамических процессов распространения примеси загрязняющих веществ в водоеме на примере озера Камлупс (Британская Колумбия, Канада). В качестве индикатора загрязнения рассматривается показатель мутности воды. В негидростатической модели учитывается почасовое изменение коротко- и длинноволновой радиации, потоков скрытого и чувствительного тепла. Получено качественное соответствие результатов моделирования данным натурных измерений. Показано влияние атмосферных параметров на пространственно-временное распределение мутности воды в озере.
Ключевые слова: лимнологическое моделирование, численный эксперимент, загрязнение озер, мутность воды, озеро Камлупс.
Введение
В связи с резким увеличением в последние десятилетия антропогенной нагрузки на окружающую среду, ведущей к загрязнению крупных озер, начинает возрастать роль чистой пресной воды. По самым пессимистичным прогнозам [1], глобальный водный кризис ожидается уже в 2025 г. Большая часть поверхностных пресных вод сосредоточена в озерах. Понимание природных механизмов озерной гидродинамики необходимо для правильного выбора стратегии предотвращения загрязнения воды. Следовательно, технологии прогнозирования и мониторинга качества воды станут весьма востребованными в ближайшем будущем.
Современные математические модели гидродинамики глубокого водоема [2 - 7] включают систему нелинейных уравнений в частных производных, представляющих собой математическую запись законов сохранения. Однако в существующих моделях, как правило, не учитывается важнейший фактор формирования течений — реальная суточная изменчивость потока тепла на поверхности озера, а задается постоянное среднемесячное значение солнечной радиации. Следует также заметить, что такие сезонные факторы, как колебание расхода воды притока, изменение термического и гидрохимического режима играют ключевую роль в формировании структуры течений и определяют качественные характеристики переноса загрязнения в районе впадения реки в озеро. Поэтому для более реалистичного воспроизведения конвективно-диффузионных процессов методами математического моделирования необходимо задавать дополнительные параметры на границе вычислительной области.
© ИВТ СО РАН, 2017
Со стороны атмосферы теплообмен с водоемом происходит в основном в результате взаимодействия таких факторов, как:
• радиационные процессы, включающие исходящую от Солнца коротковолновую радиацию (длина волн от 0.1 до 4 мкм) и длинноволновую радиацию , излучаемую атмосферой и водной поверхностью (длина волн от 4 до 120 мкм);
• турбулентный перенос тепла, состоящий из потока скрытого тепла Н^, связанного с испарением или конденсацией влаги, а также потока чувствительного тепла Дз, образующегося при перепаде температур между водой и приповерхностным воздухом [8].
Следует заметить, что солнечная радиация, достигающая поверхности Земли, на 99 % является коротковолновой, поэтому эти два названия часто отождествляются.
Для оценки качества воды удобно использовать показатель мутности, характеризующий зависимость прозрачности воды от наличия в ней мелкодисперсных взвешенных примесей различного происхождения.
Целью работы является разработка численной модели для исследования процессов загрязнения в пресном водоеме на основе распределения мутности воды с учетом суточной изменчивости параметров атмосферы.
1. Математическая модель и численный метод 1.1. Основные уравнения модели
Негидростатическая 2.5В-модель для воспроизведения динамики распространения загрязнения на основе пространственно-временного распределения мутности в водоеме, которая учитывает эффект вращения Земли и записана в приближении Буссинеска, включает в себя следующие уравнения: • уравнение мутности воды
дФ диФ дтФ дЪ дх дх
дх \ ж дх ) д х * "
£);
• уравнения количества движения
ди ди2 дит дЪ дх дх
дь диь дЪ дх дх
1 дт> д ( ди\ д „
:--/ + ^ + — ( ) + ь -
р0 дх дх \ ' 1 '
дwv
дх) дх
(«• £)
д (ду\ д {тдь\ ^ ^
7Г + ТГ + - 2П,и,
дх \ дх) дх \ дх )
дт дит дчи2 дЪ дх дх
1 др д
р0 дх дх
дх у
дт\ д
дх у дх
- — + — ( Кх— ) + — ^^ - — + 2ПУи - 2Пхь;
91 Ро
• уравнение неразрывности
• уравнение энергии
ди дт
дх дх '
дТ дпГ дчуТ дЬ дх дх
д / дТ\ д / дТ\ 1 дН.
ш{в'ш) + ж) + ^ -г.
во1
уравнения баланса солености в озере
№ + диЯ + dw<S = д_ iD dS_N + д_ (D dS_\ dt dx dz дх \ ж дх J dz\ z dz J
Здесь и, V — горизонтальные компоненты скорости; т — вертикальная компонента скорости; 0,х, и — компоненты вектора угловой скорости вращения Земли; д — ускорение свободного падения; ср — удельная теплоемкость; Т — температура; Б — соленость; р — давление; р0 — плотность воды при стандартном атмосферном давлении, температуре Т^ и солености Для учета влияния силы Кориолиса вдоль однородного направления (Оу) в модель включена вторая горизонтальная составляющая скорости V, которая может в некоторых областях доминировать над другими компонентами [2].
Поглощение коротковолновой радиации Н301 рассчитывается по закону Бугера — Ламберта — Бэра:
Н301 = Н3зо1,0 (1 - Г3) ехр (-баЬзЛ) ,
где т ^ ~ 0.2 — коэффициент отражения воды; ^ 0.3 м 1 — коэффициент поглощения солнечной радиации в воде; й = \Ьг — х\ — глубина, м. Приток солнечной радиации на поверхность озера Н$ЗО1,0 определяется следующим соотношением [9]:
\ S0 (ад — aw) cos ( [а (С) + b (С) ln (cos ()] , cos(() > 0; HSsolfl = <
10, cos(C) < 0.
Здесь ¿"о = 1367 Вт/м2 — солнечная постоянная; а (С) и Ь (С) — эмпирические коэффициенты, зависящие от степени покрытия небесного свода облачностью С [10]; £ — зенитный угол Солнца; эмпирические функции
ад = 0.485 + 0.515 (l.014 — 0.16/^cos
aw = 0.039 (rw/ cos ()0'3
представляют соответственно молекулярное рассеяние и поглощение перманентных газов; rw — содержание водяного пара в атмосфере, кг/м2.
Для определения коэффициентов диффузии Dx, Dz, Кх, Kz используется k-ш двух-параметрическая модель турбулентности Уилкокса [11], состоящая из уравнений для кинетической энергии и частоты турбулентных пульсаций [12].
В качестве уравнения состояния выбрано уравнение Чена — Миллеро [13], связывающее плотность воды р с температурой Т, соленостью S, давлением р и справедливое в диапазоне 0 < Т < 30 °C, 0 < S < 0.6 г/кг, 0 < р < 180 бар.
1.2. Начальные и граничные условия
Начальные условия для уравнений математической модели (2)-(7) задаются в виде и = 0, v = 0, w = 0, Ф = Фь, Т = TL, S = SL
при Ь = 0, где Ф^, Тъ, вь — мутность, температура, соленость в озере соответственно; £ — время. Начальное поле давления определяется из решения уравнений состояния
и гидростатики с граничным условием на поверхности р = ра (ра — атмосферное давление) методом Рунге — Кутты четвертого порядка точности. Уравнение гидростатики выводится из уравнения (4) при и = V = = 0:
др
дх
Граничные условия для уравнений (2)-(7) имеют вид а) на свободной поверхности
= к* = ^ит!, ш = о, ?Ф = 0. Л? = ^ = о.
' (Л 7 ^ ^ 7 ^ ^ 7 Г\ ^ 7 — £ С\ 'О
ох р0 ох р0 ОХ ОХ р0 Ср ох
Сдвиговое напряжение ветра на поверхности озера описывается законом
и 2 2 V 2 2
Г8иг! = СЮРауЩо + У2о • Що. Г8игГ = СюРа^ «2о + Що • УЮ.
где ра — плотность воздуха у поверхности воды; и1о, ь1о — составляющие скорости ветра на высоте 10 м над поверхностью озера; с10 = 1.3 • 10-3. Тепловой поток Нп^ складывается из длинноволновой радиации Н^ и скрытого Н^ и чувствительного Нз тепла:
Нп et = Нг,ш + Нь + Н8.
Параметризация компонентов тепловых потоков осуществляется по модели 3 из [14]:
• Щт = (1 + 0.17С2) Т\ - ьшаТ4. где Та — температура воздуха, К; Т — температура воды, К; а = 5.669 • 10-8 Вт/м2/К4 — коэффициент Стефана— Больцмана; ~ 0.96 — коэффициент излучения воды. Коэффициент излучения атмосферы рассчитывается по формуле
£а = СеТА.
здесь С€ = 9.37 • 10-6°С-2 — излучательная способность воздуха;
( 17.67 Т?
• Нь = ¡и (еА - е.ш). и = 6.9 + 0.345 и%0. = 6.112 ехр ( у? + ^ 5
где еш — давление насыщенного водяного пара, ГПа; еа = 0.01КН еш — давление водяного пара в атмосфере, ГПа; КН — относительная влажность, %; ¡и — коэффициент массоотдачи, Вт/м2/ГПа; и1о = л/и\о + — скорость ветра, м/с; Т? — температура воздуха, °С;
• Н3 = $ и (ТА -Т). где [ = 0.62 ГПа/К;
б) на твердых границах (на дне)
<9Ф дТ Нкео дБ
и = 0. г; = 0. ю = 0. — = 0. — =--^. — = 0.
Оп Оп Ро Ср Оп
где Нёео — геотермальное тепло; п — направление внешней нормали к области; о = 0.05 м;
в) на границе втекания реки (на левой границе)
и = ии. V = 0. т = 0. Ф = Фд. Т = Тд. Б = Бя.
где ии — скорость речного притока; Фд, Тд, ви — мутность, температура, соленость реки соответственно;
г) на открытой границе задаются условия радиационного типа [15]
, п М Ъ гр с
~дI + °фдх = 0, Ф = и,у, Ф,1,^, и простые градиентные условия
д'ш дх
1.3. Численный метод решения уравнений модели
Решение задачи основано на методе конечного объема, согласно которому скалярные величины (мутность, температура, минерализация воды) определяются в центре сеточной ячейки, а компоненты вектора скорости — в средних точках на границах ячеек. В целях приближения расчетной области к прибрежному профилю озера применяется метод блокировки фиктивных областей [16]: приравниваются нулю компоненты скорости в выключенной зоне за счет использования в ней больших значений коэффициентов вязкости.
Численный алгоритм нахождения поля течения и температуры опирается на разностную схему Кранка — Николсон. Конвективные слагаемые в уравнениях аппроксимируются по противопотоковой схеме QUICK [17] второго порядка. Для согласования рассчитываемых полей скорости и давления разработана процедура SIMPLED (Semi-Implicit Method for Pressure Linked Equations with Density correction) для течений с плавучестью [18], представляющая собой модификацию метода SIMPLE Патанкара [16]. Алгоритм SIMPLED корректирует поля скорости и давления с учетом вариации плотности в гравитационном члене уравнения количества движения для вертикальной составляющей. Системы разностных уравнений на каждом шаге по времени решаются методом релаксации. Результаты тестирования численного алгоритма на примере квадратной каверны представлены в [19].
2. Результаты и обсуждение
В качестве исследуемой области выбрано вертикальное сечение озера Камлупс [20], соответствующее направлению впадения реки Томпсон. Озеро Камлупс находится на юго-западе Канады (провинция Британская Колумбия) в 340 км северо-восточнее Ванкувера и расположено между 50°26' -50°45' с.ш. и 120°03'-120°32' з.д. по течению реки Томпсон, имеет вытянутую форму (рис. 1, а). Расчетная область имеет протяженность 10 км и глубину 138 м (рис. 1, б), начало системы координат совпадает с устьем реки (рис. 1, а). Глубина участка на границе раздела река-озеро составляет 15 м.
В вычислительном эксперименте начальное распределение температуры в озере Кам-лупс имеет постоянное значение, равное 2.4 °С, в то время как температура воды в реке соответствует 3.6 °С и нагревается на 0.2 °С в сутки [2]. Река Томпсон впадает в озеро со скоростью 0.01 м/с, минерализация воды в озере и реке равна 0.1 г/кг. Начальное значение мутности озерной воды составляет 2.0 ЛТИ, изменение мутности в устье реки представлено на рис. 2. [20]. Расчетная область (рис. 1, б) покрывается равномерной ортогональной сеткой с шагами Ъх = 25 ми кх = 3 м. Шаг по времени 60 с.
В качестве атмосферных данных выступает информация о температуре воздуха, относительной влажности, атмосферном давлении, облачности из архива погодных усло-
О '1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Длина, м
Рис. 1. Морфометрия оз. Камлупс: а — батиметрия, б — расчетная область (продольное сечение)
Рис. 2. Мутность воды в устье р. Томпсон с 01.04.1975 по 10.05.1975 г. [20]
а
вий метеостанции г. Камлупс в период с 01.04.2015 по 10.05.2015 г.1. На основе этих сведений рассчитаны значения коротко- и длинноволновой радиации, потоков скрытого и чувствительного тепла (рис. 3). В вычислительном эксперименте значения т,^
V
и Тдиг£ приняты равными нулю.
Пространственно-временное распределение мутности воды в озере Камлупс (рис. 4) показывает, что в начальной стадии моделирования примесь, поступающая из реки Томпсон, распространяется в глубоководную часть озера за счет возникающего под действием силы тяжести вдольсклонового течения (рис. 4, а). Затем по мере усиления темпов теплонакопления и прогрева речного притока менее плотная водная масса распространяется по поверхности озера, что способствует концентрации максимальных
1http://climate.weather.gc.ca
Рис. 3. Компоненты тепловых потоков, рассчитанные на основе атмосферных данных метеорологической станции г. Камлупс в период с 01.04.2015 по 10.05.2015 г. (местное стандартное время)
4000
6000
Расстояние, м
10000
Рис. 4. Распределение мутности воды (ЛХИ) в оз. Камлупс на 5-е (а) и 35-е (б) сутки моделирования (устье р. Томпсон слева)
а
значений мутности в верхних слоях водоема (рис. 4, б). Полученные результаты численного эксперимента (рис. 4) качественно соответствуют измеренному в апреле 1975 г. распределению мутности воды в районе впадения реки Томпсон в озеро Камлупс [20] (рис. 5).
Для определения влияния приходящего на водное зеркало переменного теплового излучения, наряду с базовым (далее эксперимент 1), проведен дополнительный расчет с постоянным потоком 127 Вт/м2 (далее эксперимент 2), соответствующим осредненно-
му значению в период с 01.04.2015 по 10.05.2015 г. Результаты динамической картины вертикального распределения мутности (рис. 6) демонстрируют, что переменный тепловой поток по сравнению с постоянным дает эффект замедления распространения за-
Рис. 5. Распределение мутности воды (ЛТИ) в оз. Камлупс 1 апреля 1975 г. [20] (устье р. Томпсон справа)
Рис. 6. Хроноизоплеты распределения мутности воды (ЛТИ) в оз. Камлупс при х = 2500 м с постоянным (а) (127 Вт/м2) и переменным (б) потоками тепла
Рис. 7. Хроноизоплеты распределения мутности воды (ЛХИ) при х = 2500 м с переменным потоком тепла на расчетной сетке с шагами Нх = 25 м, Нг = 3 м (а) и Нх = 12.5 м, Нг = 1.5 м (б)
грязняющих веществ на поверхности озера: если в эксперименте 2 речная водная масса с мутностью 2.5 ЛТИ достигает расстояния 2500 м от устья реки Томпсон на шестые сутки (рис. 6, а), то в эксперименте 1 — на седьмые сутки (рис. 6, б). Причиной этого служат низкие значения суммарной тепловой энергии в ночное время, которые приводят к уменьшению скорости поверхностного течения [14]. Однако при учете почасового изменения атмосферных параметров примеси загрязняющих веществ быстрее достигают придонной области: в эксперименте 1 водная масса с мутностью 4 ЛТИ приходит на 10-11-е сутки (рис. 6, б), в то время как в эксперименте 2 — только на 17-18-е сутки (рис. 6, а). Расчеты на более подробной сетке с шагами кх = 12.5 ми Нг = 1.5 м (с переменным потоком тепла) не выявили существенных различий результатов (рис. 7).
Задержка в распространении примеси (см. рис. 6) связана, вероятнее всего, не столько с амплитудой колебаний теплового потока (летом амплитуда может быть существенно выше благодаря солнечной радиации), сколько с сезонными особенностями значений теплового потока: в весенний период в ночное время происходит выхолаживание водной поверхности, приводящее к увеличению плотности воды, что замедляет скорость течения водных масс и может инициировать локальные опускания вод при температуре максимальной плотности (4 °С).
Заключение
Разработана математическая 2.5В-модель для воспроизведения динамической картины распространения примеси загрязняющих веществ в пресном озере на основе распределения мутности воды. Проведена апробация модели на примере канадского озера Камлупс. Результаты расчетов продемонстрировали, что на начальном этапе моделирования загрязнение, поступающее вместе со стоком реки Томпсон, распространяется в глубинную область озера, а затем, когда происходит интенсивный прогрев вод притока и поверхности водоема, концентрация максимальных значений мутности наблюдается в верхних слоях озера.
Анализ влияния почасовой вариации метеорологических параметров на динамику распределения загрязняющих веществ показал, что переменный тепловой поток, в отличие от постоянного, дает эффект более медленного распространения примесей в поверхностных слоях и более быстрого продвижения в глубоководную часть озера.
Благодарности. Работа выполнена при финансовой поддержке Президента РФ (грант № МК-4790.2016.5).
Список литературы / References
[1] Данилов-Данильян В.И. Водные ресурсы — стратегический фактор долгосрочного развития экономики России // Вестн. РАН. 2009. Т. 79, № 9. С. 789-798. Danilov-Danilyan, V.I. Water — a strategic factor for the long-term development of the Russian economy // Vestn. RAN. 2009. Vol. 79, No. 9. P. 789-798. (In Russ.)
[2] Holland, P.R., Kay, A., Botte, V. Numerical modelling of the thermal bar and its ecological consequences in a river-dominated lake //J. Mar. Syst. 2003. Vol. 43, No. 1-2. P. 61-81.
[3] Lawrence, S.P., Hogeboom, K., Le Core, H.L. A three-dimensional general circulation model of the surface layers of Lake Baikal // Hydrobiologia. 2002. Vol. 487. P. 95-110.
[4] Vasiliev, O.F., Bocharov, O.B., Kvon, V.I., Ovchinnikova, T.E., Kvon, D.V. Numerical modeling of thermal bars in deep lakes // Proc. of the 3-rd Intern. Conf. on HydroScience and Engineering. Cottbus; Berlin, 1998. 20 p. (CD-ROM).
[5] Цветова Е.А. Особенности конвекции подо льдом в глубоком озере // Интерэкспо ГеоСибирь. 2016. Т. 4, № 1. С. 110-114.
Tsvetova, E.A. Features of convection under ice in a deep lake // Interexpo GEO-Siberia. 2016. Vol. 4, No. 1. P. 110-114. (In Russ.)
[6] Блохина Н.С. Влияние ветра на развитие термобара и течений в водоемах различной глубины в период таяния ледового покрова // Вестн. Мос. ун-та. Сер. 3: Физика. Астрономия. 2015. № 4. С. 102-108.
Blokhina, N.S. The influence of wind on the development of a thermal bar and currents in reservoirs of different depths during ice cover melting // Moscow Univ. Physics Bulletin. 2015. Vol. 70, No. 4. P. 319-325.
[7] Farrow, D.E. A model for the evolution of the thermal bar system // EJAM. 2013. Vol. 24, No. 2. P. 161-177.
[8] Цыденов Б.О. Численное моделирование эффекта весеннего термобара в глубоком озере: Дис. ... канд. физ.-мат. наук. Томск, ТГУ, 2013. 145 с.
Tsydenov, B.O. Numerical modeling of the effect of the spring thermal bar in a deep lake: Dis. ... kand. fiz.-mat. nauk. Tomsk, TGU, 2013. 145 p. (In Russ.)
[9] Hurley, P. The air pollution model (TAPM) Version 2. Pt 1: Technical description // CSIRO Atmospheric Research Technical Paper No. 55. Australia: CSIRO, 2002. 51 p.
[10] Александрова М.П., Гулев С.К., Синицын А.В. Уточнение параметризации коротковолновой радиации на поверхности океана на основе прямых измерений в Атлантическом океане // Метеорология и гидрология. 2007. № 4. С. 45-54.
Aleksandrova, M.P., Gulev, S.K., Sinitsyn, A.V. An improvement of parametrization of short-wave radiation at the sea surface on the basis of direct measurements in the Atlantic // Russian Meteorology and Hydrology. 2007. Vol. 32, No. 4. P. 245-251.
[11] Wilcox, D.C. Reassessment of the scale-determining equation for advanced turbulence models // AIAA J. 1988. Vol. 26, No. 11. P. 1299-1310.
[12] Цыденов Б.О., Старченко А.В. Применение двухпараметрической к - ш модели турбулентности для исследования явления термобара // Вестн. Том. гос. ун-та. Математика и механика. 2014. № 5(31). С. 104-113.
Tsydenov, B.O., Starchenko, A.V. Application of the two-parametric к -ш turbulence model for studying the thermal bar phenomenon // Tomsk State Univ. J. Math. Mech. 2014. No. 5(31). P. 104-113. (In Russ.)
[13] Chen, C.T., Millero, F.G. Precise thermodynamic properties for natural waters covering only limnologies range // Limnol. Oceanogr. 1986. Vol. 31, No. 3. P. 657-662.
[14] Tsydenov, B.O., Starchenko, A.V. To the selection of heat flux parameterization models at the water-air interface for the study of the spring thermal bar in a deep lake // Proc. SPIE. 2015. Vol. 9680. P. 1-8.
[15] Orlanski, I. A simple boundary condition for unbounded hyperbolic flows //J. Comput. Phys. 1976. Vol. 21, No. 3. P. 251-269.
[16] Patankar, S. Numerical heat transfer and fluid flow. N.Y.: CRC Press, 1980. 214 p.
[17] Leonard, B. A stable and accurate convective modeling procedure based on quadratic upstream interpolation // Computer Methods in Applied Mechanics and Engineering. 1979. Vol. 19, No. 1. P. 59-98.
[18] Tsydenov, B.O., Kay, A., Starchenko, A.V. Numerical modeling of the spring thermal bar and pollutant transport in a large lake // Ocean Modelling. 2016. Vol. 104. P. 73-83.
[19] Цыденов Б.О. Численное моделирование конвективных течений в каверне // Тр. Меж-дунар. конф. "Перспективы развития фундаментальных наук". Томск: ТПУ, 2009. Т. 2. С. 673-676.
Tsydenov, B.O. Numerical modelling convective flows in a cavity // Proc. of the Intern. Conf. "Prospects of Fundamental Sciences Development". Tomsk: TPU, 2009. Vol. 2. P. 673676. (In Russ.)
[20] The limnology of Kamloops Lake, B.C. / B.E. St. John, E.C. Carmack, R.J. Daley, C.B.J. Gray, C.H. Pharo. Vancouver, 1976. 167 p.
Поступила в 'редакцию 24 января 2017 г.
A numerical study of impurity propagation in a freshwater lake on the basis of water turbidity distribution
tsydenov, bair o.
Tomsk State University, Tomsk, 634050, Russia
Corresponding author: Tsydenov, Bair O., e-mail: [email protected]
A mathematical 2.5D model for simulating the hydrodynamic processes of pollutant propagation on an example of Kamloops Lake (British Columbia, Canada) is presented. Water turbidity is considered as an indicator of pollution. The non-hydrostatic model includes the continuity, momentum, energy, salinity, and turbidity equations and takes into account the diurnal variability of the fluxes of shortwave and longwave radiation and latent and sensible heat.
Turbulence closure of the system is performed with a two-parameter k - w Wilcox-type model and algebraic relations for the coefficients of turbulent diffusion. The convection-diffusion equations are solved by a finite volume method to satisfy the integral conservation laws. The numerical algorithm for finding the flow and temperature fields is based on a Crank-Nicholson difference scheme. The convective terms in the equations are approximated with a second-order upstream scheme, QUICK. The systems of grid equations at each time step are solved by an under-relaxation method.
Numerical modelling has demonstrated that the river water with high values of turbidity first moves along the slope to the deeper parts, then spreads over the lake
© ICT SB RAS, 2017
124
B. O. ЦbIценов
surface. The results of simulation have shown qualitative agreement with in-situ measurement described by John et al. (1976). It was also found that the variable heat flux affects the space-time distribution of turbidity.
Keywords: lake modelling, numerical simulation, water bodies contamination, water turbidity, Kamloops Lake.
Acknowledgements. This research was supported by the Grant of the President of Russian Federation (No. MK-4790.2016.5).
Received 24 January 2017
