Improvements in the Land and Crop Modeling Over Flooded Rice Fields by Incorporating the Shallow Paddy Water

Abstract

Flooded rice paddies are important for modifying land surface energy and water budgets, especially in Asian countries. This study incorporated shallow paddy water into the Noah with Multi-Parameterization (Noah-MP) model to enhance its performance in capturing the distinct features of small Bowen ratios over flooded rice fields. The paddy surface water was parameterized as one integrated layer along with the top soil layer, and meteorological measurements from two crop sites in Japan, that is, SAITO (early rice) and SAGA (late rice), were employed for model evaluation at the field scale. The simulation results show that the model performance was significantly improved by combining the incorporation of paddy water and the calibration of rice crop parameters, particularly at SAGA. Compared with the reference run using the original version of Noah-MP for SAGA, the underestimation in latent heat and the overestimation in sensible heat during daytime were decreased by ∼74 W m⁻² (∼67%) and ∼92 W m⁻² (∼55%), respectively. Approximately 60%–70% of this improvement was contributed by using calibrated rice crop parameters, while the rest of 30%–40% was from further incorporating paddy water. The decreased ground surface resistance owing to the presence of paddy water was crucial for capturing the features of small Bowen ratios. The observed water depth might help mitigate the underestimation of latent heat nonlinearly. This work may benefit the study of land-atmosphere interactions and local and regional weather and climate in Asia with the widely used coupled Weather Research and Forecasting/Noah-MP model.

Introduction

Agriculture is by far the largest water-consuming sector and agricultural irrigation constitutes about 70% of freshwater withdrawals worldwide (Food and Agriculture Organization, 2010). Irrigated agriculture plays a crucial role in global food production and security, accounting for ∼20% of total cropland but producing approximately 40% of the world's food (Siebert and Döll, 2010). The past century has witnessed a rapid expansion of irrigation, with the global area equipped for irrigation increasing from ∼53 million hectares (Mha) in 1901 (Wisser et al., 2010) to ∼310 Mha circa 2005 (Siebert et al., 2013). With the growing population and increasing demand for food, global irrigated areas could further increase to as high as 1,800 Mha by 2050 (Puy et al., 2020). To properly address the challenges in food security and water scarcity in the context of global warming, it is essential to advance the understanding of the nexus among agriculture, water, and climate.

Extensive studies have demonstrated the significance of irrigation in altering land surface water and energy budget and in modifying weather and climate at local and regional scales. By transferring water from surface or subsurface water bodies (e.g., rivers, lakes, reservoirs and aquifers) to the cropland surface, irrigation tends to cool and moisten the air over irrigated lands via the repartitioning of net radiation into sensible and latent heat fluxes (Bonfils & Lobell, 2007; Chen et al., 2018; Kueppers et al., 2007; Leng et al., 2013; Ozdogan et al., 2010; Sacks et al., 2009; Xu et al., 2019). In some cases, irrigation-induced abundance of water vapor in the atmosphere can enhance cloud cover, convection, and downstream precipitation (Pei et al., 2016; Qian et al., 2013; Sacks et al., 2009). Additionally, the temperature and/or moisture contrasts between irrigated croplands and nearby non-irrigated areas may cause changes in the mesoscale or regional circulation patterns (Douglas et al., 2009; Guimberteau et al., 2012; Saeed et al., 2009; Shukla et al., 2014; Tuinenburg et al., 2014). Due to the limited observations of water use for irrigation, most previous studies on the regional and global estimates of irrigation effects were based on numerical simulations, and their main concern was the representation of irrigation-induced changes in soil moisture conditions rather than crop species (Boucher et al., 2004; Kueppers et al., 2007; Leng et al., 2013; Ozdogan et al., 2010; Sacks et al., 2009). Some recent modeling studies have focused on the high dependence of irrigation implementation on crop species (McDermid et al., 2017; Xu et al., 2019). For example, Xu et al. (2019) indicated that it was imperative to use crop-specific triggering thresholds to model irrigation in maize and soybean fields. However, few studies have distinguished paddy rice from other crops when investigating the climatic responses to irrigation.

As a staple food for ∼50% of the global population, rice is typically grown in flooded paddy fields and requires a large amount of water during its course of growth. Irrigated rice is estimated to consume approximately 34%–43% of the total worldwide freshwater withdrawals for irrigation (Bouman et al., 2007). Generally, the growing-season water input to rice fields is about 1,300–1,500 mm (Bouman et al., 2007) in Asia, where over 80% of global rice production and ∼75% of global consumption occur (Sekhar, 2018). Using climatic datasets as model inputs, some offline models, such as CROPWAT (Smith, 1992) and ORYZA2000 (Bouman & van Laar, 2006; Bouman et al., 2001) were widely used to estimate irrigation water requirements for rice fields as well as paddy rice crop performance in a changing climate (e.g., Chung & Nkomozepi, 2012; Devkota et al., 2013; Kuo et al., 2006; Wang et al., 2014). How do weather and climate respond to heavily irrigated rice paddy fields? To address this issue, a land surface model (LSM) coupled with a crop growth model for paddy rice would be very useful. Masutomi, Ono, Mano, et al. (2016) and Masutomi, Ono, Takimoto, et al. (2016) developed the MATCRO-Rice model based on the MATSIRO LSM (Takata et al., 2003), and a new term for heat flux stored in the paddy surface water was incorporated by following the model of Maruyama et al. (2017) and Maruyama and Kuwagata (2010) to consider the case of flooded rice paddies. In addition to representing paddy water effects on surface energy cycle as in Masutomi, Ono, Mano, et al. (2016) and Devanand et al. (2019) also introduced paddy field hydrology (Xie & Cui, 2011) into coupled WRF-CLM4 simulations, and further investigated the irrigation effects on the Indian Summer Monsoon precipitation. Joseph and Ghosh (2023) developed an improved irrigation module for the Variable Infiltration Capacity model by adopting similar methodology to parameterize the surface water and energy balance processes modified by ponded paddy fields. But neither of them considered the dynamic crop-growth simulation for paddy rice, and their model validation part neglected surface heat fluxes which are essential indices of the land-atmosphere interaction.

Noah-MP is a well-known land surface scheme that was augmented from the legacy Noah LSM with improved physics and multi-parameterization options (Niu et al., 2011; Yang et al., 2011). Liu et al. (2016) introduced dynamic corn and soybean simulations into Noah-MP and the enhanced model (i.e., Noah-MP-Crop used in this study) showed improved performance in simulating leaf area index (LAI), crop biomass and surface heat fluxes. A dynamic irrigation scheme was incorporated into Noah-MP to represent crop-specific irrigation applications over corn and soybean fields (Xu et al., 2019). Zhang et al. (2020) employed Noah-MP with the above-mentioned dynamic crop and irrigation schemes, and suggested that joint modeling of crop growth and irrigation improved the simulation of crop yield over irrigated regions in the central United States. However, neither of the dynamic crop-growth and irrigation schemes implemented in Noah-MP was applicable for paddy rice. Additionally, we selected Noah-MP to investigate the influence of paddy water over flooded rice fields because it has been implemented in Weather Research and Forecasting (WRF), which provides multioptions for atmospheric physical processes and has the largest number of users in the atmospheric science field. Therefore, by utilizing field observations from two paddy rice sites in Japan, the objectives of this study were to: (a) incorporate the effect of shallow paddy water into the widely used Noah-MP LSM coupled with a dynamic crop scheme (i.e., Noah-MP-Crop; Liu et al., 2016); (b) evaluate the performance of the improved Noah-MP-Crop model with paddy water against in situ observations; and (c) identify the roles of some essential parameters in capturing surface energy exchanges over rice paddies by a series of sensitivity tests. This work aimed to improve the representation of land-atmosphere interactions over flooded rice paddies in Noah-MP with the ultimate goal of coupling it with WRF. It is necessary for future modeling studies with WRF to assess how flooded rice paddies affect the atmosphere and to improve weather and climate predictions via enhanced paddy rice representation, especially in Asian countries.

The remainder of this paper was organized as follows. Section 2 presented field observations for model initialization and verification and depicted the model modifications due to the presence of paddy water in detail, as well as a set of numerical experiments conducted in this study. The calibration of some crop parameters using field measurements was also presented in Section 2. Section 3.1 evaluated the improved model performance in capturing the features of small Bowen ratios over flooded rice fields against the observations and Section 3.2 investigated the roles of some essential parameters such as ground surface resistance, water depth, and specific leaf area (SLA) in land and crop modeling. Finally, conclusions and some discussions were provided in Section 4.

Methodology and Data Model Description

LSMs are the part of climate models that simulate processes occurring at the Earth's surface, and they solve the coupled fluxes of water, energy and carbon between the land surface and atmosphere in the context of direct and indirect human forcings and ecological dynamics (Fisher & Koven, 2020). As an augmented version of the legacy Noah LSM (Chen & Dudhia, 2001; Chen et al., 1996), Noah-MP utilizes multiple options for key land-atmosphere interaction processes to represent seasonal and annual cycles of vegetation, hydrology and snow (Liu et al., 2016; Niu et al., 2011; Yang et al., 2011). The validation against in situ and satellite data showed that compared to the legacy Noah LSM, Noah-MP significantly improved the simulation of surface heat fluxes, soil moisture, land skin temperature, runoff and snow (Cai et al., 2014; Niu et al., 2011; Yang et al., 2011). Some options for runoff and groundwater schemes are provided in Noah-MP, but they were not activated in this study due to the absence of observation for model evaluation. All numerical experiments in this study were conducted at the field scale using the offline version of Noah-MP (Niu et al., 2011) running within the High-Resolution Land Data Assimilation System (HRLDAS v3.9; Chen et al., 2007). HRLDAS was developed to provide coupled WRF simulations with finescale land-state initialization such as soil temperature and soil moisture, the latter of which plays an important role in the development of deep-convection but remains very difficult to obtain. The heart of the HRLDAS infrastructure employed here is the Noah-MP land model. Additionally, HRLDAS is a highly efficient and parallelized LSM driver, which can be executed at both field and regional scales and has been widely employed to test and evaluate LSMs (e.g., Gao et al., 2015; Li et al., 2018; Zhang et al., 2016).

This study utilized the following combination of the model physics options: dynamic crop model (Liu et al., 2016), Ball-Berry-type canopy stomatal resistance (Ball et al., 1987) and Noah-type factor for stomatal resistance (Niu et al., 2011), and Monin-Obukhov scheme for the surface exchange coefficient (Brutsaert, 1982). Croplands are important in land-atmosphere interactions and the Noah-MP model is capable of dynamically simulating crop growth in two ways. The first is for generic dynamic vegetation (Niu et al., 2011; Yang et al., 2011), which does not distinguish growth characteristics such as planting/harvesting date, growing season, and dry matter partitioning among different crop species. The second is the dynamic crop scheme (Liu et al., 2016) used in this study, which can differentiate corn and soybean in peak LAI values as well as the growing season length by crop-specific parameters prescribed in the look-up table. A brief description of the Noah-MP-Crop model was provided in Supporting Information S1; more details were documented in Niu et al. (2011) and Liu et al. (2016).

The Original Surface Energy Balance

The Noah-MP model separates the vegetation canopy layer from the ground surface and utilizes a “semitile” subgrid scheme to represent the land surface heterogeneity (Niu et al., 2011). The Supporting Information S1 provided a brief description of the semitile scheme. Here, we presented the general equation for the surface energy balance (see Equation 1 below) in the original version of Noah-MP, and the detailed calculation of each component can be found in Niu et al. (2011). [Image Omitted. See PDF]where $$\mathrm{R}\mathrm{n}$$, $$\mathrm{S}\mathrm{H}$$, $$\mathrm{L}\mathrm{H}$$, and $$\mathrm{G}\mathrm{s}$$ represent the net radiation (unit: W m⁻²), sensible heat (unit: W m⁻²), latent heat (unit: W m⁻²), and the heat flux stored in soil (unit: W m⁻²), respectively. In the presence of paddy water, one main difference in the surface energy balance was the introduction of the heat flux stored in paddy water.

The Modified Surface Energy Balance in the Presence of Paddy Water

To incorporate the effect of flooded rice paddies on land-atmosphere exchanges, shallow paddy water was parameterized as one integrated layer along with the top soil layer in Noah-MP. The key modification in the surface energy balance was indicated by Equations 2 and 3 below. [Image Omitted. See PDF] [Image Omitted. See PDF]where $$\mathrm{R}\mathrm{n}$$, $$\mathrm{S}\mathrm{H}$$, and $$\mathrm{L}\mathrm{H}$$ are the net radiation (unit: W m⁻²), sensible heat (unit: W m⁻²), and latent heat (unit: W m⁻²), respectively. The total ground heat flux (i.e., $$\mathrm{G}\mathrm{t}\mathrm{o}\mathrm{t}$$ in Equation 2, which was equal to the heat flux stored in soil when no paddy water was considered in the original model) was divided into two terms (see Equation 3): the heat flux stored in soil ($$\mathrm{G}\mathrm{s}$$) and paddy water ($$\mathrm{G}\mathrm{w}\mathrm{a}\mathrm{t}$$). The two new terms of heat flux storage, that is, $$\mathrm{G}\mathrm{s}$$ and $$\mathrm{G}\mathrm{w}\mathrm{a}\mathrm{t}$$, were calculated as follows: [Image Omitted. See PDF] [Image Omitted. See PDF]

Here, $${\lambda }_{\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{k}}$$ is the bulk thermal conductivity of the top layer of paddy soil and was computed as a weighted average of the thermal conductivity for paddy water and the top soil. The model employs four-layer soil structure, and $${\increment}z$$ indicates the distance between ground surface and the first model level of soil columns, which is half of the depth of the top soil layer ($${\increment}z$$ = 0.05 m here). $${T}_{g}$$, $${T}_{\mathrm{s}\mathrm{l}1}$$, and $${T}_{w}$$ represent the ground surface, top-layer soil, and paddy water temperatures, respectively. $${C}_{w}$$ is the volumetric specific capacity of water ($${C}_{w}=4.188\,\mathrm{M}\mathrm{J}\,{\mathrm{m}}^{-3}{\mathrm{K}}^{-1}$$), D_w is paddy water depth, and $$dt$$ is the model time step. In this study, it was assumed that $${T}_{w}$$ is equal to the ground surface temperature, that is, $${T}_{w}={T}_{g}$$, as in the paddy rice model of Maruyama et al. (2017). To simplify the $$\mathrm{G}\mathrm{w}\mathrm{a}\mathrm{t}$$ calculation, this study ignored the temperature differences among the inflow irrigation water, paddy water, and outflow drainage water. The heat fluxes associated with the downward movement of paddy water into the soil layer were also ignored.

One unique aspect of the incorporation of paddy water in this study is that the water depth varies with the developmental stages of paddy rice, which can be read from observations or prescribed in the look-up table. This variable water depth affects the new term of heat flux stored in paddy water (i.e., $$\mathrm{G}\mathrm{w}\mathrm{a}\mathrm{t}$$ in Equation 5), and consequently, the surface energy budgets. Additionally, the surface albedo over flooded rice paddies was modified to 0.1 and a roughness length of 0.001 m (Maruyama & Kuwagata, 2010) was used to consider the presence of the paddy water surface. The variation in roughness length with paddy water depth (which was considered in Maruyama et al., 2017) was ignored in this study to simplify the model. The decrease in ground surface resistance to evaporation or sublimation due to the presence of paddy water (see Section 3.2) was considered by using a tenth of the adjusted value for wet soil based on Sellers et al. (1992).

Data for Model Input and Evaluation

HRLDAS is essentially an uncoupled land modeling system that utilizes a combination of observed or analyzed meteorological forcing to drive the Noah-MP LSM to simulate long-term evolution of land-state variables (Chen et al., 2007). To initialize and evaluate the original and modified versions of the HRLDAS/Noah-MP model, field observations from two different paddy rice sites (i.e., SAITO and SAGA) for the growing season of 2002 (i.e., from transplanting to harvesting; see Table 1) were employed in this study. Note that, the cropping seasons at the two sites were very different from each other, with paddy rice planted and matured about 3 months earlier at SAITO than at SAGA (Table 1). Both the SAITO (Miyazaki Plain) and SAGA (Chikushi Plain) sites were located in representative rice-growing areas that experienced humid temperate climates in Japan (Maruyama & Kuwagata, 2010). Table 1 shows the locations and dates of rice phenology at the two sites. The soil type for both rice paddy fields was gray soil. More details on the field experiments were documented by Maruyama and Kuwagata (2008).

Table 1 Locations and Dates of Rice Phenology for Two Crop Sites

Crop site	Location	Elevation (m)	Date of rice phenology
			Transplanting	Heading	Maturity	Harvesting
SAITO (early rice)	(32°06.0′N, 131°22.7′E)	11	March/18	June/15	July/25	July/20
SAGA (late rice)	(33°12.2′N, 130°16.8′E)	10	June/15	August/27	October/8	October/8

The forcing atmospheric variables required for HRLDAS/Noah-MP simulations were surface air temperature, air pressure, specific humidity, wind, downward short- and long-wave radiation, and precipitation. Some variables in the growing season, that is, hourly surface air temperature, relative humidity, wind speed at 2 m above ground level (interpolated from the observed level), and downward short- and long-wave radiation fluxes above the canopy, were available from field measurements. In the absence of field measurements during non-growing season, all of the forcing variables were employed from the nearest Automated Meteorological Data Acquisition System (AMeDAS) site or the nearest meteorological observatory of the Japan Meteorological Agency (JMA) to construct forcing files for a complete year for model spin-up. Table 2 lists the locations of the nearest stations and their data availabilities. According to a comparison among field measurements, AMeDAS and observatory data in the growing season (Table S3 in Supporting Information S1), surface air temperature and wind speed were from the nearest AMeDAS site for SAITO to fill the forcing data during non-growing season, while for the SAGA site, 2-m air temperature from the nearest AMeDAS site was employed. Additionally, surface air pressure was from the nearest observatory site, which was the only source for it. Precipitation was from the nearest AMeDAS site, and this choice was somewhat arbitrary since no field measurements were available for data evaluation.

Table 2 Nearest Site Information and Data Availability

Crop site	Nearest AMeDAS or observatory site	Location	Elevation (m)	Data availability
SAITO (early rice)	Saito (AMeDAS)	(32°05.4′N, 131°23.5′E)	11	Temperature, Wind, and Precipitation
	Miyazaki (Observatory)	(31°56.1′N, 131°25.0′E)	9	All forcing variables
SAGA (late rice)	Shiraishi (AMeDAS)	(33°10.8′N, 130°08.3′E)	4	Temperature, Wind, and Precipitation
	Saga (Observatory)	(33°15.8′N, 130°18.4′E)	6	All forcing variables

Crop parameters in the dynamic crop module (as listed in Table 3) were generally from previous literature (e.g., Bouman et al., 2001; Masutomi, Ono, Mano, et al., 2016; Masutomi, Ono, Takimoto, et al., 2016) and from the calibration by using observed biomass records (i.e., dry matter weight of leaf, stem, root, and panicle) at the SAGA site. Notably, these parameter values regarding phenology at the two sites were different from each other (the reason was discussed later in Section 4). Phenological development plays a critical role in simulating crop growth. First, we showed the estimation of phenological parameters that determine rice growth stages. Temperature is usually utilized by many crop models to indicate growth stages (e.g., Bouman et al., 2001; Li et al., 2017), as it is the main driver of phenological development (van Keulen et al., 1982). Similarly, accumulated growing degree days (GDD) was calculated in Noah-MP-Crop as a heat unit to differentiate phenology stages along with crop-specific thresholds (i.e., GDDS1-GDDS5 in Table 3). The minimum and maximum temperatures for GDD accumulation in rice (i.e., GDDTBASE and GDDTCUT in Table 3) were 8°C and 40°C, respectively (Bouman et al., 2001). The growing season in the model was generally determined by planting and harvesting dates (i.e., PLTDAY and HSDAY in Table 3) from site observations. Note that the growing season of paddy rice in this study started from transplanting because of the absence of field measurements before it. The key developmental stages for rice are emergence (not considered in this study), transplanting, panicle initiation, heading or flowering, and maturity (Bouman et al., 2001). To distinguish among them, the GDD thresholds (i.e., GDDS1-GDDS5 in Table 3) were estimated by following the values of a normalized heat unit (i.e., the ratio of the accumulated GDD at the current time step to total GDD accumulation until maturity) in MATCRO-Rice (Masutomi, Ono, Mano, et al., 2016; Masutomi, Ono, Takimoto, et al., 2016). In this study, GDD accumulated after transplanting and reached 1488°C ·d and 1904°C ·d at harvesting (i.e., GDDS5 estimation; see Figure S1 in Supporting Information S1) for SAITO and SAGA, respectively. According to the normalized ratios at transplanting, heading, and maturity, as well as some critical changes in carbohydrate partitioning in Masutomi, Ono, Takimoto, et al. (2016), GDDS2-GDDS4 were further estimated for SAITO and SAGA (Table 3). Carbohydrate loss in leaves, stems, and roots due to turnover and senescence began after flowering, and the rates (i.e., DILE_FC, LF_OVRC, ST_OVRC, and RT_OVRC in Table 3) were estimated using the biomass data for SAGA in 2012 (see Figure 1a). Next, the partitioning ratios of dry matter to various organs (i.e., leaf, stem, root, and grain), which are defined as functions of the development stages, were then parameterized by biomass data from SAGA (see Figure 1). At the beginning of the growing season, all the carbohydrates produced were partitioned into leaves (∼35%), stems (∼35%), and roots (∼30%). As the rice plant grew, more carbohydrates entered the leaves and stems, with partitioning ratios up to 0.45 (see Table 3). After heading, the dry weights of the leaves, stems, and roots started to decrease and most carbohydrates were partitioned into panicles. The partition ratios estimated here are close to those parameterized by Masutomi, Ono, Takimoto, et al. (2016).

Table 3 Some Crop Parameters Used for SAITO and SAGA

Parameter name	Description	Parameter value
		SAITO	SAGA
PLTDAY	(Trans)planting date (in Julian day format)	77	166
HSDAY	Harvesting date (in Julian day format)	201	281
GDDS1	Growing degree days (GDD) from transplanting to initial vegetative (°C·d)	15	19
GDDS2	GDD from transplanting to normal vegetative (°C·d)	542	318
GDDS3	GDD from transplanting to panicle initiation (°C·d)	697	1618
GDDS4	GDD from transplanting to flowering (°C·d)	1080	1808
GDDS5	GDD from transplanting to physical maturity (°C·d)	1488	1904
GDDTBASE	Minimum temperature for GDD accumulation (°C)	8
GDDTCUT	Maximum temperature for GDD accumulation (°C)	40
BIO2LAI	Leaf area per living leaf biomass (m2 g−1)	0.032
DILE_FC(PGSa)	Coefficient for leaf temperature stress death (s−1)	PGS = 6: 0.15
LF_OVRC(PGS)	Leaf turnover coefficient (s−1)	PGS = 6: 0.15
ST_OVRC(PGS)	Stem turnover coefficient (s−1)	PGS = 6: 0.2
RT_OVRC(PGS)	Root turnover coefficient (s−1)	PGS = 6: 0.1
LFPT(PGS)	Fraction of carbohydrate flux to leaf	PGS = 3: 0.35
		PGS = 4: 0.45
		PGS = 5: 0.2
STPT(PGS)	Fraction of carbohydrate flux to stem	PGS = 3: 0.35
		PGS = 4: 0.45
		PGS = 5: 0.4
RTPT(PGS)	Fraction of carbohydrate flux to root	PGS = 3: 0.3
		PGS = 4: 0.1
		PGS = 5: 0.1
GRAINPT(PGS)	Fraction of carbohydrate flux to grain	PGS = 5: 0.3
		PGS = 6: 1.0

^a PGS indicates the plant growth stage (see more details in Supporting Information S1 and Liu et al., 2016), and some parameters, such as the partitioning ratios of carbohydrates, are defined as functions of PGS.

LAI is an essential variable in cropland-atmosphere interactions and affects photosynthesis, leaf biomass, albedo, transpiration, and surface energy budgets, including sensible and latent heat fluxes (Liu et al., 2016). In Noah-MP-Crop, LAI is computed as the product of green leaf biomass and SLA, which varies with crop species and growth stage. However, owing to the lack of detailed observation data, a constant SLA value (i.e., BIO2LAI in Table 3) was prescribed throughout the entire growing season in Noah-MP-Crop for simplicity (Liu et al., 2016). The value of 0.032 m² g⁻¹ for both SAITO and SAGA in this study agreed with the 0.017–0.045 m² g⁻¹ range for rice reported in previous studies (Bouman et al., 2001; Masutomi, Ono, Mano, et al., 2016; Masutomi, Ono, Takimoto, et al., 2016). To further evaluate the role of SLA, two more sensitivity tests were performed using a minimum value of 0.017 m² g⁻¹ and a maximum value of 0.045 m² g⁻¹ (see Section 3.2).

Model verification involved surface energy fluxes (i.e., net radiation, ground heat flux, and sensible and latent heat fluxes) at hourly intervals. The main focus was the evaluation of sensible and latent heat fluxes, as one major purpose of this work was to improve the capture of small Bowen ratios over flooded rice fields. Additionally, the hourly water temperature was used to evaluate the performance of the incorporated paddy water module. Hourly soil temperatures at 2, 6, and 25 cm below the soil surface were also employed. A linear interpolation of the observed soil temperatures at 2 and 6 cm was utilized to estimate the top-layer soil temperature (i.e., 5 cm below the soil surface) for model evaluation. The growth of paddy rice modeled by the modified Noah-MP-Crop model was validated against the biological records of leaf and stem area indices (LAI and SAI) at intervals of 10–20 days. All model evaluation data were obtained from field measurements at both sites during the growing season.

Experimental Design

Table 4 documents a series of numerical simulations respectively conducted for SAITO and SAGA. All experiments here were performed from 1 March 2002 to 28 February 2003, with a 10-year spin-up period which ensured the modeled soil temperature and moisture reach an equilibrium state. Simulation results for the 2002 growing season with spun-up soil temperature and moisture were then employed in the following analysis.

Table 4 Experimental Design

Experiment name	Description	Key setting difference
		Paddy water and its depth	Crop parameters	Surface resistance
DEFLT	Control run with paddy water not considered and using default crop parameters.	No	Default	Not decreased
NOWAT	Same as DEFLT, but with calibrated crop parameters.	No	Calibrated	Not decreased
OBSWAT	Same as NOWAT, but with paddy water considered.	Yes/Observed	Calibrated	Decreased
RSWAT	Same as OBSWAT, but without surface resistance decreased.	Yes/Observed	Calibrated	Not decreased
FIXWAT	Same as OBSWAT, but with the paddy water depth fixed at 10 cm in the growing season.	Yes/Fixed	Calibrated	Decreased
SLAMAX (or SLAMIN)	Same as OBSWAT, but with the maximum value of 0.045 m2 g−1 (or the minimum value of 0.017 m2 g−1) for SLA.	Yes/Observed	Calibrated	Decreased

A brief description of the numerical experiments was provided as follows. The DEFLT run, which served as a reference run and utilized the default parameter configuration for the dynamic crop model in Noah-MP, was the starting point of this study. The NOWAT run was conducted using calibrated crop parameters for paddy rice from field measurements (see Table 3 and Section 2.2). In addition to the employment of calibrated crop parameters, the shallow paddy surface water layer was carefully considered (see Section 2.1) in the OBSWAT run, where the paddy water depth was updated daily by observations (see black dots in Figures 2a and 2b). The decreased ground surface resistance owing to the presence of paddy water in the OBSWAT run was set as a tenth of the adjusted value for wet soil based on Sellers et al. (1992). The difference in the modeled surface heat fluxes between the two simulations (i.e., NOWAT and OBSWAT) indicated quantitative improvements by reasonably incorporating the paddy water module into Noah-MP-Crop. Furthermore, a set of sensitivity experiments was performed to investigate the roles of some essential parameters, such as ground surface resistance and paddy water depth, in land and crop modeling over rice paddies. Instead of decreasing the surface resistance in the case of OBSWAT, the RSWAT run was conducted using option 1 for the surface resistance scheme (Sakaguchi and Zeng, 2009; as shown by dark green lines in Figures 2c and 2d), which was also employed in the cases of DEFLT and NOWAT. A constant value of 10 cm for paddy water depth was fixed throughout the growing season in the FIXWAT case. Additionally, the role of SLA (i.e., BIO2LAI in Table 3) was investigated by conducting two more sensitivity simulations (i.e., SLAMAX and SLAMIN), where a maximum value of 0.045 m² g⁻¹ and a minimum value of 0.017 m² g⁻¹ were employed.

Results

The predominant characteristic of submerged or even flooded rice fields is a small Bowen ratio (i.e., the ratio of sensible to latent heat flux) owing to the presence of shallow paddy water, as indicated by field observations for both SAITO and SAGA (see the black boxes in Figures 3e–3h). In the 2002 growing season, the observed daytime average sensible and latent heat fluxes at SAITO were ∼27 W m⁻² and ∼178 W m⁻², respectively. The Bowen ratio was even smaller at SAGA, where daytime averages of ∼3 W m⁻² and ∼293 W m⁻² were observed for sensible and latent heat fluxes, respectively. One major purpose of this work was to capture such a distinguished feature by incorporating paddy water into Noah-MP-Crop and employing calibrated crop parameters for rice. Section 3.1 presented the model performance evaluation against field measurements, and the roles of some key parameters were further investigated in Section 3.2.

Model Performance

We first presented the performances of the two reference simulations (i.e., DEFLT runs for SAITO and SAGA), which were the starting points of this study. The results show overestimated Bowen ratios at both sites with biases in surface sensible and latent heat fluxes that were more significant for the SAGA site (see Figures 3–5, and Figure S2, Tables S1 and S2 in Supporting Information S1). Compared with the observed average diurnal cycles (see the black boxes in Figures 3e–3h), both reference simulations generated overestimations in sensible heat and/or underestimations in latent heat, especially during daytime. Additionally, the SAGA biases were more remarkable than those at SAITO, where the average sensible heat at peak hours was overestimated by ∼30 W m⁻², while the modeled diurnal curve of latent heat agreed well with observations (see the blue lines in Figures 3e and 3f and Table S2 in Supporting Information S1). Similar biases were also shown in the modeled seasonal patterns, and the overestimation in sensible heat and the underestimation in latent heat were more significant at SAGA as well (see Figure S2 in Supporting Information S1). The root-mean-square errors (RMSEs) of the sensible and latent heat fluxes in the SAITO reference run were 39 and 49 W m⁻², respectively (Table S1 in Supporting Information S1). At SAGA, the DEFLT run produced a greater RMSE of 97 W m⁻² for sensible heat and up to 119 W m⁻² for latent heat (see Figure 4 and Table S1 in Supporting Information S1). Net radiation is the driver for surface heat and water exchanges between paddy rice fields and the atmosphere, and its diurnal cycles and seasonal patterns were both generally well reproduced by the reference runs (see the blue lines in Figures 3a and 3b; Figures S2a and S2b in Supporting Information S1) despite a ∼33 W m⁻² average overestimation during the daytime for SAITO (see Figure 5 and Table S2 in Supporting Information S1). Given a fairly good approximation of the ground heat flux at SAITO (Figure 3c), the sum of sensible and latent heat fluxes (i.e., $$\mathrm{S}\mathrm{H}+\mathrm{L}\mathrm{H}=\mathrm{R}\mathrm{n}-\mathrm{G}\mathrm{t}\mathrm{o}\mathrm{t}$$) was larger than that observed, which could partly explain the overestimation of latent heat with the incorporation of paddy water (see the red line in Figure 3g). Figure 6 further demonstrate that these model biases in surface heat fluxes mostly concentrated on daytime hours at both sites. On average, the reference run for SAGA produced an overestimation of ∼107 W m⁻² for sensible heat and an underestimation of ∼127 W m⁻² for latent heat during daytime throughout the growing season (see the black and blue boxes in Figures 6f and 6h, and Table S2 in Supporting Information S1). Specifically, these larger model biases in the daytime Bowen ratio for SAGA tended to be even more remarkable in the middle and late periods of the growing season (i.e., August–October; see Figures 6f and 6h), while greater overestimation occurred at the beginning of the growing season for SAITO (i.e., March–April; see Figures 6e and 6g).

Using calibrated crop parameters for paddy rice by field measurements greatly improved the capture of small Bowen ratios over flooded rice fields, especially at SAGA. For SAITO, an improved simulation of sensible heat flux was generated by the NOWAT run (see Figures 3e and 6e), with the RMSE decreasing by ∼9 W m⁻² (∼23% of that in the DEFLT run; see Figure 4 and Table S1 in Supporting Information S1). However, the daytime latent heat flux was overestimated by ∼12 W m⁻² (see Figures 3 and 6g, Figure 5 and Table S2 in Supporting Information S1), which could be associated with the greater net radiation than that observed (see Figures 3a and 6a, Figure 5, and Table S2 in Supporting Information S1). Compared with SAITO, the employment of calibrated crop parameters in the NOWAT run for SAGA showed more significant improvements in surface heat fluxes (see Figures 3f, 3h, 6f, and 6h and Figure S2 in Supporting Information S1). The RMSEs of both sensible and latent heat fluxes decreased by approximately 35%–40%, from 97 to 119 W m⁻² in the reference run to 59 and 79 W m⁻² in the NOWAT run, respectively (see Figure 4 and Table S1 in Supporting Information S1). During daytime hours, both the overestimation in sensible heat and the underestimation in latent heat decreased by ∼55 W m⁻² (approximately a 40%–50% decrease compared with the DEFLT run) at SAGA (see Figure 5 and Table S2 in Supporting Information S1). However, non-negligible gaps between the observed and modeled surface heat fluxes remained. For instance, the NOWAT run at SAGA still overestimated the daytime sensible heat by ∼52 W m⁻² while underestimating the daytime latent heat by ∼69 W m⁻² (see the red boxes in Figures 6f and 6h, Figure 5 and Table S2 in Supporting Information S1). For SAITO, the sensible heat during the daytime was also still overestimated by ∼19 W m⁻² (∼65% of the overestimation in the DEFLT run; see Figures 3e and 6e, Figure 5 and Table S2 in Supporting Information S1).

The Noah-MP-Crop performance in the surface heat flux simulation was most improved by combining the reasonable incorporation of paddy water and the calibration of crop parameters for rice, especially for SAGA (see Figures 3–6). Compared with the NOWAT run (or the reference run) for SAGA, the RMSEs of the sensible and latent heat fluxes in the OBSWAT run decreased by ∼4 W m⁻² (or ∼43 W m⁻²) and by ∼11 W m⁻² (or ∼50 W m⁻²; see Figure 4 and Table S1 in Supporting Information S1). For SAITO, simulation of the sensible heat flux was also slightly improved by the additional incorporation of paddy water, with the RMSE decreasing from ∼30 W m⁻² in the case of NOWAT to ∼26 W m⁻² in the OBSWAT run (Figure 4 and Table S1 in Supporting Information S1). Although the overestimation of latent heat flux worsened at SAITO by the OBSWAT run (Figure 4 and Table S1 in Supporting Information S1), the model performance in surface energy fluxes was better than that of the reference run in terms of overall RMSEs for both sites, especially sensible and latent heat fluxes (Figure 4 and Table S1 in Supporting Information S1). In terms of daytime averages, model improvements in surface heat fluxes were more significant at the SAGA site. Compared with the NOWAT run for SAGA, the additional incorporation of paddy water in the OBSWAT run further decreased the overestimation of sensible heat and underestimation of latent heat by ∼19 W m⁻² and ∼35 W m⁻² (∼18% and ∼27% of the bias in the reference run), respectively (see the right panels in Figures 3, 5 and 6, and Table S2 in Supporting Information S1).

In addition to remarkably improving the capture of small Bowen ratios over flooded rice fields, the enhanced Noah-MP-Crop model was also capable of simulating the growth characteristics of paddy rice. For both SAITO and SAGA, LAI and SAI were well reproduced by employing calibrated crop parameters for paddy rice in the NOWAT runs, whereas they were greatly underestimated in the reference cases (Figure 7). The further incorporation of paddy water in the OBSWAT runs changed them slightly (see the red lines in Figure 7), indicating the significance and necessity of crop parameter calibration in both crop and land modeling. Despite some cold biases, the general trend of water and soil temperatures throughout the growing season was captured by the OBSWAT runs (Figure 8), which demonstrated the validity of the paddy water module. For SAGA, the overestimation of soil temperatures by the reference run (see Figures 8d and 8f) was significantly mitigated by the calibration of crop parameters for rice, as well as the incorporation of paddy water, especially in the middle and late periods of the growing season.

These results show that by using realistic crop parameters for paddy rice and by incorporating shallow paddy water, the Noah-MP-Crop model performed much better at capturing the distinct features of small Bowen ratios over flooded rice fields; however, the model improvements at the two sites slightly differed from each other. For SAITO, despite an overestimation in latent heat, the combination of calibrating rice crop parameters and incorporating paddy water helped mitigate the overestimation in sensible heat, with the RMSE decreasing from ∼39 W m⁻² in the reference run to ∼26 W m⁻² in the OBSWAT run. The model performance was more significantly improved at SAGA, and the RMSEs of the sensible and latent heat fluxes in the OBSWAT run decreased by ∼43 and ∼51 W m⁻², respectively, indicating an improvement of more than 40% compared with the reference case. Additionally, the growth characteristics of paddy rice were better reproduced at both sites as well as soil and water temperatures.

Roles of Some Key Parameters

Some sensitivity tests were performed in this section to further identify the respective roles of some key parameters, including ground surface resistance, water depth, and SLA (i.e., BIO2LAI in Table 3).

Decreased ground surface resistance, which was due to the presence of paddy surface water, played an essential role in better capturing the features of small Bowen ratios, especially at SAGA. Compared with the RSWAT runs, the decrease in surface resistance in the OBSWAT runs (see the dark green lines in Figures 2c and 2d) helped mitigate the overestimation of sensible heat and/or the underestimation of latent heat (see Figures 5 and 9, and Table S2 in Supporting Information S1). For instance, the daytime average sensible heat at SAGA was overestimated by ∼53 W m⁻² in the RSWAT run, which was ∼20 W m⁻² (∼60%) higher than that in the OBSWAT case (see Figures 5b and 9f, and Table S2 in Supporting Information S1). Meanwhile, the underestimation of the daytime latent heat in the RSWAT run reached up to ∼70 W m⁻², which was ∼35 W m⁻² (∼100%) larger than that of the OBSWAT run (see Figures 5b and 9h, and Table S2 in Supporting Information S1). For SAITO, the unchanged surface resistance (see the dark green line in Figure 2c) also generated a larger overestimation in sensible heat during the daytime, from ∼4 W m⁻² in the OBSWAT run to ∼22 W m⁻² in the RSWAT run (see Figures 5a and 9e, and Table S2 in Supporting Information S1).

Water depth, a newly incorporated variable in the paddy water module, probably had a nonlinear influence on the surface heat flux simulations (see modifications in surface energy budgets in Section 2.1.2). Although the observed water depth was shallower than the 10 cm fixed value in the FIXWAT run throughout the growing season (Figures 2a and 2b), a similarly shallow layer of paddy water in the OBSWAT case did not indicate a decrease in latent heat flux (see Figures 5, 9g, and 9h; Table S2 in Supporting Information S1). Conversely, the daytime average latent heat modeled by the FIXWAT run was ∼9 W m⁻² smaller than that in the OBSWAT case for both SAITO and SAGA (see Figure 5 and Table S2 in Supporting Information S1). This might be due to an increased heat capacity of paddy water at deeper depths, which resulted in cooler water temperatures during the daytime (not shown) and reduced evaporation from the water surface. Additionally, as the water depth changed from the OBSWAT to the FIXWAT runs, the first affected variable was the ground heat flux, which, increased by ∼12 W m⁻² on average during the daytime (see Figures 5, Figures 9c and 9d and Table S2 in Supporting Information S1). Considering the minimal changes in net radiation (see Figure 5 and Table S2 in Supporting Information S1), the sensible and latent heat fluxes were both decreased in the FIXWAT run, resulting in a more significant underestimation of the latent heat at SAGA (see Figure 5 and Table S2 in Supporting Information S1).

These two parameters have one thing in common: they modified the surface energy fluxes with little change in crop growth (not shown). In fact, some parameters associated with the dynamic crop module (e.g., SLA investigated here) also played non-negligible roles in the simulation of surface energy fluxes by affecting paddy rice growth. The results of the sensitivity tests (i.e., SLAMAX and SLAMIN) showed that the 0.017–0.045 m² g⁻¹ range for SLA resulted in a peak difference of ∼6 m² m⁻² in LAI and ∼1.5 m² m⁻² in SAI for both SAITO and SAGA (see Figure 10). The modeled surface heat fluxes changed accordingly. For instance, the overestimation of daytime sensible heat flux (or the underestimation of daytime latent heat flux) at SAGA was increased by ∼27 W m⁻² (or ∼20 W m⁻²) when the SLA changed from 0.045 to 0.017 m² g⁻¹ (see Figures 5, Figures 9f and 9h, and Table S2 in Supporting Information S1). Additionally, such a decrease in SLA also produced an increase of ∼15 W m⁻² in daytime net radiation, which was non-negligible compared to the cases of ground surface resistance and water depth.

Generally, the parameters influencing the simulation of surface heat fluxes over flooded rice fields could be divided into two categories: without (e.g., surface resistance and water depth) or with (e.g., SLA) remarkable changes in the growth characteristics of paddy rice. Decreased ground surface resistance owing to the presence of paddy water was the most crucial factor for better capturing the features of small Bowen ratios over flooded rice fields. Compared with the RSWAT run for SAGA, the decreased ground surface resistance in the OBSWAT case helped mitigate the overestimation of sensible heat by ∼20 W m⁻² and the underestimation of latent heat by ∼35 W m⁻² during the daytime. Water depth possibly had a nonlinear influence on the surface heat flux simulations and more realistic data helped capture the features of small Bowen ratios over flooded rice fields. Crop parameters, such as SLA, produced significant changes in LAI, and consequently, surface energy fluxes.

Conclusions and Discussion

Rice, the major crop type in Asia, is typically grown in flooded paddy fields and requires large amounts of irrigation water throughout the growing season. The presence of shallow paddy water in rice fields can modify surface water and energy budgets and affect weather and climate via changes in surface properties. However, the effect of flooded rice paddies has not been considered in many land surface schemes in climate models such as Noah-MP, which is coupled into the WRF model and is widely used in simulating land-atmosphere interactions. To improve the model performance of Noah-MP-Crop over rice fields, shallow paddy surface water was incorporated in this study as an integrated layer along with the top-layer soil; for simplicity, the water temperature was assumed to equal the ground surface temperature as in Maruyama et al. (2017). Meteorological and biological measurements from two paddy rice sites in Japan, SAITO (early rice) and SAGA (late rice), were employed to initialize and evaluate the enhanced Noah-MP-Crop model. One unique aspect of this study is the identification of the roles of some key parameters for capturing the surface heat fluxes over flooded rice fields.

Compared with other dryland crops such as corn and soybean (Chen et al., 2018; Liu et al., 2016; Xu et al., 2019), the predominant characteristic of submerged or even flooded rice fields is a small Bowen ratio owing to the presence of shallow paddy water. For instance, the daytime averages of ∼27 and ∼178 W m⁻² were respectively observed for sensible and latent heat fluxes at SAITO. Generally, the performance of Noah-MP-Crop in reproducing the distinct features of small Bowen ratios over flooded rice fields was significantly improved by combining the incorporation of paddy water and the calibration of crop parameters for rice; however, the model improvements at the two sites slightly differed from each other. For SAITO, the overestimation in sensible heat was mitigated and the RMSE decreased from ∼39 W m⁻² in the reference run to ∼26 W m⁻² in the OBSWAT run. The original Noah-MP-Crop model with default crop parameters produced more remarkable biases at SAGA, and the RMSEs of the sensible and latent heat fluxes decreased from ∼97 to ∼119 W m⁻² in the reference run to ∼54 and ∼68 W m⁻², respectively. Because of using hourly data for evaluation, these RMSEs in the OBSWAT runs were larger than those in Masutomi, Ono, Takimoto, et al. (2016) and Maruyama and Kuwagata (2010), where daily values were employed and average RMSEs of less than 20 W m⁻² for both sensible and latent heat fluxes were reported. Importantly, this study agreed with Masutomi, Ono, Takimoto, et al. (2016) in that the flooded surfaces and irrigation largely affected the simulation of sensible and latent heat fluxes over rice paddies.

Crop parameter calibration was essential for improving the capture of small Bowen ratios at both sites. Using calibrated crop parameters for rice alone decreased the RMSEs of the sensible and latent heat fluxes at SAITO by ∼9 W m⁻² (∼20% of that in the reference run). The model improvements were more remarkable at SAGA, where the RMSEs of both sensible and latent heat fluxes decreased by ∼40 W m⁻² (approximately 34%–40% of that in the reference run). During the daytime, the overestimation and underestimation in sensible and latent heat, respectively, decreased by ∼55 W m⁻² (an approximate decrease of 35%–50%) in the SAGA reference run. The realistic incorporation of paddy water in the OBSWAT run further improved the simulation of surface heat fluxes at SAGA, generating an additional decrease of ∼19 W m⁻² in the overestimation of sensible heat and ∼35 W m⁻² in the underestimation of latent heat during the daytime. One thing to note here is that those calibrated parameter values regarding phenology (i.e., GDD thresholds utilized to distinguish the developmental stages) at SAITO and SAGA were different from each other (see Table 3). This could be associated with the daylength effect on rice growth rates not yet considered in Noah-MP (Bouman et al., 2001; Maruyama & Kuwagata, 2010). Rice is a short-day plant, and daylength affects the induction of flowering along with temperature. Suboptimal photoperiods less or more than optimum photoperiod could result in a longer photoperiod-sensitive phase (Bouman et al., 2001; Maruyama & Kuwagata, 2010), which was from the end of the basic vegetative phase to panicle initiation (approximately corresponding to the period of GDDS2-GDDS3 in this study). Compared with SAITO, paddy rice at SAGA was grown in hot months, during which daylength was very likely to be suboptimal rather than optimum. As a result, the photoperiod-sensitive phase at the SAGA site was longer, which was indicated by a bigger spread between GDDS2 and GDDS3.

The roles of some key parameters in surface heat flux simulations were further investigated using sensitivity tests, and the results show that decreased ground surface resistance played a major role in helping the paddy water module better capture the features of small Bowen ratios. For instance, compared with the OBSWAT simulation for SAGA, the unchanged surface resistance in the RSWAT case resulted in an overestimation of sensible heat by ∼20 W m⁻² and an underestimation of latent heat by ∼34 W m⁻² during the daytime. Paddy water depth possibly affected the simulation of surface heat fluxes in a nonlinear manner. For SAGA, the employment of the observed water depth, which was shallower than the 10 cm fixed value in the FIXWAT run, generated a decrease of ∼11 W m⁻² in the underestimation of the daytime latent heat flux. Crop parameters such as SLA were different from ground surface resistance and water depth in that SLA modified the surface heat fluxes via remarkable changes in growth of paddy rice growth (e.g., LAI). The variation in SLA from 0.017 to 0.045 m² g⁻¹ at SAGA caused a decrease of ∼27 W m⁻² in sensible heat and an increase of ∼20 W m⁻² in latent heat during the daytime, indicating the necessity of crop parameter calibration.

The results of this study are very useful and beneficial for improving the regional-scale simulations of surface heat fluxes and enhancing the understanding of land-atmosphere interactions over flooded rice fields, especially in Asian countries. Also, the paddy water module in this study can serve as a useful reference for the incorporation of paddy water into other land and climate models. However, there are some limitations or sources of model uncertainty in this study. We incorporated paddy water as one integrated layer along with the top soil layer, and the paddy water depth was prescribed by observation or in the look-up table. In this case, the influence of paddy water was represented without the implementation of a dynamic irrigation scheme employed in Devanand et al. (2019) and Joseph and Ghosh (2023), and the paddy water depth in the above two studies was updated every time step based the water balance of precipitation reaching the ground, irrigation applied, evaporation from the ponded water and infiltration into the soil. Hence, we ignored some sink terms in the water balance such as runoff and infiltration as well as the source of irrigation water, which is very important for regional-scale hydroclimate simulations. Future efforts should be directed to the incorporation of the modified water balance process when transitioning the land and crop modeling over flooded rice paddies from the field to regional scale. Additionally, due to the absence of field measurements, the atmospheric forcing conditions were constructed by using observations from the nearest observatory or AMeDAS site during non-growing season, which may also contribute to the model uncertainty. For instance, despite the smaller errors to field measurements at SAGA, surface air temperature from the nearest AMeDAS site was on average ∼1°C higher than observation in the growing season of 2002 (see Table S3 in Supporting Information S1). Further, those GDD thresholds used to distinguish the SAITO and SAGA developmental stages (Table 3) were different from each other, which should be normalized by incorporating the daylength effect on GDD accumulation (Bouman et al., 2001). Most importantly, future work will focus on the influence of flooded rice fields on regional weather and climate using the enhanced WRF/Noah-MP-Crop model.

Acknowledgments

The authors would like to thank three anonymous reviewers for extremely helpful feedback on this manuscript. This work was funded by JSPS KAKENHI (Grant 19H03084). The author Xiaoyu Xu also would like to acknowledge the support from the Scientific Research Foundation for New Teachers of Nanjing University of Aeronautics and Astronautics (Grant YAT21001) and the National Natural Science Foundation of China (Grant 42105170).

Data Availability Statement

The numerical simulations were mostly performed on the computing facilities in the Center for Computational Sciences at the University of Tsukuba. Field measurements for the SAITO and SAGA sites were employed in the creation of this manuscript (Maruyama, 2021). The major code modifications and model outputs were attributed to Xu et al. (2022). Additionally, surface observations for the nearest AMeDAS or meteorological observatory stations were obtained from the Japan Meteorological Agency website and were also included in the Zenodo repository (Xu et al., 2022). All the figures in this study were made with the NCAR Command Language (NCL) version 6.6.2 (UCAR/NCAR/CISL/TDD, 2019).