Estimating Grassland Production in Central Europe Using Multi-Source Remote Sensing Data and a Novel Compilation of Field Observations

1. Introduction

Grasslands have large economic and societal importance given their role in ecosystem services provisioning including food production, water regulation, soil stabilization and biodiversity [1,2,3]. The role of grasslands in human well-being is also the consequence of their vast spatial extent. Depending on the definition of grasslands they cover about 1/3 of the global land area, or even more [4,5]. Besides the widely known benefits, grasslands are also important components of the global carbon cycle [6]. According to the current knowledge, grasslands can be more effective “carbon pumps” than other biomes in the sense that the accumulation of organic matter in soils is more effective compared to other vegetation types [7,8,9,10]. It was reported that grasslands can also be more efficient regarding the protection of already sequestered carbon compared to e.g., forests [7]. Grasslands support the buildup and stabilization of soil organic carbon through sustained root derived inputs, microbial necromass formation and the development of mineral associated organic matter (MAOM; [11,12]). Given the role of herbaceous vegetation in the global carbon cycle and in animal husbandry, understanding the regional patterns of grassland productivity, their interannual variability and their stress resilience is of high importance [13]. Along with point or farm scale studies of grassland functioning, large scale, spatially explicit investigations are also highly needed [14,15].

Quantification of grassland productivity in large spatial scales (even in the form of aboveground biomass (AGB), harvested aboveground biomass (HAB) or total biomass) is relatively rare in the literature. This is a major scientific gap that is in large contrast with e.g., the cropland related datasets and yield monitoring. Quantification of large-scale biomass or net primary production (NPP) is hampered by lack of observation data (including infrequent and unevenly distributed field campaigns, and persistent under sampling of belowground biomass), and also due to the highly contrasting spatial scales of the observations and the target areas (e.g., compare biomass sampling for ~1 m² versus the area of countries on the order of 100 000 km² or so [16]).

The definition of grasslands, and even the estimation of the spatial extent of grasslands in a given region is also challenging due to lack of consensus on the terminology. For example, the European CORINE Land Cover system differentiates Pastures from Natural grasslands, whereas agronomic definitions often rely on management history (e.g., permanent vs. temporary grassland). In addition, shrub-grass mosaics, fallows and vegetated buffer strips may or may not be classified as grassland depending on the context, resulting in substantial discrepancies in reported grassland area [17]. Clearly, any attempt that tries to upscale in situ point or other estimates based on proxy data of grassland biomass is of high relevance [15].

One obvious approach to estimate grassland areas and biomass production for large spatial scales is the exploitation of remote sensing information. State-of-the-art sensors deployed at polar orbiting satellites provide a variety of data on the land surface including widely used, reflectance based spectral indices like normalized difference vegetation index (NDVI), enhanced vegetation index (EVI) and plant phenology index (PPI) and also derived (i.e., partly observation, partly model based) products like leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR) and other biophysical parameters [18,19]. Clearly, recent technological advances largely support the identification of grassland areas and the indirect estimation of the production and phenological phases of grasses. These datasets could potentially support productivity estimations and upscaling once observation data with sufficient amount and quality is available.

The present study focuses on Hungary, that is located in the drought-prone Central Europe, surrounded by the Alps and the Carpathians [20]. Climate change signal is very strong here during the past two decades with strong overall warming and shift in the precipitation patterns (less rainy periods with larger rain events per episode). This kind of destabilization of the climate challenges all ecosystems including the sensitive grasslands [21]. Severe drought events in the past few years (2022, 2024, 2025) emphasize the need for the quantification of the resilience of grasslands as they are more likely exposed to water shortage due to the shallow rooting zone (as compared e.g., to forests). In this sense, understanding grassland production and its interannual variability is a key factor. As grasslands can be considered as “local climate sentinels”, their evaluation can provide information for other regions of Europe considering climate resilience.

Hungarian grasslands represent the high and low end of grassland productivity along the pedo-climatic gradient ranging from highly productive closed grasslands in nutrient rich forest soils to very low productivity, open sandy grasslands with large bare soil coverage and <1 LAI [22,23]. This property makes the target area an exciting testbed for biomass estimation and investigation of upscaling of point or parcel level observations and field data to the country scale.

The primary aim of this study was to construct a unique, observation-based reference dataset for Hungarian grassland production. We introduce a novel identification method to resolve current ambiguities in grassland spatial coverage and use multi-source Earth Observation data to upscale point level AGB and NPP estimates. Although country-specific, this work provides a transferable framework and an ecological benchmark for assessing resilience of European grasslands, particularly in those regions facing similar pedo-climatic transitions and increasing drought frequency. The upscaling model (point-to-satellite) can be applied to other EU member states that are facing the same “scale mismatch” between ~1 m² plot samples and national statistics.

The study contributes to the assessment of grassland related ecosystem services, and provides a reference for further studies under a transient climate.

2. Materials and Methods

2.1. Target Area

The study focuses on Hungary (with a total area of ~93,030 km²) in Central Europe (Figure 1a). The climate of Hungary is affected by continental, oceanic and mediterranean air masses with long term mean temperature of 11.1 °C (1991-2020), and mean annual precipitation of 616 mm [24]; Figure 1b). In each year of the study period (2017-2024), due to the overall warming, annual average temperatures were higher than the long term mean and annual precipitation sums were typically lower than 616 mm (exception is 2017 and 2019; see Table S1 in the Supplementary material). 2022 and 2024 was characterized by extreme drought, especially in the eastern part of the country.

The land use in Hungary is overwhelmingly dominated by agricultural activity. According to the CORINE database (see Section 2.4.1.) arable lands occupy approximately 50% (over 4.7 million hectares) of the total country area. The second most significant category is forestry; forests cover about 19% of the territory. Natural grasslands and pastures are also prominent, covering roughly 10% of the land. Note that this study examined several different grassland databases in order to compare them in terms of grassland area (see Section 2.4.), which means that this number is not exact but a widely used estimation. Artificial surfaces cover approximately 5%, while specialized cultivation, such as vineyards and fruit plantations represent a smaller but distinct portion of the landscape. Finally, water surfaces make up roughly 2% of the total area.

2.2. Reference Datasets

First, we performed an exhaustive literature review and research (including personal communication with experimentalists) to identify possible data sources in Hungary from a variety of research teams. The target datasets included direct biomass samplings with destructive methods, and we also used forage yield data from Hungarian National Parks along with orthophoto-based bale number estimation and appropriate conversion factors. The collected grassland related production data is mapped in Figure 1c and detailed in Table 1.

2.2.1. Harvested Aboveground Biomass Reference Data

We have data on harvested biomass from three areas within the Kiskunság National Park (Figure 1c with pink, cyan and purple). In Bugac (46.69° N, 19.6° E, data from [25,26]), the surface coverage (%) of each plant species was recorded three times a year within 78 quadrats, 63 arranged as 10-m resolution grid and 15 additional random positions. In 2019 and 2020, measurements were also carried out in a 30-m grid, which means 78 more measurement points. After that, the above-ground biomass removed from the plots from a 0.1 m diameter area, then it was oven-dried for 48 h and weighed.

Although the database, which was collected in a small area, with a dense grid and high temporal resolution, has great research potential, the comparison with satellite images taken at 10-300 meters is better supported by the average biomass value representative of the area. Thus, the mean value of the first samples of each year is examined in reference to a 16-hectare square (Figure 2).

Near Fülöpháza (46.87° N, 19.42° E, data from [27]) within the study area (approximately 1 km²) 16 sites were selected in 1999, eight dominated by Festuca vaginata and another eight by Stipa borysthenica, the two dominant species of sandy grasslands in the region. The percentage cover and density changes of Festuca and Stipa was visually estimated in late June or early July each year (following peak growth but before potential drought-induced dieback). Above-ground biomass was estimated from the surface areal cover (%), based on species specific conversion factors determined through linear regression between visually estimated cover and measured biomass in a separate study in the same site [28]. Unfortunately, the sites are only characterized by one coordinate, so we use the average of the 16 sites as a center point of a 16-hectare square around the point (Figure 2).

Near Orgovány (46.78° N, 19.47° E, data from [29]) homogeneous grassland patches (ca. 5 m in diameter) were chosen that covered the variation in grassland productivity within a 1 km² area. Then one, randomly located 0.5×0.5 m plot in each patch was sampled. Above-ground vascular plant biomass was harvested in each sampling plot at the peak of the growing season. Biomass was sorted by species, and only living material was considered in the analyses. Biomass samples were dried at 60 °C for 48 h and then weighed. Although the aim of the survey was to assess grasslands with different productivity levels, from the perspective of satellite data, they are too close to each other and thus represent a larger, inhomogeneous area covered with trees and bushes. In order to obtain a value characteristic of the landscape, we defined two 16-hectare squares around two transects of samples and averaged the measured biomass values within them (Figure 2). Note that the presence of bushes and trees within squares can be a significant source of error.

2.2.2. Yield Data from Management Registers

We used mowing data originally provided in bales ha^-1 units from four national parks and from the plot around the eddy covariance tower at the Hegyhátsál experimental site [30,31]. Parcels from Bükk and Hortobágy National Park are in Eastern Hungary, near Lake Tisza (Figure 1c with magenta and grey; data source is the Directorate of Bükk and Hortobágy National Parks). Fertő-Hanság National Park is in North-Western Hungary (Figure 1c with yellow; data source is the Directorate of Fertő-Hanság National Park) just like Hegyhátsál (46.96° N, 16.65° E; Figure 1c with red). Agricultural parcel boundaries were taken from the parcel-level farmers’ claim data submitted to the Common Agricultural Policy (CAP) subsidy system in Hungary (referred as Anonymised Database of Subsidies; data courtesy of the Hungarian State Treasury). Mowing dates and mowed bale counts came from management registers. It is common for multiple plots to be recorded with the same yield, which we tried to filter out (e.g., by merging neighboring plots with identical yields; Figure S1a in the Supplementary material) but this may still cause some systematic errors. Unfortunately, we were not able to fine-tune the boundaries of the plots and the unmowed areas based on satellite data everywhere, but in our experience, this can only have a minor impact on the results (Figure S1b in the Supplementary material). Yield by weight calculated using the standard value of 300 kg bale^-1, as proposed by the National Park or by the farmer, except for Fertő-Hanság, where data is based on precisely measured bale weight.

2.2.3. Biomass Based on Bale Digitalization Using Orthophotos

Bale numbers have been determined based on orthophotos to nearly 120 areas by the authors. The orthophotos for the year 2023 cover the entire country with a spatial resolution of 40 cm (data provided by Lechner Non-profit Ltd.; Figure 3a). We determined the parcel boundaries mainly based on orthophotos, but later refined them using Planet satellite images, taking into account simultaneous mowing (Figure 3b, c). The size of bales was measured on stereo photos (Figure 3d), then two types of bales were identified based on the dimensions. Weight was calculated based on estimates made by farmers: 300 kg bale^-1 for 125 cm diameter and 400 kg bale^-1 for 150 cm diameter. It should be noted that counting bales, the delineation of parcel boundaries, the determination of bale height, and the uniform calculation of bale weight may all contain errors that accumulate.

2.2.4. Yield Data from the Hungarian Central Statistical Office (HCSO)

The territorial statistics of the HCSO on agricultural production are based on the collection of data on actual use provided by farm organizations and individual farms according to its own methodology [32], and annual grass hay yields are reported by county (NUTS3 level) in tons per hectare. The definition of grassland used by the HCSO and the methodology of data collection are explained in detail below (see Section 2.4.3.). Note that the data was obtained through personal communication.

2.3. Remote Sensing Information

2.3.1. Copernicus LAI

The Copernicus Global Land Service (CGLS) is earmarked as a component of the Copernicus Land Monitoring Service (CLMS) to operate “a multi-purpose service component” that provides a series of bio-geophysical products on the status and evolution of land surface at global scale [33]. The CGLS provides a collection of 300 m LAI, FPAR and Fraction of Green Vegetation Cover (FCover). These products were produced and delivered from Proba-V Vegetation (VGT) data from January 2014 to June 2020, and from the Ocean and Land Colour Instrument (OLCI) on the Sentinel-3 platform since July 2020, ensuring continuity of the 300 m products. The variables are calculated globally on a 10-day basis, provided in netCDF4-CF format containing the variables values (LAI, FPAR or FCover), associated with some quantitative (QA) and qualitative (QC) quality indicators.

LAI is an integrative measure of plant development stage as it quantifies the total leaf area per unit surface area. Higher leaf area is typically associated with larger biomass in herbaceous vegetation, and can indicate larger photosynthetic capacity and biomass accumulation. This is an intrinsic canopy variable that should not depend on observation conditions. Note that vegetation LAI as estimated from remote sensing includes all the green contributors such as the understory when existing under forests canopies [33].

In our research, the observation-based mowing data always represents the first mowing/sampling during the year, while the highest value of LAI, representing the peak biomass of the grassland, is highly variable. In some cases, there is only one green-up period (referred here as season 1 associated with the first half of the year), and then LAI drops and remains low until the end of the year (Figure 4a).

In the majority of the cases, we observe a second green-up period in the second half of the year (referred here as season 2; Figure 4b-f). Sometimes even the highest value of the entire annual LAI time series is reached during season 2 (Figure 4e). Therefore, in this study, we examined the maximum values of the LAI time series divided into two seasons: before and after August 1 in each year. We refer to this data as MAX LAI (represented by 2 values).

Due to data availability and the applied smoothing of the LAI time series, the LAI value does not decrease significantly after mowing (Figure 4d-f). We had to keep this in mind during the post-processing of observation data (detailed in Section 2.5.).

2.3.2. High-Resolution Vegetation Phenology and Productivity (HR-VPP) Product

The CLMS disseminates the so-called High-Resolution Vegetation Phenology and Productivity (HR-VPP) product for Europe from 1 January, 2017 onwards [34]. This phenology related product is derived from Sentinel-2 images using information about the photosynthetically active green leaf development stage. HR-VPP is available at 10 m resolution over the entire EEA39 (32 member countries, the UK and 6 cooperating countries in the Western Balkans). The layers of vegetation phenology and productivity parameters (VPPs) are derived from seasonal trajectories of PPI [35], since PPI was found to have the strongest relationships with photosynthetic phenology as represented by annual gross primary production (GPP), and it has the highest robustness against background noise in snowy conditions [36]. In their research, Tian et al. [36] examined the performance of NDVI, EVI2, and PPI derived from Sentinel-2 data for Europe-wide phenology mapping, and among the three parameters, they recommended the PPI index for VPP derivation.

The VPPs and corresponding quality layers are generated using the Fortran version of TIMESAT 4 software package, developed by Lund and Malmö University, for analysing time-series of high-spatial resolution satellite sensor data [37,38]. Since all satellite images from the year are used for processing, the derived VPPs are published annually, after the end of the growing season(s).

Of the published quantities, the maximum PPI index [39] and total productivity [40] (referred as MAX PPI and TPROD, respectively) were used. TPROD is the integral of PPI that is the sum of all daily PPI during the growing season. The latter (being the growing season integral computed as the sum of all daily values minus their base level value [34]) can be considered as proxy of potential biomass accumulation. The quality flag layer provided with the VPPs were applied to preliminary filtering: only pixels with high or medium reliability (QFLAG >6) were used.

Note that in HR-VPP the number of output seasons per year was limited to a maximum of two per year according to the TIMESAT logic, of which both were taken into account in this study. According to our experience the presence and significance of the second season in our country is relatively small compared to the first season in most of the years (Figure 5 and Figure 6).

2.3.3. MODIS NPP

Yearly NPP estimates derived from Terra/MODIS (MODerate Resolution Imaging Spectroradiometer) measurements were used as supplementary information, based on the MOD17 gross and net primary production (GPP/NPP) model [41]. The official MODIS MOD17A3HGF version 6.1 products [42] were downloaded for Hungary for the period 2001-2024, providing annual NPP estimates in gC m⁻² year⁻¹ at a spatial resolution of 500 m. These gap-filled annual NPP values are generated from the 8-day GPP.

When calculating annual, country-mean NPP, only pixels with at least a 90% grassland cover were used in the analysis, based on National Ecosystem Map of Hungary [43] (detailed later in Section 2.4.5.). Resampling this dataset to the MODIS grid resulted in values representing the proportional coverage (%) of each ecosystem category within individual 500 m MODIS pixels [44]. The summed proportion of all grassland categories was used for pixel selection.

2.4. Estimation of the Spatial Extent of Grasslands in Hungary

The estimation of the grassland area of Hungary was based on different data sources. According to Belényesi et al. [17] the total grassland area identified differs largely depending on the methodology and grassland definitions. We would like to emphasize that it was out of scope of the present research to decide which database is more reliable when it comes to grassland areas, so below we present statistics from several different sources. Note that quantification remains a major source of uncertainty in the study (that will be discussed later).

2.4.1. CORINE Land Cover (CLC)

The CORINE (Coordination of Information on the Environment) program was launched by the European Commission, and now the inventory is one of the main products of the European Environment Agency’s CLMS [45]. The aim of the program was to develop a standardized methodology for producing continent-scale land cover, biotope, and air quality maps. The well-known and widely used CORINE Land Cover (CLC) product offers a pan-European land cover and land use inventory with 44 thematic classes. The dataset has a minimum mapping unit (MMU) of 25 hectares for areal phenomena and a minimum mapping width (MMW) of 100 m for linear phenomena and is available as vector and as 100 m raster data. The product is updated with new status and change layers every six years – with the most recent update made in 2018, so here we computed statistics from the status layer of 2018. Due to MMU and MMW, polygons which consist solely of grassland are rare. Comparative analysis of CORINE Land Cover categories and other Copernicus products (such as High Resolution Layer Grasslands [46], Herbaceous cover [47] and CLCplus Backbone [48]) was carried out (Table S2 in the Supplementary material). This analysis showed that categories Natural grassland (321) and Pastures (231) contain the “purest” grasslands, but other categories like Airports (124) and Sparsely vegetated areas (333) also contain grassland areas. Here we considered only the 321 and 231 categories as grasslands.

2.4.2. Copernicus High Resolution Layer (HRL) Grasslands

The CLMS HRL Grasslands map the location and size of permanent and temporary grasslands, in addition to providing information on mowing events across Europe [46]. Key application areas of the products range from implementation or monitoring of various European Union policies such as the Habitats Directive or the CAP, to contributing to research on how climate change is impacting grassland ecosystems. The grassland status layers provide a basic land cover classification with two thematic classes (grassland / non-grassland), produced yearly, and has spatial resolution of 20 and 100 m for 2015, and 10 and 100 m from 2017 onwards. The change layers that provide information on changes in grassland vegetation cover between years are published every 3 years.

2.4.3. Grassland Area Data from HCSO

The regional reports on agricultural activities published by the HCSO are based on data collected using its own methodology from agricultural organizations, individual farms, and selected individual farms [32]. Data on actual use are collected using questionnaires and supplemented with data from other statistical surveys.

The methodological descriptions of data collection clearly define the concept of “grassland” in the HCSO nomenclature: only grassland that has been used for at least five years for grazing animals (pasture) or mowing (meadow), or grassland that is not used but has received support for being kept in good environmental and ecological condition, is considered grassland. Certain areas are not recorded as grassland but as arable land. These include areas covered with herbaceous vegetation less than five years old (but at least one year old) that are part of crop rotations, as well as areas sown with grasses for mowing (hay) or silage, which are considered temporary grassland. Plant mixtures planted for bee pasture, plants sown for green manure, and green fallow are also classified as arable land. All this information must be taken into account when comparing the HCSO’s land use data with data from other databases [17].

In 2020, methodological changes were made to the HCSO’s data collection, including the setting of a new economic threshold for land use [49]. To ensure comparability, the data for 2010, 2013, and 2016 were also recalculated. In this work, we present data calculated using the new methodology.

2.4.4. National Agricultural Map (NAGRI)

In Hungary, the Space Remote Sensing Department of the Lechner Knowledge Centre performs a reliable mapping of the entire territory of Hungary which annually carry information on grasslands and agricultural areas (called National Agricultural Maps; [50,51]). One of the main objectives of mapping is the delineation of permanent grassland (i.e., areas with confirmed presence through a period of 5 years). These thematic maps for several years are based on a standardised methodology using cloudless Sentinel-2 images [52], spectral indices generated from them and temporal integrals derived from Sentinel-1 radar images [52,53]. The core of the methodology is a standard Random Forest algorithm [54] implemented in a Python environment [55]. Most of our training and test data is from the Anonymised Database of Subsidies (data courtesy of the Hungarian State Treasury) indicating crop type and submitted to CAP system. The thematic resolution of our maps is increased to 14-23 categories by integrating polygons for forests, water and detailed grassland categories (i.e., non-cultivated, saline, wet, herbaceous or woody grasslands affected by weeds) into the reference data system. After classification, irrelevant areas such as buildings and roads are masked based on ancillary GIS data. The satellite data were resampled uniformly to a 20×20 meters EOV (EPSG:23700) grid, so the spatial resolution of our result maps is also 20 meters. The frequency of mapping is once a year, in October, at the end of growing season. The accuracy of the results was measured with the built-in overall accuracy (OA), cohen kappa (Kappa), precision, recall and F1-score metrics of the Scikit-learn Python package [55]. The average value of OA averaged in space is 95.20% (standard deviation = 1.52), while the mean values of the Kappa are 93.97% (standard deviation = 0.89), and the mean precision, recall and F1 score values of the grassland category are 85.94%, 94.71% and 90.03%, respectively [51]. Although the maps produced each year are highly accurate in themselves, the margin of error is comparable to the changes that occur between years, so the difference between the two layers cannot be considered an unambiguous change in surface coverage.

2.4.5. National Ecosystem Map (NECOMAP)

The National Ecosystem Map of Hungary [43,56] was prepared to fulfill the obligations related to the European Union Biodiversity Strategy, Mapping and Assessment of Ecosystems and their Services (MAES). The map shows the spatial distribution of ecosystems in Hungary using a three-level category system in the form of a thematic raster with a spatial resolution of 20×20 meters for the base year 2015/16. The information is supplemented by four additional layers, thus helping (i) to distinguish between trees and shrubs, (ii) to distinguish vegetation-covered areas in urban areas, (iii) to identify saline lakes, and (iv) further enhancing the base map’s Level 3 forest categories with Level 4 categories [57]. The map is primarily based on the processing of remote sensing data from 2016 and, to a lesser extent, from 2015 and 2017 using the Random Forest algorithm [54] and other regularly updated databases. During processing, time series Sentinel-1 and Sentinel-2 satellite images [58], spectral indices and radar characteristics derived from them were used. This was followed by a rule-based integration of official databases such as Land Parcel Identification System (LPIS), road, railway and surface water coverage, national soil and forest databases. The results underwent multi-level quality control [59,60]. The map currently is a status layer, with the future goal of extending it to monitoring between years. Here we present the sum of Grasslands and other herbaceous vegetation (3) and Fens and mesotrophic wet meadows, grasslands with periodic water effect (5120) classes.

2.5. Post-Processing of Observation Data

Due to the well-recognized ambiguity of productivity and NPP [61], first the nomenclature must be set, and the methodology must be properly described that is used to approximate ANPP and NPP. In case of heterogeneous data sources that are used in the study, the different biomass metrics must be harmonized.

Most of the data sources detailed in Table 1 provide estimation on harvested aboveground biomass (HAB), and it is used to approximate the total aboveground biomass (AGB, including standing dead biomass; [15]). In grasslands, HAB and AGB, estimated after the intensive spring growth, are widely used to estimate ANPP, that is aboveground net primary production. ANPP, by definition, is the growth (i.e., biomass increment) of the aboveground plant parts within a year. In our study HAB is considered to be equal to AGB in case of complete clipping (when all biomass is removed from the target plot). If the HAB is the result of mowing, additional considerations are needed (e.g., [62]). Grass mowing with machinery always leaves a fraction of the total biomass at the site to support regrowth (the cut part is left at the site for a few days to dry out, then it is packed into bales and removed from the site). The residue is referred to as stubble in some studies. In this sense HAB is not equal to AGB, thus NPP is not approximated well.

In order to use a harmonized dataset, correction was applied to the HAB data observed in situ or estimated from other methods. Using the available LAI dataset, we estimated the change in LAI (that is considered as a proxy for biomass) before and after harvest. LAI before harvest was approximated by the maximum LAI in the first half of the year. As the relatively smooth LAI dataset distorts (masks) the immediate LAI after mowing, we uniformly used LAI=1 as an approximation for the post-harvest LAI. This is supported by ceptometer based measurements from the Hegyhátsál grassland experimental site in Hungary (unpublished data; see [30,31] for site description). The change in LAI was associated with HAB and the total LAI was used to estimate ABG. This procedure ensured the comparability of the data and the homogeneity of the entire dataset (meaning as homogeneous as possible in spite of the uncertainties) that is used for the model construction and upscaling. In this sense, in the Results section HAB always implicitly means corrected biomass observation.

As peak biomass observations are available in Hungary for the first mowing (that is the only cut in some regions), ANPP was estimated first for this time period (typically around the beginning of June). After mowing the regrowth can fail in some regions during summer due to lack of precipitation (Figure 4). However, in autumn, after the occasional summer heat and drought, there is a secondary growth with typically smaller biomass production. As this secondary growth contributes to the annual ANPP, we used the relationship established based on the first cut, and extrapolated it to the second one using the remote sensing data and the selected VI/biophysical variable (MAX LAI, MAX PPI or TPROD) in the second half of the year. Annual ANPP was estimated as the sum of the two simulated AGB.

It is important to note that as most of the observation-based data reported air-dry biomass instead of dry matter (DM; the latter is called oven-dry, that is the biomass after its water content is removed typically by heating), the results provide HAB and ANPP with some water content. In the interpretation and application of the data this should be kept in mind. In order to keep consistency, we intentionally provide the results in gBM unit (grams of biomass; and also in MtBM unit for the country totals that is megatons of biomass), since we refer to biomass that can be considered as air-dry. In this sense the data provided here in gBM unit differs from the data in gDM units that are widely used in the literature. It was out of scope of the present study to provide precise conversion factors that can be used to estimate dry matter, but the conversion can be done with approximations (e.g., assuming 14% water content).

2.6. Statistical Analysis

Simple linear models were used to establish relationship between the observation-based field evidence (HAB) and the remote sensing based biophysical data (MAX LAI, MAX PPI and TPROD). The plot- and county-wise datasets were combined with three different logics to support probabilistic ensemble estimation. Least absolute residuals (LAR) fitting was used as it is less sensitive to the presence of outliers and it provides a more robust fit. The calculations were made in MATLAB [63].

During the calculations, critical consistency check was performed to set the baseline for the model construction. The minimum expectation for a given dataset was that clear (but not necessarily strong) relationship was sought between the in situ data and the remote sensing proxy. In other words, we considered a dataset as usable (or informative) if small values of the observation data correspond to small values of the remote sensing data and vica versa. If this condition is not fulfilled, then the dataset can be dropped. In case of linear regression-based model evaluation normalization was also done since the predictors introduced above have different magnitudes. This normalization enabled the direct comparison of the sensitivity of the linear models (i.e., comparison of the slopes) that could be used to further perform the consistency check.

After the calculation of the regression equations, the established relationships were used to estimate total aboveground biomass production for the whole country, using an ensemble of equations and vegetation indices. In this sense we quantify the uncertainty related to the reference dataset selection, and to the uncertainty that is arising from the grassland area identification. Interannual variability of grassland production is also addressed.

Assumptions on belowground productivity (belowground net primary production, BNPP) supplements our work, which means that total NPP is estimated based on ANPP. As NPP is one of the most important measures of the carbon cycle of grasslands, its quantification is essential given the fact that in grasslands root production is larger than shoot production [64,65]. The probabilistic approach for HAB and ANPP estimations supports the ensemble (more robust) estimation of the NPP of Hungarian grasslands.

We used five different conversion factors to estimate BNPP/ANPP (and also BNPP/NPP and ANPP/NPP for easier interpretation) for the NPP calculations (Table 2). The annual GPP sums from the Bugac [30,66] and Hegyhátsál [30,31,67] eddy covariance (EC) sites were used to estimate NPP as ~GPP/2, that is a robust estimation used in the literature ([68]; note that using alternative conversion factors was out of scope of the study). Using the biomass observations at Bugac and Hegyhátsál (Table 1) the carbon equivalent of the biomass was calculated, then ANPP was calculated, and the NPP component ratios were estimated (Table 2). We used the Mokany et al. [69] study’s mean estimation for the root:shoot ratio for temperate grasslands (that is BNPP/ANPP=4.22) as an alternative conversion factor (the other ratios were also calculated). Hui and Jackson [70] provided explicit equations for the BNPP/(BNPP+ANPP) ratios based on a synthesis of field experiments for global grasslands. They proposed equations for estimating the ratio from long term mean of precipitation and temperature (see Figure 1 and Figure 2 in [70]). We used those equations driven by the long-term climate data for Hungary to estimate the BNPP/(BNPP+ANPP) ratios provided in Table 2. This ratio was converted to the ratios presented in Table 2. Using the ANPP estimations based on the presented linear models the ANPP/NPP fraction was used to estimate total grassland NPP.

Given the fact that NPP is calculated from ANPP by using the presented constants, the spatial pattern of NPP clearly follows that of ANPP, thus spatially explicit visualization is not performed. Instead, we present NPP per unit area, that, of course, can be upscaled to the whole country with a selected area estimation (e.g., NAGRI). Calculation of NPP per unit area establishes a way to compare the findings with NPP taken from previous studies.

3. Results

3.1. Grassland Area Estimation

Figure 7 summarizes the grassland area estimations for the whole country, calculated from the spatially explicit datasets detailed in the Material and methods section. The figure clearly shows that there is a large variability between the estimations.

The highest grassland area was identified by Copernicus HRL, covering ~1.35 million hectares in 2023, which is roughly 14.5% of the country’s territory. The variability between years shows that the total area of grasslands is decreasing. 2015 deviates from this pattern, which can be explained by a methodological change in mapping. The decrease between 2017 and 2023 is 198,550 hectares, which represents nearly 13% of the grasslands in 2017. NAGRI also shows a steady decline in total grassland area (2019 is an exception), as it decreased by 15.8% (222,000 hectares) between 2015 and 2024. The difference between the years might be associated with classification errors to some extent. The NECOMAP mapped 25% less grassland than the corresponding NAGRI data. This is probably associated with the fact that this dataset has a more detailed categorization. The CORINE Land Cover map for 2018 provides even smaller area, with less than 1 million hectares. This can be explained by the fact that we only used two categories from CORINE. In this dataset the spatial resolution might also restrict some additional grassland mapping. The smallest grassland area (about 50% of Copernicus HRL or NAGRI) is estimated by HCSO statistics and the trend shows increasing of the grassland territory. The reason for this is likely the applied definition of grassland, as explained in detail in the Metrials and methods. Supplementary material Figure S2 shows an example to demonstrate the differences between the grassland categories in the databases we used.

Based on multi-year averages where data was available, the total grassland area varies significantly across sources: Copernicus HRL reports the highest coverage at 1,426,229 ha, followed by NAGRI (1,317,885 ha), NECOMAP (1,041,958 ha), CORINE (923,340 ha), and HCSO (754,720 ha). The difference between the minimum and maximum multiannual average is 671,509 ha.

3.2. Relationship Between Observation-Based and Remote Sensing Data

We examined the relationship between the biomass observations and the proxies derived from satellite images using three data selection approaches. Selection A (SELA) contains all the biomass data at our disposal, including point, plot, and county-level HCSO data, while Selection B (SELB) only includes data collected by mowing in parcels, while Selection C (SELC) examines HCSO data alone as it has large spatial representativeness (Figure 8, Figure 9 and Figure 10). We also examined the data sampled using the cutting method at specific points (Fülöpháza, Bugac, and Orgovány) separately, but since these samples were taken in areas with extremely low yields and thus the national extrapolation of the models fitted to them cannot be considered representative, we do not present them here as a separate selection group.

Biomass data from all three data selections were plotted against both Sentinel-3 LAI data (Figure 8) and parameters from the HR-VPP database (MAX PPI and TPROD; Figure 9 and Figure 10). The vertical axis of the graphs in Figure 8, Figure 9 and Figure 10 shows HAB, while the horizontal axis shows the parcel-mean (or NUTS3 level mean in case of HCSO reference data) of the satellite proxy. Since there is no uniform number of samples available for each year, it is not feasible to examine each year separately; the data for each year are shown together on the plots.

It is clear from the plots that the HAB values remain below 1200 gBM m^-2 in all cases, and the data derived from orthophoto digitization (blue dots) have the largest variability. The data from the point measurements mentioned above, as well as a significant portion of the data from the Hortobágy and Bükk National Parks, are associated with the lowest values, below 200 gBM m^-2 (grey and magenta dots).

The consistency of the MAX LAI-based relationship is surprisingly low and it is rather noisy. In the case of SELA and SELB, they do not clearly align with the fitted line, and there is also a noticeable plateau between LAI 5 and 6 m² m^-2. In the case of SELC, i.e., when using the HCSO data, the average MAX LAI value of some higher-yield counties is around 2 m² m^-2, which results in a very small slope of the fitted line. In the HR-VPP-based comparison, the TPROD graphs show better consistency, and the R² values are also the highest here (0.41 in the case of SELA).

The fitted regression equations are summarized in Table 3. The slope and intercept values are not comparable due to the markedly different magnitudes of the satellite proxies (LAI varies from 0 to 7, MAX PPI ranges between 0 and 3, while TPROD is between 0 and 400).

When the fitted lines are plotted together after normalization (Figure 11; see Section 2.6.) the slopes of the lines turned out to be comparable, from which the model based on HCSO (SELC) and MAX LAI values deviated significantly. Considering the consistency criteria described in Section 2.6. and also the large deviation from the other slopes, we excluded this equation (marked by red in Table 3) from the analysis and did not use it in the national scale extension. In summary, we used 8 equations to upscale the observations to the country level (black, blue and green equations in Table 3).

3.3. Country-Level HAB and ANPP

The eight equations presented above were used to derive a national HAB for each year, but only one of them (the one based on TPROD and SELA) is presented here for two extreme years (Figure 12), as well as the average and standard deviation between years (Figure 13). Supplementary material Figure S4 shows maps for all examined years from 2017 to 2024.

The results were produced at a spatial resolution of 10 meters, and integration into hexagons with diagonal length of 10 km was performed solely for visualization purposes. Note that we colored all hexagons where the proportion of grassland was at least 10%, but obviously only grassland pixels were taken into account in the spatial aggregation. When deriving average yields (Table 4 and Figure 14) and later country totals, we always used the highest available spatial resolution: 300 m for LAI and 10 meters for HR-VPP. In all cases, the maps represent the sum of the two seasons, i.e., the yield (ANPP) for the entire year.

The aggregated maps for 2022 and 2023 show a difference in yield between an extremely dry year and a year with higher precipitation (Figure 12). In both years, yields in the driest areas were less than 200 gBM m^-2 (Figure 12, yellow color). The maximum value of ANPP in 2022 is 128 gBM m^-2, and there is no hexagon where it exceeds 500 gBM m^-2. In 2023 the minimum of the average ANPP values are still below 200 gBM m^-2, but in many cases exceed even 500 gBM m^-2.

The long-term average map shows that the central, lowland, and sand ridge areas of the country have lower, while the western and northeastern hilly and mountainous regions have higher average yields (Figure 13).

The long-term standard deviation of ANPP characterizes each area in terms of stability of annual production. It is interesting to note that the spatial correlation between the mean and the standard deviation is non-significant, meaning that there is no relationship between mean biomass production and its temporal stability.

The ANPP calculated from the models (Table 3) is averaged annually for each model and also averaged for the whole time series (Table 4 and Figure 14). It is striking that models based on MAX LAI estimate significantly higher yields than those based on HR-VPP values. The large difference between 2022 and 2023 is most evident in the models based on TPROD, with a difference of 95 gBM m^-2 between the average yields for the two years. Using the estimates of all models, the multiannual mean grass yield in Hungary is 325.7±56.6 gBM m^-2 (median is 316 gBM m^-2).

The individual ANPP time series (from 2017 to 2024; see Figure 14) estimated by the different proxy-data selection combinations were compared with the time series of the annul means and medians to propose possibly one single method that approaches the ensemble relatively well. Considering the annual means (last row in Table 4; upper values) the MAX PPI-SELA combination was associated with the lowest bias (-1.6 Mt year^-1) with relatively high explained variance (R²=0.85). The highest explained variance was associated with the MAX PPI-SELB combination (R²=0.91). TPROD-SELB had the second smallest bias (-13.2 Mt year^-1). Considering the median time series, again the TPROD-SELB combination seems to be suitable to represent the robust estimation if bias is the main criteria (bias=0.4 Mt year^-1, R²=0.74). The TPROD-SELA combination is associated with the highest explained variance (86%, bias= 4.9 Mt year^-1). Depending on the planned application, the combination with the highest R² or with the smallest bias (in absolute sense) can be this selected.

Using the average annual yields (Figure 14) and the grassland area calculated based on various grassland databases (Figure 7), national grassland yield totals were obtained, which are shown in Figure 15. Although grassland area showed a declining tendency in several databases (Copernicus HRL and NAGRI; Figure 7), this clear trend is not evident in the national grassland yield, which can be explained by changes in average yields between years. If we use each data selection-grassland area database combinations as individual estimation, the annual averages vary between 31.9 and 38.3 Mt year^-1, with a multiannual mean of 35.1±9.2 Mt year^-1 (median is 35.9 Mt year^-1) biomass.

We examined whether it is the grassland area estimation or the model construction that has larger influence on the country totals. Using multiannual mean ANPP data, Table 5 shows that the standard deviation of the results that are caused by the grassland area estimation varies between 7.6 and 9 Mt year^-1, while it ranges between 3.1 and 5.5 Mt year^-1 caused by the model selection method (considering only a single grassland area database). Results focusing on individual years are presented in Table S5 in the Supplementary material, which corroborate these results (with some interannual variability in the standard deviations). The results indicate that the grassland area database selection has larger influence (with almost a factor of 2) on the country totals than the remote sensing based linear model selection.

The country total time series, using all 40 individual estimations for the different proxy-data selection-grassland area combinations (Figure 15) were compared with the annul means and medians to propose possible simpler method that is closest to the ensemble mean and median. Considering the mean time series, the NECOMAP-MAX PPI-SELB combination is the closest to the robust estimation (bias=-1.3 Mt year^-1, R²=0.82). Using the median of the 40 individual estimations for the country totals per year, the NECOMAP-TPROD-SELB (bias=0.9 Mt year^-1, R²=0.82) combination and the NECOMAP-TPROD-SELA (bias=0.5 Mt year^-1, R²=0.81) combinations seem to be a good choice.

If we combine the results from the country-mean ANPP and the country total biomass production, and selecting the median as the proposed approach in both cases, the NECOMAP-TPROD-SELA combination seems to be an optimal choice.

3.4. Estimated Total NPP for Hungary

Table 6 summarizes the mean grassland NPP (2017-2024) for Hungary based on the different models for ANPP and also based on the different conversion factors between ANPP and NPP (described in Section 2.6.). The estimations span a wide range from 313 gC m^-2 year^-1 to 773.5 gC m^-2 year^-1 where the high end is associated with the M2005 conversion factor. Across all predictors and data selection methods the HJ2006TEMP conversion factor provides the smallest estimates (366 gC m^-2 year^-1). Based on the MAX LAI, MAX PPI and TPROD proxies the overall NPP mean is 551, 403 and 421 gC m^-2 year^-1, respectively. The variability of NPP (expressed as the standard deviation of the estimates) is larger due to the selection of the conversion method than due to the selection of the predictor dataset.

Supplementary material Figure S5 shows the frequency distribution of the estimated NPP values from Table 6. The distribution is skewed which means that the mean is not the best predictor for robust estimation. Nevertheless, the long term mean of grassland NPP in Hungary is 447 gC m^-2 year^-1 (with standard deviation of 110 gC m^-2 year^-1). The median value is 421 gC m^-2 year^-1.

Figure 16 shows the annual sums of NPP (per unit grassland area) for Hungary based on the ensemble mean and median using the NPP data calculated from all combination of the predictors (Table 3) and the NPP conversion factors for the entire study period per year.

The figure also shows the estimated NPP for the Bugac and Hegyhátsál eddy covariance sites for comparison (Figure 16). Note that the Hegyhátsál eddy covariance system was disassembled in 2019, so no further data is available for GPP. The NPP calculated in this study shows relatively stable performance, at least compared to the point level eddy covariance based NPP. The Hegyhátsál data shows consistently higher values than the best estimate from the presented modelling exercise.

Interestingly, up to 2020 the Bugac NPP was similar to the one observed in Hegyhátsál in spite of the fact that it is located in a sandy, drier region. In the same time, the NPP derived from this study was considerably lower. However, after 2020 NPP dropped at Bugac, and shows very similar NPP sums than the presented country-mean results. It suggests that the Bugac site might be considered as a sentinel for the whole country, at least in the 2021-2024 time period. Data from additional years are clearly needed to support the statement.

4. Discussion

4.1. Grassland Area Estimations

Grasslands are highly relevant ecosystems both globally and locally due to their ecosystem services, which include biodiversity conservation, water purification, and carbon sequestration, just to name a few. Another important role is their agricultural use and contribution to food security [2]. Planning effective grassland management, setting the right mowing intensity, and estimating expected yields are all important issues in terms of milk, meat, and wool production worldwide [5]. Both the preservation of natural values and consistent agricultural use require knowledge of the spatial extent and location of grasslands.

The first and perhaps the most complex issue of grassland mapping is the definition of grassland itself. There are so many heterogeneous and diverse grassland types worldwide that it is not obvious to come up with a general definition [19]. This difficulty is further complicated by the dynamic variability of grasslands within a year, as well as the transition between grasslands and croplands or wooded/forested areas or wetlands between years. This variability can be caused by natural factors (climate, water supply, erosion, wildfires) and/or human intervention (grazing, mowing, irrigation, fertilization, sowing, cultivation). These factors demonstrate how complex it is to create a spatially explicit database that can characterize grasslands from all aspects.

In many European countries, including Hungary, there is no database focusing exclusively on grasslands that could serve both nature conservation and other grassland management and utilization planning needs and considerations at the national level, especially those related to forage production. The compilation of European best practices is currently underway as part of the LIFE IP GRASSLAND-HU project [71], the full results of which will be published later by researchers working on the topic. We present here a summary of the first results based on [17].

According to the finding of the experts, there are few mapping efforts specifically aimed at assessing the area and condition of grasslands; rather, we generally find studies covering several habitat or land cover types, much of which has been/is being carried out only in sample areas at the level of individual projects. In Slovakia, between 1999 and 2006, field sampling was carried out at 20,000 sites under the supervision of the DAPHNE nature conservation association [72,73]. Later, a detailed ecosystem map of Slovakia was also completed [74]. The Romanian Grassland Database is not the result of a systematic mapping effort, but rather a collection of surveys conducted in 21,685 quadrats across the country [75]. The Baltic states (Estonia, Latvia, and Lithuania) launched a project in 2014 (called LIFE Viva Grass) to study and support the maintenance of biodiversity and ecosystem services provided by grasslands [76]. This project was implemented in three sample areas per country covering a total of 465,600 hectares.

Exceptions to project-level surveys are the Czech Republic, where an annually updated national database is available for all areas of environmental value [77,78], and Ireland, where a national, high-resolution land cover map was prepared under the coordination of an interdisciplinary working group [79,80]. The use of remote sensing in mapping methodologies was diverse, but only Ireland used advanced, rule-based and machine learning classification methods, including artificial intelligence, for the semi-automated processing of aerial and satellite imagery (Sentinel-2). This system also integrated several other databases available at the national level into the processing workflow.

In Hungary, currently there is no national level, spatially explicit database providing comprehensive information on the location, extent, type, ecological status, and utilization of grasslands and grassy habitats [17]. This dataset would be highly useful for ecologists and nature conservation experts, as well as for professionals in the agricultural sector involved in animal farming and feeding [81].

The mapping method based on the General National Habitat Classification System (ÁNÉR; [82]) would be suitable for the high-precision mapping and classification of Hungarian grasslands, but it involves intensive fieldwork, so its nationwide application would be extremely time-consuming and costly. This mapping method is used to produce the habitat maps of the National Biodiversity Monitoring System in 5×5 km squares every 10 years (NBmR; [83]), as well as the mapping of Natura 2000 sites, which cover approximately 21% of the country’s area. This method was also applied between 2003 and 2006 to create the Hungarian Habitat Map Database (MÉTA; [84,85]), which covers the entire country; however, the mapping units were not based on actual land cover or habitat boundaries but based on 35-hectare hexagons referred as “MÉTA” hexagons.

In 2022, the Grassland Management Working Group of the Hungarian Livestock Breeders Association and the National Chamber of Agriculture conducted a sample-based field survey in several locations across the country on the condition of grasslands, which was structured around livestock feeding considerations [86,87]. The experts carried out their work on an approximately 20,000-hectare area in various parts of the country. In addition, information on a further 88,404 hectares was collected using questionnaires, which asked about the size and yield of grasslands and bulk fodder production areas used by farmers, the types and application of grassland management, the number of mowing, grazing methods, the number of rotations, the length of the grazing period, the species and number of animals grazed, and the use of Ivermectin. These grasslands represent 11.5% of the total grassland area in the country [87].

Among the available spatial databases at the national level, [17,81] present several in detail with the aim of highlighting their heterogeneity in terms of information content, mapping objectives and methodology, and accessibility. The differences between the various maps lie in the varying (i) mapping objectives, (ii) definitions of grassland, (iii) grassland typology, (iv) data production methodology, (v) spatial resolution of maps, and (vi) temporal resolution and update frequency of data. The comparative study of Hungarian grassland extent covers seven databases for the year 2016 (Figure 6 in [17]).

In the current study, we used some of the datasets presented in [17], and performed the comparison for every year between 2015 and 2024 (depending on the data availability), and included two additional European datasets (CORINE and Copernicus HRL Grasslands) in the analysis like [81] did in their study. As shown in the results section, the area of mapped grasslands shows high variability (see Section 3.1. ad Figure 7), which in our case can also be attributed to the reasons listed above.

Clearly, standardization is a necessary next step in grassland mapping also taking into account the cropland mapping initiatives. For example, set aside can be considered as grassland or cropland as well, but if the definitions and methodologies differ, they can be categorized for both land cover types, or for none of them. Clearly, the total area of the country has to be mapped, and unmapped or double-counted areas are avoidable. This is a challenge for the forthcoming years.

4.2. Reference Productivity Datasets for Model Construction

Communicating grassland productivity is not straightforward despite the apparent simplicity of the term “biomass” or “biomass production”. Several studies discussed the terminologies and measurement techniques that aimed to estimate some aspects of grassland productivity [15,88]. The simplest and most popular method is the destructive sampling of smaller or larger plots (with clipping) and quantifying the dry biomass and/or its carbon equivalent. HAB can approximate the total AGB (including standing dead biomass) or can only approximate some percent of it due to the remaining residue. Some methods differentiate between recently dead biomass (standing dead plants) and old dead biomass which complicates the quantification [15]. HAB and AGB, estimated as peak biomass, are widely used to estimate ANPP. ANPP, by definition, is the growth (i.e., biomass increment) of the aboveground plant parts within a year. Herbivory, exudation and other processes can modulate ANPP, but generally HAB or AGB is a good proxy for ANPP, especially in temperate steppe grasslands [89,90]. ANPP can be expressed as dry matter or carbon equivalent, depending on the research aim. Total NPP, that is ANPP plus BNPP is an important variable that the carbon cycle community uses [91]. Quantification of BNPP is especially problematic since destructive methods are not effective to capture e.g., total fine root biomass. For this reason, assumptions were needed to quantify BNPP from ANPP.

Note that destructive sampling obviously modifies productivity and grassland functioning, as disturbance (i.e., the removal of the cut biomass) changes allocation patterns and results in NPP that is different from the undisturbed grassland. It is fair to state that the productivity of mowed grasslands represents mowed grasslands only, and not grasslands with other management like grazing or intact ecosystems. It would be highly beneficial to standardize the sampling protocols and initiate consistent and representative sampling of biomass, supported by remote sensing data and studies like the present one.

In this study a novel compilation of reference data was presented to support country-average ANPP and NPP. Given the heterogeneous quality of the data, uncertainties are obviously present. Nevertheless, the dataset can be used for further studies to build improved models and to clarify spatial and temporal differences in biomass production. The collected biomass data are comparable to those presented in [81] where the mean ANPP ranged between 112 and 376 gBM m^-2 between 2010 and 2014 (this study reported data between 235 and 470 gBM m^-2).

4.3. Remote Sensing-Based Upscaling of Grassland Production

Several studies have examined the relationship between multispectral satellite data and field-based measurements of grassland biomass production (see [19] for a comprehensive review). The majority of these studies rely on vegetation indices or satellite-derived products specifically designed for vegetation monitoring. With respect to modeling this relationship, both mathematical approaches (using linear or non-linear formulations) and machine learning methods have been shown to yield satisfactory results.

Zhao et al. [92] used MODIS PSNnet to estimate AGB over a large (192,512 km²) Chinese arid and semi-arid temperate continental monsoon region. They used significant amount of in situ measurements from 2005 to 2012, establishing their models using 975 sites and validating them using 230 sites. This study examined multiple mathematical approaches, but the unitary linear method was found to be the most accurate one with the R² of 0.55, RMSE of 26.67 g m^-2 and the precision of 0.69%.

NDVI and derived biophysical variables (LAI and FCover) were used to estimate biomass in the work of Dusseaux et al. [93] on a western French study area, comparing SPOT products to field measurements from 37 sites. According to their results, using linear regression, R² was 0.68 in case of LAI, while field data showed weaker correlation with NDVI (R²=0.3) and FCover (R²=0.5).

In [94], five MODIS vegetation indices (NDVI, EVI2, SAVI, MSAVI, OSAVI) and two raw multispectral bands (RED and NIR) were used to estimate biomass production of Irish grasslands. Their method was applied in two study areas with the total area of 171.3 ha, in two temporal ranges (2001-2012, 2001-2005). Besides multiple linear regression, an artificial neural network and adaptive neuro-fuzzy inference system (ANFIS) were also applied to investigate the efficiency of the machine learning approach. That ANFIS algorithm performed the highest R² (0.85 and 0.76) and the lowest RMSE (11.07 and 15.35 kg DM ha^-1 day^-1) on both study areas.

Quan et al. [95] used a radiative transfer model called PROSAILH to calculate LAI and DM content of grasslands from Landsat 8 data. The result of LAI×DM was regarded as the estimated aboveground biomass. The method was tested on a study area covering 29,600 km² located in Qinghai Province, China. Considering ground truth data, a total of 135 30×30 m plots were sampled in 2014 and 2015. The results of the linear regression showed an R² of 0.64 and RMSE of 42.67 g m^-2.

The study of Guerini Filho et al. [96] was conducted on a natural grassland with the area of 23 ha, located in the Brazilian Pampa. The results showed the maximum R² of 0.61 for linear regression between Sentinel-2 reflectance bands and the 57 samples collected on the field. In this study simple linear models were constructed to approximate HAB and ANPP at large spatial scales, with 10 m spatial resolution, based on MAX LAI, MAX PPI and TPROD. The R² values were moderate (spanning from 0.11 to 0.41) that can be explained by the diversity of the datasets, grassland spectral properties affected by biodiversity and the large spatial scale that the study aimed to explore. Also, as the remote sensing proxy data behaved differently, as their biophysical meaning and limited information content also contributed to the uncertainty.

The comparison of the works [92,94] shows that using the same remote sensing data (MODIS) for smaller regions, better models can be constructed. However, collecting adequate amount of field data with a balanced spatial distribution can also lead to satisfactory results on larger scale [95].

Obviously, more complicated methods are available in the literature for the estimation of biomass and ANPP using remote sensing and in situ observations [19]. In our case the selection of the simplest method is justified by the fact that we aimed to understand the effect of the selection of the driving remote sensing dataset on the applicability as a predictor. Since there is no similar previous study for Hungary up to the knowledge of the authors, this first yet simple method is fully supported. Forthcoming studies can focus on more complex methods like machine learning, XGBoost and others that could possibly provide refined estimates, but this was out of scope of the present study.

The results clearly indicate that now the remote sensing-based grassland mapping is the largest source of the uncertainty in the production estimations. This clearly claims for a harmonized definition of grasslands in Hungary and in the European Union, and the construction of grassland cadastre/register database that support policy, environmental protection initiatives and animal husbandry [17,81].

Similarly to croplands, there is a major need to collect and synthesize grassland biomass data for the whole country. This can be done with coordinated efforts aiming at the construction of spatially explicit, freely available datasets. Methodology for grassland biomass sampling must be harmonized to avoid complications that are related to the conversion factors, residual stubble and other methods [97].

In a wider context, other variables needed to be included in the analysis. Clearly, remote sensing data is affected by the soil conditions, especially in open grasslands that occupy relatively large areas in Hungary. Topography, soil water content, species distribution and other information can supplement future studies.

4.4. Estimation of Total NPP

Explicit NPP estimations for grasslands are relatively rare in the literature for countries or for larger regions. It is not surprising given the methodological difficulties in the observations and the modelling approaches. Note that grassland models are still highly uncertain in terms of their overall performance and their sensitivity to climate change factors that makes the estimations even more difficult [98,99].

In an early yet pioneering study Vuichard et al. [100] estimated the NPP of European grasslands with the PaSim model considering fertilization and management effects. The estimated values ranged between 600 and 1000 gC m^-2 year^-1. The study investigated two conditions in the NPP simulations, the “cut grassland” and the “grazed grassland” scenarios. The “cut grassland” scenario showed high NPP values in the western and southern part of Europe (almost 1100 gC m^-2 year^-1), 900 gC m^-2 year^-1 in the middle of Europe, while between 600 and 800 gC m^-2 year^-1 in the east part of Europe. In the mountainous grassland areas in Central Europe the NPP values were less than 600 gC m^-2 year^-1. The “grazed grassland” scenario showed consistently higher overall NPP values (with difference around 300 gC m^-2 year^-1) compared to the “cut grassland” scenario in all parts of Europe.

Schulze et al. [8] published the main findings of the CARBOEUROPE-IP project considering several components of the European carbon and greenhouse gas balance, representing the summarized knowledge in the early 2000s. According to the results based on atmospheric observations and ground-based measurements across Europe, grassland NPP was ~750 ±150 gC m^-2 year^-1. No region-specific information was presented in that study.

In the Ciais et al. [101] study, grassland NPP was estimated from different data sources including the PaSim model. European mean NPP varies considerably between the data sources (see Table 1 in [101]) ranging between 1237 and 558 gC m^-2 year^-1. The study proposed the 750 ± 150 gC m^-2 year^-1 for NPP that was the base information for Schulze et al. [8].

Chang et al. [102] used the ORCHIDEE-GM (Organizing Carbon and Hydrology in Dynamic Ecosystems with Grassland Management) process-based ecosystem model to simulate different grassland conditions in Europe in spatially explicit manner. ORCHIDEE-GM was constructed by the widely used ORCHIDEE dynamic global vegetation model and the PaSim model that was mentioned above [103], which means that PaSim is rather relevant in the grassland related literature in Europe. Chang et al. [102] used management and fertilization information with land cover change data from the HILDA database in order to perform realistic simulations. The study found that overall NPP increased between 1991 and 2010 in Europe. The changes are caused by climate change, rising atmospheric CO₂ concentration and nitrogen addition. Based on the Chang et al. [102] study the European mean NPP was 600 gC m^-2 year^-1 around 2010.

Focusing specifically in Central Europe, Barcza et al. [104] used the Biome-BGC model to simulate (among others) NPP for Hungarian grasslands. For the 2002 and 2009 time period the mean NPP was 958 gC m^-2 year^-1 which is considerably higher than the result from the present study. Barcza et al. [104] also used the MOD17 NPP product (see Section 2.3.3.) for Hungarian grasslands, with NPP of 578 gC m^-2 year^-1 (2002-2009 mean) that is closer to the presented results, although the target time period is different.

In summary, the updated, data-driven estimate for the NPP of Hungarian grasslands is lower than the European average, and even lower than the values published earlier. It can be attributed to the dry years after 2020, but it can also indicate that the NPP of Hungarian grasslands was overestimated in the past as well. Clearly, more advanced modeling studies are needed in the future to constrain the NPP for the drought-prone grasslands of Central Europe.

4.5. Limitations of the Study

There are several uncertainties and limitations associated with the present study that have to be mentioned. Those uncertainties obviously affect all three main topics addressed here, namely the grassland area estimation, the in situ observations and the remote sensing-based upscaling. In terms of grassland mapping Section 4.1. provided a detailed overview about the uncertainties. Here we would like to point out that spatial and statistical accuracy depends significantly on the original purpose of the mapping and the spatial resolution. Input reference data, satellite data, classification, post-processing, masking, visual interpretation can all be burdened with errors [17]. The EU level initiatives might resolve at least part of the issues thus grassland are estimations can be constrained supporting the scientific community.

Considering the reference data, as mentioned above, biomass sampling was done using diverse methods. One obvious limitation of the study is the typically undocumented destructive sampling methodology. In case of Fülöpháza, Orgovány and Bugac, the measurements were consistently planned and accurately documented, and the measured quantity is in the order of gBM m^-2. In contrast, in case of data from the national parks, the calculated or estimated number of bales and the uniformly defined bale weight of 300 kg may also be a source of error. The HCSO estimates yield using ton per hectare unit and assigns completely homogeneous yields to the NUTS3 level regions, which can lead to a significantly large error when projected onto square meters.

In situ data has limited representativeness, and is inevitably affected by other factors like microtopography, soil texture, soil organic matter content, soil hydrology, occasional presence of elevated groundwater, herbivory and other undetectable processes. In certain cases, we average several measurements and assign the average to an arbitrarily defined square that is considered to be representative to the landscape. This homogenizes small-scale differences, thus reducing the accuracy of the field data. Obviously, those uncertainties hold information that can be exploited in future studies.

In Hungary, due to the widespread mowing at least once in a year minimizes the risk of leaving “old” standing dead biomass in the field [15]. The presence of litter is another question. Litter in grasslands has its own decomposition dynamics, and it is hard to provide an estimate for the age of litter. Dynamic models like Biome-BGCMuSo, PaSim, ORCHIDEE or ModVege are needed to properly quantify litter dynamics.

In most cases, the moisture content of the measured biomass is not recorded. Some of the cut biomass measurements are oven-dried, while in the case of baled data, the plant can only dry naturally (leading to air-dry moisture content that is not exactly known) during the period between mowing and baling. Data collection should strictly accompanied by detailed documentation of the observation data.

Remote sensing-based information is another source of uncertainty. In mowing data, the inaccuracy of plot boundaries and the number of mixed edge pixels can have a slight but measurable impact on the plot averages calculated from satellite proxies. A more significant source of error may be the inhomogeneity within the parcel, the size of unmowed areas or strips, the presence of individual trees or groups of bushes, and mosaics with wetlands.

Perhaps the most critical source of error is the impact of grazing on the combination of observations and satellite data, and thus on the biomass estimated from it. The Fülöpháza dataset is from a long-untouched area, undisturbed by mowing or grazing. Among the national parks, the Fertő-Hanság plots are mowed, while in Hortobágy, in most cases, grazing also takes place alongside mowing. Even if it is known whether a given plot has been grazed, the interval and intensity of grazing are usually not known, so it was not possible to focus on this in the present study.

Last, but not least, the remote sensing data time series are heavily filtered and smoothed, which effectively masks the disruptive effect of clouds, but also introduces uncertainty regarding biomass. Obviously, improvements are needed in terms of the quality of remote sensing data, and the biomass estimation methods. Uncertainty can be quantified and error propagation strategies can be applied. In our study the ensemble method was used to quantify uncertainty, but e.g., the linear models can be replaced by nonlinear equations, also considering observation data uncertainty. Additional data like meteorology and soil texture can be included to further understand the variability of the proxy-observation diagrams. Nevertheless, our study provides a solid foundation for further studies with many directions for improvements.

5. Concluding Remarks

The study provided pragmatic solutions for three major scientific gaps in our understanding of Central European, drought-prone grasslands. First of all, the area covered by grassland is highly uncertain but essential for the calculation of the total biomass production, being a huge economic value. Our results indicated large variability of total grassland area per terminology/scientific use. The study demonstrated the importance of proper grassland area estimation for practical applications. Second, there is a need to obtain country-wise, consistent biomass sampling that covers different grassland types under variable soil and climatic conditions. The presented dataset, that is the collection of diverse data sources, is unique and contributes to the understanding of grassland to climate fluctuations. Third, there is a strong need to estimate grassland production in a spatially explicit manner. This is solved in the study within an ensemble framework using remote sensing data and the collected observations. The constructed dataset is unprecedented in Central Europe and provides a solid basis for other studies as well.

Focusing specifically on carbon cycle and greenhouse gas budget, the calculated ANPP and NPP can act as a reference for upscaled biogeochemical models aiming to quantify the present-day carbon budget of Hungarian grasslands. It can be used to benchmark other dataset and provides guidelines even at 10 m spatial resolution.

The main message of the study is the importance of the grassland categorization and spatial extent estimation. Even if we can constrain the biomass yield for unit grassland area using increasingly sophisticated methods/models, the upscaling is still a sensitive issue. The national or continental scale totals can only be quantified if a consistent and consensus-based estimation exists for the grassland cover. This land use/land cover issue is of course not standalone, since some of the grasslands are converted to croplands or even affected by reforestation or afforestation. Set-aside in croplands can be considered as temporal grassland, but this might not be suitable for e.g., the CAP logic. In any case, there is still a lot to do with grassland categorization and spatially explicit productivity estimates not only in Central Europe, but globally as well.

Hungary’s current climate and geography make it a proxy for the future of many other European regions. The present study provides a baseline in two ways. First, it is a “drought frontier baseline”, as Southern and Central Europe face some aridification, the behavior of Hungarian grasslands (which already exist on the edge of water-limiting conditions) acts as a “baseline” for what other European grasslands may look like if aridification intensifies in 20-30 years. Second, it is a “pedo-climatic gradient baseline”, because the study covers many ecosystems from nutrient-rich forest soils to open sandy grasslands, thereby capturing a range of conditions found across the continent. The presented results provide valuable insight into the future of European grasslands in general.