Central Taiwan’s hydroclimate in response to land use/cover change

Land use/cover change (LUCC) has taken place since the 1990s in central Taiwan; however, its impacts on the local and regional hydroclimatology are not understood thoroughly. This study is grounded in a numerical experiment using the Weather Research and Forecasting (WRF) model and statistical assessments of continuous land cover and gridded precipitation data derived for central Taiwan. We incorporate survey-based land use data in 1995 and 2007 in driving WRF to simulate selective non-rainy and rainy (dry and wet) cases under weak synoptic forcings in July and August (JA). The two land-use conditions reveal changes in simulation fields on account of increased urban and built-up lands. Results averaged over the dry cases show increased (diminished) sensible heat fluxes and 2 m temperatures (latent heat fluxes and 2 m specific humidity) in 2007 compared to that in 1995. The wet-case simulation further identifies intensified precipitation over the downwind areas of urban and built-up lands, strongly subject to local topography and prevailing winds. Statistical assessments of the Landsat land cover and gridded precipitation data verify significant increasing trends in urbanization and the JA rainfall. Regression-based analysis that scales the effect of the LUCC on the change in precipitation corroborates the WRF simulation: LUCC has induced eastward, downwind association with the JA rainfall.


Introduction
Land use/cover change (LUCC) can induce significant variations in land surface fluxes and water budget partitioning, thereby impacting local and regional climate [1][2][3][4][5][6][7]. The influences of LUCC in the Earth system have been a major concern for the global change assessment. For example, extensive anthropogenic irrigation and deforestations result in changes in canopy resistance, surface albedo, soil moisture, and soil temperature that can potentially modify surface heat fluxes [8,9]. Many studies have shown irrigation's impacts on local, regional, and global climate, and even the atmospheric circulation changes (e.g. [10][11][12][13][14][15]). While deforestation usually causes less evapotranspiration and higher temperature, leading to reduced precipitation [16][17][18], precipitation in deforested areas could also be more extensive than that in high-density forests [19]. Reference [20] further indicated that deforestation could increase local convection and precipitation by introducing more lateral water vapor transport from the surrounding oceans over the Maritime Continent regions.
Variations of heavy rainfall events subject to regional LUCC have been discussed in previous studies (e.g. [1,21,22]). Some other studies showed that the rainfall cycles through the controls (e.g. transpiration and the partitioning of surface net radiation) of different land surface conditions can further alter the feedback mechanisms between the land surface and consequential strong convective events [23][24][25][26]. While copious studies have investigated how the dynamics of urban expansion (as one type of LUCC) impacts local atmospheric circulations, precipitation quantities, and microenvironment (temperature or humidity) through analyzing observed data or modeling results Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Over the past few decades, the world has been experiencing significant LUCC, and Taiwan is no exception. Several major administrative areas in Taiwan have undergone rapid development, and central Taiwan ranks among the top. To simulate the impact of urbanization in central Taiwan on precipitation, [32] conducted a numerical experiment using the MM5 model. Their findings based on a hypothetical urban expansion showed that the warmer and drier conditions of built-up lands create a more unstable atmospheric condition that further enhances local convection. In association with the effect of the landsea breeze and topography, their simulations suggested increased precipitation over downwind areas. [30] thereafter used the Weather Research and Forecasting (WRF) model to examine how meteorological simulation responds to different land use data in Taiwan, including the default USGS and MODIS and a proprietary product derived from SPOT satellite images. They found the respective overestimation of irrigated cropland and urban areas in the USGS and MODIS land use data, resulting in erroneous simulations (e.g. overpredicting the surface wind speed and daytime temperatures). Similar studies are abundant, and most of them imply the importance of accounting for more realistic and accurate LUCC conditions to avoid significant biases in numerical simulations.
To the best of our knowledge, there exists a gap in literature emphasizing the assessment of LUCC on island-scale climate subject to complex terrain. This research gap could be attributed to the lack of extended observations for both climate and land use/cover data at a resolution fine enough to clarify such intricate land-atmosphere relationships.

Data and methods
2.1. Land use/cover data: sources and derivation To understand the impact of LUCC on central Taiwan (i.e. Miaoli, Taichung, Changhua, Yunlin, and Nantou counties/cities in this study) climate, we acquired and derived two types of land use/cover data from different sources. The first is survey-based data acquired from the National Land Surveying and Mapping Center (NLSC, www.nlsc.gov.tw/LUI), who applied cadastral survey supplemented by aerial images, GIS analysis, and field survey to the development of the land use data at a spatial resolution of 50 m in 1995 and 2007 (NLSC95 and NLSC07 hereafter). The NLSC adopts a three-tier hierarchical classification system that shows different number of classes in the two years. To prepare a compatible land cover data in the numerical model (i.e. WRF) environment, we performed reclassification for the NLSC data in both years according to the 24 land-use categories retrieved from the default USGS data. A series of GIS practices (e.g. reproject and resample) were used to prepare the 1 km NLSC data for Taiwan with 10 common USGS land-use classes in 1995 and 2007. We plotted the processed NLSC95 and NLSC07 data (figure 1) and recognised the predominant LUCC, which is from agricultural land use types to urban and built-up lands. The processed NLSC data were regarded as the true land use/cover states used for driving WRF; the same data were also used for calibrating satellite-based land cover from 1995 to 2007 described next.
Owing to the two-year NLSC data that fell short of providing a continuous correspondence between LUCC and local circulations, we derived the second type of land cover data from Landsat 4-5 Thematic Mappter (TM) satellite images (available from the EarthExplorer, USGS, https://earthexplorer.usgs. gov/) at a 30 m resolution according to the following procedure. We scrutinised all Landsat images in wintertime from 1995 to 2007, identifying one image per year with the least cloud cover over central Taiwan. Data collection category level is L1TP, indicating the data have been radiometrically calibrated and orthorectified using digital elevation models. Afterward, we applied a hybrid unsupervised-supervised classification scheme based on ISODATA (Iterative Self-Organizing Data Analysis Techniques Algorithm) and MLE (Maximum Likelihood Estimate) to classify satellite imagery in ERDAS IMAGINE. Similar concepts of using a hybrid classification scheme for the investigation of LUCC have been employed by other studies (e.g. [33]). In our case, ISODATA first recognised various clusters from the original image in any given year, which was then merged into five dominant classes in the region (i.e. urban, cropland, mixed forest, water body, and barren) with the guidance of the unchanged areas between the NLSC95 and NLSC07 data. Signature files generated from ISODATA were used by MLE to map the five land cover classes in each year, and each land cover map was then superimposed with water bodies in 1995 and unchanged areas between the NLSC95 and NLSC07 data, assuming these areas remain unchanged throughout the years. The final land cover data were upscaled to 5 km spacing, and the percentage of each land cover class in a coarser grid can thus be derived. The classification results perform to a satisfactory level, with high spatial correlation coefficients and accuracy (both>0.8) for urban and overall classes (see figures S1 and S2 is available online at stacks.iop.org/ERL/15/034015/mmedia).

Precipitation data and case selection
One of the essential elements of the case selection for numerical simulations and the characterization of spatiotemporal patterns is precipitation data. We used both station and gridded precipitation data in this study: while the former was acquired from 216 weather stations over central Taiwan managed by the Central Weather Bureau (avaliable from the Data Bank of Atmosphere & Hydrologic Research, https://dbar. pccu.edu.tw/Default.aspx), the latter product was recently developed by the Taiwan Climate Change projection and Information Platform (TCCIP, https://tccip.ncdr.nat.gov.tw/). It is worth noting that the TCCIP data is the best available gridded data for the whole Taiwan island with very fine spatiotemporal resolution (5 km, daily) and extended records , developed using all available station data and other gridded products (e.g. APHRODITE, Yatagai et al [34]). Rigorous QA/QC, filling of missing data, and the temperature-rainfall-topography relationship have all been accounted for in the generation of the TCCIP gridded data [35].
The temporal resolution of the station data is hourly, which is finer than the daily gridded data. We wished to identify cases that exhibit more intense convective rainfall in the afternoon (1200-1800 local standard time, LST), so hourly rainfall records available from the station data were used for case selection for better identifying rainfall events in the designated time frame. Non-rainy (dry) and rainy (wet) events in July and August (JA) were identified from the station data from 1999 to 2016, and events affected by typhoons were excluded from our selection based on historical typhoon warnings. On the other hand, we derived the total cumulative rainfall and maximum daily rainfall in JA from the gridded data for further statistical assessments.
Case selection was conducted to illustrate the differences between energy budget driven by the NLSC95 and that by NLSC07 land use, in compliance with a specific procedure summarized herein. The specific procedure was designed to identify the most representative 10 dry and 10 wet events in JA by systematically scrutinising all the events from 1999 to 2016. The 10 dry cases were selected based on (1) 24 h accumulated rainfall less than 5 mm for all the stations lasting at least two consecutive days, and (2)<10 mm rainfall occurring at two stations at most for the third consecutive day (before or after). In contrast, the 10 wet cases were selected from those rainy events in which at least 40 stations were exhibiting maximum hourly rainfall greater than 15 mm during 1200-1800 LST. Table S1 lists the dates of all the dry and wet cases, and figure S3 shows the observed 6 h rainfall for the wet cases using station data. The aforementioned systematic search was implemented by means of two criteria based on the number of stations and rainfall thresholds. In fact, sensitivity analysis of the criteria has been performed: decreasing (increasing) the number of stations and/or the rainfall threshold would result in an increase (decrease) in the number of cases; for example, if decreasing 40 stations to 20, or 15 mm to 10, we would obtain approximately triple or double amount of wet cases. Along with examining the While dry-case simulation is primarily to explore changes in energy budget in generally sunny conditions, wet-case simulation can reveal how LUCC influences thunderstorm characteristics (e.g. intensity and spatial patterns). According to the case selection procedure, major weather systems (e.g. typhoons and weather fronts) and their impacts should be largely circumvented, enhancing the sensitivity of the WRF simulation to specified LUCC.

Model configuration and simulations
We used the WRF version 3.9.1 [36] with the Noah land surface model [37] to conduct numerical simulations. The simulation domain is nested and shown in figure 1. The outer coarsest domain has a resolution of 18 km (D1), ranging down to 6 and 2 km (D2 and D3). The vertical layer is composed of 45 sigma levels, and vertically stretched with finer resolution (from 50 to 350 m) below 2 km and coarser resolution (about 625 m) above 2 km. Regarding the cumulus parameterization scheme, the most common one applied to East Asian (Taiwan) regions is either the Kain-Fritsch or the Betts-Miller-Janjic (BMJ) scheme [38,39]. A preliminary assessment of the two schemes revealed that the BMJ scheme performed better in estimating precipitation amounts and patterns in our cases, so the BMJ scheme [40][41][42][43] was adopted and used in D1 and D2. No cumulus parameterizations were used in D3. Other physical options set to initialize the WRF simulation include the RRTM [44], Goddard [45], WSM5 [46], YSU [47], and the revised MM5 Monin-Obukhov [48] schemes for the longwave radiation, shortwave radiation, cloud microphysics, PBL, and surface layer schemes, respectively. The selection of these physical parameterizations (listed in Tale S2) was based on recommendations from previous studies [38,49,50].
We acquired the initial and boundary conditions from the National Centers for Environmental Prediction's (NCEP) Final Analysis (FNL) of global data at 1°×1°spatial and 6 h temporal resolutions. The total simulation time of each case was 36 h, beginning at 1200 UTC the day before and ending at 0000 UTC the day after. According to previous numerical studies (e.g. [49,51,52]), it is known that the time required by the regional model (e.g. WRF) to develop realistic cloud and rainfall patterns from the large-scale initial condition is generally 6-12 h, and the model shows the highest skill over 12-36 h. Thus we used the first 12 simulation hours as the model spin-up period, and then obtained major outputs for analysis from the remaining 24 h. The designed numerical experiment intends to capture the fine-scale meteorological process (e.g. precipitation) and to partition the energy budget for the aforementioned cases under weak synoptic forcing.
In order to discriminate central Taiwan's climate in response to two land use data in 1995 and 2007, we conducted two sets of identical WRF runs except that one was driven by the NLSC95 data and the other by the NLSC07 data. Areas outside central Taiwan were configured with the NLSC07 data for both sets of runs. In other words, we performed 2×10 runs for both dry and wet cases, and then obtained two sets of dryand wet-case results for further analysis. In both sets of simulations, we assured the consistency among the physical options (i.e. two land-use conditions, but consistent physical parameterizations).

Statistical analysis of LUCC-precipitation relationship
To corroborate findings from numerical simulations, we conducted statistical assessments of observed datasets. First, trend analysis was applied to the Landsat land cover as well as the gridded precipitation data using Sen's slope estimator [53]. Significant trends estimated by Sen's slope, while irrespective of the cause, can help identify sensitive areas bearing possible resemblance to other observed data and simulated patterns (e.g. intensified or weakened precipitation at specific locations resulting from the WRF simulation).
Another regression-based approach was adopted to scale the effect of LUCC on the interannual variations of precipitation, Assuming the change of precipitation (DP) is is linearly dependent on the change of urban/built-up percentage (DX ) in a particular grid, which can be expressed as where dP dX can be approximated as the coefficient of simple regression of P against X. Similar discussions were made by [54] and [55], while in their case, different correspondences were made (e.g. between hydrological quantities and temperature or precipitation). In our case, the above analysis was performed both on-grid and off-grid: whereas the former is to identify any significant correlations between P and X at the same grid, the latter is to seek any remote connections with certain displacement of the P grid. Technically, we followed the equation below for the off-grid analysis: where i and j denote the grid location, m and n the total number of zonal and meridional grids, E and SN the eastward and southward/northward displacement of the P grid. If E and SN are both zeros, the equation is identical the on-grid analysis. In addition, urban/ built-up percentage X was derived for 5 km grids in consistent with the alignment of P grids for this analysis, and regression is not performed for precipitation data outside central Taiwan or over the ocean (i.e. no westward displacement).

WRF simulation in response to LUCC
To better illustrate the differences in energy budget and rainfall patterns caused by the two sets of land use data, we calculated the average of simulated fields (e.g. temperature, fluxes, and precipitation) over the 10 dry (or wet) cases driven by the NLSC95 and NLSC07 land use, and then subtracted the averaged 1995 fields from the 2007 fields to generate the composite plots. In the dry-case simulation, we found that increased urban and built-up lands enhance sensible heat fluxes and diminish latent heat fluxes, as shown in figure 2. Changes in energy budget further lead to an increase in 2 m temperatures and reduction in 2 m specific humidity (q), illustrating an urban heat island in general. In the wet-case simulation, similar but less prominent patterns of energy fluxes and temperatures were obtained. The increase (reduction) in sensible heat (latent heat) fluxes and 2 m temperatures (2 m q) can still be observed. LUCC through perturbing land surface fluxes can influence local circulations, which can be further examined in the wet-case simulation. The thickness of the planetary boundary layer (PBL) indicates the mixing layer height; a warmer surface can enhance the mixing, thereby increasing the thickness of PBL ( figure 3). The PBL and temperature patterns thus bear a high resemblance. A warmer surface can also reduce air density and lead to a decrease in sea level pressure (SLP); lower SLP anomalies can be found over areas with increased urban and built-up lands. In accordance with the SLP anomalies, 10 m winds indicate local convergence taking place at the low SLP anomalies, and intensified onshore wind vectors can be located at certain areas (e.g. along the coast of Changhua County). The vertical profile (from 1000 to 800 mb) of the areal mean upward velocity in Taichung city also indicates a stronger upward motion in 2007 compared to that in 1995. Consequently, precipitation patterns have been altered: Major positive precipitation anomalies were identified not right above but in proximity to those areas where more urban and built-up lands were located and extended (figures 1 and 4(a)). These areas are mostly over the hilly areas, suggesting the underlying topography also plays an important role in altering precipitation patterns.
To manifest the combined effect of LUCC, local topography, and prevailing winds on precipitation, we overlaid the simulated precipitation anomalies with the wind field, along with the topography as background images, as shown in figure 4(b). We found intensified precipitation over the downwind areas of Figure 2. In dry-case simulations, the difference between the average of simulated fields driven by NLSC07 and that by NLSC95 land use, with respect to (a) sensible heat, (b) latent heat, (c) 2 m temperature, and (d) 2 m specific humidity. more urban and built-up lands, consistent with the major finding of [32] but having more details disclosed. The intensified precipitation patterns basically align over the windward side of the hills and specific basins, such as the Jialishan mountain in Miaoli, Taichung and Puli (Nantou) Basins, and Bagua Plateau in Changhua.

Observed trends in precipitation and relationship with urban expansion
Findings from the WRF simulation were further supported with the statistical assessments of observed datasets. As described in section 2.4, we first applied trend analysis to the Landsat land cover and gridded precipitation. Figure 5(a) indicates the observed trends in the percentage of urban areas derived from the Landsat data. In line with the NLSC data, several areas were identified with significant urbanization trends over central Taiwan, especially in the Taichung metropolitan area. In figures 5(b) and (c), central Taiwan exhibits mostly upward trends in the total cumulative as well as maximum daily rainfall in JA, based on the long-term precipitation records from 1960 to 2015. Certain areas marked with significant upward trends in the JA rainfall coincide with our simulated patterns, such as the Jialishan mountain and Taichung and Puli Basins previously mentioned. We also ensured that the field trend is significant using the false discovery rate approach [56][57][58]. Nevertheless, some coastal enhancement was also found, which should be a derivative of some other mechanisms than LUCC.
Next we performed the regression-based analysis to associate the interannual variations of precipitation    (2)), shown over grids where urban percentage greater than 10%. (c) As in (b) but P replaced from total to maximal rain. with LUCC. Following equation (2) for the off-grid association analysis that sifted through all E and SN grids can actually yield many significant correlations between P (total or max rain) and X (percentage of urban areas) over different grids. In figure 6, we presented selective results of dP dX, of which the average is highest over heavily urbanized areas (i.e. grids of urban areas >50% in figure 6 (a)). In both cases where the total and max rain as variable P, significant associations were identified for the eastward displacement of the P grid (E=4, ∼20 km). Thus, in agreement with the WRF simulation, the statistical assessment also suggests urbanization in central Taiwan has induced eastward, downwind intensification of the JA rainfall in the past decades.

Conclusions and recommendations
Central Taiwan has experienced rapid urbanization since the 1990s. Through numerical and statistical assessments, we concluded that changes in land surface conditions were responsible for the altered islandscale circulations subject to complex terrain. Our demonstration factored impacts of LUCC derived from realistic, survey-based land use data in the WRF simulation, producing significant differences in dry and wet cases under weak synoptic forcing in JA. The dry-case simulation suggested an urban heat island effect, resulting from further expansion of existing metropolitan areas. On the other hand, the wet-case simulation indicated that intensified precipitation could be associated with not only LUCC but also local topography and prevailing winds. Findings from the numerical experiment were corroborated by the statistical assessments of the Landsat land cover and gridded precipitation data: Significant increasing trends in urbanization and the JA rainfall and the off-grid LUCC-rainfall association were identified.
A more recent version of the NLSC land use data in 2015 (NLSC15) was just released and shown along with the NLSC95 and NLSC07 data in figure S4; the percentage change from the NLSC95 to NLSC07, and from the NLSC07 to NLSC15 were listed in table S3. Whereas there was a significant increase (16.5%) in irrigated cropland areas from 2007 to 2015, the change in urban and built-up land was minute (−0.8%). In contrast, the increase in urban and built-up land areas was much more dramatic from 1995 to 2007 (66.79%). This update suggests that the impacts of LUCC on central Taiwan's hydroclimate should have been made a decade ago, and we suspect that such impacts are still lasting due to similar land surface conditions. Urbanization along with other types of LUCC (e.g. irrigation and deforestation) can exert effects on land surface energy and water budget, further increasing the complexity in the Earth system under climate change. While our findings only verified the influences of LUCC in a particular region, various LUCC taking place in different regions suggests that global hydroclimatic conditions have been experiencing changes more rapidly than expected. More studies are required for the discussion of the consequence of LUCC associated with climate change in both regional and global contexts.