Model Based Spatial Distribution of Oxygen-18 Isotopes in Precipitation Across Canada

Abstract

Challenges inherent in mesoscale (large domain) hydrologic modelling of remote ungauged basins include validation of model results and quantification of uncertainty in the predictions. Isotope equipped hydrological models, such as isoWATFLOOD, have the ability to simulate both quantity and isotopic composition of streamflow and runoff generation processes providing more options for model validation, but first require information about isotopic composition of precipitation across the model domain. The stable water isotopic composition of precipitation (18OPPT) is available from stations sampling monthly composite precipitation, but the data are spatially and temporally discontinuous, particularly in northern areas of Canada. Here we use new data from the Canadian Network for Isotopes in Precipitation (CNIP) to create and evaluate a variety of empirical relationships to develop improved interpolations of the spatial distribution of 18OPPT across Canada. The goal of this research is to develop 18OPPT prediction models that can be incorporated directly within the isoWATFLOOD hydrological modelling framework to provide spatially variant 18OPPT patterns as forcing data for the iso-hydrological model. Comparison of model results has identified models capable of simulating annual 18OPPT distributions, but also identified seasons and areas where the geographical and climatological parameters included in this analysis were not able to simulate the measured distributions. Spring, summer and fall model results were satisfactory; however, winter model results were more variable, indicating increased complexity in the driving forces of 18OPPT patterns during this season. Overall, model results improve with the addition of time-variant climate parameters, this finding being especially significant during the winter. Improving the precipitation input fields within isotope-equipped hydrological models will provide a valuable tool for water use management within large, remote, and often ungauged Canadian rivers and will facilitate studies of both climate variability and surface hydrology in remote regions.

Les dfis inhrents la modlisation hydrologique chelle moyenne (mso-chelle) des bassins non jaugs loigns englobent la validation des rsultats de modle et la quantification de l'incertitude lie aux prdictions. Les modles hydrologiques et isotopiques, par exemple le modle isoWATFLOOD, ont la capacit de simuler la fois la quantit et la composition isotopique des processus de gnration du dbit et de l'coulement, ce qui offre davantage d'options pour la validation de modle. Cependant, ils exigent d'abord des donnes l'gard de la composition isotopique des prcipitations dans l'ensemble du domaine modlis. La composition en isotopes stables des prcipitations (d18OPPT) est disponible partir de stations procdant l'chantillonnage de composites mensuelles. Toutefois, les donnes sont discontinues sur les plans spatial et temporel, en particulier dans les rgions nordiques du Canada. Ici nous utilisons de nouvelles donnes tires du CNIP, c'est--dire le Canadian Network for Isotopes in Precipitation (Rseau canadien pour les isotopes dans les prcipitations), afin de crer et d'valuer diverses relations empiriques dans le but de concevoir de meilleures interpolations de la distribution spatiale des d18OPPT l'chelle du Canada. L'objectif de la recherche consiste laborer des modles de prdiction des d18OPPT pouvant tre intgrs directement au cadre de modlisation hydrologique isoWATFLOOD afin de fournir des tendances des d18OPPT spatialement variables en tant que donnes de forage pour le modle d'observation hydro-isotopique. La comparaison des rsultats de modle a permis d'identifier des modles capables de simuler les distributions annuelles de d18OPPT. Cependant, elle a galement permis de cerner des saisons et des rgions o les paramtres gographiques et climatologiques inclus dans l'analyse ne permettaient pas de simuler les distributions mesures. Les rsultats de modle au printemps, en t et l'automne se sont avrs satisfaisants. Toutefois, les rsultats de modle obtenus en hiver taient plus variables, ce qui indique une complexit accrue des facteurs dterminants des tendances des d18OPPT au cours de cette saison. Dans l'ensemble, les rsultats de modle s'amliorent avec l'ajout de paramtres du climat variables dans le temps, cette constatation tant d'une importance toute particulire au cours de l'hiver. L'amlioration des champs de saisie des donnes de prcipitations dans les modles d'hydrologie isotopique fournira un outil prcieux pour la gestion de l'utilisation de l'eau des rivires canadiennes de grande importance, loignes et souvent non jauges et elle facilitera en outre les tudes de la variabilit climatique et de l'hydrologie de surface dans les rgions loignes.

Introduction

Hydrological models are critical tools for the sustainable management and future developments of water resources in Canada. However, problems arise in calibrating these models over the large, often remote basins for which these predictions are required, primarily due to a lack of available data at sufficient resolutions. New tools are being sought that are capable of providing reliable, large-scale predictions of current and future water availability. One such tool to improve calibration and constrain uncertainty within hydrological models is the simulation of stable water isotopes in conjunction with streamflow. This technique allows for a second level of validation for models, whereby both streamflow and isotope simulations are used to constrain model simulations. The isoWATFLOOD model (Stadnyk-Falcone, 2008), an isotopically-enabled version of the WATFLOOD hydrological model, has created the opportunity to simulate both quantity and isotopic composition of streamflow and runoff generation processes. Accurate simulation of isotopic variability in streamflow requires the continuous definition of isotopic forcing from climate, which can be provided from CNIP or GNIP (Canadian or Global Network for Isotopes in Precipitation) stable water isotopes in precipitation (¹⁸O_PPT) records. Similar to precipitation forcing, such station values of isotopes in precipitation require distribution across the watershed domain. The need for estimates of distributed ¹⁸O_PPT data where measurements are not available is required not only for hydrological modelling, but also for regional hydrological, ecological and climatic studies (Gibson et al., 2005; Stadnyk-Falcone, 2008; Birks and Gibson, 2009). The main objective of this study is to develop empirical relationships to predict average annual and seasonal average ¹⁸O_PPT patterns spatially across Canada. In evaluating which model parameters provide the best prediction of annual and seasonal isotope compositions of precipitation we will gain insight into the geographical and climatological controls on isotopic labeling of meteoric water. The long-term goal of developing this methodology is to implement the ¹⁸O_PPT prediction model within the isoWATFLOOD hydrological modelling framework to provide spatially variant ¹⁸O_PPT patterns as forcing data for the iso-hydrological model.

The abundance of oxygen-18 in precipitation is temporally and spatially variable, especially in seasonal climates similar to those over most of Canada. Changes in the oxygen-isotope signature of an air mass can be complex, and are the result of many different variables. ¹⁸O_PPT is controlled by Rayleigh distillation of atmospheric vapour, which is predominantly determined by changes in air-mass temperature due to the strong temperature-dependence of the moisture carrying capacity of air (Dansgaard, 1964; Yurtserver and Gat, 1981; Rozanski et al., 1993). In mid- to high-latitudes, this relationship is described by the depletion of heavy isotopes from water vapour as the air mass moves towards higher altitudes (altitude effect), inland from the coast (continental effect), or the transport of vapour from lower to higher latitude locations (latitude effect). Numerous geographic and climatic parameters are correlated with changes in temperature and hence rain-out, which can include latitude, elevation, distance inland from the coast and humidity. Changes in moisture sources (e.g. in Canada: Pacific, Arctic, Gulf of Mexico and/or North Atlantic, and potentially large lake systems) and transport processes, rainout history and seasonality (Lawrence and White, 1991; Gat et al., 1994; Simpkins, 1995; Kohn and Welker, 2005; Moran et al., 2007; Birks and Edwards, 2009) will also influence the stable-isotope composition of precipitation, further complicating ¹⁸O_PPT prediction. Therefore, interpolation methods with the ability to accurately predict changes in ¹⁸O_PPT via correlation with easily obtainable proxy variables are a valuable tool for extending the station based measurements of ¹⁸O_PPT across the country. Bowen and Revenaugh (2003) suggest that these informed interpolation methods are more likely to provide an improved spatial estimation of the ¹⁸O_PPT pattern over a watershed than currently used conventional methods that do not incorporate any geographical or climatological parameters. Informed interpolation ¹⁸O_PPT fields have been conducted on both global (Birks et al., 2002; Bowen and Wilkinson, 2002; Bowen and Revenaugh, 2003; Van der Veer et al., 2009) and regional (Dutton et al., 2005; Liebminger et al., 2006; Liu et al., 2008; Lykoudis et al., 2009; Zhao et al., 2011) scales. The empirical relationships developed from regional studies provided better approximations of the spatial and temporal distributions of ¹⁸O_PPT, likely due to the use of regional isotopic data, the inclusion of climate variables, and the refinement of the scale at which the ¹⁸O_PPT relationship is constructed. The poor performance of the Bowen and Wilkinson model over north central Canada (Bowen and Wilkinson, 2002) indicates that a more regional model that better represents the controls on isotopic labeling of precipitation is needed for this region. This study builds on the previous global efforts by developing multivariate linear and stepwise regression models from Canadian data and incorporating both geographic and climate variables within the model parameterizations. It represents the first informed interpolation of the spatial distribution of ¹⁸O_PPT for Canada. Due to the intended use of the model, publically available data and simplistic parameterizations are utilized to ensure seamless implementation into isoWATFLOOD.

Data Sources

The ¹⁸O_PPT data utilized in this study is provided from a new CNIP database combined with the existing GNIP database (International Atomic Energy Agency/World Meteorological Association, 2006). Together, these databases include monthly compositions for 35 Canadian stations, with the majority of records spanning at least several years, or more than 25 monthly compositions. Previous depictions of precipitation isotope fields across Canada (Bowen and Wilkinson, 2002) were based on the limited number of Canadian stations included in the GNIP database, often composed of data from discontinuous sampling campaigns. New data from the CNIP network initiated in 1997 consists of 17 stations distributed across Canada, collecting weighted monthly precipitation samples ( Figure 1 ). This network marks the first time that stations located in the southern and northern regions of the country have been simultaneously sampled, and includes 10 new stations at a longer time-series than previously available. Data are selected based on record length and completeness so that only stations with at least one complete annual cycle of isotope compositions are used. The final dataset includes 2908 monthly composite samples of ¹⁸O_PPT from 33 stations ( Table 1 ). Records with lengths less than five to ten years may introduce bias into the models, by sampling time-periods not representative of long-term averages for a given location. These biases may include not capturing both phases of the El Nino/La Nina Southern Oscillation (ENSO), the Arctic Oscillation (AO) and climate change, and should be kept in mind during interpretation of the model results.

Figure 1. Spatial Distribution of CNIP/GNIP Stations across Canada. Smaller circles indicate locations where the CNIP station is located at an elevation lower than 200m while larger circles indicate stations located at elevations higher than 200m.

Table 1.  ¹⁸O_PPT sampling record and sampling site geographic information.

CNIP Site City	Province	Latitude ()	Longitude ()	Altitude (m)	Record Length	^18OPPT ()	n
Esther	Alberta	51.67	110.2	707	19972003	16.3	73
Whitehorse	Yukon Territory	60.72	135.07	702	19611989	21.3	95
Edmonton	Alberta	53.57	113.52	671	19611966	17.1	51
Wynyard	Saskatchewan	51.77	104.2	561	19651982	15.6	77
Mayo	Yukon Territory	63.62	134.13	504	19851989	20.8	43
Saskatoon	Saskatchewan	52.1	106.43	504	19972001	14.5	35
Ste. Agathe	Quebec	46.05	74.28	395	19751982	12.6	82
Atikokan	Ontario	48.75	91.62	393	19751982	12.7	76
Chapais	Quebec	49.82	74.95	382	19972007	14.1	122
Experimental Lakes Area	Ontario	49.67	93.72	369	19972007	12.3	123
Bonner Lake	Ontario	49.38	82.12	245	19972007	13.8	121
Snare Rapids	Northwest Territories	63.52	116	241	19972007	22.4	112
Simcoe	Ontario	42.85	80.27	240	19751982	9.4	81
Egbert	Ontario	44.23	79.77	224	19972002	10.4	65
Gimli	Manitoba	50.62	96.98	223	19751982	14.2	76
Yellowknife	Northwest Territories	62.28	114.27	205	19891993	20.5	58
Fort Smith	Northwest Territories	60.02	111.97	203	19611965	18.2	40
Bay D'Espoir	Newfoundland	47.98	55.82	190	19972007	9.2	121
Saturna Island	British Columbia	48.78	123.13	178	19972007	10.5	122
Ottawa	Ontario	45.32	75.67	114	19792007	10.8	339
Resolute Bay	Nunavut	74.72	94.98	67	19891993	24.1	59
Alert	Nunavut	82.2	62.3	62	19891993	28.9	60
Goose Bay	Newfoundland	53.32	60.42	46	19972007	14.5	122
Harcourt	New Brunswick	46.5	65.27	40	19972007	10.4	123
Truro	Nova Scotia	45.37	63.27	40	19751983	9.3	92
Churchill	Manitoba	58.45	94	28	19891993	17.2	56
Cambridge Bay	Nunavut	69.06	105.07	23	19891993	23.3	55
Victoria	British Columbia	48.65	123.43	20	19751982	10.3	72
Mould Bay	Northwest Territories	76.14	119.2	15	19891993	24.9	58
Eureka	Nunavut	80	85.56	10	19891993	26.2	55
Hall Beach	Nunavut	68.47	81.15	8	19891993	21.7	55
The Pas	Manitoba	53.97	101.10	272	19751982	16.6	73
Calgary	Alberta	51.02	114.02	1049	19922001	17.9	116

Other similar studies developing spatially distributed ¹⁸O_PPT fields have had significantly more measured data available for model construction (e.g. Bowen and Wilkinson (2002), 232 stations; Dutton et al. (2005), 79 stations; and Liu et al. (2008), 55 stations). One of the limitations of this analysis is the sparse station ¹⁸O_PPT data available for Canada. Station coverage is denser along the southern portion of Canada than in the north ( Figure 1 ). The vast size, remoteness and harsh climate of northern Canada make the collection of high resolution field measurements a challenge for most meteorological and hydrological monitoring. The lack of station data in the area is a limitation to the accuracy of interpolations of ¹⁸O_PPT, but also emphasizes the necessity for the development of empirical models to estimate where station based data are not available. As more data are collected in future years, they can be integrated into the proposed model to improve Canadian ¹⁸O_PPT interpolations.

Station annual and seasonal ¹⁸O_PPT averages are constructed from monthly composite samples over each stations period of record. Four seasons are selected to represent the Canadian annum: winter (December, January, February; or DJF), spring (March, April, May; or MAM), summer (June, July, August; or JJA) and fall (September, October, November; or SON). Monthly precipitation amounts and average temperatures are available for each CNIP/GNIP isotope record and are used to amount weight annual and seasonal average ¹⁸O_PPT according to Equation 1:

where i is the month index and P represents precipitation amount. All ¹⁸O_PPT values are reported in per mille () and utilize the notation relative to the VSMOW (Vienna Standard Mean Ocean Water) standard (Craig, 1961):

For model regressions, annual and seasonal climate normals for total monthly precipitation amounts and mean monthly temperatures are used and obtained from Canada's Historical Weather and Climate Data archive (http://climate.weatheroffice.gc.ca/).

Methods

Five model types are evaluated in this study, each representing different regression techniques and parameterizations. The first model type is based on the Bowen and Wilkinson (2002) methodology (herein referred to as the BW model), that makes use of a two-step regression technique incorporating geographic parameters resulting in the following equation:

Where: LAT represents latitude (degrees), ALT represents altitude (meters) and A through D are regression coefficients.

Model types two through four (herein referred to as SIM1, SIM2 and SIM3 models, respectively) are created using simple multivariate regression, each with a slightly differing parameterization (displayed in Equations 4, 5 and 6, respectively):

where: LONG represents longitude (degrees), and PRECIP and TEMP are the long-term climate normal total monthly precipitation (mm) and mean monthly temperature (C), respectively. Similarly, A through F are regression coefficients. SIM1 shares the same parameterization as the BW model; however, it is derived utilizing simple multivariate regression rather than the two-step regression technique, therefore yielding slightly different results. Given the dependence of ¹⁸O_PPT on continentality, moisture sources, and pathways, longitude is included as a proxy for these effects. Time variant climate parameters are included in two of the model parameterizations to provide physically-based insight into the controls on seasonal variability.

The fifth and final model employs a stepwise regression procedure for model parameter selection. Stepwise regression is carried out on the parameters included in Equation 6 to create a series of stepwise models for the analysis (herein referred to as the SW model). The stepwise regression procedure is included in this study to determine which parameters are statistically significant both annually and seasonally. This method provides an automatic procedure for selecting which variables should be entered into the regression models by conducting a series of t-tests at the 95% confidence level. For each of the five model types, an annual model and four seasonal models are created, resulting in 25 models for comparison and analysis.

Model Results

Model performance is evaluated statistically through the adjusted coefficient of determination (adj. R²) and the range and variance of the model residuals. Spatial and temporal trends in simulated ¹⁸O_PPT and model residuals are examined and discussed to gauge the various strengths and weaknesses of the individual models. For the simple and stepwise multivariate regression models, the underlying assumptions behind the regression framework are tested to ensure that violation of these assumptions does not occur. Two of the main assumptions in using a simple multivariate regression model are the normality of residuals and small amounts of multicollinearity between parameters. Large amounts of multicollinearity result in an uncertain response surface, and consequently the model parameters may no longer be given any meaningful interpretation. The normality of the model residuals are evaluated through a Lilliefors test (Lilliefors, 1967) for normality at the 95% confidence level. All annual and seasonal models derived using simple multivariate regression passed the Lilliefors test at the 0.05 alpha-level, indicating normality of residuals. The presence of multicollinearity is determined by the calculation of variance inflation factors (VIF). The VIF is describes how much the variance of an estimated regression coefficient is increased due to the presence of collinearity. VIF's are calculated for each parameter within each model. Table 2 indicates the results of the VIF calculations, denoting a yes (Y) if any of the parameters within a model have a VIF greater than 10. If the VIF is greater than 10, high collinearity is present and action should be taken to reduce the collinearity and caution should be used when placing physical meaning on the regression coefficients.

Table 2. Summary of variance inflation factor results for simple linear regression (SIM) and stepwise regression (SW) models. A VIF greater than 10 is indicated with a (Y) while a VIF less than 10 is indicated by a (N).

	Models with VIF Exceeding 10
Season	SIM1	SIM2	SIM3	SW
Annual	Y	N	Y	N
DJF	Y	N	N	N
MAM	Y	N	Y	N
JJA	Y	N	N	N
SON	Y	N	Y	Y

VIF's indicate SIM1 shows high amounts of multicollinearity. This is expected, and is due to the inclusion of both LAT and LAT² within the regression model. Between the annual and seasonal SIM3 and the SW models, four instances of high multicollinearity occur, with VIF's ranging between 15 and 30. These high VIF's are present in the annual, MAM and SON models and in all instances relate to the LAT parameter, indicating that latitude and temperature are highly correlated during the spring and fall seasons. No other predictor variables show VIF's greater than 10, indicating no additional instances of multicollinearity.

Tables 3 and 4 summarize the model parameterizations and performance for the annual and seasonal models, respectively. Stepwise regression results are indicated with an (R), denoting the parameters retained within each model type. Model results are analyzed and interpreted within the Discussion Section.

Table 3. Annual model parameterization and summary statistics. Bowen and Wilkinson and simple regression parameterizations indicated by an (X) Stepwise regression retained parameters indicated by an (R).

		Parameterization
Season	Model Type	LAT	LAT^2	ALT	LONG	PRECIP	TEMP	Residual Range ()	Adjusted R^2	Residual Variance ()
Annual	BW	X	X	X				1.8 to 2.5	0.97	0.9
	SIM1	X	X	X				1.8 to 2.5	0.97	1.0
	SIM2	X		X	X			2.2 to 2.5	0.97	1.1
	SIM3	X		X	X	X	X	2.3 to 1.6	0.97	0.9
	SW	R		R		R		2.2 to 1.8	0.97	0.9

Discussion

For the annual stepwise regression, latitude, altitude and precipitation were retained while longitude and temperature were eliminated. The retained parameters in the seasonal SW models also indicate the importance of latitude, altitude and precipitation as predictor variables. Latitude survived the stepwise regression for all four seasonal SW models, while altitude and precipitation were both retained in three out of four seasonal models. Precipitation did not survive the regression for the JJA model, indicating that the amount of precipitation is not highly correlated to changes in ¹⁸O_PPT during the summer months. Temperature and longitude remained only in half of the seasonal stepwise regression models (JJA & SON and JFM & MAM, respectively). Longitude is more prevalent in the winter and spring model results, likely due to the enhanced role of continentality controlling the isotopic labeling of precipitation during colder seasons. The temperature difference between the vapour source (ie. ocean) and the condensation site is largest during the winter and spring seasons. The reduced role of temperature and continentality in the summer months is also consistent with the dominance of convective, localized precipitation in the continental interior of Canada during this season (Koster et al., 1993).

Table 4. Seasonal model parameterization and summary statistics. Parameters included in each of the models are indicated with an (X). Stepwise regression retained parameters indicated by an (R).

		Parameterization
Season	Model Type	LAT	LAT^2	ALT	LONG	PRECIP	TEMP	Residual Range ()	Adjusted R^2	Residual Variance ()
DJF	BW	X	X	X				6.5 to 9.0	0.70	16.4
	SIM1	X	X	X				6.0 to 6.6	0.78	12.5
	SIM2	X		X	X			6.5 to 7.2	0.78	12.2
	SIM3	X		X	X	X	X	3.1 to 5.2	0.91	4.8
	SW	R			R	R		3.5 to 5.1	0.91	4.8
MAM	BW	X	X	X				2.8 to 2.9	0.95	2.2
	SIM1	X	X	X				3.0 to 2.8	0.95	2.2
	SIM2	X		X	X			3.4 to 3.0	0.95	2.3
	SIM3	X		X	X	X	X	3.3 to 2.8	0.97	1.7
	SW	R		R	R	R		3.0 to 2.6	0.96	1.7
JJA	BW	X	X	X				4.3 to 2.7	0.95	1.4
	SIM1	X	X	X				3.8 to 2.3	0.95	1.2
	SIM2	X		X	X			4.2 to 2.0	0.94	1.5
	SIM3	X		X	X	X	X	3.5 to 1.9	0.95	1.4
	SW	R		R			R	3.5 to 2.2	0.95	1.4
SON	BW	X	X	X				3.2 to 3.0	0.94	2.5
	SIM1	X	X	X				2.8 to 3.5	0.95	2.4
	SIM2	X		X	X			2.9 to 3.4	0.95	2.5
	SIM3	X		X	X	X	X	2.5 to 2.7	0.97	1.6
	SW	R		R		R	R	2.5 to 2.7	0.97	1.5

Regression coefficients from models not showing multicollinearity are analyzed to determine if they are consistent with the expected controls on the isotopic labeling of precipitation. For all annual and seasonal model parameterizations, latitude and altitude show the expected negative relationship with ¹⁸O_PPT. Longitude, precipitation and temperature consistently show positive correlation with ¹⁸O_PPT for the annual models, however, the magnitude and direction of correlation changes seasonally. Temperature is positively correlated with ¹⁸O_PPT for the SON, DJF and MAM seasons, but becomes negative during the JJA season. Longitude shows a similar but opposing trend to temperature, where negative correlations are present in the SON, DJF and MAM seasons and positive correlations are found in the summer. Precipitation regression coefficients always remain positive.

Evaluation of model performance statistics indicates that annual models perform comparably and overall, satisfactorily ( Table 3 ). On average the models have an adjusted R² value of 0.97 and a residual variance of approximately 0.96. The addition of climate parameters (PRECIP and TEMP) did not result in a noticeable improvement due to the low temporal resolution of this analysis. Comparison of seasonal model type performance ( Table 4 ) show that models with time variant climate parameters (SIM3 and SW models) consistently outperformed models composed solely of geographic parameters (BW, SIM1, and SIM2 models); typically showing higher adjusted R² and lower residual variance. Figure 2 displays seasonal variations in simulated ¹⁸O_PPT for six strategically picked CNIP stations from various climate regions (Environment Canada, 2010) across Canada. Two stations from coastal locations are selected: Bay D'Espoir, NL in Atlantic Canada, and Saturna Island, BC located on the Pacific coast. Calgary, AB and Chapais, QC are chosen as east and west continental stations, representing prairies and northeastern forest climate regions, respectively. Snare Rapids, NWT is located within the Mackenzie District, and is therefore considered to represent continental stations above 55 latitude. Finally, Resolute Bay, NU represents the Arctic tundra CNIP stations.

Figure 2. Annual temporal variations in simulated ¹⁸O_PPT from six select CNIP stations. Measured ¹⁸O_PPT is indicated by a black line.

Analyzing the model results temporally, the MAM, JJA and SON seasonal models all perform satisfactorily; however the DJF models show surprisingly low R² values (0.70 to 0.91), high error variance (4.8 to 16.4), and large residual magnitudes; indicating increased complexity in the driving forces of ¹⁸O_PPT patterns during the winter season (e.g., enhanced kinetic effects during phase changes at sub-zero temperatures, pronounced continentality) resulting in weaker relationships between ¹⁸O_PPT and the predictor variables. Another significant finding is the large difference in performance between DJF model parameterizations solely comprised of geographic parameters and parameterizations including both geographic and climate parameters. The addition of climate parameters to the DJF models significantly improved results and increased R² values from 0.78 (SIM1, SIM2) to 0.91 (SIM3, SW), while decreasing the residual variance substantially (12.5 compared to 4.8). This improvement in performance is not as apparent in the remaining seasons, where all model types perform comparably, with typical R² values ranging from 0.94 to 0.97 and residual variances range from 1.2 to 2.5, suggesting strong relationships between ¹⁸O_PPT and the parameters sets. The probable cause of the poor DJF model performance is due to the weaker relationship between ¹⁸O_PPT and latitude during the winter, when interannual variations in isotope labelling are largest. Linear regression of ¹⁸O_PPT versus latitude for the winter season shows an R² of 0.71, in comparison to R² values ranging from 0.85 to 0.91 for identical regressions completed for the remaining seasons. To better understand model performance, spatial trends in the DJF model residuals are assessed ( Figure 3 ).

Figure 3. Regression model residuals for the DJF model series. Negative residuals (measured <modelled) are shown in black and positive residuals (measured>modelled) are shown in light gray. Residuals are calculated by subtracting modelled ¹⁸O_PPT from measured ¹⁸O_PPT. The size of the circles represents the relative magnitude of the residuals.

There are two types of patterns evident in the residual maps for the DJF seasonal models ( Figure 3 ). Despite the slight differences between the geographically based models and the combined geographical and climatological models, in general all DJF models typically show similar large-scale trends in the sign of the residuals, with the major difference between models being in the magnitude of the residuals. One consistent pattern in the DJF residual maps is the presence of high-magnitude negative residuals within the south-central continental CNIP stations, indicating the models predict more positive ¹⁸O_PPT values than measured. Birks and Edwards (2009) found that circulation-dependent factors play a significant role in the isotopic signature over central Canada, which is where a large portion of the negative residuals (ie. model over prediction) are occurring. The sets of parameters used in the winter models may not adequately capture the depletion of heavy isotopes from the longitudinal transport of air masses from the Atlantic or Pacific oceans inland or the influence of depleted Arctic air masses over central Canada. Additionally, in this region measured ¹⁸O_PPT data can show anomalously negative or positive values when associated with more zonal or meridional flow, respectively (Birks and Edwards, 2009). Based on the winter model results, it is apparent that longitude does not make a significant enough model contribution to progressively deplete ¹⁸O_PPT signatures inland from the coast and cannot account for the different isotopic labeling occurring at a single geographic location under different circulation patterns or multiple moisture sources. A possible solution to this problem is the introduction of separate models for various climate regions in future model iterations.

Another trend visible in Figure 3 is the extremely positive residuals at a group of western CNIP stations between 55 and 70 N latitude (Snare Rapids, Yellowknife, Whitehorse, Fort Smith, Mayo), indicating that the model simulations predict more negative ¹⁸O_PPT than measured. The isotope time series from Fort Smith and Mayo are both less than five years in length and there is no overlap between the time periods for these stations (19851989 and 19611965, respectively). Therefore, the large positive residuals at this location could be related to biases in the data used to develop the models. However, it may also indicate unaccounted processes or sources of moisture. As with the south central anomaly, the positive residuals could indicate a circulation-dependant control on the isotopic composition of precipitation in this location, similar to the anomalously positive ¹⁸O_PPT values that occur when meridional transport of moisture across the Rockies results in less rainout.

The anticipated trend towards more negative ¹⁸O_PPT values as you move from the coast towards the interior stations is present in all of the winter models; however, the winter models based solely on geographical parameters largely under predict at coastal locations. The models that incorporate climate parameters have small to medium residuals at these locations indicating that the addition of climate parameters to the seasonal models may help to capture some of the physical processes behind the seasonality in isotope patterns, especially at coastal locations where isotope patterns are influenced by high amounts of precipitation and moderate temperatures. Spatial maps similar to those shown in Figure 3 created for the remaining MAM, JJA and SON models also show similar findings; all models predict similar large-scale patterns and the addition of climate parameters to regression models improved the evaluation statistics and decreased the magnitude of residuals. This is further highlighted in a series of plots detailing the seasonal variations in model residuals ( Figure 4 ) for all 33 CNIP stations.

Figure 4. Residuals from seasonal regression models plotted against longitude for each of the five model types.

Figure 4 establishes that generally each model type provides similar spatial residual patterns for each respective season. Therefore, the remainder of the discussion will focus on the analysis from the SIM3 model because it consistently out-performs the other model types ( Table 3 , Table 4 ). A time-series of residual maps based on the seasonal models from the SIM3 regression is created to spatially assess the changes in residuals annually and from season to season ( Figure 5 ).

The distributions of residuals for annual, MAM, JJA, and SON show different trends than those introduced for the DJF model. Across central Canada a shift from negative residuals in the winter to a mix of negative and positive residuals in the remaining months occurs, indicating no major model bias exists at this location for these months. The SIM3 model underestimates the ¹⁸O_PPT for stations located near the Great Lakes for most of the year, until fall when there is a reversal in the sign of the residuals indicating the model overestimates ¹⁸O_PPT for this season, a trend similar to the annual models. Fall is a time of high evaporation in the vicinity of the Great Lakes and the switch in the sign of the residuals may reflect increased contributions of recycled moisture (Van der Ent et al., 2010) or the moderating effect of the lakes on local temperature. Isotopic evidence for this recycling has been noted in the higher d-excess of precipitation and river runoff in the vicinity of the North American Great Lakes (Gat et al., 1994). Generally, high Arctic residuals tend towards negative residuals in the winter and spring seasons and positive residuals in the summer and fall seasons, indicating the model is not capturing the full seasonality inherent to this area.

Figure 5. Regression model residuals for the annual, DJF, MAM, JJA and SON SIM3 models. Negative residuals (measured<modelled) are shown in black and positive residuals (measured>modelled) are shown in light gray. Residuals are calculated by subtracting modelled ¹⁸O_PPT from measured ¹⁸O_PPT. The size of the circles represents the relative magnitude of the residuals.

Spatial Distribution of Annual and Seasonal ¹⁸O_PPT over Canada

Annual and seasonal contour maps from SIM3 model simulations are created for the Canadian domain through use of Canada3D, a digital elevation model (DEM) with a horizontal grid spacing of 30 arc-seconds. Elevation at atitude and longitude pairs is obtained by sampling the DEM at a 0.250.25 resolution. Gridded temperature and precipitation annual and seasonal long-term averages are extracted from the National Centers for Environmental Prediction (NCEP) North American Regional Reanalysis (NARR) dataset. NCEP reanalysis data are provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Website at http://www.esrl.noaa.gov/psd/. NARR is a long-term high-resolution climate dataset for the North American domain at a grid resolution of approximately 0.3 degrees (32km). Annual and seasonal climatologies are derived from monthly mean climatologies, and subsequently used as input into the SIM3 model. The SIM3 model simulation of ¹⁸O_PPT across Canada annually and seasonally is presented in Figure 6 . The ¹⁸O_PPT bands visible in the figure are due to the dependence of oxygen-18 on latitude, and deviations from the banding effect are due to altitude or climate effects. The ¹⁸O_PPT pattern becomes more depleted in heavy isotopes as it travels poleward, which is expected due to the removal of moisture from air masses (Yurtsever and Gat, 1981). The Maritimes indicate more enriched ¹⁸O_PPT patterns, which is predominantly due to their close proximity to the ocean and lower latitude location.

Figure 6. SIM3 annual and seasonal ¹⁸O_PPT patterns across Canada.

The final maps of the distribution of isotopes in precipitation generated in similar studies (Bowen and Wilkinson, 2002; Dutton et al., 2005; Liu et al., 2008), are based on correcting an estimated ¹⁸O_PPT map by applying a contoured residual field. The emphasis in this study is to develop improved empirical relationships to approximate the isotopic labeling of precipitation across Canada thereby decreasing the residuals and the need for correction. For this reason we have presented our estimated ¹⁸O_PPT maps ( Figure 6 ) and the residuals ( Figure 7 ) separately. Interpolated maps of the annual and seasonal residuals ( Figure 7 ) are created as companions to Figure 6 , to give an estimate of the uncertainty of the SIM3 model results across Canada. However, it should be noted that additional uncertainty exists when interpolating within areas where there are very little data, therefore application of the residual map contours to the ¹⁸O_PPT estimates to arrive at a final isotope signature should be carried out with caution.

Conclusions

This research has developed and evaluated five predictive ¹⁸O_PPT models capable of providing estimates of the average annual and seasonal composition of Canadian ¹⁸O_PPT using a combination of geographical and climatological parameters. The statistics describing the model-data comparisons show satisfactory model results for annual, spring, summer and fall models. However, winter ¹⁸O_PPT predictions show significantly lower adjusted R² and higher residual variances than the other seasonal models. Seasonal model results improve with the addition of climate parameters (temperature and precipitation), this finding being especially significant in the winter season. There are still aspects of model development and performance that could benefit from improvement before implementation within the isoWATFLOOD framework. Our results suggest that inclusion of additional parameters selected to capture some of the influence of circulation patterns on the isotopic labeling of Canadian precipitation should be useful for developing improved predictive ¹⁸O_PPT models suitable for finer spatial and temporal scales. Additional parameterization needs to be carefully scrutinized so as not to limit the practicality of the model to hydrological applications where forcing data adds unnecessary uncertainty to model outcomes. Predictions of ¹⁸O_PPT should ideally be within limits of analytical uncertainty (0.2 for ¹⁸O_PPT) based on reported laboratory standards. Additional data or a change in model development methodology is required to refine the model to this level of accuracy. For future studies, it is suggested that the Canadian domain is separated into climate zones, and that several ¹⁸O_PPT models are created, each reflecting the conditions of their respective climate zone. This methodology may strengthen the relationships between ¹⁸O_PPT and the model parameters. Lastly, models should be created for both oxygen-18 and deuterium to allow for the creation of a simulated isotope framework, and the use of d-excess as a diagnostic tool. This study shows great promise in determination of a relationship to accurately predict ¹⁸O_PPT across Canada. The new annual and seasonal maps of the distribution of ¹⁸O_PPT in Canadian precipitation will provide an improved basemap for climatological, ecological, hydrological and hydrogeological studies requiring information on the long-term seasonal and annual labeling of precipitation. Upon completion, the models will be embedded within isoWATFLOOD to enable internal simulations of ¹⁸O_PPT at the basin scale, and facilitate improved regional simulations of water isotopes in hydrologic compartments.

Figure 7. SIM3 annual and seasonal residual contours.

Acknowledgements

The authors would like to thank all of the reviewers who provided their comments and insight, helping to improve this paper. This research has benefited from discussions with Dr. John Gibson, Dr. Yi Yi, and Dr. Tom Edwards. This research was funded by a University of Manitoba Graduate Scholarship, a Natural Sciences and Engineering Research Council of Canada (NSERC) Canada Graduate Scholarship, and Alberta Innovates-Technology Futures. Additional gratitude is expressed to the Canadian Network for Isotopes in Precipitation for providing the updated isotopes in precipitation data necessary to conduct this research.