1 Introduction

The microwave radiative transfer model (RTM) provides a link between the land surface state variables (e.g., soil moisture and temperature) and the satellite-observed brightness temperatures (T_b). It plays a key role on the satellite data inversion and land data assimilation. Several RTMs^[1-5] have been proposed to estimate T_b at the top of atmosphere over various surface conditions. However, RTMs have been reported to exhibit significant uncertainties, which are mainly arising from the formulations of soil roughness, vegetation optical thickness, soil dielectric constants^[6-7], parameter values^{[5, 8]}, land surface variables, and auxiliary fields^[9-10], and affect large-scale forward T_b simulation^{[6, 8]}, soil moisture retrieval^[11], and the performance of a satellite data assimilation system^[12].

How to reduce the uncertainties in microwave RTMs is a major challenge, especially at a large spatial scale. De Lannoy et al.^{[5, 8]} used the particle swarm optimization algorithm and the Markov chain Monte Carlo (MCMC) method to determine the optimal parameter values, respectively. However, the parameter calibration method is limited to the number of the parameters, since RTMs include a variety of parameterizations and parameter values. The multi-model ensemble strategy^[12] is seen as a good choice to improve the model simulation by exploiting the diversity of the individual model predictions with efficient weighting schemes. The weights for each ensemble member can be the same or established by certain statistical schemes. The Bayesian model averaging (BMA) is a statistical post-processing method to produce a probabilistic forecast from the original ensembles^[13]. It has gained popularity in the diverse fields such as medicine, economic, hydrology, and meteorology^[13-16]. The BMA scheme can make use of the advantages of each individual mode and provide a better description of the predictive uncertainty, which is superior to some multi-model methods.

In this study, we use the BMA method to improve the microwave T_b predictions from the community microwave emission model (CMEM)^{[4, 7]} to simulate the passive microwave T_b at low frequencies. The CMEM is coupled to a land surface model, i.e., the community land model version 4.5 (CLM4.5)^[17], which can provide the required input for the CMEM^[18]. The ensemble for the BMA scheme is generated from different vegetation opacity parameterizations (VOPs), soil dielectric constant parameterizations (SDCPs), and soil roughness parameterizations (SRPs), because the three factors play key roles in the simulation of T_b at the the top of atmosphere^[19-20].

This study aims to investigate how the BMA scheme can improve both the accuracy and the reliability of the T_b predictions. Particularly, we examine the performance of the BMA scheme with the CLM4.5-CMEM over different vegetation types. The paper is organized as follows. Section 2 gives a brief description of the CMEM and CLM4.5. Section 3 describes the BMA methodology and experiment design. The evaluation results of the BMA scheme are presented in Section 4, and the conclusions are shown in Section 5.

2 Model descriptions 2.1 Community land model

The CLM4.5 from the National Center for Atmospheric Research in the United States is a community-developed land model for simulating biogeophysical and biogeochemical processes. It contains several modifications over previous versions, e.g., reducing the biases of soil carbon in the arctic and revising canopy radiation and photosynthesis schemes^[17]. CLM4.5 represents the land surface as a hierarchy of sub-grid types, e.g., glacier, lake, wetland, urban, and vegetated land. The vegetated land is composed of multiple plant functional types (PFTs). Biogeophysical and biogeochemical processes are simulated for each sub-grid. The detailed descriptions of CLM4.5 are given in Ref. [17].

The forcing dataset from the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (hereafter ITP)^[21] is used to run the CLM4.5 in this study. It has a spatial resolution of 0.1° × 0.1° and a three-hourly temporal resolution (1979−2012), which covers mainland China (15°N−55°N, 70°E−140°E). This dataset is generated by combining the ground-based measurements, the satellite-retrieved products, and other data from the global land data assimilation system^[22-23]. The VOP in the CLM4.5 uses 16 default vegetation classes. For the CMEM simulations, these classes are further refined into the six classes defined in Ref. [24], in which the C4 crops are not included in the CLM4.5. Figure 1 shows the surface characteristics over mainland China.

2.2 CMEM

The CMEM is developed by the ECMWF and used to simulate the low frequency passive microwave T_b (1 GHz−20 GHz) of the surface^{[4, 7]}. It includes four modules to compute the contributions from soil, vegetation, snow, and atmosphere. The polarization T_b over the snow free areas at the top of atmosphere (T_{Btoa, p}) is expressed as follows:

(1)

(2)

where τ_{atm, p} represents the atmospheric optical depth. T_{Bau, p} and T_{Bad, p} are the upward and downward atmospheric emissions defined by

(3)

(4)

where T_{eff, p} is the effective temperature of the soil medium, T_c is the canopy temperature, and ω_p is the single scattering albedo. T_{Btov, p} is the T_b at the top of vegetation. T_{Bsoil, p} and T_{Bveg, p} are the soil and vegetation layer contributions, respectively. r_{r, p} is the soil reflectivity of the rough surface (one minus the emissivity e_{r, p}), and τ_{veg, p} represents the vegetation optical depth. Although the CMEM has a snow emission model, the effect of the snow is ignored in this study.

The CMEM has a modular choice of the physical parameterizations for three dielectric layers, i.e., soil, vegetation, and atmosphere, and includes 1 440 different configurations^[6]. The vegetation optical depth and soil dielectric constant are two key components of the microwave emission model, and make the major contributions to the simulation of T_b at the top of atmosphere^[19]. Moreover, soil roughness has important effects on the T_b simulation^{[20, 25]}. Therefore, 24 configurations of the CMEM are used in this study to generate the ensemble (see Tables 1 and 2) for the BMA scheme, which concerns four vegetation opacity models, three soil dielectric models, and two soil roughness models. The other parameterizations are the same as those in Ref. [6]. Here, general descriptions of the SDCPs, the SRPs, and the VOPs are presented (see Refs. [4] and [7] for the details).

2.2.1 Soil dielectric model

Three semi-empirical soil dielectric models are available in the CMEM, i.e., Dobson et al.^[1], Mironov et al.^[26], and Wang and Schmugge^[27] (hereafter Do, Mi, and Wa, respectively, see Table 1). The Wa and Mi models consider the effects of the bound water and are valid for 1 GHz-10 GHz^[6], while the Do model is available between 1 GHz and 18 GHz and is valid for a larger range of soil texture types^[26].

2.2.2 Soil roughness model

The CMEM uses a semi-empirical approach to consider the effect of the soil roughness^[4], i.e.,

(5)

where r_{s, p} is the soil reflectivity of the smooth surface, Q represents a polarization mixing factor, N is the angular dependence, h is the roughness parameter, and θ is the incidence angle. Although five parameterizations are available in the CMEM, only two are used in this study. This is because that other parameterizations include the dependencies on the correlation length and the clay fraction, which may introduce large uncertainties^[25]. The first one is the simplest parameterization proposed by Choudhury et al.^[28] (hereafter Ch), i.e.,

(6)

(7)

where σ is the standard deviation of the soil roughness, f is the frequency, and k is the wave number. Another approach proposed by Wegmüller and Mätzler^[29] (hereafter We) is

(8)

(9)

(10)

where the subscripts v and h represent the vertical polarization and the horizontal polarization, respectively.

2.2.3 Vegetation opacity model

In the CMEM, the T_b at the top of the vegetation is represented through the so-called τ-ω approach (see Eq. (4)). Jackson and Schmugge^[30] (hereafter Ja) proposed a simple parameterization to express the optical thickness as a function of the vegetation water content (V_WC) and an empirical vegetation structure parameter b as follows:

(11)

Wigneron et al.^[3] (hereafter Wi) calculated τ_{veg, p} by

(12)

(13)

where t_p represented the angular effect on τ_{veg, p}, τ_nadir was the nadir optical depth, N_LAI was the leaf area index, and b′ and b” were the vegetation structure parameters. Kirdyashev et al.^[31] (hereafter Ki) used a different method to compute τ_{veg, p} between 1 GHz and 7.5 GHz as follows:

(14)

where k was the wave number, ε_vw represented the dielectric constant of the vegetation water (the imaginary part), ρ_water was the water density, and a_geo was a vegetation structure parameter. Wegmüller et al.^[32] (hereafter We) considered k in the attenuation along the viewing path, and extended the scheme (1 GHz-100GHz).

3 Data and method 3.1 Microwave observations

The brightness temperatures simulated by the CMEM are compared with the C-band satellite observations (6.9 GHz) from the advanced microwave scanning radiometer on the Earth observing system (AMSR-E) gridded datasets (0.25°×0.25°), produced by the National Snow and Ice Data Center. The observations are twice daily, i.e., 1:30 p.m. (ascending) and 1:30 a.m. (descending), with an incidence angle of 55°^[33]. The vertical polarization is relatively insensitive to the vegetation coverage^[34]. Thus, the T_b data at the vertical polarization can be used to evaluate the microwave RTM^[18] or directly assimilated into the land surface-hydrology model to improve the soil water content estimations^{[12, 35-36]}. Previous studies^[37-38] suggest that the radio frequency interference (RFI) at the C-band may affect the use of the AMSR-E data in North America, South Africa, Europe, Middle East, and Japan, while the accuracy of the T_b observations over other areas is very good (less than 1K). In addition, the descending pass AMSR-E observations may be affected by the sun-glint effects^[6]. Therefore, the vertically polarized data at 6.9GHz covering the period from 2003 to 2011 for the descending orbit over mainland China are chosen here to calibrate the BMA weights and evaluate the CLM4.5-CMEM outputs.

3.2 BMA approach

The BMA algorithm^[13] is a statistical postprocessing method to produce the forecast probability density functions (PDFs) from ensembles, and provides a reliable description of the total modeling uncertainty. The BMA predictive PDF of any quantity of interest is a weighted average of the PDFs centered on the bias-corrected forecasts from a set of individual models. The weights are the estimated posterior model probabilities, and represent the relative predictive skill for each ensemble member. The BMA predictive PDF is expressed as follows:

(15)

where y represents the predictive variable of interest (T_b in this study), f = f₁, ···, f_N is an ensemble of predictions (24 configurations of the CMEM shown in Table 2), w_k denotes the posterior probability of the forecast k, and g_k(y|f_k) is the conditional PDF of y. The weight w_k is non-negative, satisfying

(1)

The probability distribution of the brightness temperatures is represented by a normal distribution here. Following Vrugt et al.^[39], we use the differential evolution adaptive metropolis (DREAM) MCMC algorithm to estimate the BMA weight w_k and variances.

3.3 CLM4.5-CMEM experiments

The CLM4.5 is first run 300 years for the spin-up through repeating 30-year (1979−2009) ITP dataset. Then, we use the 1979−2012 forcing data to run CLM4.5 so as to generate the simulated land surface dynamic fields between 2003 and 2012 (e.g., soil water content, soil temperature, canopy temperature, and 2-m air temperature), which are used as the inputs to the CMEM. The static fields required by the CMEM (see Fig. 1), including soil texture, vegetation types, and elevation, are from the surface dataset in the CLM4.5. It should be noted that a quality check for both the AMSR-E measurements and the CLM4.5-CMEM predictions between March and October is taken into consideration to eliminate the data affected by the snow cover (snow water equivalent > 10⁻⁴ kg·m⁻²) and temperature below 0℃.

In the present study, the period of 2003-2006 is chosen as the training period, while the period of 2007-2010 is chosen for the evaluation. We calibrate each ensemble member to obtain the BMA weights and variances for each grid cell independently. To ensure a meaningful calibration at each grid cell, we impose a minimum of 20 valid data points (N_i) per year to compute the BMA weights for each member. The requirement of N_i > 20 is used for the calculation of the evaluation statistics as well.

3.4 Statistical measures

The following two statistical indicators are defined here to evaluate the performances of each model member and BMA simulation:

(16)

(17)

where D_RMSD is the root-mean-squared difference (RMSD), r is the correlation coefficient, S is the T_b from the model simulation, O is the AMSR-E observation, and σ_S and σ_O represent the standard deviations.

4 Results 4.1 Spatial distribution of C-band brightness temperature over China

Figure 2 shows the spatial distributions of the mean T_b at the C-band averaged during the validation period of 2007-2010 at the vertical polarization for the descending orbit from the AMSR-E observations and model simulations. The CLM4.5-CMEM simulations show similar broad patterns, and capture the spatial variations of the AMSR-E observations well (spatial correlations r > 0.9). However, they show an overestimation of the observed T_b over northwest China and Tibetan Plateau and an underestimation for the densely vegetated conditions (e.g., southeast and northeast China). Moreover, different parameterizations in the CMEM differ in the magnitude with varying mean biases. The Ja (P1, P5, P9, P13, P17, and P21) and Wi (P4, P8, P12, P16, P20, and P24) vegetation opacity models significantly underestimate T_b over most of the vegetation-covered areas, even with negative biases less than -100 K for the forest areas. By contrast, the other two vegetation opacity models, i.e., Ki and We, reduce the underestimations, and capture the spatial variations better. It can be seen from Fig. 3 that the BMA weights have a close relationship with the individual model performance. The Ki and We parameterizations have larger weights than the Ja and Wi parameterizations. Moreover, the CLM4.5-CMEM simulations show larger discrepancies among different PFTs over south and northeast China, which suggests that the vegetation contribution dominates the signal of the T_b simulation over the vegetation-covered areas. For northwest China and the Tibetan Plateau, the surface emissions are mainly controlled by the soil dielectric and soil roughness models due to the sparse vegetation. Relatively, the soil roughness has a larger effect on the T_b simulation than the soil dielectric constant over the bare soil, e.g., northwest China. Compared with the We SRP (P13-P24), the Ch scheme (P1-P12) shows higher T_b values, and even leads to an overestimation for the We vegetation opacity model over southern China (P3, P7, and P11) and for the Ki model over the crop areas (P2, P6, and P10).

The SAA method, in which each model has the same weight, can improve the T_b simulations over the densely vegetated conditions compared with the Ja and Wi parameterizations, e.g., a reduction in the underestimation over southeast and northeast China. However, the simply arithmetic ensemble mean does not improve the accuracy of the T_b predictions for all conditions. For example, the Ki and We VOPs are closer to the AMSR-E observations than the SAA observations. In contrast, the BMA-based reconstructions reproduce the spatial variation of the AMSR-E observations best, with the highest spatial correlation coefficient (r > 0.99). Moreover, a single model simulation by the SAA method shows large biases in simulating the T_b fields over the Qinghai-Tibetan Plateau, the Taklimakan Desert, and southeast and northeast China, but these biases are reduced significantly by the BMA scheme.

4.2 Statistical evaluation of temporal variation

The RMSDs and the correlation coefficients between the CLM4.5-CMEM simulations, including 24 members, the BMA scheme, and the arithmetical average, and the AMSR-E observations, are presented in Figs. 4 and 5, respectively. It can be found from Fig. 4 that the Ja and Wi VOPs show higher RMSDs, e.g., more than 100 K over parts of southwest, southeast, and northeast China, compared with the Ki and We schemes, ranging from 5 K to 50 K. The averaged RMSDs are 107 K, 81 K, 22 K, and 18 K for Ja, Wi, Ki, and We, respectively. The SAA shows a middle performance, where

(2)

Moreover, the average RMSD for the Ch SRP is lower than that for the We SRP (see Fig. 4). These results indicate that higher weights (see Fig. 3) are strongly associated with lower RMSDs (see Fig. 4) and higher correlation coefficients (see Fig. 5). It indicates that the BMA weights can reflect the predictive skill of the model.

The values of the RMSD for the BMA prediction are lower than those obtained from any of the single model simulation and the SAA method. Another substantial improvement is that the correlation coefficients between the BMA-based reconstructions and the AMSR-E observations are all higher than 0.5 except for the Tibetan Plateau and the wetland in northeast China (see Fig. 5). The average RMSD and correlation coefficient for the BMA scheme in China are 4 K (see Fig. 6(a)) and 0.64 (see Fig. 6(b)), respectively. Thus, the uncertainty of the individual models in the model ensemble is substantially identified through the BMA weights and reduced by the BMA scheme, and the BMA-based reconstruction can generally reveal the spatial-temporal variation of the AMSR-E observations and better reproduce the T_b fields than the SAA method.

4.3 Effects of vegetation types

Because different schemes differ substantially in their representations of the vegetation optical thickness, the differences in the model-predicted T_b induced by the vegetation types can be quite large. The averaged RMSDs and the correlation coefficients of the ensemble simulations from the CLM4.5-CMEM against the AMSR-E observations are presented in Figs. 7 and 8, respectively, where Avg. presents the averaged value of the same type vegetation, and Std. presents the standard deviation.

In general, the CLM4.5-CMEM simulations have a lower D_RMSD and a higher r over the bare soil (see Fig. 7(a)) than those over low vegetation types (see Figs. 7(e)−7(g)). The largest D_RMSD is observed over high vegetation types (i.e., forest, see Figs. 7(b)−7(d)). The SDCP and SRP have larger ranges than the VOP. This suggests that the VOP make more contributions to the T_b simulations than the SDCP and SRP over the vegetation-covered areas, which confirms the results of Fig. 2. Compared with the SAA method, the BMA approach significantly improves the T_b predictions for all vegetation classes, showing lower DRMSD (see Fig. 7) and higher r (see Fig. 8). The RMSDs between the BMA-based reconstructions and the AMSR-E observations range from 2.6 K to 4.9 K, while the SAA ranges from 16 K to 100 K (see Fig. 7). It should be noted that the BMA-derived estimations show a high RMSD (4.8 K) over the bare soil, because the CMEM simulations do not consider the effect of desert and thus have large biases over the desert areas (e.g., the Taklimakan Desert, see Fig. 4). Besides large RMSD, the Ja VOP has negative correlations with the AMSR-E observations for the vegetation conditions (see Fig. 8). It leads to a relatively poorer performance for the SAA predictions (0.20-0.65) compared with the Ki and We schemes (0.46-0.70). In contrast, the BMA multi-model ensemble predictions have the strongest correlations for all vegetation classes, ranging from 0.52 to 0.79. This suggests that the BMA scheme for the CLM4.5-CMEM can generate more skillful T_b fields and agrees better with the observations, with lower RMSD and higher correlation coefficient, than the best individual predictions from the ensemble simulations for all vegetation types.

5 Conclusions

In this study, the BMA scheme is used to improve the T_b predictions of the CMEM coupled with the land surface model CLM4.5 through a weighting strategy based on a multi-model ensemble. The ensemble is generated from four VOPs, three SDCPs, and two SRPs. The results show that the BMA scheme can improve the accuracy of the T_b predictions at the C-band, and can reproduce the spatial variation of the AMSR-E observations with a spatial correlation coefficient more than 0.99 better than the simulations from any individual member or the SAA.

The results show that the VOPs can make more contributions to the T_b simulations over the vegetation-covered areas than the SRPs and SDCPs. The SAA method can reduce the RMSD (50K) compared with the Ja (107K) and Wi (81K) VOPs, but show worse performance than the Ki (22K) VOP and the We (18K) SRP. In contrast, the BMA multi-model ensemble predictions show the lowest RMSD (4K) and the strongest correlation (r = 0.64) for all vegetation classes.

Despite the improved accuracy of the T_b estimations compared with the individual model simulations, the results also show that the BMA predictions have large errors (5K-10K) over the deserted areas (i.e., the Taklimakan Desert) in northwest China, and are weakly correlated with the observed T_b values over the Tibetan Plateau (r < 0.4). This suggests that there is probably room for the improvement in the performances of the CMEM over these areas.

Although the AMSR-E observations are tested at the C-band in this study, the BMA scheme can also be extended to other frequencies (e.g., the L-band or X-band) and satellite observations, e.g., the soil moisture and ocean salinity and soil moisture active-passive satellites. Moreover, the improved T_b predictions of the RTM using the BMA scheme in this study will facilitate the T_b assimilation into the land surface model for improving the surface and root zone soil moisture simulations at regional and global scales.

Acknowledgements The AMSR-E data were downloaded from the National Snow and Ice Data Center (NSIDC, http://www.nsidc.org/data/amsre/index.html). We would like to thank the anonymous reviewers for their comments and suggestions about this paper.