The carbon budget of the managed grasslands of Great Britain constrained by earth observations

Grasslands cover around two thirds of the land area of Great Britain (GB) and are important reservoirs of terrestrial biological carbon (C). Outside a few well-monitored sites the quantification of C dynamics in managed grasslands is made complex by the spatio-temporal variability of weather conditions combined with grazing and cutting patterns. Earth observation (EO) missions produce high-resolution frequently-retrieved proxy data on the state of grassland canopies but synergies between EO data and biogeochemical modelling to estimate grassland C dynamics are under-explored. Here, we show the potential 5 of model-data fusion (MDF) to provide robust near-real time analyses of managed grasslands of GB (England, Wales and Scotland). We combine EO data and process-based modelling to estimate grassland C balance and to examine the role of management. We implement a MDF algorithm to (1) infer grassland management from vegetation reduction data (Proba-V), (2) optimise model parameters by assimilating leaf area index (LAI) data (Sentinel-2) and (3) simulate livestock grazing, grass cutting, and C allocation and loss to the atmosphere. The MDF algorithm was applied for 2017 and 2018 at 1855 fields 10 sampled from across GB. The algorithm was able to effectively assimilate the Sentinel-2 based LAI time series (overlap=80%, RMSE=1gCm−2, bias=0.35 gCm−2) and predict livestock densities per area that correspond with independent census-based data (r=0.68). The mean total removed biomass across all simulated fields was 6 (±1.8) tDM ha−1y−1. The simulated grassland ecosystems were on average C sinks in 2017 and 2018; the GB-average net ecosystem exchange (NEE) and net biome exchange (NBE) for 2017 was -232±94 and for 2018 was -120±103 gCm−2y−1. The 2018 summer drought reduced C sinks, with a 15 9-fold increase in the number fields that were C sources (NBE>0) in 2018 compared to 2017. We conclude that management in the form of sward condition and the timing, intensity and type of defoliation are key determinants of the C balance of managed grasslands. Nevertheless, extreme weather, such as prolonged droughts, can convert grassland C sinks to sources. 1 https://doi.org/10.5194/bg-2021-144 Preprint. Discussion started: 14 June 2021 c © Author(s) 2021. CC BY 4.0 License.

on livestock density and timing of cutting. Such data are available at regional-to-global scales from agricultural censuses, but not at field scale or in near real time. Spatial disaggregation of census data introduces significant temporal and spatial uncertainties to model predictions Vuichard et al. (2007); Rolinski et al. (2018). A solution to this is the use of specialised algorithms that detect the timing of grazing and cutting events at field scale from EO-based data (Yu et al., 2018;Reichstein 95 et al., 2019). However, EO-based biomass removal detection algorithms do not consider grassland ecosystem dynamics and are based on image and/or time-series analysis (Reinermann and Asam, 2020).
In this study, we combine EO data and process modelling to detect management operations and estimate resulting C dynamics over a large domain (GB) and at fine resolution (sub-field scale). We use a parsimonious process-model of grassland C dynamics (DALEC-Grass) that is integrated into a probabilistic model-data fusion (MDF) algorithm (CARDAMOM). CAR-100 DAMOM generates field-specific calibrations of DALEC-Grass by assimilating EO-based LAI time series. The algorithm has been validated against an extensive dataset of field measured C flux and biomass data from two managed grasslands in Scotland and further tested using 4 years of EO-based LAI data in 3 variably-managed fields in South West England (Myrgiotis et al., 2020(Myrgiotis et al., , 2021. Here, the MDF algorithm is implemented for 2017-2018 on a large sample of 1855 managed grassland fields in GB (England, Wales and Scotland). The fact that in 2018 GB was affected by a summer heat and drought wave (summer 2018 105 was ≈ 1 o C warmer than summer 2017) allows us to examine the impact of climate anomalies on grassland C balance (Kendon et al., 2018(Kendon et al., , 2019. The experimental approach uses land cover data (UK Land Cover Map) to sample one managed grassland field per 25km 2 across GB. Grazing intensity, cutting timing and yields and C pools and fluxes are simulated by DALEC-Grass for every field using a local model calibration. In order to assess the MDF approach, C cycle process estimates are compared to biomass utilisation from the relevant literature and livestock density data from recent agricultural census data. The aim of 110 this study is to answer four key questions: 1. Can we detect grassland vegetation management variability over large domains at field scale by assimilating EO information on leaf area index? 2. What is the C balance of managed grasslands and how does it vary across GB?
3. Which factors control the predicted C balance and biomass removals? and below-ground plant tissues based on the assumption that C allocation to roots increases after sufficient leaf area has been developed such that there are diminishing returns on further canopy expansion (Reyes et al., 2017). The model uses a simple scheme to describe C allocation to soil with two pools considered; a more labile litter pool to which dead plant material and manure-C are added, and a more recalcitrant soil organic carbon pool (SOC), which receives C from the litter pool only. The model's 31 parameters are presented in Table A1 in Appendix A . For detailed information on the model's concept, its struc-130 ture, its validation against measured CO 2 fluxes and its testing under variable grassland management using EO data we refer to (Myrgiotis et al., 2020(Myrgiotis et al., , 2021. Grassland management events are inferred from the combination of vegetation reduction time-series, which are provided to the model as inputs, and EO-based data on LAI, which are assimilated through a MDF algorithm. At each time step the algorithm reads the vegetation reduction information and decides whether to simulate the corresponding biomass removal as 135 a grass cutting, a grazing event or ignore it as being unlikely based on (1) the simulated aboveground biomass (AGB) at that time and (2) a set of biophysical and agronomic conditions.

Carbon Data Model Framework
The Carbon Data Model Framework (CARDAMOM) is a Bayesian MDF framework that is tailored for use in ecosystem biogeochemistry studies (Bloom et al., 2016). The algorithm is provided with prior information on the distribution of each 140 model parameter (i.e, min, max, shape). By assimilating observational data (LAI, soil C stocks etc) the algorithm updates the distribution of model parameters following the rules of Bayesian inference. CARDAMOM has been used in studies of ecosystem C fluxes in the past, including two recent studies using DALEC-Grass (Myrgiotis et al., 2020(Myrgiotis et al., , 2021. A key aspect of CARDAMOM is the use of ecological and dynamic constraints (EDC), which are conditions applied to the parameter sampling process in order to ensure the mathematical, ecological and biogeochemical sensibility (or "common sense") of the simulated 145 system. In simple terms, CARDAMOM examines whether the simulated pools and fluxes that result from implementing the model with a sampled parameter vector behave in realistic ways; i.e. do not exceed certain user-defined, widely-accepted and literature-based limits. The EDCs used in CARDAMOM for this study are presented in Table A2 in Appendix A .
Bayesian inference is performed in CARDAMOM using the root mean square error (RMSE) between the simulated and the EO-based LAI time-series to calculate and attribute likelihoods to every sampled parameter vector. The Metropolis-Hastings 150 (MH) method has been used in two previous CARDAMOM / DALEC-Grass studies. However, the duration of implementing the CARDAMOM MDF algorithm for a single field is a crucial aspect for a large-domain study such as the present one. The MH method is time-consuming and so the Simulated Annealing (SA) algorithm was used here (Kan et al., 2016). After testing the SA algorithm to identify the optimal number of repetitions for achieving acceptable chain convergence the number of parameter proposals was set equal to 5000000 per chain. The uncertainty around the assimilated (Sentinel-2) LAI data was set to 15% (relative standard deviation) in this study. However, it should be noted that the uncertainty around remote sensing-based LAI data is poorly determined but largely underestimated (Zhao et al., 2020).
A uniform distribution was used for each of the 31 DALEC-Grass parameters with the range for each parameter presented in Table A1 (Appendix A). In Myrgiotis et al. (2020) DALEC-Grass parameter priors were refined through implementing the model using known vegetation management (cutting dates, livestock density time-series) and by assimilating field-measured 160 LAI and CO 2 flux data. In Myrgiotis et al. (2021) these model priors have been further refined using EO-based vegetation anomaly time-series (for vegetation management inference) and by assimilating Sentinel-2 based LAI time-series. The limits of the uniform parameter distributions used in the present study are based on the results of Myrgiotis et al. (2021) but 4 parameters were allowed to vary more than these suggested in order to better consider the variability in management factor across GB grasslands (indicated with * in Table A1 in Appendix A ).

Location of managed grasslands
For the identification of the location and limits of managed grasslands we used the 2018 Land Cover Map plus (LCM), which is produced and published annually by the Centre of Ecology and Hydrology (CEH) of the UK (www.ceh.ac.uk/crops2015).
The LCM includes geo-referenced polygons of improved grassland fields in GB that are identified as such using a combination of reflectance data. The LCM data are validated against ground observations of land-use type.

Earth observation data
DALEC-Grass uses information on vegetation canopy reduction to infer vegetation-related management operations (i.e. grazing or cutting). To produce this time series for each simulated field we obtained LAI data from the Copernicus Global Land Service (CGLS) database for GB for 2017 and 2018. The CGLS LAI data comprise top-of-atmosphere reflectance processed products from the Proba-V satellite and have a spatial resolution of 300m and a temporal resolution of 10 days (?). Each data point 175 has an uncertainty attributed to it. In periods of high cloud coverage the uncertainty over the LAI estimates increases as the provided LAI value was estimated using a machine learning model built with time series for past years and neighbouring pixels. For each simulated field the corresponding time series has been converted, from their original 10-day time-step, to a weekly time step using linear interpolation. Thereafter, the reduction in between subsequent dates in the time series has been calculated. When the change between week n and n + 1 was positive the reduction value for week n is 0. Hereafter, we refer to 180 this time series as the "vegetation reduction" time series. The EO-based LAI data that are assimilated in CARDAMOM were extracted from Sentinel-2 images. Atmospherically-corrected images at 20m resolution (L2A product) were downloaded from the Amazon Web Services (AWS) Sentinel-2 data pool. The images were processed to remove pixels with cloud and haze and, then, used to calculate LAI (at 20m) using the sen2cor algorithm (Weiss and Baret, 2016).
Values of soil C (gCm −2 at 60cm depth) at 300m resolution were obtained from the most recent version of the SoilGrids 190 database (Hengl et al., 2014). For every simulated field the mean and standard deviation (SD) of the corresponding SoilGrids pixels are used to define the range of the model's initial SOC pool size parameter.
Agricultural census-based data on the number of sheep and cattle (beef and dairy) were obtained from the EDINA AgCensus database in 2020 (AgCensus, 2020). They are used to independently evaluate the estimates from the inversion of EO data.
The AgCensus data are produced by spatially disaggregating the numbers of cattle and sheep recorded at the level of local 195 administrative units (for each UK country), into a 5km grid of the UK. The most recently available livestock data for each constituent UK country refer to different years; 2010 for England, 2015 for Wales and 2017 for Scotland.

Sampling of grassland fields
Implementing the MDF algorithm for the thousands of fields that are classified as improved grassland in the LCM database is 200 computationally demanding and time consuming. In addition, the resolution of the vegetation reduction data, which are used in DALEC-Grass to infer grazing and cutting, is 300m (9ha). This means that vegetation reduction data for fields smaller than 9ha are increasingly noisy. For this reason, and taking into account that the average managed grassland field is 5-9ha in size, we set a minimum limit of 6ha (and a maximum of 13ha) when filtering the LCM dataset to obtain the location of fields. Moreover, the number of EO data points available for each field depends on the time of image capturing and the amount of cloud cover 205 at overpass. As a consequence, the number of dates of available EO data can vary considerably between polygons. We set a limit of having at least 30 Sentinel-2 data points (during 2017-2018) for a sampled field to be simulated. The fields that met the above-mentioned conditions were allocated to 25km 2 cells of a 5km grid of GB. One field was randomly selected from each cell, which resulted in a set of 2108 fields (Fig, ??). CARDAMOM was implemented for each of the selected fields for 2017 and 2018 by running DALEC-Grass at a weekly-time step while assimilating the corresponding EO-based data when available.

Assessment and analysis of MDF predictions
To assess the effectiveness of the LAI assimilation process we quantify the level of fit between MDF-predicted and EO-based time-series using (1) the % of overlap between EO-based data points (mean) and corresponding MDF-predicted ranges (95% confidence interval); and (2) the RMSE and (3) the bias between of simulated and observed data. To account for the possibility that some of the simulated fields may not be managed grasslands but classified as such in the LCM data, we remove from the 215 results any fields for which the estimated overlap is < 50% (see results for size of post-MDF dataset).
To answer our first question, the MDF-predicted weekly grazed biomass is converted into livestock units (LSU) per ha following the assumptions that : (1) 1 cattle is 1 LSU and one sheep is 0.11 LSU; (2) 1 LSU weighs 650kg; (3) an animal demands ≈ 2.5% of its weight in the form grass dry matter (DM) when grazing; and (4) 47.5% of DM is C (Vertès et al., 2018).
The MDF-predicted and independent census-based LSU per ha are compared using the correlation coefficient (r) and the root 220 mean square error (RMSE) as the assessment measures. The MDF-based estimates of grass biomass utilisation across GB are assessed against data from the (Qi et al., 2017) study.
To answer our second question we present and examine the annual and seasonal C balance and the cumulative annual fluxes of the simulated fields. To assess what controls the predicted C balance of the simulated grasslands (our third question) we quantify the correlation coefficient between meteorological model drivers, management-related model parameters and MDF 225 predictions of C cycling. In order to provide a more quantitative assessment of the factors that control grassland C dynamics we quantify the relative impact of management and climate on the MDF-predicted NBE. We use the model meteorological drivers and the posterior model parameters related to management and climatic controls for every simulated field to train a random forest (RF) model that estimates NBE. 75% of the data are used to train the RF model and 25% to assess its predictive ability (coefficient of determination). Thereafter, we use the Shapley Additive Explanations (SHAP) method to quantify how much 230 each RF-predictor affects the RF-predicted NBE (?). The SHAP method examines the structure of the RF model and provides the weight (SHAP value) that the model gave to each predictor. SHAP values can be seen as the machine learning equivalent of the coefficient of determination (r 2 ). The estimated SHAP values are normalised (0-1) to be comparable to r 2 . We note that RF is used in this study solely to support MDF data analysis and not for predictive purposes.
Finally, for each simulated field and model output the MDF algorithm produces a mean and 95% confidence interval. To 235 answer the fourth question of the study, we quantify the predictive uncertainty around an output by calculating its relative confidence range (RCR). RCR is equal to the size of the MDF-predicted 95% confidence interval divided by the corresponding mean, and multiplied by 100 (expressed as %). We present and examine the estimated RCRs to identify the key factors that affect them.

Assimilation of EO-based LAI data
For 12% of the initial dataset of 2108 simulated fields, our analysis failed to generate a simulated-vs-observed LAI overlap > 50% . These field were removed from the analysis, so the final reported dataset includes 1855 fields. Based on these field, three performance metrics indicated that CARDAMOM effectively assimilated the provided EO-based LAI time-series (   The analysis suggests that the 1855 simulated fields were managed with varying intensity. The majority of simulated fields were grazed-only (75%) and no cut-only fields were simulated. Grazed biomass exceeded cut biomass in 85% of the fields 260 (GCD > 0) and cut biomass exceeded grazed biomass in the remaining 15% (GCD < 0). The mean MDF-predicted annual a yield of 9.76±2.03 tDMha −1 y −1 (Qi et al., 2017(Qi et al., , 2018. Most of the MDF-predicted first grass cuts (85%) occurred between the first half of May and the second half of July. For the fields where more than one cut was identified the period between the first and last cut was ≈ 2 months. The MDF-predicted day-of-year of first cut increased northwards with the mean date of first cut in northern England and Scotland being 3-6 weeks later than in the southern half of GB. This spatial pattern is likely the combined effect of differences in the onset and duration of 270 the grass growing season and in related management decisions. Due to the small share of cut-and-grazed fields in the simulated dataset, and in order to make the spatial pattern more visible, we present the average month of first cut on a regional basis ( Fig.   A1 in Appendix A ).

Predicted C balance and dynamics
MDF-based C cycle estimates show that management affected the C balance of the simulated grassland ecosystems signifi-275 cantly. The difference between grazed and cut biomass volume (GCD) is used to present the impact that these two biomass removal methods have on C balance. The mean annual GPP across GB fields was 30% higher (1992±400 gCm −2 y −1 ) for fields where most biomass was removed via grazing (GCD>0) compared to those where cut biomass exceeded grazed biomass (GCD<0) (1518±426 gCm −2 y −1 ) ( Figure 5). REco was higher for fields dominated ny grazing also. The mean NEE across GB was -232±94 gCm −2 y −1 , the relative role of grazing compared to cutting did not have affect NEE significantly, and 95% of 280 the simulated fields were net C sinks at the ecosystem scale. When considering the role of cutting and grazing C removals and returns to the ecosystem, the impact of cutting as a biomass removal method becomes significant. The NBE of fields dominated by cutting was 38±70 gCm −2 y −1 while fields with more grazing had a NBE of -126±95 gCm −2 y −1 . On average, 60% more C was removed (grazed and cut) in mostly-cut (GCD>0) grasslands than in mostly-grazed (GCD<0) grasslands. The flux of C into the SOC pool was, on average, 66% larger in mostly-grazed than mostly-cut fields. The annual soil C sequestration rate 285 (annual change in SOC size) for mostly-cut grasslands was 116±52 gCm −2 y −1 and 36±40 gCm −2 y −1 for grazer dominated fields. The spatial distribution of MDF-predicted GPP, REco, NEE, NBE, removed biomass and C flux into SOC is presented in the cartograms of Figure A2 in Appendix A .
Seasonal NEE varied across GB, with strongest sinks in Spring and Summer, strongest sources in Autumn, and close to neutral met exchanges during Winter (Fig. 6). However, there were clear inter-annual differences between 2017 and 2018 in 290 the analysis. Across the southern third of GB (the Midlands and Southern England) many grasslands became C sources during the summer of 2018 while remaining areas were weaker sinks than in 2017 6). This pattern was driven by the 2018 European drought and heat wave, which affected GB as a whole and was particularly acute in the southern half of England (Sibley, 2019). The three-week rolling average VPD (DALEC-Grass met driver) in summer 2018 across GB was 50% higher than in summer 2017 (Fig. A4 in supplementary material). The GB-average GPP, REco, Hr, C flux to soil and soil C sequestration

Controls on C cycling in UK grasslands
Correlation coefficients generated across the 1855 fields show the links between meteorological drivers, key processes (model parameters) and model outputs (C exchanges). There are strong positive correlations between GPP, Reco, C inputs to soil 305 and root:shoot ratio, and these factors are strongly negatively correlated with NBE and NEE. The most productive fields (higher GPP) are associated with high inputs of C to soils and are the strongest C sinks (more negative NEE and NBE).
Among modelled processes, the ratio of C allocation to roots relative to stems and leaves (root/shoot ratio) is the most strongly correlated with the net C balance of the simulated fields (Fig. 7). More-frequent and higher-yielding grass cuts reduce root-toshoot ratio and, therefore, reduce the flux of C to litter and, subsequently, to the recalcitrant soil C pool (SOC  We expanded on the correlations-based analysis of the MDF results by using (1) the MDF-predicted data on NBE, and (2) the corresponding meteorological drivers and (3) model parameters describing climatic and management controls on grass 320 growth, to train a RF model that estimates NBE. The resulting RF model was able to explain 93% of the variance in MDFpredicted NBE (r 2 =0.93) using 5 predictors: VPD, mean T, solar radiation, all posterior DALEC-Grass parameters related to climatic effects on grass growth and all posterior parameters related to grassland management. The weight (normalised SHAP) attributed by the RF model to these 5 predictors suggests that management parameters (aggregated weight for all parameters in the group) were the most important factor for grassland NBE over the simulated period (Table 2. The normalised SHAP for 325 management parameters was the highest among the 5 NBE predictors in 2017 (contributed by 34% to NBE) and the second highest (contributed by 38%) in 2018. The 2018 summer heat wave caused the contribution of VPD to NBE to increase from 3% in 2017 to 40% in 2018. Overall, these results reaffirm the conclusions of our correlations-based analysis and clarify the importance of grassland management relative to climate and climatic anomalies.  The size of the uncertainty around MDF estimates is quantified using the RCR (relative confidence range) of MDF outputs.
The mean RCR is 42±9% for LAI, 21±10% for GPP, 18±6% for REco and 26±16% for grazed biomass (Figure 8). MDF predictive uncertainty is therefore a small faction of the mean estimate of these scalar variables. The GB-average RCR for LAI, GPP and grazed biomass prediction increased from 44%, 26% and 27% in fields where cut biomass did not exceed grazed biomass (GCD>0) to 54%, 40% and 52% in fields where cut biomass exceeded grazed biomass (GCD<0). The higher 335 RCR (mean and SD) for LAI and grazed biomass is caused by the spatio-temporal uncertainty in the vegetation reduction time series. This input-related uncertainty leads to the MDF algorithm being less effective in identifying cutting instances in some simulated fields; i.e. sampled parameter vectors produce varying predictions on the timing and intensity of grass cuts.
The impact of input data uncertainty on RCR is also reflected on the shape of RCR distributions for GPP and grazed biomass (violin plots in Fig. 8). The assimilated EO LAI data condition the simulated LAI, which combined with the simulated removals 340 (grazing/cutting) determine the weekly GPP at each simulated field. Reducing the uncertainty (spatial and temporal) around the vegetation reduction data will lead to less variable predictions of cutting timing and intensity and, thus, lower predictive uncertainty overall.  to the field-size limits (6-13ha) used in sampling for fields across GB, the share of less intensively managed grasslands is likely biased high in the simulated fields' dataset (Qi et al., 2018).

355
Cut-only grasslands were under-represented in our analysis. We argue that this is an artefact of (1) the noise in the vegetation reduction time-series, which led to cut-only grasslands failing to pass the 50% overlap limit and being excluded from the final dataset; and (2) the inclusion of fields 6-13ha in size in the set of simulated fields. The inclusion of fields that are 6-13ha in size could could have led to an under-representation of cut-only fields but we note that no data exist on the percentage of managed grasslands that are grazed-only, cut-and-grazed and cut-only.

The C balance of managed grasslands
The majority of managed grasslands in GB were net C sinks during 2017 (NEE=-289±76 gCm −2 y −1 ) and 2018 (NEE=-174±74 gCm −2 y −1 ) at the ecosystem level, i.e. based on CO2 gas exchanges. Numerous studies have concluded that managed temperate grasslands are, on average, C sinks but NEE estimates vary greatly between ≈ -700 gCm −2 y −1 to almost C-neutrality Skiba et al., 2013;Ciais et al., 2010;Rutledge et al., 2015;Ammann et al., 2020;Soussana et al., 2007;Ammann 365 et al., 2007;Gilmanov et al., 2007). The scale of NEE increase between 2017 and 2018 is comparable to past field-based estimates under normal and heat-wave conditions ( et al., 2011). When considering C fluxes that also included grazing/cutting removals and manure return to the soil (i.e. NBE), the simulated grasslands were net C sinks during 2017 (NBE=-191±82 gCm −2 y −1 ) and close to C neutral in 2018 (NBE=-49±70 gCm −2 y −1 ). These NBE estimates are comparable to those in the literature but the inconsistency in the variables included in NBE calculations makes comparisons less than straightforward 370 Skinner, 2008). Based on RF-based analyses of climate drivers and model parameters we argue that the increase in mean annual NEE and NBE between 2017 and 2018 was caused by the 2018 summer heatwave. The negative effect of elevated annual temperatures on grassland biomass productivity and NEE has been examined in measurements and model-based studies on European grasslands before (? Ciais et al., 2005;Jansen-Willems et al., 2016;Thompson et al., 2020).
The mechanistic understanding, and the model representation, of how plants respond to increased VPD is improving but key 375 aspects are still disputed (Grossiord et al., 2020;Massmann et al., 2019).
Our study shows that biomass removals were key for the C balance of managed grasslands. The role of cutting relative to grazing as a biomass removal method was found to be particularly important. Our MDF-predicted dataset does not allow us to compare grazed-only with cut-only grassland C balances because no cut-only fields were simulated. Grasslands in which most biomass was removed via cutting had a lower GPP and REco as opposed to grasslands in which grazing was the main 380 biomass removal method (Fig. 5). However, when GPP and REco are summed the resulting NEE did not vary significantly between mostly grazed and mostly cut grasslands. All of the simulated grasslands were grazed and underwent more or less frequent defoliation during the simulated period. When cutting occurs the leaf area of a grassland is reduced close to zero, which represents a diminution of the grassland's photosynthetic capacity. According to DALEC-Grass, in the post-cutting period the simulated grassland allocates almost all of its C to aboveground tissues (stems,leaves) in order to build up the leaf 385 area necessary to increase photosynthetic activity and sustain growth. This causes a smaller root-to-shoot ratio in grazed and cut grasslands compare to grazed-only ones as well as a reduction of root-based C inputs to the litter pool (lower Rh). Grazing that occurs during the post-cut period maintains the volume of leaves at relatively low levels. This leads to a reduced annual GPP for grasslands that are both grazed and cut and also means lower manure-C returns to the soil, which explains the weaker C sinks (higher NBE) of most grazed-and-cut grasslands compared to grazed-only fields. These findings are supported by the relevant 390 literature. For example, Skinner (2008) found that higher biomass removals increase NBE based on C flux measurements in cut-and-grazed temperate grasslands in the USA. While Koncz et al. (2017) used eddy covariance measurements of C fluxes at a cut-only and a grazed-only field in Hungary and found that the cut field had a more positive NBE than the grazed field. Senapati et al. (2014) reached the same conclusion as Koncz et al. (2017) using eddy covariance measurements from a cut-only versus grazed-only experiment in France. Soussana et al. (2010) reviewed studies on European managed grassland C balance 395 and found that grazed-only grasslands had the lowest NBE, followed by cut-only grasslands with cut-and-grazed grasslands having the highest NBE (NBE in this study included animal methane-C and C-leaching fluxes).

Factors that control managed grassland C dynamics
Management and climate have a combined effect on C dynamics and disentangling their individual impacts is a challenge (Ammann et al., 2020). Here, we used a correlation matrix of model drivers, parameters and outputs to understand how climate The conclusions that we draw in regards to which factors have more influence on grassland C dynamics are based on two assumptions. Firstly, we assume that the simulated grasslands are well-optimised for the intended use; i.e. to sustain different 410 types of livestock (e.g. dairy and/or beef cattle and/or sheep). This means that each sward is maintained in good condition and that farmers manage their fields optimally based on their long-term experience. Secondly, the fact that a large share of the simulated fields (especially in the southern half of England) experienced continuous weeks of unusually hot and dry weather conditions during one (2018) of the two simulated years is treated as a climate anomaly; i.e. climate in 2018 is not representative of normal climatic controls on C balance. Based on these assumptions, we argue that the MDF-inferred vegetation management 415 at each simulated field was adapted to weather conditions. Therefore, significant changes in ecosystem C cycling were beyond the control of human management and can, thus, be mostly attributed to weather conditions. We conclude that management is more important than climate in terms of the C balance of managed grasslands in GB. However, we note that climatic anomalies, such as heat waves and droughts, can reduce the relative importance of management as a determinant of grassland C balance.
In simple terms, human decisions can adjust grassland sink or source strength and this depends mostly on the soil's existing

Predictive uncertainty
The uncertainty of MDF predictions was overall relatively small (Fig. 8). Predictive uncertainty was noticeably higher for 425 fields that were both cut and grazed. This was caused by the fact that the MDF algorithm does not infer cutting simply by translating large reductions in vegetation as cuts. The algorithm examines each weekly vegetation reduction to decide whether to simulate it as cutting or grazing depending on the simulated amount of foliar biomass at the time, which is itself controlled by the assimilated EO-based LAI time-series. A weekly vegetation reduction will be simulated as a cut when it is reasonable both biophysically (i.e. simulated and observed LAI data before and after cut have good fit) and agronomically (i.e. it happens 430 during cutting season and will yield > 1.5 tDMha −1 ). This higher predictive uncertainty when cutting occurs suggests that the best way to obtain more accurate predictions is to improve the quality (spatial and temporal resolution) of the vegetation reduction time-series.

Limitations
This study uses a MDF algorithm that depends on EO data and process modelling of C dynamics in grasslands. The Proba-435 V-based vegetation reduction time-series that are used to drive DALEC-Grass have a resolution (9ha) that is coarse when compared to the average size of grassland fields in GB. These noisy data on vegetation reduction cause increased uncertainty in MDF predictions especially in regards to the timing of cutting events. Moreover, most areas of GB are affected by frequent cloudiness, which means that the number of Sentinel 2-based LAI data points per year and simulated field is limited compared to other parts of the world, although we ensured 30 images per field over two years in our selection process.

440
DALEC-Grass is a model that was developed and tested under UK conditions showing a very good skill in predicting C allocation and CO 2 fluxes under variable management and different soil conditions. Its parameters were refined, firstly, using field measured LAI and, subsequently, using EO-based LAI data. DALEC-Grass in its present form does not consider the effects of soil moisture and nitrogen cycling on grass growth. We believe that the MDF algorithm can be applied at locations were soil moisture and nitrogen are not limiting factors for grass growth. It should be noted that DALEC-Grass has been 445 validated against data from grasslands dominated (>90%) by perennial ryegrass (Lolium perenne) and its ability to simulate swards with larger shares of herbs and forbs has not been tested.

Future work
Our overarching aim is to produce a computational ecosystem science framework that is (1) able to utilise the swathes of EO data that are increasingly becoming available while (2) being easy to adapt and incorporate new knowledge gained from 450 field/lab experiments and ground observations. This study showed that the MDF algorithm will benefit most from (1) improving the temporal resolution and quality of EO LAI data used, and from (2) introducing processes relevant to soil moisture and N cycling aspects. We believe that by advancing on the first of these two fronts the algorithm will be able to produce more precise estimates across the UK and regions with similar agro-climatic conditions in Europe and beyond. Introducing soil moisture and N cycling-related processes to DALEC-Grass will pave the way for more detailed consideration of the effects of fertiliser 455 use and different grass mixtures, and for its application at climatically-critical rangelands and pastures across the world (e.g. tropical and dry regions). DALEC-Grass has a structure that facilitates the incorporation of modelling advances made with other DALEC-based models such those presented in Revill et al. 2021 for foliar N and Smallman and Williams 2019 for soil moisture. We also note that the quality of soil C-related data is critical for better constraining below-ground C pools and fluxes (e.g. heterotrophic respiration). The production of relevant spatial data in the future will improve the credibility for MDF 460 estimates.

Conclusions
This study presented how, by fusing earth observation data and biogeochemical modelling of managed grassland C dynamics, a MDF framework can detect biomass removals and use them to predict grassland C fluxes and balance. We showed that MDFpredicted annual yields and livestock density mirror ground based information well. In agreement with a range of studies on 465 temperate grasslands in Europe and beyond, our study reaffirms the C sink potential of managed grasslands in GB. In contrast to previous measurements and model-based studies, however, we showed how MDF can quantify and interpret C dynamics across a large domain (GB) while also resolving sub-field scale variability in vegetation management. It is widely accepted that climate change is manifesting itself, among other ways, as more frequent droughts in northern Europe (Peters et al., 2020). Our study showed how the most prolonged drought (2018) that has been recorded since 2000 affected the C balance of managed 470 grasslands. It highlights that the ability of temperate maritime grasslands to sequester C could be significantly affected by prolonged heat waves and drought. Various climatic and management-related factors affect both the annual C balance and the seasonal grassland biomass utilisation in livestock farming in GB; and northwest Europe in general. National targets for C neutrality in the agricultural sector and climate change create a challenging future for UK grassland farming. MDF represents a robust method for monitoring the biophysical state and C dynamics of any grassland field, over any domain in near-real time. 475 We argue that farmers and governments alike can use MDF approaches like this to provide needed monitoring tools for C balance, and guidance on adaptation and mitigation of climate change effects on agriculture towards meeting net-zero goals.
Appendix A * parameters for which the prior was wider than suggested by (Myrgiotis et al., 2021).