Геохимия, 2020, T. 65, № 10, стр. 1009-1024

ВВЕДЕНИЕ

Рост глобальной температуры, наблюдаемый в последние десятилетия, связывается многими исследователями со стремительным ростом содержания в атмосфере парниковых газов, и прежде всего СО₂ [Grace 2000; IPCC 2013], темп роста концентрации которого в конце 2000-х годов достиг 2.20 ± 0.01 ppm/год (или 1.7% за год) и продолжает увеличиваться (https://www.esrl.noaa.gov/ gmd/ccgg/trends/global.html#global). Максимальная концентрация CO₂, измеренная в обсерватории Мауна-Лоа (на Гавайях 19°29′ с.ш. 155°36′ з.д. на высоте 3400 м), в январе 2020 г. достигла 412.30 ppm. Не существует однозначного мнения касательно наблюдаемого роста концентрации СО₂ в глобальном масштабе. Возможно, этот процесс связан с естественными колебаниями, и происходит не только за счет возрастающего объема поступлений из антропогенных источников (Peters et al., 2012).

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

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

Данная статья является продолжением серии работ по моделированию физиологических процессов сосны обыкновенной (Pinus sylvestris). В работе (Соколов и др., 2019а) мы подробно рассмотрели набор моделей динамики транспирации и выбрали их модель, соответствующую набору имеющихся экспериментальных данных. Выбор осуществлялся на основе метода сбалансированной идентификации (Соколов, Волошинов, 2018; Sokolov, Voloshinov, 2019; Соколов и др., 2019б) путем минимизации среднеквадратичной ошибка кроссвалидации. С использованием полученных результатов (временные ряды значений устьичной проводимости и дефицита воды в растении) и дополнительных наборов данных, полученных в результате пассивных (наблюдения в естественных условиях) и активных экспериментов, мы попытаемся описать путь СО₂ в листе, оценить работу составных частей биохимического конвейера фиксации CO₂.

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

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

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

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

РОЛЬ ВНЕШНИХ ФАКТОРОВ СРЕДЫ И ВНУТРЕННИХ ПАРАМЕТРОВ В МОДЕЛЯХ ФОТОСИНТЕЗА НА УРОВНЕ ЛИСТА

Процесс поглощения CO₂ растением (и соответствующую модель) можно представить в виде потоковой диаграммы (рис. 1).

Рис. 1.

Потоки CO₂ в листе. Схематическое представление модели конвейера поглощения CO₂ (обозначения см. в табл. 1).

Углекислый газ из атмосферы (форма 1) через устьица (форма 2) поступает (диффундирует согласно закону Фика) в межклеточное пространство (форма 3), откуда молекулы СО₂ диффундируют (согласно закону диффузии растворенных газов в жидкости) сквозь стенки мезофильных клеток (форма 4) к реакционным центрам в хлоропластах (форма 5), где и вступают в биохимический цикл фотосинтеза (форма 6), результатом которого является синтез ассимилятов (форма 8). На рис. 1 отображен и противоположный процесс – дыхание (форма 7), т.е. потребление богатых химической энергией веществ на биохимические процессы в различных частях листа. Подробное изучение дыхания не является целью данной работы, поэтому для простоты предположим, что появляющийся при этом углекислый газа попадает непосредственно в хлоропласты.

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

Рассмотрим более подробно элементы диаграммы на рис. 1

В основном, в процессе фотосинтеза используется CO₂ атмосферы (форма 1). Небольшая часть углекислоты может доставляться к листьям вместе с ксилемным соком, но она обычно не превышает нескольких процентов.

Концентрация СО₂ имеет выраженную суточную динамику, определяемую процессами фотосинтеза и дыхания (Lloyd et al., 2002). Суточный максимум атмосферной концентрации CO₂ регистрируется в 6:00–7:00 ч местного времени, минимум – в вечерние часы (16:00–19:00 ч). Максимальная величина суточной амплитуды для CO₂ наблюдается в июле около поверхности земли и составляет в среднем 21 ppm (Тимохина и др., 2014). Результаты, полученные в Швеции (Perttu et al., 1980) в конце августа в 120-летнем сосновом древостое (высота – 15.8 м) показали, что с 6 до 21 ч колебания концентрации СО₂ во времени составляли 10–12 ppm от среднего уровня. Небольшое увеличение концентрации СО₂ происходило с 20 до 24 ч. В дневное время у поверхности почвы концентрация СО₂ на 10–15 ppm выше, чем в кронах.

Известно, что при концентрациях, близким к естественным, зависимость фотосинтеза сосны от величин концентрации СО₂, близких к естественным, носит линейный характер (Голомазова, 1987, Щербатюк и др., 1992) и некоторое ее уменьшение днем снижает интенсивность процесса. В работе А.Г. Молчанова (1983) выявлено воздействие изменений естественной концентрации СО₂ на световые кривые фотосинтеза. Снижение концентрации СО₂ на 7% вызывало падение интенсивности фотосинтеза на плато до 30%.

Молекулы атмосферного СО₂ через пограничный слой воздуха над листом, попадают через устьичные отверстия (форма 2) в межклеточное пространство (форма 3). Проводимость пограничного слоя для хвои сосны высока – 8.1 моль H₂O м^–2 с^–1 (Stangl et al., 2019), т.е. сопротивление пограничного слоя достаточно мало, и им в расчетах можно пренебречь или включить его в устьичное.

Закон Фика

позволяет рассчитать устьичную проводимость для воды G_sw при известных значениях транспирации E, концентрации водяных паров в атмосфере W_a и насыщающей концентрации водяного пара в межклетниках при температуре листа W_i. В предположении равенства температуры листа и воздуха разность концентраций может быть заменена в первом приближении на дефицит давления водяного пара в воздухе VPD (Барри, Даунтон, 1987).

Из множества факторов, влияющих на устьичную проводимость, выделим солнечную радиацию Q и дефицит воды в растении WD, которые оказались ведущими факторами в модели водного режима (Соколов и др., 2019а). В летнее время (а именно, к этому сезону относится моделируемый период) температура в диапазоне 14–23°С оказывает меньшее влияние.

Аналогично водяному пару, поток CO₂ через устьица (фотосинтез) P определяется разностью концентраций углекислоты в воздухе (C_а) и межклетниках (C_i) и устьичной проводимостью для углекислого газа (G_s). Согласно закону Фика для CO₂:

Поток молекул CO₂ из межклеточного пространства (форма 3) через мезофилл (форма 4) в центры карбоксилирования в хлоропластах (форма 5) описывается законом диффузии растворенных газов в жидкости:

Мезофильная проводимость характеризует диффузию CO₂ из межклеточного пространства листа до реакционных центров в хлоропластах и включение в биохимический цикл реакций темновой стадии фотосинтеза (Цельникер и др., 1993, Laisk, Loreto, 1996). В ряде работ G_m разделяют на две составляющих: физическую и химическую компоненты (Laisk, Oja, 1998). Первая зависит в основном от длины пути СО₂ в жидкой фазе клетки и суммарной поверхности хлоропласта. Химическая составляющая определяется активностью карбоксилирующих ферментов, прежде всего РДФ-карбоксилазы. Анализ компонент мезофильного сопротивления показал, что доля физической компоненты значительно меньше, чем химической (Цельникер, 1978; Laisk, Oja, 1998). Мезофильную проводимость (обе составляющие) или, как ее иногда называют, внутреннюю проводимость связывают с первоначальным наклоном кривой зависимости ассимиляции A от C_i (Moldau, Kull, 1993; Laisk, Loreto, 1996; Eichelmann, Laisk, 1999; Eichelmann et al., 2004).

Однако, в последнее время все больше авторов термин “мезофильная проводимость” связывают только с диффузией CO₂ через мезофилл листа: межклеточные воздушные пространства, клеточную стенку и межклеточную жидкость. То есть, мезофильную проводимость рассматривают как синоним диффузионной проводимости листа, а химическая часть, которая часто лимитирует фотосинтез, уже не входит в понятие мезофильной проводимости. Появление новых приборов и методик привело к значительному росту числа работ, где изучается мезофильная проводимость – за последние 15 лет число ссылок, связанных с различными аспектами мезофильной проводимости для CO₂, увеличивалось экспоненциально, и превысило тысячу публикаций (Flexas et al., 2008, Stangl et al., 2019). Ошибки, вызываемые пренебрежением G_m при создании моделей фотосинтеза, значительно ухудшают их качество (Flexas et al., 2008).

Мезофильную проводимость в полевых условиях оценить значительно труднее, чем устьичную G_s, так как измерения G_m требуют одновременных измерений газообмена и изотопной дискриминации ¹³C (∆¹³C). Дискриминация может быть выведена из сигнатуры δ¹³C продуктов фотосинтеза, то есть растворимых сахаров листа, флоэмного содержимого (Hu et al., 2010; Ubierna, Marshall, 2011). Эти методы технически трудны, особенно в полевых условиях. Лишь в сравнительно недавно удалось с приемлемой точностью провести длительные одновременные измерения газообмена на уровне побега и ¹³C дискриминацию в столетнем древостое сосны обыкновенной на севере Швеции, получить почасовые оценки G_m параллельно с G_s, Р и E (Stangl et al., 2019).

Другим способом оценить концентрацию С_с, а затем подсчитать G_m, является проведение (активных) экспериментов по определению углекислотных компенсационных пунктов (УКП). Для этого побег помещают в замкнутую систему, где СО₂ постепенно “выедается” и через некоторое время СО₂-газообмен становится равным нулю. Концентрация углекислоты в камере, при которой это происходит, получила название УКП. В ряде моделей компенсационный пункт приравнивают к С_с(Лайск, 1977; Цельникер и др., 1993). У С4 растений УКП близок к нулю. У С3-растений, к которым относится сосна, в связи с фотодыханием концентрация CO₂ на реакционных центрах (С_с) может быть значительно больше нуля (Цельникер и др., 1978, Eichelmann еt al., 2004), особенно при низкой освещенности (Голомазова, 1987).

Проводимости G_m и G_s реагируют на изменение окружающих условий по-разному. Описаны быстрые реакции G_m на некоторые переменные окружающей среды. Эти переменные включают интенсивность или качество света (Loreto et al., 2009; Campany et al., 2016), межклеточную концентрацию CO₂ (C_i) (Flexas et al., 2007) и температуру листа (Warren, 2008).

При постоянной G_s рост G_m будет увеличивать эффективность использования воды (WUE) (Flexas et al., 2010), которая определяется как отношение нетто-ассимиляции углерода (Р) к потере воды в ходе транспирации (E). Это происходит из-за того, что увеличение G_m не имеет прямого влияния на транспирацию, но может увеличить фотосинтез, что приводит к увеличению отношения Р/E. В моделях WUE G_m часто вводится как постоянная, эмпирическая поправка в зависимости между C_i/C_a или экстраполируется на основании ее корреляции с G_s, хотя корреляция с G_s не всегда высокая (Klein et al., 2015).

Доставленный в хлоропласты СО₂ в центрах карбоксилирования вступает в биохимический цикл фотосинтеза (форма 6). Здесь, с помощью накопленных (в реакциях световой фазы фотосинтеза) богатых энергией веществ, CO₂ превращается (в ходе темновой фазы фотосинтеза) в глюкозу, а затем в крахмал (форма 8). Основные факторами, определяющими скорость фиксации CO₂ являются ФАР (фотосинтетически активная солнечная радиация), температура и концентрация C_c.

ОБЪЕКТ ИССЛЕДОВАНИЯ И ЭКСПЕРИМЕНТАЛЬНЫЕ ДАННЫЕ

Исследования проводились на базе полевого стационара Института леса КНЦ РАН, расположенного в 50 км к северу от г. Петрозаводска (62°13′ с.ш. и 34°10′ в.д.). Район исследований, пробные площади и конкретный объект изучения – 55-летнее сосна обыкновенная (Pinus sylvestris L.) высотой 14 м описаны в (Болондинский, Кайбияйнен, 2003, Соколов и др., 2019а). В ходе исследования, наряду с транспирацией (Веселков, Тихов, 1984) и микроклиматическими параметрами, исследовался СО₂-газообмен побегов в разных частях кроны (Болондинский, Кайбияйнен, 2003).

Для непрерывной автоматической регистрации СО₂-газообмена побегов сосны использовалась 6-канальная автоматическая установка на базе инфракрасного газоанализатора Infralyt-4 (Junkalor, DDR). В дальнейшем ряд показателей СО₂-газообмена уточнялся с помощью портативной газометрической системы LI-6200 (LI-COR, США). Измерения СО₂-газообмена проводились с частотой 10 показаний/ч по каждому каналу. По ним рассчитывались средние получасовые значения. Интенсивность фотосинтеза рассчитывалась как на сухой вес, так и на площадь хвои (мкмоль СО₂ м^–2 с^–1). Относительная погрешность измерений составляла 10–12%.

Ассимиляционные камеры изготавливались из прозрачного жесткого полиэтилена и устанавливались с южной стороны кроны на однолетних побегах 4–5-ой мутовок. В ходе измерения обеднение по CO₂ воздуха в камере не превышало 10%. В летнее время, чтобы избежать перегрева хвои, мы использовали дополнительную систему для увеличения потока воздуха через камеры. При прокачке воздуха через камеру со скоростью 2.8 л/мин воздух сменялся в камере каждые 10 с, и хвоя не перегревалась больше, чем на 1°С.

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

Температура воздуха рядом с камерами регистрировалась термодатчиками прибора ТЭТ-2 (Гидрометприбор, Латвия), защищенными от попадания прямых солнечных лучей цилиндрическим жестяным экраном. Разницу между температурой хвои и температурой окружающего воздуха измерялась медно-константановыми термопарами (диаметр проволоки 0.1 мм). В каждом побеге закреплялись 4–5 термопар.

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

Наиболее тщательно с минимальным количеством пропусков, не превышающим 1–2 суток в месяц, СО₂-газообмен измерялся с 18 апреля по 30 октября. Для моделирования был выбран фрагмент данных с 20 июня по 23 июля. Всего 1632 измерения. Данные по E имеют пропуски (1537 точек). Отсутствие данных было вызвано, в основном, техническими проблемами. Выбранный период попадал на окончание интенсивного апикального роста ствола (27.05–20.06). После ослабления аттрагирующих свойств верхушечной меристемы часть углеводов транспортировалась в зону камбиального роста и к корням. В конце июня начался интенсивный радиальный рост, закончившийся в 3 декаду июля. В июле наблюдался также интенсивный рост хвои и ее дозревание. В результате, в исследуемый период у дерева, также, как и во время апикального роста в высоту, наблюдалась повышенная потребность в ассимилятах.

В рассматриваемый период (июнь-июль) в 50-см слое почвы содержалось 35–45 мм влаги, что несколько меньше, чем в мае – 65 мм. Последняя величина близка к полевой влагоемкости. Минимум влажности почвы (4–5%) регистрировался в зоне наибольшего распространения корней (10–30 см) обычно в июле-августе. Хорошая увлажненность почвы объясняется существованием довольно мощного пласта грунтовых вод (уровень грунтовых вод 1–1.5 м) и достаточным количеством осадков в период исследований – 92 и 392 мм за вегетацию и за год соответственно. В период, выбранный для моделирования, сосна не испытывала недостатка почвенной влаги (Кайбияйнен, 1990).

В начале июля в течение трех суток проводились измерение содержания CO₂ в воздухе на разных высотах. Погрешность измерений не превышала 1%. Несмотря на сложную картину смены восходящих и нисходящих потоков СО₂, связанных с изменением мощности истоков и стоков СО₂ (фотосинтез, дыхание ветвей и ствола, выделение СО₂ из почвы), на фиксированных высотах в определенный момент времени концентрация СО₂ изменялась слабо. Полученные в результате значения приведены в табл. 2.

Углекислотный компенсационный пункт у сосны нами определялся как с помощью портативной системы Li-6200, так и на стационарной установке. Величину УКП мы приравнивали к концентрации CO₂ на реакционных центрах (С_с). В зависимости от внешних факторов среды и периода вегетации величина УКП колебалась днем (при больших значениях радиации) в пределах 55–67 мкмоль/моль (июле-августе – 65–68 мкмоль м^–2 с^–1). Уменьшение величины освещенности, начиная с 200 мкмоль м^–2 с^–1, приводило к увеличению УКП (концентрации CO₂) (рис. 2). Примерно такая зависимость отмечалась и другими исследователями (Голомазова, 1987).

Рис. 2.

Зависимость УКП (концентрации CO₂ в хлоропластах С_c) у побегов сосны от освещенности (Q).

ПРИМЕНЕНИЕ ТЕХНОЛОГИИ СБАЛАНСИРОВАННОЙ ИДЕНТИФИКАЦИИ ДЛЯ АНАЛИЗА МОДЕЛЕЙ ФОТОСИНТЕЗА

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

Формальные постановки задач идентификации не приводятся – предполагается, что используемая технология сбалансированной идентификации позволяет для заданной модели явления и используемого набора данных определить неизвестные функции и оценить погрешность построенной модели. Подробности можно найти в (Соколов, Волошинов, 2018; Sokolov, Voloshinov, 2019).

Для рассматриваемых ниже моделей основным множеством измерений

Для оценки погрешности моделирования используется процедура перекрестного оценивания с разбиением множества измерений D на 34 подмножества: каждое содержит измерения за одни сутки (до 48 точек). Полученную для такого разбиения методом сбалансированной идентификации среднеквадратичную ошибку кроссвалидации σ_cv, которую будем называть погрешностью моделирования данных, можно трактовать как меру предсказательной точности модели. Здесь и далее все погрешности даны в процентах от величины стандартного отклонения измерений (квадратный корень из дисперсии), что обеспечивает нормировку от 0% (при полном совпадении модельных значений со всеми измерениями) до 100% (для примитивной модели в форме постоянной функции, равной среднему значению измерений).

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

Замечание. Чтобы избежать путаницы, в обозначениях функций используются дополнительные строчные буквы, соответствующие аргументам. Например, Pq(Q) или Pqvpd(Q, VPD).

Модель 1: Pq(Q)

Результаты моделирования фотосинтеза как функции солнечной радиации приводятся на рис. 3а. Погрешность моделирования σ_cv = 51.84%.

Рис. 3.

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

Модель 4: Pqt(Q, T)

Результаты моделирования фотосинтеза как функции солнечной радиации и температуры на рис. 3г. Погрешность моделирования σ_cv = 51.00%.

Результаты моделирования показывают, что основным фактором, определяющим зависимость фотосинтеза от внешних параметров, является солнечная радиация. Как и следовало ожидать, добавление к основному фактор (Q) других факторов (VPD и T) не сильно уменьшает оценку погрешности. Это объясняется тесной взаимосвязью факторов внешней среды (некоторая синхронность или когерентность).

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

$Er{{r}_{k}} = {{P}_{k}}_{{}} - Pqt({{Q}_{k}},{{T}_{k}}),$

аппроксимируем функцией, зависящей от времени суток. Результаты, приведенные на рис. 4, показывают занижение модельного фотосинтеза в утренние часа и его завышение в вечерние. Это объясняется отсутствием в модели устьичного регулирования газообмена – утром, когда воды в растении достаточно, устьица открыты полностью, а после полудня устьица прикрываются из-за дефицита воды в растении.

Рис. 4.

Аппроксимация остатков моделирования фотосинтеза (как функции радиации и температуры) функцией времени суток: точки – остаточные ошибки, кривая – аппроксимация.

Модель 5: Устьичная проводимость

Воспользуемся устьичной проводимостью, рассчитанной ранее в модели транспирации (Соколов и др., 2019а). Это позволит на основе закона Фика рассчитать концентрацию СО₂ в межклетниках Ci:

$P = (Ca - Ci)\,Gs \times 1.0816;\,\,\,\,Ca = 360$

и представить фотосинтез в виде

$P = Pqt\left( {Q,T} \right)Pci\left( {Ci} \right);\,\,\,\,Pci\left( {360} \right) = 1,$

где функции мультипликативного представления должны удовлетворять условиям монотонности

$\begin{gathered} d{\text{/}}dQ(Pqt(Q,T)) \geqslant 0, \\ d/dCi(Pci(Ci)) \geqslant 0, \\ Pqt(0,T) \leqslant 0. \\ \end{gathered} $

Последние три условия представляются оправданными для биохимического процесса фотосинтеза – при увеличении потока энергии (солнечной радиации) и концентрации реагента (СО₂) скорость реакции не должна убывать, при отсутствии света фотосинтез не идет (но имеется ночное дыхание, включенное в P). Такие предположения подтверждаются экспериментально, например, световые кривые фотосинтеза демонстрируют насыщение по радиации при постоянстве других условий (Болондинский, Виликайнен, 2014).

Расчеты показывают заметное снижение погрешности σ_cv = 50.44%.

При анализе найденных функций мультипликативного представления фотосинтеза (рис. 5) следует отметить, что функция Pci(Ci) работает на фотосинтез при Ci в диапазоне 90–250 ppm – ночная концентрация CO₂ (около 360 ppm) быстро “выедается” при появлении света.

Рис. 5.

Фотосинтез как мультипликативная функция (модель 5). Два сомножителя: (а) зависимость от радиации и температуры (линии уровня); (б) зависимость от концентрации СО₂.

Модель 6: Ассимиляция (брутто фотосинтез) минус дыхание

Разделим (нетто) фотосинтез на две (знакоопределенные) составляющие: ассимиляцию CO₂ (брутто фотосинтез) и дыхание:

$\begin{gathered} \begin{array}{*{20}{l}} {P = Aqt(Q,T)Aci(Ci) - R(h),} \\ {Aqt(Q,T) \geqslant 0;\,\,\,\,Aqt(0,T) = 0,} \end{array} \\ 0 \geqslant R(h) \leqslant 0.2;\,\,\,R\left( 0 \right) = R\left( {24} \right). \\ \end{gathered} $

Обе функции мультипликативного представления ассимиляции предполагаются неубывающими и неотрицательными, при отсутствии света ассимиляция равна 0. Дыхание предполагается периодической, неотрицательной и ограниченной функцией времени суток. Погрешность снизилась σ_cv = 50.20%. Полученная функция дневной динамики дыхания приведена на рис. 6.

Рис. 6.

Дыхания как функция времени суток.

В анализируемый период ночное дыхание у однолетней хвои, в отличие от растущей хвои, имело низкий уровень. Отчасти этому способствовали низкие ночные температуры. При температурах 14–16°С дыхание достигало 0.2 мкмоль м^–2 с^–1. В утренние часы темновое дыхание также имело невысокие значения, постепенно возрастая к вечеру. Если температура воздуха была в дневное время умеренной (18–20°С), величины дыхания не превышали 0.2–0.3 мкмоль м^–2 с^–1 Дыхание в дневные часы значительно возрастало, если имела место полуденная депрессия фотосинтеза. Однако, в рассматриваемый период последнее явление практически не наблюдалось.

Модель 7: Изменение CO₂ в окружающем воздухе

Воспользуемся табл. 2, содержащей набор измерений атмосферного СО₂ в разное время суток, для аппроксимации суточного хода атмосферного СО₂ периодической функцией (рис. 7).

Рис. 7.

Суточный ход СО₂ в атмосфере: точки – измерения, кривая – аппроксимация их периодической функцией времени суток.

Заменим в модели постоянное значение СО₂ на найденную функцию времени суток. Результат – погрешность σ_cv = 50.07%.

Модель 8: Углекислотные компенсации точки

В точке компенсации фотосинтез равен нулю:

$Aqt\left( {Q,T} \right)Aci\left( {Ci} \right) - R\left( h \right) = 0$

т.е. ассимиляция (брутто фотосинтез) равна дыханию. Воспользуемся собранным материалом по УКП (рис. 2) и добавим в задачу идентификации (в целевую функцию) специальный штраф за отклонение ассимиляции от дыхания для всех записей УКП. Таким образом, измерения УКП доопределяют ассимиляцию при малых значениях концентрации CO₂. Добавление в задачу дополнительной качественной (наличие точек компенсации) и количественной (значения параметров в точках компенсации) информации приводит к уменьшению погрешности и существенному изменению функции Aci. Теперь, как и следовало ожидать, при минимальных значениях концентрации Сi скорость реакции ассимиляции СО₂ близка к нулю (см. рис. 8).

Рис. 8.

Ассимиляция как мультипликативная функция (модель 8). Два сомножителя: (а) зависимость от радиации и температуры (линии уровня); (б) зависимость от концентрации СО₂ в межклетниках.

Результат – погрешность σ_cv = 49.91%.

Модель 10: Мезофильная проводимость как функция радиации Gm(Q)

Погрешность увеличилась: σ_cv = 50.14. Это свидетельствует о том, что модель стала слишком сложной для имеющихся данных – слишком много промежуточных искомых функций, не подкрепленных экспериментом. Тем не менее, приводим найденную функцию (рис. 9 а). Полученный результат можно объяснить – при выходе фотосинтеза на насыщение при высоких значениях Q (более 2000) нет необходимости стимулировать поток СO₂ путем обеспечения высокой мезофильной проводимости (например, при помощи повышения концентрации и скорости регенерации рибулозобисфосфаткарбоксилаза/оксигеназы).

Рис. 9.

Мезофильная проводимость: как функция солнечной радиации (а), как функция дефицита воды в растении (б).

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ И ИХ ОБСУЖДЕНИЕ

Формально можно рассматривать проведенное исследование как попытку выбрать лучшую комбинацию гипотез и экспериментальных данных и объединить их в математической модели. Модификация моделей проводилась эволюционным способом от простого к сложному. В целом, критерием такого “эволюционного отбора” является погрешность кроссвалидации – если погрешность модифицированной модели уменьшалась, она “выживала” и использовалась как основа дальнейшей модернизации. Некоторая часть такого отбора была рассмотрена выше, другая (не меньшая по объему), содержащая тупиковые ветви, была опущена. Одна из таких ветвей была связана с мультипликативным представлением:

Сначала это направление успешно эволюционировало, однако, при добавлении в модель УКП проиграло основной ветке в борьбе за минимизацию погрешности.

Анализ (табл. 3) динамики погрешности моделирования данных σ_cv позволяет выделить ключевые, существенные гипотезы, привлечение которых приводит к заметному уменьшению погрешности: учет устьичной проводимости, учет дыхания и использование точек компенсации.

Полученные модели позволили рассчитать некоторые среднесуточные показатели, приведенные на рис. 10.

Рис. 10.

Среднесуточные значения: (а) фотосинтеза (P); (б) транспирации (E). (в) Эффективности использования воды (WUE). (г) Дефицит воды в растении. Прерывистые линии – соответствующие средние за период моделирования значения: aP = 2.52, aE = 0.196, aWUE = 12.8, aWD = 2.32.

Анализ среднесуточных значений фотосинтеза и транспирации (рис. 10а, 10б) показывает значительную синхронность их изменений в течение первых 25 дней анализируемого периода. В результате и величины WUE (рис. 10в) сохраняли определенное постоянство и имели всего два небольших выброса – на 16 и 22 день (5 и 11 июля). Как в прохладные пасмурные дни начала периода, так и далее на протяжении нескольких солнечных дней с 28 июня по 2 июля среднеднесуточные транспирационные расходы составляли 0.15 ± 0.06 ммоль м^–2 с^–1, что несколько ниже обычных значений для начала июля. В какой-то мере это было обусловлено невысокой температурой воздуха (18–20°С) и VPD, не превышающего 1 кПа.

В анализируемый период наблюдался радиальный рост дерева, а также рост молодой хвои. Наличие мощных аттрагирующих центров обуславливала большую потребность в ассимилятах, результатом чего был высокий среднесуточный уровень фотосинтеза (до 3.5 мкмоль м^–2 с^–1). Учитывая, что устьичная проводимость была ниже среднего уровня и составляла в первую половину дня 0.04–0.06 моль м^–2 с^–1, фотосинтез не снижался благодаря высокой мезофильной проводимости. Расчеты показали (Соколов и др., 2019а), что днем наблюдалась дегидратация ствола, запасы влаги в котором уменьшались. Обезвоживание дерева в солнечные дни постепенно нарастало (рис. 10г). Далее, в ходе кратковременного потепления на 14-й день (до 25°С) и увеличения VPD до 2 кПа наблюдался максимальный за все время уровень транспирации (1.2 ммоль м^–2 с^–1), который был обусловлен большим VPD. При этом устьичная проводимость снизилась на 10–15%. В целом ряде исследований последних десятилетий установлена обратная связь устьичной проводимости как с VPD, так и с WD (Warren, 2008; Campany et al., 2016; Stangl et al., 2019). Это подтверждают и наши предыдущие исследования (Соколов и др., 2019а)

В следующие дни дерево использовало низкий уровень VPD в утренние часы для максимального раскрытия устьиц. При этом диффузия CO₂ через устьица была значительной при минимальных потерях влаги. Это отразилось на 16-й день в величине WUE, возросшей на 80%.

А прохладную, преимущественно пасмурную погоду начала исследуемого периода дефицит воды в дереве (WD), был в целом ниже среднего уровня (1–2 л). В последующий солнечный период с более высокими температурами воздуха WD возрос до 4 л. Слабое снижение величин WUE было обусловлено последовательно возрастающим WPD, приводящим к небольшому росту транспирации, несмотря на уменьшение устьичной проводимости. Далее после нескольких дней с дождями водный баланс дерева пришел в норму (WD < 0.3 л) и величины WUE сохраняли относительное постоянство около среднего уровня.

Нужно отметить, что дегидратация ствола на протяжении всего июля оказывала значительное влияние на устьичную проводимость, что обусловило, как правило, ее невысокие значения на протяжении дня. Лишь в дождливые периоды при минимальных значениях VPD устьичная проводимости была высокой, что обеспечивало относительно высокий уровень поглощения CO₂ даже при низкой инсоляции. Именно этим объясняются очень большие выбросы WUE на 26–28-й день – до 85 мкмоль CO₂ (ммоль Н₂О)^–1.

Используя средние за период моделирования среднесуточные значения фотосинтеза – 2.52 мкмоль м^–2 с^–1, а также транспирации – 0.196 ммоль м^–2 с^–1 можно оценить величины поглощения CO₂ и испарения Н₂О на 1 га древостоя сосняка черничного.

Измерения, проведенные на модельном дереве 45–55-летнего возраста (высота 14.7 м, диаметр 16 см, протяженность кроны 8 м) показали, что площадь всей хвои составила 45 ± 5 м². При этом площадь световой хвои и хвои 3-го года и выше составляла 33% и 35% соответственно. Используя скорости транспирации (данные Б.М. Веселкова), связанные с освещенностью и возрастом, можно привести площадь всей хвои кроны к площади хвои одного типа, например, освещенной 1–2-летнего возраста. То есть, если бы вся крона дерева состояла из освещенной хвои 1–2-летнего возраста, то именно данная площадь обеспечила бы такую скорость транспирационного расхода влаги, как и реальное наше дерево. Проведенный по нескольким методикам расчет (Кайбияйнен, 1990) показал, что приведенная площадь хвои всей кроны к площади 1–2-летней световой хвои имела величину 32 ± 5 м². Все эти количественные расчеты не очень точны, есть ошибки в интегрировании расхода влаги побегами на всю крону.

Для пересчета на гектар воспользуемся данными измерений параметров модельных деревьев, полученных при изучении таксационных характеристик нашего древостоя (Казмиров и др., 1977). В сосняке черничном свежем вся хвоя составляет 4% от общей фитомассы всех деревьев – 47.6 ц/га. С учетом соотношений, полученных для нашего модельного дерева, поглощать углекислоту как 1–2-летняя хвоя будут примерно 33.3 ц/га. Достаточно простой расчет показал, что за выбранный период 1 гектар леса поглотил 14.6 т, при этом расход воды составил 464 т. Используя наши данные о распределении месячных сумм СО₂-газообмена и транспирации побегов сосны (рис. 11), можно для гектара леса за вегетацию дать приблизительную оценку поглощения CO₂ – 34 т, и испарения Н₂О – 1351 т. Погрешность этих оценок на уровне 20–25%.

Рис. 11.

Распределение месячных сумм транспирации (1) и фотосинтеза (2) побегов в верхней части кроны экспериментальной сосны.

За месяцы вегетации выпало 92 мм осадков (Кайбияйнен, 1990), на 1 га – 920 м³ воды. По нашим расчетам гектар леса испарил несколько больше. Но нужно учесть, что к моменту начала вегетации в почве содержались значительные запасы талых вод – около 1000 тонн/га. Таким образом, величина испарения с 1 га в целом соответствует литературным данным.

ЗАКЛЮЧЕНИЕ

Объединение знаний о процессах фиксации СО₂ растением, (раздел “Роль внешних факторов среды и внутренних параметров в моделях фотосинтеза на уровне листа”), с экспериментальными данными (раздел “Объект исследования и экспериментальные данные”) на основе методов математического моделирования (раздел “Применение технологии сбалансированной идентификации для анализа моделей фотосинтеза”) позволило шаг за шагом модифицировать (усложнять) описание объекта и в результате выбрать модель, сложность которой соответствует количеству и качеству экспериментального материала. Использование УКП позволило оценить динамику концентрацию CO₂ в хлоропластах и зависимость от нее ассимиляции СО₂. Однако надежно рассчитать проводимость мезофилла, как функции внешних условий не удалось. Возможно, проведение дополнительных активных экспериментов с целью получения УКП в более широком диапазоне переменных и привлечение дополнительных знаний о дневном дыхании (фотодыхании) позволит продвинуться дальше.

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

Работа была выполнена с использованием оборудования центра коллективного пользования “Комплекс моделирования и обработки данных исследовательских установок мега-класса” НИЦ “Курчатовский институт”, http://ckp.nrcki.ru/.

В расчетах использовался сервис решения наборов независимых задач оптимизации https://optmod. distcomp.org/apps/vladimirv/solve-set-opt-probs. Это один из нескольких сервисов оптимизационного моделирования https://optmod.distcomp.org на платформе Эверест, http://everest.distcomp.org/ (Sukhoroslov O. et al., 2015).

Работа выполнена сотрудниками ИППИ РАН при финансовой поддержке РНФ проект 16-11-10352 (разд. “Применение технологии сбалансированной идентификации для анализа моделей фотосинтеза” и “Результаты моделирования и их обсуждение”). Финансовое обеспечение исследований сотрудниками Института леса КарНЦ РАН осуществлялось из средств федерального бюджета на выполнение государственного задания КарНЦ РАН (0220-2014-0010) (разд. “Роль внешних факторов среды и внутренних параметров в моделях фотосинтеза на уровне листа” и “Объект исследования и экспериментальные данные”).

ПРИЛОЖЕНИЕ: ФРАГМЕНТ ФАЙЛА ЗАДАНИЯ ДЛЯ МОДЕЛИ МЕЗОФИЛЬНОЙ ПРОВОДИМОСТИ

Ключевое слово “Var:” означает неизвестную функцию-переменную, которую необходимо идентифицировать (определить), “Param:” – известную (заданную) функцию, “EQ:” – уравнение. Знак “#” означает комментарий

Var:   P(t)   # фотосинтез

Param: Q(t)   # радиация

  T(t)   # температура

  h(t)   # время суток

Param: Ca(h)   # концентрация CO₂ в атмосфере от времени суток h

Var:   Ci(t ≥ 64  # концентрация CO₂ в межклетниках

  Cc(t) ≥ 64  # концентрация CO₂ в хлоропластах

   A(t) ≥ 0   # ассимиляция CO₂

   Aqt(Q, T) ≥ 0 # зависимость ассимиляции от Q и T

  Acc(Cc) ≥ 0 # зависимость ассимиляции от C_c

   R(h) ∈ [0, 0.2]; R(0) = R(24) # дыхание от времени суток

Param: Gs(t)    # устьичная проводимость

Var:  Gm ∊ [0, 0.3] # проводимость мезофилла

EQ:  P = (Ca(h) – Ci) Gs × 1.0816 # закон Фика для устьичной проводимости

  P = (Ci – Cc)Gm × 1.0816 # закон Фика для мезофилла

  P = A – R (h) # фотосинтез = ассимиляция – дыхание

  A = Aqt(Q, T)Acc(Cc) # ассимиляция

   Aqt(0, T) = 0 # при Q = 0 ассимиляция не идет

  Acc (360) = 1 # нормировка

  d/dgCc(Acc(Cc)) ≥ 0 # Acc не убывает при росте концентрации C_c

  d/dgQ(Aqt(Q, T)) ≥ 0 # Aqt не убывает при росте концентрации Q