## 2 Materials and Methods
### Key Input Data for BEPS Simulation  
Sentinel-2 Level-1C top-of-atmosphere (TOA) reflectance products are used to derive LAI. These products are organized into tiles of 100 by 100 km\({}^{2}\) in UTM projection and WGS4 ellipsoid and are available from the Copernicus Open Access Hub ([[https://scihub.copernicus.eu/dhus/#/home](https://scihub.copernicus.eu/dhus/#/home)]([https://scihub.copernicus.eu/dhus/#/home](https://scihub.copernicus.eu/dhus/#/home)) (accessed on 1 June 2019)). Two tiles (14 UNA and 14 UNV) covering the SMAPVEX16-MB region were selected and mosaicked. For BEPS simulation, a region of approximately 40 km \(\times\) 60 km is subset to match the region shown in Figure 2.  
The Sen2 Cor atmospheric correction processor [60] is used to convert the TOA reflectance to top-of-canopy (TOC) reflectance, i.e., Sentinel-2 Level-2A product. A Simple Ratio index (SR), defined as a ratio of near-infrared band to red band reflectances (or band 8A to band 4 reflectances from Sentinel-2) was used to derive LAI at 20 m resolution, while taking into account the effect of bidirectional reflectance distribution [61]. We chose band 8A in 20 m resolution rather than band 8 in 10 m resolution because band 8A is more sensitive to LAI. In [61], the LAI is modeled as  
\[L=f\left(\text{SR}*f_{\text{BEDF}}\left(\theta_{s},\theta_{v},\varphi\right)\right) \tag{6}\]where \(\theta_{\text{s}}\) is the solar zenith angle, \(\theta_{\text{v}}\) is the view zenith angle, and \(\varphi\) is the relative azimuth angle between the sun and the viewer. The bidirectional reflectance distribution function (BRDF) kernel \(f_{\text{BRDF}}\) is simulated using a physically based geometrical optical model [62]. In this way, the effect of angular changes to SR is minimized. Although we used the terminology \"LAI\" in this study, our retrieved LAI from Sentinel-2 includes the area of all green elements. Although Sentinel-2 has a high revisit frequency, the availability of cloud-free and cloud shadow-free images is limited (Table 1). A locally adjusted cubic-spline capping (LACC) method is used to smooth the time-series of LAI and discretize the smoothed LAI curve into daily values for BEPS modeling [63]. In LACC, a variable local smoothing parameter, which controls the local smoothness of the fitted curve, is automatically determined according to the local curvature of the original seasonal variation pattern. An iteration procedure is designed to produce a seasonal capping curve by progressively replacing abnormally low values with fitted values. Normally, less than 5 iterations are required to produce a smoothed capping curve. The lower-right corner in our study area is excluded from GPP simulation owing to too few Sentinel-2 images available for LACC inputs. Six meteorological variables including air temperature, relative humidity, wind speed, and atmosphere pressure at 2 m above the surface, and incoming solar shortwave flux and total precipitation at the surface level are used to drive BEPS. These variables are acquired from a global reanalysis dataset, Modern-Era Retrospective Analysis for research and Applications, Version 2 (MERRA-2) [64]. In BEPS, the soil up to a 2 m depth is stratified into five layers for the simulation of soil water and heat fluxes. Four soil attributes including texture, bulk density, organic content and reference depth are extracted from the Soil Grids ([[https://soilgrids.org](https://soilgrids.org)]([https://soilgrids.org](https://soilgrids.org)) (accessed on 15 June 2019)) that provide global soil attributes at 250 m resolution [65]. Clay and fine loamy soils dominate in the eastern and southern parts of the study area and coarse loamy and sandy soils are prominent in the western part (Figure S1). To match the pixel resolution of Sentinel-2 data, the meteorological variables are resampled to 20 m resolution using a bilinear method, and the soil attributes are resampled to 20 m using the nearest neighbor method.