2.1 Snow-related processes in the Noah land surface model

In this study, we employ the Noah LSM (Chen et al., 1996; Koren et al., 1999; Ek et al., 2003) to simulate the single-site land surface processes (Mitchell, 2005), including the surface energy and water flux, and to verify energy and water budgets in the near-surface atmospheric layer by simulating the soil moisture, soil temperature, and snowpack. The Noah LSM is a stand-alone and one-dimensional column model, developed through multi-institutional cooperation. In the soil, to simulate soil moisture and soil temperature, we selected four layers with thicknesses of 10, 30, 60, and 100 cm, respectively, from top to bottom, for a total depth of 2 m. The model also evaluates various other variables, including skin temperature, snow depth, snow water equivalent, snow density, and canopy water content (Mitchell, 2005). The energy and water fluxes are calculated through the surface energy and water balance equations, respectively. Due to its adequate complexity and computational efficiency (Mitchell et al., 2004), the Noah LSM has been coupled to the operational NWP model of the Korea Meteorological Administration (KMA), named the Korean Integrated Model (KIM; Hong et al., 2018) – see Koo et al. (2017) for details of the coupled KIM–Noah LSM system.

The current Noah LSM (version 3.4.1) uses a single-layer representation of the snow processes considering a bulk snow–soil canopy layer (Sultana et al., 2014). If air temperature is less than 0 ^∘C, the resulting precipitation is considered snow. The fractional snow cover is determined as a function of snow water equivalent (SWE) using a generalized snow depletion curve. Snow albedo is calculated based on the fractional snow cover and the maximum snow albedo (Ek et al., 2003). Snow depth is represented by SWE and the bulk snow density (Jonas et al., 2009). The equations in the Noah LSM describe the heat exchanges at the snow–atmosphere and snow–soil interfaces as well as snow accumulation, sublimation, and melting (Suzuki and Zupanski, 2018). The abovementioned snow processes contain certain estimated coefficients or constants, known as parameters, which employ typical, empirical, or a priori values. The parameters are provided as lookup tables based on their samples in the field or lab. Traditionally, they are tuned by trial and error to calibrate the model against historical observations in a specific location; however, a systematic and objective procedure is essentially required for a large number of stations (Duan et al., 2006; Rosolem et al., 2013). We explain below the details of the snow-related parameters to be optimized for various stations in SK.

2.1.1 Fractional snow cover (FSC)

The FSC (σ_s) is important for the accumulation and ablation processes (Livneh et al., 2010). As a function of SWE (W_s; in meters) extracted by the atmospheric input values (Livneh et al., 2010), σ_s varies nonlinearly as in Eq. (1), following the empirical snow depletion curves of Anderson (1973):

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{\mathrm{s}}=\mathrm{1}-e^{-P_{\mathrm{s}}W}+W{e}^{-P_{\mathrm{s}}}.\end{array}$$

Here, P_s is the distribution shape parameter, and $W=W_{\mathrm{s}}/W_{\max }$, where W_max is the threshold of W_s above which σ_s is 100 %. Note that, from Eq. (1), σ_s is a function of P_s and W_max – these two parameters are to be optimized.

Figure 1 represents the responses of the snow variables to the variations in the snow-related parameters for given ranges. It is noteworthy that P_s has a positive correlation with snow cover (Fig. 1a). For example, σ_s increases as P_s increases, resulting in relatively slow snowmelt. In Eq. (1), the value of P_s usually ranges between 2 and 4 (e.g., Anderson, 1973; Koren et al., 1999), and its default value in the Noah LSM is 2.6. We seek the optimal value of P_s, which lies between 2 and 4 and is suited to SK.

Figure 1Responses of the snow variables to the variations in the snow-related parameters for given ranges. (a, b) Responses of FSC, for W_s=0.02 m, to variations in P_s (with W_max=0.08 m) and in W_max (with P_s=2.6), respectively. (c, d) Responses of SA, for ${\mathit{\alpha }}_{\max ,\mathrm{sat}}=\mathrm{0.2}$ and t=10 d, to variations in α_max,CofE (with C=0.5) and in C (with ${\mathit{\alpha }}_{\max ,\mathrm{CofE}}=\mathrm{0.85}$), respectively. (e, f) Responses of SD (cm), for W_s=0.02 m and $T_{\mathrm{air}}=-\mathrm{5}$ ^∘C, to variations in P₁ (with P₂=0.0017 g cm^−3 ∘C⁻¹) and in P₂ (with P₁=0.05 g cm⁻³).

Download

The SWE threshold, W_max, has a negative correlation with snow cover, as shown in Eq. (1), and it is more sensitive compared to P_s within a given parameter's range (Fig. 1b). In the Noah LSM, the values of W_max are prespecified in a table (VEGPARM.TBL), varying with the land cover types (LCTs). W_max has the largest value over forest, reflecting the irregular geometry of forest cover (Livneh et al., 2010). Previous studies suggest the uncertainty range in the values of W_max; for instance, Livneh et al. (2010) used 0.04 m for forest and 0.02 m for non-forest, respectively, whereas Wang and Zeng (2010) used 0.2 m for tall vegetation and 0.01 m for short vegetation. The default values in the Noah LSM are 0.08 m for forest and 0.04 m for non-forest. We estimate the optimal W_max values, suited to SK, in the range between 0.01 and 2 m.

2.2 Optimization tool: micro-genetic algorithm

The genetic algorithm (GA) is a global optimization algorithm developed by John Holland in the 1970s (e.g., Holland, 1973, 1975) and is based on Darwinian principles of natural selection (Golberg, 1989). It uses reproduction selection, crossover, and mutation to operate a set of potential solutions, i.e., population or individuals, which are expressed by a string, called a chromosome: its binary form is called a gene (Krishnakumar, 1990; Rudnaya and Santosa, 2000). The selection operator first selects good solutions or eliminates bad solutions based on the fitness value; then, the crossover operator exchanges the genetic information between the solutions using the single-point or uniform types. The mutation operator modifies the value of each gene of the chromosomes by replacing it with the opposite value, e.g., 0 with 1, which prevents premature convergence. When a new generation is created, the above processes are repeated until the convergence condition or the prescribed number of iterations is satisfied.

The micro-GA is an advanced and simplified GA with smaller generation sizes, thus requiring less computational time than the conventional GA (Krishnakumar, 1990; Wang et al., 2010). It has been used in meteorology for optimal parameter estimation (e.g., Yu et al., 2013) or scheme-based optimization (e.g., Hong et al., 2014, 2015; Park and Park, 2021; Yoon et al., 2021). Its main difference from the conventional GA is the population size; for example, the micro-GA uses 5 individuals, while the conventional GA uses more than 30 individuals. Note that the conventional GA with a small population quickly converges to non-optimal solutions due to insufficient information; however, the micro-GA solves this problem by using elitism, which assigns the best individual among the five individuals based on the fitness evaluation and carries it to the next generation – this guarantees the preservation of the good solutions during the generations. Furthermore, the micro-GA does not take mutation to achieve diversity; instead, it uses the re-initialization which starts with a new individual whenever the diversity is lost.

2.2.1 Coupling the micro-GA with the Noah LSM and parallelization

Figure 2 describes the process of parameter optimization in the micro-GA–Noah LSM coupled system. (1) The micro-GA initializes the snow parameter combinations represented by the binary encoding through the random samples of the individual. (2) The micro-GA controls the Noah LSM by editing the parameter-related files, such as GENPARM.TBL, VEGPARM.TBL, and the Fortran code (module_sf_noahlsm.F), and prepares the forcing data for each station. (3) As recommended in Carroll (1996), the five individuals configured with the different snow parameters execute the ensemble runs of the Noah LSM in parallel. (4) The performance of each Noah LSM is evaluated in comparison with the observation through a given fitness function. (5) The micro-GA selects the highest fitness by comparing a number of individuals through the tournament selection. (6) New combinations for the next generation are produced through the crossover using the selected ones in the previous step. (7) When the convergence is satisfied, the other four individuals except the best individual marked by elitism are randomly regenerated. (8) The micro-GA repeats these processes until the prescribed entire iteration converges into a global maximum of the fitness function.

Although the micro-GA is computationally more efficient than the conventional GA, it still demands substantial computing time because each individual serially executes the model. Therefore, we have developed a parallel processing system in the micro-GA–Noah LSM coupled system. Instead of sequentially performing each individual and calculating the fitness within a generation, we run the model simultaneously for all populations to obtain the fitness and select the best individual when all tasks are finished (see the dashed box in Fig. 2). This new parallel system linearly reduces the execution time, which is proportional to the number of individuals. In addition, since the coupling system was created in a shell script, it is possible to assign multiple cores for model execution for various stations. The new parallel processing system, created by reflecting these two main points, improves the computation time – making it different from the non-parallel processing of a coupled system, e.g., the coupled micro-GA and the Noah LSM with multiple physics options (Noah-MP) model (see Hong et al., 2014).

Figure 2A flowchart of parameter optimization from the micro-GA–Noah LSM coupled system. The dashed box depicts the parallel system for the Noah LSM, running for each individual.

Download

2.3 Data

The land surface processes were forced by six meteorological fields from ASOS (https://data.kma.go.kr, last access: 24 October 2022): wind speed (m s⁻¹), wind direction (degrees), temperature (K), relative humidity (%), surface pressure (hPa), and precipitation rate (kg m⁻² s⁻¹). When missing data exist in less than 72 h, linear interpolation was performed except for precipitation. Stations with a missing rate greater than 1 %, during the entire experimental period, have been excluded. For the initial and boundary conditions, downward shortwave/longwave radiation (W m⁻²), precipitation rate (kg m⁻² s⁻¹), soil temperature (K), soil moisture (m³ m⁻³), and surface temperature (K) have been obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) – the fifth-generation ECMWF reanalysis-Land (ERA5L) hourly data (Muñoz-Sabater, 2019) – having a spatial resolution of 9 km and four soil layers with thicknesses of 7, 21, 72, and 189 cm, respectively, from top to bottom for a total depth of 2.89 m. We have used the data at the ERA5L grid's nearest point to the ASOS station.

The snow observations (i.e., SD, FSC, and SA) are used for the model verification and the fitness function calculation. For SD, the hourly model outputs are evaluated using the hourly ASOS data. To confirm the snow season, we have excluded the SD observations lower than 0.1 cm. For FSC and SA, we have no ASOS observations over SK; thus, we have used the MODIS/Terra Snow Cover Daily L3 Global 500 m SIN Grid radiance data (Hall and Riggs, 2021). They are generated from the MODIS/Terra Snow Cover 5-Min L2 Swath 500 m data (Hall et al., 2006) by selecting the best observation based on a scoring algorithm when they are closest to nadir with maximum coverage of the cell (Hall and Riggs, 2007). In particular, FSC is generated by the normalized difference snow index (NDSI). The MODIS snow data at the points nearest to the ASOS locations were extracted and used for verification of the model-generated FSC and SA. Being a polar-orbiting satellite, MODIS contains only one observation per day; thus, we have extracted the model output for verification at 02:00 UTC, when the satellite (Terra) passes over SK. For the calculations, we have converted the percent values of FSC and SA to the decimal values; then, we excluded observational data with values below 0.05 (i.e., 5 %) for both FSC and SA.

For the optimization experiment, we have selected some stations that represent different land covers in SK, aiming at having a representative combination of snow-related parameters over SK. We have defined a representative set of LCTs within a 2.5 km radius from the ASOS observations, excluding the water body. The LCTs have been taken from the MODIS (on board Terra and Aqua) Land Cover Type Yearly Climate Modeling Grid (CMG) Version 6 (Friedl and Sulla-Menashe, 2015), in which maps are provided from the land cover classification schemes of the International Geosphere-Biosphere Programme (IGBP), the University of Maryland (UMD), and the leaf area index (LAI), all at a 0.05^∘ spatial resolution in geographic latitude/longitude projection (see Sulla-Menashe and Friedl, 2018), for the entire globe from 2001 to 2019. Finally, we have compiled a set of five representative stations for each different LCT – deciduous broadleaf forest (DBF), mixed forest (MF), woody savanna (WS), cropland (CL), and urban and built-up lands (UB) – as shown in Table 1.

Table 1Five representative LCTs over SK, following the IGBP classification – DBF, MF, WS, CL, and UB. For each LCT, five selected stations are shown with the station name (abbreviation in parentheses), location in latitude (^∘ N) and longitude (^∘ E), ratio of LCT in a 2.5 km buffer (%), soil type, and missing ratio (%). The OPT_5 experiment employs only the stations highlighted in bold, while the other experiments use all the stations.

Download Print Version | Download XLSX

4.1 Spinup analysis

Numerical prediction models generally require spinup to reach a statistical equilibrium state where the initial conditions under a forcing are adjusted to the model's own physics/dynamics and numerics (Bonekamp et al., 2018). Without sufficient spinup, the LSMs can generate severe bias of initial conditions (Cosgrove et al., 2003). Prior to the optimization experiments, we have conducted a spinup experiment in one of the stations, Seoul, to check the appropriate spinup time. It was carried out in two ways: (1) using a spinup period recursive in 9 years (e.g., Jun et al., 2020) and (2) using a spinup period that was not included in the analysis.

First, the Noah LSM has been repeatedly executed using the atmospheric forcing for 9 years. This recursive simulation has been conducted from 1 May 2009 to 30 April 2018 to see whether the model was able to reach an equilibrium by setting the repetition loop to 0, 300, 600, and 1000. Our results indicated no significant differences; thus, we concluded that repetition was not required. Second, we have performed sensitivity tests to identify the spinup period due to changes in the initial conditions by adding biases (±0.1 m³ m⁻³ for soil moisture and ±3 K for soil temperature) to the ERA5L data. As a result, we found that the adequate spinup periods were about 3 months and 1 year for soil moisture and soil temperature, respectively; however, the snow variables were insensitive to the initial condition changes, thus requiring no spinup period. Although the spinup is not necessary for this study that focuses on the snow processes, we have performed the optimization experiments starting from May when snow is absent.

4.2 Optimal estimation of snow parameters

To optimize snow parameters specialized in SK, we have employed the micro-GA–Noah LSM coupled system using the observations over SK. Figure 4a shows the evolution of the fitness function for OPT_5 in a total of 200 generations and Fig. 4b for OPT_W in a total of 100 generations. Since OPT_W optimizes solely the W_max parameter, it has smaller generations. In OPT_5, the fitness function converges at the 160th generation, while the fitness function of OPT_W quickly converges in all LCTs (Fig. 4b). The convergence occurs at the 3rd generation for DBF, 70th generation for MF, 7th generation for both WS and CL, and 12th generation for UB.

Figure 4The fitness function for generations during (a) five-snow-parameter optimization (OPT_5) and (b) W_max optimization (OPT_W) for DBF (black), MF (blue), WS (green), CL (yellow), and UB (red) LCTs.

Download

As a result, we have obtained the optimized six snow parameters over SK (Table 3). OPT_5 simultaneously generates the optimized five snow parameters (P_s, α_max,CofE, C, P₁, and P₂) associated with the FSC, SA, and SD, while OPT_W, depending on the LCTs, generates the optimized W_max associated with the FSC. The first snow parameter, P_s, is optimized from its standard value of 2.6 to 2.7097, which results in an increase in the FSC. The second snow parameter, W_max, is optimized depending on each LCT. In detail, the W_max in DBF and WS increases from 0.08 to 0.1632 m and from 0.03 to 0.0406 m, respectively. They lead to a decrease in the FSC due to a negative correlation. On the other hand, the W_max in MF and UB decreases from 0.08 to 0.0529 m and from 0.04 to 0.0284 m, respectively, thus increasing the FSC. The optimized CL shows a similar value from 0.04 to 0.0406 m, which means that the current value was correct for SK. The third snow parameter related to the SA, α_max,CofE, decreases from 0.85 to 0.7387, inducing a decrease in SA. The fourth snow parameter, C, also shows a similar value from 0.5 to 0.5355, and thus this value was correct for SK. The fifth snow parameter, P₁, increases from 0.05 to 0.0698 g cm⁻³, resulting in a decrease in SD. The last snow parameter, P₂, decreases from 0.0017 to 0.0002 g cm^−3 ∘C⁻¹, leading to an increase in SD.

Table 3Summary of optimized snow parameters related to snow variables. Minimum (min), default, and maximum (max) are the ranges used in the optimization process. Default is the empirical value used in the Noah LSM.

Download Print Version | Download XLSX

We have investigated the mean bias (MB) using the box plot expressing the quartile and the distribution of extreme values: it explains how much the bias of the CNTL is reduced in optimization experiments by comparing the model with the observations. Before optimization, the CNTL showed underestimated FSC and SD and overestimated SA (−0.133, −4.39 cm, and 0.0408, respectively; see Fig. 5). However, the bias patterns in FSC and SA vary for each station owing to the lower spatial and temporal resolution of satellite observations. On the other hand, the SD shows an underestimation at all stations; the increase in the SD due to fresh snow was underestimated, and snowmelt was proceeding faster than the observation.

Figure 5Box plots of (a) FSC bias, (b) SA bias, and (c) SD bias (cm) for CNTL, VRF_5, and VRF_6. The maximum differences are indicated with the black star symbol (e.g., 0.637 (CNTL), 0.643 (VRF_5), 0.570 (VRF_6) for FSC, 0.605 (CNTL), 0.563 (VRF_5), and 0.525 (VRF_6) for SA and 34.1 cm (CNTL), 45.1 cm (VRF_5), and 46.3 cm (VRF_6) for SD). Each mean of snow variables is indicated as a black circle (e.g., −0.133 (CNTL), −0.145 (OPT_5), and −0.149 (VRF_6) for FSC, 0.0408 (CNTL), 0.0298 (VRF_5), and 0.0281 (VRF_6) for SA and −4.39 cm (CNTL), −2.81 cm (VRF_5), and −2.45 cm (VRF_6) for SD).

Download

Figure 6Scatter plots of observations (OBS) and model results (LSM) for snow variables FSC (left panels), SA (middle panels), and SD (cm; right panels) from the verification experiments – CNTL (black dots), VRF_5 (blue dots), and VRF_6 (orange dots), which are evaluated over different LCTs. (a–c) DBF represented by station UL, (d–f) MF by GM, (g–i) WS by NG, (j–l) CL by BR, and (m–o) UB by SL.

Download

The performance has been evaluated using the improvement ratio, which indicates how much the RMSE, MB, and coefficient of determination (R²) of experiments using optimized parameters (i.e., VRF_5 and VRF_6) are improved compared to CNTL, as shown in Eq. (7) (Table 4). In VRF_5, new parameter values – P_s, α_max,CofE, C, P₁, and P₂ – optimized by the micro-GA result in an improvement in the RMSEs for FSC, SA, and SD of 0.7 %, 5.4 %, and 13.7 %, respectively. However, the RMSE of FSC relatively weakly improved by about 0.7 % because the other parameter, W_max, is not yet optimized. In terms of MB, we anticipate that the increase in P_s will overcome the underestimated FSC. However, VRF_5 strengthens the underestimation of FSC from −0.133 to −0.145, and thus it deteriorates the MB by about 9.1 % (Table 4 and Fig. 5a). Regarding the SA, the optimized α_max,CofE decreases the SA to solve the overestimation in CNTL. The other parameter, C, has been optimized to its default value, 0.5355, which means that this was an appropriate constant for SK snowfall prediction. Therefore, the MB of the SA is improved by 26.9 % by reducing the SA from 0.0408 to 0.0298 (Table 4 and Fig. 5b). Next, SD shows the greatest RMSE improvement of 13.7 %. In fact, the Noah LSM suffers from a negative bias for SWE, especially in early spring (Sheffield et al., 2003; Ek et al., 2003; Pan et al., 2003; Mitchell et al., 2004; Jin and Miller, 2007; Livneh et al., 2010). Because SD is proportional to SWE, the underestimation can be exhibited due to negative bias of the SWE. However, the optimized P₁ leads to a decrease in SD, and thus it intensifies the underestimation for SD. On the other hand, the optimized P₂ increases the SD as follows: when the air temperature is warmer than −15 ^∘C, the fresh snow density slowly increases, which quickly induces an increase in SD following Eq. (5). Therefore, the optimization of P₂ solves the underestimated SD by about 35.9 % due to the increased SD from −4.39 to −2.81 cm within most of the temperature ranges (Table 4 and Fig. 5c). We also investigated R², which measures the proportion of variation for a dependent variable that can be explained by an independent variable. Although the R² values are low in FSC and SA, the difference between CNTL and the verification experiment (e.g., VRF_5) has 95 % statistical significance, as evaluated with a two-tailed t test. After optimization, the R² values in VRF_5 improve by 3.3 % and 1.5 % for FSC and SD, respectively. However, these changes are insignificant compared to the other statistics such as RMSE and MB.

Table 4The RMSE, MB, and R² of snow variables and improvement ratios (%) in parentheses from CNTL to VRF_5, and VRF_6 over the 25 representative stations. The difference between CNTL and verification experiments (i.e., VRF_5 and VRF_6) has 95 % statistical significance, as evaluated with a two-tailed t test.

Download Print Version | Download XLSX

To supplement insufficient improvement in the FSC, we have additionally optimized W_max as a function of LCT (OPT_W) using the optimized values of five parameters from OPT_5. Here, we have only used the FSC to define the fitness function, not considering SA and SD; thus, the fitness function is defined using Eq. (8), where the FSC is the only element of x, and the normalized process with Eq. (9) is not necessary. As a result, OPT_W further improves the RMSE of the FSC in VRF_6 compared to VRF_5 in most stations: the significant decreases in W_max over MF and UB lead to an increase in the FSC, possibly alleviating the underestimation problem of the FSC in VRF_5.

Figure 7Time series of the snow variables for DBF (e.g., UL) from May 2009 to April 2018: (a) FSC, (b) SA, and (c) SD (cm). Observations are in gray dots, and model results are in black dots for CNTL and in orange dots for VRF_6.

Download

Finally, all six parameters related to the snow variables have been verified in VRF_6 with the same 25 stations used in CNTL. When the optimized five parameters are used except W_max (VRF_5), SA and SD are improved, and FSC shows a weak improvement in RMSE performance (Table 4). However, when the optimized W_max depending on the LCTs from OPT_W is used (VRF_6), the FSC appears in a larger positive impact with other variables. As a result, an improvement in RMSE for FSC, SA, and SD is 3.3 %, 6.2 %, and 17.0 %, respectively. However, the MB for the FSC strengthens from 9.1 % to 11.9 % in VRF_6 (Table 4 and Fig. 5a) due to larger negative bias, especially in the DBF. On the other hand, SA and SD reduce the MB against the CNTL and enhance the improvement ratio from 26.9 % to 31.0 % and from 35.9 % to 44.2 %, respectively (Table 4 and Fig. 5b–c). Like the RMSE, the R² of FSC and SD also improved in VRF_5 and VRF_6. The SA worsened in VRF_5 and was a bit more severe in VRF_6. However, they are still small impacts compared to RMSE and MB.

To understand more details of the improvements due to the optimization, we analyzed the scatter plots that compare the observations and the model results in Fig. 6 and listed their RMSE and R² in Table 5. Since the observation patterns are different for different stations, we selected the representative station for each LCT. For FSC, it is relatively hard to recognize the explicit bias patterns, as shown in Fig. 6 (left panels); however, compared to CNTL, the RMSE decreased in VRF_5 and further reduced in VRF_6 (see Table 5). The VRF_6 revealed the largest R² values over most LCTs, except WS (station NG) and CL (station BR). In particular, VRF_6 produced the highest FSC over MF (station GM) (see Fig. 6d), with the smallest RMSE and the largest R², which significantly alleviated the underestimation problem. For SA, its overestimation in CNTL has been prominently reduced in both VRF_5 and VRF_6 – see Fig. 6 (middle panels). For instance, SA decreased over DBF (station UL) in both VRF_5 and VRF_6, with a larger decrease for VRF_6 (Fig. 6b). The performance statistics of both VRF_5 and VRF_6 demonstrated improvements over most LCTs except UB (station SL) (see Table 5). For SD, the parameter optimization brought about remarkable improvement compared to FSC and SA – see Fig. 6 (right panels). Note that SD is optimized using the hourly in situ observations (i.e., larger amount of data), while both FSC and SA are optimized using the daily satellite observations. For example, VRF_6 with DBF produced notably large SD values (Fig. 6c) with the lowest RMSE and the highest R² (Table 5), diminishing the underestimation problem in CNTL. It is hard to say which verification experiment gives the best results (i.e., VRF_5 versus VRF_6), but the performance with optimized parameters is usually better than CNTL in terms of RMSE (e.g., for most LCTs such as DBF, MF, WS, and UB) and R² (e.g., for LCTs including DBF, MF, and CL). Overall, both VRF_5 and VRF_6 produced snow variables that are closer to observations than CNTL for most LCTs (i.e., stations), and VRF_6 generally showed the lowest RMSE and the highest R² in all the snow variables.

Table 5Statistics of model performance using non-optimized parameters (CNTL) and optimized parameters (VRF_5 and VRF_6) over different LCTs represented by different stations – DBF represented by UL, MF by GM, WS by NG, CL by BR, and UB by SL. The RMSEs and R² values are shown for three snow variables – FSC, SA, and SD.

Download Print Version | Download XLSX

Figure 7 shows temporal changes in the snow variables after parameter optimization by comparing their time series of the observations and the model simulations (CNTL versus VRF_6) for DBF represented by UL. The CNTL shows positive or negative biases in FSC, positive bias (overestimation) in SA, and negative bias (underestimation) in SD: these biases are all reduced down in VRF_6. The bias patterns in Fig. 7 are consistent with those in Fig. 6a–c.

Lastly, we have investigated how the optimized snow parameters can affect the other variables in the LSM. Figure 8 depicts the time series of the differences of LSM variables (soil temperature, sensible heat flux, and soil moisture) between VRF_6 and CNTL (i.e., VRF_6 minus CNTL) following the changes in SD. Although the LSM variables here are not directly optimized, they respond to the optimized snow parameters through associated physical processes. Note that the underestimation of SD in CNTL has been alleviated in VRF_6 by using the optimized snow parameters (see Figs. 7c and 8a). Next, soil temperature in the first soil layer (7 cm) increases as SD increases after optimization, which consequently increases sensible heat flux. The residual of the surface energy balance is close to zero, implying that the surface energy balance is well conserved even after optimization. Soil moisture depends on snowmelt, following the trend of increased snowfall in the previous winter. Extreme fluctuations sometimes appear in the time series analyses due to nonlinear effects, but we can understand the overall tendency according to the increased SD on the land surface.

Figure 8Time series of the difference between CNTL and VRF_6 for the UL in DBF from May 2009 to April 2018: (a) SD (cm), (b) soil temperature in the top soil layer (7 cm) (ST; K), (c) sensible heat flux (SH; W m⁻²), and (d) soil moisture in the top soil layer (7 cm) (SM; m³ m⁻³).

Download