Numerical Evaluation of JULES Surface Tiling Scheme with High-Resolution Atmospheric Forcing and Land Cover Data

Land surface heterogeneity exists at all spatial scales and has many important effects on energy, momentum and mass exchange between land and atmosphere. Land surface models (LSMs) partially consider surface subgrid heterogeneity (SSGH) effects through surface tiling methods. In this study, a series of numerical experiments were conducted to evaluate the performance of the Joint UK Land Environment Simulator (JULES) LSM’s surface tiling scheme by combining atmospheric forcing data and landcover fraction data at different horizontal resolutions. Our tests quantitatively show that the surface tiling scheme can have a significant impact on model simulated fields, but cannot reflect SSGH effects adequately. Soil is considered to be homogeneous across a single grid cell in the current surface tiling scheme, which results in soil-moisture dependent anomalies in simulated carbon flux. Our numerical results indicate that the current JULES LSM needs additional updates to consider subgrid scale heterogeneities of subsurface processes and soil-vegetation interactions. (Citation: Park, J., H.-S. Kim, S.-J. Lee, and T. Ha, 2018: Numerical evaluation of JULES surface tiling scheme with highresolution atmospheric forcing and land cover data. SOLA, 14, 19−24, doi:10.2151/sola.2018-004.)


Introduction
The mass, energy and momentum fluxes between land surface and the atmosphere are primarily controlled by the land surface characteristics including the topography, the land cover and use, and the types of vegetation present. These characteristics exhibit distinct geophysical and physiological properties and can vary nonlinearly regardless of the spatial scale. The wide degree of land surface heterogeneity and the dependences of fluxes on surface characteristics result in spatially nonlinear and heterogeneous fluxes. Thus, a critical issue in land surface modeling is the reliable simulation of the effects of land surface heterogeneity, including its interconnection with biogeochemical exchanges between land surfaces and the atmosphere (Bonan et al. 1993; de Vrese and Hagemann 2016;Essery et al. 2003;Giorgi et al. 2003;Shrestha et al. 2016;Tian and Xie 2008).
The spatial resolutions of current atmospheric models range from a few tens of kilometers (for regional scale models) to a few hundred kilometers (for global scale models). Because of the rising demand for hyper-resolution modeling (Beven and Cloke 2012;Bierkens et al. 2015;Wood et al. 2011) and computational development, researchers are attempting to increase these models' horizontal resolutions by about an order of magnitude. However, doing so requires extensive computational resources, and several obstacles remain to be overcome, including achieving computational stability and physical consistency. Consequently, significant land surface heterogeneity is currently inevitable within each grid cell of a coarse-resolution climate model, and land surface models (LSMs) have to consider the effects of surface subgrid heterogeneity (SSGH).
Two main methods have been used to simulate SSGH in LSMs. The first method is surface tiling, also known as the mosaic method, which utilizes multiple fine-resolution tiles with different characteristics within a single grid cell, calculates each tile's flux individually, and averages the areal-fraction weighted flux (Koster and Suarez 1992). The surface fluxes originating from each tile are assumed to have blended completely horizontally. The second method is the statistical-dynamical method, which either parameterizes how the subgrid variability in the model state variables affects the grid-averaged fluxes or develops new parameters to aggregate the effects of fine-scale SSGH on larger-scale fluxes (Giorgi and Avissar 1997).
Significant differences exist among those methods. The tiling approach (or flux-aggregation method) is generally reported to improve the model performance more than the parameter aggregation method (Ament and Simmer 2006;Heinemann and Kerschgens, 2005). However, the tiling method's assumption of complete horizontal atmospheric mixing is valid only under specific environmental conditions. Many studies have reported that this assumption becomes less valid with increasing SSGH (Albertson and Parlange 1999;Ma et al. 2008). Thus, the objectives of this paper are (1) to test the LSM surface-tiling scheme on a local scale, (2) to quantify the effects of SSGH on errors in flux estimation, and (3) to identify the causes of those errors by combining the atmospheric forcing data and fractional land cover data at various spatial resolutions.

JULES LSM and its surface tiling scheme
The Joint UK Land Environment Simulator (JULES) Clark et al. 2011) is a community LSM that was designed to be coupled with other weather and climate models. The number of offline applications of the model, such as drought impact simulation, is increasing (Van den Hoof and Lambert 2016). In the original version of JULES, each grid cell was composed of up to nine surface tiles (including five vegetated types and four non-vegetated types). The vegetated surface type included different plant functional types (PFTs) with different structural and physiological characteristics, and all of the surface tiles within the same grid cell shared the same meteorological forcing conditions and underlying soil types.
A recent version of JULES (vn4.8) is adopted for testing and evaluation. JULES calculates the moisture, carbon, momentum,

Numerical Evaluation of JULES Surface Tiling Scheme with High-Resolution Atmospheric Forcing and Land Cover Data
resolutions and subgrid heterogeneities): (1) reference run (REFR run): higher resolution land surface data and higher resolution forcing data; (2) lower resolution forcing run (LOWR run): higher resolution land surface data and lower resolution forcing data; (3) SSGH run: lower resolution land surface data with consideration of SSGH and lower resolution forcing data; and (4) NOSSGH run: lower resolution land surface data with no consideration of SSGH and lower resolution forcing data. The initial soil wetness condition for each run of the model is derived after a sufficient number of spin-ups with the same atmospheric forcing data period have been performed over a maximum of 100 iterations, in order to allow soil moisture to reach equilibrium. This results in differences in initial soil wetness condition for each model run. In contrast to soil moisture, all of the other initial states are spatially homogeneous. The effects of meteorological forcing data resolution are tested by comparing the LOWR run and the REFR run, and SSGH effects are tested by comparing the SSGH and NOSSGH runs with the reference run. Since observational evaluation of the REFR run is currently not feasible due to a spatiotemporal lack of measurement data, the REFR run was assumed to be closer to actual conditions than lower-resolution model runs.
The estimated carbon (gross primary productivity (GPP)) and energy fluxes (latent heat flux (LE) and sensible heat flux (H )) of each model run are averaged both spatially and temporally. All of the anomalies present in LOWR, SSGH, and NOSSGH runs are calculated based on differences with the REFR run. To analyze the effects of the SSGH scheme, the anomalies in the SSGH and NOSSGH runs are compared. To determine the causes of the anomalies in the SSGH, the anomaly in soil water content is compared with the anomaly in GPP.

Performance of JULES surface tiling scheme
The effects of the meteorological forcing data resolution can be ignored, because the differences in flux estimation between the LOWR and REFR run are less than 1% for all fluxes. The JULES surface tiling scheme improves the major surface flux estimation when compared with the higher resolution model run (Fig. 3). Overall, the magnitude of the error reduction that results from considering SSGH varies by flux variable. The carbon flux (GPP) is affected more than the energy fluxes (LE and H ) when SSGH is considered (Figs. 3a, 3b, and 3c). The temporal variations in the GPP estimation errors exhibit distinct patterns related to the surface tiling scheme (Fig. 3d). When SSGH is not considered, the anomaly demonstrates a clear seasonal pattern: it reaches its highest value at the beginning of the growing season, decreases and energy fluxes for each tile under given conditions of incoming short-and long-wave radiation, air temperature, specific humidity, wind speed, precipitation, and air pressure. Both the leaf-level net photosynthesis (A) and the tree-water uptake are assumed to be constrained by the soil water stress factor ( β) as in the following equation : where θ, θ c , θ w (all with units of m 3 m −3 ), and A P (mol CO 2 m −2 s −1 ) are the mean soil moisture concentration in the root zone, the critical and wilting point concentrations, and the potential photosynthesis rate, respectively. The leaf-level photosynthesis rate A is scaled up to the canopy-level rate using the leaf area index (LAI) and the LAI-related absorbed photosynthetically active radiation. In this study, the LAI is estimated by the phenology module. The tree-water uptake, or transpiration, is estimated by considering both available soil water and the fraction of the roots reaching each layer of soil (i.e., at depths of 0.1 m, 0.25 m, 0.65 m, and 2.0 m). The rooting depths and root distribution patterns differ for each PFT. While the surface fluxes are calculated independently for each tile, both the soil temperature and moisture are assumed to be homogeneous across each grid cell, meaning that all of the PFTs within a single cell share the same soil water availability. The PFTs and their related parameters can be modified for research purposes. For example, both increasing the number of PFTs and using trait-based parameterization can improve the model's performance on a global scale (Harper et al. 2016). In addition, the optimal PFT parameters can be estimated using data assimilation techniques (Raoult et al. 2016). In this study, we follow both of these approaches (Harper et al. 2016;Raoult et al. 2016) and increase the number of tree PFTs from two (broadleaf (BL) and needleleaf (NL)) to five (three BL PFTs: Quercus spp., Betulaceae spp., and other broadleaf species; and two NL PFTs: Pinus spp. and other needleleaf species) for more realistic local land-cover representation. Parameter values for quantum use efficiency, nitrogen concentration at canopy top, scaling factor relating maximum carboxylation and leaf nitrogen concentration, canopy structure, maximum ratio of internal to external carbon dioxide concentration, and dark respiration scaling factor are derived using an optimization procedure based on eddy-covariance measured flux data assimilation for two PFTs, and others are obtained from the TRY database (Kattge et al. 2011). It should be noted that these modifications are under verification by comparison with eddy-flux data and watershed-scale water yield data.

Design of numerical experiments
The hyper-resolution meteorological dataset (90 m × 90 m and 270 m × 270 m) are derived using the WRF multi-nesting technique (Skamarock et al. 2008;Lee et al. 2016) without statistical downscaling for an entire growing season (March to October 2014) with an hourly time resolution over forested catchments (Fig. 1). The time step of JULES is one hour as forcing meteorological data, and the grid size of JULES ancillary data is set to be identical to that of atmospheric forcing data. The types of land cover are classified using the forest cover map provided by the Korean Forest Service (Kim 2017), as most of the target area is located on forested land. The map provides vectored forest area that is classified by the dominant overstory tree species. Each 90 m × 90 m grid cell is assigned a single dominant PFT. For the 270 m × 270 m grid cells, the land cover types are assigned based on either a single dominant land cover type (no SSGH) or the fractions of the PFTs in the nine subgrid cells (SSGH).
To test the JULES surface tiling scheme, four sets of simulations are conducted (Fig. 2). Each simulation combines one of the three meteorological driving datasets with one of the three kinds of land surface fraction data (which each have different spatial   until August, and then increases again until the end of the growing season. On the other hand, the anomalies in energy flux exhibit weaker seasonal trends regardless of whether SSGH is considered or not (Figs. 3e and 3f). These numerical results support those of previous studies that reported significant improvements in model performance by incorporating the surface tiling, or mosaic, approach (Essery et al. 2003;Giorgi et al. 2003;Shrestha et al. 2016). The seasonality of the errors in GPP estimation might be caused by differences in either the initial soil moisture or the phenological response. When SSGH is not considered, small NL patches are disregarded in the BL-dominated grid cells and are not represented in the land surface (L2 and L3 in Fig. 2). While BLs show clear phenological changes, NLs tend to maintain substantial leaf areas during the dormant season, which causes differences in LAI mostly in leaf onset period when SSGH is considered (Fig. 4a). This results in seasonal carbon uptake differences between the BL and NL PFTs. BLs have reduced carbon uptakes during their leaf onset period, and disregarding some of the NLs causes greater differences at the beginning of the growing season (Fig. 4b). Because plant transpiration contributes partially to the LE, these phenological effects also result in seasonal LE flux errors (Fig. 3e). However, note that despite the improvements in flux estimation, the surface tiling scheme still results in errors of considerable magnitude (greater than 7%) (Figs. 3d, 3e, and 3f). These errors in SSGH run are significantly correlated with degree of mixture between BL and NL PFTs, which results in higher errors on grid with increasing SSGH (figure not shown). Figure 5 shows the relationships between the mean anomalies in the GPP and the mean anomalies in the soil water limiting factor ( β) for each grid cell. Regardless of whether SSGH is considered or not, the anomalies in the GPP increase with the anomalies in the β. When SSGH is not considered, this correlation is much weaker, and 42% of the GPP anomalies can be explained by the changes in the β conditions. On the other hand, the surface tiling approach results in both smaller anomalies in the β and a stronger correlation between the β and GPP anomalies. The SWC anomalies can explain 90% of the GPP anomalies.

Effects of sub-surface soil water mixing on carbon flux estimation
This correlation in GPP anomaly and SWC anomaly is caused by disregarding spatial configuration of PFTs and their interactions with soil. The JULES LSM assumes that each surface tile in a grid cell shares the same meteorological and soil conditions . This results in at least one of the PFTs experiencing changes in soil water availability due to the surface tiling approach. The PFTs with higher water demands benefit more when they share soil water with PFTs that have lower demands. BLs exhibit phenological changes that generally result in less tree water use in both the spring and the autumn due to their leaf area reductions. On the other hand, they have higher photosynthetic and transpiration rates than NLs after their leaf areas have fully expanded under abundant water conditions. Consequently, sharing underlying soil water benefits NLs in both the spring and the autumn. BLs can be exposed less to soil water limiting conditions during the middle of growing season when sharing soil water with Fig. 5. Correlation between grid-mean soil water stress factor (β) anomaly and gross primary productivity (GPP) anomaly in surface subgrid heterogeneity (SSGH) and NOSSGH runs.  NLs, which results in higher photosynthetic rates in those periods. Field studies have also reported on the mixture effects caused by BL and NL PFTs on both the SWC and productivity (Falge et al. 2002;Pretzsch et al. 2015;Tang et al. 2014). However, the assumption of complete soil water mixing can result in erroneous flux estimation (Hartmann 2016;Rihani et al. 2010;Wagener and Gupta 2005). Sub-surface characteristics such as the soil depth, hydraulic conductance, and field capacity demonstrate significant heterogeneity across various scales and are strongly influenced by the adjoining vegetation. Soil pore geometry, macroscopic networks, and the amounts of organic compounds are all affected by the vegetation types and vertical root distribution patterns (Dunne et al. 1991;Jobbágy and Jackson 2000;Li and Shao 2006). Although many attempts have been made to include these effects in LSMs (Melton et al. 2017), most of the current LSMs do not stochastically consider them, either wholly or in part. In addition, most of the soil processes, such as soil moisture content calculation, are treated as grid-mean values, thus failing to consider sub-surface heterogeneity. Therefore, including both the subsurface heterogeneity and the interconnections between the vegetation and the soil in each grid cell can produce more accurate estimations of land surface fluxes (Maxwell and Miller 2005).
In addition, JULES LSM surface tiling relies on a statistical approach in which each grid cell is divided into fractional cover without any spatial information. However, the spatial configuration of the landscape has significant effects on both the surface energy fluxes and the multi-scale circulations (Oliphant et al. 2003;Zhou et al. 2011). Vertically discretized land surface modeling enables better estimations of the surface fluxes as well as superior reproductions of their temporal and spatial patterns (Naudts et al. 2015).

Conclusions and future work
This study numerically tests the JULES LSM surface tiling scheme and analyzes the impact of subgrid heterogeneity on the model-estimated surface fluxes using strategically designed numerical experiments with real land surface heterogeneity data. The results indicate that the current JULES LSM surface tiling scheme does not fully reflect the effects of SSGH, partly because it fails to consider the effects of subsurface heterogeneity. This may result in much greater errors under conditions of both limited soil water and the coexistence of different PFTs with distinct characteristics (tree-grass or tree-shrub mixtures) (Shrestha et al. 2016). The JULES community is attempting to introduce consideration of subsurface heterogeneity by incorporating the soil tiling scheme; We suggest that, in addition to the soil tiling scheme, more precise representation of interconnections with subsurface processes should be introduced in LSMs.