Using SAR-based products to calculate potato carbon uptake in a tropical Andean region

ABSTRACT

Gross primary productivity (GPP) is an essential parameter to estimate the efficiency of carbon transfer in terrestrial ecosystems. The daily GPP has been monitored in the past using mainly optical satellite imagery. However, GPP has never been monitored using satellite Synthetic Aperture Radar (SAR) imagery. We evaluate the possibility of using Sentinel-1 SAR data to estimate GPP and compare with field GPP measurements using eddy covariance (EC) systems in three commercial potato fields, located in the Andean region of Colombia, using three different water irrigation levels. Raw and processed radar data from Sentinel-1 were compared with the daily and accumulated GPP from the EC. Results indicate that SAR data has a high correlation with the accumulated GPP measured by the EC in the field. The normalized radar backscatter and radar brilliance coefficients with VH polarization show a good correlation with the EC accumulated GPP in irrigated potato fields (R²: 0.77–0.81), while radar vegetation indices show a good correlation with the EC accumulated GPP in potato fields with no irrigation (R²≈0.82). In particular, the accumulated GPP during the whole potato crop cycle was estimated with good accuracy (~5% error with irrigation and ~10% error with no irrigation).

KEYWORDS:

• Synthetic aperture radar
• eddy covariance
• vegetation indices
• gross primary productivity
• solanum tuberosum L
• Sentinel-1

Introduction

Monitoring the total carbon cycle is important to determine if an ecosystem is a sink or a source of greenhouse gas (GHG) emissions, which means that the ecosystem is mitigating GHG emissions (sink) or increasing GHG emissions (source) leading to global warming and climate change. A critical parameter in the carbon cycle is the vegetation’s gross uptake of CO₂ called Gross Primary Productivity (GPP) (Anav et al., 2015). The accumulated GPP is used as a parameter to determine the total amount of carbon dioxide fixed by vegetation through photosynthesis and is the basis for quantifying the efficiency of the ecosystem carbon transfer (C. Liu et al., 2019; Woodwell & Whittaker, 1968). It has been found that the global GPP exhibits spatial and temporal changes due to nutrient availability and climatic conditions (Davi et al., 2006; Woodwell & Whittaker, 1968), requiring continuous monitoring. The GPP can be estimated in the field using the Eddy Covariance (EC) technique (Tagesson et al., 2015). However, EC measurements are limited to the tower footprint area (Tagesson et al., 2015). Optical satellite remote sensing platforms have been used in the past to monitor the daily GPP based on vegetation indexes (VIs) (empirical models) (Badgley et al., 2017; Box et al., 1989; Castaño-Marin et al., 2023; Gitelson et al., 2008; Madugundu et al., 2017; Sims et al., 2006; G. Wu et al., 2020; Xin et al., 2017; Q. Zhang et al., 2015), VIs combined with meteorological variables such as temperature and photosynthetically active radiation (PAR) (semi-empirical models) (Badgley et al., 2017; Gao et al., 2014; Jiang et al., 2021; Sims et al., 2008; C. Wu et al., 2009, 2010; Q. Zhang et al., 2015), and light use efficiency models (LUE) (Badgley et al., 2017; Box et al., 1989; Castaño-Marin et al., 2023; Gitelson et al., 2008; Madugundu et al., 2017; Running et al., 2000; Sims et al., 2006; C. Wu et al., 2010; Xin et al., 2017; Q. Zhang et al., 2015). LUE models estimate the PAR absorbed by the plant canopy (APAR) and the efficiency of converting APAR into fixed carbon. Nevertheless, satellite optical imagery is limited by high cloud cover in humid areas such as the tropics. Cloud-free multispectral imagery using unmanned aerial vehicles (UAVs) has been recently proposed by Castaño-Marin et al. (2023) to monitor the daily GPP in the Andean region.

Satellite-based active sensors such as synthetic aperture radars (SAR) can capture day and night data under all weather conditions and penetrate cloud cover (Mandal et al., 2021; Steele-Dunne et al., 2017). In particular, SAR-based satellite remote sensing has been used in the past to monitor crops in terms of inventory (classification) (Georgiana et al., 2019; Harfenmeister et al., 2019; Kaplan et al., 2021; Khabbazan et al., 2019; C. Liu et al., 2019; McNairn & Brisco, 2004; Shang et al., 2022; Sivasankar et al., 2018) and in terms of crops’ biophysical parameters such as canopy height (Georgiana et al., 2019; Harfenmeister et al., 2019; Kaplan et al., 2021; Khabbazan et al., 2019; C. Liu et al., 2019; McNairn & Brisco, 2004; Shang et al., 2022; Sivasankar et al., 2018), leaf area index (LAI) (Harfenmeister et al., 2019; Kaplan et al., 2021; C. Liu et al., 2019; Mandal et al., 2021; McNairn & Brisco, 2004; Shang et al., 2022; Sivasankar et al., 2018; Steele-Dunne et al., 2017), vegetation water content (VWC) (Harfenmeister et al., 2019; Kushwaha et al., 2022; Mandal et al., 2020; McNairn & Brisco, 2004; Shang et al., 2022; Steele-Dunne et al., 2017; Xie et al., 2019), water stress (McNairn & Brisco, 2004; Steele-Dunne et al., 2017), soil moisture (Georgiana et al., 2019; C. Liu et al., 2019; McNairn & Brisco, 2004; Shang et al., 2022; Steele-Dunne et al., 2017), phenological stage (Khabbazan et al., 2019; C. Liu et al., 2019; Mandal et al., 2020, 2021; Shang et al., 2022), biomass (Georgiana et al., 2019; Harfenmeister et al., 2019; Kushwaha et al., 2022; C. Liu et al., 2019; Mandal et al., 2020, 2021; Shang et al., 2022; Sivasankar et al., 2018; Steele-Dunne et al., 2017; Yang et al., 2019), yield (Barbouchi et al., 2022; Georgiana et al., 2019; C. Liu et al., 2019; McNairn & Brisco, 2004; Shang et al., 2022; Sivasankar et al., 2018), and Kc (Kaplan et al., 2021, 2023); mostly using the C and K bands that interact more with the canopy (Georgiana et al., 2019; C. Liu et al., 2019; Mandal et al., 2021; McNairn & Brisco, 2004; Shang et al., 2022; Sivasankar et al., 2018; Steele-Dunne et al., 2017). However, no studies have been done to monitor GPP with SAR, specifically, carbon uptake. GPP is different from biomass because GPP is an indicator of productivity and the rate of carbon fixation in the ecosystem, while biomass provides a static measure of the amount of carbon accumulated in the form of organic matter (Chi et al., 2016; Fu et al., 2020). GPP allows the assessment of the ecosystem’s impact on the global carbon cycle (Pahari et al., 2018).

In this study, we use preprocessed Sentinel-1 Ground Range Detected (GRD) SAR Level 1 (L1) remote sensing imagery and find correlations with GPP measurements taken with an EC system in the field to estimate carbon sequestration in potato crops under different kinds of water irrigation levels, in a tropical Andean region. We investigate the relationship between Sentinel-1, L1 GRD preprocessed SAR products such as the normalized radar backscatter (σ, with polarizations VV and VH), the radar brightness coefficient (β, with polarizations VV and VH), several vegetation radar indexes proposed in the literature, and the daily and the accumulated GPP. We derived empirical models relating Sentinel-1 GRD preprocessed SAR products with the daily and accumulated GPP. We also investigate the effect of different water irrigation levels on the GPP accumulated during the crop life cycle and its effect on the GRD SAR products for each irrigation level. The region is known for its high potato production levels, which makes it an ideal place to study the impact of agricultural practices on carbon uptake. Potato is now the world’s third most important food crop after wheat and rice in terms of human consumption, and it is a staple food in the Andean highlands (Devaux et al., 2021).

Materials and methods

Test sites and field measurements

This study was conducted using three different water irrigation levels in three commercial potato fields in Colombia’s Andean region, the western savannah of Cundinamarca (Figure 1). The first potato field was in the municipality of Facatativa, with an area of 9.15 hectares (ha) (4.802348° N, 74.28834° W; ∼2568 m above sea level), this plot was planted on 27 February 2020, and harvested on 30 July 2020 (153 days). This plot was irrigated by sprinklers according to a fixed schedule without considering the technical criteria of water balance. The total irrigated water was 100 mm. The second potato field was in the municipality of Tenjo, with an area of 6.0 ha (4.87033° N, 74.1294° W; ∼2572 m above sea level), this plot was planted on 1 August 2020, and harvested on 9 December 2020 (130 days), without irrigation (rainfed). The third plot was in the municipality of Subachoque, on a 3.11 ha area (4.888668 N, 74.18668 W; ∼2609 m above sea level), this plot was planted on 10 February 2022, and harvested on 9 August 2022 (179 days), and was irrigated by sprinklers, according to the crop’s water needs detected through a water balance, to maintain the soil under optimal water conditions near field capacity through the crop cycle and the total irrigated water was 59.6 mm. The Diacol Capiro variety was used in all three production systems, with a planting density of 33,333 plants ha⁻¹.

Figure 1. Study area. Base maps obtained from ESRI® geographic service. Imagery base maps obtained from ©Maxar Technologies on Facatativa (date 11 July 2021), and Subachoque and Tenjo (date 26 January 2021). Up: Geographical location of the three sites considered. Bottom: Images of the three potato plots, plot boundaries are shown in yellow, the EC tower is shown as a red star, and the EC footprint climatology of each site is shown in green.

The periods between June and August and between December and February have the lowest rainfall due to the double passage of the Intertropical Convergence Zone (ITCZ). The accumulated precipitation in each site was as follows: Facatativa 271 mm, Tenjo 229 mm, and Subachoque 681.4 mm. All plots are in a fluvial-lacustrine plain, with flat slopes (0–1%), terraces between low and medium levels, and soils derived from volcanic ash. The Tenjo plot is part of the Helena soil unit and the Subachoque plot is part of the Toscana soil unit. Both belong to the Andisol order (Pachic Melanudands in Tenjo and Typic Hapludands in Subachoque). Otherwise, the Facatativa plot is on an Inceptisol (Andic Humudepts) within the Vuelta soil unit, characterized by deep soils, well-drained, and medium to moderately fine textures (IGAC, 2013). Common practices employed in these production systems include sowing vegetative or asexual seeds, foliar and soil fertilization, weed control, hilling, insecticide and fungicide treatment, and chemical dehaulming or burning. Burning involves the application of a non-selective post-emergent contact herbicide during the senescence phase of plants. The herbicide causes immediate foliage death, allowing the tubers in the soil to develop their epidermis within a period of two to 3 weeks after application. This is done to accelerate the crop’s harvest maturity and harvesting (Martínez-Maldonado et al., 2021).

EC system measurements

The EC technique was used to measure carbon dioxide (CO₂) fluxes between the soil – crop system and the atmosphere. Data quality filters were applied, and gap-filling methods were used to complete the missing data. Accurate net ecosystem exchange (NEE) partition into GPP and ecosystem respiration (RECO) was done by fitting a light response curve (Tagesson et al., 2015). In Facatativa (irrigated production system), the EC station was installed on 20 March 2020, 22 days after sowing (DAS), In Tenjo (rainfed production system), the EC station was installed on 13 August 2020, (12 DAS), in Subachoque (optimal irrigated production system), was installed on 12 February 2022 (2 DAS). The EC tower consisted of an IRGASON system, which comprised an open-path gas analyzer (EC 150, Campbell Scientific, Inc., Logan, Utah, USA), and a 3D sonic anemometer (CSAT3A, Campbell Scientific, Inc., Logan, Utah, USA).

The height and sampling frequency of the sonic anemometer, as well as the description and configuration of sensors for monitoring PAR, global solar radiation (Rg), air and soil humidity (HR, WSC), air and soil temperature (Tair, Ts), soil heat flux (G), and precipitation, are described by Martínez-Maldonado et al. (2021).

The area of influence for surface-atmosphere exchange in an EC system, called the flux footprint or simply the footprint, contains the spatial extent and position of the surface area contributing to the turbulent flux at the EC receptor position, at a specific point in time (Kljun et al., 2004, 2015). Several footprint models have been developed in the past, see for instance (Chen et al., 2009; Kljun et al., 2004, 2015; Kormann & Meixner, 2001; H. Zhang & Wen, 2015) and the references therein. The EC footprint is typically computed every 30 minutes (Chen et al., 2009; Kljun et al., 2015; H. Zhang & Wen, 2015). Here, we used the Easy Flux DL software from Campbell Scientific (Campbell Scientific, 2023) that averages the data fields stored in the data logger from the EC system measurements over a period of 30 minutes (Campbell Scientific, 2022) using the Kljun et al. (2004) footprint model. EC flux measurements can provide estimates of components of the natural carbon cycle in the area by aggregating (Chen et al., 2009) the time series of computed footprints (every 30 minutes) into long-term flux footprint patterns called the footprint climatology of the area, for selected hours, days, months, seasons, or years (Chen et al., 2009; Kljun et al., 2015; H. Zhang & Wen, 2015). The footprint climatology provides a quantitative estimate of the changes in the spatial representativeness of the EC system over time (H. Zhang & Wen, 2015). Combining footprint climatology with remotely sensed data provides spatially explicit information on vegetation structure, topography, and possible source/sink influences on the measured fluxes (Kljun et al., 2015). Here, we derived a footprint climatology for each site (see Figure 1) by accumulating the estimated time-series footprint (every 30 minutes) during each crop life cycle (see the total number of days for each crop cycle, in the previous test sites and field measurements section). We use the footprint climatology R version provided by Kljun et al. (2015) (although, there are also Python and Matlab versions of this algorithm), to aggregate the time series of footprints into a footprint climatology.

SAR images and processing

SAR imagery was obtained from Sentinel-1 (S1), a mission developed by the European Space Agency (ESA) and its Copernicus initiative. S1 encompasses a constellation of two polar-orbiting satellites (S1A and S1B), operating day and night and capturing C-band SAR imagery, with a nominal frequency range from 4 to 8 GHz (center frequency of 5.405 GHz), and 3.75 to 7.5 cm wavelength with an incident angle that ranges from 20° to 50°. S1 allows imagery acquisition regardless of weather and illumination conditions. However, on 23 December 2021, S1B suffered an anomaly related to the instrument electronics power supply provided by the satellite platform, leaving it unable to deliver radar data (European Space Agency, 2022).

S1 can operate a wide swath of 250 km, acquiring data in four modes: Stripmap – SM, Interferometric Wide swath – IW, Extra Wide swath – EW, and Wave – WV; and distributes products in four categories: raw Level-0, processed Level-1 Single Look Complex – SLC, processed Level-1 Ground Range Detected – GRD, and processed Level-2 Ocean – OCN. Also, S1 supports operation in single polarization (HH or VV) and dual polarization (HH+HV or VV+VH). The two-satellite constellation offers a 6-day exact repeat cycle at the equator; nevertheless, a single satellite can map in the IW mode once every 12 days, in a single flight direction (ascending or descending). Level-1 GRD products are available in three-pixel spacings: full in mode SM (3.5 × 3.5 m), high in modes SM (10 × 10 m), IW (10 × 10 m), and EW (25 × 25 m); medium in modes SM (40 × 40 m), IW (40 × 40 m), and EW (40 × 40 m), and mode WV (25 × 25 m). The IW swath mode is the main acquisition mode over land (European Space Agency, 2023b).

This study used high-resolution IW GRD products consisting of SAR imagery that have been detected, multi-looked, and projected to ground range using an Earth ellipsoid model (Filipponi, 2019; Veloso et al., 2017). Forty-five S1 images were acquired from the Copernicus Open Access Hub (European Space Agency, 2023a) (Table 1) covering the three plots of potato crops described in the test sites section, on all the dates available within the period starting with the first date of EC data acquisition and ending with the date before burning the crops or chemical dehaulming (Table 1).

Table 1. Sentinel-1 imagery used in this study.

Site	Period	Mission	Flight direction	Number of images
Facatativa	19-Mar-2020 06-Jul-2020	S1B	Ascending	9
			Descending	9
Tenjo	11-Aug-2020	S1B	Ascending	9
	23-Nov-2020		Descending	8
Subachoque	11-Feb-2022 22-Jul-2022	S1A	Descending	10

S1 IW SAR GRD L1 imagery was preprocessed using the standardized workflow for S1 GRD L1 imagery (Bauer-Marschallinger et al., 2021; Filipponi, 2019; Google Earth Engine, 2023; D. Kumar, 2021; Sentinel Hub, 2022; Wagner et al., 2021) designed to be performed using the Sentinel Applications Platform (SNAP), open-source software with executable tools and Application Programming Interfaces (APIs) which have been developed by ESA to facilitate the utilization, viewing, and processing of a variety of remotely sensed data (European Space Agency, 2023c). We made some small modifications on the standard preprocessing workflow, consisting of six main steps (Figure 2). The first step consists of selecting a subset of the study area to avoid processing the entire SAR image. The second step applies the orbit file, available for each SAR scene in its product metadata providing accurate information about the satellite position and velocity. The third step performs thermal noise removal to reduce noise effects in image intensity in both polarizations (VV and VH). The fourth step is calibration, a procedure that converts digital pixel values to radiometrically calibrated SAR backscatter, in this case, sigma naught (${\sigma }_0$) or normalized radar backscatter and beta naught (${\beta }_0)$ or radar brightness coefficient. The fifth step applies a range Doppler terrain correction for the distortion caused by topography based on a digital elevation model (DEM) implementing an orthorectification method (Small & Schubert, 2008) for geocoding SAR scenes from images in radar geometry. The DEM used was obtained from NASA’s Shuttle Radar Topography Mission (SRTM) with 1 arc-second resolution (Farr et al., 2007). The output of this standardized preprocessing workflow of Sentinel-1 SAR GRD L1 data consists of georeferenced images of σ and β (VV and VH polarizations) with a spatial resolution of 10 m × 10 m i.e. identical to the high-resolution spectral bands in Sentinel-2 (S2) imagery, enabling registration of these two modalities (Filipponi, 2019). We did not perform speckle filtering as in Filipponi (2019) since it averages SAR values using a moving window, reducing resolution and the native values on each pixel as indicated in Kaplan et al. (2023). Because of the strong along-track Doppler centroid variation of IW in S1 (Wegmuller et al., 2016), a co-registration operation is needed to spatially align the processed time-series images using SNAP (the sixth step). With the workflow described before and illustrated in Figure 2, five layers and/or bands were obtained in each SAR scene: ${\sigma }_0$ with polarization VH (${\sigma }_{VH}^0$), ${\sigma }_0$ with polarization VV (${\sigma }_{VV}^0$), ${\beta }_0$ with polarization VH (${\beta }_{VH}^0$), ${\beta }_0$ with polarization VV (${\beta }_{VV}^0$), and the local incidence angle ($\mathit{\theta }$) in decimal degrees.

Figure 2. SAR-imagery (Sentinel-1) Ground Range Detected (GRD) L1 preprocessing workflow. Sigma naught (${\sigma }_0$), beta naught (${\beta }_0)$, ${\sigma }_0$ with polarization VH (${\sigma }_{VH}^0$), ${\sigma }_0$ with polarization VV (${\sigma }_{VV}^0$), ${\beta }_0$ with polarization VH (${\beta }_{VH}^0$), ${\beta }_0$ with polarization VV (${\beta }_{VV}^0$), and the local incidence angle ($\mathit{\theta }$).

Using the incidence angle, a normalization of $\mathit{\sigma }$ and $\mathit{\beta }$ was performed as indicated in Kaplan et al (2021, 2023):.

(1) ${\sigma }^0\mathit{corr}={\sigma }^0\times \mathit{\theta }$(1)

(2) ${\beta }^0\mathit{corr}=\frac{{\beta }^0}{\mathit{sin}\left({\left(\mathit{Radians}\left(90-\mathit{\theta }\right)\right)}^3\right)}$(2)

Hence, new images were added to the dataset: ${\sigma }_{VH}^0\mathit{corr}$, ${\sigma }_{VV}^0\mathit{corr}$, ${\beta }_{VH}^0\mathit{corr}$, and ${\beta }_{VV}^0\mathit{corr}$ using Equations (1)(1) ${\sigma }^0\mathit{corr}={\sigma }^0\times \mathit{\theta }$(1) and (2(2) ${\beta }^0\mathit{corr}=\frac{{\beta }^0}{\mathit{sin}\left({\left(\mathit{Radians}\left(90-\mathit{\theta }\right)\right)}^3\right)}$(2) ). In addition, four radar vegetation indexes were computed: the cross-polarization ratio (Copol) (Mandal et al., 2020), the dual-polarized Ratio Vegetation Index (RVI) (Kaplan et al., 2021; Mandal et al., 2020), the Dual Polarization SAR Vegetation Index (DPSVI) (dos Santos et al., 2021; Mandal et al., 2020; Periasamy, 2018), and the modified DPSVI (DPSVIm) (dos Santos et al., 2021):

(3) $\mathit{Copol}=\frac{{\sigma }_{VH}^0}{{\sigma }_{VV}^0}$(3)

(4) $\mathit{RVI}=\frac{4{\sigma }_{VH}^0}{{\sigma }_{VH}^0+{\sigma }_{VV}^0}$(4)

(5) $\mathit{DPSVI}=\left(\frac{\left(max\left\{{\delta }_{VV}^0\right\}-{\delta }_{VV}^0\right)+{\sigma }_{VH}^0}{\sqrt{2}}\right)\text{break}\left(\frac{{\delta }_{VV}^0+{\delta }_{VH}^0}{{\delta }_{VV}^0}\right){\delta }_{VH}^0$(5)

(6) $\mathit{DPSVIm}=\frac{{\left({\delta }_{VV}^0\right)}^2+{\delta }_{VV}^0{\delta }_{VH}^0}{\sqrt{2}}$(6)

where $\mathit{max}\left\{{\delta }_{VV}^0\right\}$ is the maximum value of ${\delta }_{VV}^0$ in the image. Notice that Mandal et al. (2020) has a wrong definition of the DPSVI, the correct definition of DPSVI and DPSVIm were taken from dos Santos et al. (2021), which is simpler and more intuitive than in the original version of Periasamy (2018).

Model calibration

The integration of the accumulated GPP, measured within the footprint climatology of the EC tower and the preprocessed Sentinel-1 GRD imagery is identical to the integration between optical remote sensing imagery and the GPP measured with the EC, given that preprocessed Sentinel-1 GRD imagery has the same spatial resolution than the high-resolution bands in Sentinel 2. Hence, a representative measure of preprocessed S1 GRD imagery (σ, β, and radar vegetation indices) within the pixels of the EC footprint climatology area (Figure 1) can be taken such as the mean value (Castaño-Marin et al., 2023; Madugundu et al., 2017; Noumonvi et al., 2019; Spinosa et al., 2023). We found that the accumulated GPP measured by the EC on the same dates of S1 GRD image acquisition (Table 1) has a linear relationship with some of the 12-radar products (Figures 3–5), so a linear model between the EC accumulated GPP and radar data was fitted using the Numpy Python package and a Polynomial model of degree one, which uses least-squares minimization. The Normalized Root Mean Square Error (NRMSE) and the coefficient of determination (R²) were computed for each linear model:

(7) $\mathit{RMSE}=\sqrt{\frac{{\sum }_{i=1}^N{\left(x_i-{\hat{x}}_{\mathit{i}}\right)}^2}{\mathit{N}}}$(7)

(8) $\mathit{NRMSE}=\frac{\mathit{RMSE}}{\max \left\{x\right\}\min \left\{x\right\}}$(8)

(9) $R^2=1-\frac{{\sum }_{i=1}^N{\left(x_i-{\hat{x}}_{\mathit{i}}\right)}^2}{{\sum }_{i=1}^N{\left(x_i-\bar{x}\right)}^2}$(9)

Figure 3. Linear regression between the GPP accumulated in Facatativa and (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, (d) ${\beta }_{VH}^0\mathit{corr}$, (e) ${\sigma }_{VV}^0$, (f) ${\sigma }_{VV}^0\mathit{corr}$, (g) ${\beta }_{VV}^0$, (h) ${\beta }_{VV}^0\mathit{corr}$, (i) Copol, (j) RVI, (k) DPSVI, and (l) DPSVIm SAR data.

Figure 4. Linear regression between the GPP accumulated in Subachoque and (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, (d) ${\beta }_{VH}^0\mathit{corr}$, (e) ${\sigma }_{VV}^0$, (f) ${\sigma }_{VV}^0\mathit{corr}$, (g) ${\beta }_{VV}^0$, (h) ${\beta }_{VV}^0\mathit{corr}$, (i) Copol, (j) RVI, (k) DPSVI, and (l) DPSVIm SAR data.

Figure 5. Linear regression between the GPP accumulated in Tenjo and (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, (d) ${\beta }_{VH}^0\mathit{corr}$, (e) ${\sigma }_{VV}^0$, (f) ${\sigma }_{VV}^0\mathit{corr}$, (g) ${\beta }_{VV}^0$, (h) ${\beta }_{VV}^0\mathit{corr}$, (i) Copol, (j) RVI, (k) DPSVI, and (l) DPSVIm SAR data.

where, N is the number of points (see Table 1), $x_i$ the ith observed accumulated GPP value, ${\hat{x}}_{\mathit{i}}$ the accumulated GPP predicted by the model, $\max \left\{\mathit{x}\right\}$ is the maximum value of all observed points, and $\min \left\{\mathit{x}\right\}$ is the minimum of all observed points. Hence, $\max \left\{\mathit{x}\right\}-\min \left\{\mathit{x}\right\}$ corresponds to the range of accumulated GPP. We use NRMSE rather than RMSE since the accumulated GPP is quite different on each crop site (Facatativa, Subachoque, and Tenjo), given the different water crop management on each site (Section Test sites and field measurements) affecting the RMSE so the different RMSE on each site cannot be compared to the others. The NRMSE allows a direct comparation between crop sites.

Leave-one-out cross-validation (LOOC) was used to test the best models found since the number of data points on each site is relatively low and it has been shown that it is an almost unbiased estimator of the true model generalization performance (Cawley & Talbot, 2010). LOOC provides a measure of the variability of the model parameters, in addition to the classical NRMSE and R² statistics. The mean, standard error, and p-values of each model parameter can be computed using t-student statistics. LOOC is performed by eliminating in date order each data point to train the linear model on the remaining points and validate the point eliminated for each data point available, see for instance Castaño-Marin et al. (2023).

Results

The resulting dataset on each acquisition date and site consists of 12 IW S1 GRD L1 preprocessed raster images: 4 from the native SAR products (${\sigma }_{\mathit{VV}},{\sigma }_{\mathit{VH}},{\beta }_{\mathit{VV}},{\beta }_{\mathit{VH}}$), 4 with incidence angle correction (${\sigma }_{VV}^0\mathit{corr},{\sigma }_{VH}^0\mathit{corr},{\beta }_{VV}^0\mathit{corr},{\beta }_{VH}^0\mathit{corr}$), and 4 from radar vegetation indexes (Copol, RVI, DPSVI, DPSVIm). We use the georeferenced footprint climatology polygons (Figure 1) to filter out pixels outside the EC footprint climatology. Slightly better correlation results were found using the median rather than the mean values within the EC footprint on Facatativa and Subachoque. Slightly better results were found on Tenjo using the mean rather than the median values. This might be due to the presence of a small speckle noise on the Facatativa and Subachoque sites that is filtered using the median.

Good GPP-SAR correlations were found using the ascending S1B orbit (Facatativa and Tenjo), while there was no GPP-SAR correlation in the descending S1B orbit. The GPP-SAR correlation on the whole dataset (ascending and descending) is poor, due to the ascending data. Some studies also show a preference for either ascending or descending SAR orbits when monitoring crops (Harfenmeister et al., 2019; Kaplan et al., 2021; Khabbazan et al., 2019; V. Kumar et al., 2022; Kushwaha et al., 2022; Mandal et al., 2020, 2021; Nasirzadehdizaji et al., 2019; Sonobe, 2019; Valcarce-Diñeiro et al., 2019; Xie et al., 2019; Yang et al., 2019). Other SAR studies on crops combine both orbits (Kaplan et al., 2023; S. Liu et al., 2021; Mandal et al., 2021). Harfenmeister et al. (2019) explain that images acquired with the same orbit are more similar to each other, but also descending orbit images are more unstable; probably due to moisture in the morning (see also Kaplan et al. (2021)). Good GPP-SAR correlations were found in the S1A descending orbit (the only data available due to the S1B mission failure in 2021).

A correlation was also found between some of the 12 parameters radar data and the daily GPP, but this relationship is highly non-linear, it quickly saturates, and it is very sensitive to the first few values (before saturation); hence, we do not recommend using this highly nonlinear daily GPP-SAR data relationship (see Appendix), since the accumulated GPP is much more useful than the daily GPP.

Figures 3–5 show the linear regression between the accumulated GPP on each of the radar data (12 parameters) acquired between the sowing date and burn of the potato crop on three sites: Facatativa (Figure 3), Subachoque (Figure 4), and Tenjo (Figure 5). As can be seen from Figures 3 and 4, the native ${\sigma }_{VH}^0\mathit{\hspace{0.17em}}$and ${\beta }_{VH}^0$ have the best NRMSE and R² results correlating with the accumulated GPP, their incidence angle corrected versions ${\sigma }_{VH}^0\mathit{corr}\mathit{\hspace{0.17em}}$and ${\beta }_{VH}^0\mathit{corr}$ also have good NRMSE and R² results correlating with the accumulated GPP, but slightly lower than the native values. The native ${\sigma }_{VV}^0\mathit{\hspace{0.17em}}$and ${\beta }_{VV}^0$ and their incidence angle corrected versions ${\sigma }_{VV}^0\mathit{corr}\mathit{\hspace{0.17em}}$and ${\beta }_{VV}^0\mathit{corr}$ seems to correlate well with the GPP accumulated in the Facatativa site but do not correlate with the GPP accumulated in the Subachoque site. A closer look at ${\sigma }_{VV}^0\mathit{\hspace{0.17em}}$and ${\beta }_{VV}^0$ in Figure 3 shows that the last five points out of nine fall close to each other in terms of their SAR values, while the accumulated GPP varies strongly. In addition, there is a gap in the SAR values between the two groups of points (the first four and the last five). In addition, the correlation between Copol, RVI, DPSVI, DPSVIm and the accumulated GPP is weak. On the other hand, in the Tenjo site, the correlation between ${\sigma }_{VH}^0,$ ${\beta }_{VH}^0$ and their incidence angle corrected versions with the accumulated GPP is weak, also no correlation between ${\sigma }_{VV}^0,$ ${\beta }_{VV}^0$, ${\sigma }_{VV}^0\mathit{corr},\mathit{\hspace{0.17em}}$and ${\beta }_{VV}^0\mathit{corr}$ and the accumulated GPP was found. In contrast, the best NRMSE and R² results correlating with the GPP accumulated by Tenjo correspond to the Copol, RVI, and DPSVI radar vegetation indexes. This could be because both Facatativa and Subachoque have good crop management (irrigation, and other agronomic practices) leading to a high accumulated GPP; hence, their results coincide. Tenjo had poor crop management (no irrigation and pests) leading to a very low accumulated GPP, so it does not have the same behavior as Facatativa and Subachoque. Figure 6 summarizes all the results.

Figure 6. Accumulated GPP from radar data in terms of (a) the NRMSE and (b) the R² score, on the three sites.

Based on the results of Figures 3–6, the best SAR parameters for high GPP accumulated by potato crops (good and excellent crop management) were: ${\sigma }_{VH}^0,$ ${\beta }_{VH}^0$, ${\sigma }_{VH}^0\mathit{corr},$ ${\beta }_{VH}^0\mathit{corr}$, while the best parameters for low accumulated GPP (no crop management) were as follows: Copol and RVI radar vegetation indexes. Note that DPSVI had a not so good overall performance on the three test sites and the DPSVIm had a very good performance in Facatativa only, but worse than the DPSVI in the two other sites. This may be because the DPSVI is the product of two indexes and ${\sigma }_{VH}^0$ (see Equation (6)(6) $\mathit{DPSVIm}=\frac{{\left({\delta }_{VV}^0\right)}^2+{\delta }_{VV}^0{\delta }_{VH}^0}{\sqrt{2}}$(6) ), we know that ${\delta }_{VH}^0$ correlates well with the accumulated GPP in Facatativa and Subachoque

but not with Tenjo and $\left(\frac{{\delta }_{VV}^0+{\delta }_{VH}^0}{{\delta }_{VV}^0}\right)=1+\mathit{Copol}$ correlates well only on Tenjo, so their product is worse than each one of them separately on the three test sites. On the other hand, DPSVIm (see Equation (6)(6) $\mathit{DPSVIm}=\frac{{\left({\delta }_{VV}^0\right)}^2+{\delta }_{VV}^0{\delta }_{VH}^0}{\sqrt{2}}$(6) ), relies heavily on ${\delta }_{VV}^0$, which only correlates well with the accumulated GPP (Figure 3) in Facatativa. Using LOOC, the estimated RMSE and R² statistics on unseen data, and the averaged linear model parameters, the estimated standard error of those parameters and their corresponding p-values (t-test) can be computed. Table 2 summarizes the results of LOOC on the best linear models found. In addition, Facatativa models can be used to predict the accumulated GPP in Subachoque with good accuracy, see Table 3. The reason is Facatativa SAR parameters ranges are lower than those in Subachoque, accounting for the larger GPP accumulated in Subachoque (see the ranges of the accumulated GPP in Figures 3 and 4). Facatativa models cannot predict the GPP accumulated in Tenjo since GPP in Facatativa is much higher than in Tenjo. The Subachoque model cannot be used to predict Facatativa accumulated GPP, since the Subachoque model predicts a larger GPP accumulated.

Table 2. Cross-validated best linear models.

Parameter	Model	NRMSE	R^2 score	p-values
FACATATIVA
${\sigma }_{VH}^0$	GPP= −380(69) + 68727(3475)$\times {\sigma }_{VH}^0$	0.179	0.72	4.5×10^−4	1.1×10^−7
${\sigma }_{VH}^0\mathit{corr}$	GPP= −366(84) + 1746(84)$\times {\sigma }_{VH}^0\mathit{corr}$	0.183	0.71	4.2×10^−4	7.7×10^−7
${\beta }_{VH}^0$	GPP= −380(69) +44378(2244)$\times {\beta }_{VH}^0$	0.179	0.72	4.5×10^−4	1.1×10^−7
${\beta }_{VH}^0\mathit{corr}$	GPP= −366(84) + 1746(84)$\times {\beta }_{VH}^0\mathit{corr}$	0.193	0.67	2.8×10^−4	4.6×10^−8
SUBACHOQUE
${\sigma }_{VH}^0$	GPP= −1185(196) + 87151(5483)$\times {\sigma }_{VH}^0$	0.174	0.71	1.5×10^−4	1.2×10^−7
${\sigma }_{VH}^0\mathit{corr}$	GPP= −1124(164) + 2282(126)$\times {\sigma }_{VH}^0\mathit{corr}$	0.172	0.71	6.5×10^−5	4.4×10^−8
${\beta }_{VH}^0$	GPP= −1185(195) +53887(3390)$\times {\beta }_{VH}^0$	0.174	0.71	1.5×10^−4	1.2×10^−7
${\beta }_{VH}^0\mathit{corr}$	GPP= −1117(159) + 36738(1933)$\times {\beta }_{VH}^0\mathit{corr}$	0.171	0.72	5.4×10^−5	3.0×10^−8
TENJO
Copol	GPP= −308(29) + 4286(241) × Copol	0.177	0.73	7.5×10^−6	2.2×10^−7
RVI	GPP= −403(39) + 1488(88) × RVI	0.177	0.72	8.7×10^−6	3.2×10^−7

Table 3. Validation of Facatativa best linear models on Subachoque.

Parameter	Model	NRMSE	R^2 score
${\sigma }_{VH}^0$	FACATATIVA	0.164	0.74
${\sigma }_{VH}^0\mathit{corr}$		0.160	0.75
${\beta }_{VH}^0$		0.173	0.71
${\beta }_{VH}^0\mathit{corr}$		0.161	0.75

Figure 7 shows the accumulated GPP measured with the EC (blue) on the same dates the S1 images were acquired and the corresponding GPP accumulated in Facatativa using the model for ${\sigma }_{VH}^0$ (Table 2) as a function of the Days After Sowing (DAS). It is important to extrapolate the GPP accumulated on the day before chemically dehaulming the crop (130 DAS) to estimate the total GPP accumulated during the crop cycle. Hence, we perform a simple fitting of the accumulated GPP as a function of DAS to be able to interpolate the accumulated GPP on any DAS, based on the shape of the accumulated GPP measured by the EC.

Figure 7. GPP accumulated in Facatativa measured (EC) vs estimated using ${\sigma }_{VH}^0$ (Table 2) as a function of DAS.

As can be seen in Figure 7, the GPP measured with the EC is slightly non-linear on the first DAS, then becomes linear. We propose the following curve, to interpolate the accumulated GPP as a function of the DAS corresponding to the dates the images were taken:

(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10)

(11) $\mathit{X}=5\frac{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}}{\mathit{max}\left\{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}\right\}}$(11)

Hence, the DAS is normalized to be in the 0.0 to 5.0 range and the regression consists only of finding $\mathit{C}$, the slope of the linear part of the curve that fits the accumulated GPP estimated with the best cross-validated linear models (Table 2) on the dates available. We use the curve fit function in the scipy. optimize Python module to find $\mathit{C}$. Notice that $\mathit{GPP}$ in Equation 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) starts at 0 when $\mathit{X}=0$ and continuously increases (as the accumulated GPP should). Leave one out cross-validation can be used also here to estimate the NRMSE and R² statistic on unseen data. Figure 8 shows the NRMSE and R² statistics of the LOOC using Equations 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) and 11(11) $\mathit{X}=5\frac{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}}{\mathit{max}\left\{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}\right\}}$(11) to fit the accumulated GPP as a function of the DAS.

Figure 8. Cross-validation results of the accumulated GPP fitted to DAS in terms of (a) the NRMSE and (b) the R² score, between the estimated accumulated GPP using the best linear models and the EC accumulated GPP.

We report here LOOC results using the interpolating curve proposed in Equations 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) and 11(11) $\mathit{X}=5\frac{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}}{\mathit{max}\left\{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}\right\}}$(11) to estimate the curve fit on unseen data. However, in practice this interpolating curve fitting function will be done using all the data points available to estimate $\mathit{C}$. Then the daily accumulated GPP can be found from Equations 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) and 11(11) $\mathit{X}=5\frac{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}}{\mathit{max}\left\{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}\right\}}$(11) by using the continuous daily EC DAS, starting with the DAS of the first image and ending with the DAS of the day before burning the canopy for harvesting. Note that the curve cannot consider dates before the first image acquisition (negative $\mathit{X}$), but it can consider any DAS posterior to that date. Figures 9–11 compare the estimated daily accumulated GPP (Equations 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) and 11(11) $\mathit{X}=5\frac{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}}{\mathit{max}\left\{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}\right\}}$(11) ) with the daily accumulated GPP measured from the EC, in terms of NRMSE and R² statistics for the three sites: Facatativa, Subachoque, and Tenjo. The initial estimate of the accumulated GPP using the best linear models (Table 2) on the dates of image acquisition is also shown as a reference. As can be seen from these figures, the fitting of the initial estimate of the accumulated GPP (Table 2) to daily DAS using a single parameter provides a very good approximation of the daily accumulated GPP measured with the EC.

Figure 9. Facatativa: accumulated GPP fitted to daily DAS (red) compared to the daily accumulated GPP from the EC (green) and the initial estimate of the accumulated GPP (blue) from the best linear models (Table 2) using (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, and (d) ${\beta }_{VH}^0\mathit{corr}$ SAR data.

Figure 10. Subachoque: accumulated GPP fitted to daily DAS (red) compared to the daily accumulated GPP from the EC (green) and the initial estimate of the accumulated GPP (blue) from the best linear models (Table 2) using (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, and (d) ${\beta }_{VH}^0\mathit{corr}$ SAR data.

Figure 11. Tenjo: accumulated GPP fitted to daily DAS (red) compared to the daily accumulated GPP from the EC (green) and the initial estimate of the accumulated GPP (blue) from the best linear models (Table 2) using (a) Copol and (b) RVI.

Using Equations 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) and 11(11) $\mathit{X}=5\frac{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}}{\mathit{max}\left\{\mathit{DAS}-\mathit{min}\left\{\mathit{DAS}\right\}\right\}}$(11) , the total accumulated GPP can be extrapolated to the last DAS date before burning the canopy for harvesting the potatoes. Figure 12 shows the percentage of relative error between the accumulated GPP using the best SAR parameters and curve fitting to DAS compared to the measured accumulated GPP on the date before burning the canopy on the three sites. The percentage of relative error is computed as,

(12) ${\epsilon }_r\left(\mathrm{\%}\right)=\frac{\left|\mathit{GP}{\mathit{P}}_{\mathit{EC}}-\widehat{\mathit{GPP}}\right|}{\mathit{GP}{P}_{\mathit{EC}}}\times 100$(12)

where, $\mathit{GP}{P}_{\mathit{EC}}$ is the total GPP accumulated at the last day before burning the crop measured with the EC and $\widehat{\mathit{GPP}}$ is the total GPP accumulated at the same day using the best linear models based on SAR parameters (Table 2) and fitting to DAS (Equations 10(10) $\mathit{GPP}=\mathit{C}\left(\mathit{X}+\frac{1}{\left(1+\mathit{X}\right)}-1\right)$(10) -12(12) ${\epsilon }_r\left(\mathrm{\%}\right)=\frac{\left|\mathit{GP}{\mathit{P}}_{\mathit{EC}}-\widehat{\mathit{GPP}}\right|}{\mathit{GP}{P}_{\mathit{EC}}}\times 100$(12) ). In the case of Facatativa, the GPP accumulated during the crop cycle, measured with the EC is 2778.84 g m⁻² at 130 DAS and using the two best SAR parameters and DAS fitting is 2638.51 g m⁻² using ${\sigma }_{VH}^0$ and 2638.47 g m⁻² using ${\beta }_{VH}^0$ that is a relative error of 5%. In the case of Subachoque, the GPP accumulated during the crop cycle measured with the EC is 3779.0 g m⁻² at 162 DAS, and using the two best SAR parameters and DAS fitting is 844.47 g m⁻² using ${\sigma }_{VH}^0$ and 3844.07 g m⁻² using ${\beta }_{VH}^0$ that is a relative error of 1.7%. In the case of Tenjo, the GPP accumulated during the crop cycle measured with the EC is 879.06 g m⁻² at 114 DAS and using SAR parameters and DAS fitting is 787.26 g m⁻² using Copol and 787.54 g m⁻² using RVI for a relative error of 10.4%.

Figure 12. Percentage of relative error in the total GPP accumulated during the crop cycle using best SAR parameters and DAS fitting compared to EC measurements.

Discussion

Results show a simple linear relationship between ${\sigma }_{VH}^0$ and ${\beta }_{VH}^0$ radar raw products and the GPP accumulated in potato crops (var. Diacol Capiro) with agronomic management ranging from good to optimal, consisting of water irrigation, pest, and disease control, among other agronomic practices. Radar vegetation indexes (Copol, RVI, DPSVI, DPSVIm) do not seem to have a good correlation with the accumulated GPP in this case. On the other hand, ${\sigma }_{VH}^0$ and ${\beta }_{VH}^0$ do not have a sufficiently good correlation with the GPP accumulated in potato crops with no irrigation. This could be attributed to the water stress experienced by the non-irrigated plants in Tenjo, which caused early senescence, increased pest, and disease damage, and consequently, a reduction in the crop’s leaf area GPP compared to the irrigated crops (Rambal et al., 2014). However, some radar vegetation indexes (Copol and RVI) show a good linear correlation with the accumulated GPP in the case of no water irrigation of the crop. Incidence angle correction (Equations 1(1) ${\sigma }^0\mathit{corr}={\sigma }^0\times \mathit{\theta }$(1) and 2(2) ${\beta }^0\mathit{corr}=\frac{{\beta }^0}{\mathit{sin}\left({\left(\mathit{Radians}\left(90-\mathit{\theta }\right)\right)}^3\right)}$(2) ) does not seem to improve the relationship between radar products and the accumulated GPP (Figures 3–12), in general. This might be because the EC footprint area is small (on average 42 pixels of size 10 m × 10 m) compared for instance to Spinosa et al. (2023) (225 pixels of size 30 m × 30 m); hence, the incidence angles vary very little. Further studies should be performed on areas with a larger EC footprint climatology area, where the radar incidence angle has a larger variation. We found no correlation between the accumulated GPP and radar products with VV polarization, except through radar vegetation indexes. Our results are better than some previous results estimating biophysical parameters from radar imagery (Harfenmeister et al., 2019; Kushwaha et al., 2022; Mandal et al., 2020) and comparable to dos Santos et al. (2021), Kaplan et al. (2021, 2023). and Mandal et al. (2021). Previous studies monitoring the GPP with remote sensing imagery are entirely based on optical sensors, using empirical, semi-empirical, and LUE models. These studies are limited by cloud cover, especially in humid regions such as the tropics, and their estimations are in terms of the daily GPP, requiring daily estimates over the crop cycle to obtain the accumulated GPP, indirectly (Castaño-Marin et al., 2023; Madugundu et al., 2017). Our results are comparable in accuracy to daily GPP estimates using VIs.

The relative error of the GPP accumulated during the crop’s life cycle seems to deteriorate (Figure 12) as the water availability decreases from optimal water irrigation (1.7% in Subachoque) to non-technical water irrigation (5% in Facatativa), and finally rainfed (10.4% in Tenjo). Hence, it seems the relative error in the GPP accumulated during the crop cycle increases as the accumulated GPP reduces, which also affects plant structure due to reduced carbon fixation and early senescence. The accumulated GPP with optimal water irrigation (Subachoque) was 3779.0 g m⁻², under non-technical water irrigation (Facatativa) was 2778.84 g m⁻², and rainfed irrigation only (Tenjo) was 879.06 g m⁻² indicating lower carbon sequestration reduces as water availability reduces to the point of generating water stress and diseases, as indicated before in the case of Tenjo. At this point, considering each plot’s crop life cycle is essential. While Subachoque had a life cycle of 179 days (~6 months), Facatativa and Tenjo had a shorter life cycle (153 and 130 days, respectively). These results indicate that in Tenjo and Facatativa the senescence of the crop occurred much earlier than in Subachoque; hence, the potato crop in Subachoque had a longer life cycle and accumulated a larger amount of GPP (~1000 g m⁻² more than in Facatativa and ~2900 g m⁻² more than in Tenjo). Further work should be made in the future to validate and generalize models for the whole potato production region in Colombia as well as to validate the hypothesis on how the level of water irrigation affects plant structure, SAR measurements, and the error in the estimated accumulated GPP during a full crop cycle. In addition, further studies should be conducted to derive relationships between SAR imagery products and GPP on other kinds of crops (wheat, soybeans, rice, corn, etc.) and other types of vegetation (forests, grasslands, shrublands, savannas, etc.) worldwide by the scientific community. We explored here Sentinel-1 IW GRD products and their relationships with GPP; however, other SAR products such as Single Look Complex (SLC), see for instance Mandal et al. (2020), and their relationship with GPP could also be explored in the future as well as using radar vegetation models (Mandal et al., 2021), such as the semiempirical Water Cloud Model (WCM) that can be paired with machine learning (ML) algorithms for inversion of crop biophysical parameters. These semiempirical models might help reduce the estimation error in the case of very low water availability, by training specialized ML models for this kind of situation. Finally, it has been found that soil water depletion (drought) significantly affects the relationship between GPP and VIs, degrading their prediction accuracy (Huang et al., 2019; Maleki et al., 2022). Hence, low water availability in crops could be an issue for both optical and radar remote sensing of the GPP. In the case of optical remote sensing, this issue can be alleviated by including the soil moisture information within the GPP models (Maleki et al., 2022), this could be a research area in the future for both optical and radar remote sensing of the GPP. In addition, further research work should be done to study how the accumulated GPP changes on different phenological stages of the crops (or on each image acquisition date) and how it relates to the increase in aerial and root (tubers such as potatoes, cassava, carrots, beets, etc.) biomass (carbon fixation).

Currently, radar remote sensing of crops focuses on crop classification and estimation of biophysical parameters (height, LAI, VWC, soil moisture, water stress, biomass, yield, Kc), but not the GPP. The accumulated GPP is a fundamental parameter for monitoring carbon exchange in agroecosystems and plant ecosystems in general. It provides crucial information about the amount of carbon absorbed through photosynthesis and converted into biomass, as well as insights into their productivity and functioning. Additionally, it helps determine their capacity to act as carbon sinks or sources. This information is essential for understanding and addressing the challenges related to climate change and the environmental and productive sustainability of crops, such as potatoes, which are one of the most important crops worldwide (CIP, 2023) and a staple food in Colombia.

Conclusions

This work shows that it is possible to estimate the accumulated GPP of different agro-ecosystems such as potato crops, using cloud-free SAR data from remote sensing platforms and simple linear regression models, with good accuracy compared to EC measurements. The results indicate that the level of water supply applied to potato crops significantly affects radar measurements as well as the accumulated GPP of the crop, due to changes in the crop canopy structure that are visible in the field. Hence, proper water irrigation is needed to improve crop carbon uptake. Cloud-free remotely sensed SAR data and field EC systems could be used in the future to derive robust models to monitor carbon uptake in potatoes, other agroecosystems, and other kinds of vegetation such as forests, grasslands, shrublands, savannas, etc. Since the accumulated GPP provides a measure of the total vegetation carbon uptake, while the current state of the art can only provide daily GPP estimates using optical remote sensing platforms, it is possible that radar remote sensing could be used in the future to monitor the accumulated GPP even in cloud-free areas, as well as the possibility of combining both optical and radar data modalities with machine learning algorithms. Monitoring the accumulated GPP is of prime importance to ensure the efficient use of resources (water, carbon, and nutrients) and to understand the carbon balance of the crop, which can indicate whether the crop is a sink or a source of greenhouse gas emissions. This can contribute to identifying inappropriate farming practices and taking corrective measures aimed at increasing the crop’s productivity potential and achieving greater agro-environmental efficiency in support of climate-smart agriculture.

Acknowledgments

This work is part of a larger project in Corporación Colombiana de Investigación Agropecuaria (AGROSAVIA) named Sistema de Información Agroclimática del cultivo de la papa en la región de Cundinamarca, Colombia (SIAP). We thank, Fabio Ernesto Martínez Maldonado for his contribution to data analysis, Zahara Lasso, Douglas Gómez, Jose Alfredo Molina Varón, Pablo Edgar Jiménez, Óscar Dubán Ocampo Páez, and Jhon Alexander Martínez Morales, for their contribution in the equipment installation process, and farmers Santiago Forero, Wilson Forero, and Alejandro Forero for providing a suitable lot for the study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data presented in this study are available upon request from the Agrosavia Intellectual Property Department. The data are not publicly available due to Agrosavia’s copyright. Your request can be sent to [email protected].

Additional information

Funding

This research was funded by the Fondo de Ciencia, Tecnología e Innovación del Sistema General de Regalías, administered by the Fondo Nacional de Financiación para Ciencia, Tecnología e Innovación—Francisco José de Caldas, Programa Colombia BIO, Gobernación de Cundinamarca and Ministerio de Ciencia, Tecnología e Innovación [MINCIENCIAS] funding number [66153], and Corporación Colombiana de Investigación Agropecuaria [AGROSAVIA] funding number TV19 1000911.

Appendix

The relationship between S1 measurements and the daily GPP is highly nonlinear. Best results were found using a hyperbolic tangent (tanh) activation function between all the raw data: ${\sigma }_{VH}^0,$ ${\beta }_{VH}^0$, ${\sigma }_{VH}^0\mathit{corr},$ ${\beta }_{VH}^0\mathit{corr}$, ${\sigma }_{VV}^0$, ${\beta }_{VH}^0$, ${\sigma }_{VV}^0\mathit{corr}$, ${\beta }_{VV}^0\mathit{corr}$ and the daily GPP:

(A1) $\mathit{GPP}=\mathit{Ctanh}\left(\mathit{ax}\right)$(A1)

where $\mathit{C}$ and $\mathit{a}$ are model parameters to be fit and $\mathit{x}$ is one of the SAR measurements. There is one exception for the case of Subachoque and the ${\sigma }_{VV}^0$, ${\beta }_{VH}^0$, ${\sigma }_{VV}^0\mathit{corr}$, ${\beta }_{VV}^0\mathit{corr}$ measurements, where a linear relationship is better. Best results for the radar vegetation indexes (Copol, RVI, DPSVI, DPSVIm) were found using a nonlinear function based on the lognormal distribution with zero mean and standard deviation $\mathit{s}$. The lognormal probability distribution function in given by:

(A2) $\mathit{pdf}\left(\mathit{x}\right)=\frac{1}{\mathit{sx}\sqrt{2\mathit{\pi }}}{e}^{-\frac{l{n}^2\left(\mathit{x}\right)}{2s^2}}$(A2)

However, this curve first goes up due to the exponential term and then goes down as increases. We do not want that behavior, so the following nonlinear function is used to model the GPP as a function of radar vegetation indexes:

(A3) $\mathit{GPP}=\mathit{C}{e}^{-\frac{l{n}^2\left(\frac{x}{\mathit{max}\left\{\mathit{x}\right\}}\right)}{2s^2\frac{x}{\mathit{max}\left\{\mathit{x}\right\}}}}$(A3)

where $\mathit{C}$ and $\mathit{s}$ are model parameters to be fit and $\mathit{x}$ corresponds to the radar vegetation indexes. Notice that $\mathit{x}$ is normalized by dividing the index by its maximum value so that it ranges in the (0,1] range so that the curve mostly goes up. The $2s^2\mathit{x}$ term in the denominator of Equation A3(A3) $\mathit{GPP}=\mathit{C}{e}^{-\frac{l{n}^2\left(\frac{x}{\mathit{max}\left\{\mathit{x}\right\}}\right)}{2s^2\frac{x}{\mathit{max}\left\{\mathit{x}\right\}}}}$(A3) guarantees that the curve does not go down, but it then saturates.

Figures A1-A3 show the NRMSE and R² statistics of the nonlinear regression curves (Equations A1(A1) $\mathit{GPP}=\mathit{Ctanh}\left(\mathit{ax}\right)$(A1) and A3) fitted to the relationship between SAR data and the daily GPP on the three sites: Facatativa (Figure A1), Subachoque (Figure A2), and Tenjo (Figure A3).

Figure A1. Nonlinear regression between the daily GPP and (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, (d) ${\beta }_{VH}^0\mathit{corr}$, (e) ${\sigma }_{VV}^0$, (f) ${\sigma }_{VV}^0\mathit{corr}$, (g) ${\beta }_{VV}^0$, (h) ${\beta }_{VV}^0\mathit{corr}$, (i) Copol, (j) RVI, (k) DPSVI, and (l) DPSVIm SAR data on Facatativa.

Figure A2. Nonlinear regression between the daily GPP and (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, (d) ${\beta }_{VH}^0\mathit{corr}$, (e) ${\sigma }_{VV}^0$, (f) ${\sigma }_{VV}^0\mathit{corr}$, (g) ${\beta }_{VV}^0$, (h) ${\beta }_{VV}^0\mathit{corr}$, (i) Copol, (j) RVI, (k) DPSVI, and (l) DPSVIm SAR data on Subachoque.

Figure A3. Nonlinear regression between the daily GPP and (a) ${\sigma }_{VH}^0$, (b) ${\sigma }_{VH}^0\mathit{corr}$, (c) ${\beta }_{VH}^0$, (d) ${\beta }_{VH}^0\mathit{corr}$, (e) ${\sigma }_{VV}^0$, (f) ${\sigma }_{VV}^0\mathit{corr}$, (g) ${\beta }_{VV}^0$, (h) ${\beta }_{VV}^0\mathit{corr}$, (i) Copol, (j) RVI, (k) DPSVI, and (l) DPSVIm SAR data on Tenjo.

Notice that despite some good correlations on Facatativa (Figure A1), the data points are mainly located on the saturation region of the nonlinear equations and the first few points determine the ascending part of the curve, hence the curve is very sensitive to those few points. In the case of Subachoque (Figure A2), no good correlations were found between the daily GPP and SAR parameters. In the case of Tenjo (Figure A3), no good correlations were found between the daily GPP and SAR data, except for the radar Copol and RVI vegetation indexes, and even in this case the points are not uniformly distributed, they form patches of points. From these results and the fact that the accumulated GPP is a more desirable measure of GPP than the daily GPP, we do not recommend deriving the daily GPP from SAR data.