The Influence of Synoptic Weather Types and Moisture Transport Pathways on Precipitation Isotopes in Southern Patagonia

We analyzed 28-year-long monthly oxygen isotope composition of precipitation (δ18Op) data from Punta Arenas (Chile) on the leeward side of the Andes to understand how different synoptic weather types and moisture transport pathways influence δ18Op variability in this region. Combining weather station 6 h precipitation data and atmospheric back trajectories, we found that in such a region where the atmospheric circulation pattern is dominated by very strong westerlies, an increased monthly proportion of easterly-delivered precipitation—with the air-mass trajectory path evading the influence of Andean “isotopic rain shadow” and having less rainout en route—would increase δ18Op. These synoptic easterlies are a result of quasi-stationary blocking-like flow that are an important but underappreciated part of regional circulation patterns and climate. In addition, synoptic easterlies are more often associated with heavy precipitation events as shown by weather station data and higher deuterium excess that indicates weaker post-condensation raindrop re-evaporation. Therefore, our analysis demonstrated the process link between the frequency of synoptic weather types characterized by blocking-like flow and temporal variations in δ18Op in Southern Patagonia. We conclude that isotope proxy paleo-records in this region could provide unique insights into the behaviors and dynamics of the large-scale Southern Hemisphere Westerly Winds over long timescales.


Introduction
Stable isotopes in water (δ 18 O and δ 2 H) have been widely used as geochemical tracers for the hydrological cycle. Isotopic fractionation occurs in water molecules during phase transitions modulated by transport and also biotic processes, leading to distinct variations in isotopic compositions across the hydrosphere [1]. Of particular importance is the isotopic composition in precipitation (most of studies focus on oxygen isotopes in precipitation, δ 18 O p ) with a large number of data collected since the 1950s [2]. The δ 18 O p data showed empirical relationships with climatic variables such as air temperature in the high-latitude regions and precipitation amount in the tropical regions [3]. The Rayleigh distillation model that depicts isotopic fractionations in air mass rainout processes could explain much of δ 18 O p variability at various temporal and spatial scales in both observations and models [3][4][5]. These well-established relationships between δ 18 O p and climatic variables led to applying stable isotope proxies in terrestrial climate archives to reconstruct past hydrological  [47].
The influence of the SHWWs on Punta Arenas is year-round. During austral summer, the SHWWs form a condensed and strong belt between 45°S and 55°S and maximum wind speeds are reached in Punta Arenas. During austral winter, the SHWW belt expands northward as far as 30°S while the wind speed only relatively weakens in the south, including at Punta Arenas [48]. Despite the seasonal variations in the SHWWs, the seasonality of precipitation in Punta Arenas is very weak with only slightly more precipitation during austral autumn [49]. The mean annual air temperature in Punta Arenas airport was 6.2°C, while the mean warmest month (January) air temperature was 10.8°C and the mean coldest month (July) air temperature was 1.6°C during the period 1990-2017 based on data from GHCN Version 3 [43]. Based on the Köppen−Geiger climate classification map [50], Punta Arenas is located at the boundary between Tundra climate (ET) and Subpolar oceanic climate (Cfc) zone.
Inter-annual or long-term climate variability in Southern Patagonia has been tied to the large-scale Southern Annular Mode (SAM), which is defined as the difference in zonal mean sea level pressure (SLP) between 40° S and 65° S [51]. This pressure gradient controls hemisphere-wide strength and latitudinal position of the SHWWs, and consequently, their associated storm tracks. Instrumental and Figure 1. Digital elevation map of Southern Patagonia. Also shown are major cities (cyan triangles) and field sites (red dots) discussed in this paper, and the elevation contour of 1000 m asl. The digital elevation model is from the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) at 30-arc-second resolution [47].
The influence of the SHWWs on Punta Arenas is year-round. During austral summer, the SHWWs form a condensed and strong belt between 45 • S and 55 • S and maximum wind speeds are reached in Punta Arenas. During austral winter, the SHWW belt expands northward as far as 30 • S while the wind speed only relatively weakens in the south, including at Punta Arenas [48]. Despite the seasonal variations in the SHWWs, the seasonality of precipitation in Punta Arenas is very weak with only slightly more precipitation during austral autumn [49]. The mean annual air temperature in Punta Arenas airport was 6.2 • C, while the mean warmest month (January) air temperature was 10.8 • C and the mean coldest month (July) air temperature was 1.6 • C during the period 1990-2017 based on data Atmosphere 2020, 11, 514 4 of 26 from GHCN Version 3 [43]. Based on the Köppen−Geiger climate classification map [50], Punta Arenas is located at the boundary between Tundra climate (ET) and Subpolar oceanic climate (Cfc) zone.
Inter-annual or long-term climate variability in Southern Patagonia has been tied to the large-scale Southern Annular Mode (SAM), which is defined as the difference in zonal mean sea level pressure (SLP) between 40 • S and 65 • S [51]. This pressure gradient controls hemisphere-wide strength and latitudinal position of the SHWWs, and consequently, their associated storm tracks. Instrumental and reanalysis data suggest the leading role of the SAM on both temperature and precipitation patterns in broad regions of the extratropical Southern Hemisphere, including Southern Patagonia [48,52]. Specifically, a positive phase of the SAM, through the stronger and southwardly positioned SHWWs, increases air temperature but only marginally decreases precipitation in the domain of Punta Arenas, and vice versa [48,52,53]. The weak influence of the SAM on precipitation in Southern Patagonia is partly due to the fact that the SAM is a hemisphere-wide climate mode, while regional precipitation variability is mainly controlled by synoptic-scale processes.

GNIP Isotope Data
Monthly precipitation isotope data from the GNIP Punta Arenas (airport) station are available from 1990 to 2017 with some short data gaps (Table S1) [34]. In total, there are 298 individual δ 18 O p measurements. The dataset contains some very high or even positive δ 18 O p values after 2002 (Figure 2b), which have been pointed out previously [49]. These samples likely experienced evaporation for unknown reasons, and the reason why these very high δ 18 O p values only appeared after 2002 is unknown. Since 2010, GNIP Punta Arenas station has started to label isotope data that showed a hint of sample evaporation if the calculated deuterium excess (d-excess, defined as δ 2 H − 8 × δ 18 O) value was below zero, although some of these potential problematic data were not necessarily impacted by sample evaporation, as d-excess usually has a relatively large scatter, and a value lower than zero is not unusual in maritime climates [1] such as in Punta Arenas. For the following HYSPLIT analysis, we filtered out GNIP data by removing those δ 18 O p data that were above −3% or that had corresponding d-excess below −10% (Figure 2a). These criteria were subjectively chosen but reasonable as δ 18 O p and d-excess values falling in these ranges are indeed too extreme to occur in monthly-composited precipitation at this latitude. Another GNIP station in Ushuaia (Argentina), located~250 km southeast of Punta Arenas (see Figure 1), does not have any δ 18 O p data above −6.8% during its monitoring period 1981-2002. These two stations have a similar climate, yet Punta Arenas has a weighted mean annual δ 18 O p value 2% -3% higher than Ushuaia; the exact value depends on how the potential outliers in both datasets were treated [49]. Considering this 2% -3% offset would justify our conservative cut-off value (−3% ) for δ 18 O p . There are a few data with d-excess below −10% in Ushuaia likely due to sample evaporation as well, but these data only occurred before 1988 and were not an intermittent and common feature. This observation would also justify our cut-off value (−10% ) for d-excess. Nonetheless, we cannot rule out the small possibility that a few of these anomalously high δ 18 O p or low d-excess values are naturally occurring. If they were, the collected precipitation samples might be caused by a dominant contribution of local convective precipitation [18,33] or by a significant role of raindrop post-condensation re-evaporation [16,28]. Both mechanisms would underplay the role of moisture transport pathway in shaping the δ 18 O p value, thus removing a few of such values does not affect our analysis on robust data. We also removed two very low monthly δ 18 O p values. These two data points appear out of bound of typical observations ( Figure 2a) and likely suggest additional factors that are not considered in our HYSPLIT framework; otherwise, these two extremes would have skewed the analysis. After the step of data filtering (all filtered out data are listed in Table S1), we obtained a final dataset with 265 δ 18 O p measurements. Each δ 18 O p datum has a corresponding d-excess value calculated from its paired δ 2 H p (hydrogen isotope composition of precipitation) datum.

of 26
May; MAM) (r = 0.38, p < 0.01) and spring (September-November; SON) (r = 0.38, p < 0.01). We thus define the δ 18 which allows us to explore other drivers controlling monthly δ 18 Op variability that are unrelated to δ 18 Op seasonality and/or air temperature. This approach is similar to but not the same as the multivariate regression approach by Krklec and Domínguez-Villar [54], who included both air temperature and monthly percentages of different moisture sources into their multivariate models.  Table S1 for details). (c) Scatter plot between monthly δ 18 Op and Global Historical Climatology Network (GHCN) monthly air temperature and their linear regression line. The use of GHCN monthly air temperature data [44] rather than GNIP monthly air temperature data is due to several potential erroneous data recordings in the GNIP database (see Figure A2d for details). One data gap in GHCN (February 2011) is filled by GNIP data. (d) Scatter plot between monthly δ 18 Op and Punta Arenas airport monthly cumulative precipitation (calculated from the sum of all 6 h precipitation data for each month; see Section 3.2 for details). There are 6 individual months during which precipitation data have temporal gaps, and these data are replaced by GHCN monthly precipitation amount data (see Figure A2a). In (c) and (d), symbols are color-coded to represent different austral seasons: summer (December-February), autumn (March-May), winter (June-August), and spring (September-November).  Table S1 for details). (c) Scatter plot between monthly δ 18 O p and Global Historical Climatology Network (GHCN) monthly air temperature and their linear regression line. The use of GHCN monthly air temperature data [43] rather than GNIP monthly air temperature data is due to several potential erroneous data recordings in the GNIP database (see Figure A2d for details). One data gap in GHCN (February 2011) is filled by GNIP data. (d) Scatter plot between monthly δ 18 O p and Punta Arenas airport monthly cumulative precipitation (calculated from the sum of all 6 h precipitation data for each month; see Section 3.2 for details). There are 6 individual months during which precipitation data have temporal gaps, and these data are replaced by GHCN monthly precipitation amount data (see Figure A2a). In (c) and (d), symbols are color-coded to represent different austral seasons: summer (December-February), autumn (March-May), winter (June-August), and spring (September-November).
Among the selected data, it is clear that the monthly δ 18 O p value has a significant, though weak, linear correlation (p < 0.001) with monthly air temperature (Figure 2c), but not monthly precipitation amount (Figure 2d). A linear regression yields the relationship between monthly δ 18 O p and monthly air temperature (T) as: where data for T are from GHCN version 3 [43] (Figure 2c). The observed seasonality in δ 18 O p is a prevalent phenomenon across the mid-and high-latitude regions even in maritime climate [1]. This δ 18 O p seasonality is driven by multiple processes, including changing moisture source region and condition as well as changing transport path and condensation history of water vapor among different seasons via the Rayleigh mechanism [1,4,32]. These net effects on δ 18 Atmosphere 2020, 11, 514 6 of 26 which allows us to explore other drivers controlling monthly δ 18 O p variability that are unrelated to δ 18 O p seasonality and/or air temperature. This approach is similar to but not the same as the multivariate regression approach by Krklec and Domínguez-Villar [54], who included both air temperature and monthly percentages of different moisture sources into their multivariate models.

Weather Station Precipitation Data
Punta Arenas airport cumulative precipitation amount data at 6 h intervals (at UTC time 00:00, 06:00, 12:00, and 18:00) are provided by the Dirección Meteorológica de Chile and are available for public download on [55]. The data were collected by Punta Arenas airport weather station. This same weather station also reports monthly climatology data to the GHCN (as presented in Section 3.1) and operates isotope sample collection for the GNIP [56]. The quality of 6 h cumulative precipitation data is robust as we find the monthly cumulative precipitation amount data, calculated from the sum of all 6 h precipitation data for each month, are clustering on the 1:1 line when plotted against either reported GHCN or GNIP monthly precipitation amount data ( Figure A2a,b). Precipitation events (at least 0.1 mm) were recorded during 17% of all 6 h intervals. With 6 h cumulative precipitation data, we can explicitly know the time and precipitation amount of each precipitation event that together constitute the monthly precipitation sample that GNIP measured for isotopic composition.

Methods
We ran the HYSPLIT back-trajectory model to reconstruct the origin and path of air mass that resulted in precipitation event recorded by the weather station described above. We employed a Python-based package for HYSPLIT (PySPLIT) [57] to achieve a fast and flexible workflow in trajectory generation and analysis. Five-day (120 h) air mass backward trajectory was generated four times (at UTC time 00:00, 06:00, 12:00, and 18:00) at initial back-trajectory heights at 500, 1000, 1500, and 2000 m above ground level (agl) for each day from 1990 to 2017 to match with the 6 h cumulative precipitation data from the weather station. The duration of 120 h of back-trajectory modeling was chosen as it was long enough to allow trajectory clustering as introduced below and short enough to save computation capacity. For gridded meteorological datasets, we used the Global Data Assimilation System (GDAS) reanalysis dataset (1.0 • spatial resolution) for the period 2005-2017. For the period 1990-2004 during which the GDAS dataset was not available, we used the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis dataset (2.5 • spatial resolution) [58]. Both datasets are archived on [59]. Based on the 6 h cumulative precipitation data from the weather station, trajectory files that are not associated with precipitation events were deleted. Then, all remaining trajectory files are associated with documented precipitation events. There are four trajectories with different initial back-trajectory heights for each event. For GDAS-based trajectories, only one trajectory at the height level that has the highest initial relative humidity was chosen for each event to better capture the level of precipitating air mass. Overall, 16%, 26%, 24%, and 34% of precipitation events have the highest relative humidity at the level of 500 m, 1000 m, 1500 m, and 2000 m, respectively. For NCEP/NCAR-based trajectories, this selection criterion was not reliable, as 64% of all events have the highest relative humidity at 500 m due to the limitation of the coarser-resolution reanalysis dataset; thus, we considered all cases of different initial back-trajectory heights specifically in the following analysis. We employed this simple but still realistic representation of precipitating air masses widely employed in literature [16,27,54,60,61], while a more complex approach has been applied using remote sensing data [32] or in an idealized framework [62]. To visualize the general trajectory pathway, we generated trajectory frequency contours by binning all precipitation-producing trajectory end point data in each 1.0 • grid.
We grouped all generated trajectories into four clusters based on where trajectories entered into Southern Patagonian landmass. Although HYSPLIT GUI has provided a k-means clustering algorithm, it suggests that only a large number of clusters are appropriate to represent all trajectories in an acceptably low total spatial variance. To simplify the solution, we grouped trajectories based on the locations of back-trajectory end points when trajectories came across the eastern, northern, western, and southern boundaries of Southern Patagonian landmass ( Figure A3). These trajectories were then clustered as the distal easterlies, northerlies, westerlies, and southerlies, respectively ( Figure 3). These four trajectory pathway clusters represent major synoptic weather types in this region (see [45]). Here the word "distal" is to discern our definition of air mass transport direction that is based on modeling its back-trajectory over a course of time from transport direction defined based on "local" weather station wind direction (e.g., Figure A1 introduced in Section 2 and [14]). Using four simplified distal trajectory clusters reflects our goal that we are mainly interested in how Patagonian "upwind" terrains influence δ 18 O p in relatively inland Punta Arenas. Costal terrains on western (the Andes) and southern (the Cordillera Darwin) boundaries would produce the isotopic rain shadow effect [11,63] on their leeward side-air masses would have to experience orographic uplift, condensation, and rainout before reaching Punta Arenas ( Figure 1). By contrast, terrains on northern and eastern boundaries are relatively flat without main topographic barriers; thus, air masses would experience different degrees of fractionation [28].
Atmosphere 2020, 11, x FOR PEER REVIEW 7 of 25 Using four simplified distal trajectory clusters reflects our goal that we are mainly interested in how Patagonian "upwind" terrains influence δ 18 Op in relatively inland Punta Arenas. Costal terrains on western (the Andes) and southern (the Cordillera Darwin) boundaries would produce the isotopic rain shadow effect [11,63] on their leeward side-air masses would have to experience orographic uplift, condensation, and rainout before reaching Punta Arenas ( Figure 1). By contrast, terrains on northern and eastern boundaries are relatively flat without main topographic barriers; thus, air masses would experience different degrees of fractionation [28]. There are 13 precipitation events based on the 6 h precipitation data in this month. The total monthly precipitation amount is 26 mm. Trajectory paths that enter eastern, northern, western, and southern boundaries of Southern Patagonian landmass are grouped into clusters of distal easterlies (brown), northerlies (green), westerlies (blue), and southerlies (orange), respectively. The trajectory line width is proportional to the 6 h precipitation amount associated with each trajectory. For reference, the thickest and thinnest paths on this graph represent 4.6 mm and 0.2 mm precipitation, respectively. Yellow triangle marks Punta Arenas. The proportions of precipitation contributed by different trajectory clusters (easterly-derived-dE, northerly-derived-dN, westerlyderived-dW, and southerly-derived-dS) in that month can be calculated. This provides an additional metric that allows us to investigate if the proportion of precipitation derived from a certain trajectory cluster could influence monthly-composited δ 18 Op and d-excess in GNIP data.
Because each trajectory is corresponding to a 6 h precipitation event, then the fraction of each trajectory cluster contributing to monthly total precipitation is calculated as (using the cluster of distal easterlies as an example): where dE is the monthly percentage of distal easterly-derived precipitation, ∑ p E is the sum of all distal easterly-derived 6 h precipitation, and P is the sum of all 6 h precipitation for a particular month. Following Equation (3), the percentages of distal northerly-derived (dN), westerly-derived (dW), and southerly-derived (dS) precipitation can be also calculated ( Figure 3). As such, we could quantify the contributions of major synoptic weather types to total monthly precipitation.
To extract the influence of synoptic weather types and moisture transport pathways on δ 18 Op in Punta Arenas, we directly looked into whether the δ 18 O p e metric is significantly correlated with dE, dN, dW, and dS, using the Pearson correlation coefficient. The statistical analysis was run for four different austral seasons. Correlation is considered significant if p-value is less than 0.1. If so, it would indicate that precipitation derived from that trajectory cluster has on average a distinct high or low δ 18 Op value relative to others, leading to statistical correlation between its proportion and monthly δ 18 O p e in our framework. Here a less conservative p-value is chosen because as presented below, the Figure 3. Visualization of GDAS-based precipitation-producing back-trajectories for an example month (August 2011). There are 13 precipitation events based on the 6 h precipitation data in this month. The total monthly precipitation amount is 26 mm. Trajectory paths that enter eastern, northern, western, and southern boundaries of Southern Patagonian landmass are grouped into clusters of distal easterlies (brown), northerlies (green), westerlies (blue), and southerlies (orange), respectively. The trajectory line width is proportional to the 6 h precipitation amount associated with each trajectory. For reference, the thickest and thinnest paths on this graph represent 4.6 mm and 0.2 mm precipitation, respectively. Yellow triangle marks Punta Arenas. The proportions of precipitation contributed by different trajectory clusters (easterly-derived-dE, northerly-derived-dN, westerly-derived-dW, and southerly-derived-dS) in that month can be calculated. This provides an additional metric that allows us to investigate if the proportion of precipitation derived from a certain trajectory cluster could influence monthly-composited δ 18 O p and d-excess in GNIP data.
Because each trajectory is corresponding to a 6 h precipitation event, then the fraction of each trajectory cluster contributing to monthly total precipitation is calculated as (using the cluster of distal easterlies as an example): where dE is the monthly percentage of distal easterly-derived precipitation, p E is the sum of all distal easterly-derived 6 h precipitation, and P is the sum of all 6 h precipitation for a particular month. Following Equation (3), the percentages of distal northerly-derived (dN), westerly-derived (dW), and southerly-derived (dS) precipitation can be also calculated ( Figure 3). As such, we could quantify the contributions of major synoptic weather types to total monthly precipitation.
To extract the influence of synoptic weather types and moisture transport pathways on δ 18 O p in Punta Arenas, we directly looked into whether the δ 18 O e p metric is significantly correlated with dE, dN, dW, and dS, using the Pearson correlation coefficient. The statistical analysis was run for four different austral seasons. Correlation is considered significant if p-value is less than 0.1. If so, it would indicate that precipitation derived from that trajectory cluster has on average a distinct high or low δ 18 O p value relative to others, leading to statistical correlation between its proportion and monthly δ 18 O e p in our framework. Here a less conservative p-value is chosen because as presented below, the monthly proportions of precipitation contributed by different trajectory clusters have skewed distributions towards 100% or 0%. If skewed data were input into bivariate correlation, it often causes an "attenuation bias" [64] in the Pearson correlation coefficient, weakening its p-value. If significant correlation was found between dW and δ 18 O e p within a specific season, we calculated the regression coefficient to determine the sensitivity of δ 18 O p to the proportion of westerly-derived precipitation.
To investigate the cause for the distinct δ 18 O p value in a certain trajectory cluster, we calculated median, lower quartile, and upper quartile values for specific humidity (normalized to the initial air mass value) and precipitation rate (both are meteorological output variables in HYSPLIT) in each hour along modeled trajectories to characterize the rainout processes for different trajectory clusters. The meteorological output variables have been used to characterize air mass rainout history in some event-based studies (e.g., [16,30]). In reality, each trajectory has its own air mass history, not only in terms of trajectory path but also along-trajectory meteorological data. Thus, the aim of this analysis is to characterize the general evolutions of air masses in different trajectory clusters rather than to quantitively model the degree of rainout and isotopic fractionation in air masses. We also analyzed the correlation between d-excess data and the proportion of precipitation derived from different trajectory clusters to diagnose if moisture recycling or post-condensation raindrop re-evaporation affected precipitation isotopes for different moisture transport pathways.

Precipitation Origin from the HYSPLIT Framework
Mean monthly proportions of precipitation derived from the HYSPLIT framework (Table 1) and precipitation-producing trajectory frequency contour plots ( Figure 4) are consistent with empirical observations that the strong SHWW belt dominates over Punta Arenas all year long. The distal westerlies on average contributed 82% (GDAS-based) or at least 85% (NCEP/NCAR-based) of precipitation despite the fact that Punta Arenas is on the leeward side of the Andes under the influence of the SHWWs. In general, the dW histogram is highly skewed towards the range of 95%-100% ( Figure A4), meaning that for around half of all months the dW are higher than 95% while other trajectory cluster proportions are less than 5%. The dW have clear seasonal variations that are also consistent with the observed SHWW seasonal migration (see Section 2). The dW are higher and lower during summer and winter, respectively, while the shoulder seasons (autumn and spring) have intermediate dW ( Table 1).
The derived monthly proportions of precipitation for each trajectory cluster are sensitive to the meteorological reanalysis dataset used for modeling back-trajectories. Although back-trajectories were modeled over different periods, using the GDAS reanalysis dataset would generate lower dW than using the NCEP/NCAR reanalysis dataset, especially during winter by at least 10%, but not for summer (Table 1). In addition, using GDAS would also generate higher dE than dN, while using NCEP/NCAR would generate the opposite ( Table 1). Considering that the NCEP/NCAR reanalysis dataset has a lower resolution, some of their modeled back-trajectory paths likely would be deflected compared to the GDAS-based back-trajectory paths, leading to biased mean estimates on the origin of precipitation. Therefore, back-trajectory clustering results for periods 1990-2004 and 2005-2017 should be treated separately in the following discussion because these were derived from different meteorological reanalysis datasets. 100% ( Figure A4), meaning that for around half of all months the dW are higher than 95% while other trajectory cluster proportions are less than 5%. The dW have clear seasonal variations that are also consistent with the observed SHWW seasonal migration (see Section 2). The dW are higher and lower during summer and winter, respectively, while the shoulder seasons (autumn and spring) have intermediate dW (Table 1). The derived monthly proportions of precipitation for each trajectory cluster are sensitive to the meteorological reanalysis dataset used for modeling back-trajectories. Although back-trajectories were modeled over different periods, using the GDAS reanalysis dataset would generate lower dW than using the NCEP/NCAR reanalysis dataset, especially during winter by at least 10%, but not for summer (Table 1). In addition, using GDAS would also generate higher dE than dN, while using NCEP/NCAR would generate the opposite ( Table 1). Considering that the NCEP/NCAR reanalysis dataset has a lower resolution, some of their modeled back-trajectory paths likely would be deflected compared to the GDAS-based back-trajectory paths, leading to biased mean estimates on the origin of precipitation. Therefore, back-trajectory clustering results for periods 1990-2004 and 2005-2017 should be treated separately in the following discussion because these were derived from different meteorological reanalysis datasets.
Furthermore, the derived monthly proportions of precipitation for each trajectory cluster are also sensitive to the choice of initial back-trajectory height for the NCEP/NCAR dataset. Although the four selected initial back-trajectory heights represent the levels at which most of atmospheric moisture is entrained [65], higher initial back-trajectory heights would generate higher dW (Table 1). Furthermore, the derived monthly proportions of precipitation for each trajectory cluster are also sensitive to the choice of initial back-trajectory height for the NCEP/NCAR dataset. Although the four selected initial back-trajectory heights represent the levels at which most of atmospheric moisture is entrained [65], higher initial back-trajectory heights would generate higher dW (Table 1). These results demonstrate the role of boundary-layer wind shear on modifying back-trajectory path, although many studies have suggested that different initial back-trajectory heights would not affect the general path of back-trajectory [54,60,65,66]. The influence of initial back-trajectory height is also illustrated in the plots of precipitation-producing trajectory frequency contours that encompassed the Atlantic Ocean and Northern Patagonia underscoring distal easterlies and northerlies when back-trajectories were initiated at 500 m agl (Figure 4b), but not when back-trajectories were initiated at higher levels ( Figure 4c).

The Effect of Synoptic Weather Types on δ 18 O p
In the correlation coefficient heat map between the monthly proportions of precipitation derived from different trajectory clusters and monthly δ 18 O p (Figure 5a), the dW and dE are negatively and positively correlated with δ 18 O e p for the shoulder seasons from the GDAS reanalysis dataset at 0.1 level. From the NCEP/NCAR reanalysis dataset, the dW and dN are negatively and positively correlated with δ 18 O e p in autumn and winter while the significance level depends on the initial back-trajectory height. No correlation was found if the initial back-trajectory height was 500 m. This suggests this level might to be too low to represent the precipitating air masses, consistent with the GDAS-based initial relative humidity data (see Section 4).
The above results present evidence that an increasing proportion of distal westerly-derived precipitation decreases δ 18 O p in Punta Arenas, although this moisture transport pathway effect is not manifested for summer during which the SHWW strength reaches maximum at this latitude and causes the distal westerlies almost the sole moisture transport pathway (Table 1). These findings are consistent with the Rayleigh distillation model that predicts lower δ 18 O p on the leeward side of the Andes due to the "upwind" orographic and rainout effects of westerlies [63]. This Andean isotopic rain shadow effect may not be as strong as in other parts of Patagonia due to the lower Andean elevation at the latitude of Punta Arenas (Figure 1), but our analysis mining the monthly δ 18 O p data still demonstrated this dynamical link. Along-trajectory meteorological data in GDAS-based distal westerlies indeed show orographic rainout during air mass transport. The specific humidity starts to slightly decrease about 5 h prior to reaching Punta Arenas (Figure 6e). The precipitation rate increases and decreases about 12 h and 2 h before reaching Punta Arenas, respectively, consistent with the orographic effect as "wet windward and dry leeward" conditions ( Figure 6f). The similar process mechanism should also apply to distal southerlies, whereas in general distal southerlies contribute very little to total precipitation (Table 1) and the along-trajectory data do not show similar changes in meteorological data as distal westerlies (Figure 6g,h).
The decrease in the proportion of distal westerly-derived precipitation implies that precipitation derived from secondary transport pathways become more important, in particular, the distal easterlies and northerlies (Table 1). However, because the NCEP/NCAR and GDAS reanalysis datasets have discrepancy in resolving the trajectory paths for distal easterlies and northerlies, caution should be taken in our interpretations. For the higher-resolution GDAS dataset, an increasing proportion of distal easterly-derived precipitation would increase the δ 18 O p . Indeed, distal easterlies have a less orographic rainout during moisture transport due to the relatively flat and lower-elevation terrain in the eastern side of Patagonia (Figure 1). The along-trajectory meteorological data show that distal easterlies have a strong moisture uptake from the Atlantic Ocean with increasing specific humidity until 4 h before reaching Punta Arenas (Figure 6a), after which the specific humidity decreases and precipitation rate increases (Figure 6a,b). Precipitation mainly occurs on site rather than en route for distal easterlies, unlike distal westerlies that show clear "up-and-down" pattern in precipitation rate. As a result, distal easterly-derived precipitation attains a higher δ 18 O p signal relative to the distal westerlies. The distal northerlies, by contrast, show no appreciable change in specific humidity, and the along-trajectory rainout is very low (Figure 6c,d). In reality, the distal northerlies in our HYSPLIT framework have a longer transport length on continent (Figure 3), similar to the setting with a stronger isotopic "continental effect" that should have increased the degree of rainout and decreased δ 18 O p [3,4]. However, the flat along-trajectory specific humidity profile might indicate an efficient moisture recycling, likely from plant transpiration flux from Patagonian grasslands (rather than evaporation), which could return water vapor efficiently [67][68][69] without kinetic fractionation and lead to relatively higher δ 18 O p [61]. Nevertheless, distal northerlies also contribute little to total precipitation (Table 1) and do not play a large enough role in modifying the δ 18 O p .  Table S2. However, from the NCEP/NCAR reanalysis dataset, the dN are positively correlated with δ 18 O e p . It seems inconsistent with the GDAS-based results described above and requires some explanations. As mentioned in Section 5.1, the trajectory paths modeled from the NCEP/NCAR dataset might be deflected, particularly for distal northerlies that seemed overrepresented at the cost of distal easterlies. Instead, we examined the spatial relationships between dN and the independent European Center for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis fields to focus on the general synoptic-scale atmospheric circulation patterns that have caused high dN. We found that the NCEP/NCAR-based dN are correlated with cyclonic anomalies (negative SLP anomalies) over Northwestern Patagonia and its adjacent Pacific Ocean, and anticyclonic anomalies (positive SLP anomalies) centered over the Eastern Antarctic Peninsula (Figure 7a). Such the blocking-like flow dynamically favor the occurrence of quasi-stationary synoptic easterlies [39,40,71], resulting in an overall reduced average zonal wind speed of prevailing westerlies (Figure 7b). In addition, the NCEP/NCAR-based dE also show similar patterns in spatial correlations with SLP and zonal wind speed fields (Figure 7c,d). Therefore, although the trajectory paths for distal easterlies and northerlies modeled from the NCEP/NCAR reanalysis dataset are uncertain compared to the GDAS, their synoptic-scale features are robust. We propose that a majority of distal northerlies modeled from the NCEP/NCAR reanalysis dataset were essentially distal easterlies if modeled from the higher-resolution GDAS reanalysis dataset, which was however not available before 2005. This could explain why the dN were overestimated over the dE from the NCEP/NCAR dataset. As such, the significantly positive correlation between dN and δ 18 O e p is still consistent with our previous finding from the GDAS-based results above.
Atmosphere 2020, 11, x FOR PEER REVIEW 11 of 25 rain shadow effect may not be as strong as in other parts of Patagonia due to the lower Andean elevation at the latitude of Punta Arenas (Figure 1), but our analysis mining the monthly δ 18 Op data still demonstrated this dynamical link. Along-trajectory meteorological data in GDAS-based distal westerlies indeed show orographic rainout during air mass transport. The specific humidity starts to slightly decrease about 5 h prior to reaching Punta Arenas (Figure 6e). The precipitation rate increases and decreases about 12 h and 2 h before reaching Punta Arenas, respectively, consistent with the orographic effect as "wet windward and dry leeward" conditions ( Figure 6f). The similar process mechanism should also apply to distal southerlies, whereas in general distal southerlies contribute very little to total precipitation (Table 1) and the along-trajectory data do not show similar changes in meteorological data as distal westerlies (Figure 6g,h). The decrease in the proportion of distal westerly-derived precipitation implies that precipitation derived from secondary transport pathways become more important, in particular, the distal

The Effect of Synoptic Weather Types on Deuterium Excess
There are multiple factors that could control the d-excess in precipitation, while major studies suggest evaporative source conditions, including sea surface temperature and boundary layer relative humidity, are the drivers of precipitation d-excess on global scale and have controlled the d-excess variations preserved in polar ice cores over millennial and orbital timescales [72,73]. However, at a specific location or region on continents, the short-term (monthly and seasonal) variations in precipitation d-excess are additionally modified by a multitude of processes unrelated to oceanic conditions, including moisture recycling (evaporation and transpiration) [20,69,74] and post-condensation raindrop re-evaporation [16,20,75].
From the similar correlation coefficient heat map analysis, we found that the dW and dE are negatively and positively correlated with d-excess, respectively (Figure 5b). We speculate that the degree of raindrop post-condensation re-evaporation is the main mechanism controlling monthly d-excess variability in Punta Arenas. From the 6 h weather station data, distal westerlies have much higher frequency of precipitation amount less than 1.2 mm (for the GDAS) and 1.5 mm (for the NCEP/NCAR) in 6 h (Figure 8a,c). As a result, distal westerly-derived precipitation consists of~25% of this "drizzle" precipitation, which is more than double that of distal easterlies (Figure 8b,d). By contrast, distal easterlies have at least 50% of precipitation that has a precipitation rate of more than 5 mm within 6 h. Therefore, distal westerlies have such as higher proportion of drizzle precipitation component in which raindrop re-evaporation is significant in decreasing d-excess [28] but not for distal easterlies. Precipitation event samples collected at Río Gallegos and Laguna Potrok Aike in Eastern Patagonia (see Figure 1 for locations) [28] and other semi-arid regions such as Southeast Australia [16] also showed that drizzle precipitation is subjective to a higher degree of raindrop re-evaporation and support our speculation.
Atmosphere 2020, 11, x FOR PEER REVIEW 13 of 25 explanations. As mentioned in Section 5.1, the trajectory paths modeled from the NCEP/NCAR dataset might be deflected, particularly for distal northerlies that seemed overrepresented at the cost of distal easterlies. Instead, we examined the spatial relationships between dN and the independent European Center for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis fields to focus on the general synoptic-scale atmospheric circulation patterns that have caused high dN. We found that the NCEP/NCAR-based dN are correlated with cyclonic anomalies (negative SLP anomalies) over Northwestern Patagonia and its adjacent Pacific Ocean, and anticyclonic anomalies (positive SLP anomalies) centered over the Eastern Antarctic Peninsula (Figure 7a). Such the blocking-like flow dynamically favor the occurrence of quasi-stationary synoptic easterlies [39,40,71], resulting in an overall reduced average zonal wind speed of prevailing westerlies (Figure 7b). In addition, the NCEP/NCAR-based dE also show similar patterns in spatial correlations with SLP and zonal wind speed fields (Figures 7c and d). Therefore, although the trajectory paths for distal easterlies and northerlies modeled from the NCEP/NCAR reanalysis dataset are uncertain compared to the GDAS, their synoptic-scale features are robust. We propose that a majority of distal northerlies modeled from the NCEP/NCAR reanalysis dataset were essentially distal easterlies if modeled from the higherresolution GDAS reanalysis dataset, which was however not available before 2005. This could explain why the dN were overestimated over the dE from the NCEP/NCAR dataset. As such, the significantly positive correlation between dN and δ 18 O p e is still consistent with our previous finding from the GDAS-based results above.

Implications for SHWW Reconstruction
The response of δ 18 O e p to dW in Punta Arenas demonstrated in our analysis indicates that δ 18 O p variations preserved in terrestrial climate archives on the leeward side of the Andes have the potential to reconstruct the long-term SHWW dynamics. When the SHWWs were strengthened at this latitude, regional precipitation was more governed by the prevailing westerlies and their associated storm tracks [36,76]. When the SHWWs were weakened with an overall reduced SLP gradient over the Drake Passage similar to the regional negative phase of the SAM, other synoptic weather types such as easterly flows that are a result of quasi-stationary blocking of the prevailing westerlies [40,71] might become more frequent in regional climate. A progressively positive shift in the SAM and an overall strengthening and poleward shift in the SHWWs during recent decades have been observed and was driven by ozone depletion and CO 2 concentration rise [51,77]. The current dominance of the SHWWs is, however, not always the case in the past. For instance, temperature proxy data have shown that the SAM was in an extreme negative phase during the fifteenth century, meaning that the SHWW strength was much reduced during that time [77]. Therefore, the increasing frequency of synoptic weather types in which moisture is transported by non-westerlies during the waning of SHWWs is important for understanding the connection between δ 18 O p and SHWW behaviors in the past (Figure 9). Furthermore, the frequency and strength of blocking-like episodes in Patagonia have been linked with the El Niño-Southern Oscillation (ENSO) and would increase during the El Niño phase through the Rossby wave train [39,71]. We propose that δ 18 O p variations reflect the synoptic-scale processes related to the interactions between ENSO and SHWWs/SAM more than the simplistic strength of the SHWWs (Figure 9). Long-term temperature changes could also affect δ 18 O p in this region, in particular over orbital timescales (e.g., before the Holocene), but their roles are less prominent on multi-decadal and centennial timescales over which temperature variations are not likely to cause large shifts in δ 18 O p that if observed would instead imply a dynamical driver in the circulation pattern [39]. In addition, temperature effect might be also small during the Holocene. Synthesis data and models show that temperature variations were less than 0.8 • C over the past 8000 years [78].
Sensitivity analysis indicates that a 20% less westerly-derived precipitation could cause a 0.6-1.3% higher δ 18 O p in Punta Arenas (Table 2). A previous study close to the Atlantic coast (Río Gallegos and Laguna Potrok Aike; see Figure 1) showed that precipitation event samples associated with westerlies have δ 18 O p signals~7% higher than easterlies and northerlies [28], which from a simple isotopic and mass balance calculation could be converted to an effect of 1.4% in δ 18 O p for their sites was there 20% less westerly-derived precipitation. The different sensitivities may reflect the fact that Punta Arenas is in a central locus of Patagonian landmass,~170 km away from the Atlantic coast, thus easterlies experience a higher degree of rainout during moisture transport compared to the sites in the Atlantic coast. In addition, the Andean divide at the latitude of Punta Arenas is not as high as other parts of Patagonia. This may further reduce the "upwind" orographic and rainout-related fractionation that has occurred in air masses delivered by westerlies. Therefore, terrestrial climate archives such as lakes and peat bogs closer to the Atlantic Ocean would be more sensitive in detecting the long-term variability in the frequency of synoptic easterlies and the SHWW dynamics from the δ 18 O p paleo-records. For example, peat moss cellulose data from Ariel Peatland (see Figure 1), which was close to the Atlantic coast (~80 km), showed that the δ 18 O p was~2% higher than during the early twentieth century than the current period [39]. Taking the sensitivity into consideration, this shift might indicate that westerly-derived precipitation was~30% less during the early twentieth century compared to today as a recent analogue for the period of weakened SHWWs. Table 2. Effect of a 20% less westerly-derived precipitation on δ 18 O p based on regression coefficients between δ 18 O e p and dW if within a specific season they are statistically significant from Figure 5a. Results are compared with another study [28] in which the authors based on event-based precipitation samples concluded that precipitation derived from westerlies are on average 7% lower than easterlies and northerlies at Río Gallegos and Laguna Potrok Aike (Figure 1). From a simple isotopic and mass balance calculation, it could be converted that a 20% less westerly-derived precipitation would have an effect of +1.4% in δ 18 O p for their sites. The SHWW dynamics also regulate hydrological regimes in Southern Patagonia. Although today, the SHWWs are delivering a large amount of precipitation even to the leeward side of Southern Patagonia, they are not as efficient as other wind directions (Figures 8 and A1). Some historic extreme precipitation events along the Atlantic coast of Patagonia have been linked with enhanced moisture flux derived from easterlies [40]. Our HYSPLIT back-trajectory modeling also showed that the extreme precipitation event on 11-12 March 2012, which caused historic flooding of river Río Las Minas in the main city of Punta Arenas [79], occurred under bursts of easterlies (Figure 10a). This event was caused by a transient migratory weather system with an anomalous anticyclonic center over the Antarctic Peninsula resulting in a blocking-like flow (Figure 10b). For this month, the dE is as high as 76%. The monthly precipitation amount is 63 mm (102 mm from conventional weather station Jorge C. Schythe in the main city of Punta Arenas [79]), which is around twice as high as the monthly average ( Figure 2d). The GNIP δ 18 O p value for this month is −4.75% , which corresponds to a +3.4% of δ 18 O e p . The d-excess for this month is 1.6% , suggesting that the high δ 18 O p value is not an artifact from sample evaporation. This month is a good example showing that changes in synoptic-scale patterns in the past were able to influence both δ 18 O p and hydroclimate. If the proportion of precipitation derived from easterlies increased, it might also indicate an increase in moisture availability for grasslands and desserts in Eastern Patagonia [46]. Therefore, from terrestrial climate archives, the δ 18 O p paleo-records can be tied to other hydrological proxies to understand the coupled evolution of the synoptic-scale patterns and hydrological regimes in Southern Patagonia over long timescales (Figure 9).   [39,71]. This Rossby wave train is responsible for the apparent ENSO-SAM phase teleconnection ( 2 ) [77,80], as the circulation anomalies of blocking-like flow is similar to the negative phase of the SAM (reduced mean SLP gradient across the Drake Passage) ( 3 ). The increased occurrence of quasi-stationary blocking-like flow would increase the occurrence of synoptic easterlies in Southern Patagonia ( 4 ). As a result, the overall strength of the prevailing westerlies would be reduced ( 5 ), although the strength of the hemisphere-wide SHWWs is directly related to the SLP gradient between mid-and high-latitude characterized by the phase of the SAM ( 6 ). As the synoptic easterlies became more frequent, their representation in regional precipitation in Eastern Patagonia would he higher, increasing the mean δ 18 O p over long timescales ( 7 ). In addition, the synoptic easterlies are more often associated with larger precipitation events that would increase the overall precipitation amount and weaken the Andean rain shadow in Eastern Patagonia ( 8 ) [46].

Implications for Paleoaltimetry Reconstruction
Our results also have implications for isotope-based paleoaltimetry reconstruction in Southern Patagonia and, by extension, Northern Patagonia and other regions with similar settings such as the South Island in New Zealand (leeward side with a single dominant moisture transport pathway) [11,81]. Traditionally, these studies have assumed that a single dominant moisture transport pathway such as the strong SHWWs in Patagonia has also persisted throughout geological timescales. Under this assumption, any negative temporal shift in δ 18 O p in the modern leeward side could be interpreted as the development of isotopic rain shadow and used to estimate the surface elevation in upwind locations using an empirical isotopic lapse rate [12,37].
This study, however, shows that an additional complexity is involved in properly interpreting the signal of surface uplift from leeward isotope data. Specifically, the shifts in δ 18 O p can be alternatively caused by changes in moisture transport pathways under different atmospheric circulation patterns even in a region where today there seems to only have a single dominant moisture transport pathway. As such, shifts in Antarctic sea ice extent, Southern Hemisphere meridional sea surface temperature gradient, and tropical ocean-atmosphere system can also impose changes in SHWW strength and position [82][83][84] and, consequently, the synoptic-scale patterns and δ 18 O p signals preserved in terrestrial climate archives in the leeward side of Patagonia. Therefore, any leeward isotope-based paleoaltimetry studies should consider how atmospheric circulation was operated in the past under different paleogeography and climate boundary conditions. Because our sensitivity analysis (Table 2) indicated that δ 18 O p at sites farther away from the Atlantic coast is less sensitive to the SHWW variability, we would expect that sites from the nearest leeward locations are most suitable for detecting the major uplift signal of the Patagonian Andes from the δ 18 O p proxy data. Although the effect of changing moisture transport pathway on δ 18 O p might be not as large as the effect of mountain uplift, it could result in a biased estimate on uplift history if proxy data are collected from sites close to the Atlantic coast where a larger fraction of moisture could be advected from easterlies that break the rule of leeward paleoaltimetry.

Conclusions
Other than air temperature and precipitation amount, it has been well recognized that moisture source and transport pathway could also influence δ 18 O p for individual precipitation event on a regional scale. However, only a few paleoclimate studies are able to link δ 18 O p variations over time with synoptic-scale shifts in atmospheric circulation patterns. We here provide a new framework for extracting the relationship between the synoptic conditions and δ 18 O p hidden in decades-long GNIP monthly δ 18 O p record at a specific station. In such framework, weather station precipitation data and back-trajectory modeling are combined to quantify the proportions of precipitation delivered from different trajectory clusters that are characterized by different synoptic weather types on monthly basis. We apply this framework to GNIP monthly δ 18 O p record in Punta Arenas located in Southern Patagonia that today is under year-round influence of the strong SHWWs. We show that distal westerlies indeed contribute to a majority of precipitation in Punta Arena even though it is located in the Andean rain shadow, but not all of the precipitation. Instead, there is a considerable amount of precipitation originated from the distal easterlies, particularly during non-summer seasons. We further show that a decreasing proportion of distal westerly-derived precipitation would increase δ 18 O p , caused by a concurrently increasing proportion of distal easterly-derived precipitation that has a relatively higher δ 18 O p signal. We then link the proportions of precipitation origin with synoptic-scale atmospheric circulation patterns to provide an observational basis and a process-based mechanism for reconstructing the long-term SHWW dynamics from δ 18 O p paleo-records preserved in terrestrial climate archives. Our analysis highlights the importance of site location in preserving the most sensitive SHWW record inferred from δ 18 O p data, as sensitivity analysis indicates that the Atlantic coast might have the largest contrast in δ 18 O p between the westerly and easterly flows. Our results have implications for isotope-based leeward paleoaltimetry not only in Patagonia but also for other similar regions. Leeward paleoaltimetry studies often assume a single moisture transport pathway or source region has persisted throughout geological timescales, while our study suggests that different atmospheric circulation patterns in the past were also able to cause shifts in δ 18 O p , which sometimes were, however, interpreted as the sign of surface uplift or evaporative enrichment in soil water.
Finally, the framework introduced in this study provides new opportunities in linking large dataset of GNIP monthly precipitation isotopes with synoptic-scale patterns of moisture transport. Such framework utilizes station-based hourly or daily precipitation data to quantify the contributions of different synoptic weather types on each monthly-composited δ 18 O p as well as d-excess datum. These monthly-composited δ 18 O p and d-excess data readily available from multiple sites and over a long monitoring period, if put in a reasonable context of moisture transport, could provide insights into hydrological processes on terrestrial landscape, with implications for partitioning terrestrial water fluxes and reconstructing atmospheric circulation patterns using isotope data.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/11/5/514/s1, Table S1: Summary of all missing data or filtered dataset from GNIP Punta Arenas, Table S2: The number of δ 18 O p and d-excess data that are used in correlation heat map analysis.
Author Contributions: Conceptualization, Z.X. and Z.Y.; methodology, Z.X.; software, Z.X.; validation, Z.X.; formal analysis, Z.X.; investigation, Z.X.; resources, Z.X. and N.B.; data curation, Z.X. and N.B.; writing-original draft preparation, Z.X.; writing-review and editing, N.B. and Z.Y.; visualization, Z.X.; supervision, Z.Y.; project administration, Z.Y.; funding acquisition, Z.Y. All authors have read and agreed to the published version of the manuscript.  Six-hour cumulative precipitation amount (n = 40,769) was divided by 6 to match with each hourly wind direction data prior to each precipitation measurement, assuming that wind in each of that 6 h contributes to that 6 h cumulative precipitation equally. There are 140 observational gaps (0.3% of all hours). (c) The quotient between precipitation percentage and time percentage for each respective wind direction.
The weather station data are available for public download on [55].  data have temporal  gaps (27%, 16%, 22%, 14%, 11%, and 6% of all 6 h intervals, respectively), which are responsible for their monthly cumulative precipitations amount underestimated compared to GHCN monthly precipitation amount data. It is still likely that a few GHCN monthly precipitation amount data are erroneous as well, which might explain the deviations from the monthly cumulative precipitation amount. (b) Scatter plot between reported GNIP monthly precipitation amount and Punta Arenas monthly cumulative precipitation amount similar as (a). (c) Scatter plot between GHCN Version 2 monthly precipitation amount and GNIP monthly precipitation amount. Clearly there are also several important erroneous data recordings in GNIP database. (d) Scatter plot between GHCN Version 3 monthly air temperature and GNIP monthly air temperature. Most data are plotted along the 1:1 line, but there are still several important erroneous GNIP data recordings. These data are mainly found from January 2015 to March 2016 (shown as open dots). These GNIP monthly air temperature values are much higher than normal, while weather station data from nearby city do not show any corresponding high air temperature, confirming that these GNIP data are erroneous. Figure A3. Map showing the delineations of eastern (brown), northern (green), western (blue), and southern (orange) moisture transport boundaries of Southern Patagonia landmass used in this study. If a back-trajectory came across one of the above boundaries, that trajectory was clustered into distal easterlies, northerlies, westerlies, and southerly trajectories, respectively. Here the western boundary of Southern Patagonia landmass is set on the east of the Andes so that any back-trajectory crossing over the Andes will be classified as the distal westerlies.