Calibration for Improving the Medium‐Range Soil Forecast over Central Tibet: Effects of Objective Metrics’ Diversity

1. Introduction

Soil moisture (SM) and soil temperature (ST) are crucial variables modulating land-atmosphere fluxes [1,2,3,4]. However, due to the complexity of ST modelling over semi-arid regions, the ST simulations directly produced by land surface model (LSM) exhibit spatiotemporal deficiencies, posing challenges in their regional weather and climate applications [5,6]. Research efforts in improving ST simulations have suggested that the manually corrected high-sensitivity land parameters could benefit the greater scale ST modeling physics [7,8], and the auto-calibrated LSM’s parameter table could benefit the joint SM-ST modeling configuration [9]. Given the great challenges in solving the highly non-linearity in joint SM-ST modeling, e.g., the high-dimensional land parameters and nonlinear physics, the composited objectives evaluating the distance between simulations and observations is proposed to enhance calibration performance (e.g., Kling-Gupta efficiency) [10], whose various internal multi-metrics’ credits need more endeavor to meet with a robust real-world application [11,12]. Therefore, evaluating the effects of the objective metrics’ diversity on calibration performance in solving the spatial complexities of surface simulations is of great significance for improving ST modeling and forecasting over semiarid land.

The LSM parameters optimization or identification has been evolving for decades with the regional application of auto-calibration techniques, primarily achieved by utilizing global search algorithms (GSA, e.g., particle swarm optimization and shuffled complex evolution; PSO [13,14,15,16,17,18] and SCE [19,20,21,22]) to seek optima against specific model objective. As LSM parameter number usually decreases GSA’ efficiency and effectiveness, especially in high-dimensional cases, early research advocated for dimensionality reduction through generalized land parameter sensitivity analysis, such as focusing on reducing insensitive parameters for specific objectives, based on land surface parameter categorization (e.g., soil, vegetation, general, and initial types) to enhance calibration optima [23,24,25]. Moreover, given the intensifying diversity of land surface model applications (e.g., runoff, fluxes), globally applicable land surface parameter estimation has garnered great attention [26,27,28,29], but this has been challenged by the largely varied sensitivities of the distinguished LSM parameters in arid and semi-arid land (ASAL) [30,31].

The well-known Noah LSM though has been widely employed in finer numerical studies [32], but faces increasingly prominent issues related to the representativeness of parameters in complex ASAL applications, such as varying sensitivities of vegetation and general parameters to thermal flux respectively [25,33,34]. This poses continuous challenges for refined land surface applications in Tibet, a region with diverse climatic zones, e.g., LSM parameter diverse advantages in different regions of a similar surface [18,28,29]. Despite the establishment of a refined SM-ST observation network under a semiarid climate over central east Tibet (i.e., northwest Naqu) [35], which features grassland as the primary land cover and clay as the dominant soil texture, more comprehensive calibration objective against the LSM parameter uncertainties reduction and surface enhancement is still required for the robust ST modeling [9,36,37,38].

In fact, with the development of land remote sensing, given the diversity of GSA’ strategies and application objectives (such as SM, ST, runoff, and fluxes), the objective metric designs of auto-calibration have greatly developed to enhance LSM modelling performances. For instance, LSM calibration using multi-source remote sensing data, the multi-objective design concentrates on the comprehensive inversion characteristics of remote sensing SM and/or ST observations are essential [39,40,41]. Similarly, in calibrating applications aimed at improving the spatial accuracy of surface state predictions, a multi-objective design that considers horizontal variations [42,43,44,45] and/or vertical stratification [46,47,48,49,50,51,52,53,54,55,56,57,58,59] of states and observations is crucial. Generally, the multi-objective metrics can, to a certain extent, address the issues of observational data fusion and multi-state complex error measurement in specific calibration applications, emphasizing the enhanced role of spatial dimensions of single or multiple land surface state errors as holistic objectives.

Moreover, as the inherent scale uncertainties of land surface state (e.g., the distance between simulation and observation could fall into the non-Euclidean space) lead to challenges in assessing and ranking LSM’ performances in high-dimensional searching space, the holistic objectives with differentiated metrics have been widely developed to simplify and enhance the calibration [50,51,52]. And they can be primarily categorized into flexible (e.g., Pareto front [53,54,55,56], and dominated Pareto [57,58,59]) and deterministic approaches. The Pareto front adjusts the cumulative distribution of metrics based on external algorithm storage and aims at general LSM modelling with globally applicable parameters, which can be an expensive evaluator that independent from GSA. While within the GSA’s evaluator, the dominated Pareto compares relationships among metrics and the deterministic approach combines various metrics, and then aims at determining the optimally combined solution of various simulations.

However, though metrics offer the deterministic reference for estimating diversity of model performances, research indicates that the application of computational methods employed by these metrics often exhibits blindness against same datasets. For instance, integrated metrics such as Nash–Sutcliffe efficiency and Kling–Gupta efficiency can be utilized for algorithm comparison, yet for actual model evaluation, direct metrics are still necessary to indicate [60,61,62,63]. The optimal applicability of direct metrics like root mean square error and mean absolute error in describing data is premised on their distributions conforming to normal and Laplacian distributions, respectively [64,65]. Furthermore, correlated coefficient (CC) is susceptible to the monotonicity and nonlinearity of two types of data [66,67]. Consequently, the performance of these metrics, when combined across different dimensions, often necessitates a comprehensive evaluation of their calibration suitability tailored to varied spatiotemporal requirements of LSM simulations [68,69].

Overall, due to the high short-term surface simulation biases and the high dimensionality of parameter space of Noah LSM, there is a lack of comparisons for the multi-objective metric methods targeted at SM-ST joint calibration to validate their ability in reducing spatial complexities of regional land simulations over ASAL. Therefore, to fill this gap, this study, based on the Noah LSM and its full parameter table, utilizes the PSO/SCE method and various SM-ST objective metrics, combined with intensive regional soil site observations, to explore the impacts of metric objective differences on the spatial heterogeneity and uncertainty of regional land surface parameters, calibration efficiency and effectiveness, as well as temporal and spatial errors in surface forecasting. And suggestions are provided for the regional ST modeling configuration, aiming to improve medium-range ST forecasting of semi-arid regions.

2. Methods

2.1. Calibration Schemes

2.1.1. Evolution Algorithms

Group and individual social behaviors are incorporated into the core PSO algorithm process [13,14,15,16,17,18]. The algorithm first randomly selects the scaled (or normalized) parameters ($\mathrm{x}$) to generate the initial population including the individual position and the speed of position change, i.e., $\left({\mathrm{x}}_{\mathrm{i}}^0,{\mathrm{v}}_{\mathrm{i}}^0\right),\mathrm{i}\in (1,\dots ,\mathrm{n}\mathrm{p})$, and further obtain the local and global optimal position (${\mathrm{P}}_{\mathrm{b}}^0$ and ${\mathrm{g}}_{\mathrm{b}}^0$) through evaluating and sorting, where the superscript represents the time slice (Figure 1). Then proceed with following steps repeatedly till the stop criteria is met. Comparing individuals’ current and previous evaluation values to obtain the current local optima (${\mathrm{P}}_{\mathrm{b}}^{\mathrm{t}}$); then sorting local optima to find the current global optima (${\mathrm{g}}_{\mathrm{b}}^{\mathrm{t}}$); each particle’s speed and position are updated using ${\mathrm{v}}_{\mathrm{i}}^{\mathrm{t}+1}=[\mathrm{w}{\mathrm{v}}_{\mathrm{i}}^{\mathrm{t}}+{\mathrm{c}}_1{\mathrm{r}}_1({\mathrm{p}}_{{\mathrm{b}}_{\mathrm{i}}}^{\mathrm{t}}-{\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}})+{\mathrm{c}}_2{\mathrm{r}}_2({\mathrm{g}}_{\mathrm{b}}^{\mathrm{t}}-{\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}}\left)\right]\times \sqrt{\mathrm{n}}{\mathrm{v}}_{\mathrm{m}}\times \left[\right({\mathrm{r}}_3-1){\mathrm{v}}_{\mathrm{r}}+1]+{\mathrm{c}}_3{\mathrm{v}}_{\mathrm{r}}$ and ${\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}+1}={\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}}+{\mathrm{v}}_{\mathrm{i}}^{\mathrm{t}+1}$ respectively, where ${\mathrm{v}}_{\mathrm{m}}$ and ${\mathrm{v}}_{\mathrm{r}}$ equal to 0.5 and 0.15, ${\mathrm{r}}_1$ and ${\mathrm{r}}_2$ equal to 0.5 and 0.15, $\mathrm{w}$ equal to 0.9, ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, and ${\mathrm{c}}_3$ equal to 2.0, 2.0, and 10-7 [9]. Note that except $\mathrm{n}\mathrm{p}$, all the other parameters that can affect the generality of PSO were vaguely related to the dimension of the parameter space.

Community and induvial social behaviors are incorporated into the core SCE algorithm process [19,20,21,22]. The algorithm first randomly selects the scaled (or normalized) parameters ($\mathrm{x}$) to generate the initial population, then further obtain initial local and global optima (${\mathrm{D}}^0$ and ${\mathrm{g}}_{\mathrm{b}}^0$) through evaluating and sorting with marked orders. Then proceed with following steps repeatedly till the stop criteria is met. Evaluating and marking the individuals to reorganize the original population into $\mathrm{n}\mathrm{c}$ communities (each community has $\mathrm{m}$ points) ; through complex competitive evaluation (CCE) of each point where the triangular probability distribution ${\mathrm{P}}_{\mathrm{k}}={\displaystyle \frac{2\times \left(\mathrm{m}+1-\mathrm{k}\right)}{\mathrm{m}\times \left(\mathrm{m}+1\right)}},\mathrm{k}=1,\dots ,\mathrm{m}$, to determine the previous generation, and the new individuals from each community are mixed to form a new population; then reorder the individuals to form $\mathrm{n}\mathrm{c}$ new communities and obtain the current local and global optima (${\mathrm{D}}^{\mathrm{l}}$ and ${\mathrm{g}}_{\mathrm{b}}^{\mathrm{t}}$) (Figure 1). Also, $\mathrm{n}\mathrm{c}$ and $\mathrm{m}$ equal to 2 and $2\mathrm{n}+1$, while the outer cycling number ($\mathrm{n}\mathrm{e}$), the internal and external iteration number of CCE ($\mathsf{\alpha }$ and $\mathsf{\beta }$) equal to $\mathrm{n}+1$, $\mathrm{n}\mathrm{c}$, and $\mathrm{m}$, respectively. Note that except $\mathsf{\alpha }$ and $\mathsf{\beta }$, all five parameters that can affect the generality of SCE were only related to the dimension of the parameter space.

Note that $\mathrm{x}$ is scaled with equation as ${\displaystyle \frac{{\mathrm{x}}_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}-{\mathrm{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{\mathrm{x}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}$ for both PSO and SCE. The totals individual number ($\mathrm{n}\mathrm{c}\times \mathrm{m}$) of SCE and total particle number ($\mathrm{n}\mathrm{p}$) of PSO both equal to $2\left(2\mathrm{n}+1\right)$. And for both SCE and PSO, the GSA stops when the evaluation contour ($\mathrm{i}\mathrm{e}$) is greater than 105 Noah runs [9]. The equitable population size and stop criteria intend to reduce that the objective metrics’ impact could be affected by the algorithms themselves and ensure the relatively equitable investigation.

2.2.2. Optional Evaluator

The evaluator of the above-mentioned GSA algorithms used in our study is shown in Figure 2, which includes a fixed physical constraint and an optional objective function. The physical constraint formula (${\mathrm{f}}_{\mathrm{c}}$) represents the soil moisture of the first two surface soil layers ($\mathrm{S}\mathrm{M}\mathrm{C}1$ and $\mathrm{S}\mathrm{M}\mathrm{C}2$) only varies between the wilting point ($\mathrm{W}\mathrm{L}\mathrm{T}\mathrm{S}\mathrm{M}\mathrm{C}$) and the soil moisture where transpiration stress begins ($\mathrm{R}\mathrm{E}\mathrm{F}\mathrm{S}\mathrm{M}\mathrm{C}$) [9,24,39].

Under this constraint ($\mathrm{r}<0$), the unscaled parameters will drive the target model to run simulation once for evaluation based on the corresponding objective value, which can be selected based on the objective type ($\mathrm{Z}$), and $\mathrm{Z}$ is a constant that varies between 1 and 8. Note that once $\mathrm{t}$ is determined, the predefined corresponding objective metric or function (${\mathrm{f}}_{\mathrm{o}}$) that measures the distances between simulation ($\mathrm{s}$) and observation ($\mathrm{o}$) is also determined at the very beginning.

2.2. Composited Metrics

For calibration schemes based on GSA, the parameter simulation problem in LSM is addressed by searching for the optimal parameters and/or simulations that minimizes or maximums the objective function (${\mathrm{f}}_{\mathrm{o}}$). Especially, ${\mathrm{f}}_{\mathrm{o}}$ has been extended into multi-dimensions, i.e., layers ($\mathrm{n}\mathrm{l}$) and variables $\left(\mathrm{n}\mathrm{e}\right)$, to meet the multiple dimensional SM and ST objectives, and eight different metrics are investigated during present study (Table 1).

This study mainly examines fixed metrics composed of fundamental measures such as linear correlation coefficient ($\mathrm{c}\mathrm{c}$) [66], Kling-Gupta efficiency ($\mathrm{k}\mathrm{g}\mathrm{e}$) [10], absolute error ($\mathrm{a}\mathrm{e}$) [65], Nash-Sutcliffe efficiency ($\mathrm{n}\mathrm{s}\mathrm{e}$) [63], and root mean square error ($\mathrm{r}\mathrm{m}\mathrm{s}\mathrm{e}$) [64] across different dimensions. Specifically, correlation coefficients (CCS), enhanced Kling-Gupta efficiencies (EKGE), mean absolute errors (MAES), Nash-Sutcliffe efficiencies (NSES), and root mean square errors (RMSES) represent the average of $\mathrm{c}\mathrm{c}$, $\mathrm{k}\mathrm{g}\mathrm{e}$, $\mathrm{a}\mathrm{e}$, $\mathrm{n}\mathrm{s}\mathrm{e}$, and $\mathrm{r}\mathrm{m}\mathrm{s}\mathrm{e}$ across both the variable and layer dimensions. Additionally, the enhanced multiple objectives (EMO) integrate the average values of the measure that combines $\mathrm{c}\mathrm{c}$, $\mathrm{a}\mathrm{e}$, $\mathrm{n}\mathrm{s}\mathrm{e}$, and $\mathrm{r}\mathrm{m}\mathrm{s}\mathrm{e}$ in the variable and layer dimensions.

Furthermore, since the surface variations of the topsoil layer could often determine the sublayer’s variations according to infiltration [32], therefore, the dimensional target of the upper layer is assumed to be the dominated Pareto solution, and the top layer’s objective that is larger (or smaller) than sublayers’ is taken as the current optimally maximum (or minimum) solution [57—59]. Consequently, the Pareto-dominant KGE (PKGE) and the Pareto-dominant multi objectives (PMO) indicates the dominated top layer’s value of EKGE and EMO respectively.

All the above-mentioned multi-objective metrics’ variable and layer dimensional number are 2 (e.g., SM and ST) and 4 respectively. CCS varies in [－1, 1]. EKGE and NSES both vary in (－∞, 1], EMO, MAES, and RMSES all vary in [0, ＋∞). PKGE varies in (－∞, 1], and PMO varies in [0, ＋∞). Therefore, the value of the metric determines the performance of the evaluator, and the direction of the metric determines the direction of the search, that is, continuously approaching the optimal value of the metric (i.e. the final ideal termination condition) towards the calibrated optimal solution.

2.3. Performance Evaluation

2.3.1. Parameter

During this study, parameter heterogeneity was defined as variations or sensitivities of land parameters across sites. Due to the immense dimensionality of parameter—site sensitivities, parameter relative sensitivities based on the two predefined limits of the parameter space are suggested, e.g., if more (fewer) sites met (failed to meet) a parameter's limit compared to others, indicating sites' relative sensitivity to that parameter within the limit's confidence [9,24]. Since the parameter relative sensitivities (or heterogeneity) are usually large while their homogeneity could be small (and thus be easily observed), here to qualify this and simplify metrics’ diversity investigation, we further propose the parameter numbers with low site sensitivities as homogeneity (H). Consequently, low H (>0) of this study indicates high heterogeneity of one metric quantitatively. Note that when all and no sites cross the parameter’s limits, H equals 0 and $\tilde{0}$ respectively.

The parameter’ spatial uncertainty is defined as the land parameter range and outlier against the sites, e.g., one parameter’s inter quartile range (IQR, >0), smaller parameter ranges and outlier numbers indicated fewer uncertainties and fewer unaccountable factors respectively [9]. Especially, to simplify different metrics’ effects on parameter uncertainty, the whole parameter space’s uncertainty is defined as the inter quartile range of all the IQR ensembles of different parameters (or IQRD) in the parameter space. Consequently, the IQRD’s inter quartile range and outlier indicate the quantitative parameter—uncertainties among metrics.

Especially, compared to SCE, the parameter number with less parameter uncertainties (PNL, >0) and the outlier number reduction of parameter uncertainties (ONR) in PSO are summarized in this study to qualify the justice that if the metrics’ parameter uncertainty is affected by GSA itself or not. As the heterogeneity and uncertainty differences of different metrics could account for the metric-informed method’s performance in solving parameter spatial complexities during SM-ST calibration, thereby, the metric with less parameter uncertainties and heterogeneity could meet the preferable LSM configuration demand in surface forecasting.

2.3.2. Objective

As the population position of on generation, e.g., the best (P_b) or medium (P_m, if non solution is met) locations that known as fitness against the number of LSM runs (or the convergence speed), could indicate the method’s performance in calibration efficiency, therefore, the better fitness values (e.g., larger EKGE values or smaller EMO values) with fewer LSM runs indicate more efficiency, where the success rate exploring evolution abilities is usually put alongside with.

Moreover, as the optimal objectives (e.g., the final EKGE or EMO values) could indicate method’s performance in calibration effectiveness, the larger or smaller optimal objectives that depend on the direction of predefined metrics indicate more effectiveness. Furthermore, since the kernel density distribution of optimal values across different sites demonstrates their spatial enrichment characteristics, the variation in enrichment between different algorithms (such as PSO or SCE) to a certain extent reflects their capacity to address the spatial disparity in SM-ST simulation.

2.3.3. Simulation

To simplify the spatial complexity among regional datasets, linear fitting between the observations (OBS) and simulations (SIM) for all sites is conducted [68]. The linear fitting’s slope (s) demonstrates the sensitivity of SIM to OBS, while its coefficient of determination (r2) or the goodness of fit demonstrates if the sensitivity or linear model is robust or not. Moreover, under the assumption of the normal distribution of the errors between SIM and OBS (${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$) of all sites, Gaussian fit of ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ that are resampled with 100 bins 100 is conducted to generate at most two signals determining the main distribution characteristic., e.g., the amplitude (or frequency, f) and center (c) [69]. Here the compound feature of f and c that is closer to the normal distribution indicates the better performance or more consistent with the assumption.

The method’s performance in optimal simulation and forecast is qualified using the spatial differences and similarities of surface conditions among different datasets, e.g., ST and SM simulations or reanalysis, and observations, by the following equation:

$${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}=\sqrt{{\displaystyle \frac{\sum _{\mathrm{j}=1}^{\mathrm{n}\mathrm{s}}{({\mathrm{s}}_{\mathrm{j}}-{\mathrm{o}}_{\mathrm{j}})}^2}{\mathrm{n}\mathrm{s}}}},{\mathrm{C}\mathrm{C}}_{\mathrm{S}}={\displaystyle \frac{\sum _{\mathrm{j}=1}^{\mathrm{n}\mathrm{s}}\left[\right({\mathrm{s}}_{\mathrm{j}}-\overline{{\mathrm{s}}_{\mathrm{n}\mathrm{s}}}\left)\right({\mathrm{o}}_{\mathrm{j}}-\overline{{\mathrm{o}}_{\mathrm{n}\mathrm{s}}})]}{\sum _{\mathrm{j}=1}^{\mathrm{n}\mathrm{s}}{({\mathrm{s}}_{\mathrm{j}}-\overline{{\mathrm{s}}_{\mathrm{n}\mathrm{s}}})}^2\sum _{\mathrm{j}=1}^{\mathrm{n}\mathrm{s}}{({\mathrm{o}}_{\mathrm{j}}-\overline{{\mathrm{o}}_{\mathrm{n}\mathrm{s}}})}^2}})$$

(1)

where $\mathrm{i}$ and $\mathrm{j}$ represent the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ time and the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ site, respectively, and $\mathrm{n}\mathrm{s}$ represent for the total number of stations. And smaller ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and/or high ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ indicate better performance.

Meanwhile, the Taylor diagram [60,61,62,63] that could assemble the comprehensive statistics (i.e., standard deviation, root-mean-square difference, and correlation) in a temporal sequence between SIM and OBS was also created for comparison with the method’s skills. Usually, a smaller distance away from the reference location (OBS) indicated more skills. Note that the SIM datasets (30 min) were linearly interpolated into 3 h for a broad comparison with the land reanalysis.

In addition, aside from the uncertainties and heterogeneous requirements in surface prediction parameters (manifested as variations in calibration performance), addressing the precision demands inherent to surface prediction (evidenced by differences in calibration robustness), this study employs indicators such as KGE increment, ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reduction, and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ increment to clarify the performance of various SM-ST objectives in the parameter—simulation and/or calibration—prediction frameworks, aiming to explore the target metrics for both the optimal configuration of LSM and the maximum benefits of surface predictions.

3. Experiments

3.1. Model and Data

The Unified Noah LSM is created to better predict the effects of land surface processes on regional weather, climate, and hydrology. It is intended to comprehend the intricate biophysical, hydrological, and biogeochemical interactions between the land surface and the atmosphere at micro- and mesoscales (Figure 3A) [32]. The simple driver Noah LSM (version 3.4.1, https://ral.ucar.edu/model/unified-noah-lsm, accessed on 31 July 2024) has been recently extended into the muti-point applications over central Tibet [9].

The SM-ST observations that are firstly derived from the highest altitude soil moisture network in the world (Figure 3B, whose elevations are above 4470 m), which is constructed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS) with four soil depths (i.e., 0–5, 10, 20, and 40 cm) [35], and are further assembled into the multi-site (i.e., 12) observations of the local warm season (i.e., covering from 1 April to 31 July 2014) over northwest Naqu city that has a typical semiarid climate by using simple quality control based time continuity correction (detailed described in Ref [9]). Also, the global land data assimilation system (GLDAS) [70] grid soil reanalysis data with resolutions of 3 h/0.25°is collected for broader comparison with the surface simulations during this study.

The gridded meteorological surface datasets that merging a variety of data sources are firstly developed by ITPCAS, with a 3 h interval (3 h) and a resolution of 0.1° × 0.1°, were produced by [71], and are further reassembled into the multi-site LSM forcing dataset by using the inverse distance-weighted quadratic spline interpolation method to drive the Noah LSM.

According to the observational soil and surface characteristics, the multi-site Noah LSM is configured with a 4-layer depth and 30-min runtime step, and the soil and vegetation types are mainly silt and grassland, while the slop type is assumed to be flat (e.g., 1). And the forcing time step (3 h) and screen height (10 m for winds and 2 m for temperature) for the LSM are the same as the input forcing data [9].

3.2. Experimental Description

Three month long warm-up run (covering the period from April 1 to July 1 of 2014) of the multi-site LSM, that initialized with the unobserved default parameters (i.e., the “General”, “Vegetation”, “Soil”, and partial “Initial” types) [32] and partially observational “Initial” parameters (i.e. SMC1-4 and STC1-4), is firstly conducted to obtain the default multi-site parameter tables including spatially distinguished “Initial” parameters for the following experimental runs [9]. Based on this, one-month long run ranged from July 1 to July 31 (or the control run briefed as CTR hereafter) is conducted as the referenced surface conditions resulted from the default LSM parameter table configuration.

The multi-objective metrics varied calibration runs that ranged from July 1 to July 15 are conducted to obtain the calibrated multi-site LSM models with metrics informed parameter tables and further investigate the metrics’ impact on calibration’s abilities in solving the spatial complexities of Noah LSM. Then the abovementioned various objective informed LSM models run from July 15 to July 31 to obtain the hopefully improved surface forecasts and further investigate metrics’ impact on surface forecast.

Therefore, the difference between CTR and calibration could account for the calibration performance, and the difference among different calibration runs could account for the metric’s impact on the calibration. Meanwhile, the difference between CTR and calibrated forecast runs could account for the calibrated models’ performances, and the difference among different calibrated forecast runs should account for the metric’s impact on surface forecast. Note that all objective metrics within both PSO and SCE algorithms are conducted to explore if the potentially improved surface forecast could be highly affected by the calibration algorithms themselves or not.

4. Results

4.1. Case Perspective

As land surface models utilize parameters and forcing inputs to prepare land surface forecasts, the issues of surface simulation and local application in typical semi-arid regions (i.e., rapidly applying calibrated parameters to surface forecasts) are exemplified here. To this end, a review of the spatiotemporal characteristics of the default forcing, initial parameters, and their overall simulation status across different periods, including control, simulation calibration, and forecast verification, is conducted to clarify the fundamental manifestations of the issues involved in this study.

4.1.1. Model Configure

The site averaged 3-h meteorological forcing values against time during the study period are shown in Figure 4 a~d. During July 2014, the diurnal variation in temperature (${\mathrm{T}}_{2\mathrm{m}}$) mostly ranged between 5 and 15 °C, with an extremely dry atmosphere whose relative humidity (${\mathrm{R}\mathrm{H}}_{2\mathrm{m}}$) values were mostly below 1%. The relatively low wind speed (${\mathrm{W}\mathrm{S}}_{10\mathrm{m}}$) generally varied between 0 and 6 m s^–1, and the wind direction (${\mathrm{W}\mathrm{D}}_{10\mathrm{m}}$) was mostly dominated by southern flow (between 180° and 270°) from July 1 to July 10 and from July 16 to July 21, respectively, but the opposite in other periods. The incoming shortwave radiation ($\mathrm{S}\mathrm{W}$) exhibited strong diurnal variation between 0 and 600 W m^–2, and the incoming longwave radiation (L$\mathrm{W}$) varied between 250 and 350 W m^–2. The pressure ($\mathrm{P}$) was generally around 586 hPa and the maximum hourly precipitation (${\mathrm{R}}_{1\mathrm{h}}$) was about 5 mm h^-1 on July 10.

The initial land parameters of all the sites (or default parameters) that needs calibration has shown great variety for the “Initial” types (Figure 4 e), and this is especially pronounced for the moisture-related parameters, i.e., the SMC and the SH2O. This should attribute to the differences in pre-experiment 3-month simulation of different sites. Furthermore, due to the lack of direct observations for other types of land surface parameters, they are configured using statistical data from optimization experiments based on limited benchmarks from the previous study (serving only as reference inputs for control experiments, consistent with conventional numerical model configurations) [9]. Therefore, in numerical operations, parameter variations among different stations mainly exist in the initial types. Consideration of multi-site calibration can account for the differences among stations with unobserved parameters under existing observational constraints, namely, the spatial heterogeneity and uncertainty of the parameter space. Consequently, the spatial heterogeneity of parameters (i.e., the sensitivity of commonly used parameters to different stations or the number of intersections of the same parameters at different stations in the parameter space) and the characteristics of uncertainty (such as the inter quartile range and outlier features of parameters at different stations) in relation to the differences in various optimization objectives are the key areas of focus for further investigation in this study.

4.1.2. Forecast Problem

The CTR simulations and observation datasets (OBS) for the surface layer are compared in Figure 4 e~f. For the whole experimental period, the linear fit for the surface soil moisture (${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$) exhibited a small increasing slope (about 0.21) with weak consistency, and the surface soil temperature (${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$) had a larger decreasing slope (about –0.45) with strong differences. Moreover, the linear fits of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ for the calibration and forecast periods were 0.22 and 0.15 respectively, and the linear fits of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ for calibration and forecast periods were -0.48 and -0.4 respectively. This indicates that the surface conditions of the forecast period were slightly better than those for calibration period. Generally, ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ fits better than ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$.

In addition, the Gaussian fits of the errors between ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ observation and simulation (${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$) for whole experimental period had a sharp and narrow distribution, which was centered around 0.15 ${\mathrm{m}}^3·{\mathrm{m}}^{-3}$ with a frequency of around 800, while the ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ distributions of the calibration and forecast periods had centered around 0.16 and 0.13 ${\mathrm{m}}^3·{\mathrm{m}}^{-3}$, with the frequency of around 500 and 300 respectively. This indicates ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ were mostly underestimated for all periods and this is more pronounced at the calibration period. Nevertheless, the ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ for different periods had shown bimodal distributions (Figure 4i), whose centers were located around -4 and 9 K (whole period), -5 and 9 K (calibration period), and -4 and 8 K (forecast period) respectively. This indicates ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ were both under- and over-estimated, and the later were more pronounced. Generally, ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ were both underestimated.

In general, though ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ in CTR exhibited better consistency with OBS than ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, the overall surface simulation underestimation of Noah LSM could be great for regional surface forecast applications. Note that either the ITPCAS forcing data sets or the improved heat-sensitive parameter Z0h (also known as CZIL) to improve ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ with the Noah LSM over a surface near our study area [7], and the non-negligible biased ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ and the spatially diversed parameter space in CTR indicated a more effective calibration in present study. Since multi-objective calibration can reduce these spatiotemporal errors through parameter identification to improve subsequent forecasts [9], the next focus is on how different target metrics affect the performance of calibration and forecasting.

4.2. Effects on Calibration

4.2.1. Optimal Parameters

Due to the significant spatial heterogeneity exhibited by most optimal parameters in PSO and SCE, this study conducted a statistical analysis of spatially homogeneous (or heterogeneous) parameters based on parameter type classifications (see Figure S1-1, briefed in Table 2). For the "Vegetation" type, except for the SCE scenario considering CCS, the number of homogeneity parameters in other scenarios is zero, indicating heterogeneity. Regarding the "Soil" type, the counts (H_p, H_s) of homogeneity parameters for PSO and SCE calibration schemes based on EKGE, EMO, MAES, and RMSES metrics are (1, 3), (2, 2), (1, 1), and (2, 1), respectively. For the "General" type, the (l1, l2) based on CCS, EKGE, EMO, MAES, and RMSES metrics are (1, A), (2, 2), (2, 2), (1, 1), and (1, 1), respectively. For the "Initial" type, the (l1, l2) based on EKGE, EMO, and MAES metrics are (4, 2), (3, 2), and (2, 1), respectively. Evidently, among all pairs, the spatial homogeneity of optimal parameters for all "Vegetation" types in PSO and SCE is relatively minimal, suggesting the strongest heterogeneity. Conversely, "Soil" and "General" types exhibit minimal spatial heterogeneity, while "Initial" types fall in the middle. Notably, QTZ and SBETA parameters consistently demonstrate homogeneity, below the parameter space threshold (0.03), across PSO and SCE schemes based on EKGE, EMO, MAES, and RMSES metrics.

Regarding the counts of homogeneity parameters in PSO and SCE schemes, when considering the disparities among metrics, we observe the following: for CCS, the counts are 1 and 40, respectively; for EKGE, both schemes yield 7; for EMO, the counts are 7 and 6; for MAES, 4 and 3; for NSES, 2 and none; both PKGE and PMO register none; and for RMSES, the counts are 5 and 2. Evidently, there exist substantial variations in the homogeneity or heterogeneity of parameters among calibration schemes based on different metrics. Notably, CCS exhibits the lowest parameter heterogeneity, followed by EKGE, then EMO, and subsequently MAES and RMSES. NSES displays relatively poor parameter heterogeneity, whereas PKGE and PMO manifest the highest degree of parameter heterogeneity.

In addition, the inter quartile ranges (IQR) of various parameters and the entire parameter space of PSO and SCE contributed by different metrics are shown in Figure 5. For PSO, the maximum IQR, which has about 4.8 of the parameter SNP in the “Vegetation” type, had the largest uncertainties, while the EMO made the largest contribution. However, the IQR, which is about 1.2 of the SBETA parameter in the “General” type, behaves oppositely, while EKGE, EMO and RMSES make the smallest contributions (Figure 5a). For SCE, the maximum IQR around 1.82 of the CZIL parameter in the “General” type has the largest uncertainties, while EMO has the largest contribution. Nevertheless, the IQR that is around 0.61 of the parameter CSOIL in the “General” type behaves conversely, while EKGE, EMO and RMSES make the smallest contributions (Figure 5b). In general, PSOs have achieved higher IQRs than SCEs on most metrics. And PSO and SCE achieved the lowest uncertainties of the parameters SBETA and CSOIL in the “General” type.

The Inter quartile Range (IQR) distribution of the global optimal parameter space for the PSO schemes across various metrics exhibits a broader and more scattered range compared to that of the SCE (Figure 5c). For PSO, the median sizes of the IQR distributions (IQRD) within the global optimal parameter space, ranked from highest to lowest, are PKGE > PMO > NSES > CCS > EMO > EKGE > RMSES > MAES. In contrast, for SCE, the order is PMO > PKGE > EKGE > EMO > MAES > RMSES > CCS. Furthermore, for PSO, the number of outliers in the IQRD is highest for EKGE with 3, followed by EMO and CCS with 2, while the rest of the metrics have 0 outliers. For SCE, EKGE and RMSES share the highest number of outliers at 2, followed by EKGE and MAES with 1 outlier each, and the rest are 0 (Figure 5d). In summary, significant differences exist in the IQRD of the global optimal parameter spaces across different metrics, with SCE exhibiting smaller IQRD but relatively more outliers. Notably, EKGE and EMO exhibit relatively large numbers of outliers in both PSO and SCE.

Furthermore, the PSO’s parameter spatial IQRs are compared with SCE in different types in Table 3, e.g., the parameter number with less uncertainties (PNL) and the outlier number reduction of parameters’ uncertainties (ONR) in PSO when compared to SCE. For the "Vegetation" type, all metrics are null except for the PNL value of EKGE, which is 2, while the ONR of all metrics is non-positive. For the "Soil" type, the PNL values are positive for all metrics except for CCS, NSES, and PKGE, which are null. The ONR values are positive for EKGE, EMO, PMO, and RMSES, while negative for the rest. Regarding the "General" type, all metrics exhibit positive PNL values except for NSES, PKGE, and PMO, whose PNL values are null. The ONR values are positive for EKGE, EMO, MAES, and RMSES, and negative for the others. In the "Initial" type, only EKGE and EMO have positive PNL values, with the rest being null. The ONR values are positive for all metrics except for CCS, PKGE, and PMO, which are non-positive. In summary, summing the PNL values across types, EKGE has the highest total (8), followed by EMO and RMSES (7), then MAES (5), with PMO and CCS having the lowest totals (1). PKGE has no PNL value. For the ONR values, EMO has the highest total (9), followed by EKGE (3), then RMSES (3), while PMO has the lowest (2). The rest of the metrics have negative ONR values.

In summary, for the SM-ST calibration of the same metric, SCE consistently achieves lower parameter uncertainty than PSO, albeit at the cost of relatively higher spatial heterogeneity. Specifically, in terms of parameter uncertainty, MAES in PSO and CCS in SCE exhibit the smallest metrics. As for parameter spatial heterogeneity, EKGE and EMO in PSO yield the smallest metrics, while SCE solely displays the smallest EKGE.

4.2.2. Effectiveness and Efficiency

Figure 6 shows the different metrics’ fitness (i.e., the best position of one population, Pb) curves of calibration, the median convergency position, and the median converged Noah run numbers of PSO (${\mathrm{C}}_{\mathrm{P}}^{\mathrm{P}},{\mathrm{C}}_{\mathrm{P}}^{\mathrm{N}}$) and SCE (${\mathrm{C}}_{\mathrm{S}}^{\mathrm{P}},{\mathrm{C}}_{\mathrm{S}}^{\mathrm{N}}$) for all sites. For CCS, both PSO and SCE had both sharply increased before 3,000 Noah runs, and both converged to 1 but at around 79,475 and 66,663 runs respectively. For EKGE, PSO and SCE have both sharply increased before 10,000 Noah runs but converge to 0.56 at 99,017 runs and 0.53 at 90,731 runs respectively. For EMO, PSO and SCE both decrease to 1 before 8,000 Noah runs but converge to 1 at 99,297 runs and 1.08 at 82,709 runs respectively. For MAES, PSO and SCE both quickly decrease to the range of 0.7-1.1 before 10,000 Noah runs but converge to 0.79 at 99,765 runs and 0.81 at 94,795 runs respectively. For NSES, PSO and SCE have both instantly reaching 1 at 187 runs, indicating the most rapid convergence among all metrics. However, for PKGE and PMO, since volatile finesses (e.g., who vary within $(-\mathrm{\infty },1]$ and $[1,+\mathrm{\infty })$ respectively) are found for all sites in each generation, nonstrict solutions can be observed. For RMSES, PSO and SCE both sharply decrease to 1 before 5,000 Noah runs but converge to 0.97 at 99,391 runs and 0.98 at 94,029 runs respectively.

Generally, except PKGE and PMO, other metrics of PSO have achieved better effectiveness as indicated their better fitness values, but with relatively worse efficiency as indicated their larger converged runs compared to those of SCE. The non solution performance for the metric PKGE and PMO of both PSO and SCE have indicated their requirements of more Noah runs in achieving convergence, or the potential failure of the Paetro dominated logic (i.e., that surface improvement likely improve the subsurface). For MAES, NSES, and RMSES, fitness curve of site C4 is found to be notably biased from (or worse than) that of other sites. Nevertheless, for all the metrics’ convergences, MAES has the largest range, and this could indicate the divergent convergence domain.

Figure 7 presents the success rate curves for calibration across various metrics. For CCS, PSO experiences a decline from 70% to 20% during the first 10,000 Noah runs, followed by a gradual decrease to near zero. In the case of EKGE, PSO initially shows a decline from 80% within the first 5,000 Noah runs, subsequently exhibiting two distinct patterns: fluctuations around 40% and 20%, respectively. For EMO, PSO drops from 80% to nearly 0% within the initial 25,000 Noah runs, with some stations subsequently exhibiting strong fluctuations between 0% and 80%. MAES follows a similar trend, with PSO declining from 80% to near 0% within the first 15,000 Noah runs, and subsequent intense fluctuations between 0% and 80% at certain stations. For NSES, PSO gradually decreases from 80% to 20% within the first 35,000 Noah runs and remains stable thereafter. PKGE and PMO exhibit similar behavior, with PSO slowly declining from 80% to 20% within the first 20,000 Noah runs and fluctuating slightly around 20% thereafter. SCE's performance in PKGE resembles that of CCS. In contrast, RMSES displays a fluctuating decline from 80% to 0% within the initial 20,000 Noah runs for PSO, followed by drastic fluctuations between 20% and 80%. However, SCE consistently demonstrates a rapid initial decrease from 80% to 20% across nearly all metrics, maintaining this level thereafter.

For all metrics, the search domain of SCE exhibits a consistent pattern, characterized by an L-shaped thin linear region. In contrast, PSO's search domain displays significant fluctuations and notable variations across different metrics (e.g., EKGE, EMO, MAES, RMSES), albeit with an overall larger area than SCE. This suggests that for most metrics, PSO demonstrates stronger evolutionary capabilities compared to SCE, which primarily contributes to PSO's slightly slower convergence rate compared to SCE.

Figure 8 presents the statistical performance of the optimal objectives across all stations for various metrics. For CCS, both PSO and SCE exhibit a concentrated distribution near 1, with PSO displaying a tighter clustering and an outlier at 0.973. In the case of EKGE, PSO and SCE concentrate around 0.58 and 0.53, respectively, with PSO showing a more focused distribution and an outlier at 0.34. For EMO, PSO and SCE are centered near 1 and 1.1, respectively, with PSO displaying a relatively dispersed distribution and an outlier at 1.5. MAES values for PSO and SCE are centered around 0.79 and 0.81, respectively, demonstrating similar distributions. For NSES, PKGE, and PMO, both PSO and SCE have concentrated distributions near 1, with NSES exhibiting a more tightly clustered distribution compared to the other two metrics. Finally, for RMSES, PSO and SCE are centered around 0.9 and 1.1, respectively, with SCE displaying a more focused distribution, and both having outliers at around 2.4.

It is evident that for the optimal solutions of PKGE and PMO, both PSO and SCE yield values of 1, indicating the absence of optimal solutions or the need for more time to locate them. In contrast, numerical optimal solutions were achieved for other metrics. Furthermore, while PSO consistently outperformed SCE in attaining better optimal solutions across almost all metrics, significant variations were observed in the enrichment levels of optimal solutions between PSO and SCE under different metrics. For instance, PSO surpassed SCE in CCS and EKGE, whereas SCE surpassed PSO in EMO, MAES, and RMSES. Notably, PSO and SCE exhibited similar performance in NSES. This underscores the disparate spatial variability characteristics of optimal solutions influenced by distinct metrics (whereby the enrichment levels of optimal solutions at different sites reflect the extent of spatial variability). Additionally, notable outliers were identified in PSO's performance within CCS, EKGE, and EMO metrics, while both PSO and SCE exhibited outliers in the RMSES metric. This indicates that for RMSES, unquantifiable factors within the spatial variability of optimal solutions are more pronounced, whereas for other metrics, PSO's performance relative to SCE is more significantly influenced.

In summary, apart from PKGE and PMO, for other metrics, PSO typically exhibits better optimal solutions, i.e., enhanced effectiveness, compared to SCE, albeit at the cost of relatively lower efficiency. Notably, for CCS, EKGE, and RMSES, the optimal solutions obtained by PSO demonstrate higher kernel densities than those by SCE. Conversely, for EMO and MAES, the performance trend is reversed.

4.2.3. Optimal Simulation

Figure S2-1 presents linear fitting (s, r²) between simulations and observations of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ under varying metrics. For ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, PSO s (in descending order) are EMO, EKGE, RMSES, MAES, PMO, NSES, CCS, PKGE, with r² values also descending from EMO to PKGE. In contrast, SCE slope (s) are EMO, PMO, EKGE, MAES, PKGE, NSES, RMSES, CCS, with r² following a similar but slightly different descending order. For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$'s linear fitting (Figure S2-1-2), PSO s are EKGE, EMO, MAES, RMSES, CCS, NSES, PKGE, PMO, while r² values show a distinct ordering: PMO, followed closely by EMO/PKGE, then MAES/RMSES/NSES, EKGE, and finally CCS. SCE fitting for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ exhibits a different ordering for s (EKGE, EMO, CCS, RMSES, MAES, NSES, PMO, PKGE) and r² values (PKGE, PMO, NSES, EKGE, EMO, RMSES, with CCS and MAES closely grouped).

Generally, for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, except for NSES, PKGE, and PMO metrics, both PSO and SCE exhibit negative s values, while the rest are positive (Table 4). This indicates that most linear relationships between calibrated simulations and observations are positively correlated, which aligns with the improvement objectives of this study. Specifically, for EMO and EKGE, the s values of PSO (SCE) in the calibration of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ are 0.96 (0.83) and 0.18 (0.23), respectively, showcasing the optimal calibration performance (Figure 9). Furthermore, it is noteworthy that for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, the highest r² value of 0.11 is comparable to the lowest r² value observed in ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ (PKGE), implicitly suggesting a greater challenge in modeling ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$.

The Gaussian fitting, i.e., center (frequency) as c(f) with units of ${\mathrm{m}}^3·{\mathrm{m}}^{-3}$(1), of ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ in Figure S2-2-1 reveals: CTR’s ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ are widely distributed, peaking at ~0.15 (f=297). CCS, PSO, SCE errors span widely around -0.04, 0.11 (f≈350, 295). EKGE's PSO, SCE errors narrowly center at 0 (f≈1276, 608). EMO's PSO, SCE errors narrowly peak at 0, 0.01 (f≈1178, 700). MAES's PSO, SCE errors widen slightly at 0.01, 0.02 (f≈344, 416). NSES's PSO, SCE errors are wide at 0.05 (f≈274, 230). PKGE's PSO, SCE errors widely center at 0.08, 0.11 (f≈322, 325). PMO's PSO, SCE errors narrowly peak at 0.02, 0.03 (f≈480, 444). RMSES's PSO, SCE errors narrowly center at -0.02, 0 (f≈426, 296). Moreover, The Gaussian fitting, i.e., center (frequency) as c (f) with units of K (1), of OBS-SIM for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ (Figure S2-2-2) shows: CTR errors have a wide bimodal dist. centered at ~7.1, -3.8 (f≈192, 134). CCS, PSO, SCE errors widely center at ~2.3, 1.1 (f≈216, 167). EKGE's PSO, SCE errors widely center at ~1.3, 2.5 (f≈200, 203). EMO's PSO, SCE errors center at ~0.85, 1.23 (f≈170, 207). MAES's PSO, SCE errors center at ~-0.06, 0.88 (f≈200, 230). NSES's PSO, SCE errors widely center at ~5.86, 5.03 (f≈169, 213). PKGE's PSO, SCE errors widely center at ~4.91, 5.01 (f≈237, 152). PMO's PSO, SCE errors widely center at ~6.1, 5.19 (f≈300, 224). RMSES's PSO, SCE errors center at ~0.16, 1.29 (f≈200, 206).

In summary, for ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EKGE's performance in both PSO and SCE is closest to a normal distribution, whereas for that of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, MAES exhibits the closest resemblance to normality (Figure 10), with EKGE performing relatively poorly (Table 5). This underscores the significant influence of metric discrepancies on calibration simulation errors, contingent upon distinct calibration objectives. Furthermore, excessively wide peaks with low frequencies in unimodal distributions (e.g., CCS, NSES, PKGE, and PMO) indicate the dispersed fitting distribution, potentially necessitating the multimodal (e.g., more than two peaks) fitting. Conversely, bimodal distributions characterized by narrower peaks may call for a single-peak fitting centered around the modes.

Figure 11a depicts temporal ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ (${\mathrm{m}}^3·{\mathrm{m}}^{-3}$) variations for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. CTR's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ is generally largest, 0.15 (decreasing during July 5th and 10th rainfalls), with a slight upward trend. For CCS, PSO's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ 0.15 increases slightly, while SCE's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ fluctuates around 0.07. EKGE and EMO show PSO(SCE) RMSEs of 0.01(0.03) and 0.01(0.02), respectively, both trending downward. MAES's PSO/SCE ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ 0.04, both declining. NSES's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ around 0.07, up trending. PKGE's PSO(SCE) ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ 0.12(0.1), up trending. PMO's PSO ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ decreases from 0.1 to 0.05, SCE's 0.07, slightly down. RMSES's PSO(SCE) ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ 0.04(0.05), both declining. Moreover, Figure 11b illustrates the overall RMSES distribution for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. Median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ ranking from highest to lowest for PSO: CCS (0.13) > PKGE (0.12) > NSES (0.08) > PMO (0.07) > RMSES (0.039) > MAES (0.038) > EMO (0.018) > EKGE (0.017); for SCE: PKGE (0.1) > CCS (0.09) > NSES (0.085) > PMO (0.065) > RMSES (0.056) > MAES (0.039) > EKGE (0.03) > EMO (0.02). Notably, EKGE and EMO exhibit the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ for PSO and SCE, respectively, whereas CCS and PKGE have the highest.

Figure 11c depicts temporal variations in spatial correlation coefficients (${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$) for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ simulations. CTR's ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ significantly drops Jul 5-6 (0.5 to -0.6), fluctuating at ~0.2 otherwise. CCS: PSO's ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ stable at -0.5, SCE increases Jul 5 (-0.4 to 0.5), fluctuating -0.2. EKGE: PSO 1, SCE 0.8. EMO: Both are ~1. MAES: PSO increases (0.2 to 1), SCE initially declines (0.4 to 0), then ~0.2. NSES: PSO ~0.7, drops post-Jul 10 to ~0.5; SCE ~0.2. PKGE: PSO ~0.45, sharp drop Jul 5 to ~-0.4; SCE -0.4 to unspecified, sharp drop, ~-0.3. PMO: PSO 0.6, sharp drop Jul 6 to -0.5, rises to 0.2; SCE 0.8, drops to 0, rises to 0.4. RMSES: PSO increases (0.5 to 0.8), stabilizes ~0.8 post-Jul 4; SCE ~0.45, declines Jul 4-6, rises to ~0.26. EMO and EKGE consistently outperform MAES for PSO and SCE in ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$, with other metrics displaying varied trends. Moreover, Figure 11d illustrates the overall ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ distribution for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. Median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ ranking from highest to lowest for PSO: EMO > EKGE > MAES > RMSES > NSES > PMO > PKGE > CCS; for SCE: EKGE > EMO > MAES > PMO > PKGE > NSES > RMSES > CCS. Notably, EMO and EKGE exhibit the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for PSO and SCE, respectively, whereas CCS have the lowest.

For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, CTR's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ shows a marked diurnal variation, averaging 8K fluctuations (Figure S2-3-1). Due to overlapping diurnal error ranges, its performance complexity surpasses ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. Notably, NSES, PKGE, PMO peak ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ >14K (CTR's max), indicating inferiority (Figure 11e). Conversely, MAES and RMSES peak at 8K, surpassing CTR. EKGE and EMO, excluding initial days, also peak near 8K, outperforming CTR. Median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ (K) ranking from highest to lowest yields the following order for PSO: PMO (7.5) > PKGE (6) > NSES (5.8) > CCS (4) > EKGE (3.5) > MAES (2.8) > RMSES (2.5) > EMO (2.48); and for SCE: PKGE (6.1) > PMO (5.8) > NSES (5.6) > CCS (3.6) > EKGE (3.3) > MAES (2.9) > RMSES (2.7) > EMO (2.5) (Figure 11f). In both PSO and SCE, EMO exhibits the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$, whereas PMO and PKGE respectively possess the highest.

Furthermore, for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, CTR's ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ varies from -0.5 to 0.7, showing distinct diurnal patterns (Figure 11g). Overlapping diurnal error ranges complicate performance compared to ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ (Figure S2-3-2). CCS and EKGE's max ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ < 0.7 (CTR's max), indicating inferiority. NSE, PKGE, PMO max ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ rival CTR, but min ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ > -0.5, outperforming CTR. EMO, MAES, RMSE max ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ ~0.8, exceeding CTR. Hence, for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performance ranks EMO, MAES, RMSE best, followed by NSE, PKGE, PMO; CCS, EKGE perform bad. Moreover, Figure 11h illustrates the overall ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ distribution for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$. Median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ ranking from highest to lowest for PSO: EMO > MAES > RMSES > NSES > PMO > PKGE > EKGE > CCS; for SCE: EMO > RMSES > CCS > EKGE > MAES > PKGE > PMO > NSES. Notably, EMO exhibit the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for both PSO and SCE, whereas CCS and NSES have the lowest.

In summary, for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EKGE and EMO exhibit the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for PSO and SCE, respectively, whereas CCS and PKGE have the highest ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ for PSO and SCE, respectively, and CCS have the lowest ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for both. For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, in both PSO and SCE, EMO exhibits the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$, whereas PMO and PKGE respectively possess the highest; EMO exhibit the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for both PSO and SCE, whereas CCS and NSES have the lowest. Generally, EKGE and EMO have the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ for PSO and SCE respectively, while EMO has the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ for both.

4.3. Effects on Forecast

4.3.1. Linear and Gaussian Fitting

Figure S3-1 illustrates disparities in linear fitting (s, r2) between SIM and OBS for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ across metrics. For ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$'s linear fit (Figure S3-1-1), PSO s (descending): EKGE > EMO > MAES > RMSES > NSES > PMO > CCS > PKGE; r² order matches. For SCE, s: EMO > EKGE > PMO > MAES > NSES > RMSES > CCS > PKGE; r² differs: EKGE > EMO > MAES > PMO > NSES > PKGE > RMSES > CCS. For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$'s fit (Figure S3-1-2), PSO s: MAES > RMSES > EMO ≥ EKGE/CCS > NSES > PKGE > PMO; r²: PMO > PKGE > EMO > RMSES ≥ MAES/NSES > EKGE > CCS. SCE's s: MAES/CCS > RMSES > EMO/EKGE > NSES > PMO > PKGE; r²: PKGE > RMSES/PMO > MAES/EMO > NSES > CCS/EKGE.

Generally, in addition to NSES, PKGE, PMO in ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, and PKGE in ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, both PSO and SCE exhibit positive s values (Table 6). This indicates that most linear relationships between validation forecasts and observations are positively correlated, which aligns with the improvement objectives of this study. Specifically, for EKGE and MAES, the s values of PSO (SCE) in the validation of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ are 0.98 (0.84) and 0.14 (0.15), respectively, showcasing the best forecast performance (Figure 12). Furthermore, it is noteworthy that for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, the highest r² value of 0.1 is much smaller than the highest r² value observed in ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ (EKGE), implicitly suggesting a greater challenge in forecasting ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$.

The Gaussian fitting (c (f)) of ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ (Figure S3-2-1) reveal: CTR centered at 0.19 (f=272); CCS, PSO, SCE at 0.15, 0.07 (f=189, 225); EKGE, PSO, SCE narrowly at 0 (f=383, 363); EMO, PSO, SCE at 0 (f=416, 359); MAES, PSO, SCE at -0.01, 0 (f=359, 284); NSES, PSO, SCE at 0.06, 0.05 (f=343, 322); PKGE, PSO bimodal at 0.13, 0.04 (f=234, 220), SCE bimodal at 0.16, 0.06 (f=199, 365); PMO, PSO, SCE widely at 0.01, 0.04 (f=367, 323); RMSES, PSO, SCE at -0.02, 0.01 (f=293, 326). Similarly, for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ (Figure S3-2-2): CTR bimodal at 7.28, -3.57 (f=211, 160); CCS, PSO, SCE widely at 3.2, -0.38 (f=187, 181); EKGE, PSO, SCE at -0.09, 3.39 (f=143, 189); EMO, PSO, SCE at -1.41, -0.98 (f=175, 148); MAES, PSO, SCE at 0.49, 0.29 (f=181, 206); NSES, PSO, SCE widely at 5.81, 4.56 (f=204, 210); PKGE, PSO, SCE widely at 4.9, 5.7 (f=214, 217); PMO, PSO, SCE widely at 6.17, 5.47 (f=221, 187); RMSES, PSO, SCE at 0.55, 0.32 (f=194, 198).

Generally, for ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EMO’s and EKGE's performances in both PSO and SCE are closest to the normal distribution, whereas for that of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, EKGE in PSO and MAES in SCE exhibit the closest resemblance to normality (Figure 13), with EKGE performing relatively poorly (Table 7). This underscores the significant influence of metric discrepancies on forecast errors. Furthermore, excessively wide peaks with low frequencies in unimodal distributions (e.g., CCS, NSES, PKGE, and PMO) indicate the dispersed fitting distribution, potentially necessitating the multimodal (more than two peaks) fitting.

4.3.2. Spatial Difference and Similarity

Figure 14a depicts the temporal ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ (${\mathrm{m}}^3·{\mathrm{m}}^{-3}$) variations for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. CTR's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ is largest (~0.15), fluctuating with a dip on Jul 24. CCS's PSO RMSES ranges around 0.13, trending up, while SCE's remains stable at 0.1. EKGE’s and EMO's PSO/SCE ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ (~0.02 / 0.01) and (~0.02), both trends slightly up. MAES's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ are ~0.04. NSES's hover at 0.07, declining slightly. PKGE's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ is ~0.12. PMO's decline from ~0.07 to 0.05. RMSES's PSO/SCE ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ are ~0.04/0.06. Furthermore, Figure 14b illustrates the overall ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ distribution for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. Median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ ranking from highest to lowest for PSO, CCS(0.12) > PKGE(0.11) > NSES(0.07) > PMO(0.05) > MAES(0.04) > RMSES(0.036) > EMO(0.028) > EKGE(0.02); for SCE, PKGE(0.11) > CCS(0.1) > NSES(0.07) > RMSES(0.052) > PMO(0.05) > MAES(0.04) > EMO(0.02) > EKGE(0.019). Notably, EKGE has the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ for both PSO and SCE, whereas CCS and PKGE have the highest, respectively.

Figure 14c shows temporal ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ variations for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. CTR's ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ significantly drops from July 20^th to 21^st (0.1 to -0.7), stable at ~0.2 otherwise. PSO and SCE ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$in CCS hover around -0.3. EKGE's PSO ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ remains ~0.8, SCE jitters ~0.3. EMO's ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ ~1 for both methods. MAES's PSO ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ ~0.8, declining gradually; SCE's drops from 0.4 to 0 (July 5th), rises slightly to ~0.2. NSES's PSO ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ starts at 0.7, dropping to ~0.5 post-July 10th; SCE jitters ~0.2. PKGE's PSO CCS ~0.45, sharply drops to ~-0.4 post-July 5th; SCE initially jitters ~-0.4, spikes, then jitters ~-0.3. PMO's PSO CCS starts ~0.6, sharply drops to ~-0.5 (July 6th), rises to ~0.2; SCE similar, starts ~0.8, drops to 0, rises to ~0.4. RMSES's PSO ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ jitters, rises (0.5 to 0.8), stabilizing ~0.8 post-July 4th; SCE starts ~0.45, drops to ~0 (July 4th-6th), rises, fluctuating ~0.26. For ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$'s ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$, PSO and SCE consistently rank EMO > EKGE > MAES; others exhibit unstable/inferior performance. Figure 14d depicts the overall ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ distribution for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. Ranking metrics by the median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$, from highest to lowest for PSO: EKGE > EMO > MAES > RMSES > NSES > PMO > PKGE > CCS; for SCE: EKGE > EMO > MAES > PMO > NSES > RMSES > PKGE > CCS. Notably, EKGE has the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for both PSO and SCE, while CCS and PKGE have the lowest, respectively.

Figure 14e displays temporal ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ variations for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$. CTR's ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ exhibits pronounced diurnal fluctuations 10K. Metric performances are intricate due to overlapping diurnal error amplitudes (Figure S3-3-1). NSES, PKGE, PMO peak ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ > 15K (CTR's max), indicating inferior performance. EMO, MAES, RMSES max ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ < 7K, superior to CTR. CCS, EKGE extreme ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ 8K (except July 1-2), also outperform CTR. For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ hierarchically show EMO, MAES, RMSES best, followed by CCS, EKGE, with NSES, PKGE, PMO worst. Moreover, Figure 14f shows the ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ distribution for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, ranked by median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$. For PSO: PMO (7.5) > PKGE (6.9) > NSES (6.6) > CCS (4) > EKGE (3) > RMSES (2.9) > MAES (2.7) > EMO (2); for SCE: PKGE (7) > PMO (6.8) > NSES (6.4) > CCS (3.8) > EKGE (3.3) > MAES (2.7) > RMSES (2.5) > EMO (2.3). EMO has the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ for both methods, whereas PMO and PKGE have the highest for PSO and SCE, respectively.

Figure 14g presents temporal ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ variations for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, with CTR's ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ displaying strong diurnal fluctuations between -0.7 and 0.7. Overlapping diurnal error amplitudes complicate performance compared to ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ (Figure S3-3-2). Notably, EMO, MAES, RMSE metrics exceed CTR's extremes, demonstrating superior performance. Moreover, Figure 14h depicts the overall CCS distribution for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$. Ranking metrics by the median CCS, from highest to lowest for PSO: EMO > RMSES > MAES > NSES > PMO > PKGE > CCS > EKGE; for SCE: RMSES > CCS > EMO > MAES > EKGE > NSES > PKGE > PMO. Notably, EMO and RMSES has the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for PSO and SCE, respectively, while CCS and PMO have the lowest, respectively.

In summary, for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EKGE has the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for both PSO and SCE, whereas CCS and PKGE behave oppositely, respectively. For ST05cm, EMO has the lowest median ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ for both methods, whereas PMO and PKGE have the highest for PSO and SCE, respectively; EMO and RMSES has the highest median ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ for PSO and SCE, respectively, while CCS and PMO have the lowest, respectively. Generally, EKGE have the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ for both PSO and SCE, while EMO and RMSES has the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ for PSO and SCE, respectively.

4.3.3. Surface States Intercomparison

Figure 15 presents Taylor Diagram plots of calibrated and CTR simulations of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ during the forecast period, compared with observations and/or GLDAS data, across various metrics. For the comparison of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ simulations with observations (Figure 15a), CTR exhibits a root mean square difference ($\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}$) greater than 0.02 ${\mathrm{m}}^3·{\mathrm{m}}^{-3}$, surpassing other simulated metrics and GLDAS. However, the correlation coefficient ($\mathrm{C}\mathrm{C}$) between CTR and observations is above 0.5, outperforming other simulations and GLDAS except for EKGE and EMO metrics. Additionally, CTR's standard deviation ($\mathrm{S}\mathrm{T}\mathrm{D}$) reaches approximately 0.03 ${\mathrm{m}}^3·{\mathrm{m}}^{-3}$, significantly higher than that of other simulated metrics and GLDAS. Thus, EKGE and EMO metrics, when applied in PSO and/or SCE, effectively improve the simulation of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. In the comparison of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ with observations (Figure 15b), CTR, GLDAS, and multiple simulations demonstrate no skill. Nevertheless, like ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, simulations using EKGE and EMO metrics consistently yield the lowest $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}$ and $\mathrm{S}\mathrm{T}\mathrm{D}$, as well as the highest $\mathrm{C}\mathrm{C}$ among all evaluated metrics.

Furthermore, for the comparison of sensible heat flux (HFX) with GLDAS (Figure 15c), CTR displays higher $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}$ and lower $\mathrm{C}\mathrm{C}$ than other simulated metrics, albeit with a relatively low $\mathrm{S}\mathrm{T}\mathrm{D}$. This suggests that while most other metrics' HFX simulations outperform CTR in terms of $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}$ and $\mathrm{C}\mathrm{C}$, their STD values are relatively increased, with EKGE and EMO ranking top two in both PSO and SCE for $\mathrm{S}\mathrm{T}\mathrm{D}$. In contrast, for the comparison of latent heat flux (LH) with GLDAS (Figure 15d), CTR exhibits lower RMSD and higher $\mathrm{C}\mathrm{C}$ than other simulated metrics, but with a relatively high $\mathrm{S}\mathrm{T}\mathrm{D}$. Notably, CTR's LH simulation surpasses other metrics in both $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}$ and $\mathrm{C}\mathrm{C}$. Specifically, EKGE and EMO rank top two for both $\mathrm{S}\mathrm{T}\mathrm{D}$ and $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}$ in both PSO and SCE, which is a notable contrast to the findings for HFX.

In summary, compared with observations, the ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ simulations of EKGE and EMO exhibit higher Taylor diagram skill (TDS) in both PSO and SCE, significantly outperforming CTR. In contrast, when compared with GLDAS reanalysis, the TDS of HFX simulations for all metrics in both PSO and SCE are superior to CTR, whereas the performance of LH simulations is the opposite. Evidently, the enhancement of surface SM and ST simulations often yields more divergent surface flux simulation results, indicating the high complexity of modeling surface states and fluxes in arid regions.

4.4. Configure and Benefit

Figure 16 compares the parameter ranges of the "best metric’s simulations" between PSO and SCE, alongside the KGE values of various metrics for surface soil moisture simulations against observations. It is observable that in PSO, the optimal parameter range of EMO is larger than that of EKGE, whereas the opposite holds true for SCE, where EMO's optimal parameter range is smaller than EKGE's (Figure 16A). The KGE values of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ from optimal simulations of different metrics indicate that in PSO, EKGE achieves the highest KGE value, whereas in SCE, EMO attains the peak (Figure 16a). For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, however, EKGE's optimal simulation yields the highest KGE value in both PSO and SCE (Figure 16b). In terms of forecasted ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EKGE consistently produces the highest KGE values in both PSO and SCE. Conversely, for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, CCS achieves the highest KGE values in both PSO and SCE, with EKGE following closely.

Figure 17 illustrates the changes in ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reductions and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ increases of the simulations for various metrics and CTR during the calibration and validation periods. For ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, during the calibration period, most metrics, except CCS, exhibit a reduction in ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ compared to CTR, with EMO and EKGE showing the most significant improvements (Figure 17a), which is also reflected in their highest ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ (Figure 17e). During the validation period, EKGE and EMO stand out among the metrics, excluding CCS and PKGE, in terms of ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reduction (Figure 17b), again accompanied by the highest ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ (Figure 17e). For ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, during the calibration period, ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reductions relative to CTR are observed for most metrics except PKGE, PMO, and NSES, with MAES, RMSES, and EMO demonstrating the most pronounced improvements (Figure 17c), which also correspond to the highest ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ (Figure 17g). Similarly, during the validation period, ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reductions are observed for most metrics except PKGE, PMO, and NSES, with MAES, RMSES, and EMO continuing to show the most significant improvements (Figure 17d), accompanied by the highest ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ (Figure 17h).

In summary, the parameter uncertainty range of EMO is slightly smaller than that of EKGE in both PSO and SCE, but the two metrics exhibit a trade-off in terms of ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reduction and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ increase during the forecast and calibration periods for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$. Specifically, EKGE shows the greatest ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reduction for both calibration and forecast periods of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, and the largest ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ increase during the forecast period of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$. Notably, EMO demonstrates the largest ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reduction and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ increase for both forecast and calibration of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, while EKGE performs poorest. This is notably different from the clear advantage of EKGE observed in our previous studies [9], which can be attributed to the use of four layers in all objective metrics in this study. This suggests that the EKGE metric with different vertical dimensions (number of layers) can significantly impact the improvement capability of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ forecasts. Additionally, the failure of PKGE and PMO during the forecast and calibration of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ (e.g., even inferior to CTR) indicates the ineffectiveness of using surface-layer-dominated Pareto objectives and highlights the limitations of adjusting subsurface simulations through improvements in surface simulations within Noah LSM.

5. Discussion

Though a comparative analysis of eight kinds of the introduced multi-objective calibration has effectively portrayed the multifaceted impacts of metric differences on the joint SM—ST calibration under the month—long calibration—prediction framework, offering insights for ST modeling and prediction in semi-arid regions, this study is nevertheless subject to limitations: 1) the imperfect datasets such as the unavailability of the site scale forcing data [35,71], which may lay certain spatial effects on the simulations; 2) the absence of solutions for calibration schemes based on the dominant Pareto metrics (e.g., PKGE and PMO) (Figure 6), which might indicate a need for extended search time or fewer physical constraints.

Significant heterogeneities in optimal parameters are observed across different objective metrics (Figure S1-1 and Table 2), while within the same metric, the heterogeneities are relatively closer across different algorithms (e.g., PSO and SCE), indicating that the parameter heterogeneities are largely determined by the choice of metrics rather than the algorithms themselves. Furthermore, for all metrics, the uncertainty of optimal parameters associated with PSO is higher than that of SCE (Figure 5), consistent with our previous findings [9]. Additionally, the "Vegetation" type exhibits a general pattern of null values (e.g., NA) in PNL and non-positive values in ONR compared to other types (Table 3), suggesting higher relative uncertainty and more unexplainable factors within the "Vegetation" type parameters. In particular, the EKGE metric comprehensively performs best in reducing spatial heterogeneities and uncertainties in LSM parameters compared to other metrics.

Moreover, significant disparities exist in the fitness curves of identical algorithms with different metrics, whereas the differences are relatively minor in the fitness curves of distinct algorithms with the same metric (Figure 6), indicating that metric variations exert a more profound influence on the convergence efficiency of calibration than the algorithm itself. Notably, the rapid convergence observed in CCS and NSES, as well as the non-convergence in PKGE and PMO, likely signify the presence of locality and sub-optimality in the numerical solutions [50]. Furthermore, substantial variations in success rates across different metrics are evident in PSO, whereas minimal changes are observed in SCE (Figure 7), suggesting that the evolutionary capability of the PSO algorithm is constrained by metric differences. Analogous to convergence efficiency, the degree of enrichment in numerical solutions also exhibits greater disparities between identical algorithms with different metrics than between distinct algorithms with the same metric (Figure 8), highlighting that metric variations have a more profound impact on the validity of numerical solutions compared to the algorithm itself.

There are considerable discrepancies among different metrics in addressing the spatial complexity of land surface modeling (e.g., fitting, errors, and similarities of land surface states during calibration and forecasting), and they are sensitive to algorithmic variations (Section 4.3). Specifically, the EKGE metric and the EMO metric exhibit the best overall performance for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, respectively (Figure 9-14). However, the EKGE metric's performance for ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ is inferior to our previous study [9], which can be attributed to the consideration of subsurface soil in the SM-ST target dimension, also demonstrating significant vertical variability in regional land surface modeling targets [36,49]. Although the EMO metric is not the overall best for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, its relatively balanced performance with no notable weaknesses can mitigate this issue to some extent (Figure 16-17).

Overall, the selection of calibration objectives must be carefully considered due to the profound impact of metric differences on the spatial heterogeneity of parameters in calibration, the efficiency and effectiveness of calibration, and the spatial complexity of surface conditions. Especially, the establishment of an automated soil observation network at regional stations, with the entire SM—ST sates serving as the joint calibration objective, can enhance operational land surface applications. Considering the less restrictive nature of EMO that combine multiple metrics in various applications, the benefits are relatively more robust. Future research should strengthen the regional application of EMO to improve the representation of surface characteristics in regional weather and climate numerical forecasting.

6. Conclusions

The surface conditions are crucial for both regional hydrology and weather. Using ITPCAS dataset from April 1 to July 31, 2014, the present study investigates the performance of various multi-objective metrics that combined with the multi-parameter tables as varied criteria of GSA on enhancing the Noah LSM calibration and forecasting. Comprehensive comparisons are conducted among these enhancements such as the optimal land parameters, objectives, and simulations, and the objective-informed forecasts that brought by these different metrics, to identify the effect of metric’s diversity on SM-ST calibration and surface forecast. Results have shown that:

The case study presented herein can be succinctly characterized as a configuration-forecasting problem. Initially, in terms of model configuration, the forcing manifestations encompass locally elevated surface temperatures (>5°C), low relative humidity (<1%), feeble wind speeds (<5m s^-1), a shift in wind direction from south to north, low atmospheric pressure (586hPa), and minor hourly rainfall intensities (<5mm h^-1). Subsequently, within the default model parameter configuration, significant spatial disparities emerge in the parameter space due to the static parameters being set as the globally optimal defaults, while the initial parameters are derived from forecasts spanning the preceding three months. The surface forecasting challenge manifests in the form of consistent ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ simulations but poor ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ simulations. Considering the variations in default model parameters across the calibration and validation (or forecast) periods, these periods are analyzed separately. Specifically, during the calibration and forecasting periods, the slope (s) and goodness of fit (r²) for the ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ simulations under the default parameter configuration are 0.22/0.14 and 0.15/0.05, respectively, with Gaussian fits of their errors exhibiting positive skew distributions centered at 0.16 and 0.13 m³.m-³. In contrast, the s/r² values for the ST05cm simulations are -0.48/0.17 and -0.4/0.15, with their errors displaying broader bi-modal distributions.

Firstly, for the optimal parameters of SM-ST calibration of the same metric, SCE consistently achieves lower parameter uncertainty than PSO, albeit at the cost of relatively higher spatial heterogeneity; specifically, in terms of parameter uncertainty, MAES in PSO and CCS in SCE exhibit the smallest; as for parameter spatial heterogeneity, EKGE and EMO in PSO yield the smallest, while EKGE in SCE displays smallest. Moreover, apart from PKGE and PMO, for other metrics, PSO typically exhibits better optimal solutions, i.e., enhanced effectiveness, compared to SCE, albeit at the cost of relatively lower efficiency; notably, for CCS, EKGE, and RMSES, the optimal solutions obtained by PSO demonstrate higher kernel densities than those by SCE; conversely, for EMO and MAES, the performance trend is reversed. Furthermore, EMO’s and EKGE’s PSO (SCE) calibration of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ with the maximum upward slope as 0.96 (0.83) and 0.18 (0.23), respectively, showcase the optimal linear fitting (Figure 9); for ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EMO’s and EKGE's performances in both PSO and SCE are closest to the normal distribution, whereas for that of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, EKGE in PSO and MAES in SCE exhibit the closest resemblance to normality (Figure 13); EKGE and EMO have the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ for PSO and SCE respectively, while EMO has the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ for both.

EKGE’s and MAES’s PSO (SCE) ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ forecasts with the maximum upward slope as 0.98 (0.84) and 0.14 (0.15), respectively, showcase the best linear fitting (Figure 12). For ${\mathrm{E}}_{\mathrm{O}-\mathrm{S}}$ of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$, EMO’s and EKGE's performances in both PSO and SCE are closest to the normal distribution, whereas for that of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, EKGE (in PSO) and MAES (in SCE) exhibit the closest resemblance to normality (Figure 13). EKGE have the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ for both PSO and SCE, while EMO and RMSES has the best ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ performances of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ for PSO and SCE, respectively. Furthermore, compared with observations, the ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ simulations of EKGE and EMO exhibit higher Taylor diagram skill (TDS) in both PSO and SCE, significantly outperforming CTR. In contrast, when compared with GLDAS reanalysis, the TDS of HFX simulations for all metrics in both PSO and SCE are superior to CTR, whereas the performance of LH simulations is the opposite.

The parameter uncertainty range of EMO is slightly smaller than that of EKGE in both PSO and SCE, but the two metrics exhibit a trade-off in terms of ${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{S}}$ reduction and ${\mathrm{C}\mathrm{C}}_{\mathrm{S}}$ increase during the forecast and calibration periods for ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$. However, due to the failure of vertical dimension expansion, EKGE performs poor in ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ simulation improvement. Evently, EMO showcases the greatest benefit in surface forecast improvement among all metrics and with the hopeful low parameter uncertainties, which shows the most promising application performance.

Specifically, in the SM-ST calibration of Noah LSM, for optimal parameters, MAES (in PSO) and CCS (in SCE) exhibit the lowest levels of uncertainty; EKGE/EMO (in PSO) and EKGE (in SCE) yield the smallest spatial heterogeneity, while other metrics demonstrate nearly irregular or non-discriminatory patterns. For optimal solutions, apart from Pareto dominance-based metrics (e.g., PKGE and PMO), other metrics do not alter the generality of GSA algorithms (such as effectiveness and convergence domain). Notably, although CCS and NSES can accelerate GSA convergence, their impacts on model calibration and prediction remain highly uncertain or even negative. Furthermore, regarding optimal modeling performance for calibration and forecast compared to observations, substantial variations exist among different metrics, among them, EMO and EKGE yield the best ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ modeling abilities, while EKGE and MAES exhibit the best ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ modeling abilities, respectively. In terms of observation-simulation error fitting, EMO and EKGE in ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ both perform optimally, while EKGE (in PSO) and MAES (in SCE) in ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ demonstrate the best performance. For optimal spatial error and similarity performance in calibration and forecast periods, EKGE’s ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ performs best, while EMO’s ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ excels.

Overall, the metrics, apart from their impact on GSA itself (e.g., convergence), could significantly influence the performances (including parameters, numerical solutions, and simulations) in SM-ST calibration and prediction. Furthermore, the vertical dimensionality of the objective metrics in this study notably affects the modeling of ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$, indicating that the improvement of surface states through metrics based on subsurface soil conditions is not absolute. Additionally, since the optimal performance of different metrics in individual or joint modeling of ${\mathrm{S}\mathrm{M}}_{05\mathrm{c}\mathrm{m}}$ and ${\mathrm{S}\mathrm{T}}_{05\mathrm{c}\mathrm{m}}$ is not entirely consistent, the selection criteria for metrics in GSA applications are not unique. Specifically, EMO outperforms others in calibrating and predicting the surface layer states. These findings could enhance our understanding of the spatial complexity of parameters and simulations in surface forecasting in semi-arid regions, thereby facilitating the improvement of regional surface forecasting.