Performance of Singular Spectrum Analysis in Separating Seasonal and Fast Physiological Dynamics of Solar‐Induced Chlorophyll Fluorescence and PRI Optical Signals

High temporal resolution measurements of solar‐induced chlorophyll fluorescence (F) and the Photochemical Reflectance Index (PRI) encode vegetation functioning. However, these signals are modulated by time‐dependent processes. We tested the applicability of the Singular Spectrum Analysis (SSA) for disentangling fast components (physiology‐driven) and slow components (controlled by structural and biochemical properties) from PRI, far‐red F (F760), and far‐red apparent fluorescence yield (Fy∗760). The proof of concept was developed on spectral and flux time series simulated with the Soil Canopy Observation of Photochemistry and Energy fluxes (SCOPE) model. This allowed the evaluation of SSA decomposition against variables that are independent of physiology or are modified by it. Slow SSA‐decomposed components of PRI and Fy∗760 showed high correlations with the reference variables (R2 = 0.97 and 0.96, respectively). Fast SSA‐decomposed components of PRI and Fy∗760 were better related to the physiological reference variables than the original signals during periods when leaf area index (LAI) was above 1 m2 m−2. The method was also successfully applied to predict light‐use efficiency (LUE) from the fast SSA‐decomposed components of PRI (R2 = 0.70) and Fy∗760 (R2 = 0.68) when discarding data modeled with LAI < 1 m2 m−2 and short‐wave radiation Rin < 250 W m−2. The method was then tested on data acquired in a Mediterranean grassland. In this case, the fast SSA‐decomposed component of apparent LUE∗ showed a stronger correlation with the fast SSA‐decomposed component of Fy∗760 (R2 = 0.42) than with original Fy∗760 (R2 = 0.01). SSA‐based approach is a promising tool for decoupling physiological information from measurements acquired with automated proximal sensing systems.


Introduction
Remote and proximal sensing of vegetation are powerful tools for the exploitation of subtle signals related to plant physiology and photosynthesis. During the last two decades, the interest of the remote sensing community toward solar-induced chlorophyll fluorescence (F) and the Photochemical Reflectance Index (PRI) (Gamon et al., 1992) has increased due to evidence of the close relationship between plant physiological properties and these optical signals (Garbulsky et al., 2011;Meroni et al., 2009;Mohammed et al., 2019). A combination of different processes maintains plants' ability to deal with varying environmental conditions and stress factors. Under optimal light conditions and adequate water and nutrients supply, photochemical reactions, including CO 2 assimilation and electron transport, occur at high efficiency. However, when absorbed solar energy exceeds photosynthetic capacity, it risks damaging the reaction centers of the photosystems and must be emitted as F in the 650-850 nm spectral range or dissipated as heat as part of the non-photochemical quenching (NPQ) (Demmig-Adams & Adams, 1992). Since all three mechanisms (carbon-fixation, F and NPQ) compete for the same absorbed energy, the characterization of both F and NPQ might enable accurate inference of photosynthesis from optical signals (Frankenberg & Berry, 2018;Porcar-Castell et al., 2014). One of the NPQ thermal dissipation mechanisms is the de-epoxidation of the xanthophyll cycle pigments (Demmig-Adams, 1990;Niyogi et al., 1997). The excess energy leads to the interconversion of the xanthophyll cycle pigments, violaxanthin to antheraxanthin and then to zeaxanthin, providing a sink for this energy (Demmig-Adams & Adams, 1996;Pfündel & Bilger, 1994;Yamamoto, 1979). This process is quickly reversible, and zeaxanthin is converted back to violaxanthin under low light conditions and during the night. These conversions result in changes of leaf absorptance around 531 nm, which are detectable with PRI (Gamon et al., 1992;Garbulsky et al., 2011).
Based on the conceptual light-use efficiency framework introduced in Monteith (1972), gross primary production (GPP) is the product of the photosynthetically active radiation (PAR) absorbed by chlorophyll (aPAR Cab ) and the efficiency with which absorbed light can be used to fix atmospheric CO 2 (light-use efficiency, LUE): Cab aPAR is often expressed as a fraction of incoming PAR (faPAR Cab ), which represents a vegetation property independent of the magnitude of the incoming PAR.
Numerous studies have shown that F and PRI measured from a variety of tower-based (e.g., Cogliati et al., 2015;Kim et al., 2021;Wieneke et al., 2018;Xu et al., 2021), airborne (e.g., Middleton et al., 2017;Rascher et al., 2015;Rossini et al., 2015;Siegmann et al., 2021;Tagliabue et al., 2019), and spaceborne platforms (e.g., Middleton et al., 2016;Sun et al., 2018;Wang et al., 2020;Zhang et al., 2020) can successfully track variations in GPP and/or LUE. However, canopy-scale F and PRI are not exclusively driven by plant physiology. Structural and biochemical properties of vegetation varying at different timescales also affect F and PRI signals remotely observed. Phenological changes in pigment pools and canopy structure (e.g., leaf area index, LAI) induce a slow variation of these signals. Faster changes take place at diurnal scales driven by directional effects (sun position and illumination conditions) and physiological responses to environmental and stress conditions (e.g., high vapor pressure deficit, VPD).
At seasonal scale, PRI is sensitive to leaf properties, especially to slow pigment pool modifications responding to environmental factors such as sun exposure, aging, or chronic stress (Filella et al., 2009). These irreversible changes are termed constitutive components of PRI (Gamon & Berry, 2012). Facultative PRI components lead to fast and reversible physiological variations related to the de-epoxidation state (DEPS) of xanthophylls in response to illumination intensity (Demmig-Adams & Adams, 1992).
Following the adaption of the LUE model proposed by Lee et al. (2013) top of the canopy (TOC) F signal can be defined as: where ϕ f is the physiological F emission yield and f esc is the fraction of all F photons that escape from the canopy. Both faPAR Cab and f esc are controlled by canopy structure and leaf biochemical properties (Martini et al., 2019;Migliavacca et al., 2017), which usually vary at seasonal timescales. Fluorescence yield ϕ f , in BIRIUKOVA ET AL.
10.1029/2020JG006158 2 of 25 turn, directly responds to the energy partitioning in the photosynthetic machinery. Therefore, it represents a short-term variability (diurnal and sub-diurnal) driven by the modulation of the physiological status of plants.
The recent development of coupled photosynthesis and radiative transfer models (RTMs) (e.g., the Soil Canopy Observation of Photochemistry and Energy fluxes, SCOPE (van der Tol et al., 2009)), able to simulate physiologically driven F and PRI (Vilfan et al., 2018), allowed a better characterization of the physiological and biophysical components of these optical signals. However, the separation of the physiological and radiative sources of variability in long-term measurements of reflectance and F is still a challenge. This knowledge gap hampers the full exploitation of long-term data series collected by unattended spectroradiometric systems, such as those described in Porcar-Castell et al. (2015).
Alternatively, statistical approaches have been successfully used to decompose overlapping signals. For example, the Singular Spectrum Analysis (SSA) is a comprehensive methodology originally established by Broomhead and King (1986) and Fraedrich (1986) and later developed by Ghil et al. (2002) and Golyandina et al. (2001). SSA is a powerful tool for decomposition, reconstruction, and forecasting of climatic time series (Ghil et al., 2002;Plaut et al., 1995;Yiou et al., 1996). Also, SSA was successfully used for the characterization of the dynamics of eddy covariance (EC) ecosystem-atmosphere fluxes (Mahecha et al., 2007;Mahecha, Reichstein, Jung, et al., 2010;Wang et al., 2012) and for the evaluation of terrestrial biosphere and semi-empirical model performances at different timescales (Mahecha, Reichstein, Carvalhais, et al., 2010;Migliavacca et al., 2015). The main idea behind SSA is that time series can be described as a sum of superimposed subsignals, which can be extracted based on their characteristic scales of variability (Mahecha, Reichstein, Jung, et al., 2010). Unlike other methods, one-dimensional SSA is non-parametric, and therefore, it does not require prior information about the number and/or frequencies of periodicities, nor a model for trend (Golyandina et al., 2018).
In this work, we extend the use of SSA to decompose time series of optical signals related to the photosynthetic activity of plants, such as PRI, far-red fluorescence (F 760 ), and far-red apparent fluorescence yield (Fy * 760 = F 760 /PAR). We hypothesize that SSA can coherently decouple fast and slow components of the variability of these signals. We test this hypothesis by evaluating to what extent SSA fast components can be attributed to vegetation's physiological responses and SSA slow components to seasonal variations of structural and biochemical properties. We also evaluate the potential of SSA fast components to predict NPQ and LUE. First, we tested SSA on time series of half-hourly data simulated with the state-of-the-art model SCOPE. Then, we applied SSA to automated proximal sensing and EC time series acquired in a Mediterranean grassland.

SCOPE Simulations
SCOPE version 1.73 was used to simulate one year of TOC reflectance factors (R), F, and fluxes at half-hourly time steps (Biriukova et al., 2021). Variations of leaf absorptance between 500 and 570 nm induced by the inter-conversions of the xanthophyll cycle pigments were simulated by the leaf RTM Fluspect (Fluspect-CX) (Vilfan et al., 2018). Fluorescence radiance was simulated by SCOPE using the fluorescence emission spectra characterized from the FluoWat leaf clip measurements (fluorescence of photosystems I and II are not separated, SCOPE parameter "calc_PSI" = 0) (Vilfan et al., 2016) and an empirical fluorescence model (van der Tol et al., 2014) (SCOPE parameter "Fluorescence_model" = 0).
The distribution of absorbed light into competitive pathways is controlled by the rate coefficients (K), which express the probability of the different fates of the excitations . The rate constant for constitutive thermal dissipation that is present in dark-adapted plants (K d ) was defined as a function of leaf temperature (T) as K d = max (0.8738, 0.0301 · (T − 273.15) + 0.0773). The rate constant for fluorescence (K f ) was set to 0.05, and the rate constant for heat dissipation as part of NPQ (K n ) was defined as K n = K no · (1 + β) · x α /(β + x α ), where K no , α, and β are parameters of the empirical model (equal to 5.01, 1.93, 10, respectively). The rate constant of photochemistry (K p ) was set to 4. The degree of light saturation (x, used for computation of K n ), the steady-state fluorescence yield (F s ), and the fluorescence efficiency amplification factor (ϕ' f ) are the output parameters of SCOPE's fluorescence model. NPQ was computed as K n /(K f + K d ).

Journal of Geophysical Research: Biogeosciences
The steady-state fluorescence yield was computed as F s = F m * (1−ϕ p ), where light-adapted fluorescence yield is F m = K f /(K f + K d + K n ), and the photochemical yield is ϕ p = ϕ 0 p · J a /J e ; where ϕ 0 p is the photochemical yield under dark-adapted conditions, J a is the actual electron transport rate, and J e is the potential electron transport rate .
The simulations included temperature correction of the maximum carboxylation rate (V cmax , [μmol m −2 s −1 ]) (SCOPE parameter "apply_T_corr" = 1). Soil heat flux was defined as a constant fraction of soil net radiation (SCOPE parameter "soil_heat_method" = 2). SCOPE was parameterized to reproduce the spectral behavior of a Mediterranean grassland using vegetation properties derived from field measurements, as well as meteorological data from the research station of Majadas de Tiétar (39°56′24.68″N, 5°45′50.27″W) (Cáceres, Spain). The station is located in a typical Mediterranean savanna ecosystem dominated by herbaceous stratum constituted by grasses, forbs, and legumes. The site is characterized by a mean annual temperature of 16°C but a strong seasonality encompassing a wet season from November to May and dry summers (Perez-Priego et al., 2015). A complete description of the study site can be found in El-Madany et al. (2018) and Perez-Priego et al. (2017).
In order to realistically represent seasonal and intra-daily meteorological conditions, we used half-hourly observations of forcing meteorological variables measured in 2016  Soil moisture content (SM p , [%]) averaged from 4 sensors at 5 cm depth was used to modulate soil R in the brightness-shape-moisture (BSM) sub-model of SCOPE . The parameterization of SCOPE was defined according to Pacheco-Labrador et al. (2019) in the same site. Soil resistance for evaporation from the pore space (r ss , [s m −1 ]) was estimated from SM p as in Pacheco-Labrador et al. (2019).
Seasonal variability of LAI [m 2 m −2 ], leaf chlorophyll (C ab , [μg cm −2 ]) and carotenoid contents (C ca , [μg cm −2 ]) was simulated from time series of midday Normalized Difference Vegetation Index (NDVI) (Tucker, 1979) measured in 2016 by Decagon SRS sensor (Decagon Devices, Pullman, WA) at the study site (Luo et al., 2018). Therefore, these parameters varied daily. LAI was derived from an empirical relationship with NDVI (Martín et al., 2020) (roughly three times NDVI). C ab was predicted using a model fit from field spectral measurements and pigment content determined from destructive samples of 25 × 25 cm grass patches, sampled in several campaigns between 2017 and 2019 (Martín et al., 2020;Melendo-Vega et al., 2018). C ab was estimated as C ab = (0.007 − (0.0001/NDVI) · log (1 + (NDVI/0.0001))) · 4,443, while C ca was predicted as a function of C ab , according to the linear model C ca = 0.24 · C ab + 0.67 using field information from the same data set. The ratio of C ab to C ca of the simulated data set ranged from 3.06 in autumn to 3.78 in spring. Other parameters were kept constant during the simulation. Leaf angle distribution was assumed spherical, V cmax was set to 80 μmol m −2 s −1 , and the slope (m, [-]) of the Ball-Berry model (Collatz et al., 1991) was set to 10. These parameters were fixed to simplify the simulations. A constant diffuse to global radiation ratio was set to 20%. Moreover, the variables from the biochemical model of SCOPE (e.g., NPQ, ϕ' f , ϕ p ), which are not part of the default output of SCOPE, were also extracted. The same way as the model fluxes (e.g., GPP), these variables were computed as a weighted mean of all leaves inside the canopy, considering leaf angle distribution and the relative depth of each leaf, which determine the amount of radiation that each leaf receives (van der Tol et al., 2009). However, unlike the fluxes, biochemical variables were not scaled by LAI since they are leaf rather than canopy scale parameters.
PRI was computed as follows (Gamon et al., 1992): where R 531 is the reflectance factor of the xanthophyll-sensitive band at 531 nm and R 570 is the reflectance factor of the reference band at 570 nm. With this formulation, PRI values can vary between −1 and 1 and are directly proportional to NPQ.

Journal of Geophysical Research: Biogeosciences
To evaluate whether extracted SSA components were related either to the constitutive or the facultative drivers of variability, SCOPE model was run in two different modes featuring and excluding physiological effects on PRI and F (i.e., the effect of the xanthophyll cycle de-epoxidation on leaf absorption and fluorescence efficiency amplification factor (ϕ' f ) on F).  . For the definition of the variables simulated with SCOPE and modeled with SSA refer to Table 1.
Only daytime data were simulated. Any situation where R in ≤ 10 W m −2 or SZA ≥ 85° was considered night. Night-time values of different model outputs were linearly interpolated between the sunrise and sunset values. However, R could strongly change due to small variations in SZA when it is large. Therefore, sunrise and sunset were simulated differently to provide a smooth baseline. For these time steps, we set SZA = 85° and used the forcing meteorological variables of the following or the previous time step, respectively. R was calculated for R in = 10 W m −2 using the optical radiative transfer module of SCOPE only so that no physio-BIRIUKOVA ET AL.  PRI were identical. The rest of the model outputs were computed for R in = 0 W m −2 using the full model. At such low radiation levels, the effects of the xanthophyll cycle on R resulted negligible, and it was assumed that these simulations could be representative of dark night-time. Night-time interpolation is necessary since SSA requires continuous data. This gap-filling is reasonable for parameters that are not expected to strongly vary during the night, such as F or R-based variables, as is the case for PRI. These assumptions might not completely hold due, for example, to physiological recovery. Since PAR equals 0 at night, the night-time gaps in the series of  760 Fy were filled with the maximum daytime values selected in the moving window of 1 day. We simulated observational uncertainty by adding Gaussian noise to the time series of PRI, F 760 , and  760 Fy with 0 mean (μ), standard deviation (σ) equal to the 95% quantile of daytime data multiplied by 0.01. LUE [µmol CO 2 /µmol photons absorbed] was computed as the ratio of simulated GPP to aPAR Cab .

Singular Spectrum Analysis
SSA is one of several potential time series decomposition techniques. It was chosen here because it is highly data-adaptive and allows for decomposing phase-modulated signals. The method can be described in four steps: embedding, decomposition, grouping, and reconstruction (Golyandina & Korobeynikov, 2014) ( Figure 2).
Step 1 Embedding. The original time series Step 2 Decomposition. Singular value decomposition (SVD) leads to elementary matrices of rank 1: where d is a rank of X. Each elementary matrix i X is defined by the eigentriple: The eigentriple consists of a singular value  i , the left eigenvector i U and the right eigenvector i V . The singular values of eigentriples are proportional to the fraction of explained variance corresponding to each eigentriple.
Step 3 Grouping. The grouping is performed by choosing the sets of eigentriples (eigentriple grouping) so that each set corresponds to an identifiable series component. The grouping procedure partitions the set of indices   1,...,d into m disjoint subsets 1 ,..., m I I . The result of this step is the grouped matrix decomposition of the expansion (6): Step 4 Reconstruction. In the last step, each matrix of the grouped decomposition (7) is transformed into a new series of length N by diagonal averaging. As a result, the initial time series   1 ,..., N f f is decomposed into a sum of m reconstructed series: Time series decomposition was implemented using R-packages Rssa (Golyandina et al., 2018;Korobeynikov et al., 2017) and spectral.methods (Buttlar et al., 2014). There are two parameters in SSA, which the analyst must set: the window length (L) and grouping of the eigentriples. The choice of L is dependent on the characteristics of the subsignals to be extracted. In general, L ≤ N/2, and the higher the L, the more detailed the BIRIUKOVA ET AL.

10.1029/2020JG006158
6 of 25 decomposition is. For the identification of the trend, L should be large enough to be separable from periodic oscillations and noise. For the extraction of a periodic component with a period T, it is advisable to have L proportional to T. The periods of the harmonic components of the time series can be identified with the periodogram.
The grouping can be done manually by analyzing the graphs of eigenvectors and their frequencies or automatically (Golyandina & Zhigljavsky, 2013). In this study, we used the automatic grouping implemented in spectral.methods R-package (Buttlar, 2015), which groups SSA components based on their common features. This method measures the commonality of components by means of the weighted correlations between the components: if the weighted correlation is high, then the corresponding components have similar behaviors and should be included in one group (Golyandina et al., 2018).
The whole algorithm was run stepwise for each frequency interval (which are specified by "borders.wl" argument, Table S1). This allows adapting L for a particular frequency bin to be extracted. The choice of frequency bands is subjective and was determined here based on the time series length and temporal resolution (30 min ). In particular, we defined the following classes: long-term or seasonal (2 weeks-1 year), diurnal (7 h-2 weeks), sub-diurnal (30 min-7 h). For each frequency bin, the chosen window length L was: 2 months, 1 week, and 1 day. The choice of window length was supported by the SSA theory (i.e., for extracting the long-term component, L was chosen large enough to be separable from the periodic component, while for the extraction of diurnal oscillations, L was proportional to the period of 1 day). The data associated with the output of SSA analysis can be found in Biriukova et al. (2021 Table S1.

Field Data Acquisition and Processing
Spectral measurements (hereafter denoted with subscript "m") were acquired in Majadas de Tiétar with the high-resolution fluorescence box (FloX) device (JB Hyperspectral Devices UG, Germany), specifically designed for retrieving F and vegetation indices in visible (VIS) and near-infrared (NIR) domains. The FloX system contains two spectrometers -QE Pro (wavelength range of 650-813 nm, spectral sampling interval (SSI) of 0.15 nm and full width at half maximum (FWHM) of 0.3 nm), and Flame (wavelength range of 340-1,020 nm, SSI = 0.65 nm, FWHM = 1.5 nm) (Ocean Optics, USA). The up-welling radiance was measured with the 25° field of view (FOV) of fiber optics. Down-welling irradiance was measured using hemispherical cosine receptors. The FloX was installed on a radiometric tower at 10 m height, which resulted in a radiometric footprint of approximately 4.4 m in diameter. An automatic rotating arm allows sequential measurements over a tree and herbaceous stratum with an interval of 15 min. In this work, we use measurements acquired over the herbaceous stratum only. Measurements were collected continuously from 6 a.m. to 7 p.m. (UTC time) with an average interval of 2 min. The data set covers the period from March 5, 2017 to October 24, 2018, with gaps associated with the instrument's technical issues from mid-June 2017 to mid-August 2017 and from mid-January 2018 to mid-March 2018.
CO 2 fluxes between vegetation and atmosphere were measured with the EC technique (e.g., Baldocchi et al., 1996) (hereafter denoted with subscript "m"). EC system consists of an infrared gas analyzer LI-COR LI7200 (LI-COR Inc, Lincoln NE) and a three-dimensional sonic anemometer Gill R3-50 (Gill Instruments Ltd., Lymington under SZA > 70° were discarded. Second, we applied linear interpolation to fill missing data points for each day that included more than five valid daytime observations. Third, large gaps in the time series were filled using the SSA gap-filling approach implemented in the function gapfillSSA of the spectral.methods R-package (Buttlar et al., 2014) Forth, the night-time baseline was recomputed for each variable by filling it with the maximum or minimum (depending on the variable) daytime values selected in the moving window of 1 day. Last, Gaussian noise was added to night-time baselines with μ = 0 and σ equal to the 95% quantile of daytime data multiplied by 0.01.
Since aPAR Cab was not available from field observations, the apparent  m LUE was estimated as the ratio of m GPP to m PAR . Light-use efficiency represents a combination of different components -green canopy structure, light absorption, and physiology (Gitelson & Gamon, 2015). For this reason, we applied the SSA decomposition on  m LUE to decouple the influence of canopy biochemical and structural properties (i.e., faPAR Cab ). In this case, we consider the SSA fast varying component extracted from , and  m LUE were decomposed into the same frequency classes as the SCOPE-simulated variables (Section 2.2): long-term or seasonal (2 weeks-1 year), diurnal (7 h-2 weeks), sub-diurnal (30 min-7 h). Parameters of the filterTSeriesSSA function used for the decomposition are also identical and reported in Table S1.

Seasonal Cycles of SCOPE-Simulated Variables
The Mediterranean climate is characterized by a strong seasonality, driven mainly by radiation and precipitations (e.g., El-Madany et al., 2018). The rainy period occurs during late fall and early spring, and the dry season in summer extends to early fall (Figure 3b). Simulated biophysical parameters, fluxes, and spectral variables are coherent with the typical phenology of the grassland at the site (Luo et al., 2020). According to the models based on NDVI observations, during the green-up period from fall to winter, simulated LAI and C ab increased from 0.5 to 2.5 m 2 m −2 and from 7.5 to 25 µg cm −2 , respectively ( Figure 3d). This variability is coherent with expected phenology and with the variability of SM p (Figure 3b). The peak of the growing season in Majadas de Tiétar occurs in spring (Luo et al., 2020). This is reproduced by simulated GPP and scope 760

F
, which featured maximum values at the beginning of May (Figures 3e, 3i and 3j). Early summer is characterized by a dry-down period, a transition to the hottest and driest season where photosynthesis is strongly inhibited by water stress. Summer T air reaches 40°C and VPD values of 75 hPa (Figures 3a and 3b), which is also represented in the simulations. at short temporal scales when canopy structural parameters (LAI) do not vary significantly (Figures 4a and 4b) is scaled with the fluorescence efficiency factor predicted by the biochemical model according to the way aPAR Cab is split in the photosynthetic machinery. Therefore,  scope 760 F can be attributed to the physiological regulation of fluorescence efficiency (Figures 4c and 4d).

Journal of Geophysical Research: Biogeosciences
For the decomposition analysis, we assume that the scaling of canopy-level parameters  scope PRI and  scope 760 F to leaf-level NPQ and ϕ' f is unnecessary if the relationship between the decomposed SSA components and NPQ or ϕ' f is evaluated for different LAI classes.

PRI
(R 2 = 0.97) (Figure 6a). This confirms that SSA was able to distinguish the long-term variability of   (Figures 7d-7f). The same results are found for noiseless simulations, which suggest that SSA cannot extract  scope (diurnal + sub-diurnal) has a large potential to track  scope 760 / F PAR better than scope 760 Fy (R 2 > 0.75), but that it is reduced by the presence of noise in the signal (Figure 7i). To ascertain that ssa 760

Fy
(diurnal + sub-diurnal) is related to vegetation physiology, we applied the same SSA decomposition to scope 0,760

Fy
. In this case, the components of the same frequency showed a lower correlation with  scope 760 / F PAR (R 2 =0.48), but the relationship was still significant ( Figure S1).

ssa and Light-Use Efficiency -Model-Based Scenario
To evaluate whether SSA-decomposed fast components of optical signals are better related to the physiological response of vegetation than the original undecomposed variables, we compared their relationships with LUE ( Figure 8). Also, according to the results presented in Figures 7c, 7f and 7i, we considered the limitations of SSA to predict fast components of different spectral variables at low LAI values and evaluated the relationships for data with LAI above 1 m 2 m −2 (Figure 9).   (Figure 8h). Additionally, these relationships are different for high and low light conditions (and, therefore, low NPQ), which makes unsuitable the evaluation of the linear models. When NPQ < 0.1 (which in the simulation occurs when R in < 250 W m −2 ), the relationship between ssa 760 Fy (diurnal + sub-diurnal) and  scope 760 / F PAR and LUE becomes negative ( Figure S2). In order to properly evaluate these relationships, we limited the analysis to the cases where  scope     (Figures 9b and 9f).

Performance of SSA in Separating Seasonal and Fast Dynamics of F and PRI
With the increasing availability of high temporal resolution optical data collected simultaneously with CO 2 fluxes at EC stations (Balzarolo et al., 2011;Gamon, et al., 2006), it is pivotal to interpret the information provided by these data sets adequately. PRI and F are related to photosynthetic activity but confounded by biophysical parameters which vary at different timescales (Gamon & Berry, 2012). To our knowledge, there were no previous attempts to use time series decomposition as a tool to disentangle physiological information from these signals.
Our approach shows that SSA decomposition of scope xan

PRI
, scope 760 F , and scope 760 Fy simulated with SCOPE for a specific case study of Mediterranean grassland ecosystem allowed to separate slow and fast varying components with different levels of accuracy. Following the two-components concept of PRI variability -constitutive and facultative, introduced in Gamon and Berry (2012), we showed that these components could be successfully distinguished using the highly data-adaptive SSA technique. The decomposed slow variability of the total scope xan PRI showed a high correlation with constitutive scope 0 PRI (Figure 6a), which is highly correlated with modified red-edge normalized difference index (mNDI), sensitive to chlorophyll content (R 2 = 0.94) (Sims & Gamon, 2002). This result is similar to the strong relationships emerging between PRI of perfectly dark-adapted leaves and mNDI obtained in Hmimina et al. (2014 and Merlier et al. (2015) at both leaf and canopy scales. The facultative variability  scope PRI was well predicted by the fast component ssa xan PRI extracted with SSA and varying at diurnal and sub-diurnal timescales (Figure 7a).
The potential to capture fast dynamics of spectral signals related to vegetation physiology was assessed by relating ssa xan PRI (diurnal + sub-diurnal) to LUE. The overall high R 2 (Figure 7c) gradually increased from low to high LAI values. In the modeled data set, the summer season combined low LAI and reduced water availability with highly variable VPD (daily range of 0-75 hPA) ( Figure 3b) and NPQ values (daily range of 0-2.6). Under these severe conditions, there is a lack of PRI response to increasing PAR, which might explain the low correlation between ssa xan PRI (diurnal + sub-diurnal) and LUE found at this period of the year.
The analyses of noiseless data discarded that these weak correlations in these circumstances were induced by noise ( Figure 7c).
As shown in Figure 8a, total ssa xan PRI and LUE exhibit a non-linear relationship with variability in intercepts between LAI classes. The fast component of ssa xan PRI (diurnal + sub-diurnal), in turn, linearly correlated with LUE when the relationships were considered for the periods of similar LAI (Figures 8g and 9g). Previous works of Hmimina et al. (2014, evaluated the relationships of total and pigment-corrected PRI with LUE for short periods, excluding sources of significant LAI changes. We assessed the SSA decomposition BIRIUKOVA ET AL.

10.1029/2020JG006158
16 of 25 Despite the robustness to noise and the capability of SSA to capture the facultative component of PRI (Figures 7a-7c (Figures 10a and 10d). This might be explained by the fact that under the strong seasonal variability of the Mediterranean ecosystem (in terms of pigment pool and LAI), the PRI signal itself covaries with LUE in a way that, after all the uncertainties, SSA decomposition cannot offer BIRIUKOVA ET AL.

10.1029/2020JG006158
18 of 25  (diurnal + sub-diurnal); (i) R 2 for linear relationships (c), (f) aggregated by NDVI classes of equal size excluding and including noise. The data presented in the figure correspond to daytime data, SZA ≤ 80°, NDVI > 0. 4 (a, b, d, e, g, and h), and R in < 250 W m −2 (c, f, and i). significant advantage respect to the original signal. However, this could be a valuable tool in less dynamic ecosystems, such as evergreen forests, where the estimation of ΔPRI has provided additional access to LUE . Also, it must be considered that ssa m PRI (diurnal + sub-diurnal) may capture additional sources of variation. On the one hand, ssa xan PRI (diurnal + sub-diurnal) responds to diurnal and sub-diurnal variations of NPQ driven by instantaneous changes in light intensity (Krause & Weis, 1991) and modulated by VPD and T air (Demmig-Adams & Adams, 1992); on the other hand, these are mixed with directional effects (i.e., sun position and illumination conditions) imposed on diurnal cycles of PRI (Biriukova et al., 2020;Hall et al., 2008;Hilker et al., 2008). Both processes vary in the same frequency bin and only change in amplitude, and SSA might be unable to separate these two components. with PAR we removed the part of the variability attributed to irradiance and reduced the number of unknowns in the LUE model (Equation 2).
We hypothesized that by decomposing scope Fy was more sensitive to noise than the decomposition of PRI; however, in all the cases, the method resulted inaccurate at low LAI. This might imply the existence of a minimum  scope 760 F threshold for SSA to capture its fast variability. Discarding data modeled with LAI < 1 m 2 m −2 improved the estimation of  scope 760 / F PAR with ssa 760 Fy (diurnal + sub-diurnal) ( Figure S3). Despite the sensitivity to noise, when limited to a sufficient LAI range, SSA successfully decoupled physiology-related information from apparent fluorescence yield both in the modeled data set ssa (diurnal+sub-diurnal) (Figures 10f and 10i). The evaluation of this tool on simulated data provided a valuable insight into the limitations and strengths of SSA in terms of noise and variable space, both for the estimation of physiologically driven fast variability and the relationship of these estimates to variables of physiological interest (LUE).
The choice of SSA decomposition parameters plays an important role in the reconstruction of time series. The decomposition can be potentially improved by varying the "cut-offs" of the frequency bins (Table S1) in which we aim to extract the subsignals. Here we tested different frequency bins of varying length and borders and heuristically chose broader bins covering three timescales interesting from an ecological point of view -seasonal, diurnal, and sub-diurnal. Depending on the time series length and temporal resolution, the borders of the frequency bins of interest might be expanded (e.g., if sampling frequency is lower, there is no need to detect sub-diurnal components), or narrowed (e.g., in case of sampling frequency higher than 30 min would be interesting to have several bins per day to explore the dynamics of the fast physiological response at different timescales during a day).

Limitations and Applicability to Field Data
This study is a model-based proof of concept for decoupling slow dynamics related to biochemical composition and structural changes and fast physiological dynamics in TOC time series of F and PRI. We chose a BIRIUKOVA ET AL.

10.1029/2020JG006158
19 of 25 model-based study because the variables required to evaluate the proposed decomposition technique (i.e., a combination of canopy scale passive and leaf-level active measurements) are only sparsely available and prone to uncertainty related to the scaling from leaf-to canopy-level processes such as NPQ. In this context, the use of synthetic data generated by a state-of-the-art process-based model allowed evaluating whether SSA can disentangle the processes acting at different timescales as encoded in the model. This kind of model-based evaluation has been successfully used for various problems lacking suitable data sets (e.g., Nelson et al., 2018). However, it should be noted that SCOPE is not a dynamic model. This implies that the model does not capture some of the physiological responses (i.e., sustained NPQ) that can be observed in the field. When analyzing the links between decomposed PRI and NPQ, it should be considered that NPQ involves mechanisms operating at different temporal scales, such as reversible energy-dependent NPQ with overnight relaxation and sustained NPQ operating at longer timescales (Porcar-Castell, 2011). However, predicting NPQ based on PRI is only possible when reversible NPQ is dominating and sustained NPQ is insignificant (Alonso et al., 2017). In SCOPE simulations, we kept a constant rate of sustained NPQ, so that modeled NPQ only includes the effect of the xanthophyll cycle activation. Therefore, when applying SSA to field data, one should consider that decomposed fast components of PRI may not represent the total variability of NPQ.
Additional considerations are associated with the parameterization of SCOPE fluorescence module with data sets limited to few species (i.e., cotton), which may also affect the accuracy of the representation of the physiological response of grassland ecosystem. Nevertheless, we assume that for the purpose of the decomposition of slow and fast temporal dynamics of F and PRI with SSA SCOPE provided reasonable variabilities controlled by physiological and non-physiological vegetation properties.
Another limitation of the method arises in the case of structurally complex canopies. SCOPE model relies on the approach of SAIL, which represents a homogeneous canopy of randomly distributed leaves. This approach cannot reliably represent canopies with a strong geometrical component due to a heterogeneous 3D structure. PRI and fluorescence models have been implemented in 3D RTMs such as FluorFLIGHT (Hernández-Clemente et al., 2017) or DART (Gastellu-Etchegorry et al., 2017); however, these implementations mainly limit to the radiative transfer of these signals, and lack of a coherent coupling with photosynthetic efficiency at a comparable level of detail. Further research is needed in the case of structurally complex canopies, but current models do not yet allow such detailed analyses.
SSA decomposition requires high temporal resolution data, acquired with at least 1 h-30 min interval to track the fast physiological response of F and activation of the xanthophyll cycle, occurring at timescales of minutes after changes in light intensity (Müller et al., 2001). In general, the higher the temporal resolution of the time series, the more accurate extraction of fast varying physiological components can be achieved. With the expanding network of automated proximal sensing systems (Aasen et al., 2019;Cogliati et al., 2015), continuous and high-resolution time series of F and R (as measured in Majadas de Tiétar, Figure 10) become increasingly available.
These systems can be used to acquire time series of PAR, PRI, and F, as well as normalized vegetation indices informing on vegetation structure. Often, these systems are installed close to EC sites, which provide measurements of NEE, partitioned on GPP and respiration, and auxiliary abiotic variables (Rebmann et al., 2018), which would be helpful for the evaluation of SSA decomposition performance under different environmental conditions. Some EC sites also provide measurements of aPAR, but often this term has to be modeled from remote sensing observations, which is still challenging. Apparent LUE * (GPP/PAR) and fluorescence yield Fy * (F/PAR) are usually easier to obtain in comparison to metrics derived using aPAR or aPAR Cab . Gitelson and Gamon (2015) showed that LUE computed as GPP/PAR is the LUE estimate the most confounded by canopy structure among other LUE formulations. Contrarily, LUE computed as GPP/ aPAR Cab mostly depends on the physiological status of vegetation. aPAR can be estimated using several downward and upward-facing quantum sensors installed above and below the canopy (Inoue et al., 2008;Jenkins et al., 2007) or using automated observation systems based on LED sensors (Kim et al., 2019). aPAR Cab can be effectively retrieved from field spectral reflectance and transmittance measurements (Serrano et al., 2000) or approximated using vegetation greenness indices (e.g., NDVI). However, there is still no standard procedure to accurately estimate aPAR Cab from optical reflectance factors. In this study, we showed the applicability of the SSA decomposition on commonly available  (diurnal + sub-diurnal) to assess whether the fast component of 

760
Fy is related to physiological status of vegetation.
In addition to time series of optical variables, leaf-level pulse amplitude-modulated (PAM) measurements of NPQ, ϕ f , and ϕ p would be greatly beneficial for validating the decomposition results. For example, simultaneous installation of automated high spectral resolution devices (e.g., FloX system, JB Hyperspectral Devices UG, Germany) with micro-PAM (Atherton et al., 2016;Magney et al., 2017;Porcar-Castell et al., 2008) can provide a data set for the validation of the method. However, the availability of these coupled datasets is still very limited.
Since SSA requires data continuity, gap-filling should be applied to time series with missing values. For this purpose, SSA has also been successfully used as a gap-filling tool (Buttlar et al., 2014;Mahecha et al., 2007). The classical SSA algorithm was modified so that SSA components are estimated based on non-missing values only, and the values of the reconstructions are imputed to missing values (Golyandina & Osipov, 2007). For night-time data, a noisy baseline should be provided as well, for example, by linearly interpolating the last daytime observation of a day and the first daytime observation of the following day, or by applying a moving window to smooth a baseline. One of the reasons why SSA-decomposition of measured data presented here did not yield robust results in the case of m PRI (Figures 10a, 10d and 10g) might be long gaps in the time series (up to two months of missing data). In the future, the method should be tested on longer and more continuous time series.

Conclusions
Automated proximal sensing is a powerful complement of gas exchange measurements in ecosystem stations monitoring water and carbon fluxes. It provides spectral signals (F and PRI) encompassing information on light-use in the photosynthetic machinery. Therefore, proximal sensing can improve our understanding of ecosystem function variability in time. However, these signals are affected by additional confounding factors operating at different temporal scales. We demonstrated the capability of SSA to separate the components related to canopy structural/biochemical properties and physiology of these signals from a simulated realistic time series of spectral and physiological variables. This decomposition was especially successful in the case of PRI and  760 Fy , whose relationship with light-use efficiency was strong but still dependent on LAI. Application of the method on measured data showed promising results in the case of  760 Fy . We expect that applying SSA to automated continuous measurements will improve the contribution of proximal sensing to the characterization of ecosystem functional properties.