Constructing a 1km Gridded Multi-Scalar Drought Index Dataset ( 1960-2012 ) in Taiwan Based on the Standardized Precipitation Evapotranspiration Index-SPEI

An index sensitive to global warming, the standardized precipitation evapotranspiration index (SPEI), is employed in this study to construct a 1-km gridded multi-scalar drought index data bank in Taiwan. A siteand scale-dependent posterior fitness assessment procedure regarding determination of the most appropriate statistical distribution is used to standardize the station’s water deficit/surplus series, thereby SPEI at various time scales. Model uncertainty at different scales is evaluated and the results show that the uncertainty is higher for the shorter 10-day scale. Contrasts in the climatic means between SPEI and popular standardized precipitation index (SPI) are compared. The climatic pre-summer monsoon heating and midsummer drought phenomena are better captured by the SPEI. The established data bank, one of the data integration tasks in the Taiwan Climate Change Projection and Information Platform (TCCIP) project, is helpful in historical drought diagnostics. It also serves as the ground truth to correct model biases while using the regional model to project hydro-climate change and impact assessment, and historical baseline in the statistical downscaling while developing the statistical relationship between the general circulation model outputs and fine-scale observations. As an example of applications, the gridded SPEI at 3-month time scale is used to study the interannual variability in springtime drought. The results show that Taiwan’s southwestern plains region is most vulnerable to the risk of droughts. Composite analysis reveals that possible causes of island-wide drought link to the turnabout of cold El Niño-Southern Oscillation (ENSO) phase but also to the interannual variability in Pacific Decadal Oscillation when the analysis is initiated from the regional perspective.


InTroDuCTIon
One of the central missions for the Taiwan Climate Change Projection and Information Platform (TCCIP, hereinafter) project is to integrate raw station data distributed from various sources (see section 2) by which complete and high resolution gridded data can be produced via imputation and interpolation methods.Within the context of exploring regional climate change under global warming, the primary motivation for initiating the TCCIP project in 2009 involves three reasons for this request.Firstly, the importance of validating the regional climate model (RCM) output is increasing (Déqué et al. 2007;Fowler et al. 2007;Giorgi et al. 2009) due to the needs for impact assessments (IPCC 2014).The gridded data generated by the RCM represents the area averaged state rather than the point process.Constructing a mesh with each gridded value representing the best estimate average of the observations inside the grid is the most appropriate way for validating the model output rather than directly comparing it with the point observations (Osborn and Hulme 1998).Secondly, the downstream impact models often require complete gridded data to facilitate implementation with no further need to consider how to handle missing values (Haylock et al. 2008).Finally, the gridded data fill estimates at locations away from the point observations, thereby allowing local climate in data-sparse areas to be studied to some extent.Amongst others, precipitation and surface air temperature data are the foci of data collection due to their pivotal roles in determining the surface hydro climate, and thereby human activities.In this regard, drought due to precipitation deficit has been one of the major concerns for the TCCIP working group, especially when our world is warming (Dai 2011(Dai , 2013)).Drought affects more people than any other natural hazard (Obasi 1994) and is of great importance in the planning and management of water resources (Wilhite 2005).A long-lasting drought has dire socio-economic consequences.Different from the permanent aridity in arid regions, drought is a temporarily recurring extreme climate event (Dai 2011) and can occur in virtually all climatic zones, including monsoonal Asia.Since the mid-1990s, prolonged and widespread droughts have occurred in India and the occurrence frequency of droughts also increases as global warming progresses (World Bank 2003).Vulnerable drought-prone countries also include China.Severe droughts occurred in 1997 and 1999 to 2002 over northern China and resulted in large economic and societal losses (Zhang 2003).As in India, China has suffered an increased risk for droughts after 1980s as a warming world generates higher surface temperatures thereby increasing drying (Dai et al. 2004;Zou et al. 2005).Although the annual precipitation amount is high, Taiwan is also prone to droughts due to its uneven temporal precipitation distribution.Severe droughts occurred in 1963 to 1964, 1991, 2002 to 2003, and the first half of 2015.Based on the secular-long station records, it has also been detected that the long-term drought occurrence trend has increased in central and southern Taiwan since the 1960s (Chen et al. 2009).
Indexing is one among many ways for describing drought severity.The first comprehensive effort to evaluate the drought severity is the Palmer drought severity index (PDSI; Palmer 1965), based on a soil water balance equation.Utilizing precipitation and temperature records to calculate moisture supply and demand, PDSI and its ensuing variants have been widely used for monitoring (e.g., Karl and Quayle 1981;Jones et al. 1996) and forecasting (e.g., Kim and Valdés 2003;Özger et al. 2009) regional droughts, and assessing the global warming effect on the areal extent and severity of global droughts (Dai 2011).The primary limitation of PDSI is its inherent time scale making it more suitable for monitoring agricultural droughts rather than hydrological droughts (Guttman 1998).To mediate this drawback, McKee et al. (1993) suggested using the standardized precipitation index (SPI hereinafter) based on an appropriate probability distribution to monitor meteorological and hydrological droughts.Since then SPI has been popular worldwide (e.g., Edwards and McKee 1997 in USA;Lana et al. 2001 in Spain; Wu et al. 2001 in China; Mishra and Desai 2005 in India) and in local hydrological and meteorological societies (Shiau 2006;Chen et al. 2009Chen et al. , 2013)).
Based on the precipitation record length, SPI can be calculated for a desired period.Such versatility allows SPI to monitor changes in both short-term water supplies such as soil moisture and long-term water resources such as groundwater and reservoir levels.The major concern of SPI is that it relies on precipitation alone.Lots of zero precipitation values in a particular season are common for dry climatic regimes such as southern Taiwan during the northeasterly monsoon season.For shorter time scales, the calculated SPI series may not be normally distributed because of the highly-skewed underlying precipitation distribution.This may cause large biases while modeling precipitation distributions from limited data samples.With regard to the climate change issues, another major concern of using SPI is that it does not explicitly take the rising temperature factor into account.The local water demand, especially during the hot season, is critically dependent on the temperature via evapotranspiration.
Recently, Vicente-Serrano et al. (2010) proposed a new drought index, the standardized precipitation evapotranspiration index (SPEI hereinafter), to combine the multi-scalar character advantages of SPI with the PDSI capacity for comprising the temperature factor in the drought assessment.Because the generally rising temperature factor is explicitly included in its calculation (detailed in section 3), SPEI is sensitive to the possible effects of global warming on regional droughts at various time scales.Furthermore, it can also be used for detecting and attributing drought severity while projecting the future climate change scenarios.Although only recently developed, SPEI has been used to study droughts worldwide, for example: Australia (Allen et al. 2011); Europe (Potop et al. 2012); Africa (Harari and Ferrara 2013); and China (Yu et al. 2014).However, its use has not been propagated into the local meteorological as well as hydrological society yet, to the best of the author's knowledge.
To fulfill the aforementioned scientific missions assigned to the TCCIP project, the goal of this study is to construct a high-resolution gridded SPEI data bank at multi time scales over Taiwan and adjacent islands to facilitate the downstream drought related studies.After illustrating the sources of acquired precipitation and temperature data and data preprocessing in section 2, section 3 discusses the methodology of how to attain the goal.Specifically, a three-step procedure called posterior fitness assessment is proposed in this study to select the most appropriate statistical distribution for standardizing the water deficit/surplus series.Model uncertainty due to the chosen statistical models while constructing the station SPEI (and SPI) and the estimates of climates at different time scales are discussed in section 4. There, the contrasts in the climatic features between SPEI and SPI maps are highlighted.As an example of applying the established SPEI gridded data to the downstream applications, the interannual variability in springtime droughts in Taiwan is investigated and preliminary results are reported in section 5. Section 6 summarizes this study.

DATA SourCES AnD PrEProCESSInG
Through the data exchange platform of TCCIP at National Science and Technology Center for Disaster Reduction (NCDR), large efforts have been made to collect the rainfall and temperature data archived at various governmental agencies.These include the Central Weather Bureau (CWB) of Ministry of Transportation and Communications (MOTC) providing the rainfall and temperature data from synoptic stations (SYNOP) and auto rainfall observing network (ARON); the Water Resources Agency (WRA) of Ministry of Economic Affairs (MEA) providing rainfall observations from the hydro stations, rainfall and temperature data from the agronomy stations of Irrigation Association of Taiwan Province before 1998 (agronomy station data became void after 1999 due to the virtualization of Taiwan Province); the Civil Aeronautics Administration (CAA) of MOTC providing rainfall and temperature observations from civil airports, and the Air Force providing rainfall data at military airports.Few rainfall and temperature records in the mountainous areas are also obtained from the state-owned Taiwan Power Company when the company explored routes to construct the island-wide electric transportation network in the 1960s and 1970s.We discarded a few stations having record length less than two years.Using the 2-yr criterion does not seem to be rigid statistically.However, such loose criterion is based on the consideration that early data in mountainous areas are invaluable and worth being kept to facilitate the subsequent missing data imputation.
A serial sifting process such as deleting/combining repetition data, handling the inconsistent station coding of some agronomy stations, and manually correcting the obvious errors in the original handwritten records has been laboriously taken to assure data quality.Interested readers in the quality assurance and control of station raw data can refer to Hung (2012) for details.One thousand seven hundred and fifteen rainfall and 559 temperature stations were used in gridded SPEI/SPI dataset construction.See Fig. 1a for the rainfall station distributions and their geographical locations in Taiwan.As shown, although the density in the western plains of Taiwan is relatively high, stations are still scarce in the remote mountainous areas of the Central Mountain Range (CMR) when the elevation is higher than 2000 m.A similar data richness problem is fairly true for the temperature stations (not shown).
Missing data is unavoidable.The incomplete time series in both rainfall and temperature data found in both daily and monthly-mean records need to be imputed.The technique described in Simolo et al. (2010) (and earlier works: Shepard 1968(and earlier works: Shepard , 1984;;Willmott et al. 1985) is simple but has been proved to be efficient in high-resolution data analysis (e.g., Brugnara et al. 2012).This method is employed to perform the imputation of missing data in station time series.Briefly, the missing records at the target sta-tion are filled using the weighted averages of the products of three weights calculated at n neighboring reference stations having the concurrent records.Number n is set to 20 in the current study.The product of weights multiplies the relative radial, elevation, and angular weights of reference stations together.The introduction of angular weight, the angular separation of any two reference stations i and j with the vertex of the angle defined at target station, avoids the undesired over-weighting of areas with high station density (Willmott et al. 1985;Efthymiadis et al. 2006).We chose the top 10 reference stations with the highest products of three weights to calculate their weighted average.Interested readers in the data imputation technique can refer to Simolo et al. (2010; also see Weng and Yang 2012 for its recent application in Taiwan) for details.
As mentioned, number of the temperature stations is far less than the number of rainfall stations.However, we need temperature data to coexist with the rainfall data at the same station to calculate the SPEI value.Spatiotemporally, temperature is in general a more homogeneous climate variable than rainfall.The error caused by the interpolation scheme is presumably smaller.Therefore, for those rainfall stations lacking temperature records, we first apply the aforementioned missing data imputation technique to fill the temperature gaps before calculating the station SPEI.Complete rainfall and temperature series (from 1 January 1960 to 31 December 2012) were therefore imputated for 1715 stations.

SPEI/SPI Calculations
The first step in calculating the SPEI is to evaluate the potential evapotranspiration (PET hereinafter).This is difficult because PET is in fact determined using many parameters such as humidity, soil wetness, insolation, vapor pressure, surface-atmosphere latent heating, and sensible heating, in addition to the ground temperature (Shuttleworth 1993;Allen et al. 1998).The popular Penman-Monteith method, recommended by the Food and Agriculture Organization (FAO) of the United Nations as the standard method for computing PET, cannot be used in the current study because the required solar radiation data is not available at most stations.Instead, we followed the simple method described in Thornthwaite (1948) to calculate PET.This method, which only needs the surface temperature and station's latitude as input, has been shown able to provide similar results as more complex methods when a drought index is calculated (Maurer 2007).Although Hobbins et al. (2008) found that the PET estimated using the Thornthwaite method could lead to errors in energy-limited regions; this should not be a critical problem in subtropical Taiwan.Since rainfall and surface air temperature are the only two climate variables with long historical records and the latter is the main theme in a warming world, the SPEI can make full use of these data to test drought sensitivity to global warming.Vicente-Serrano et al. (2010) provided the step-by-step procedure using Thornthwaite method to calculate PET.
With the monthly (and daily) PET data available, the monthly (or daily) series D i k consisting of the aggregation of differences d i between the precipitation P i and PET i , d i = P i -PET i , can be formulated as where i is the monthly (or daily) index, k is the chosen time generated by the TCCIP working group.In this report, we will mainly discuss the datasets with 10-day, 1-, and 3-month time scales.The standardization for the D i k (or D * i k ) series is comprised of a transformation of a statistical distribution appropriately describing the long-term series of observations to standard Gaussian distribution.McKee et al. (1993) used the two-parameter Gamma distribution to develop SPI.Recently, the two-parameter Gamma distribution was also adopted by Chen et al. (2013) to construct SPI and develop probabilistic drought forecasting in southern Taiwan.Guttman (1999) recommended the three-parameter Pearson type III distribution (sometimes referred to as the 3-parameter Gamma distribution) while calculating SPI.Since negative values (i.e., water deficit) exist in the D i k series during the dry seasons, the two-parameter Gamma and three-parameter Pearson type III distributions commonly used in standardizing the D * i k series (having zero lower bound) may not be adequate to standardize the D i k series to obtain the SPEI.Vicente-Serrano et al. ( 2010) suggested the three-parameter log-logistic distribution for standardizing the D i k series worldwide.However, the log-logistic distribution drawback is the existence of statistical moments for a very limited range of distribution parameters (Rowiński et al. 2002).Therefore, it seems that the universally "best" distribution does not exist.
This study proposes a posterior fitness assessment procedure to select the most appropriate statistical distribution for standardizing D i k and D * i k series.Nine theoretical distributions including the 2-parameter Gamma (GAM), 2-parameter Gumbel (GUM), 3-parameter Pearson type III (PE3), 3-parameter log-normal (LN3), 3-parameter generalized logistic (GLO), 3-parameter generalized extreme-value (GEV), 3-parameter generalized Pareto (GPA), four-parameter Kappa (KAP), and five-parameter Wakeby (WAK) distributions, are tested first to model the D i k values on the chosen time scale and at the specific station.We use the Lmoment method (Hosking 1990; Fortran routines are available at the Internet location http://lib.stat.cmu.edu/general/lmoments) to estimate the parameters for the above distributions.The L-moments, though analogous to the conventional central moments, are able to characterize a wider range of distributions and are more robust to the outliers (Hosking 1992;Ulrych et al. 2000).
After the parameters are estimated the corresponding and determined cumulative distribution function F(x) of a specific distribution, where random variable x represents the observed D i k ( D * i k ) series, is inversely mapped to the quantile function (i.e., the probit function called by statisticians) of standard Gaussian distribution z p , ( ) to get the SPEI (SPI).The erf -1 is the inverse function of error function.The z p , namely the desired SPEI (SPI) series, has no analytical expression.We use the popular formula of Abramowitz and Stegun [1965, Eq. (26.2.23)] to numerically approximate the z p as follows: where ln t p 1 2 = ^h, when 0 < p ≤ 0.5.When p > 0.5, the @ .The coefficients c 0 = 2.515517, Finally, two L-moment ratios: L skewness 3 x and L kurtosis 4 x , of the estimated SPEI/SPI (i.e., the transformed z p series) can be calculated as 3 3 2 x m m = , and x m m = , where 2 m , 3 m , and 4 m are the L-moments obtained from the probability-weighted moments (PWMs; cf.Greenwood et al. 1979;Wallis 1982;Haktanir 1991).The relationships between the L-moments and PWMs are as follows: The PWMs of order s, w s , s = 0, 1, 2, 3, are calculated as where i is the range of observations arranged in ascending order, N is the sample size, SPEI i (or SPI i ) is the transformed z p series, and f i is the distribution-independent plotting-position estimator (cf.Hosking and Wallis 1995).Hosking and Wallis (1997) suggested the Landwehr formula (Landwehr et al. 1979) , to get the reasonable results.
If the observed D i k series is well described by the selected statistical distribution, the transformed z p series should then approximate the standard Gaussian distribution which has L skewness 0 x x Y Y are the estimated ( , ) 3 4 x x by the specific distributions.As an example, Fig. 1b shows the selected statistical distributions at individual stations for the SPEI on the 1-month time scale.Clearly, a unified distribution does not exist.Instead, stations at high mountainous regions prefer PE3 (orchid dots) and KAP (red dots) distributions, whereas many stations in the western plains choose GLO (silver dots) and WAK (blue dots) distributions.The LN3 (aqua dots) and GEV (green dots) distributions are mainly found at stations in the foothills of CMR.
Choosing a distribution with many parameters seems to be odd as the basic precept in statistical modeling is to keep the number of parameters as low as possible.Nevertheless, since the aim of the current study is to provide a high quality SPEI databank facilitating the use of well-established Gaussian process methods, and the estimated D i k ( ) D * i k series, if necessary, can be easily obtained through the analytic quantile functions of these theoretical distributions, the concern should not be critical.
In principle, the D * i k series at stations can have its own statistical distribution while constructing the station SPI.However, we deliberately used the same distribution decided by the D i k series.In this way, any differences between SPEI and SPI series are attributed to the temperature variations (through PET) alone.But, how good (or bad) is such an approach?Figures 2a to d respectively show the Lmoment ratio diagrams (abscissa: L skewness 3 x ; ordinate: L kurtosis 4 x ) for station D, D * , SPEI, and SPI series calculated at 1-month time scale.The estimated ( , ) x x Y Y of D and D * series are obtained by employing the ascending-ordered D and D * series in Eq. ( 5).The theoretical L-moment ratios for two 2-parameter (exponential and Gaussian distributions) and five 3-parameter (PE3, LN3, GLO, GEV, GPA) distributions are also drawn as points and colored lines, respectively.Note that the black point (0.0, 0.1226) represents the standard Gaussian normal distribution.Before standardization, large diversity in the estimated ( , ) and D * series is evident (blue dots), especially in the D series (Fig. 2a).This is understandable due to the add-on complex of temperature variations.Note that many blue points in Fig. 2b cluster around the PE3 curve (dashed olive line).
It endorses the popularity of the PE3 distribution adopted by the local hydrological society in the SPI calculation.After standardization, the estimated ( , ) x x Y Y of SPEI (Fig. 2c) and SPI (Fig. 2d) series converge into the surroundings of Gaussian point, indicating the normality.As shown, the degree of convergence in the SPI series is comparable with that in the SPEI series.While a distribution with many parameters is requested by the more complex D series on the 1-month time scale, it should also satisfy the need of simpler D * series, though possibly redundant.Panels in Fig. 3 show the similar L-moment ratios diagrams for the 10-day time scale.At this shorter time scale the clustered points run away from the Gaussian point more, indicating that both D (Fig. 3a) and D * (Fig. 3b) series become more non-Gaussian.Notice that the ratio diagrams recommend the GLO (PE3) distribution for a large portion of the D (D * ) series.Using a distribution suitable for the D series on the D * series would cause some inconveniences in the degree of normality of SPI series (Fig. 3d).Nevertheless, the Gaussian requirement for the SPEI series at the 10-day time scale is also less satisfied (Fig. 3c).The problem is thus not merely attributed to the chosen statistical distribution.Rather, it likely points out the limitation of using the L-moment method for parameter estimation.The robustness of L-moments to the outliers also indicates their lesser sensitivity to the heavy tailed distributions and thus modifications on the L-moments have been suggested (e.g., Elamir and Seheult 2003;Hosking 2007).The analysis of heavy-tailed data is a critical issue while moving towards the higher frequency end.Therefore, the use of SPEI/SPI, currently based on the L-moments method, at 10-day time scale should be judged with cautions.
In contrast, the aforementioned problem largely disappears for the 3-month time scale (Fig. 4).Both D and D * series become more Gaussian.Here, many D * series show their preference towards the GPA distribution (Fig. 4b), whereas the D series show no such preferences (Fig. 4a).Nevertheless, both SPEI and SPI series achieve normality, as shown in Figs.4c and d.

Griddization
A built-in NCL (NCAR Command Language; NCAR: National Center for Atmospheric Research) package, Natgrid, is used in this study to interpolate station SPEI/SPI data from an irregularly-spaced distribution to a rectilinear 1-km resolution horizontal mesh.Natgrid implements a natural neighbor interpolation method (Sibson 1981) based on Watson's package nngridr (Watson 1992(Watson , 1994)).One distinguishing quality of the natural neighbor interpolation method is the way in which a set of neighboring points (i.e., the natural neighbors) is selected for interpolating at a given point.The selection process avoids the problems common to methods based on choosing a fixed number of neighboring points or all points within a fixed distance.Natgrid uses a weighted average method that is more sophisticated than the usual inverse distance weighted average method.Another distinguishing quality of this method is the way that the weights are calculated for the functional values at the natural neighbor coordinates.These weights are based on proportionate areas rather than distances.Note that the current package does not allow any missing values in the input data.Readers interested in the theoretical/technique details should refer to Watson's books (1992Watson's books ( , 1994)).
The above griddization method is implemented in the empirical orthogonal function (EOF) space rather than directly in the physical space.The orthogonality property of EOF provides a complete basis where the time-varying SPEI (or SPI) field can be decomposed into Where T k and S k , the k'th principal component and the k'th Eigen pattern, are time-dependent and space-dependent alone, respectively.The above interpolation scheme is first implemented on each S k pattern to convert it into G k mesh (G denotes gridded level) and then summed up through Eq. ( 6) where S k is replaced by G k to obtain the gridded SPEI (or SPI).To alleviate/remove the local noise, the summation is terminated when the leading EOFs [denoted as symbol N in Eq. ( 6)] together explain the 99% of the total variance amount.

uncertainty
If the calculated SPEI/SPI series were truly Gaussian, the mean ( n ) of the series is then zero and the standard deviation (v) is 1, i.e., SPEI/SPI ~ N(0, 1).Note that a SPEI (SPI) of zero indicates a value corresponding to 50% of the cumulative probability of D (D * ) value in the selected distribution.To examine the expected normality, Figs. 5 and 6 respectively show the spatial distributions of the estimated means ( n X ) and standard deviations ( v X ) of gridded SPEI/ SPI at 10-day, 1-, and 3-month time scales.The n X of SPEI at 10-day scale shows a slight dry-bias to the southwest of Taiwan, whereas it expresses a slight wet-bias in the southern tip of CMR at both 1-and 3-month scales.The former one (latter ones) is (are) associated with slightly over-(un-der-) estimated v X in the surrounding areas.Similar slight wet-bias associated with underestimated variability is also found to the southwest of Taiwan for SPI at 1-and 3-month scales.At 10-day scale, however, the n X ( v X ) of SPI express- es a larger spread of wet-bias (underestimated variability) in the western plains interfered with "bulls eye" of dry-bias (overestimated variability) in its southwestern corner.In addition, a slightly underestimated v X appears in the northeast- ern coast for SPEI at 3-month scale and for SPI at all three scales.Conclusively, larger uncertainty is mainly found at the 10-day SPEI scale as well as SPI, and over the west/ southwest regions.Note that these regions have distinguishing dry and wet seasons around the year.
To further examine the model uncertainty at the seasonal time scale, panels in the left, middle, and right columns of Fig. 7 show the mean monthly-mean differences averaged in January, April, July, and October (from-top-to-bottom) between the estimated D (i.e., D X ) and observed D series ( ) D D + X at 10-day, 1-, and 3-month scales, respectively.The D X series are obtained from the quantile functions of site- dependent statistical distributions.Regardless of which time scale, larger model errors occur in the warming months and generally over the remote mountainous areas.This implies that either the selected statistical distribution is still inappropriate or the imputed precipitation and/or temperature exists larger uncertainties therein (cf.Weng and Yang 2012).Furthermore, the employed distributions tend to overestimate the D series (i.e., wetter).In April, +0.5 ~ +1.0 mm day -1 wet bias can be found in the southern portion of CMR at 1-and 3-month time scales.In July, the wet-bias reaches +2.0 mm day -1 over the western plains (southwestern plains and East Rift Valley) at 10-day (1-month) time scale.Note that at 3-month scale a dry-bias with -3.0 mm day -1 error also appears in the foothills of southern part of CMR.Different panels in Fig. 8 shows similar results but examines the mean monthly-mean differences between the estimated D * (i.To explore the effects of interannual temperature variability on the regional water balance, Fig. 9 shows the mean seasonal correlations between the SPEI and SPI series at the 1-month scale.Lower (higher) correlations indicate larger (smaller) interfering effects of temperature changes.To highlight the regions sensitive to the temperature factor, the correlations were subtracted from 1.0.The sensitive regions extend northward from the southern tip of Taiwan in spring (March-to-May; MAM) to the entire western plains in fall (September-to-November; SON).The sensitivity largely increases over the southwestern plains in winter (Decemberto-February; DJF).In DJF, water balance in the southeastern regions is also interfered with by the temperature variability.The decreasing (increasing) temperature effect in the warming MAM and JJA (cooling SON and DJF) seasons implies the more active (passive) role of rainfall processes

Climatology
Figures 10 and 11 respectively show the spatial distributions of the climatological-means in the SPEI and SPI at 10-day time scale.Only the means in the 1 st and middle days of each month are presented.Since both SPEI and SPI series at specific station are generated using the same statistical distribution, any differences between them are caused by the temperature factor alone (via PET).Although resembling each other, the evolutionary wet/dry climate described by the SPEI and SPI patterns indicates that the degree of wetness (dryness) in the period of 15 April -15 October (1 November -1 April) decreases after taking the temperature factor into account.The reduction in wetness in the warming period due to the increased PET is found mainly in the southwestern regions facing the prevailing southwesterly winds, whereas the regions with reducing dryness in the cooling period due to the decreased PET are located in the northeast/north coasts in early winter facing the prevailing northeasterlies and later in the northwest/west coasts facing the prevailing northwesterlies.
The SPEI patterns feature drier southwestern plains in the 1 May and 15 May phases before the South China Sea (SCS hereinafter) monsoon (i.e., Meiyu-Baiu) onset, the socalled pre-summer monsoon heating phenomenon, and the subsequent fast transition into the wet phase.Furthermore, it also highlights the mid-summer drought phenomenon during the monsoon break (Asai and Fukui 1977;Ho and Wang 2002;Karnauskas et al. 2013) in July.Both climatic phenomena cannot be clearly addressed by the SPI patterns.
One of the distinguishing SPEI (and SPI) features is that it can be calculated at different time scales to monitor drought and moist conditions in different hydrological subsystems (Vicente-Serrano et al. 2010).At a shorter (longer) time scale, the drought and moist periods have shorter (longer) duration but a higher (lower) occurrence frequency.It is also important to understand how the temperature factor, coming through PET, is involved with different time scales for soil wetness, river discharge, reservoir storage, and even the long-term groundwater variability.Figure 12 shows the climatic SPEI and SPI patterns in January, April, July, and October at the 1-month time scale.It is immediately recognized that the degree of water surplus/deficit and dry-to-wet contrast described by the monthly SPEI pattern is more severe than what we expect compared to SPI as well as SPEI at the 10-day time scale.For example, the mid-summer drought in July (cf.Fig. 10) is in fact more durable in northern Taiwan.In addition, contrasts between SPEI and SPI patterns show that using the SPEI degree of local wetness described by SPI in July decreases, whereas the degree of local dryness in April and October increases over the southwestern and western parts of Taiwan.In January, however, the degree of dryness decreases island-wide in the mid-winter.
Similar contrasts between SPEI and SPI patterns are also evident in the 3-month time scale (Fig. 13).Noticeably, the dry-to-wet contrast between the western (southwestern) and eastern (northwestern) Taiwan in November (May), which represents the overall condition during the fall (spring) season, is emphasized by the SPEI pattern in this seasonal time scale.In February, the SPEI pattern also highlights the "Dongyu" (freezing wintertime rainfall in Chinese) phenomenon in the northeast horn of Taiwan comprising the Northeast/Yilan National Scenic Area and Lanyang Plain, facing the prevailing northeasterly monsoon winds.

APPlICATIon: InTErAnnuAl vArIABIlITy of SPrInGTIME DrouGhTS
With long-term (1960 -2012), multi-scalar, and 1-km gridded SPEI/SPI datasets becoming available, researchers in various disciplines are ready to conduct various applications related to the high-resolution analyses such as verification of regional hydrological aspects associated with historical typhoon events (e.g., Hsiao et al. 2013;Shih et al. 2014) and interactions between vegetation and climate system (e.g., Hickler et al. 2005;Heumann et al. 2007;Jain et al. 2009).Utilizing the established datasets, tasks are underway inside the TCCIP working group.As an example, preliminary results regarding the springtime drought problem are briefly reported in this section.

regional Aspects
The interannual variability in spring rainfall in Taiwan is largely associated with El Niño-Southern Oscillation (ENSO) (Wang and Hwu 1994;Chen et al. 2003Chen et al. , 2008Chen et al. , 2013)).Analyses conducted in previous studies still relied on the station-level data and often attempted to identify the top-down causal relationship between the large-scale lowfrequency variability as the cause and the regional hydroclimate response.To provide a holistic spatial distribution perspective regarding the regional variability in dry/wet conditions during the springtime, anomalies (respect to the 1960-to-2012 mean) of gridded SPEI with the 3-month time scale in May, which represents the overall MAM condition, are subject to EOF analysis (EOFA).To test the high-resolution analysis robustness, we also used the SPEI series at 28 CWB/SYNOP stations (cf.black dots in Figs.14b or 15b for their geographical locations) to conduct the same EOFA.The raw rainfall and temperature records at these operational stations are relatively complete.In both high and low resolution analyses, only the first two EOF modes are statistically significant according to North's rule of thumb (North  et al. 1982) and individually explain 69.8 and 11.4% (60.6 and 12.2%) of the total variability of the gridded (station) anomalies, respectively.
Figure 14 shows the first EOF (EOF1) patterns for the gridded (Fig. 14a) and station (Fig. 14b; after spatial interpolation) anomalies and the corresponding time series principal components (PC1s; Fig. 14c).The correlation coefficient of these two PC1 series between high (histogram) and low (red line) resolution analysis is 0.98, indicating a highly coherent time evolution.The EOF1 pattern of the low-resolution analysis shows a smooth structure with larger (smaller) weights located in the southwestern plains (northeast horn) of Taiwan.The central portion of mountainous Taiwan also has a smaller weight (cool colors), whereas the northwestern region shows a larger weight (warm colors).While the gross feature is similar, the high-resolution analysis expresses a much finer spatial structure.Centering in the middle between Kaohsiung and Tainan district (cf.Fig. 1c), EOF1 pattern puts more weights on the entire southwestern and western plains suggesting that these regions are most vulnerable to the threat of springtime drought/flood in some years.The larger vulnerability also covers the northwestern region including the Shei-Pa national park to its southeast.Note that the smaller weight previously in Central Taiwan now shifts to the southern part of the CMR but the one previously in the northeast horn keeps stagnant.
The second EOF (EOF2) patterns for the high-and lowresolution analyses and the corresponding PC2 are shown in three panels of Fig. 15.The correlation coefficient of two PC2 series between the high and low-resolution analysis reaches 0.87, significant at 0.5% level.Unlike the islandwide coherence in the EOF1, spatial distribution of EOF2 pattern in both high (Fig. 15a) and low (Fig. 15b) resolution analyses expresses a north-to-south dipole structure of negative-to-positive SPEI anomalies indicating the coexistence between a wetter south and a drier north.As shown in Fig. 15c, the SPEI magnitude in some particular years is smaller (larger) than -0.84 (0.84), reaching the moderate drought (flood) severity according to Agnew's classification (Agnew 2000) and thereby implying a 'southern flood (drought) vs. northern drought (flood)' pattern.For the highresolution analysis, the negative SPEI anomalies (shaded in warm colors) centering in the Taipei Valley (cf.Fig. 1c) is accompanied by the positive anomalies (shaded in cool colors) centering in the southwestern Pingtung Plains.The former stretches southeastward to the upper part of the East Rift Valley, whereas the latter stretches northwestward towards the western plains.Such detailed stretching is almost absent in the low-resolution analysis.In addition, there is a mismatch of the positive anomalous center between the low and high-resolution analysis.
We can use the PC series to assess the influence of temperature on the dry (or wet) severity.This is achieved by correlating the anomalous gridded SPI series with the PC series obtained from the EOFA of gridded SPEI anomalies.The resultant correlation maps are shown in Fig. 16.Since both SPEI and SPI share the same statistical distribution, the regions with lower (higher) absolute values of correlations imply a higher (lower) influence of temperature on the relative dry/wet severity.The resultant correlation map of EOF1 (Fig. 16a) indicates that the influences of temperature are mainly in the plains areas including the southwestern Pingtung Plains, East Rift Valley, Taipei Valley, and northeastern Lanyang Plains, consistent with the contrasts in climatic settings between the SPEI and SPI in May (cf.Fig. 13).In contrast, the relative dry/wet conditions over the mountainous areas and western plains north of Tainan are determined mainly by the absence/presence of precipitation but less affected by the temperature variability.
In contrast, the correlation map of EOF2 (Fig. 16b) shows a coherent relationship between the positive (negative) anomalies to the southern (northern) Taiwan and the superimposed positive (negative) correlations.Near the positive/negative action centers, larger anomalies are associated with higher correlations (in absolute values), suggesting that the north-south diploe structure is primarily attributed to the rainfall processes and the temperature variability only plays a minor role.

large-Scale Associations
To shed some light on seeking the potential sources for the variability of regional wet/dry conditions in spring, we conducted composite analysis on the seasonal-mean sea surface temperature [SST; data source: Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST); Rayner et al. 2003], low-tropospheric (850 hPa) winds, and associated stream function anomalies using data reanalyzed by the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) (Kalnay et al. 1996).The composited difference is calculated between the ensemble averages of the anomalies in the selected 7 driest springs and in the 7 wettest springs (according to the PC series) before being divided by two.The Monte Carlo white noise test (Overland and Preisendorfer 1982) with 5000 resampling among non-repeated 53 years (1960 -2012) is used to test the local SST significance and vector wind anomalies at the 0.1 confidence level.
The resultant composite corresponding to the EOF1 of regional SPEI anomalies is shown in Fig. 17a.When an island-wide drying event occurs (cf.Fig. 14a), Taiwan as well as southeastern China is well-situated under the influence of an anomalous anticyclonic circulation.Furthermore, the associated east-northeast wind anomalies with downslope subsidence are prone to the formation of dry condition over the southwestern-to-western plains of Taiwan after crossing the CMR.As shown, the above regional anomalies are in fact only a part of large-scale circulation anomalies posited in the vast western and northern Pacific oceans.To the northeast direction, an anticyclonic gyre in the North Pacific Ocean is associated with warm SST anomaly (SSTA) in the central Pacific and cold SSTA on the west coasts of the North American Continent, resembling a low-frequency Pacific Decadal Oscillation pattern (PDO; Mantua et al. 1997; see Mantua and Hare 2002 for a review) in its negative phase.
The atmospheric-oceanic anomalies configuration herein signifies the weakened Aleutian low, southwestward displacement of western North Pacific Subtropical High (WNPSH), and weakened East Asian jet aloft (not shown).The rainfall activity associated with the frequent westerly disturbances within the mid-latitude frontal system is thus retarded over the south-to-southeast China, Taiwan, and north Philippines,   but enhanced over the Yangtze River-to-Yellow River region in east China, Korea, south Japan, and south Philippines, according to the composite analysis using APHRODITE (stands for Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation of Water Resources, cf.Yatagai et al. 2012) data (not shown).
To the southeast direction, an anomalous cyclonic circulation in the Philippine Sea-to-tropical western Pacific is accompanied by a wedge-shaped warm SSTA pattern with equatorial symmetric maximum appearing near the Niño 4 region (averaged SSTA in the 5°N -5°S, 160°E -150°W domain).Meanwhile, cold SSTA appear in the far eastern tropical Pacific south of equator.Notice that the Philippine Sea low-level cyclonic anomaly has a 'little brother' over the Solomon Islands in the southern hemisphere and is associated with a upper-level anticyclonic anomaly (not shown), suggesting the Matsuno-Gill-type pattern (Matsuno 1966;Gill 1980) in response to the equatorial warm SSTA forcing west of the dateline.The low-latitude atmospheric and oceanic scenarios just described are reminiscent of the tropical Pacific-East Asian teleconnection pattern during the turnabout of cold ENSO phase (Wang et al. 2000).However, a classic quadrupole response of the Gill's atmosphere to the equatorial heating is obscure east of the dateline and is asymmetric with respect to the Equator west of the dateline.Meanwhile, the presumed cold (warm) SSTA to the northwestern flank of Philippine (anti) cyclonic anomaly due to the local air-sea interactions (Wang et al. 2000) is statistically insignificant.
As a result, when we initiate the analysis from the regional perspective, a hybrid large-scale teleconnection pattern is likely formed from the combination of tropical and mid-latitudinal anomalous circulations.In other words, the general notion: the excess (deficient) spring rain in Taiwan can be attributed to the strengthened southwesterly (northeasterly) winds induced by the anomalous Philippine anticyclone (cyclone) which is maintained by the local air-sea interactions during the decaying stage of warm (cold) ENSO phase, needs to be scrutinized (cf.Chen et al. 2008).
As mentioned previously, the composited SSTA distribution in the North Pacific resembles the negative phase of the low-frequency PDO teleconnection pattern.The seasonal PDO index, which is defined as the first PC of the EOFA of seasonal SSTA in the North Pacific north of 20°N (Zhang et al. 1997), is calculated and then correlated with the SPEI-PC1 series.While the correlation in the lag-1 season is -0.32, it reaches -0.43 in the lag-0 season which is significant at the 0.1-level.It thus suggests that except for the tropical influences, which seem to be case dependent, the interannual variability in circulation and SST in the North Pacific regime in the preceding seasons could also play some roles in affecting the springtime dry/wet conditions in Taiwan, though its detailed process needs to be further explored.
Figure 17b shows the composited large-scale SST and low-level circulation anomalies related to the 'southern flood vs. northern drought' regional pattern (cf.Fig. 15a) described by the EOF2.To the southwest (northwest) of Taiwan, the low-level west-northwest (west-southwest) wind anomalies appear in the northern part of the SCS (Yangtze River Valley).The associated upper-level trough deepens in the eastern coasts of the Indochina Peninsula, the coasts of southern China, and the northern part of the SCS (not shown), thereby enhancing the rainfall activity in the southwestern corner of Taiwan.Increased rainfall also appears in the adjacent regions of Vietnam and northern part of the Philippines based on the APHRODITE analysis (not shown).The cold SSTA therein are in response to the atmospheric forcing via the Ekman pumping.To the northeast (southeast) of Taiwan, an anomalous anticyclonic (cyclonic) gyre is found in the western subtropical (tropical) Pacific.Between them, the bifurcate easterly anomalies with the ridge line stretching westward towards the east coasts of Taiwan are well-posited near the northeast horn of Taiwan, resulting in the rainfall deficits there.The decreased rainfall signal also extends westward into Fujian and Zhejiang provinces over southeast China (not shown).
In the western Pacific sector, a conspicuous wave train pattern with north-south orientation emanates from the equatorial Pacific near the dateline.Note that another wave train pattern with the opposite anomalous cyclone and anticyclone phase is juxtaposed in the eastern Pacific north of 20°N.The above anomalous circulations in the pan-North Pacific basin vertically express a barotropic (baroclinic) structure north of 20°N (within the deep tropics), whereas the Philippine Sea cyclonic anomaly tilts northwest accompanying the upper-level cyclonic circulation over the northern part of the SCS (not shown).Within the low-latitudes equatorial westerly anomalies which suppress the upwelling and thicken the thermocline depth could cause the warm SSTA (≥ 0.5°C) in the central to western tropical Pacific.The cold SSTA strip near the northern Marianas may be induced by various processes or their combination.The cyclonic anomaly may raise the off-equator thermocline via Ekman pumping leading to the decrease in SST.The enhanced northerlies in the northeastern flank in this season also reduce the SST via wind-evaporation feedback.It is also likely that such cold SSTA, which extends westward to the east coasts of the Indochina Peninsula and eastward to the north of the Hawaiian Islands, is ignited in response to the atmospheric forcing and signifies the equator ward shift in mid-latitude westerlies in East Asia.
It is intriguing that the western Pacific wave train attains a spatial pattern similar to the positive summertime Pacific-Japan (PJ) teleconnection pattern phase (Nitta 1987).The positive (negative) active centers of height anomalies are near Japan and the Aleutian Islands (over the Philippines and Okhotsk Sea).The summertime PJ pattern has been interpreted as a remote response to the basin-wide warming SSTA over the Indian Ocean after El Niño, which suppresses the convection over the Philippines and enhances the Meiyu/ Baiu rain band in East Asia (Xie et al. 2009).However, our springtime version of the PJ-like pattern shows its association with the SSTA in the central tropical Pacific Niño 4 region rather than in the Indian Ocean.The associated rainfall activity shows a sandwiched pattern, namely the suppressed rain in southeast China (including the northeast Taiwan) is accompanied by the enhanced rain in the Vietnam-Philippines to its south and in the Yangtze river basin-Korea-Japan to its north (not shown), which is not fully consistent with the typical summertime PJ related rainfall pattern.Although four out of the chosen seven years (1966, 1969, 1977, 2010;cf. Fig. 15c) in the positive phase of composite analysis are indeed after weak-to-strong El Niños, two out of seven (1972,2002) are actually before moderate-to-strong El Niños.The correlation between PC2 series and concurrent (lag-1 season) Niño 4 index is only 0.24 (0.10).The correlation between the PC2 series and summertime PJ index, which is commonly defined as the difference in 850-hPa geopotential height anomalies between two grid points (155°E, 35°N: region east of Japan and 125°E, 22.5°N: region east of Taiwan) (Wakabayashi and Kawamura 2004), is only 0.13.None of them are statistically significant.Whether the detected springtime PJ-like pattern represents a physical mode and its possible linkage with other known climate modes is worth further exploration.

SuMMAry
A complete (without missing values), long-term (1960 -2012), and high-resolution (1-km) gridded data bank including both SPEI and SPI at various time scales (from 10-day to 2-yr) has been constructed as one of the data integration tasks in the TCCIP project.Such tasks will not be materialized without the pivotal role played by the NCDR/ MOST in collaboration with several governmental agencies during the data collection phase (cf.Fig. 1a).Only the results at 10-day to 3-month scale are shown in this study.Apart from facilitating the diagnostics of climatic processes related to the regional hydro-climate variability, the established quasi-observed gridded data, which use gauge records alone, serve as the ground truth to correct the model biases while using state-of-art dynamical RCMs to project future hydro-climate change and impact assessment in the regional scale.It also serves as the historical baseline in the statistical downscaling while developing the statistical relationship between the large-scale GCM outputs and fine-scale observational data.Moreover, the introduced SPEI uses both rainfall and surface air temperature as input, which is in contrast to the conventional SPI that is based on rainfall alone.This would allow the SPEI to account for the increased warming almost certainly in the coming decades (Solomon et al. 2007;IPCC 2014) on the devastating risk of drought in the East Asia (Dai 2011;2013).The established gridded SPEI/SPI datasets are accessible for interested researchers through the data service platform at NCDR: http:// tccip.ncdr.nat.gov.tw/ds/.
After acquiring the complete station rainfall and temperature series via the missing data imputation technique (Simolo et al. 2010) and calculating PET via the Thornthwaite method (1948) to obtain the D series, this study proposes a site-and scale-dependent posterior fitness assessment procedure to standardize the D series and thereby SPEI.The procedure is comprised of 3 steps to attain the optimal standardization.Firstly, the specific D series is modeled using nine commonly used statistical distributions (GAM, GUM, PE3, LN3, GLO, GEV, GPA, KAP, and WKA).There, the L-moment method is employed to estimate the parameters associated with individual distributions.Secondly, the estimated parameters are used to obtain distribution-dependent SPEI after inversely mapping the associated CDF to the probit function of standard Gaussian distribution.Thirdly, the distribution-dependent SPEI series are used to calculate the corresponding L skewness 3 x and L kurtosis 4 x via the Lmoments based on the PWMs.The distribution that has the smallest Euclidean distance between the estimated , ( ) x x Y Y and (0, 0.1226), the theoretical ( , ) 3 4 x x of the standard Gaussian distribution, is then chosen as the most appropriate one to standardize the D series (cf.Fig. 1b).Such a strategy intends to provide a near-normal SPEI to facilitate using well-established methods (such as EOFA) for the Gaussian process.Note that in this study we used the same statistical distribution to construct both the SPEI and SPI at the same stations and scales.Through direct comparison between them, possible temperature variability effects on the dry/wet condition at different time scales can be highlighted.
The L-moment ratio diagrams  show that the uncertainty of the employed statistical models is higher at the shorter 10-day time scale in both SPEI and SPI series which is in response to the larger non-Gaussian D and D * series.There, only the marginal normality can be attained.The cause is likely due to L-moment method use in the parameter estimations of heavy-tailed distributions and future refinement is needed (Hosking 2007).Spatially, SPI (SPEI) at this time scale has the wet (dry) bias with underestimated (overestimated) variance mainly over the western (southwestern) plains .Seasonally, SPEI at different time scales all express the overall wet bias with maximum around 2 -3 mm day -1 in spring and summer seasons (Fig. 7), whereas SPI mainly shows larger dry bias in summer with maximum around 4 mm day -1 over the southern part of CMR (Fig. 8).In the worst scenario, this is roughly equivalent to 15 -20% relative error with respect to the mean seasonal accumulation over the data scarce regions.Figure 9 shows that the southwestern-to-western plains, especially in SON and DJF seasons, are sensitive to the temperature effect while considering the interannual variability in regional water balance.This is inconsistent with the finding in section 5 when studying the springtime variability (cf.Figs. 14 and 16a).
The observed mean wet/dry climate contrasts between the SPEI and SPI patterns at the 10-day time scale (Fig. 10 vs. Fig.11), namely the degree of wetness (dryness) decreases in the warming (cooling) period, are mainly found in the windward regions facing the prevailing seasonal winds and is accounted for by the increased (decreased) PET during the southwesterly (northeasterly) monsoon.At this short time scale SPEI is able to highlight the pre-summer monsoon heating and mid-summer drought phenomena climatologically occurring in the first half of May and July, respectively.Similar contrasts are also observed at 1-month (Fig. 12) and 3-month (Fig. 13) seasonal time scales and suggest a more passive role of wintertime low temperature (thereby a lower PET) caused by the passage of dynamic rainfall processes and a more active role of summertime high temperature (thereby a higher PET) which may result in the thermodynamic rainfall process.Therefore, different temperature factor roles occur through the PET that affect the local dry/wet conditions between warm and cold months/seasons can be identified easier using the SPEI.
As an application example, the gridded SPEI-3m data in May are used to study the interannual variability in springtime drought problem in Taiwan.The EOFA identifies two dominant modes, one with an island-wide coherent dry (or wet) pattern emphasizing the southwestern/western plains (Fig. 14) and the other with a north-south orientated dry-to-wet dipole pattern (Fig. 15).The correlation pattern between the PC1 series and gridded SPI anomalies (Fig. 16a) suggests the temperature effect importance on the regional water balance over the Pingtung Plains, East Rift Valley, Taipei Valley, and Lanyang Plains.In contrast, the correlation pattern between the PC2 series and gridded SPI anomalies (Fig. 16b) suggests that the dipole pattern formation is attributed mainly to the rainfall variability and temperature effect is likely negligible.
The composite analysis results (Fig. 17) indicate that, from the regional perspective, the island-wide EOF1 drying pattern in spring is in part attributable to the interannual PDO fluctuations in the North Pacific and is in some cases associated with the ENSO turnabout in the tropical Pacific during their negative phases.Both likely coordinate to form an anomalous anticyclone over Taiwan with downslope subsidence strengthening the drought condition in the western and southwestern plains.The large-scale environment associated with the 'southern flood vs. northern drought' dipole pattern of EOF2 expresses a low-level convergent zone extending from the northern part of SCS eastward to the northern Marianas and a low-level divergence with bifurcated easterly anomalies near the east coasts of Taiwan.The dipole pattern with regional dry/wet variability is consistent with the enhanced (suppressed) rainfall activity in Vietnam and the northern Philippines (southeastern China).The depicted springtime PJ-like wave train pattern is intriguing.The physical mechanism behind it and possible connection with the summertime PJ pattern are worth further exploration (Weng and Yang, in preparation).
Fig. 1.(a) The topography of Taiwan and the geographic locations of rainfall stations operated by various governmental agencies: the automatic rainfall observing network (ARON; red dots) by the Central Weather Bureau (CWB), the hydro stations (blue dots) by the Water Resources Agency (WRA), Ministry of Economic Affairs (MEA), and synoptic stations (SYNOP) by the CWB plus airport observations by the Civil Aeronautics Administration (CAA), Ministry of Transportation and Communications (yellow dots); (b) the topography of Taiwan and specific statistical distributions: Pearson type III (PE3; orchid dots), log-normal (LN3; aqua dots), generalized logistic (GLO; silver dots), generalized extreme-value (GEV; green dots), generalized Pareto (GPA; yellow dots), Gamma (GAM; dark orange dots), Gumbel (GUM; dark orange dots), four-parameter Kappa (KAP; red dots), and five-parameter Wakeby (WAK; blue dots) distributions, employed at individual stations while calculating the SPEI on the 1-month time scale; (c) the administrative districts of Taiwan.

(
Hosking and Wallis 1997).The most appropriate distribution is thus the one with the smallest Euclidean distance d,

Fig. 2 .Fig. 3 .
Fig. 2. L-moment ratio diagrams for station (a) D (i.e., P minus PET), (b) D * , (c) SPEI, and (d) SPI time series calculated at the time scale of 1-month are shown as tiny blue points.The theoretical L-moment ratios for Pearson type III (PE3; dashed olive line), log-normal (LN3; dashed black line), generalized logistic (GLO; dashed red line), generalized extreme-value (GEV; dashed blue line), generalized Pareto (GPA; solid black line), exponential (Exp; olive point), and normal (large black point) distributions are also shown as lines and points with L skewness ( ) 3 x and L kurtosis ( ) 4 x as the abscissa and ordinate coordinates, respectively.

Fig. 4 .
Fig. 4. The configurations are the same as in Fig. 2, except for the 3-month time scale.
e., D * [ ) and observed D * series ( ) model uncertainty in D * [ series is in general small- er than the D X series.Regardless of which month and which time scale, some sporadic wet-biases (+0.3 ~ 0.5 mm day -1 ) can be found in the eastern/northeastern areas of Taiwan.In July, however, large dry-biases reaching -4.0 mm day -1 are found in the foothills of southern CMR at 3-month time scale.Similar dry-bias (-1.0 ~ -2.0 mm day -1 ) also appears in the East Rift Valley regions in October.Is a 3.0 -4.0 mm day -1 error in the mountainous areas acceptable?Mean rainfall accumulation at Alisan (the nearest ref-erence station) is 608.4 mm inJuly (1960July (  -2012 averaged) averaged), which is equivalent to 19.6 mm day -1 .The largest relative error is thus about 15.3 -20.4%, which reflects the statistical stability data scarcity problem in the mountainous areas.

Fig. 9 .
Fig. 9.The mean seasonal correlations between the SPEI and SPI time series on the 1-month time scale.Note that the correlations are subtracted from 1.0 while drawing the results.

Fig. 10 .
Fig. 10.The climatological-means spatial distributions in the SPEI at 10-day time scale in the 1 st and the middle days of each month.

Fig. 11 .
Fig. 11.The climatological-means spatial distributions in the SPI at 10-day time scale in the 1 st and the middle days of each month.

Fig. 12 .
Fig. 12.The climatological-means spatial distributions in the SPEI (left column) and SPI (right column) at 1-month time scale in January, April, July, and October (from-top-to-bottom).

Fig. 13 .
Fig. 13.The climatological-means spatial distributions in the SPEI (left column) and SPI (right column) at 3-month time scale in February, May, August, and November (from-top-to-bottom).
Fig. 16.The correlation maps between the anomalous gridded SPI series and (a) the PC1 series and (b) the PC2 series obtained from the EOFA of gridded SPEI anomalies.
Fig. 17.The composited SST spatial distribution (color shadings), 850-hPa vector wind, and stream function (contours; contour interval is 2.0 × 10 6 m 2 sec -1 ) anomalies in MAM related to (a) the first and (b) the second EOF of the gridded SPEI-3m anomalies in May.Only locally significant SST and wind vectors at 1.0% level are shown.