In this time of a changing climate, it is important to know whether lake levels will rise, potentially causing flooding, or river flows will dry up during abnormally dry weather. The Great Lakes region is the largest freshwater lake system in the world. Moreover, agriculture, industry, commerce, and shipping are active in this densely populated region. Environment and Climate Change Canada (ECCC) recently implemented the Water Cycle Prediction System (WCPS) over the Great Lakes and St. Lawrence River watershed (WCPS-GLS version 1.0) following a decade of research and development. WCPS, a network of linked models, simulates the complete water cycle, following water as it moves from the atmosphere to the surface, through the river network and into lakes, and back to the atmosphere. Information concerning the water cycle is passed between the models. WCPS is the first short-to-medium-range prediction system of the complete water cycle to be run on an operational basis anywhere. It currently produces two forecasts per day for the next three days. WCPS generally provides reliable results throughout the length of the forecast. The transmission of errors between the component models is reduced by data assimilation. Interactions between the environmental compartments are active. This ongoing intercommunication is valuable for extreme events such as rapid ice freeze-up and flooding or drought caused by abnormal amounts of precipitation. Products include precipitation; evaporation; lake water levels, temperatures, and currents; ice cover; and river flows. These products are of interest to a wide variety of governmental, commercial, and industrial groups, as well as the public.
Environment and Climate Change Canada’s (ECCC) operational implementation of a numerical forecast system linking atmospheric, surface, river and lake/ocean models is described.
The Water Cycle Prediction System (WCPS) was recently implemented over the Great Lakes and St. Lawrence River watershed (WCPS-GLS version 1.0) by Environment and Climate Change Canada’s (ECCC) Canadian Centre for Meteorological and Environmental Prediction (CCMEP, formerly CMC). WCPS simulates the complete water cycle, following water as it moves from the atmosphere to the surface, through the river network and into lakes, and back to the atmosphere. The implementation of WCPS-GLS, following a decade of research and development, supports Canada’s obligations under the Boundary Waters Treaty (International Joint Commission 2016). The Great Lakes region is the largest freshwater lake system in the world. It is a densely populated region that is important for agriculture, industry, and commerce. Commercial shipping and recreational boating are active, in particular because of the Great Lakes–St. Lawrence Seaway System that connects all the Great Lakes—through the St. Lawrence River—to the North Atlantic Ocean.
WCPS is a chain of interconnected models (Fig. 1). The individual models represent processes in the atmosphere, at the land surface and in the soil (at 10-km resolution), in large bodies of water and marine ice (at 2-km resolution), and in rivers (at 1-km resolution). WCPS-GLS produces two forecasts per day for the next three days. Products include spatially varying precipitation, evaporation, river discharge, water level anomalies, surface water temperatures and currents, and ice concentration and thickness (see “Products generated by WCPS-GLS” section for further information). This information is useful to water resources and management authorities, flood forecasters, hydroelectricity producers, navigation, environmental disaster managers, search and rescue teams, agriculture, and the public.
To our knowledge, WCPS is the first system to predict the complete water cycle at a short to medium range on an operational basis anywhere. Representing the full cycle enables us to assess whether the water balance at the surface is closed. Previously, Nagai et al. (2011) described a numerical forecast system of the complete water cycle. This system linked atmosphere, ocean, waves, hydrology, and land surface models. However, there is no indication that this system was ever run on an operational basis (i.e., on an ongoing, scheduled basis in real time with 7 days per week support).
The Global Flood Awareness System (GloFAS; Alfieri et al. 2013; Dottori et al. 2016; European Commission 2016) predicts the water cycle daily on a quasi-operational basis. GloFAS was jointly developed by the European Commission and the European Centre for Medium-Range Weather Forecasts (ECMWF). This system couples ECMWF’s Variable Resolution Ensemble Prediction System (VarEPS; Miller et al. 2010) and its land surface module, Hydrology Tiled ECMWF Scheme of Surface Exchanges over Land (HTESSEL; Balsamo et al. 2009, 2011), with the LISFLOOD hydrology model (van der Knijff et al. 2010). However, the river discharges are not used as boundary conditions for the ocean model. Consequently, the water cycle represented is incomplete.
Similarly, the recently introduced UCK2 regional coupled environmental prediction system (Lewis et al. 2018) predicts the water cycle over the United Kingdom. UCK2, developed by the Joint Weather and Climate Research Programme and the Met Office, couples models for the atmosphere (Met Office Unified Model; Brown et al. 2012), land surface [Joint UK Land Environment Simulator (JULES); Best et al. 2011; Clark et al. 2011], ocean [Nucleus for European Modelling of the Ocean (NEMO); Madec et al. 1998, 2008], and ocean waves (WAVEWATCH III; Tolman 2009). The land surface and ocean models are not coupled via riverine discharges. This system is expected to run operationally in the near future.
In both France (Coustau et al. 2015) and the United States (National Oceanic and Atmospheric Administration 2016), forecasts of river discharges are produced daily on an operational basis. Météo-France produces national medium-range hydrological forecasts by coupling ECMWF’s ensemble forecasts with a land surface model and a hydrologic model. In the United States, the Weather Research and Forecasting (WRF) Model hydrological modeling system (WRF-Hydro), driven by an assortment of observations and model output, produces river forecasts of varying lengths. However, to our knowledge, in neither country does a short- or medium-range ocean/lake model ingest the forecast discharges as boundary conditions.
The aim here is not to provide a detailed evaluation of all model components (Fortin et al. 2017), but rather to answer three particular questions: What is the overall quality of WCPS-GLS? Does the quality of the forecast degrade significantly as a function of lead time? How does error propagate through the system? To address these questions, the next section describes the component models of WCPS and their linkages, the third section outlines the observations used to verify the system’s performance, and the fourth section verifies the performance of the individual component models. Products generated by WCPS-GLS are presented in the fifth section. The sixth section outlines future development of WCPS while the seventh section presents our conclusions.
THE COMPONENT MODELS AND THEIR LINKAGES.
WCPS is a chain of interconnected numerical models (Fig. 1). A detailed description is provided in the appendix. Additional information and a thorough evaluation are also available in the technical documentation (Fortin et al. 2017). Briefly, WCPS couples ECCC’s atmospheric Global Environmental Multiscale (GEM) model (Côté et al. 1998a,b; Girard et al. 2014); the surface and soil model Interactions between Soil, Biosphere, and Atmosphere (ISBA; Bélair et al. 2003a,b); the lake/ocean model NEMO (Madec et al. 1998, 2008); the Los Alamos Sea Ice Model (CICE; Hunke and Lipscomb 2010); the river-routing model WATROUTE (Kouwen 2010) for flow through the river network; and the Coordinated Great Lakes Regulation and Routing Model (CGLRRM; Tolson 2009) for flow through the lakes’ interconnecting channels. Together, these models simulate the complete water cycle, albeit in varying levels of detail for different aspects of the cycle.
The component models of WCPS pass information to each other both offline and online during model integration (see appendix). GEM–ISBA passes estimates of surface runoff, or excess surface water, to WATROUTE. WATROUTE also requires estimates of recharge, or water that enters the river from the base of the soil column. Currently, climatological (2004–09; Deacu et al. 2012) recharge is used since it is overestimated by the operational surface data assimilation system (see appendix). WATROUTE provides NEMO–CICE with time-varying riverine discharges. GEM and NEMO–CICE interactively exchange fluxes of heat, mass, and momentum across the marine boundary layer at every common model time step. This tight coupling allows a coherent response to rapidly changing conditions in the atmosphere or lakes (e.g., due to the formation of lake ice).
A variety of products are used to verify the forecast skill of the various model components (Fortin et al. 2017). This includes observations from meteorological and hydrometric stations, satellite retrievals of sea surface temperature and lake ice, surface flux measurements from flux towers, and vertical profiles of water temperature from moorings. Comparisons against estimates from independent models are also used. Satellite-derived lake levels are currently at too coarse a horizontal resolution to be of interest for the Great Lakes (National Aeronautics and Space Administration 2001).
Surface currents are important for a variety of applications, such as search and rescue and oil spill response. Since it is difficult to evaluate the quality of surface currents directly, it is assessed indirectly by verification of winds, water level, and vertical stratification in the lakes. In the future, surface currents will be verified using drifters and observations at fixed locations. Satellite-derived estimates of surface currents are of little interest in the dynamically active Great Lakes as they reflect mean currents rather than eddies (Heuzé et al. 2017).
For many parameters only a single type of verifying observation is available. In contrast, for surface water temperature, both in situ observations and satellite retrievals are available. Simulated temperatures have been verified separately against these two data types. A simultaneous comparison would be of interest. In situ observations are useful since they provide precise information, albeit with gaps in coverage. Satellite-derived data are useful in that they cover an extended area, albeit with less precision.
The verification presented in the current study is limited to parameters directly related to the hydrological cycle. We verify model output against analyses of a variable’s spatial distribution and against individual observations. The observations used to evaluate precipitation, sea ice cover, and lake levels are described below. Specifically, we compare simulated precipitation (Fig. 2) to the regional deterministic configuration of the Canadian Precipitation Analysis (CaPA-RDPA; Lespinas et al. 2015). We used version 4.0 (v4.0) of CaPA; in addition to the gauge networks identified in Lespinas et al. (2015), CaPA v4.0 assimilates Canadian and U.S. radar data following Fortin et al. (2015). Figures 6 and 7 present analyses from Regional Deterministic Prediction System (RDPS). See the appendix for the generation of these analyses. We verify the simulated lake ice concentration of Fig. 9 against daily ice charts from the Canadian Ice Service (Canadian Ice Service 2016).
The verification of the component models’ products in Figs. 3–11 uses the individual observations identified in Table 1. For lake levels, only international gauging stations were used for verification. These monitoring stations are used to establish water levels in a body of water shared by Canada and the United States and to determine transboundary streamflows. Their designation and operating procedures are coordinated by the two countries. During the period of this study, Port Stanley was offline due to waterfront construction.
Two corrections were applied to simulated lake levels prior to verification. Since WCPS-GLS does not currently represent the varying effects of ongoing glacial isostatic adjustment on Great Lakes water levels, estimated adjustments (Bruxer and Southam 2008) were applied to the simulated level at each observation station for verification purposes. Additionally, simulated lake levels had drifted, likely as a result of using climatological (2004–09) riverine inflows as boundary conditions for NEMO’s simulations from 2004 to 15. Therefore, the mean bias in lake level for 15 June–15 October 2016 over all international gauging stations on a given lake was removed from all simulated levels at all gauging stations on that lake. Note that, although a long-term bias can be removed in this way, analyzed riverine inflows are a preferred source of boundary conditions as their temporal variations affect lake dynamics.
NEMO–CICE does not assimilate directly the ice charts of the Canadian Ice Service (Fig. 9). However, NEMO–CICE assimilates and the Canadian Ice Service incorporates the same estimates of ice concentration and thickness derived from RadarSat-2 satellite images. Thus, the simulated and verifying estimates of ice cover are not completely independent. Of the various products used for the verification of all other parameters, only surface-level atmospheric temperature and riverine discharge are also assimilated by WCPS-GLS (see appendix).
Warm season case study: 25 June–2 July 2016.
As an introduction to WCPS-GLS and its representation of the water cycle, we present a case study of 25 June–2 July 2016. (A cold-season case study and a verification of the long-term performance of WCPS-GLS follow in the next two subsections.) During this period, two waves of precipitation passed through the Great Lakes’ watershed from west to east (Fig. 2). We expect the rivers and lakes to respond to the precipitation’s passage. Furthermore, since lake levels during this period were in no way unusual for the season (J. P. Smith et al. 2016), hydrological drivers were likely well balanced. Thus, the response to the precipitation may be considered typical, rendering the period a suitable test of the prediction system’s performance.
The accuracy of the simulated precipitation varies within the Great Lakes’ basin. During the period’s first wave of precipitation (25–27 June), simulated 24-h accumulated precipitation is underestimated in the western part of the domain in comparison to the CaPA analyses (see the “Verifying observations” section) and overestimated in the east. During the second wave (29 June–1 July), precipitation is initially underestimated in the west, overdeveloped in the central-west region, and then underestimated in the eastern half of the watershed. The online supplemental material (https://doi.org/10.1175/BAMS-D-16-0155.2) provides animated forecasts coinciding with the second wave of precipitation, riverine discharge, surface water temperature, and changes in lake level.
Differences between simulated and observed riverine discharges are shown in Fig. 3 for the most-downstream and most-upstream stations. The most-upstream stations are closest to a river’s headwaters, though they are not necessarily in the headwaters. No discharge observations are assimilated upstream of the most-upstream stations. Similarly, although there are no other stations between the most-downstream stations and the lakes, these stations are not necessarily at the river’s mouth. Although the most-downstream stations are the stations that benefit most from upstream data assimilation, 50%–79% of the most-downstream stations are also that river’s most-upstream station (Table 2); there is only one observation station on the river.
Differences between riverine discharges as simulated and observed respond to the passage of the precipitation waves. For example, in response to the first wave of precipitation, the discharge differences between Lakes Huron, Erie, and Ontario increase rapidly then diminish (Fig. 3; see Fig. 11 for lake locations). The differences reflect local errors in precipitation. As expected, differences at stations not benefiting from upstream data assimilation are noticeably greater than at those that do benefit; data assimilation reduces discharge errors caused by precipitation errors and/or issues within the land surface and/or river-routing model.
Deacu et al. (2012) investigated the impact of precipitation errors on simulated streamflows in the Great Lakes region from 2004 to 2009. They found that driving the model with an early experimental version of CaPA (Mahfouf et al. 2007) versus precipitation simulated by a regional version of GEM with a 15-km horizontal resolution (Mailhot et al. 2006) degraded the Nash–Sutcliffe model efficiency coefficient for Lake Erie but improved it by 0.06 for Lake Ontario. Deficiencies found in this early version of the analysis (Carrera et al. 2010) have since been corrected.
More recently, experiments were performed in the Canadian Rockies for October 2014–September 2015 (P.-M. Farah 2016, personal communication) where WATROUTE was forced alternately by runoff driven by CaPA analyses and by precipitation estimated by the High Resolution Deterministic Prediction System (HRDPS; Milbrandt et al. 2016). Version 3.0 of CaPA replaced version 2.4 in mid-November 2014. Version 2.4 assimilated surface synoptic observations (SYNOP), aviation routine weather reports (METARs), Réseau Météorologique Coopératif du Québec (RMCQ), and American meteorological cooperative network [Standard Hydrometeorological Exchange Format (SHEF)] precipitation observations (Lespinas et al. 2015). In version 3.0, the assimilation of radar data was added (Fortin et al. 2015). Averaged over 32 observation stations, the temporal variability improved with CaPA precipitation: the Nash–Sutcliffe model efficiency coefficient and the correlation coefficient increased by 19.3 and 0.6, respectively. The volumes of the simulated discharge also improved: the mean absolute error decreased by 48% and the percent bias improved by 94%.
The underestimation of the simulated precipitation over Lake Superior’s watershed on 25 June is reflected in the underestimation on 25–27 June of the inflow to the lake from the group of most-downstream stations (brown lines in Fig. 4). Interestingly, the forecasts ending from 1800 UTC 27 June to 1800 UTC 28 June, where the event falls approximately within the forecast’s day 2, missed the event completely. In contrast, forecasts launched earlier and ending from 0600 UTC 26 June to 0600 UTC 27 June, where the event represents the forecast’s day 3, agree as a group better with the observations. Given the absence of a single discharge observation on some of the watershed’s rivers, we are unable to formally evaluate the accuracy of the forecast total inflows (blue lines in Fig. 4). However, since the earlier forecasts of inflow from the group of most-downstream stations agree best with the observed event, it seems likely that the higher estimated total inflows of the earlier forecasts are also closer to reality. In contrast, the quality of the forecast inflow to Lake St. Clair from the most-downstream stations deteriorates with lead time; many brown forecast lines end at values considerably higher than the observed inflows. Interestingly, forecasts tend to be far more consistent for Lakes Michigan and Huron (not shown) than Lake St. Clair; localized errors are likely smoothed out in larger watersheds.
Mean lake levels are strongly affected by inflows, which drive multiday variations in lake levels. Spatial variations in levels are mostly linked to wind and atmospheric pressure changes, which manifest themselves as intraday variations in lake levels (Hamblin 1987). Note that, while the version of NEMO used in WCPS-GLS assumes constant surface pressure, the surface of the lake does respond to the wind stress, creating seiche events. The lack of inverse barometer effects can be expected to reduce the amplitude of storm surge events. Simulated lake levels, debiased by removing each station’s mean bias over the period shown (see the “Verifying observations” section), are presented in Fig. 5 (see Fig. 11 for station locations). NEMO is unable to capture the high-frequency (subdiurnal) variations seen on Lakes Superior and Ontario. This may be a forcing issue (10-km winds driving a 2-km lake model), be due to NEMO model physics (e.g., it does not represent nonhydrostatic seiches, is possibly overly diffusive, or is lacking in resolution), or reflect the lack of an explicit coupling to a wave model (e.g., Powers and Stoelinga 2000). NEMO, however, does reproduce the low-frequency (multiday) changes in lake levels on Lakes Michigan and Huron. NEMO also represents the medium-frequency (diurnal) waves on Lake Erie. The alternating 25+-cm fluctuations on 1–2 July likely represent a wind-driven seiche. NEMO does not capture the full magnitude of the initial peak at Toledo, Ohio, but does forecast subsequent peak levels, particularly at the longer lead times.
The lakes’ forecast response to the passage of the second wave of precipitation is examined in Figs. 6–8. During this period, surface (10 m) analyzed winds are highly variable (Fig. 6, right column). As expected, given their smoother surface, winds are noticeably stronger over the lakes.
Forecast surface currents in the Great Lakes (Fig. 6, left column) are primarily wind driven (Harrington 1894). However, marine circulations driven by density variations, bathymetry, and riverine and lake discharges also impact surface currents. The current from the south shore of the eastern half of Lake Superior toward the lake’s center at 0600 UTC 30 June is wind driven. In contrast, the offshore-oriented surface current in Lake Superior’s northwest corner is likely driven by the significant (300–400 m3 s–1) inflow from the Nipigon River. This river is supplied by Lake Nipigon and its watershed and by water diverted from Lake Ogoki into Lake Nipigon. The inflow’s impact may be heightened by the relative shallowness of the nearshore area; the surface current stops abruptly and the surface water cools noticeably where the lake deepens offshore (Figs. 6 and 7; National Centers for Environmental Information 1999). Similarly, the steep bathymetric ridge in central Lake Michigan likely disrupts the region’s wind-driven surface currents at 0600 UTC 1 July.
In general, the direction of the simulated surface currents is likely realistic, given the generally good agreement between observed (right column) and forecast (left column) changes in lake level at individual locations (Fig. 6). However, the magnitude of the change in lake level is often underestimated. Thus, current strength may be underestimated. The lack of an explicit coupling to a wave model likely contributes to errors in turbulent fluxes between the atmosphere and the lakes with effects on the mixed layer depths and surface currents (e.g., Powers and Stoelinga 2000).
Surface water temperatures can vary greatly over the Great Lakes in response to the watershed’s 10° range in latitude, differences in atmospheric surface-level temperature, and the 128-m range in mean lake depths (Environmental Protection Agency 2016). In Fig. 7, observed and forecast surface water temperatures are warmest on Lake Erie, the most southerly lake and the shallowest, and coolest on Lake Superior, the most northerly and deepest. Simulated surface water temperatures cool offshore as depth increases.
Wind action, riverine inflow, and lake discharge also impact surface water temperature by mixing or exchanging warmer surface waters and colder deep waters. On Lake Michigan, in the presence of strong surface winds and currents, observed and forecast coastal waters cool rapidly on the southeast side from 1 to 2 July while remaining warm on the western side. This suggests an overturning circulation with cold waters rising on the east and warm waters sinking on the west. Similarly, wind-driven upwelling and downwelling occur along the north and south shores of Lake Ontario, respectively, on 2–3 July (not shown) forming a pool of cool water that spreads south.
A cold bias in the forecast surface water temperature at the western tip of Lake Superior is caused in part by a greater model depth at the observation location. Since the simulated water is deeper, it is colder and less likely to warm as rapidly as observed. Similarly, the warm bias of the two northernmost stations on Lakes Michigan and Huron are caused in part by the lesser model depth at the observation location. The 2-km NEMO configuration is not able to resolve the finescale structure of coastal areas. In the future, we plan to add a 2D hydraulic model (Matte 2014; Matte et al. 2017a,b) to WCPS-GLS to simulate processes in the lake’s coastal regions.
Evaporation from 30 June to 3 July 2016 was forecast to be strongest over Lakes Michigan and Erie (Fig. 8). At this time, surface-level air was very dry (not shown). In contrast, low-level atmospheric humidity was high over Lake Superior, where condensation was forecast. Moreover, Lake Superior’s cooler surface waters (Fig. 7) likely cooled the low-level air mass, promoting condensation. Forecast and observed (Blanken et al. 2011; Spence et al. 2011, 2013) evaporation/condensation from the lakes are in good general agreement. In Lake Superior, the near-zero observed values are at most times collocated with or adjacent to forecast regions of near-zero evaporation/condensation. Similarly, in Lake Ontario, the strong observed evaporation agrees well at most times with forecast values either at the observation location or nearby. Evaporation is at times underestimated on Lake Erie.
Cold season cold snap: 11–18 February 2016.
To illustrate the capacity of the system to capture important changes in lake ice cover, we now present results from a rapid freeze-up event (Fig. 9). Such events are fairly typical in the Great Lakes and occur most years. A full evaluation of ice forecast skill is presented in Fortin et al. (2017) and shows that WCPS-GLS captures well the seasonal cycle of ice conditions with forecast errors typically of a similar magnitude to uncertainties in ice charts. In Fig. 9, simulated lake ice concentration on Lake Erie at 1800 UTC over a series of days are presented alongside daily ice charts from the Canadian Ice Service (Canadian Ice Service 2016; see the “Verifying observations” section). Additionally, the 48-h forecast of the concentration of lake ice from 1200 UTC 11 February 2016 is shown in Fig. DE1.
The forecasts and charts of Fig. 9 show an initial rapid buildup in the western and central parts of the lake followed by a gradual reduction in ice cover along the north shore. However, the model apparently overestimates ice concentrations in the lake’s central region. The available RadarSat observations indicate a slight overestimation in simulated ice concentration on 13 February, gradually transitioning to a slight underestimation on 15 February. Thus, actual concentrations may be between the forecasts and charts of Fig. 9. (Note that simulated and retrieved estimates on Lakes Superior, Michigan–Huron, and Ontario differed by no more than 3% on average over the period.) The model’s overestimation of ice concentrations on Lake Erie is likely explained by simulated surface-level atmospheric temperatures that were too cold during this week by up to 3°C on the western and southwestern edges of Lake Erie.
Long-term performance: 15 June–15 October 2016.
To provide an overview of the ability of WCPS-GLS to forecast the state of the atmosphere, rivers, and lakes in the Great Lakes region, we verify the forecasts from 15 June to 15 October 2016 against the observations listed in Table 1. We consider the bias, or the mean difference between the simulated and observed values (simulated minus observed) over the verification period. We also consider the standard error, which we define as the standard deviation of the difference between the simulated and observed values. Figure 10 presents the bias and standard error of two basic aspects of each modeling compartment. Figure 11 presents the spatial distribution of the differences shown in Fig. 10. In Fig. 10, the mean error over all observation stations is given as a function of forecast hour. Errors for forecast hour 0 are calculated from the analysis used to initialize the model. In Fig. 11, the mean error over the full forecast (hours 0–84) is given for each station. Figures 11d and 11e show the locations of the lakes and the international gauging stations, respectively. For the precipitation, we provide only the standard deviation of differences for the entire Great Lakes watershed given this score’s high degree of noise.
Several main themes emerge from the long-term verification of WCPS-GLS, including how error propagates between model compartments. For instance, both atmospheric surface-level temperatures and surface water temperatures are too cold. Since the atmospheric and lake models are tightly coupled (see appendix), these cold biases are likely connected. The surface water temperature’s bias is nonzero at forecast hour 0. In contrast, the bias of the 2-m atmospheric temperature is initially small and increases with forecast hour. Thus, the error likely originates with NEMO’s pseudoanalysis of surface water temperatures and subsequently propagates into the surface-level atmosphere during the forecast.
Errors in estimates of water volumes also propagate within WCPS-GLS between component models. Precipitation, riverine discharge, inflow to the lakes, and lake levels tend to be predicted more accurately in the western part of the watershed (Lakes Superior and Michigan–Huron) than the eastern part (Lakes St. Clair, Erie, Ontario). In the western region, the precipitation received during the verification period was approximately normal versus the 15-yr (2002–16) average precipitation as estimated by CaPA analyses (M. Beauchemin 2016, personal communication). Furthermore, ∼80% and ∼40% of the discharge within the watersheds of Lakes Superior and Michigan–Huron, respectively, represent the persistence below a major dam of the most recent observed discharge (Table 3). In contrast, in the eastern region, the amount of precipitation received in August 2016 was generally normal to ∼50% above normal. However, in all other months from May to October 2016, the eastern region received ∼50% less than the normal precipitation. Additionally, no discharge in the eastern part of the watershed represents the persistence of an observation. Hence, in the eastern region we have an abnormally small amount of precipitation during the verification period combined with a less tightly constrained river-routing model.
The abnormally dry conditions of the verification period are likely responsible for the tendency in the eastern region of precipitation, riverine discharge, and inflow to the lake to be overestimated. Errors in precipitation are likely propagating into the river discharges. The current use of climatological recharge in WCPS-GLS combined with an unusually dry season likely further increases the positive discharge biases. The discharge errors then propagate into the inflows. Finally, errors in inflow likely propagate into the lake levels both directly from the lake’s immediate watershed and indirectly through the lakes’ connecting channels.
Two additional sources of error in lake levels must also be considered. It is highly likely that part of the error in the estimated flow through the connecting channels between lakes is attributable to errors in the precipitation that falls directly onto the surface of the lakes. The magnitude of errors in precipitation onto lake surfaces is difficult to evaluate given the lack of midlake precipitation gauges. However, the CaPA precipitation analyses of Fig. 2 incorporate overlake radar observations. There is no evidence in Fig. 2 of a difference in error characteristics between overlake and overland precipitation. Furthermore, the changes in lake levels of Fig. 5, which are forced partly by overlake precipitation minus evaporation (P − E), are reasonable.
The interlake flows, which also drive lake levels, are represented by the Coordinated Great Lakes Regulation and Routing Model (CGLRRM; see appendix). We have not evaluated the performance of CGLRRM. This will be addressed during the update of CGLRRM (see “Future work” section).
Despite the evidence that erroneous volumes of water are being transferred between model compartments, not all such errors propagate unchecked through the modeling system; the assimilation of riverine discharge observations reduces discharge errors (“Warm season case study” subsection). Interestingly, even though the watersheds of Lakes St. Clair, Erie, and Ontario have the lowest fractions of the most-downstream stations that are also the most-upstream stations (Table 2) and the most-downstream stations benefit most from data assimilation upstream (“Warm season case study” subsection), these lakes’ discharges and inflows are characterized by a strong error growth with lead time from forecast hour 48. This apparent contradiction is likely explained by the fact that the assimilated discharge observations have no impact on upstream values; the benefit of the assimilated discharge propagates downstream with the observed discharge (see appendix). Consequently, the fraction of verified discharges that are purely natural and unaffected by data assimilation increases with forecast hour. Thus, the errors at the longer lead times likely represent the increasingly unchecked propagation of erroneous water volumes in an abnormally dry region.
Finally, error growth with forecast hour is low for precipitation, lake levels, and surface water temperatures. In the western part of the watershed (Lakes Superior and, Michigan–Huron), which received more normal precipitation amounts than the eastern part and where discharges are more tightly tied to observations (Table 3), errors also grow slowly for riverine discharges and inflows to the lakes. The persistence of discharge observations downstream of major dams in the watersheds of Lakes Superior and Michigan–Huron undoubtedly reduces error growth with lead time for discharge and inflow (Table 3). Nonetheless, the low error growth with lead time of multiple variables suggests that 84-h forecasts of the hydrological cycle have the potential to be useful through day 3.
Statistical postprocessing might improve forecasts given that all variables verified in Fig. 10 are characterized by a systematic bias that changes with lead time and varies from lake to lake. As experience is gained with the system, it will be possible to implement a model output statistics (MOS) postprocessor to reduce this bias and possibly also reduce the standard error. An adaptable method will need to be used as the system is expected to be upgraded on a yearly basis to benefit from improvements in numerical weather prediction models used operationally at ECCC. This approach is expected to work well for atmospheric and surface variables, as it is already used operationally for weather forecasting (Wilson and Vallée 2002). However, streamflow forecasts are notoriously difficult to bias correct in cold regions because streamflow exhibits a strong seasonality, with most of the large flow events occurring during the spring freshet. For this reason, it is expected that a reforecast experiment—in which the current modeling system is run back in time for many years—will be required. Such experiments have already been conducted for individual components of the forecasting system (Deacu et al. 2012; Dupont et al. 2012) but not yet for the coupled system.
PRODUCTS GENERATED BY WCPS-GLS.
Table 4 provides a list of currently available analyses and forecasts. The list includes products from the hydraulics model Simulation Hydrodynamique OPérationnelle (SHOP) (Matte 2014; Matte et al. 2017a,b). SHOP was not described above as it is not currently running on an operational basis within the Great Lakes region. However, SHOP’s analyses of St. Lawrence River flows from Montréal to Trois Rivières, produced by CCMEP on an operational basis since 2013, are considered to be of interest to the readers since its analyses of flows on the St. Clair and Detroit Rivers will be available on an operational basis in the near future.
Products can be viewed on the web map service (WMS) of ECCC’s Recherche en Prévision Numérique (RPN-WMS). WMS is a protocol to serve georeferenced map images over the Internet (www.opengeospatial.org/standards/wms). Various viewers such as Adaguc (Royal Netherlands Meteorological Institute 2009) or the European Commission’s GLoFAS (European Commission 2016) can be used to display layers of interest. ECCC’s implementation of the WMS viewer Adaguc can be accessed online (http://collaboration.cmc.ec.gc.ca/cmc/hydro/www/adaguc).
Clicking on “cancel” in the login window redirects users to the registration page. The website’s username and password are sent to the e-mail address provided. The following http request retrieves a list of available fields: collaboration.cmc.ec.gc.ca/rpn-wms/?service=wms&version= 1.1.1&request=GetCapabilities.
In the near future, discharge observations from the U.S. Geological Survey (USGS) and Ontario Power Generation will be assimilated by WATROUTE and used for model verification. Additionally, forecast recharge from GEM-LAM (see Fig. 1 and the appendix) is expected to soon replace the climatological values currently used by WATROUTE. These upgrades should improve the riverine discharge and inflow forecasts. For the lakes, base lake levels will be corrected, removing the current lakewide bias. Uncertainties in surface currents will be assessed by comparison with surface drifters and observations at fixed locations. To simulate processes in the lakes’ coastal regions at a resolution finer than NEMO’s 2-km resolution, we plan to add a 2D hydraulic model (Matte 2014; Matte et al. 2017a,b) in these regions; the hydraulic model’s unstructured grid will provide more details of the nearshore circulation. Finally, a surface wave forecasting system for the Great Lakes is in development using WAVEWATCH III (Tolman 2009) and will be coupled to WCPS-GLS in future.
In the longer term, the techniques used to assimilate riverine discharge and lake ice observations will be upgraded. In the future, the value of an assimilated observed discharge will modify the values of WATROUTE’s simulated discharges upstream and/or forward in time. It is expected that the updated assimilation method will improve the quality of the forecast discharges for longer periods. Similarly, instead of the current method of inserting satellite-derived estimates of lake ice concentration and thickness into the CICE fields directly, the observed and simulated values will be blended horizontally to produce a smooth transition between purely observed and purely simulated values.
Additionally, to estimate water temperatures of larger rivers, we will add the river basin model (RBM; Yearsley 2009, 2012) riverine temperature model to WCPS. The skill of RBM, coupled to GEM and WATROUTE, has already been evaluated for the Kitimat River in British Columbia (Benyahya et al. 2017, manuscript submitted to Atmos.–Ocean). Elsewhere, riverine temperatures will be estimated from soil and surface-level atmospheric temperatures. Where available, observed riverine temperature will be assimilated. WATROUTE will pass the temperature of riverine inflows to NEMO.
Further planned upgrades for NEMO include the assimilation of observations of lake water level and temperature. Given the apparent propagation of biased surface water temperatures in the pseudoanalysis to the atmosphere (see “Long-term performance” subsection), the assimilation of water temperatures by NEMO–CICE will optimally be linked with the 2-m atmospheric temperatures analyzed by RDPS (see appendix). Similarly, the assimilation of observed lake water levels, which will include a representation of the glacial isostatic adjustment, will optimally simultaneously consider analyzed riverine inflows. Additionally, we plan to feed riverine discharges back into the process generating the soil moisture analysis. Thus, if discharges forced by the simulated surface runoff and drainage are overestimated, analyzed soil moisture will be reduced.
In the future, the 2.5-km HRDPS (Milbrandt et al. 2016) will provide boundary conditions for WCPS-GLS. The system’s GEM-LAM component will also be updated to resemble HRDPS. At this point, the 2.5-km horizontal resolution of the atmosphere will almost perfectly match the 2.0-km resolution of NEMO–CICE.
Finally, to improve the quality of WCPS as a whole, a statistical postprocessing system may be incorporated. Such a system, by reducing bias and possibly error in variables transferred between component models, could improve the quality of operational forecasts. Simultaneously, the diagnosed biases would inform researchers of weaknesses in the component models.
WCPS-GLS is the first operational modeling system of the complete water cycle at a short to medium range on an operational basis anywhere. WCPS-GLS generally produces reliable results throughout the 84-h forecast. As a rule, products from the atmospheric and lake water models are more mature and reliable than those of the river-routing model. Unfortunately, errors in initializing surface water temperatures appear to propagate into the atmosphere. Similarly, errors in precipitation appear to propagate into riverine discharges and inflows and subsequently into forecast lake water levels. The assimilation of riverine discharge observations limits the propagation of the precipitation errors.
The products generated (see “Products generated by WCPS-GLS” section) provide a new capacity at ECCC and have already generated considerable interest. Flood forecasters, generators of hydroelectricity, and educators regularly check the riverine discharge forecasts. The lake ice forecasts are used by the Canadian Ice Service to address gaps in coverage by synthetic aperture radar observations from satellite. ECCC’s National Environmental Emergencies Centre consults both surface currents and winds when responding to environmental emergencies such as oil spills. Similarly, the Canadian Coast Guard can use the simulated surface water currents for search-and-rescue efforts. High-water-level warnings for eastern Lake Erie are issued by the Ontario Storm Prediction Centre in advance of significant surge or seiche action that could lead to property damage. The advance warning enables the authorities to mitigate property damage where possible. Conservationists and ecologists will likely be interested in the forecasts of both lake water levels and temperatures. Finally, the interactive coupling of the atmosphere and ice–ocean models reduces errors in forecasts of near-surface conditions and is under subjective evaluation by weather forecasters.
D. Durnford, V. Fortin, and G. C. Smith are co–first authors. We are deeply grateful for funding in support of this project provided by the International Joint Commission (IJC), Public Safety Canada’s Search and Rescue New Initiatives Fund (SAR-NIF), and the National Science Foundation. The IJC funded the International Upper Great Lakes Study. We have also greatly benefited from many discussions with the Canadian and U.S. members of the Coordinating Committee on Great Lakes Basic Hydraulic and Hydrologic data and its Hydraulic, Hydrology, and Vertical Levels subcommittees. The breadth of knowledge encompassed by this group and their willingness to collaborate has greatly facilitated our work. Our thanks to the reviewers, whose comments deepened the manuscript.
APPENDIX: A DETAILED DESCRIPTION OF THE MODELS FORMING WCPS AND THEIR LINKAGES.
WCPS is a chain of interconnected numerical models. The first model in the chain is the Global Deterministic Prediction System v5.1 (GDPS; Qaddouri et al. 2015). GDPS is a numerical weather prediction system based on the GEM 3D atmospheric model (Côté et al. 1998a,b; Girard et al. 2014). GEM is ECCC’s operational weather forecasting model. It represents processes in the atmosphere. Processes on the surface and in the soil as well as exchanges of momentum, heat, and moisture between the atmosphere and the surface are represented by the land surface model ISBA (Bélair et al. 2003a,b). GDPS runs on a global domain with a 25-km horizontal resolution, 80 vertical levels, and a 12-min time step.
The second model in the WCPS chain is version 5.0 of RDPS (Caron et al. 2016). RDPS is a limited-area version of GEM over North America. RDPS uses a 10-km horizontal resolution, 80 vertical levels, and a 5-min time step. GDPS updates lateral boundary conditions hourly. The lateral boundary conditions inform RDPS on the state of the atmosphere on all sides of its domain at all levels.
Both GDPS and RDPS analyze the state of the atmosphere at the start of the forecast by assimilating observations. Numerous types of observations from multiple satellites are assimilated. Radiosonde, aircraft, ship, and buoy data are also ingested. These observations provide information about winds, temperatures, water vapor, clouds, rain, and snowpack. In all, some 2.2 and 12 million observations are assimilated each day, respectively, by RDPS (M. Reszka 2016, personal communication) and GDPS (J. Morneau and J. St.-James 2016, personal communication). An analysis of the state of the land surface and the soil column is generated by a sequential assimilation scheme (Bélair et al. 2003b). The land surface model ISBA uses this analysis as the starting point for the forecast.
The third model in the WCPS chain is a second limited-area version of GEM (GEM-LAM). GEM-LAM’s domain covers the Great Lakes–St. Lawrence River watershed. Like RDPS, GEM-LAM uses a 10-km horizontal resolution and 80 vertical levels. However, GEM-LAM’s time step is shorter at 3.75 min to permit two GEM-LAM time steps per 7.5-min NEMO–CICE time step (see below). GEM-LAM starts its forecasts from the analysis of the state of the atmosphere produced by RDPS. RDPS also provides hourly lateral boundary conditions for GEM-LAM throughout the forecast.
GEM-LAM and RDPS are very similar applications of GEM. The primary difference between these two applications is their treatment of surface water temperature and ice cover on the Great Lakes. RDPS relies on a 10-km analysis of surface water temperature and lake ice cover while GEM-LAM relies on a 2-km pseudoanalysis provided by the lake and ice model (see below). Additionally, whereas the surface water temperature and ice cover remain unchanged throughout the RDPS forecast, they evolve during the GEM-LAM forecast.
GEM-LAM is closely linked with the three remaining models in the WCPS chain of models. These are the 3D hydrodynamic lake/ocean model NEMO (Madec et al. 1998, 2008), the 2D dynamic and thermodynamic CICE (Hunke and Lipscomb 2010), and the 2D hydrological river-routing model WATROUTE (Kouwen 2010). These four models together form a closed circle of models. They simulate the complete water cycle, albeit in varying levels of detail for different aspects of the cycle. To facilitate the exchange of information between the four models, GEM-LAM, NEMO, and CICE are launched each day at 0000 and 1200 UTC while WATROUTE is launched at 0600 and 1800 UTC.
The NEMO–CICE ice–ocean model (Madec et al. 2008; Hunke and Lipscomb 2010; Dupont et al. 2015; G. C. Smith et al. 2016) represents processes occurring within and at the surface of the Great Lakes [Dupont et al. 2012; however, with a different sea ice component at that time, namely, Louvain-la-Neuve Sea Ice Model (LIM2)]. Currents, water temperatures, and marine ice evolve dynamically. NEMO coupled with a marine ice model has been run by CCMEP in the Gulf of St. Lawrence [initially with Modèle du Golfe du St. Laurent (MoGSL); Saucier et al. 2003; Smith et al. 2013; replaced in 2014 by CICE] on an operational basis in real time since 2011. NEMO and CICE are set up on the Great Lakes at a 2-km horizontal resolution. NEMO represents 35 vertical levels. NEMO and CICE both use a 7.5-min time step. Thus, GEM-LAM processes two time steps for every time step of NEMO and CICE.
Initial conditions for NEMO–CICE are generated by assimilating observations of the concentration and thickness of marine ice derived from RadarSat satellite images. The observed values are inserted directly into the model. Although this data assimilation method is unsophisticated and no observations of either water temperature or salinity are ingested, the quality of the marine ice and surface water temperature analyses matches or surpasses the operational CCMEP analyses of these variables (Smith et al. 2013). Furthermore, this procedure provides initial conditions on NEMO’s 2-km grid, which is at a higher resolution than the 10-km operational CCMEP analyses.
Atmospheric forcing fields at the surface of the water, such as winds, temperature, moisture, and solar radiation (Table A1), are provided by GEM-LAM to NEMO–CICE every NEMO–CICE time step. During the following NEMO–CICE time step, or 7.5 min later, these forcing fields drive currents in the surface waters and modify the surface water temperature and the ice cover. NEMO–CICE calculates the fluxes of momentum, heat, and moisture between the atmosphere and the surface since the horizontal resolution of NEMO–CICE is higher than that of GEM-LAM. Moreover, NEMO–CICE is better informed on the characteristics of the marine ice and its snow cover. After calculating fluxes for the different categories of water, ice, and snow, NEMO–CICE aggregates the fluxes and passes them to GEM-LAM. The data are transmitted between GEM-LAM and NEMO–CICE via CCMEP’s server Globally Organized System for Simulation Information Passing (GOSSIP; Smith et al. 2018). GOSSIP merely transmits the data. GEM-LAM and NEMO perform all needed interpolation and regridding of data.
Thus, the P − E forcing for water levels is determined interactively by GEM-LAM and NEMO–CICE based on respective surface properties (winds, temperature, humidity) from the two models. GEM-LAM determines the precipitation and NEMO–CICE calculates the evaporation. Changes in total water levels are then determined by the resulting P − E plus river runoff (P − E + R) from the coupled model.
As mentioned, the exchange of data between GEM-LAM and NEMO–CICE occurs at every common time step (7.5 min). This enables WCPS to capture the effects of rapid changes (e.g., the rapid ice freeze-up event shown in Fig. 9 and discussed in the “Cold season cold snap” subsection). While reducing this coupling frequency to 1 h or longer (as used in many other coupled systems) may not have a large impact on the overall forecast skill, the additional computation cost of the added exchanges is seen as small compared to the potential benefits for extreme events.
WATROUTE provides boundary conditions to NEMO–CICE in the form of time-varying riverine discharges (Table A1). The average simulated outflow from each river over the 48 h ending 6 h prior to the start of the analysis-generating process is calculated by NEMO–CICE at the start of that process and is ingested throughout the process. The same mean discharges are ingested during the following 84-h forecast.
WATROUTE simulates the flow of water (Kouwen 2010; Deacu et al. 2012; Gaborit et al. 2017) through a river network. WATROUTE’s domain covers the Great Lakes and St. Lawrence River watershed to Tadoussac at a ∼1-km resolution. The basic river network is derived from the USGS hydrological data and maps based on shuttle elevation derivatives at multiple scales (HydroSHEDS) dataset (Lehner et al. 2008; Gaborit et al. 2017). It is then modified locally to promote agreement between observed and simulated subwatersheds. The initial 7.5-min model time step diminishes as needed to permit the routing of large pulses of streamflow generated by significant rainfall or snowmelt.
Mathematically, WATROUTE solves the continuity equation in each grid cell, starting with the cell at the highest elevation and progressing downstream: a grid cell’s change in river channel storage over a time step equals the difference of its inflow and outflow during the time step. Outflow calculated for a cell upstream constitutes inflow to the grid cell downstream. ISBA, the land surface scheme in GEM-LAM and RDPS, combines with WATROUTE to constitute a 1D distributed hydrological model.
GEM-LAM provides WATROUTE with hourly accumulations of surface runoff (Table A1). The hourly accumulations are ingested by WATROUTE at the start of its forecast and applied at the appropriate time. WATROUTE’s data assimilation cycle uses runoff purely from GEM-LAM. However, since GEM-LAM’s forecast is launched 6 h before WATROUTE’s, WATROUTE’s forecast uses runoff from GEM-LAM for all but the last 6 h, at which point RDPS provides the runoff. Surface runoff is calculated in ISBA following Eq. (1):
where P is precipitation reaching the ground directly, R is runoff from vegetation, M is snowmelt, E is evaporation, and F is infiltration into soil. The variability of surface runoff causes frequent rapid changes in riverine streamflow.
WATROUTE also requires hourly accumulations of recharge, or the water that percolates through the soil column and enters subterranean reservoirs. Recharge is released into the river network at a slowly evolving rate, simulating the groundwater contribution to river flows. Currently, climatological (2004–09; Deacu et al. 2012) recharge produced by Global Environmental Multiscale model-Hydro (GEM-Hydro) (Gaborit et al. 2017) is used owing to its overestimation by the sequential surface data assimilation system (see above). This decision to ignore recharge from the GEM model will be revisited in 2017 when the sequential surface data assimilation system will be replaced by the Canadian Land Data Assimilation System (CaLDAS; Carrera et al. 2015). CaLDAS relies on an ensemble Kalman filter procedure. It will include the land surface model Soil, Vegetation, and Snow (SVS; Alavi et al. 2016; Husain et al. 2016). CaLDAS assimilates remotely sensed information on soil moisture in addition to surface observations of temperature and dewpoint.
Initial conditions for WATROUTE’s forecasts are generated by assimilating hourly observations of discharge. The observed values are inserted directly into the model. Currently, observations from 213 ECCC stations are assimilated in the Great Lakes area. Downstream of major dams where the flow can fluctuate unpredictably (see Table A2), the most recent observed discharge value persists throughout the 84-h forecast. Specifying persistence of the most recent observed discharge on the Nipigon River is the method currently used to represent the diversion of water from Lake Ogoki and its watershed into Lake Superior’s watershed.
Flow through the interconnecting channels between the lakes is represented by CGLRRM (Tolson 2009). CGLRRM uses a stage–discharge curve to calculate the volume of water discharged from each lake. The calculation considers the mean water level of either the lake discharging the water (Lakes Superior, Erie, Ontario) or of the lake at both ends of the connecting channel (Lakes Michigan–Huron and St. Clair). CGLRRM is embedded within NEMO–CICE.
A supplement to this article is available online (10.1175/BAMS-D-16-0155.2)