A Simple Harmonic Model for FAPAR Temporal Dynamics in the Wetlands of the Volga-Akhtuba Floodplain

Abstract

The paper reports a technique used to construct a reference time series for the fraction of absorbed photosynthetically-active radiation (FAPAR) based on remotely-sensed data in the largest Russian arid wetland territory. For the arid Volga-Akhtuba wetlands, FAPAR appears to be an informative spectral index for estimating plant cover health and its seasonal and annual dynamics. Since FAPAR algorithms are developed for multiple satellite sensors, all FAPAR-based models are suggested to be universal and useful for future studies and long-term monitoring of plant cover, particularly in wetlands. The model developed in the present work for FAPAR temporal dynamics clearly reflects the field-observed seasonal and annual changes of plant cover in the Volga-Akhtuba floodplain wetlands. Various types of wetland plant communities were categorized by the specific parameters of the model seasonal vegetation curve. In addition, the values derived from the model function allow quantitative estimation of wetland plant cover health. This information is particularly important for the Volga-Akhtuba floodplain, because its hydrological regime is regulated by the Volzhskaya hydropower plant. The ecosystem is extremely fragile and sensitive to human impact, and wetland plant cover health is a key indicator of regulatory efficiency. The present study is another step towards developing a methodology focused on arid wetland vegetation monitoring and conservation of its biodiversity and natural conditions.

1. Introduction

1.1. Background

The Volga-Akhtuba floodplain extends for more than 400 km along the Volga River, from the city of Volgograd and southeast to the Volga delta and the Caspian Sea. Its unique wetland ecosystems are the foundation of local economies, including electricity generation, agriculture, livestock production, fishery, chemical industry, recreational tourism, and so forth. Surrounded by desert and dry steppes, the floodplain is a huge oasis of human activity within three administrative territories of Russia: the Volgograd region, the Astrakhan region and the Republic of Kalmykia. Ecosystem monitoring along with the regulation of the hydrological regime in the area allows understanding, predicting and mitigating various natural and human impacts on ecosystems, particularly those associated with hydropower dam water discharge. These actions are crucial to preserving the development of the region [1]. In recent years, accelerating desertification processes have been observed in the Volga-Akhtuba floodplain [2]. Local wetlands are extremely sensitive to the hydrological regime in the arid environment, and their monitoring is essential for biodiversity conservation [3,4,5,6,7,8,9,10,11]. Because of the enormous size of the area under study (over $9000{\mathrm{km}}^2$), remote sensing techniques represent the only source of an unbiased and balanced broad view of the problem. In this framework, the vegetation activity is not only a significant factor in the local economy (e.g., in agriculture, livestock production), but it is also a good indicator of the integral condition of ecosystems [12,13,14,15,16,17,18,19,20,21].

Since the 1930s, studies of plant cover in the Volga-Akhtuba floodplain have been performed locally in the field using geobotanical methods [17,18,19,20] and hydrological research [5,6,22]. The use of remotely-sensed data only began in recent years. Research on the vegetation index dynamics using Landsat data has been done in the test areas of the Volga-Akhtuba floodplain mostly for the summer low-water season [2,23]. However, Landsat time series applicability is methodologically limited because of the temporal resolution of Landsat data, which does not allow direct comparison of the behavior of plant cover for different years or for different sources of remotely-sensed data.

1.2. FAPAR Index

The fraction of absorbed photosynthetically-active radiation (FAPAR) is one of the fundamental terrestrial state variables in the context of global change sciences. FAPAR is recognized as an essential climate variable by the Global Climate Observing System [24]. FAPAR represents the fraction of photosynthetically-active radiation absorbed by plant cover for photosynthesis. This index includes reflectance values in the red and infrared bands because green plants strongly absorb solar radiation in the red spectral region owing to photosynthesis and strongly reflect and scatter electromagnetic waves in the near infrared. Thus, the FAPAR index corresponds only to the living elements of the canopy that are green. Reflectance in the blue band is used to decontaminate the red and the near-infrared bands from atmospheric effects [21,25,26,27,28].

In this framework, FAPAR is an informative index for estimating Volga-Akhtuba wetland plant cover health and its seasonal and annual dynamics. Since FAPAR algorithms are developed for multiple sensors, we suggest that FAPAR-based models should be universal and useful for future studies and long-term monitoring of plant cover, particularly in wetlands [12,21,29]. Many studies on FAPAR dynamics in various plant ecosystems, including estimation of plant cover health, have been done in different areas and types of plant ecosystems worldwide. This index has been shown to be a reliable indicator of the state of a plant ecosystem [12,25,26,29]. However, FAPAR has never been applied in the Volga-Akhtuba floodplain.

1.3. FAPAR Temporal Dynamics

When analyzing the temporal dynamics of FAPAR for a particular year, a researcher often needs a reference model that serves as a “normal”, “regular” curve for examining any deviations. This reference model should include a time-continuous deterministic formula that provides a reference value of FAPAR for any point in time. The current work is aimed at presenting a simple version of such a model that requires no additional information beyond some FAPAR time series (possibly partial) for one or more reference growing seasons for calibration.

The model is relevant for geographic areas with distinct non-growing and growing seasons and particularly for arid wetlands. It is based on several physical and mathematical considerations that seem to be very well-suited for monitoring vegetation activity. These considerations were adjusted and tested, using remote-sensing and field data for several land sites in the Volga-Akhtuba floodplain. The model also provides several values that can be used to classify plant cover based on its yearly vegetation activity pattern. These indicator values may also be useful for comparing land sites with different environmental conditions and vegetation types. Other studies have applied similar approaches [30].

1.4. Test Sites

A detailed overview of test sites used in this study is provided in

Table 1. List of test sites.

Ref. No.	Site Name/Land Type	Location (N Latitude, E Longitude)	Prevailing Plant Cover Type
1	Pahotniy runnel vicinage	$48°{42}^{\prime }44.{9498}^{″}$ $44°{43}^{\prime }12.{6910}^{″}$	Woodlands (oakeries, black poplar and ash communities), abandoned agricultural lands
2	Zamora/dried lake vicinage	$48°{28}^{\prime }29.{4686}^{″}$ $44°{59}^{\prime }1.{2918}^{″}$	Woodlands (black poplar and ash communities), secondary feather grass-sagebrush steppe
3	Chichera lake vicinage	$48°{33}^{\prime }39.{7494}^{″}$ $45°8^{\prime }23.{6545}^{″}$	Steppe meadows and cereal meadows
4	Kaloshi lake vicinage	$47°{59}^{\prime }9.{0350}^{″}$ $46°{21}^{\prime }7.{2180}^{″}$	Cereal meadows
5	Tarpan runnel vicinage	$47°{45}^{\prime }0.{1975}^{″}$ $46°{36}^{\prime }37.{2794}^{″}$	Herb and cereal meadows temporally varying in proportions, partially overgrazed
6	Osochniy runnel vicinage	$47°{48}^{\prime }9.{5740}^{″}$ $46°{41}^{\prime }16.{6068}^{″}$	Cereal and Glycyrrhiza glabra meadows, partially overgrazed
7	Karasyachiy Ilmen vicinage	$47°{42}^{\prime }24.{3036}^{″}$ $46°{48}^{\prime }44.{6429}^{″}$	Cereal with Glycyrrhiza glabra meadows and steppe meadows
8	Maiorskoe lake vicinage	$47°{36}^{\prime }6.{7570}^{″}$ $46°{44}^{\prime }54.{1360}^{″}$	Herb and cereal meadows temporally varying in proportions
9	Akhtuba-Kriusha interfluve	$47°{14}^{\prime }46.{0665}^{″}$ $47°{17}^{\prime }53.{5291}^{″}$	Herb and cereal meadows with Glycyrrhiza glabra temporally varying in proportions
10	Sukhaya Akhtuba runnel vicinage	$47°{02}^{\prime }38.{5398}^{″}$ $47°{43}^{\prime }47.{7512}^{″}$	Herb and cereal meadows temporally varying in proportions
11	Bolkhuny desert	$48°7^{\prime }52.{5642}^{″}$ $46°{32}^{\prime }41.{8272}^{″}$	Desert, overgrazed
12	Leninsk steppe	$48°{44}^{\prime }30.{3269}^{″}$ $45°{18}^{\prime }30.{5118}^{″}$	Dry steppe
13	Sukhodol crop fields	$48°{36}^{\prime }42.{285}^{″}$ $44°{55}^{\prime }14.{4912}^{″}$	Woodlands (ash and oleaster communities), settlements and abandoned agricultural lands
14	Seroglazovo desert	$47°{08}^{\prime }11.{2606}^{″}$ $48°{06}^{\prime }8.{1519}^{″}$	Desert
15	Proran runnel vicinage	$48°{39}^{\prime }44.{0826}^{″}$ $45°9^{\prime }8.{7671}^{″}$	Marshes

, and their spatial distribution is shown in Figure 1. These 15 land sites have different types of plant cover (from desert to wetland), and each is approximately $6\times 6{\mathrm{km}}^2$ in size. Site details are listed in Table 1.

2. Methodology and Data

2.1. General Considerations

For further discussion, we will use dimensionless time t, normalized within the growing season:

$$0\le t=\frac{{[DOY-G_B]}_{mod365}}{{[G_E-G_B]}_{mod365}}\le 1$$

(1)

where $DOY$ is the current day of the year and $G_B$ and $G_E$ are the days at the beginning and end of the growing season, respectively. Modulo M operation ${[·]}_{modM}$ is reasonable for regions where $G_E$ is less than $G_B$, for example, in the Southern Hemisphere.

Because we are focused on the “normal” annual dynamics of vegetation activity, which are intended to serve as a reference, we assume that the state of the activity is close to some steady year-to-year life cycle, and nothing happens during the non-growing season. Mathematically speaking, this assumption means that the FAPAR function is periodic over normalized time with a period of 1. According to the Fourier series theory, any finite continuous periodic function can be represented as a converging sum (possibly with an infinite number of summands) over j of elementary harmonic terms $b_j\cos 2\pi jt+c_j\sin 2\pi jt$. A partial sum omits the tail part of the sum with term numbers over some n. This partial sum yields an approximate result with some predictable accuracy. The greater the number of terms, the more rapid oscillations are included in the model. In our experience, vegetation activity does not change too quickly, so oscillations shorter than $1/6$ of a growing season can be neglected in the model. Thus, considering n between 4 and 6, we can write the general form of the FAPAR model function:

$$f\left(t\right)=a_0+\sum _{j=1}^n\left(b_j\cos 2\pi jt+c_j\sin 2\pi jt\right)$$

(2)

where coefficients $a_0$, $b_1,\cdots ,b_n$, $c_1,\cdots ,c_n$ are determined in the following sections.

2.2. Additional Constraints

Presuming that the transitions between non-growing and growing seasons proceed normally, we consider vegetation activity to have a smooth fade-in at the beginning of the growing season and to fade out smoothly at the end. This consideration means a derivative of zero as the boundary conditions:

$${\left.\frac{\partial f}{\partial t}\right|}_{t=0}={\left.\frac{\partial f}{\partial t}\right|}_{t=1}=0,\frac{\partial f}{\partial t}=2\pi \sum _{j=1}^{n}j\left(b_j\sin 2\pi jt-c_j\cos 2\pi jt\right)$$

(3)

which give:

$$\sum _{j=1}^{n}j{c}_j=0\Rightarrow c_1=-\sum _{j=2}^{n}j{c}_j$$

(4)

Finally, from Equations (2) and (4), we obtain:

$$f\left(t\right)=a_0+\sum _{j=1}^n{b}_j\cos 2\pi jt+\sum _{j=2}^n{c}_j\left(\sin 2\pi jt-j\sin 2\pi t\right)$$

(5)

where $a_0$, $b_1,\cdots ,b_n$, $c_2,\cdots ,c_n$ are the best coefficients to approximate the reference FAPAR data, given a particular number of harmonics n. The number n is a trade-off between formula complexity, reference data approximation accuracy and reference data availability (because the number of reference points can be insufficient for determining a large number of coefficients).

2.3. Computing Harmonic Coefficients Using Reference Data and the Least Squares Method

Let some N reference FAPAR values $F_i$ be known for time instances $t_i$. The reference value is usually accompanied by some estimated standard deviation ${\sigma }_i$ that indicates the accuracy level of a particular measured FAPAR value $F_i$. The described dataset can either be available from special databases as ready-to-use quantities obtained from remotely-sensed data [31], or it can be manually derived from multispectral satellite imagery (e.g., using special software like ArcGIS, BEAM, GRASS GIS) [21,26,27,32]. If no estimated standard deviations are available for the whole dataset, all ${\sigma }_i$ are typically set equal (e.g., ${\sigma }_i=1$ for all i). If only a part of ${\sigma }_i$ is unknown, then unknown values can usually be set to the mean value of known ones.

Thus, we have a set of triplets $\{t_i,F_i,{\sigma }_i\}$, which are assumed to comply (to a certain extent) with the model given by Equation (5), so that for every i (from 1 to N):

$$F_i=a_0+\sum _{j=1}^n{b}_j\cos 2\pi j{t}_i+\sum _{j=2}^n{c}_j\left(\sin 2\pi j{t}_i-j\sin 2\pi t_i\right)+r_i,\mathbb{E}\left[r_i\right]=0\mathbb{E}\left[r_i^2\right]={\sigma }_i^2,$$

(6)

with $r_i$ being a (desirably small and generally unbiased) unknown deviation of an observed FAPAR value $F_i$ from the model, arising from both measurement uncertainties and the model inaccuracy and $\mathbb{E}[·]$ standing for mathematical expectation. The principle criterion for choosing $a_0$, $b_1,\cdots ,b_n$, $c_2,\cdots ,c_n$ is now the minimal sum of residuals squared, with respect to its accuracy indicators ${\sigma }_i$:

$$\sum _{i=1}^N\frac{r_i^2}{{\sigma }_i}=\sum _{i=1}^N\frac{1}{{\sigma }_i}{\left(a_0+\sum _{j=1}^n{b}_j\cos 2\pi j{t}_i+\sum _{j=2}^n{c}_j\left[\sin 2\pi j{t}_i-j\sin 2\pi t_i\right]-F_i\right)}^2\underset{a_0,b_j,c_j\in \mathbb{R}}{\stackrel{\phantom{a_0,b_j,c_j\in \mathbb{R}}}{\to }}min$$

(7)

Strictly speaking, the function minimized in the criterion given in Equation (7) and the particular weights, which are chosen as $1/{\sigma }_i$ in our case, both depend on the probability distribution of $r_i$. The distribution is usually unknown, and thus, it is conventionally assumed to be close to Gaussian and independent for each $r_i$. This assumption yields the sum of squares in Equation (7). With regard to weights, we should keep in mind that the estimated values of standard deviation ${\sigma }_i$ are usually not very accurate, so this quantity should not affect the solution too much. In fact, using weights in the form of $1/{\sigma }_i^2$ depends more strongly on the accuracy of ${\sigma }_i$ and yields apparently worse results with real data.

Note that the model on the right side of the Equation (6) is linear over unknowns $a_0$, $b_j$, $c_j$; therefore, it can be rewritten as a matrix equation in the following form:

$$\begin{array}{ccc}\hfill z& =& Hx+r,z={[F_1,\cdots ,F_N]}^T,x={[a_0,b_1,\cdots ,b_n,c_2,\cdots ,c_n]}^T,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill H& =& \left[\begin{array}{ccccccc}1& \cos 2\pi t_1& \cdots & \cos 2\pi n{t}_1& \sin 4\pi t_1-2\sin 2\pi t_1& \cdots & \sin 2\pi n{t}_1-n\sin 2\pi t_1\\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 1& \cos 2\pi t_N& \cdots & \cos 2\pi n{t}_N& \sin 4\pi t_N-2\sin 2\pi t_N& \cdots & \sin 2\pi n{t}_N-n\sin 2\pi t_N\end{array}\right]\hfill \end{array}$$

(9)

with ${[·]}^T$ meaning matrix transposition. In this notation for obtaining the optimal values of $a_0$, $b_j$, $c_j$, let us first introduce the goal function J:

$$J(a_0,b_1,\cdots ,b_n,c_2,\cdots ,c_n)={(Hx-z)}^{T}W(Hx-z),$$

(10)

$$W=\left[\begin{array}{ccccc}1/{\sigma }_1& 0& \cdots & 0& 0\\ 0& 1/{\sigma }_2& & 0& 0\\ ⋮& & \ddots & & ⋮\\ 0& 0& & 1/{\sigma }_{N-1}& 0\\ 0& 0& \cdots & 0& 1/{\sigma }_N\end{array}\right]$$

(11)

that is, of course, the same as in Equation (7). The goal function J then has its optimal value at point $\tilde{x}$ where:

$${\left.{\nabla }_{x}J\right|}_{x=\tilde{x}}={\left.\left[\frac{\partial J}{\partial a_0},\frac{\partial J}{\partial b_1},\cdots ,\frac{\partial J}{\partial b_n},\frac{\partial J}{\partial c_2},\cdots ,\frac{\partial J}{\partial c_n}\right]\right|}_{x=\tilde{x}}=0$$

(12)

The latter is a system of linear equations over x, which is resolved by:

$$\tilde{x}={\left(H^{T}WH\right)}^{-1}{H}^{T}Wz$$

(13)

The inverse matrix in Equation (13) exists if and only if there are $2n$ different $t_i$ values in the reference data. Therefore, reference data must contain FAPAR values for at least $2n$ different dates within the growing season. For $n=6$, there are 12 values, which are normally available. Typical FAPAR data series usually contain values for more than 30 different days in each year (10-day composites), and most of them fall within the growing season.

The desired coefficients for the harmonic model in Equation (6) are just the corresponding components of $\tilde{x}$ as noted in Equation (8), for example, $a_0={\tilde{x}}_1$, $b_1={\tilde{x}}_2$, $c_2={\tilde{x}}_{n+2}$, and so on. The coefficient $c_1$ is then calculated according to Equation (4). All coefficients for the harmonic sum in Equation (2) are now determined.

2.4. Accuracy Assessment

To validate the model quantitatively, three values were considered. Since the accuracy of reference FAPAR data varies significantly, direct comparison with the model is not always correct. Consequently, some weighting should be made when calculating mean values for differences between the model and data.

The first accuracy measure is a weighted mean of residuals between the model and actual data ($RWM$, residuals’ weighted mean):

$$RWM=\frac{1}{N}\sum _{i=1}^N\frac{w_i}{\overline{w}}\left[F_i-f\left(t_i\right)\right],w_i=\frac{1}{{\sigma }_i},\overline{w}=\frac{1}{N}\sum _{i=1}^N{w}_i$$

(14)

In real terms, this value represents the overall deviation of the model from the reference data, and it should be as small as possible. The least squares procedure guarantees that it is close to zero.

The second quantity is a weighted root-mean-square deviation of residuals ($RWD$, residuals’ weighted RMS deviation):

$$RWD=\sqrt{\frac{1}{N}\sum _{i=1}^N\frac{w_i^2}{{\overline{w}}^2}{\left[F_i-f\left(t_i\right)\right]}^2}$$

(15)

Smaller values of $RWD$ show that the appropriate model has been chosen to approximate the reference data series.

The last accuracy measure is the theoretically-estimated standard deviation ($ESD$) of the model predictions derived from a priori estimates of reference data and observability issues in the approximation problem. Since the uncertainty of the model arises from uncertainties in the estimation of harmonic coefficients and all of these values are multiplied by factors between 0 and 1, the upper bound for the standard deviation of the model is estimated to be:

$$ESD=\sqrt{{\tilde{\sigma }}_{a0}^2+\sum _{j=1}^n\left({\tilde{\sigma }}_{bj}^2+{\tilde{\sigma }}_{cj}^2\right)}$$

(16)

where ${\tilde{\sigma }}_{a0}$, ${\tilde{\sigma }}_{bj}$ and ${\tilde{\sigma }}_{cj}$ stand for estimated standard deviations of harmonic coefficients. They are defined from the least squares procedure as square roots of the diagonal elements of estimation error covariance matrix P, which in our case is the following:

$$P={\left(H^{T}WH\right)}^{-1}{H}^{T}H{\left(H^{T}WH\right)}^{-1}$$

(17)

Thereafter, ${\tilde{\sigma }}_{a0}=\sqrt{P_{11}}$, ${\tilde{\sigma }}_{b1}=\sqrt{P_{22}}$, ${\tilde{\sigma }}_{c2}=\sqrt{P_{n+2n+2}}$, and so on, and ${\tilde{\sigma }}_{c1}$ is computed according to Equation (4) as:

$${\tilde{\sigma }}_{c1}=\sqrt{\sum _{j=2}^n{\left(j{\tilde{\sigma }}_{cj}\right)}^2}$$

(18)

Thus, although $ESD$ is a purely theoretical value that highly depends on a priori quantities, it also reflects the uncertainty of harmonic coefficients based on the number of reference data points and their diversity in time within the growing season and on the estimated errors in reference FAPAR data.

2.5. Reference Data

FAPAR products used in the present study are based on SeaWiFS and MERIS satellite sensor data [26,27] downloaded from the European Commission Joint Research Centre (JRC) website [31]. The data are FAPAR JRC products, spatially averaged over $3\times 3$ pixel fragments (i.e., about $6\times 6{\mathrm{km}}^2$), 10-day time composites created according to the time composite algorithm developed by Pinty et al. [32]. Please refer to [21] for more details.

Test areas under study are listed in Table 1 and shown in Figure 1 on the map. They were chosen to have different types of relatively homogeneous intrinsic plant cover. Detailed descriptions of plant communities were made in the field at reference points within all test sites. The time span of the data used extends from 1997–2006. The years 1997–2005 encompassed a period of relatively high Volga runoff featuring healthy plant communities, so FAPAR measurements for this period were used as reference data to calibrate the model. In contrast, the year 2006 was extremely dry, and the vegetation activity and FAPAR dynamics for 2006 differ significantly from those of 1997–2005.

3. Results

3.1. FAPAR Yearly Temporal Dynamics Model

The complete methodology underlying the derivation of the model and its rationale are given in Section 2. Let F be the value of the fraction of absorbed photosynthetically-active radiation (FAPAR index). Its model (reference) value $F^{*}$ for a given day of the year is then calculated as follows:

$$F^{*}\left(DOY\right)\equiv f\left(t\right)=a_0+\sum _{j=1}^n\left(b_j\cos 2\pi jt+c_j\sin 2\pi jt\right),t=\frac{{[DOY-G_B]}_{mod365}}{{[G_E-G_B]}_{mod365}}\in [0,1]$$

(19)

where:

• $DOY$ is Day Of Year;
• $G_B$ is the beginning of a growing season; for the Volga-Akhtuba floodplain, we assume $G_B=60$ (i.e., the beginning of March);
• $G_E$ is the end of a growing season; for the Volga-Akhtuba floodplain, we accept $G_E=330$ (i.e., the end of November);
• t is the relative (dimensionless) time within the growing season;
• n is the total number of harmonics; currently, we suppose n is in the range from 4–6; for all 15 test sites, the sixth harmonic has an amplitude of less than $0.01$, which is definitely beyond the accuracy of the determination of the FAPAR index itself [26];
• j is the number of the harmonic terms;
• $a_0$, $b_j$, $c_j$ are the harmonic coefficients computed using an almost standard (except for some additional constraints and weighting described in Section 2) least squares methodology for the reference FAPAR time series;
• small f is used instead of $F^{*}$ just to clearly distinguish between the function of $DOY$ and the function of relative normalized time t;
• ${[·]}_{modM}$ means modulo M operation, which is reasonable for regions where $G_B$ is greater than $G_E$ (e.g., in the Southern Hemisphere); $M=365$ in our case because 365 is the longest time interval between two different days within a year.

3.2. Model-Based Indicators

The model-derived quantities listed below may be used for the classification of a land site based on its yearly vegetation activity pattern.

• $a_0$, a yearly mean FAPAR value.
• $AMP$, amplitude of the main harmonic and peak-to-peak amplitude (main harmonic overall span) $PP$:

  $$AMP=\sqrt{b_1^2+c_1^2},PP=2AMP$$

  (20)
• $MAXF$, maximum value of the FAPAR model function:

  $$MAXF=\underset{t\in [0,1]}{max}f\left(t\right)$$

  (21)
• $WA{V}_{FL-FH}$, weighted active vegetation period (expressed in days), which is defined as the integral of a special weighting function, $w_{FL-FH}\left(f\right)$, over time. This function has two parameters, $FL$ and $FH$, which are physically interpreted as the “FAPAR value of definitely low vegetation activity” and the “FAPAR value of definitely high vegetation activity,” respectively.

  $$WA{V}_{FL-FH}=(G_E-G_B)\underset{0}{\overset{1}{\int }}{w}_{FL-FH}\left(f\left(t\right)\right)dt,w_{FL-FH}\left(f\right)=\left\{\begin{array}{cc}1,\hfill & f\ge FH,\hfill \\ \frac{f-FL}{FH-FL},\hfill & FL\le f<FH,\hfill \\ 0,\hfill & f<FL.\hfill \end{array}\right.$$

  (22)
  For all days with a FAPAR value higher than $FH$, that is, the “high vegetation activity” threshold, $WA{V}_{FL-FH}$ just counts these days (since weighting function equals one). All days having a FAPAR value lower than $FL$, that is, the “low vegetation activity” threshold, are ignored (because the weighting function is zero). Finally, days with transition values of FAPAR between $FL$ and $FH$ are counted with weights proportional to the difference between the day’s FAPAR value and the threshold. This approach allows avoiding a leap between vegetation period estimates for FAPAR time series that are close to the threshold, but slightly lower or higher. Without weighting, given only one threshold parameter $FH$, a vegetation activity with FAPAR of $0.31$ on a 100-day period, for example, will be considered high for $FH=0.3$. However, the same vegetation activity with FAPAR $0.29$ will be considered low. Meanwhile, the accuracy of the FAPAR estimation is often off by several hundredths; thus, the practical difference between $0.31$ and $0.29$ is completely negligible. Therefore, weighting is required to obtain adequate estimates without leaps like those shown above. For the test sites under study, $FL$ and $FH$ are accepted as $0.2$ and $0.3$, respectively.

If $AMP>AM{P}_0$, then FAPAR dynamics are high enough to compute additional characteristics. If the main amplitude $AMP$ is not high enough, the rest of the indicators below are not informative and are poorly conditioned. Currently, $AM{P}_0=0.05$ is accepted.

• $PHASE$, the phase of the main harmonic,

  $$\phi =-\arctan \frac{c_1}{b_1}\in [0,2\pi ],PHASE={\left[G_B+\frac{\phi }{2\pi }{[G_E-G_B]}_{mod365}\right]}_{mod365}$$

  (23)
• $SHIR$, the secondary harmonics’ intensity ratio (i.e., the ratio between the cumulative amplitude of all secondary harmonics and the amplitude of the main harmonic). It represents how “wobbly” the reference FAPAR curve is.

  $$SHIR=\frac{\sqrt{{\displaystyle \sum _{j=2}^n}{b}_j^2+c_j^2}}{\sqrt{b_1^2+c_1^2}}$$

  (24)
• $DOYMax$, the day when the FAPAR model function reaches its maximum:

  $$DOYMax={\left[G_B+[arg\underset{t\in [0,1]}{max}f\left(t\right)]·{[G_E-G_B]}_{mod365}\right]}_{mod365}$$

  (25)

3.3. Application to Test Sites

The following charts (see Figure 2) depict several comparisons between the derived model and actual FAPAR data. FAPAR data from 1997–2005 correspond to the high Volga runoff period [23]. These values served as reference information. Summarizing the interpretation of the numerical results for test sites as compared to in situ data given in Section 4, from the ecological point of view, the model-derived quantities listed above seem to represent the ecological traits of particular land sites. Given the consistency between the model, actual FAPAR data and vegetation activity observed in the field, use of the model as a reference for analyzing and interpreting deviations from the model is reasonable. The complete set of charts and numerical results for all 15 test sites is given in Figures S1, S2 and S3 and Table S1 in supplementary materials.

4. Discussion

4.1. A Comparison of the Model with Reference Data and Vegetation Activity Observed in the Field

First of all, the plots clearly show that the model function in Equation (19) matches the actual reference data well. Therefore, the chosen model function is adequate, and a sufficient number of harmonics n is accepted. The order of Fourier partial series ($n=6$) might seem relatively high, but that is justified by the particular location of the area under study. The hydrological regime throughout the Volga-Akhtuba floodplain is regulated by a hydropower plant dam located near Volgograd (see the map in Figure 1). Spring water discharge from the power station dam leads to partial inundation of the Volga-Akhtuba floodplain. The length of inundation differs greatly by location [5,22], and in some places, the inundation can induce a local peak in vegetation activity in April or May, particularly in meadow areas. In addition, the local climate usually features a hot and dry summer, during which vegetation activity gradually declines. A relatively wet fall season follows, giving another potential local peak in vegetation activity in September or October. A clear example of these kinds of variations is shown in Figure 2c. These fluctuations imply that the model FAPAR function should have several harmonics.

By analyzing the model-derived indicators, one can categorize the type of plant cover. In particular, the marshes at Test Site #15 (Figure 2a) are characterized by the highest $AMP$, $MAXF$ and $DOYMax$ values and the lowest $PHASE$ among all test sites considered. This pattern of values indicates that over moistened areas, like the one under consideration, they have the highest vegetation activity. Since they are completely inundated during the special hydropower plant discharge (generally occurring in May–June [22]), $MAXF$ cannot be reached before the flooding ends, which results in a high $DOYMax$.

Field data also suggest that marsh vegetation begins its active growth only after the flooding ends because of the long local inundation period. The data also show that marsh vegetation sustains active growth because of the wet conditions throughout the summer. Composed mostly of sedges or common cane, the vegetation develops high projective coverage (up to 100%) and high biomass [16,18,19] in the middle of the growing season, which is maintained until fall. The model is consistent with these data.

The woodlands alternating with abandoned agricultural lands at Test Site #1 (Figure 2b) provide a different indicator pattern: relatively low $AMP$, $MAXF$ and $DOYMax$ values and high $WAV$ and $PHASE$ values. A remarkably long vegetation period ($WAV$) can be explained as follows. Weed vegetation on abandoned agricultural lands is relatively active in the spring and gradually weakens in the summer, whereas woodland vegetation activity increases at the beginning of summer. A relatively long growing season is the result. The beginning of active plant growth at Site #1 is explained by rising groundwater levels, which initiates tree leaf emergence and a growth burst in grass weed communities, which nonetheless never reach high projective coverage and biomass [16,18,19]. Woodlands occupy approximately 15% of Test Site #1 (which is quite high coverage for the Volga-Akhtuba floodplain), and because weed community vegetation ends by the middle of June or earlier, the overall vegetation activity never reaches high values. Because areas of that kind are not inundated, the maximum vegetation activity is observed during the flood period and is then kept at the average level.

Cereal meadows at Test Site #4 (Figure 2c), along with herb meadows at Sites #5, 8, 9 and 10 (Figures S1e, S2b, S2c and S2d, respectively, in supplementary materials), are the most common plant communities in the Volga-Akhtuba floodplain. These communities develop in inundated areas with a medium inundation period from about two weeks to a month [18,19,20,23]. In terms of FAPAR dynamics, cereal meadows provide the highest $SHIR$ values and low $AMP$ and $WAV$ values, but intermediate $MAXF$ values, which are shown in the chart. This pattern reflects the behavior of the plant communities observed in the field. The first secondary wave represents ephemeral spring vegetation that is quite active in areas of that kind. The main harmonic represents the primary peak in annual vegetation. Given cereal grass phenology, however, maximum vegetation activity never reaches high values and never lasts for very long, which is normal for such communities. However, again, the peak is reached no earlier than when flooding ends because of inundation. The autumn peak (third secondary wave) occurs during a relatively cool period with a traditional increase of precipitation. Herb meadows at Sites #5, 8, 9 and 10 represent a similar vegetation pattern, but typically with the third autumn local vegetation wave being significantly smoothed or even eliminated [20,23], which is also reflected by FAPAR model curves. Field observations prove that the cereal meadow vegetation curve has three waves, which correspond to spring pre-flood vegetation with ephemeral plant communities; a summer maximum (post-flood), when the highest coverage (up to 80%) and biomass are developed; and an autumn intensification, which occurs at the end of the growing season after the temperature decreases and the precipitation increases.

As expected, at most wetland sites, the FAPAR curve for the extremely dry growing season of 2006 (shown using a distinct line style) deviates significantly from the model, as in Figure 2c. This deviation, however, is insignificant for marshes near Leninsk (#15) (Figure 2a), which are always well inundated, including during the low hydropower plant water discharge in spring 2006. The other example comprises non-inundated sites in the upstream segment of the floodplain, such as Pahotniy runnel vicinage (#1) (Figure 2b), representing woodlands and abandoned agricultural lands, which are not typical for floodplain territory and depend on inundation to a lesser degree.

During the growing season of 2006, significant decreases in projective coverage and biomass were observed in water meadow plant communities [16,18,19,20]. These decreases may be explained by the drastic reduction in inundated areas and the shorter duration of inundation [22]. This explanation is also reflected by FAPAR dynamics if compared to the model function.

The test sites representing dry steppe and deserts outside the floodplain, including Site #14 (Figure 2d), are non-sensitive to the 2006 anomaly in terms of FAPAR dynamics. All values are minimal, with extremely low $AMP$ and $MAXF$, which differ greatly from the sites located inside the Volga-Akhtuba floodplain. The dry sites have plant communities with very low projective coverage [23]. Such communities are likely to increase their activity in response to precipitation (mostly in spring), and there is not any relation to the flood period. This pattern is also reflected by the FAPAR model function.

Comparing model-derived indicators like $AMP$, $SHIR$, $DOYMax$, and so forth, for different yearly vegetation activity patterns, one can see that the values differ considerably for various types of vegetation cover. Consequently, they can be used as objective indicators. For example, a marsh plant community in Figure 2a has high $AMP$ and $MAXF$ values. The more dynamic vegetation activity of cereal meadows in Figure 2c has a higher value of $SHIR$ than the plant communities in Figure 2b, which have weaker dynamics. Very low $MAXF$ values characterize the desert ecosystem in Figure 2d. This observable distinction between the indicators for significantly different vegetation activity patterns means that these indicators can be used for the classification of wetland plant cover and for the objective, quantitative estimation of its health.

It is important to note that all relations described above are rarely clear on an actual FAPAR curve for a particular year, because it is affected by many factors being in play during a particular period of time. The model curve has an advantage from averaging, while still retaining important features of the annual dynamics of vegetation activity.

4.2. Model Accuracy and Limitations

The main advantage of the derived model is that no additional data beyond the FAPAR time series (which is possibly incomplete) are needed, but this advantage is also its main limitation. If any factors like anomalous temperature, precipitation, runoff, and so forth, significantly affect photosynthesis, they can bias the model (if data are used as the reference) or divert FAPAR values from it. Thus, the model is not expected to match actual measurements exactly because real FAPAR values are always affected by deviations of the environment from its “normal” state; that is, even an ideal model would not precisely fit the data in our case. Nevertheless, the only way to assess the practical usefulness of the model is to compare it to the actual data. With this information in mind, three quantities were selected to serve as indicators of accuracy (please see Section 2.4 for their explicit definitions).

• $RWM$, weighted mean of residuals between the model and actual data, indicating whether the model is biased from the actual data. For all test sites, this value appeared to be effectively zero, as expected.
• $RWD$, weighted root-mean-square deviation of residuals between the model and actual data, showing the degree of its consistency. $RWD$ appeared to be $0.04$–$0.05$ in all cases, which corresponds well to the uncertainty of FAPAR observations themselves. Some evident examples of errors in reference FAPAR observations with the magnitude of $0.05$–$0.1$ are seen in Figure 2a,b,d in winter (between Days of Year 0 and 50).
• $ESD$, theoretically-estimated standard deviation of the model, based on the accuracy of reference data and observability properties of the estimation problem for harmonic coefficients. This measures the theoretical predictive power of the model for situations in which no significant anomalies are present, as compared to the reference dataset. Anyway, $ESD$ strongly depends on a priori estimated values for the uncertainty of FAPAR observations, which are typically not very accurate in the reference data. Thus, $ESD$ should be viewed a bit like an order-of-magnitude measure for what it is intended to indicate. These quantities turned out to be around $0.02$–$0.05$. In Figure 2, the estimated 2-σ corridors around the model are shown in yellow. They mean 95% probability to fit for the Gaussian distribution of errors.

In theory, because the model does not account for annual variations of significant factors like temperature, precipitation, flooding, and so forth, the straightforward comparison of $RWD$ and $ESD$ makes no sense. However, for a useful model, $RWD$ and $ESD$ should agree on the order of magnitude; that is, the model itself should adequately estimate its predictive power. This estimation was always the case for reference data.

It is also important to note that estimated standard deviations of the reference data come from the spatial root-mean-square deviation of FAPAR values throughout the specific territory. Thus, they mostly reflect not the real accuracy of FAPAR assessment, but its spatial heterogeneity (i.e., the plant cover variation inside the test site). These quantities directly affect the theoretically-estimated standard deviation of the model ($ESD$) and consequently have the same effect in terms of accuracy.

Vegetation dynamics, and particularly FAPAR, are obviously highly correlated with physical parameters of the environment (e.g., temperature, humidity, river runoff). This correlation explains some significant deviations in the actual data from the reference model. Despite their requirement for additional data on the physical state of the ecosystem, the models for the relation between vegetation dynamics and these physical parameters are of great interest. This research is planned for the future.

5. Conclusions

Applying conventional interpretation techniques that use remotely-sensed data to wetlands is often challenging because of the intrinsic heterogeneity and temporal variability of the area. These factors require researchers to develop new methods of data processing, and the described simple model for the fraction of absorbed photosynthetically-active radiation (FAPAR) yearly dynamics is one of them. The model is designed as an approximation by harmonic time series of a priori chosen order (currently $n=6$) with some natural restrictions, along with model-derived values to identify important traits of local ecosystems.

The developed reference model and model-derived quantities have predictive power, which is limited to years with environmental conditions (including hydroelectric power plant dam water discharge) similar to reference years. The model goes together with several accuracy and quality measures. It is spatially and temporally scalable; that is, it can be used from pixel-level to larger local areas with relatively homogenous ecosystems, and from yearly series to decades in time. The results of its application to 15 test sites in the Volga-Akhtuba floodplain for a period of nine years from 1997–2006 show the following.

• The designed model function matches well with actual data. All quality measures are consistent and acceptable. The predicted and validated accuracy of the model is at a level of several hundredths (typically $0.02$–$0.05$, compared to a full FAPAR range of one).
• Model-derived indicators can be used for plant cover and ecosystem type classification.
• Indicators can also serve as an objective measure of plant cover health and its seasonal and interannual dynamics. They allow quantitative description of anomalies, which is illustrated for the degradation of ecosystems in 2006, which was abnormally dry. The efficiency of runoff regulations can be estimated by comparison of model-derived quantities to reference values.