Heat-Mitigation Effects of Irrigated Rice-Paddy Fields Under Changing Atmospheric Carbon Dioxide Based on a Coupled Atmosphere and Crop Energy-Balance Model

Known as the heat-mitigation effect, irrigated rice-paddy fields distribute a large fraction of their received energy to the latent heat during the growing season. The present hypothesis is that increased atmospheric CO2 concentration decreases the stomatal conductance of rice plants and increases the air temperature by means of an increased sensible heat flux. To test this hypothesis, a coupled regional atmospheric and crop energy-balance model is developed and applied to a 300 × 300 km2 region in Japan. Downscaling meteorological variables from grid-mean values of mixed land use (3 × 3 km2) generates realistic typical diurnal cycles of air temperature in rice paddies and adjacent residential areas. The model simulation shows that, on a typical sunny day in summer, doubling the CO2 concentration increases the daily maximum grid-mean air temperature, particularly where rice paddies are present, by up to 0.7 °C. This CO2 effect on the grid-mean air temperature is approximately half the effect of the reduction in rice-paddy area that is postulated to occur on a time scale similar to that of the atmospheric CO2 change. However, within the internal atmospheric boundary layer of the rice paddies, the CO2 effect on the air temperature (+ 0.44 °C) still exceeds the effects of the land-use change (+ 0.11 °C). These results show a potentially important interplay of plant physiological responses regarding atmospheric CO2 in the heat-mitigation effect of rice-paddy fields under a changing climate.


Introduction
The heat-mitigation effect of irrigated crop fields has been recognized for a long time (Holmes 1970). In response to ongoing climate changes, an interest has emerged regarding whether and how this effect enhances or mitigates environmental changes (Dong et al. 2016;Nocco et al. 2019;Du et al. 2019;Yang et al. 2020). In particular, rice paddies comprise the typical landscapes in many parts of eastern Asia, and their role in heat mitigation may be of interest.
Previous studies have reported that the vicinity of rice paddies is characterized by cooler air compared with the other surrounding land that does not have rice paddies (Oue et al. 1994;Wakiyama 2007;Dong et al. 2016;Kuwagata et al. 2018); this is because rice paddies distribute a large fraction of their received energy to latent heat (Oue 2005;Ikawa et al. 2017). Recent studies have reported that changing the areal extents of rice paddies affects the regional climate (e.g., Dong et al. 2016). However, to the best of our knowledge, no study has investigated how the heat-mitigation effect of rice paddies is altered by the physiological response of rice plants under a changing climate.
Increases in atmospheric carbon dioxide [hereinafter, this concentration in the background atmosphere is denoted as [CO 2 ] decrease the latent heat flux and increase the sensible heat flux via the stomatal closure of rice plants (Yoshimoto et al. 2005). Vegetation responses to increasing [CO 2 ] likely alter the energy balance and subsequently influence the state of the atmosphere, which in turn affects land-surface processes (Vilà-Guerau de Arellano et al. 2012Arellano et al. , 2014. In an earlier study (Ikawa et al. 2018), we used a simple one-dimensional clear-sky atmospheric-boundary-layer (ABL) model to show that the air temperature and humidity in the atmosphere over a rice paddy are influenced by different [CO 2 ] and different cultivars. However, meteorological variables over lowland rice-paddy areas in Japan are also influenced by various other processes at different spatial and temporal scales, such as advection from the ocean (Kuwagata et al. 1994), thereby making it difficult to perform realistic simulations using a one-dimensional model. It is important to quantify temperature environments precisely based on a realistic diurnal variation for assessing the heat stress on rice grains because the susceptibility of crops to high temperature may vary throughout a day (Matsui and Hasegawa 2019).
We suggest that a combination of a regional atmospheric model and rice-crop energybalance model could realistically simulate the processes relevant to the heat-mitigation effects of rice paddies on different spatial scales. The regional atmospheric model would simulate the spatial variations of air temperature due to different land use and topography (Watanabe and Shimoyama 2015), whereas the crop energy-balance model would simulate ecophysiological responses to an in situ environment (Maruyama and Kuwagata 2010;Ikawa et al. 2018). While community models such as the Weather Research and Forecasting model (e.g., Liu et al. 2016) may also be useful tools, the heat-mitigation effect would be specific to crop species , and we are interested in quantifying the effects of ecophysiological processes that occur in rice paddies on ABL conditions. The present study addresses the following specific research questions: (1) To what extent do the stomatal responses of rice plants to different levels of [CO 2 ] affect the air temperature of rice-paddy areas in Japan? (2) How is the change in air temperature due to the stomatal response compared with the impact of land-use change? The objectives corresponding to these questions are (i) to couple a regional atmospheric model and a rice-crop energy-balance model and simulate the air temperature at different levels of [CO 2 ], and (ii) to run the model with different areal fractions of rice paddies under current [CO 2 ] and compare the results with those obtained in objective (i).

Methodology
The rice-crop energy model is coupled with the atmospheric model as follows: the atmospheric model uses the surface fluxes estimated by the rice-crop energy-balance model for the rice paddies, and the crop energy-balance model uses the meteorological variables estimated by the atmospheric model. Each model is described below.

Regional Atmospheric Model
A hydrostatic mesoscale atmospheric model with the Boussinesq approximation is used to simulate the flow and scalar transport within the ABL. The model is based on that of Kimura and Arakawa (1983) but modified to allow air density to change vertically (Watanabe and Shimoyama 2015). The simulation domain is assumed to be static with no geostrophic force; the other settings are listed in Table 1. Using Einstein's summation notation (i, j = 1, 2, and 3 unless stated otherwise), we define x i as a longitudinal coordinate (x, positive towards the east), a lateral coordinate (y, positive towards the north), and a terrainfollowing height ( z * ), respectively. The basic equations are the continuity equation and the equations of motion. All variables are Reynolds averaged unless otherwise noted.
The continuity equation is where h is the vertical distance between the atmosphere top ( z TOP , defined at 6 km) and the terrain height ( z G ), and ρ 0 (kg m −3 ) is the time-mean air density and is assumed to be dependent on height only. The velocity components u 1 , u 2 , and u 3 are the longitudinal (u), lateral (v), and vertical ( w * ) components, respectively. The longitudinal and lateral components (i = 1, 2) are governed by where t is time, ε ij3 is the alternating unit tensor, f c is the Coriolis parameter (s −1 ), c p (J kg −1 K −1 ) is the specific heat of dry air at constant pressure, Θ (K) is the reference potential air temperature of the calculation domain, π′ is the spatial departure of the Exner function (π = T Θ −1 , where we define T as the grid-mean air temperature) from its reference Π, which varies with height only (dΠ/dz = − g[c p Θ] −1 ), g (m s −2 ) is the acceleration due to gravity, θ′ (K) is the spatial departure of the potential air temperature from Θ, and K M (m 2 s −1 ) is the vertical eddy viscosity obtained from a level-2 turbulence-closure model Yamada 1974, 1982). The turbulence-closure model is also used to evaluate the vertical eddy diffusivity K H (m 2 s −1 ) for scalars, as used in Eq. 4. The advection terms are approximated by a third-order finite-difference scheme, and a sixth-order horizontal filter is applied to reduce shortwave numerical oscillations. Using the hydrostatic assumption, the vertical velocity component ( u 3 = w * ) is diagnosed from the continuity equation (Eq. 1), and the Exner function is calculated from the hydrostatic equation with � = 0 at the top boundary. Note that if c p is assumed to be spatially constant, then it does not need to be specified because it is cancelled upon being substituted into the first term on the right-hand side of Eq. 2. Finally, the equation governing a scalar χ [the gridmean potential air temperature θ (K), specific humidity q (g kg −1 ), and CO 2 mixing ratio (g kg −1 )] transport is The CO 2 mixing ratio is converted to atmospheric CO 2 concentration, C (μmol mol −1 ), using the molar mass of CO 2 and dry air.

Land-Surface Model
The subgrid parametrization for energy fluxes over mixed land surfaces is realized using the area-weighted average of the fluxes calculated for six land types (rice paddy, city, other crops, grass, forest, and water surface) within a single grid (Fig. 1). Unlike a parametrization for a bulk land surface, the present method allows the use of parameters determined specifically for each land-use category (Kimura 1989;Avissar and Pielke 1989). Downscaling to estimate the meteorological variables within the internal boundary layer (IBL) over each land-use category is performed by means of postprocessing, as explained in Sect. 2.5.

Rice-Crop Energy-Balance Model
The rice-crop energy-balance model is based on the model of Ikawa et al. (2018), which comprises an energy-balance submodel (Watanabe 1994;Maruyama and Kuwagata 2010;Maruyama et al. 2017) and a canopy-photosynthesis submodel (de Pury and Farquhar 1997). The energy fluxes are calculated based on a so-called double-source model. For example, the latent heat flux λE (W m −2 ) is the sum of the latent heat fluxes from the canopy (λE c ) and ground (λE g ) as follows Within IBL U(z * 2 ), θ(z * 2 ), q(z * 2 ), C(z * 2 ) Δz *2 Above IBL T s = const.

Fig. 1
Mixed land surfaces are aggregated within each surface grid of the regional atmospheric model following the procedure of Kimura (1989), after which downscaling of the meteorological variables within the internal boundary layer (IBL) is performed. The sensible heat flux (H) and latent heat flux (λE) from each land-use category are aggregated based on their areal fraction and used as the boundary condition of the atmospheric model. The CO 2 flux (F c ) is assigned only for rice paddies in this study. The ground heat flux (G) and surface temperature (T s ) are determined for each land-use category for each surface grid. The wind speed (U), potential air temperature (θ), specific humidity (q), and atmospheric CO 2 concentration (C) within each grid are assumed to represent the aggregated land surface, and these meteorological variables within the IBL over each land-use category (denoted by the subscript 'a') are calculated during postprocessing where λ (J kg −1 ) is the latent heat of vaporization, z * 1 is the height at the centre of the lowest grid layer of the atmospheric model (≈ 15.5 m), U (m s −1 ) is the wind speed, T c (°C) is the canopy temperature, T g (°C) is the ground temperature, and C Ec and C Eg are the vapour transfer coefficients for the canopy and the ground, respectively, the former being a function of the average stomatal conductance g sw (mol m −2 s −1 ) for water vapour in the canopy. Note that herein we refer to T c as the effective aerodynamic temperature for sensible heat exchange. The sensible heat flux H (W m −2 ) is calculated in a similar manner and is set proportional to a product of the temperature gradient and the transfer coefficient.
Here, we approximate the average stomatal conductance, g sw by a parallel-resistance formula, which calculates the weighted mean value from g sw of sunlit (g sw,sun ) and shaded (g sw,shd ) leaves based on the photosynthesis and stomatal conductance submodel (FvCBL model) (Farquhar et al. 1980;Leuning 1995) as follows where L c is the leaf-area index (LAI), with L c,sun and L c,shd being the LAI for sunlit and shaded leaves, respectively, Q c,sun and Q c,shd (μmol m −2 s −1 ) are the absorbed photosynthetic photon flux densities by sunlit and shaded leaves, respectively, and D c (kPa) is the vapour pressure deficit between leaf and canopy air. The specific humidity and CO 2 concentration in canopy (C ac ) are approximated as q and C at the canopy top, respectively. Meteorological variables at the canopy top are estimated from those at height z * 1 using Monin-Obukhov similarity theory as a first-order approximation.
The maximum carboxylation rates per leaf area at 25 °C (μmol m −2 s −1 ) for sunlit leaves ( v c max 25,sun ) and shaded leaves ( v c max 25,shd ) are calculated based on de Pury and Farquhar (1997), but modified using the leaf parameters for sunlit and shaded leaves for the FvCBL model (e.g., Sun et al. 2014;Masutomi et al. 2016), instead of using the bulk parameters for sunlit and shaded leaves, as follows where L is the cumulative leaf area from the canopy top, χ n is a parameter for converting from leaf-nitrogen concentration to the maximum carboxylation rate, N 0 (mmol m −2 ) is the specific leaf-nitrogen concentration at the canopy top, N b (mmol m −2 ) is the specific leaf-nitrogen concentration not associated with carboxylation, k n is the coefficient of leafnitrogen allocation in the canopy, k b is the coefficient of beam-radiation extinction in the (5d) C Ec = f g sw , canopy, and V cmax25 is the maximum carboxylation rate at 25 °C for the whole canopy calculated as Note that the coefficient k n is defined differently using a factor of L c from the original model of de Pury and Farquhar (Hikosaka et al. 2016). The stomatal conductance model based on Leuning (1995) is given as where m l is an empirical coefficient, P n (μmol m −2 s −1 ) is the net leaf photosynthesis, C s (μmol mol −1 ) is the leaf surface CO 2 concentration, Γ (μmol mol −1 ) is the CO 2 compensation point of P n , D 0 is the reference vapour pressure deficit (kPa), and g sw,min (mol m −2 s −1 ) is the minimum value of g sw . We use values of m l , Γ, D 0 , and g sw,min following our earlier study (Ikawa et al. 2018). The FvCBL model also requires the aerodynamic resistance of a leaf, which is calculated based on the leaf heat-exchange coefficient and wind speed within the canopy (Masutomi et al. 2016).
Most of the FvCBL parameters are based on leaf gas-exchange measurements for a common rice cultivar (Koshihikari) in the rice Free-Air CO 2 Enrichment (FACE) experiment (Ikawa et al. 2018(Ikawa et al. , 2019 (Table 2). The temperature function is also revised for photosynthesis parameters based on the CO 2 concentration at the carboxylation site (Bernacchi et al. 2002). Accounting for mesophyll conductance may be preferable for accurately estimating photosynthesis under changing [CO 2 ] (Sun et al. 2014). For further details concerning the rice-crop energy-balance model, see Ikawa et al. (2018) and Maruyama and Kuwagata (2010).

Energy-Flux Models for Other Land-Use Categories
The bulk formulae based on Monin-Obukhov similarity theory for the surface stress τ (kg m −1 s −2 ) and surface fluxes of sensible heat and latent heat are used for the other land-use categories. Using the surface temperature T s , the surface stress and energy fluxes are calculated as

Table 2
Model settings for the rice-crop energy-balance model based on Ikawa et al. (2017Ikawa et al. ( , 2018 *Based on Tatsumi (Dyer and Hicks 1970;Kondo et al. 1978). For each land-use category, the aerodynamic roughness length z 0 (m), the thermal roughness length z T (m), the zero-plane displacement d (m), the evaporation efficiency β (dimensionless), the albedo (dimensionless), the soil heat capacity (J m −3 K −1 ), and the soil heat conductivity (W m −1 K −1 ) are set separately (Table 3). The energy fluxes (H and λE) are solved together with the heat-conduction equation of the soil. The temperature of the open water surface is left unchanged during the simulation.

Model Simulation
The model simulations are conducted for a specific study area (ca. 300 × 300 km 2 ) in Japan (34.5-37.5 °N, 137.5-140.5 °E) (Fig. 2). The simulation domain is gridded by a 3 × 3 km 2 grid, and each square is categorized into the aforementioned land-use types based on the land-use map of 2006 obtained from the Geospatial Information Authority of Japan (GSI). Rice paddies account for 10% of the total land area, and 64% of all the grid squares has a non-zero areal fraction of rice paddies. Altitude data are also obtained from the GSI. The weather and rice-growth conditions are taken as those on a typical summer day when rice fields are irrigated (water depth of 0.05 m) and covered fully with leaves (LAI was set at 4 m 2 m −2 ). The global shortwave flux R s (W m −2 ) is determined based on extra terrestrial solar radiation, the solar zenith angle, and precipitable water to account for the effect of humidity on the attenuation (Kuwagata et al. 1990). The initial state of the potential temperature θ is set at 25 °C at the sea level and increases vertically at the rate of 0.005 K m −1 . The wind speed is set to zero and the relative humidity is set to 70% at all vertical levels as the initial condition.
The top boundary conditions for wind speeds are set to zero. Radiation conditions are applied for the velocity components normal to the horizontal boundaries. Zero-gradient conditions are applied for scalars and velocity components parallel to the horizontal boundaries. Thus, the spatial variations of the velocity fields are induced by the thermally induced atmospheric circulation due to differences in energy balance and topography.
The simulations are run with a timestep of 30 s for the atmospheric model and 60 s for the rice energy-balance model. A spin-up simulation with a control condition is run for 20 h starting at 0700 local time (LT = UTC + 9 h) to generate typical summer meteorological conditions in the area. Then, the following experimental simulations are run for 16 h (from 0300 to 1900 LT): • Simulation 1 (Control) [CO 2 ] of 400 μmol mol −1 . Table 3 Model settings for the land-surface models for land-use categories other than rice paddies Land • Simulation 4 Increased rice-paddy area within each grid square by 20% from the control. • Simulation 5 Decreased rice-paddy area within each grid square by 20% from the control.
We assume that the time period over which the rice-paddy areas change by 40% corresponds reasonably to that in which [CO 2 ] increases from 280 to 800 μmol mol −1 on a century scale. The changes in rice-paddy area are compensated by the area of city within each grid square, while the areas of the other land-use types remain unchanged in each of the five simulations.
The downwards radiation fluxes at the surface (the global shortwave flux R s , the diffuse shortwave flux, and the longwave flux) are prescribed based on the initial state of the atmosphere; however, the present study do not consider the precise radiation transfer within the atmosphere and the feedback effects of the land surface on the radiation conditions (see Appendix 1). The photosynthetic photon flux density (μmol m −2 s −1 ) above the canopy is assumed to be proportional to R s (W m −2 ) by a factor of 2, and the distribution between the absorbed photosynthetic photon flux density, Q c,sun and Q c,shd is calculated following de Pury and Farquhar (1997). Atmospheric tendencies are not considered in this study.

Heat-Budget Analysis of the Atmospheric Boundary Layer
The heat budget of the atmosphere within the ABL is calculated (e.g., Kuwagata et al. 1990) from the beginning of the simulation (0300 LT) until 1300 LT when differences in simulated meteorological variables between different runs are pronounced. The aim of this analysis is to differentiate the effects of heating via the sensible heat flux at the Herein, t total (s) is the total simulation duration (10 h) in seconds and z * h is the peak ABL height during the period of the heat-budget analysis. We set z * h at 2 km based on the simulation without the topographic effects (Online Supplement Fig. S1).

Downscaling Meteorological Variables
The simulated meteorological variables for the lowest grid layer represent the average atmospheric state over the mixture of different land-use categories in each surface grid square, regardless of the fact that there might be different meteorological variables within each IBL (Fig. 1). First, we define the grid-mean meteorological variables (i.e., U, T, q, and C) downscaled for the height of z r within the IBL for each land-use category as The meteorological variables at z r (i.e., U a , T a , q a , and C a ) are estimated from those at z * 1 using Monin-Obukhov similarity theory (e.g., Garratt 1994). For example, the variables U a and T a are calculated as We select z r to be 1.5 m above the zero-plane displacement assuming that the meteorological variables at the height z r represent the state of the specific land use on the subgrid scale of a few hundred metres.

Evaluation of Model Performance
The reason for using the regional atmospheric model rather than the slab model, which was used in our previous study (Ikawa et al. 2018), is to simulate the diurnal cycle of air temperature as realistically as possible over the complex terrain in Japan. We use hourly observed data of air temperature from a rice paddy [(x, y) = (165 km, 185 km) in Fig. 2] and from a residential area ('city') [(x, y) = (167 km, 183 km)] within the simulation domain (Kumagaya city) using the dataset of Kuwagata et al. (2014), and we compare those with the simulated air temperature T a averaged over the areas within 20 km from each site and at altitudes below 100 m (the number of grids, n = 119 for rice and n = 122 for city) under different conditions of downwards shortwave radiation at the surface, R s . We choose the dataset of Kuwagata et al. (2014) because their study area comprises a typical landscape mixed with farm areas and adjacent cities in the Kanto plain, which accounts for a large fraction of the present study domain (ca. x > 150 km, y < 250 km in Fig. 2). Note that while Kuwagata et al. (2014) reported daily and monthly data, they did not explicitly present the hourly data used here.

Comparison with Observations
The diurnal patterns of air temperature T a in the control simulation and the simulation in which the R s values are reduced by 40% are similar to the observed air temperatures under similar meteorological conditions over both the rice-paddy and residential areas (Fig. 3). As the downwards radiation in the simulations is assigned typical summer values for this region, the diurnal variations in the value of R s do not agree perfectly between the observations and simulation results (Fig. 3b). Nevertheless, it can clearly be seen that the difference in air temperature between the two land-use types is greater under the higher R s values during the daytime (R s > 0) both in the observations and simulation results. The greater difference in T a values during the periods of higher R s values is consistent with observations within and near the domain of the present study (Kuwagata et al. 2014(Kuwagata et al. , 2018. As we do not consider the urban canopies in our land-surface model, the simulated diurnal patterns of air temperature in the city areas may not be accurate (Kusaka et al. 2001;Kusaka and Kimura 2004). However, while the effect of urban canopies is generally pronounced at night, the present study focuses primarily on daytime air temperature.
(14d) T z * 1 = z * 1 Π z * 1 − 273.15, Figure 4a shows the simulated potential temperature θ(z * 1 ) together with the relative sizes of the velocity vectors for the control simulations. The simulations are able to generate typical velocity fields that are strongly influenced by the sea breeze around the coastal areas (e.g., Kuwagata et al. 1994). It is evident that the value of θ(z * 1 ) is greater in the area of Tokyo bay [(x, y) ≈ (200 km, 150 km)] than in the adjacent areas northward where rice paddies are extended over large areas, even though these areas have similar topographies (Fig. 2a). Doubling [CO 2 ] clearly increases the value of θ(z * 1 ) in an extensive area of the simulation domain (Fig. 4b). At 1300 LT when the difference in potential temperature θ(z * 1 ) among the simulations is greatest, the difference in θ(z * 1 ) is almost 0.7 °C in some areas. The temperature difference is particularly clear wherever rice paddies are aggregated, as shown in Fig. 2b.   Fig. 3 Performance of model evaluated by comparing observed air temperature in rice-paddy and city areas in July-August 2012 (61 days) (Kuwagata et al. 2014, 'Obs.' in the legend) and simulated air temperature (T a is the air temperature at 1.5 m above the zero-plane displacement) averaged over rice paddies (n = 119) and cities (n = 122) within 20 km of each site ('Model'): a comparisons by categorizing observed data into two sets with the daily maximum value of the global shortwave flux either higher (34 days, Sunny) or lower (27 days, Cloudy) than 800 W m −2 ; b control simulation and simulation with 40% of shortwave flux performed for sunny and cloudy days Both changing the areal fraction of rice paddies and changing [CO 2 ] clearly affects the simulated value of T a over rice paddies and cities (Fig. 5). We reason that the change in the value of T a in the simulations for the other land-use types reflects the changes in air temperature over a large spatial scale beyond the rice-paddy areas. We therefore choose to investigate the air temperature T a in cities and compare it with the case of rice paddies.

Near-Surface Air Temperature
For rice paddies, the changing [CO 2 ] has a greater effect on T a than did the changing land use. Doubling [CO 2 ] increases the peak value of T a by 0.44 °C ± 0.09 °C, and the increased value of T a exceeds 0.7 °C in some areas (n= 4231) (Figs. 5a and 6a). In contrast, θ (z *1 )(Sim2)−θ (z *1 )(CTRL) (K) θ(z *1 ) (K)  decreasing rice-paddy areas by 20% increases the peak value of T a by 0.11 °C ± 0.11 °C (n= 4231) (Figs. 5a and 6c). For cities, by contrast, the effect of land-use change on the air temperature T a is greater than the effect of changing [CO 2 ] (Fig. 5b). Doubling [CO 2 ] and reducing rice-paddy areas increases the peak value of T a over the cities by 0.07 ± 0.07 °C and 0.12 ± 0.13 °C, respectively (n= 4231) (Figs. 5b and 6b, d). The effects of land-use change on the value of T a are more spatially variable compared with those of the change in [CO 2 ] (Fig. 6), and some cities experience an increase in the value of T a of approximately 0.7 °C in response to decreasing rice-paddy areas (Fig. 6d). If both doubling [CO 2 ] and reducing rice-paddy areas were included, the peak value of T a of the rice paddies and cities would increase by 0.53 ± 0.15 °C and 0.18 ± 0.16 °C, respectively (data not shown).

Fig. 6
Histograms of the air-temperature change 1.5 m above the zero-plane displacement (dT a ) in simulation 2 (increased background atmospheric CO 2 concentration by 400 μmol mol −1 ) compared with the control conditions (CTRL) (a, b) and in simulation 5 (decreased rice-paddy area by 20%) compared with the control conditions (c, d) for rice paddies (a, c) and cities (b, d) at 1300 LT. The data are averaged from 4231 grid squares in which the areal fractions of city and rice paddy were non-zero. The numbers indicate the mean and standard deviation

Other Meteorological Variables and Energy Balance
The directions of changes in the specific humidity, q a in the different simulations are opposite to those of the air temperature T a (Fig. 7a). Doubling [CO 2 ] decreases the magnitude of q a in the rice paddies by 0.40 g kg −1 at 1300 LT. Decreasing the rice-paddy area decreases the value of q a , and increasing the rice-paddy area increases the value of q a . Doubling [CO 2 ] decreases ΔC a (the C a anomaly from [CO 2 ]) (Fig. 7b). The differences in ΔC a between the simulations are attributed to the differences in the canopy photosynthesis (Online Supplement Fig. S2). Similar to the air temperature T a in the rice paddies (Fig. 5a), the changes in q a and ΔC a in the rice paddies are greater in response to [CO 2 ] changes than to land-use changes (Fig. 7a, b). Increasing the rice-paddy areas increases the wind speed U a and decreasing the rice-paddy areas decreases the value of U a because of the smaller roughness length of rice paddies compared with cities (Fig. 7c), but these changes are marginal. We held R s invariant between different simulations (Fig. 7d). The differences in the values of T a and q a between the rice paddies and the cities in the simulations are attributed to the differences in the energy distribution of the fluxes λE and H. The grid-mean energy-balance components differ more in the simulations with different rice-paddy areas than in those with different [CO 2 ], whereas the energybalance components in the rice paddies are influenced more by different [CO 2 ] than by different areal fractions of rice paddies (Fig. 8). A clear change in the value of H of Fig. 7 Simulated meteorological variables over rice paddies (1.5 m above the zero-plane displacement): a specific humidity (q a ), b anomaly of atmospheric CO 2 concentration within the IBL from the background concentration (ΔC a ), c wind speed (U a ), and d the global shortwave flux (R s ). The simulations were performed under control conditions (CTRL), increased or decreased background atmospheric CO 2 concentration (+ 400 or − 120 μmol mol −1 ), and increased or decreased rice-paddy areas (+ 20% or − 20%). The data are averaged from 4231 grids in which the areal fraction of rice paddy is non-zero the rice paddies is observed in the simulation with elevated [CO 2 ], with the value of H increasing by up to 50 W m −2 (Fig. 8e). The changes in λE of the rice paddies between different simulations are attributed to different values of g sw , but the differences in , increased or decreased background atmospheric CO 2 concentration (+ 400 or − 120 μmol mol −1 ), and increased or decreased rice-paddy areas (+ 20% or − 20%). The data are averaged from 4231 grids in which the areal fraction of rice paddy is non-zero Fig. 9 a Stomatal conductance (g sw ) of rice leaves and b canopy temperature (T c ) obtained from the model simulations. The simulations were performed under control conditions (CTRL), increased or decreased background atmospheric CO 2 concentration (+ 400 or − 120 μmol mol −1 ), and increased or decreased ricepaddy areas (+ 20% or − 20%). The data are averaged from 4231 grids in which the areal fraction of rice paddy is non-zero λE between the simulations are much smaller than that for g sw on a percentage basis (Figs. 8d and 9a). Differences in the effective aerodynamic temperature, T c between different [CO 2 ] simulations are greater than the differences in the air temperature T a (Fig. 9b). On average, the elevated [CO 2 ] increases T c by 1.0 °C at 1300 LT. Slight differences in g sw between different rice-paddy areas and consequent small changes in the value of T c and energy fluxes are likely due to different meteorological conditions (i.e., T a and q a ), as shown in Figs. 5 and 7.

Temperature, Humidity, and CO 2 within the Atmospheric Boundary Layer
The differences in potential temperature θ between the simulations are evident up to 1 km above the land surface during the daytime (Fig. 10a). Similarly, the differences in the specific humidity q are clear up to 1 km above the ground (Fig. 10b). The vertical profiles of both θ and q deviate more from the control in the cases of different rice-paddy areas than in the cases of different [CO 2 ]. The clear difference of the vertical distribution of CO 2 concentration C among the simulations extends to even the upper regions of the atmosphere (≈ 2 km) (Fig. 10c). Note that the formation of a mixed layer in the area is influenced strongly by local advection, as supported by observations (Kuwagata et al. 1990;Kimura et al. 1997). Simulations without topography result in a typical vertical θ structure in a mixed layer (e.g., Stull 1988) (Online Supplement Fig. S1).

Heat-Budget Analysis of the Atmospheric Boundary Layer
The net heating of the atmospheric column (Q n ) is low near the coastal plains, moderate in inland areas with extensive rice paddies, and particularly high at the bottom of the basin Fig. 10 Vertical distributions of simulated a potential air temperature (θ), b specific humidity (q), and c deviation of atmospheric CO 2 concentration from the background (ΔC) over the terrain normalized height ( z * ) at 1300 LT averaged for rice-paddy areas (areal fraction > 0.5, n= 301) under control conditions (CTRL), increased or decreased background atmospheric CO 2 concentration (+ 400 or − 120 μmol mol −1 ), and increased or decreased rice-paddy areas (+ 20% or − 20%) (Fig. 11a). These results are consistent with the observations reported in Kuwagata et al. (1990). These spatial patterns of Q n are clearly governed by heating via advection, Q adv (Fig. 11e).
Doubling [CO 2 ] increases Q n in most areas, particularly in the areas with extensive rice paddies because of the increased heating due to the surface sensible heat, Q H (Fig. 11b, d). Conversely, the magnitude of Q adv in the areas with rice paddies tends to decrease, partially offsetting the effect of Q H (Fig. 11f).
(e) (f) Fig. 11 Heat-budget analysis averaged over 10 h (0300-1300 LT). The net heating in the atmospheric column averaged over the 10 h (Q n ) is separated into the contributions from the sensible heat flux from the land surface (Q H ) and advection (Q adv ) under control conditions (CTRL) (a, c, e). Their differences in simulations with increased background atmospheric CO 2 concentration (+ 400 μmol mol −1 ) compared with control conditions are also shown (b, d, f). The contours indicate topography

Reduced Heat-Mitigation Effect of Rice Paddies under Future CO 2 Concentration
The simulations show that doubling [CO 2 ] clearly increases the air temperature T a in both rice paddies and areas of other land-use types (Figs. 5 and 6). The clear difference in the value of T a between the different [CO 2 ] simulations near midday is attributed to the changes in the energy balance via the CO 2 response of the stomatal conductance, g sw (Figs. 8 and  9). Our model simulations do not consider other factors that are sensitive to CO 2 , such as

Sim2 -CTRL (without topography)
Sim2 -CTRL Fig. 12 The relationships between the differences in the components of the atmospheric heat budget (Q n , Q H , and Q adv ) and the differences in air temperature 1.5 m above the zero-plane displacement (dT a ) in the cities and rice paddies between the simulations of the control (CTRL), doubling background atmospheric CO 2 concentration (Sim2), and reducing the rice-paddy area by 20% (Sim5). Data are selected from the grid areas with a small rice-paddy areal fraction (< 1%, n = 636) and a large rice-paddy areal fraction (> 60%, n = 104) for panels a-f. For the comparisons between the control (CTRL) and Sim5 (g-i) simulations, data are selected from the grid areas with a small rice-paddy areal fraction (< 1%, n = 760) and a large rice-paddy areal fraction (> 40%, n = 216) based on the land-use map with reduced rice-paddy areal fraction the dynamic vegetation change in other land-use types (e.g., Notaro et al. 2007). Nonetheless, our results suggest that elevated [CO 2 ] increases air temperature in rice paddies at a greater rate than in other land uses that would be less sensitive to the [CO 2 ] change, such as city areas, provided that the atmospheric circulation pattern remains similar under elevated [CO 2 ]. The spatial variations of the T a change between the different simulations shown in Fig. 6 are influenced by the balance between the local heating and advection of the immediate atmospheric column (Fig. 12). During the daytime, the atmosphere above the areas dominated by rice paddies receives heat from the surrounding areas, resulting in a positive value of Q adv (Fig. 11e), which is positively related to the areal coverage of rice paddy within each grid (Online Supplement Fig. S3). Under the elevated [CO 2 ] condition, the heat gradient between rice-paddy dominated areas and the surrounding areas decreases, resulting in the negative change in the magnitude of Q adv (Figs. 11f and 12c, f). Despite reduced value of Q adv , the atmosphere above the rice-paddy dominated areas gains heat due to the intensified local heating (Q H ) (Figs. 11d, and 12b, e), resulting in a positive change in T a under elevated [CO 2 ] both in rice paddies and cities.
Most of the rice paddies and cities within the simulation domain experience a positive change in T a under elevated [CO 2 ] (Fig. 6a, b), although there are many grid areas with a very small rice-paddy areal fraction of less than 1%. While the increase in T a of both rice paddies and cities in the rice-paddy dominated areas is attributed to the local heating (i.e., Q H ), the positive change in T a in the areas with a small rice-paddy areal fraction is mainly attributed to the increase in Q adv or reduced negative Q adv . These relationships between the changes of T a and heat-budget components are also observed in the simulation without topography (Fig. 12d-f).
The difference in T a between the simulations with different rice-paddy areas is also influenced by the changes in heat-budget components (Fig. 12g-i). Reducing rice-paddy areas increases the value of Q H and decreases the value of Q adv in the rice-paddy dominated areas, and the net change in the atmospheric heat balance of Q n is positively related to the change in T a . The systematic bias in the T a change between rice paddies and cities that is observed in the case of different [CO 2 ]. (Fig. 12a-f) is not observed in the case of different areal fractions of rice paddies (Fig. 12g-i) because the changes in energy balance of the land surface and resultant effects on the IBL structure are small in the case of different rice-paddy areas (Fig. 8). It is noteworthy that the changes in T a due to the change of rice-paddy areas vary considerably by location (Fig. 6c, d). The large spatial variation of the change in T a associated with the land-use change is likely related to the broad range of the change in the magnitude of Q adv (Fig. 12i).
Within similar geographical conditions (e.g., the areal fraction of rice-paddy greater than 0.5 and altitude less than 100 m), the spatial variations of the T a change in rice paddies also depend on the wind speed. The change in T a under elevated [CO 2 ] is clear in those grid squares with a moderate wind speed U a (Fig. 13), which is necessary to enhance atmosphere-land surface mixing (e.g., Ikawa et al. 2018), but a further increase in the wind speed enhances the local advection and dilutes the signal from the immediate land surface. This study assumes that land uses other than rice paddies are carbon neutral because our purpose is to isolate the CO 2 effects of rice paddies on the atmosphere. This assumption enables us to use the value of ΔC a as a tracer that tracks the influence of rice paddies on the atmosphere: the greater ΔC a , the greater the impact of rice paddies on the atmosphere.
Decreasing rice-paddy areas has a greater effect on air temperature over mixed land surfaces than do [CO 2 ] effects as indicated by T a in cities and also by the vertical distributions of potential temperature θ within the ABL (Figs. 6 and 10). Nevertheless, the change in T a in the cities (Fig. 6), the surface energy fluxes (Fig. 8), and the vertical distribution of θ ( Fig. 10) among different [CO 2 ] conditions indicate that the effect of doubling [CO 2 ] still affects the air temperature and energy balance over the areas beyond the rice paddies. The [CO 2 ] effect on θ and T is approximately half that of reducing the rice-paddy areas by 20%, which is shown clearly in the grid-mean Bowen ratio (Fig. 8c). Additional model simulations indicate that the average magnitude of the changes in T a in cities caused by doubling [CO 2 ] is equivalent to a 10% decrease in the rice-paddy area within each grid square. We therefore argue that the impact of the physiological response to [CO 2 ] on the regional temperature is not negligible, at least in comparison to land-use change.

Atmospheric Effects on Rice-Plant Temperature Mediated via Stomata
Studies based on the FACE experiment have shown that the plant temperature increases via stomatal closure under elevated [CO 2 ], and that this increased temperature affects the rice grain quality, phenology (e.g., days to heading), and physiological processes (e.g., flowering time) Zhang et al. 2015;Kobayasi et al. 2019). Therefore, the question arises as to whether we need to consider an extra increase in plant temperature through the increase in atmospheric temperature via stomatal closure under elevated [CO 2 ].
We computed the canopy temperature, T c under elevated [CO 2 ] while keeping the other meteorological variables (i.e., T a , q a , and ΔC a ) of the control simulation constant (Fig. 14a). This simulation yields the value of T c under elevated [CO 2 ] without considering the feedback from rice paddies to the atmosphere. The simulation results indicate that the feedback from rice paddies under doubled [CO 2 ] affects the value of T c by 0.19 °C at 1300 LT (Fig. 14b). The weaker sensitivity of T c compared with T a is attributed to the fact that an increased temperature tends to enhance the latent heat flux compared with the sensible heat flux, and the increase of T a (+ 0.44 °C) increases the value of T c by 0.28 °C if no humidity and CO 2 feedbacks are considered (Online Supplement Fig. S4). Further consideration of the effects of decreased humidity and decreased the CO 2 concentration C due to enhanced photosynthesis result in a net increase in the value of T c by 0.19 °C.
ΔC a (CTRL) [μmol mol -1 ] U a Fig. 13 Wind speed (U a ) and difference of air temperature 1.5-m above the zero-plane displacement (T a ) over rice paddies in simulation 2 (background atmospheric CO 2 concentration increased by 400 μmol mol −1 ) compared with the control (CTRL) simulation at 1300 LT. The anomaly of atmospheric CO 2 concentration from the background (ΔC a ) is also shown in colour. The data are only for altitudes below 100 m and the areal fractions of rice paddy greater than 0.5 (n = 151) The significance of the changes in T c due to the atmospheric change via stomatal closure may deserve further investigation. Rice grains are very sensitive to temperature, especially at the reproductive stage (Usui et al. 2016;Matsui and Hasegawa 2019). Moreover, temperature has a cumulative effect over the course of the rice-growing stage (Zhang and Tao 2013;Ishigooka et al. 2017), and even a small change in temperature may have a consequence in the developmental stage over time. Additionally, note that crops generally have periods in which they are particularly sensitive to temperature (e.g., germination, flowering, and ripening), and whether we should consider vegetation feedback to the atmosphere depends on the sensitivity of a particular phenological or physiological process to temperature (Combe et al. 2015).

Limitations and Future Work
This study is limited to typical summer conditions in Japan, when the effects of plant physiological processes on the atmosphere are the greatest. However, the ecophysiological properties of rice paddies change between the growing and non-growing seasons and also among different growth stages (Ono et al. 2013(Ono et al. , 2015. Although our aim in using a regional atmospheric model is to simulate the diurnal variations of meteorological variables as realistically as possible, we do not consider the air-land surface CO 2 exchange in areas of land uses other than in rice paddies. This is because we are interested in isolating the effects of the CO 2 responses of rice plants, and we assume that considering momentum, heat, and vapour fluxes would be sufficient to create the baseline structure of the ABL (Santanello et al. 2013). The complexity of CO 2 fluxes in areas of other land-use types such as cities (Ueyama and Ando 2016) is also why we do not consider CO 2 -flux effects from other land-use types. Different simulations with topography conducted here do not result in a clear difference in velocity fields. However, changing the energy balance over extensive land-surface areas beyond rice-paddy areas may alter the patterns of atmospheric circulation. The changes in the atmospheric variables via different atmospheric flows under elevated [CO 2 ] is worth further investigation.

(b) (a)
Fig. 14 a Simulated canopy temperature (T c ) using average meteorological variables (air temperature, specific humidity, and anomaly of atmospheric CO 2 concentration from the background) simulated under control conditions (CTRL), doubled background atmospheric CO 2 concentration (simulation 2, 'Sim2'), and doubled background atmospheric CO 2 concentration with other meteorological variables of CTRL ('Sim2 without feedback'). b Difference in T c between the simulations The choice of the stomatal conductance model may be another factor in the uncertainty of the simulation results. The sensitivity of the stomatal conductance to CO 2 is reduced by introducing the approximations of Medlyn et al. (2011). For example, the peak dT c value decreases by 0.3 °C in Fig. 14b and the increase in atmospheric heating would be less. The difference is primarily attributed to the fact that Leuning-type models include the intercept Γ for the CO 2 effect, whereas the Medlyn model does not, at least in its original form. Further assessments of the stomatal conductance model are beyond the scope of this study. Because the increase in T c shown in Fig. 14 is reasonable compared with the results of the past rice FACE experiments (Yoshimoto et al. 2005;Sikma et al. 2020), we retain the current model for this study.
Regarding further realistic simulations, some other key aspects that can be considered to improve our model are an accurate temperature response and cloud formation. The effect of doubling [CO 2 ] on T a is conserved in both the rice paddies and cities under the simulations with a background air-temperature increase of 2 °C (average T a increases by 0.44 °C and 0.07 °C, respectively, at 1300 LT). However, incorporating appropriate temperature functions in the model depends on further investigation on the temperature responses of ecophysiological properties of rice paddies. Cloud formation interplays with [CO 2 ] and affects the radiation budget. Decreased humidity under elevated [CO 2 ] may suppress ABL cloud formation and affect the radiation intensity (Vilà-Guerau de Arellano et al. 2012).
The good agreement between the simulated and observed diurnal patterns in air temperature in rice paddies and cities suggests that the coupled atmosphere and land-surface model presented herein may also be useful for investigating how land-use changes affect the regional climate (Fig. 3). Our results suggest that the effects of land-use change on the atmosphere vary considerably by location, and the spatial variation may be associated with the variation of the heat advection. Land-use changes may be assessed further based on real land-use change scenarios using land-use maps from different years and by effectively introducing the key land-surface processes of other types of land use (e.g., urban canopy, CO 2 sequestration by forests) (Kusaka et al. 2001;Kusaka and Kimura 2004;Watanabe et al. 2004), depending on the research hypothesis to be tested.

Conclusion
A coupled atmosphere and crop energy-balance model is developed, and meteorological variables and energy balance near the land surface are simulated for past, current, and future [CO 2 ] and compared with cases involving different areal fractions of rice paddies. The model simulations show that doubling [CO 2 ] clearly increases the daily maximum air temperature in most of the areas when the air temperature is highest on a typical sunny day. Our model simulation also quantitatively evaluates the balance between the local heating and advection. The areas with a greater areal extent of rice paddies experience a net heating under elevated [CO 2 ], despite the negative change in advective heating received from the surrounding areas. Although the land-use change postulated herein has a greater impact than [CO 2 ] effects on the air temperature over larger spatial scales, the air temperature within the IBL (viz. T a ) over rice paddies is increased more by changing [CO 2 ] than by changing the rice-paddy areas. These results suggest an important interplay of plant physiological responses to CO 2 on the heat-mitigation effects of rice paddies under a changing climate.
where ρ 0s is the air density at the sea level (kg m −3 ) and R is the gas constant of air (J kg −1 K −1 ).
The downwards longwave flux R ld (W m −2 ) is fixed to 400 W m −2 because its diurnal variation is much smaller than that of R s . In the simulation with an increased air temperature, the value of R ld under the increased air temperature R ld+ , is approximated as where σ is the Stefan-Boltzmann constant, and h is the average atmospheric height over the calculation domain. The fraction of diffuse radiation f d used in the photosynthesis model is estimated based on Spitters et al. (1986) as follows