Machine learning accelerates parameter optimization and uncertainty assessment of a land surface model

The performance of land surface models (LSMs) strongly depends on their unknown parameter variables so that it is necessary to optimize them. Here I present a globally applicable and computationally efficient method for parameter optimization and uncertainty assessment of the LSM by combining Markov Chain Monte Carlo (MCMC) with machine learning. First, I performed the long-term ensemble simulation of the LSM, in which each ensemble member has different parameters' variables, and calculated the gap between simulation and observation, or the cost function, for each ensemble member. Second, I developed the statistical machine learning based surrogate model, which is computationally cheap but accurately mimics the relationship between parameters and the cost function, by applying the Gaussian process regression to learn the model simulation. Third, we applied MCMC by repeatedly driving the surrogate model to get the posterior probabilistic distribution of parameters. Using satellite passive microwave brightness temperature observations, both synthetic and real-data experiments were performed to optimize unknown soil and vegetation parameters of the LSM. The primary findings are (1) the proposed method is 50,000 times as fast as the direct application of MCMC to the full LSM; (2) the skill of the LSM to simulate both soil moisture and vegetation dynamics can be improved; (3) I successfully quantify the characteristics of equifinality by obtaining the full non-parametric probabilistic distribution of parameters.


Introduction
Land surface modeling is fundamental technology to understand, monitor, and predict terrestrial water, energy, and carbon cycles. As a first land surface model (LSM), Manabe (1969) developed the simple bucket model in which heterogeneity of soil, vegetation, and topography was neglected. For many years, LSMs have been evolving from the simple bucket model by including physiological processes (e.g., Sellers et al. 1986;Sellers et al. 1996;Dickinson et al. 1993), slope hydrology (e.g., Liang et al. 1994;Takata et al. 2003), groundwater processes (e.g., Koirala et al. 2014) and vegetation dynamics (e.g., Cox et al. 2000;Montaldo et al. 2005;Ivanov et al. 2008). As the complexity of LSMs increases, they have many unknown soil and vegetation parameters which cannot be directly observed. Therefore, it is necessary to infer those unknown parameters from the observation data which include the information of model's state variables by evaluating gaps between simulation and observation. This parameter optimization is the grand challenge in hydrologic data assimilation.
There are two major obstacles in the parameter optimization of LSMs. First, model parameters nonlinearly affect simulated variables so that the cost function, which maps parameters to the gaps between simulation and observation, may have a lot of local minima. In addition, model parameters mutually interact with each other (e.g, one parameter is sensitive to simulated variables only if another parameter is above a threshold). Therefore, it is not straightforward even to just find the set of sensitive parameters to observable variables (e.g., Rosero et al. 2010;Pappas et al. 2013;Sawada and Koike 2014). Second, different combinations of parameters show similar skills to reproduce observations so that many combinations of parameters can be equally acceptable, which is called equifinality in hydrologic literature (Beven and Freer 2001).
Searching a single optimal representation of parameters is clearly insufficient. A computationally efficient method to obtain the probability density function of parameters, which is called uncertainty assessment by some works in hydrology (e.g., Moradkhani et al. 2005), is crucially needed.
Despite a lot of efforts to optimize soil and vegetation parameters in LSMs, to our best knowledge, no study realized globally applicable parameter optimization and uncertainty assessment of LSMs by explicitly resolving the nonlinear effects of parameters and the equifinality. Towards a globally applicable method to optimize unknown parameters of LSMs, many previous studies have assimilated globally applicable satellite observations into LSMs. For example, Yang et al. (2007) developed the auto-calibration system, which can optimize soil parameters such as soil porosity, soil texture, and soil surface roughness, by assimilating microwave brightness temperature observed by Advanced Microwave Scanning Radiometer for Earth observation system (AMSR-E). Their system includes the Simple Biosphere model 2 (SiB2: Sellers et al. 1996) as an LSM, and a radiative transfer model to convert the LSM's state variables to microwave brightness temperature.
Although they successfully obtained a single optimal set of the parameters using the Shuffled Complex Evolution (SCE) algorithm (Duan et al. 1992), they provided no uncertainty assessment (see also Yang et al. 2005;Yang et al. 2009;Qin et al. 2009;Tian et al. 2009). Sawada and Koike (2014) extended the work of Yang et al. (2007) by incorporating a dynamic vegetation model into the LSM. Since microwave brightness temperature is sensitive to both surface soil moisture and vegetation water content, simultaneous estimation of both hydrologic and ecologic model parameters was realized and the skill of the LSM to simultaneously simulate both soil moisture and vegetation dynamics was improved. However, no uncertainty assessment was provided in their work. Bandara et al. (2015) optimized the parameters of the soil property in the Joint UK Land Environment Simulator (JULES: Best et al. 2011) using a satellite-observed soil moisture product based on Soil Moisture and Ocean Salinity (SMOS). They successfully retrieved a single optimal soil property by the particle swarm optimization (Kennedy and Eberhart 1995) without the information of uncertainty. Kato et al. (2013) optimized 24 ecological parameters using satellite-observed fraction of absorbed photosynthetically active radiation. They realized the estimation of both optimal parameters and their uncertainties by a variational method although they assumed the Gaussian probability density distribution of parameter uncertainties and the linear relationship between parameters and observations.
Many previous hydrologic studies proposed methods based on the Markov Chain Monte Carlo (MCMC) sampler and successfully applied them to parameter optimization and uncertainty assessment of lumped hydrological models (e.g., Vrugt et al. 2003;Vrugt et al. 2008;Schoups and Vrugt 2010;Vrugt et al. 2013) and distributed LSMs (e.g., Vrugt et al. 2003;Rosero et al. 2010). The limitation of those MCMC-based algorithms is that many trials of model integration ( (10 3~1 0 6 ) ) are required and they cannot be completely parallelized. Therefore, it is infeasible to directly apply the previously proposed MCMC-based algorithms to the globally distributed parameter optimization and uncertainty assessment of LSMs considering their computational cost.
Here I aim to accelerate parameter optimization and uncertainty assessment of an LSM using the technique of statistical machine learning based surrogate modeling which is theoretically investigated in the field of applied mathematics called uncertainty quantification (Sullvan 2015). Although this technique has been used for the parameter sensitivity analysis of atmospheric models (e.g., Qian et al. 2018), hydrological models (e.g., Dell'Oca et al. 2017Maina and Guadagnini 2018;Parente et al. 2019), and ecological models (e.g., Hawkins et al. 2019), few studies have applied it to parameter optimization and uncertainty assessment of LSMs with globally applicable satellite observations. In this study, a statistical surrogate model, which can mimic the relationship between model parameters and gaps between simulation and observation, is developed using machine learning. Then, using the statistical surrogate model, which is computationally cheaper than the original LSM, the probability distribution of parameters is sampled by the MCMC sampler. I demonstrate the potential of this approach in both synthetic and real-data experiments.

Problem statement
In this section, the problem of parameter optimization and uncertainty assessment is formulated. A discrete-time state-parameter dynamical system can be formulated as: where is the state variable (e.g., soil moisture, soil temperature and biomass) at time t, is the LSM, is the model parameters, is the external forcing (e.g., precipitation and incoming radiation), and is the noise process which represents the model error. Since the whole space of the state variables cannot be observed, it is useful to formulate an observation process as follows: where is the simulated observation (e.g., microwave brightness temperature) at time t, ℎ is the observation operator which is a radiative transfer model in this study, and is the noise process which represents the observation error.
where is the observation (e.g., microwave brightness temperature from a satellite) at time t. The purpose of parameter optimization and uncertainty assessment in this study is to obtain the posterior probability distribution of model parameters by the timeseries of simulation (3) and observation (4) from time 0 to T and the Bayes' rule ( | : ) ∝ ( 0: | ) ( ) where ( | : ) is the posterior probability distribution of given by the observations : , ( 0: | ) is the likelihood, and ( ) is the prior probability distribution of .
In the following section, I explain the LSM ( in equation (1)), called EcoHydro-SiB, used in this study. The targeted model parameters are also shown there. In section 2.3, I explain the observation operator (ℎ in equation (2)) used in this study. The method to efficiently obtain the posterior probability distribution of parameters ( ( | : ) in equation (5)) is explained in section 2.4.

Land surface model
EcoHydro-SiB, the improved version of SiB2 (Sellers et al. 1996), was developed as a module of a previously developed land data assimilation system, called Coupled Land and Vegetation Data Assimilation System (CLVDAS: Sawada and Koike 2014;Sawada et al 2015;Sawada 2018). In this section, the schemes related to parameter optimization and uncertainty assessment are explained. Please refer to Sawada and Koike (2014) for the complete description of this LSM. The variables of parameters which are not optimized can also be found in Sawada and Koike (2014) and references therein. In this study, Ecohydro-SiB was run with the horizontal resolution of 0.25 degree and the timestep was 1 hour.
EcoHydro-SiB solves vertical interlayer water flows from surface to 2m depth using a one-dimensional Richards equation. The depth of surface soil layer was set to 5cm and that of the other soil layers was set to 10cm. To solve the Richards equation, capillary suction ψ (m) and hydraulic conductivity (m/s) are formulated as functions of volumetric soil moisture w (m 3 /m 3 ) by the van Genuchten water retention curve (van Genuchten 1980): where is the saturated hydraulic conductivity (m/s), is the residual water content (m 3 /m 3 ), is the saturated water content or porocity (m 3 /m 3 ), and are the parameters. In this study, and n are the targeted parameters to be optimized.
EcoHydro-SiB simulates vegetation growth and senescence by solving a carbon balance equation: where , , and , are the carbon pools of leaves, stems, and roots, respectively (kg/m 2 ), , , and are the carbon allocation fractions of leaves, stems, and roots, respectively, , , and are the normal turnover rates of leaves, stems, and roots, respectively, and are the stress factor related to water stress and temperature stress, respectively. See Sawada and Koike (2014) for the equations to calculate the allocation fractions and stress factors. NPP is the net primary production (mol m -2 s -1 ) and is calculated by the photosynthesis-conductance model. In this model, NPP is the function of the maximum Rubisco capacity of the top leaf, defined as 0 in this paper. See Sawada and Koike (2014) and Sellers et al. (1996) for details of the photosynthesis-conductance model. The linear relationship between the carbon pool of leaves and leaf area index (LAI) is assumed: where is the specific leaf area (m 2 /kg). Since stem and root biomass is required to support leaves, the following constraint is formulated: where is the model parameter (see also Ivanov et al. 2008). If the relationship shown in equation (13) is violated, is set to zero and no carbon biomass is allocated to leaves. In this study, 0 and are the targeted parameters to be optimized.  Mo et al. (1980). Since this model explicitly calculates the emission and absorption by vegetation water, simulated microwave brightness temperature strongly depends on LAI calculated by equations (9)-(12). Land surface soil emissivity was calculated by the advanced integral equation model with the incorporation of a shadowing effect (Kuria et al. 2007). Since this model explicitly calculates the dielectric constant of the soil-water mixture, simulated land surface emissivity strongly depends on surface soil moisture. Please refer to Kuria et al. (2007) and Sawada and Koike (2014) for the complete description of the radiative transfer model.

MCMC with a statistical surrogate model
As shown in section 2.2, there are four LSM's parameters to be optimized in this study: saturated hydraulic conductivity , parameter of van Genuchten water retention curve , maximum Rubisco capacity of the top leaf 0 , and the factor controlling the relation between the carbon pools . The four optimized parameters were normalized by the followings: where 1 , 2 , 3 , and 4 are the optimized parameters which range from 0 to 1, the superscript def means the "default" variable for each parameter. The default variables were specified according to the global soil and land use map (see section 3). For instance, can move from 1 ( 1 = 0) to 1 ( 1 = 1) in the optimization scheme. The possible ranges of parameters such as 1 and 1 are summarized in Table 1.
The prior ( ) ( = [ 1 , 2 , 3 , 4 ] ) in equation (5) is assumed to be the uniform distribution. It should be noted that the selection of the prior ( ) , which implicitly includes the selection of the targeted parameters, the possible ranges of the parameters, and the shape of the distribution, was subjective although it was based on the previous works' experiences of parameter sensitivity analysis and parameter calibration. The selection of the prior ( ) significantly affects the results of parameter optimization and uncertainty assessment. However, in this study, I focused on the efficient parameter optimization and uncertainty assessment given the specific prior and the objective selection of the prior is not the scope of this study.
The posterior ( | : ) in equation (5) was obtained by the MCMC sampler. I used the Metropolis-Hastings algorithm (Hastings 1970). First, the initial parameter vector 0 was specified. Then, the processes shown below were iterated: 1. For each iteration i, generate a candidate parameter vector . It is sampled from the distribution q( | ).
2. Evaluate the cost function C( ) . The cost function C should be the indicator of the difference between simulation and observation. Larger C indicates that the gap between simulation and observation is smaller. If b ≥ a, reject the candidate parameter vector and + = In this study, the distribution q( | ) was set to the Gaussian distribution whose mean and standard deviation were set to and 0.1, respectively. The cost function C should be proportional to the targeted distribution ( | : ). The cost function C was defined using the root-mean-square error (RMSE): where is the observation error in the satellite-observed brightness temperature and was set to 1 (K) in this study. See also equations (2)-(4).
It is infeasible to directly apply the MCMC sampler to parameter optimization and uncertainty assessment of the LSM because whenever the cost function C is evaluated, I need to run the LSM and the radiative transfer model for the long period which includes the period of spinup to obtain the simulated observation formulated in equation (3). The Metropolis-Hastings algorithm requires the 10 4~1 0 6 iterations and these iterations cannot be parallelized. Therefore, the MCMC-based algorithms which are more efficient than the direct application of the Metropolis-Hastings algorithm have been proposed (e.g., Vrugt et al. 2003;Vrugt et al. 2008) and I used a statistical surrogate model. To develop the statistical surrogate model, g( ), RMSE was normalized: where and are the maximum and minimum RMSEs in the 400 ensemble members, respectively. The Gaussian process regressor learned the relationship between and and predicted as a function of . Then, g( ) predicted the cost function: where ( ) is the predicted normalized RMSE by the Gaussian process regression.
In summary, the proposed method is the following: 1. Generate 400 ensemble members of model parameters by the Latin hypercube sampling.
2. Drive the LSM and the radiative transfer model and calculate RMSE using those parameter ensembles in parallel.
3. Construct the statistical surrogate model g( ) from the dataset of the -C combinations by the Gaussian process regression.
4. Run the Metropolis-Hastings algorithm with 10 5 iterations in which the cost function C( ) was evaluated by the statistical surrogate model g( ) .

Data
The International Satellite Land Surface Climatology Project 2 (ISLSCP II) data (Global Soil Data Task Group, 2000) and the Food and Agricultural Organization (FAO) global dataset (Food and Agricultural Organization 2003) were used to derive the LSM's default parameters (see also section 2.4).
The Global Land Data Assimilation System (GLDAS) v2.1 meteorological forcing (Rodell et al. 2004;Sheffield et al. 2006) was used to drive the LSM. The meteorological forcing data necessary to drive the LSM are surface pressure, precipitation, surface air temperature, relative humidity, incoming solar radiation, incoming longwave radiation, and wind speed. The original temporal resolution of this data is 3-hourly and the data were linearly interpolated to hourly, which is consistent to the timestep of the LSM (see section 2.2). The original horizontal resolution of this data is 0.25 degree and is consistent to the horizontal resolution of the LSM (see section 2.2).
The observed microwave brightness temperature is from the AMSR-E L3 product provided by Japan Aerospace eXploration Agency (JAXA) (Kachi et al. 2013 v3.2 product (Dorigo et al. 2017) was used to validate the skill of the LSM to simulate surface soil moisture. It should be noted that ESA CCI SM is not completely independent to the AMSR-E since it includes the AMSR-E based soil moisture retrieval. The spatial resolution of the data is 0.25 degree which is consistent to the LSM (see Section 2.2).

Study area
The study area is a part of the Sahel region and is shown in Figure 1. This study area was chosen because there was a steep gradient of climate and landscape from the rain forest in the southern part of the study area to the savanna and desert in the northern part of the study area. I could discuss how the proposed algorithm works in the different climate and vegetation conditions. The Sahel area has been chosen as study areas in the other land data assimilation studies (e.g., Albergel et al. 2018) by the same reason.

Synthetic experiment
To deeply discuss the advantages and limitations of the proposed algorithm, an ideallized synthetic experiment was performed. In the synthetic experiment, I specified the synthetic true model parameters, shown in Table 1, and drove the LSM to generate the synthetic true land state. From the synthetic true land state, the synthetic observation (i.e. microwave brightness temperature at 6.925 and 10.65 GHz) was generated using the radiative transfer model. Using this simulated observation by the synthetic true parameters, the algorithm described in section 2.4 was performed. The timing of the synthetic observation was the same as real observations from AMSR-E (see Section 3).
The meteorological forcing and the default parameter variables were retrieved from the real data described in Section 3.
This synthetic experiment was performed in the three points, the Site I-III shown in Figure   1b

Real-data experiment
In addition to the idealized synthetic experiment, I demonstrated the potential of the proposed algorithm in the real-world application. Here I used the real AMSR-E observed brightness temperature data for parameter optimization and uncertainty assessment. The study area was shown in Figure 1b. The study period and the method to spinup the model were identical to the synthetic experiment (see section 4.2). Although the "true" parameters cannot be known in the real-world application, the performance of the   (14)) were randomly drawn from 0 to 1, but the other parameters ( 2 , 3 , and 4 in equations (15)- (17)) were fixed to the synthetic truth shown in Table 1. Because the synthetic truth of 1 is set to 0.75 (see Table 1), RMSE is small when 1 is around 0.75. As 1 deviates from 0.75, RMSE increases. Figure 2a indicates that the response of parameters to RMSE is mostly linear except for n in the van Genutchen water retention curve (orange dots in Figure 2a).  (5), with the simulation driven by the parameters from the prior, ( ) in equation (5) (i.e. random draws from the uniform distribution [0,1]). Figure 5 shows that the parameter optimization successfully improves the skill of the LSM to simulate observable microwave brightness temperature. Figure 6 indicates that the parameter optimization also improves the skill of the LSM to reproduce the synthetic truth of LAI and soil moisture. The parameter optimization substantially reduces the uncertainty in the simulation of LAI and the median of LAI simulated by the MCMC-sampled parameters was more consistent to the synthetic truth than the simulation by the completely random samples of the parameters (Figures 6a and   6b). The parameter optimization also reduces the uncertainty in the simulation of surface soil moisture especially in the dry seasons (Figures 6c and 6d). Although the root-zone soil moisture cannot be directly observed by the satellite, the parameter optimization positively impacts to the simulation of root-zone soil moisture (Figure 6e and 6f) because the parameter adjustment affects the whole state space of the LSM.
In the Site II (0.625E, 7.125N see Figure 1b), annual precipitation is approximately 1590 (mm/yr). In the synthetic truth, LAI changes from 2 to 4. In this densely vegetated condition, microwave emission from land soil surface is strongly attenuated by canopy and microwave brightness temperature is always insensitive to surface soil moisture (Sawada et al. 2017) so that the observation has the information of only vegetation dynamics in the whole period. In addition, the contribution of vegetation water content to microwave radiative transfer becomes saturated when LAI is large (see Sawada et al. 2017). Therefore, the observation is less informative in the Site II than the Site I. However, the posterior distributions of the parameters recover the synthetic truth less accurately than the Site I. Figures 7a, 7b, 7c, and 7d show that the posterior distributions of 1 , 2 , and 3 in the Site II is flatter than the Site I and their mode is not consistent to the synthetic truth (red dashed lines). Although the sampled parameters successfully minimize RMSE between simulated and observed brightness temperatures (Figure 8), the observation is not very informative to constrain the four unknown parameters in the densely vegetated area. Figure 7e shows that there are correlations between sampled parameters such as 1 -2 , 2 -3 , and 3 -4 . The advantage of the proposed method is that the equifinality due to the insufficiently informative observation can be explicitly detected by the relatively cheap computational cost.
The parameters drawn from the posterior distributions substantially reduce the uncertainty in the simulation of LAI (Figures 9a and 9b). Figure 9b shows that the simulation with the optimized parameters still has the relatively large uncertainty in the seasonal peak of LAI. This is because the contribution of vegetation water content to microwave radiative transfer becomes saturated when LAI is large (e.g., LAIs of 3.5 and 4.5 cannot be accurately distinguished by the microwave signal; see Sawada et al. 2017) .   Figures 9c, 9d, 9e, and 9f indicate that the optimization has a smaller impact on the soil moisture simulation than the Site I.
In the Site III (0.125E, 20.125N), annual precipitation is approximately 110 (mm/yr). This virtual site is located in the arid area and there is no vegetation in the synthetic truth. Figure 2c indicates that the optimization problem in the Site III is completely different from those in the Sites I and II. Parameters of 1 , 3 , and 4 are not sensitive to RMSE and the gap between simulation and observation is mainly controlled by the n of the van Genuchten water retention curve ( 2 ). Moreover, the response of 2 to RMSE is quite non-linear. RMSE suddenly increases when 2 is larger than ~0.5.
The Gaussian process regressor is less accurate to predict RMSE calculated by the LSM and the radiative transfer model than the cases of the Sites I and II (Figure 3c) although the predicted and true RMSEs are reasonably correlated (R 2 =0.89). Since the parameter range of the transition from low RMSE to high RMSE is small (Figure 2c), there should be a large number of ensembles to accurately include this transition in the training dataset and/or more complex statistical models might be needed to resolve the nonlinear response of the parameters.
The proposed method accurately samples the posterior distribution even in this extreme case. Figures 10a, 10c, and 10d show that the parameters 1 , 3 , and 4 are uniformly distributed in the range [0,1], which implies little sensitivity of these parameters to microwave brightness temperature. Figures 10b shows the parameter 2 is highly sensitive to microwave brightness temperature. This sensitivity analysis is consistent to the results shown in Figure 2c. The mode of the posterior distribution of 2 deviates from the synthetic truth due to the small sensitivity of the change in 2 from 0.2 to 0.5 (see the flat orange line in this range in Figure 2c). Figure 11 shows that this parameter optimization succesfully minimizes the gap between simulated and observed brightness temperatures. Figures 12c and 12d indicate that the parameter optimization significantly reduces the uncertainty in surface soil moisture. However, the skill to simulate root-zone soil moisture is degraded, which is not consistent to the results obtained in the Sites I and II. Although there is the sensitivity to root-zone soil moisture in the range of the sampled parameters, all of the 400 ensembles simulate the almost identical brightness temperatures ( Figure   11), LAI ( Figure 12b) and surface soil moisture (Figure 12d). Therefore, the optimization reveals that the observation has the insufficient information to constrain root-zone soil moisture in this site. Generally, in arid regions, surface variables (surface soil moisture and vegetation) are not well correlated with root-zone soil moisture (e.g., Kumar et al. 2009) so that it is difficult to infer root-zone soil moisture from surface observations.

Real data experiment
Here I demonstrate the advantages and limitations of the proposed method in the real data experiment in which the real AMSR-E brightness temperature observations are assimilated into the LSM. While in the idealized synthetic experiment I assumed that the 4 selected unknown parameters are the sole source of errors, in the real data experiment I neglected many other error sources such as the parameters which are not optimized, the meteorological forcing ( in equation (1)), the model structure ( in equation (1)), and the observation process ( in equation (2)). This point should be noted for the interpretation of the results. Figures 13a, 13b, 13c, and 13d show the medians of the 10 5 sampled parameters of 1 , 2 , 3 , and 4 , respectively. The proposed method can obtain the spatially distributed optimal parameters by the reasonably cheap computational cost. Compared with the other parameters, the saturated hydraulic conductivity prefers the default variable ( Figure  13a). Figure 13b shows that the optimal variables of the van Genuchten's n tend to be smaller than the default variables in the wetter area while they tend to be larger than the default variables in the drier area. Figures 13c and 13d shows that the optimal variables of maximum Rubisco capacity of the top leaf 0 and the factor controlling the relation between the carbon pools tend to be larger and smaller than the default variables, respectively.
As shown in the synthetic experiment of the virtual Site III, the posterior distribution of the parameters with no significant sensitivity to observation tends to be a uniform distribution (see Figure 2c and Figures 10a, 10c, and 10d). Therefore, the Kullback-Leibler divergence (KLD) (Kullback and Leibler 1951) between the uniform distribution and the posterior distribution can be a good index of the parameter sensitivity: where ( , ) is the KLD between two probabilistic distributions, p, and q. If two distributions are equal for all , ( , ) = 0. If two distributions greatly differ from each other, the large value of ( , ) is obtained. Therefore, the KLD can be used as an index to evaluate the closeness of the two probabilistic distributions. In this study, large KLD means that the sampled posterior distribution of the parameters substantially differs from the uniform distribution and the parameter is sensitive to observations. Figures 13e, 13f, 13g, and 13h show the KLDs of the parameters of 1 , 2 , 3 , and 4 , respectively. The saturated hydraulic conductivity tends to be uniformly drawn and less sensitive to observation than the other parameters (Figure 13e). Figure 13f shows that the van Genuchten's n is sensitive in the whole study area. This is because the van Genuchten's n affects not only soil water dynamics but also the soil moisture thresholds in vegetation water stress (i.e. field capacity and wilting point), which is crucially important to the simulation of vegetation growth and senescence. Figures 13g and 13h show that the maximum Rubisco capacity of the top leaf 0 and the factor controlling the relation between the carbon pools are sensitive only in the vegetated area (see also the distribution of observed LAI in Figure 1b). It is not straightforward to evaluate the ensemble simulation. Here I defined the two indices to evaluate the skill of the LSM with 400 parameter ensembles to reproduce observed LAI and surface soil moisture. First, I calculated the bias for each ensemble member and took its root-mean-square: where is the ensemble size, is the simulated observation by the ensemble , is the observation, and − ̅̅̅̅̅̅̅̅̅̅ is the temporal mean of the simulation-observation differences, which is bias for the ensemble . Second, I calculated the unbiased rootmean-square error (ubRMSE) by removing the bias between observation and ensemble mean. where E( ) is the expectation operator and ( ) is calculated as the temporally averaged ensemble mean. Finally, I calculated the improvement rate for each score (i.e., and ) as follows:

= −
where is the improvement rate, is the evaluation metrics ( or ), the subscript mcmc means the evaluation metrics from the simulation with the MCMC sampled parameters, and the subscript unif means the evaluation metrics from the simulation with the parameters sampled from the prior uniform distribution. Figures 15a and 15b show the spatial distribution of the improvement rates for and of LAI against GLASS LAI (see section 3), respectively. Although the was reduced by the parameter optimization in the large portion of the study area, there is the substantial degradation in the west part of 10N-12N. In this area, the optimized parameters 1 and 2 are unnaturally small (Figures 13a and 13b) due probably to the overfitting induced by the neglected errors such as errors in the model structure, the radiative transfer model, and the meteorological forcing. This overfitting might degrade the skill of the LSM to simulate LAI. Although the impact of the parameter optimization on of LAI is relatively small, there is the degradation in the part of the vegetated area. In this area, the seasonal cycle of the simulated LAI is slightly delayed to that of observed LAI (not shown), which can also be found in the other LSMs (e.g., Jarlan et al. 2008). This LSM's structural error cannot be fixed in the current parameter optimization framework. Figures 15c and 15d show the spatial distribution of the improvement rates for and of surface soil moisture against ESA CCI (see section 3), respectively.
The parameter optimization greatly increases in the northern arid area.
However, in this area comes mainly from the gap between simulated and observed minimum baseline surface soil moisture in the dry seasons (not shown).
Considering the error in the satellite soil moisture retrievals (approximately 0.05 m 3 /m 3 ) and the small impacts of this baseline soil moisture on the water, carbon, and energy balances, this gap should be removed for evaluation (e.g., Entekhabi et al. 2010). Figure   15d shows that the parameter optimization reduces more than 10 % in many parts of the study area.

Conclusion
In this study, the globally applicable method for parameter optimization and uncertainty assessment of the LSM was developed. The characteristics of the proposed method are the followings. First, the proposed method can retrieve the non-parametric posterior probabilistic distributions of the unknown parameters so that it is robust to the nonlinearity of the effects of the parameters. Second, from the sampled non-parametric posterior probabilistic distributions, the equifinality of the LSM's parameters can be detected. Third, I effectively used the satellite observation of microwave brightness temperature which is sensitive to the important ecohydrological variables (i.e., soil moisture and vegetation water content) and has the global and all-weather capability.
Fourth, and most importantly, the reasonably cheap computational cost was achieved by the technique of the statistical machine learning based surrogate modeling, which makes the proposed method globally applicable. The advantages described above were demonstrated in both the idealized synthetic experiment and the real data experiment.
Towards the global-scale parameter optimization and uncertainty assessment of the LSMs, several challenges remain. First, the prior should be carefully designed. The prior stated in this study implicitly includes the selection of parameters, the range of parameters, and the shape of distributions. In addition, the explicit consideration of the error in meteorological forcing and model structure is helpful to further improve the performance of the LSM's simulation. It can be done using more than one dataset of meteorological forcing and multimodel ensembles. Second, in addition to the microwave satellite observations such as AMSR-E, there are many other important satellite observations which can contribute to parameter optimization and uncertainty assessment. Extending the framework proposed in this study to explicitly consider a number of observation types (i.e. a number of targets or objectives) is an important research topic.

Figure 2.
Simulated RMSEs of brightness temperature between the LSM with the synthetic true parameters and that with the random draws of the normalized parameters of saturated hydraulic conductivity (blue), the van Genuchten's n (orange), maximum Rubisco efficiency at the leaf top (green), and the factor controlling the relation between the carbon pools (red) in (a) the Site I, (b) the Site II, and (c) the Site III. In each dot, the single targeted parameter was randomly drawn, and the other parameters were fixed to the synthetic truth (see Table 1 and section 5.1).