Introduction

Groundwater, as the largest freshwater storage component of the hydrological systems, is an essential resource necessary to sustain agricultural, industrial, and domestic activities in Canada. About one third of Canadians depends on groundwater for drinking water and up to 80% of Canada’s rural population uses groundwater for its entire water supply (Council of Canadian Academies 2009). Groundwater also plays an important role in sustaining ecosystems and the water cycle’s response to climate change at regional and global scales. Increased human withdrawals of groundwater or changes in climate have resulted in growing pressure on groundwater resources, which poses a serious threat to water security and potentially causes a decline in agricultural productivity and energy production (Frappart and Ramillien 2018). Monitoring and understanding groundwater storage changes is thus critical for maintaining sustainable economic development and healthy ecosystems, and for better understanding the hydrological cycles and climate change (Chen et al. 2016).

Monitoring and quantifying groundwater is difficult because it deals with water in a complex subsurface environment. Well monitoring is the traditional approach for estimating groundwater storage but often needs knowledge of the aquifer structure and some critical parameters such as transmissivity and storativity (or specific yield), which are difficult to obtain. In addition, well monitoring is not spatially continuous and has a high cost for a large region. There are very limited wells in some areas, especially in remote and harsh environments such as cold regions in Canada due to difficulties of access and monitoring.

Satellite remote sensing is increasingly being used in hydrological studies for its large spatial coverage and cost-effectiveness, and its ability to provide data in a timely manner. The Gravity Recovery and Climate Experiment (GRACE) satellite mission, launched in 2002, provides a unique approach to detect gravity changes globally due to redistribution of mass in the Earth system (Tapley et al. 2004). By separating the contributions to mass changes, GRACE-observed gravity changes can be used to derive terrestrial water storage (TWS) change, which includes the changes of groundwater storage (Wground), soil-water content (Wsoil), snow-water equivalent (Wsnow), and surface-water storage (Wsurf). Thus, GRACE provides a practical way to estimate Wground when Wsoil, Wsnow and Wsurf can be estimated.

The GRACE TWS data have been used to quantify variations and long-term trends of Wground at regional and global scales in various studies (e.g., Rodell et al. 2007, 2009; Swenson et al. 2006; Famiglietti et al. 2011; Richey 2014; Huang et al. 2016; Bahanja et al. 2018; Shamsudduha and Taylor 2020; Opie et al. 2020). As GRACE cannot separate the different components of TWS, the Wsoil, Wsnow and Wsurf estimated from other independent approaches such as land-surface models (LSMs), were usually used to separate Wground from GRACE TWS. Famiglietti et al. (2011) and Scanlon et al. (2012) analyzed long-term Wground change in the California Central Valley (USA) region by combining TWS from GRACE, Wsoil from GLDAS (Global Land Data Assimilation System) LSMs, Wsnow from the Snow Data Assimilation System (SNODAS) datasets, and Wsurf from in-situ surface-water reservoir measurements. Thomas and Famiglietti (2019) extended the analysis of groundwater change in USA by using Wsoil and Wsurf (surrogated by surface runoff) from GLDAS LSMs, Wsnow from SNODAS, and GRACE TWS, and identified climate-induced groundwater depletion.

While the GRACE satellite mission provides an innovative tool for groundwater estimation through its integration with other independent in-situ and model data sources, so far studies and information on the groundwater climatology for Canada’s landmass are still very limited. On the other hand, the vast majority of the existing studies, as mentioned before, have relied on the land-surface models (LSMs) included in GLDAS. Applying uncalibrated and unvalidated global-scale models to a region involves high uncertainty. This is particularly so for Canada where the hydrological processes are much more complex due to the cold-region processes dominating the landmass such as snow accumulation/melt and soil freeze/thaw. Indeed, model comparison studies have revealed large biases and uncertainties in simulating the water cycle over Canada’s landmass (Wang et al. 2015a; Xia et al. 2015; Mortimer et al. 2020). EALCO (Ecological Assimilation of Land and Climate Observations model) has been developed by Natural Resources Canada with a specific focus on cold-region mechanisms for land-surface radiation transfer, energy balance, water dynamics, and carbon and nitrogen biogeochemical cycles which play various roles in the physiological and ecohydrological processes. In particular, the freeze/thaw processes for soils, lakes and snow cover, and the remote-sensing-based vegetation and land-surface dynamics (e.g., albedo, surface-water cover), have demonstrated robustness in a number of international LSM intercomparison studies (Zhang et al. 2008; Widlowski et al. 2011; Wang et al. 2014, 2015a, b). This gave confidence in using EALCO results for this study on the Canadian landmass.

While EALCO includes lake/river surface evaporation simulations, it does not include the river routing process so it is difficult to apply its results (e.g., snowmelt) directly in this study. In many published studies, surface water is ignored (Huang et al. 2016; Chen et al. 2014; Rodell et al. 2007), while in some others LSM-simulated surface runoff is used as a proxy for surface-water storage change (Shamsudduha and Taylor 2020; Thomas and Famiglietti 2019). In this study, EALCO-simulated surface runoff is used as a proxy for surface-water storage.

The objective of this study is to assess the seasonal variations and long-term trends of the Wground for Canada’s landmass for the period 2003–2016 using the TWS from GRACE monthly spherical harmonic (SH) solutions and the Wsoil, Wsnow and Wsurf values from the Canadian LSM of EALCO. This study calculates a suite of parameters representing the spatial and temporal variations of Wground, including the occurrence time and the interannual variability of the maximum and minimum Wground in a year, the range of seasonal variation of Wground and its interannual variability, and the long-term trend of Wground over 2003–2016. The results are presented in national maps and aggregated into the hydrological units of the major drainage areas (MDAs) of Canada. The GRACE SH-derived Wground is also compared to those derived from two GRACE monthly Mascon solutions and the groundwater well measurements.

Data and methods

TWS variations in most cases are mainly contributed by the changes in Wground, Wsnow, and Wsurf. In this study, Wground is computed as

$$ {W}_{\mathrm{ground}}=\mathrm{TWS}-{W}_{\mathrm{soil}}-{W}_{\mathrm{snow}}-{W}_{\mathrm{surf}} $$
(1)

The TWS used in this study includes the GRACE monthly SH solutions (Landerer 2020), containing monthly TWS changes processed by three different data processing centers: the Center for Space Research (CSR; University of Texas, USA), GeoForschungsZentrum (GFZ; Potsdam, Germany), and Jet Propulsion Laboratory (JPL, USA). The TWS data were downloaded from the GRACE Tellus website (JPL 2020). All three TWS datasets were processed from the latest release RL-06 V03 with 1° × 1° global grids. The datasets were calibrated with standard corrections including a glacial isostatic adjustment (GIA). Post-processing filters including a de-stripping filter and a 300-km-wide Gaussian smoothing filter were applied to reduce correlated errors (Landerer and Swenson 2012). Scaling coefficients for the global land grids were applied to restore much of the energy removed by the postprocessing filtering. A detailed description of the data processing is available in Landerer and Swenson (2012). The monthly TWS data are anomalies to the baseline average over the study period of January 2003 to December 2016. Note that the baseline in the original datasets is the average over the period of January 2004 to December 2009. Because the differences among the three TWS datasets were found to be small over the study region (Wang and Li 2016), the averages of the three datasets were used in this analysis. For the result comparisons, two other TWS products derived from the GRACE monthly Mascon solutions (GRCTellus. JPL RL06M.MSCNv02 and CSR RL06M.MSCNv02) were also used in this study. More detail about the Mascon solution data can be found in Wiese et al. (2018). There are 13 months with missing GRACE TWS data during the study period, which include June 2003, January and June 2011, May 2012, March and September 2013, February 2014, May, October and November 2015, and April, September and October 2016.

To estimate GRACE-derived Wground in Eq. (1), the simulated soil moisture Wsoil, snow water equivalent Wsnow, and surface runoff (as a proxy for surface-water storage Wsurf) from the LSM EALCO were used. EALCO is driven by climate forcing at a half-hourly time step and 5-km spatial resolution. The monthly Wsoil data include water contained in the soil column of up to 4.0 m depth including seven soil layers (0–10, 10–20, 20–40, 40–80, 80–140, 140–240, 240–400 cm). When the water table is above the bottom of the seventh layer (i.e., 4.0 m), the soil layers below the water table will be exposed to groundwater and the actual number of soil (unsaturated) layers will be determined by the actual depth of the water table. The monthly Wsoil, Wsnow and Wsurf for the period 2003–2016 are calculated from the half-hourly EALCO Wsoil, Wsnow and Wsurf values based on the actual days used for each GRACE monthly solution. Since the GRACE TWS data are anomalies relative to the baseline average over the study period, the Wsoil, Wsnow and Wsurf anomalies were calculated by subtracting the time-mean baseline over the same period to make EALCO and GRACE data comparable.

Eight parameters are calculated from the derived monthly Wground data to characterize the spatiotemporal variations of Wground (Table 1). They include: the long-term trend (τ) represented by the Thiel-Sen linear regression slope of monthly Wground time series for 2003–2016; the months with maximum/minimum Wground in a year (Mmax/Mmin) and their interannual variations (σmax/σmin); and the range of seasonal Wground variation (ΔWground) and its interannual variations (represented by the standard deviation σΔground and the coefficient of variance COV of ΔWground). For the calculation of σmax and σmin, if the difference of Mmax (or Mmin) among different years is larger than 6, the difference is subtracted by 12.

Table 1 Parameters used for characterizing Wground
Full size table

The study area is shown in Fig. 1. It covers most of the Canadian landmass except areas (the grey area in Fig. 1) having no groundwater recharge due to permafrost as simulated by the EALCO model—Fig. S1 of the electronic supplementary material (ESM)—or snow/ice covering probability larger than 50% in the warm season (Fig. S2 of the ESM, Trishchenko and Ungureanu 2021). The study region includes eight MDAs of Canada. The GRACE TWS data are provided in 1° × 1° grids, but the actual resolution is much coarser (>350 km × 350 km). The Rocky Mountains with permanent snow/ice occupies a large portion of the Pacific MDA and the remaining regions (e.g., west coast) barely fits this coarse resolution. The GRACE grid often includes fractional areas of permanent snow/ice; therefore, the entire Pacific MDA is excluded in the analysis. The study also excludes the Mississippi River MDA due to its small area (~27,000 km2). All the results are presented in national maps under the Lambert Conformal Conic (LCC) projection and are also aggregated into the MDAs. The MDAs provide hydrological units that are frequently used for data collection and compilation, and for spatial analysis of environmental, economic and social statistics (Statistics Canada 2003).

Fig. 1
figure 1

Map of the study region, covering eight major drainage areas (MDAs, polygons shown in colors) over Canada’s landmass. Note that the Nelson River MDA excluded a small portion in the west, the Western and Northern Hudson Bay MDA excluded the portion in high Arctic, the Great Slave Lake MDA excluded a small portion in the north and the west, and the Arctic MDA only included its southern part. The grey areas, having permafrost or permanent snow/glacier probability larger than 50%, or areas smaller than the GRACE footprint, are excluded in this study

Full size image

Daily groundwater level data from groundwater observation wells available in two GRACE grids in the Nelson River MDA and the St. Lawrence MDA are used for Wground comparisons. The well data for groundwater level (Wlevel) can be downloaded from the Groundwater Information Network (GIN) portal (GIN 2021), which connects Alberta’s groundwater observation well network (GOWN), and the Quebec groundwater monitoring network (QGMN 2021). For the study period, there are multiple wells with either partial or complete data for both comparison grids. For each grid, the wells that have a long overlap monitoring period with the study period are selected. These include four wells located in the Nelson River MDA and six wells in the St. Lawrence MDA (Fig. 1). Monthly well data are obtained by averaging the daily data for each month, and subtracting their respective time-mean baseline. It was noted that the wells in the Nelson River MDA and the St. Lawrence MDA have observations covering the time periods of 2003–2016 and 2005–2016, separately. Finally, the mean Wlevel change of the wells in each comparison grid is calculated to represent the Wlevel change for that specific grid.

Results

Seasonal variations of W ground

Figure 2 shows the average monthly Wground over 2003–2016 and provides a general overview of the spatial and seasonal variations of Wground over the study area. In January, Wground presents small spatial variation across the study area with a value around the annual average (0 mm) over the study period. From February to June, the Wground in the western region (including the Western and Northern Hudson Bay, the Great Slave Lake, and the Arctic MDAs) and the Southern Hudson Bay area slightly increases till early spring, after which the Wground rapidly increases and reaches its yearly highest value in June due to high groundwater recharge after snowmelt and soil thaw. For example, the Southern Hudson Bay MDA has the highest Wground (~75 mm above the annual average) in June and the Western and Northern Hudson Bay MDA has the highest Wground (~ 50 mm above the annual average) in June. In contrast, the Wground in east Canada gradually decreases during this period and mostly reaches its lowest value until snowmelt and soil thaw start. South Canada, from southern Ontario to east Prairie and including the Maritime Provinces, reaches the lowest Wground in March or April—for example, the Maritime Provinces MDA and the south portion of the St. Lawrence MDA (mainly in south Ontario), reaches the lowest Wground of ~100 mm below the annual average in April. After that, the Wground rapidly increases and reaches a high value around 75 mm above the annual average in June. The long winter in the north delays the snowmelt. As a result, the north part of east Canada has 1–2 months delay to reach the lowest Wground—for example, the north portion of the St. Lawrence MDA and the Northern Quebec and Labrador MDA reach the lowest Wground (~75 mm below the annual average) mostly 1 month later in May and the north Labrador portion has the lowest Wground (over 100 mm below the annual average) 2 months later in June. In June, a large portion of the study area shows a high Wground except that the Northern Quebec and Labrador MDA still presents a low Wground with the value of around 50 mm below the annual average.

Fig. 2
figure 2

The average monthly Wground over 2003–2016, showing its spatial and seasonal variations over the landmass

Full size image

In summer, the western region and the southwestern Hudson Bay area have high evapotranspiration. As a result, they lead to relatively low groundwater recharge, resulting in net water loss in these regions. Most of these regions experience rapidly decreasing Wground in summer, after which the Wground reaches its lowest value before winter (Fig. 2)—for example, the Arctic MDA arrives at the lowest Wground (~30 mm below the annual average) mostly in October/November, while southwestern Hudson Bay indicates that it took 1 month later to reach the lowest Wground (~60 mm below the annual average). In contrast, the Northern Quebec and Labrador MDA in east Canada mainly shows a rapid increase of Wground in summer. These regions have high precipitation in fall, resulting in another high groundwater recharge before winter. As a result, the Wground in most of these regions reaches the annual peak (e.g., over 100 mm above the annual average in north Labrador) in October or November, after which the Wground rapidly decreases in January. Like in the winter and spring, south Canada also presents different seasonal patterns in Wground variations during summer and fall, whereby it has substantial precipitation after spring, resulting in rapid groundwater recharge. The highest Wground (e.g., over 100 mm above the annual average in the southern Ontario) occurs in July or August, after which the Wground rapidly drops in October and then gradually decrease throughout fall.

It is observed from Fig. 2 that the study area overall presents low Wground in early spring and high Wground in early summer. The months having the lowest/highest Wground vary by regions. Figure 3a,b shows the spatial distribution of the months having maximum/minimum Wground (Mmax/ Mmin), for 1 year for the study period. As discussed earlier, the Mmax/Mmin presents large spatial differences over the study area. Generally, the study region shows Mmax mainly in late spring and early summer, except for the Northern Quebec and Labrador MDA which has Mmax mainly in October or November. In contrast to Mmax, the Mmin generally appears in early spring (April or May) for the study region, except for south Ontario and the south part of the Nelson River MDA, which have Mmin in March and the north portion of the Southwestern Hudson Bay MDA, which has Mmin in fall (October or November).

Fig. 3
figure 3

The months with (a) maximum Wground (Mmax) and (b) minimum Wground (Mmin), and their interannual variations (c–d) represented by the standard deviations of Mmax and Mmin, respectively

Full size image

Figure 3c,d shows the interannual variations (σmax and σmin) of Mmax and Mmin, respectively, among the 14 years. Generally, the σmax and σmin both have values less than 2 months for most of the study area. As shown in Fig. 3c, the σmax shows small spatial variations without pronounced patterns. In east Canada, the southern Ontario and the southern part of the Northern Quebec and Labrador MDA have σmax of less than 1 month, whereas the remaining areas (e.g. far-north Quebec and the corridor along the St. Lawrence River) have σmax of ~2 months. For the western region, the σmax is generally larger than the east region. The σmin presents small values of mainly less than 1 month in east Canada and relatively large values of up to 3 months in the western region. In east Canada, the σmin presents very small spatial variations with abnormal values larger than 2 months in two small areas: southern Ontario and far north Quebec. For the western region, the σmin in the western part is generally 1 month larger than the eastern part.

Figure 4a shows the maximum range of seasonal Wground variation (ΔWground), averaged over the study period. The ΔWground presents large spatial differences with continuous increase from west to east in the study region. The large portion in east Canada has ΔWground above 250 mm, with some areas (e.g. Newfoundland island) having values of up to 400 mm. The large portion of the central area has ΔWground around 120 mm. For the western region, the ΔWground is low with the smallest ΔWground of <30 mm in the Prairie region.

Fig. 4
figure 4

(a) Range of seasonal variation (ΔWground) of Wground, (b) its interannual variation represented by the standard deviation (σΔWground), and (c) the coefficient of variance (COV)

Full size image

The interannual variation of ΔWground among the 14 years represented by the standard deviation (σΔWground), as shown in Fig. 4b, presents a similar spatial pattern to ΔWground (Fig. 4a). However, in relative terms, the interannual variations of ΔWground represented by the coefficient of variance (COV) shows small spatial differences and the COV is mostly less than 30% except in some small areas scattered over the study region such as Maritime Provinces MDA which has high values of 35–50%.

Long-term trends of W ground

Figure 5 shows the trends of Wground over the 14-year study period. In general, the trends present large spatial variability across the study region, mainly varying from increasing trends (i.e., water gain) of >10 mm/year to decreasing trends (i.e., water loss) of < −15 mm/year. Significant increasing trends are mainly found in east Canada (except for far north Quebec) and a central zone in northeastern Manitoba. The increasing trends of Wground in east Canada are likely associated with the increase of precipitation during the study period (Li et al. 2016). The St. Lawrence MDA experienced severe drought in 2001 and 2002 and the St. Lawrence River water level plunged to its lowest point in more than 30 years, which led to a very low Wground at the beginning of this study period (Wikipedia 2021). The consistent water gain over the 14 years is also a result of the recovery of Wground after the severe drought. Significant decreasing trends are observed mainly in the western region. The highest decreasing trend of Wground appears in some areas of the Great Slave Lake MDA. The large decreasing trends of Wground in the western regions are likely due to the decrease of precipitation (Li et al. 2016).

Fig. 5
figure 5

The trends (mm/year) of Wground over the period 2003–2016

Full size image

W ground for major drainage areas (MDAs)

The variations of Wground for the MDAs are summarized in Table 2 and Fig. 6. The Mmax for the Arctic and the Southwestern Hudson Bay MDAs, the Western and Northern Hudson Bay and the Great Slave Lake MDAs occurs in June, followed by the St. Lawrence MDA in July. The Maritime Provinces and the Nelson River MDAs have Mmax in August. The Northern Quebec and Labrador MDA demonstrates late Mmax in October. Most of the MDAs have Mmin in April except the Northern Quebec and Labrador MDA and the Arctic MDA, having Mmin in May and October, respectfully.

Table 2 The characteristics of groundwater storage (Wground) for major drainage areas (MDAs)
Full size table
Fig. 6
figure 6

Seasonal variations of Wground for the MDAs in the study area. The black lines represent the average monthly Wground over 2003–2016. The grey areas represent the Wground variation ranges in 2003–2016. The error bars show the spatial variability of the Wground in the MDAs, represented by its standard deviations

Full size image

The ΔWground varied between 54.8 and 280.5 mm among the MDAs, with the highest and lowest ΔWground in the Maritime Provinces MDA and the Western and Northern Hudson Bay MDA, respectively. The three MDAs in east Canada and the Southwestern Hudson Bay MDA demonstrate ΔWground larger than 100 mm. All the four MDAs in the western region have small ΔWground values around 50–60 mm, which are less than quarter of that for the Maritime Provinces MDA. For the entire landmass, the Wground reaches the high peak in July and the low peak in April (Table 2), with a ΔWground of 117.7 mm or 699.9 km3. It is worth noting that there is large spatial variability of Wground in each MDA especially in the Maritime Provinces MDA (Figs. 2 and 6).

The long-term trends of Wground for the MDAs (Table 2) show that the 4 MDAs in east Canada (Maritime Provinces, St. Lawrence, Northern Quebec and Labrador, and Southwestern Hudson Bay) have increasing trends from 2003 to 2016. The Maritime Provinces MDA has the highest Wground increase of 6.25 mm/year, followed by the St. Lawrence MDA having an increasing trend of 4.97 mm/year. The Northern Quebec and Labrador and the Southwestern Hudson Bay MDAs shows 1.80 and 2.45 mm/year of Wground increase, respectively. Among four MDAs in the western regions, the Great Slave Lake MDA has the largest Wground decrease of 6.67 mm/year. The Western and Northern Hudson Bay MDA and the Arctic MDA show Wground decrease of 3.93 and 4.76 mm/year, separately. For the Nelson River MDA, although the entire drainage area has a slight decreasing trend of −0.04 mm/year, it presents large spatial variability with increasing trends of up to 8 mm/year in some areas of western portion (mainly in Prairie pothole region) and decreasing trends of up to −15 mm/year in some areas of the eastern portion. For the entire study region, a slight increase of 0.002 or 0.01 km3/year of overall Wground demonstrates that there is no significant trend over the 14-year period of 2003–2016.

Comparisons of different GRACE solutions

The Wground discussed in the preceding analyses is based on the GRACE SH solutions. The results are compared with those derived from two GRACE Mascon solutions—JPL Mascon solution (JPL-MSCN) and CSR Mascon solution (CSR-MSCN)—in Fig. 7 and Table 3. Generally, the Wground derived from all three solutions presents highly consistent seasonal patterns and trends for all the MDAs, although the Maritime Provinces MDA and the Arctic MDA present relatively large differences among the three solutions. More specific comparisons, which include correlation coefficients (R) and root mean square errors (RMSE), are summarized in Table 3. The SH Wground is highly correlated to those derived from JPL Mascon solution (R range 0.75–0.97) and CSR Mascon solution (R range 0.83–0.98). The largest differences between SH and the two Mascon solutions occur in the Maritime Provinces, with RMSE of 49.4 mm for SH vs. JPL-MSCN and 28.9 mm for SH vs. CSR-MSCN, followed by the Arctic MDA having RMSE of 22.4 and 18.9 mm for the JPL-MSCN and the CSR-MSCN, respectively. The smallest differences between SH and the Mascon solutions occur in the Northern Quebec and Labrador and the Western and Northern Hudson Bay MDAs (Table 3). For the entire study region, the SH-derived Wground has a R value of 0.94 with JPL-MSCN and 0.98 with CSR-MSCN. The corresponding RMSEs are 10.5 and 5.58 mm, respectively.

Fig. 7
figure 7

Time series of monthly Wground derived from three GRACE solutions over the period 2003–2016: the GRACE spherical harmonic solution (SH), the Jet Propulsion Laboratory Mascon solution (JPL-MSCN) and the Center for Space Research Mascon solution (CSR-MSCN)

Full size image
Table 3 Comparisons of Wground derived from the spherical harmonic solution, the JPL Mascon solution, and the CSR Mascon solution
Full size table

Comparisons to groundwater level (W level) measurements

Due to the difficulty in acquiring groundwater storage observations, this study compares the SH-derived Wground to the groundwater level (Wlevel) measured in wells in this study area. Observations from a total of 10 wells, with 4 wells within a GRACE SH grid in the Nelson River MDA and 6 wells within a GRACE SH grid in the St. Lawrence MDA, are used (Fig. 1). The relationships (R values) between each two Wlevel measurements vary from 0.047 to 0.821 with a mean value of 0.46 for the grid in the Nelson River MDA, and from 0.3 to 0.94 with a mean value of 0.56 for the grid in the St. Lawrence MDA. The low correlations among the wells in each grid are partially due to the relatively large measurement uncertainties compared with the small Wlevel variations. Figure 8 compares the time series of GRACE-based Wground and Wlevel (well average) over the study period for the two GRACE SH grids. The Wground and Wlevel present similar trends in both grids. For the grid in the Nelson River MDA (Fig. 8a), both Wground and Wlevel show four subtrends during the study period, i.e., increasing trends for the period 2003–2006, decreasing trends during the period 2006 to 2009, increasing trends from 2009 to 2011, and decreasing trends in 2011–2016. For the grid in the St. Lawrence MDA (Fig. 8b), both Wground and Wlevel do not present notable long-term trends. The correlation coefficients between Wground and Wlevel are 0.3 for the Nelson River MDA, and 0.26 for the St. Lawrence MDA. It is observed that the Wground and Wlevel seasonal variations show phase differences, which partially result in the low correlation coefficients between Wground and Wlevel. The phase differences can be clearly seen from the monthly average values shown in Fig. 9. Generally, the Wground shows the seasonal pattern 1 month earlier than the Wlevel for the grid in Nelson River MDA. When moving the Wlevel data 1 month forward, the R increases to 0.37 from 0.3. For the grid in the St. Lawrence MDA, the Wground and Wlevel show different phase differences. After moving the Wground time series 1 month forward, the R between Wground and Wlevel increases to 0.51 from 0.26.

Fig. 8
figure 8

Variations of Wground and Wlevel for grids in a the Nelson River MDA and b the St. Lawrence MDA over the period 2003–2016

Full size image
Fig. 9
figure 9

Average monthly Wground and Wlevel over the study period for grids in a the Nelson River MDA and b the St. Lawrence MDA. The vertical bars show variability of Wlevel among the observation wells, represented by the standard deviations

Full size image

The range of seasonal variation for Wlevel and Wground is ~0.2 m (ΔWlevel) and ~40 mm (ΔWground), respectively, for the grid in the Nelson River MDA. The grid in the St. Lawrence MDA has ΔWlevel of ~1.0 m and ΔWground of ~110 mm. Both grids show large variabilities of Wlevel among the well measurements (standard deviation of ~0.4 m for the grid in the Nelson River MDA and ~0.6 m for the grids in the St. Lawrence MDA, see Fig. 9). The large differences can also be seen from the low correlation coefficients among the well measurements as discussed in the preceding. It is noted that the variability of Wlevel among the well measurements is twice that of the seasonal variation range (ΔWlevel) for the grid in the Nelson River MDA, and about half of the ΔWlevel for the grid in the St. Lawrence MDA, suggesting relatively large uncertainty in Wlevel measurements and small variations of the Wlevel signal. This poses large challenges for directly comparing GRACE-based groundwater estimations with in-situ well measurements, and it is a major reason for the low R values between Wground and Wlevel discussed in the preceding.

Discussion

In this study, Wground is estimated using the GRACE-observed TWS and LSM-based Wsoil, Wsnow, and Wsurf from the EALCO model. Errors or uncertainties in TWS, Wsoil, Wsnow, and Wsurf estimates would directly propagate to Wground. Mortimer et al. (2020) evaluated nine Northern Hemishphere gridded Wsnow products from three categories—reanalysis models, passive microwave remote sensing combined with surface observations, and stand-alone passive microwave retrievals. Evaluation against snow course measurements in Canada shows high RMSE and bias, and low correlation. The uncertainty for all products tends to increase with deeper snow. Results showed the errors (RMSE as a percentage of mean in-situ Wsnow) are larger than 50% for Canada. Xia et al. (2015) used soil-moisture observations from seven observational networks in USA with different biome and climate conditions to evaluate soil-moisture products simulated from four LSMs, including Noah, Mosaic, SAC (Sacramento soil moisture accounting), and VIC (Variable Infiltration Capacity model). The errors (RMSE as a percentage of mean in-situ soil moisture) show that strong seasonal variations varied by model and soil layer, having values larger than 20% for most of the observational networks. The EALCO model used in this study is developed in Canada with comprehensive algorithms for cold-region processes, and has undergone extensive calibration and validation for the Canadian landmass (e.g., Zhang et al. 2008; Widlowski et al. 2011; Wang et al. 2014, 2015a, b); therefore, the errors in the estimates for these water variables are expected to be smaller in this study than those discussed in the preceding.

GRACE-observed TWS errors, including errors from measurements, spatial and spectral leakages, post-processing, and GIA adjustment, would also bring uncertainties to the Wground estimates. The combined error from measurements and leakages depends on the drainage area/basin size (Zhang et al. 2016; Wang et al. 2014; Landerer and Swenson 2012). Generally, it would be within 15 mm at the MDA level in Canada. It is worth noting that the GIA uplift rates are mostly under 10 mm/year in Canada (Fatolazadeh and Goita 2021; Argus et al. 2021), which is much smaller than the magnitudes of seasonal variations of the TWS signal. Thus, the GIA impact on the water characterization is more likely on the long-term trend rather than the seasonal patterns.

Integrating all the errors from TWS, Wsoil, Wsnow and Wsurf, Frappart and Ramillien (2018) showed that the errors on trend estimates of Wground derived from GRACE observations and LSMs were around 10% or less in several areas/basins around the world (e.g., 9.5% in the California Central Valley (USA), Scanlon et al. 2012; 13.6% in the north of China, Feng et al. 2013; 7.14% in Colorado Basin (USA), Castle et al. 2014; and 2.2% in the Tigris and the Euphrates Basin (Asia), Voss et al. 2013). In cold regions like Canada, snow and ice are important features of the landscape. Wsnow is a primary contribution to TWS in winter time. The large uncertainty of LSM-based estimates of Wsnow might significantly impact accuracy of the GRACE-derived Wground. Directly quantifying the uncertainty of Wground is still difficult due to the lack of independent data and it needs to be addressed in future studies.

Summary

Using the GRACE-based TWS and the EALCO modelled Wsoil, Wsnow, and Wsurf, the seasonal variations and trends of the Wground over a large portion of Canada’s landmass for the period of 2003–2016 were assessed. It is observed that the overall Wground has maximum/minimum values in July/April with a seasonal variation range (ΔWground) of 118 mm (or 700 km3) for the landmass. East Canada has relatively large ΔWground mainly in the range of 200–250 mm with the largest values of up to 400 mm in the Newfoundland. The west region has relatively small ΔWground, mostly around 125 mm, with the smallest values (lower than 50 mm) appearing in the Prairie area. The months having the maximum and the minimum groundwater storage (Mmax and Mmin) vary by regions. The western region and the Southwestern Hudson Bay major drainage area (MDA) have the Mmax and Mmin mainly in spring and fall, respectively, while, in contrast, east Canada has the Mmax and Mmin mostly in fall and spring, respectively. The trend analyses over the study period indicate that significant increasing trends or water gains of up to 10 mm/year are observed in east Canada (except for far north Quebec). Significant decreasing trends or water loss are mainly observed in the central-west region with values of above −10 mm/year. The increasing trend largely offsets the decreasing trend in the study area, and the overall Wground for the entire landmass does not show a significant trend over 2003–2016. The comparison of Wground from the GRACE spherical harmonic solution and two Mascon solutions shows that they present consistent seasonal patterns and trends. The comparison of GRACE-based Wground with groundwater well measurements shows that they present similar long-term trends but with phase difference (~1 month) in seasonal variations.