Seasonal variations and long-term trends of groundwater over the Canadian landmass

Detailed knowledge of groundwater storage improves the understanding and management of water resources. Observations from Gravity Recovery and Climate Experiment (GRACE) satellites have provided data on global terrestrial-water-storage (TWS) changes since 2002. Combining GRACE-TWS and land-surface model (LSM) estimates of soil water, snow-water equivalent and surface-water storage provides a method to quantify groundwater storage (Wground). This study examines the Wground seasonal variations and trends for Canada’s landmass during the period 2003–2016 using GRACE-TWS and the Canadian LSM EALCO (Ecological Assimilation of Land and Climate Observations) model. The results show the study region has a maximum seasonal variation (ΔWground) of 118 mm (volume equivalent 700 km3), with the maximum/minimum Wground appearing in July/April. Eastern Canada has relatively large ΔWground values, up to 400 mm in Newfoundland. The Prairie region has the smallest value (<50 mm). The western and central regions show the maximum/minimum Wground mostly in spring/fall. In contrast, eastern Canada has the maximum/minimum Wground mostly in fall/spring. South Ontario and the Prairie area show the maximum/minimum Wground in summer/winter. Additionally, the Wground trends over the 14-year study period present large spatial variability, with increasing trends of up to 10 mm/year in eastern Canada and decreasing trends (similar magnitudes) in the west. The increasing trend largely offsets the decreasing trend in the study area, and the overall Wground for the region does not show a significant trend during 2003–2016. Comparison of Wground with groundwater well measurements present similar long-term trends but with a phase difference in seasonal variations.


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 (W ground ), soil-water content (W soil ), snow-water equivalent (W snow ), and surface-water storage (W surf ). Thus, GRACE provides a practical way to estimate W ground when W soil , W snow and W surf can be estimated.
The GRACE TWS data have been used to quantify variations and long-term trends of W ground at regional and global scales in various studies (e.g., Rodell et al. 2007Rodell et al. , 2009Swenson 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 W soil , W snow and W surf estimated from other independent approaches such as landsurface models (LSMs), were usually used to separate W ground from GRACE TWS. Famiglietti et al. (2011) and Scanlon et al. (2012) analyzed long-term W ground change in the California Central Valley (USA) region by combining TWS from GRACE, W soil from GLDAS (Global Land Data Assimilation System) LSMs, W snow from the Snow Data Assimilation System (SNODAS) datasets, and W surf from insitu surface-water reservoir measurements. Thomas and Famiglietti (2019) extended the analysis of groundwater change in USA by using W soil and W surf (surrogated by surface runoff) from GLDAS LSMs, W snow 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 globalscale 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. 2014Wang et al. , 2015a. 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 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 W ground for Canada's landmass for the period 2003-2016 using the TWS from GRACE monthly spherical harmonic (SH) solutions and the W soil , W snow and W surf values from the Canadian LSM of EALCO. This study calculates a suite of parameters representing the spatial and temporal variations of W ground , including the occurrence time and the interannual variability of the maximum and minimum W ground in a year, the range of seasonal variation of W ground and its interannual variability, and the longterm trend of W ground 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 W ground 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 W ground , W snow , and W surf . In this study, W ground is computed as 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 , 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 To estimate GRACE-derived W ground in Eq. (1), the simulated soil moisture W soil , snow water equivalent W snow , and surface runoff (as a proxy for surface-water storage W surf ) 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 W soil 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 W soil , W snow and W surf for the period 2003-2016 are calculated from the halfhourly EALCO W soil , W snow and W surf 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 W soil , W snow and W surf anomalies were calculated by subtracting the timemean baseline over the same period to make EALCO and GRACE data comparable.
Eight parameters are calculated from the derived monthly W ground data to characterize the spatiotemporal variations of W ground (Table 1). They include: the long-term trend (τ) represented by the Thiel-Sen linear regression slope of monthly W ground time series for 2003-2016; the months with maximum/minimum W ground in a year (M max /M min ) and their interannual variations (σ max /σ min ); and the range of seasonal W ground variation (ΔW ground ) and its interannual variations (represented by the standard deviation σ Δground and the coefficient of variance COV of ΔW ground ). For the calculation of σ max and σ min , if the difference of M max (or M min ) among different years is larger than 6, the difference is subtracted by 12.
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 km 2 ). 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).
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 W ground comparisons. The well data for groundwater level (W level ) 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 W level change of the wells in each comparison grid is calculated to represent the W level change for that specific grid.

Results
Seasonal variations of W ground Figure 2 shows the average monthly W ground over 2003-2016 and provides a general overview of the spatial and seasonal variations of W ground over the study area. In January, W ground 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 W ground 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 W ground 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 W ground (~75 mm above the annual average) in June and the Western and Northern Hudson Bay MDA has the highest W ground (~50 mm above the annual average) in June. In contrast, the W ground 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 W ground in March or Aprilfor example, the Maritime Provinces MDA and the south portion of the St. Lawrence MDA (mainly in south Ontario), reaches the lowest W ground of~100 mm below the annual average in April. After that, the W ground rapidly increases 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 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 W ground -for example, the north portion of the St. Lawrence MDA and the Northern Quebec and Labrador MDA reach the lowest W ground (~75 mm below the annual average) mostly 1 month later in May and the north Labrador portion has the lowest W ground (over 100 mm below the annual average) 2 months later in June. In June, a large portion of the study area shows a high W ground except that the Northern Quebec and Labrador MDA still presents a low W ground with the value of around 50 mm below the annual average. 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 W ground in summer, after which the W ground reaches its lowest value before winter (Fig. 2)-for example, the Arctic MDA arrives at the lowest W ground (~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 W ground (~60 mm below the annual average). In contrast, the Northern Quebec and Labrador MDA in east Canada mainly shows a rapid increase of W ground in summer. These regions have high precipitation in fall, resulting in another high groundwater recharge before winter. As a result, the W ground 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 W ground rapidly decreases in January. Like in the winter and spring, south Canada also presents different seasonal patterns in W ground variations during summer and fall, whereby it has substantial precipitation after spring, resulting in rapid groundwater recharge. The highest W ground (e.g., over 100 mm above the annual average in the southern Ontario) Fig. 2 The average monthly W ground over 2003-2016, showing its spatial and seasonal variations over the landmass occurs in July or August, after which the W ground rapidly drops in October and then gradually decrease throughout fall.
It is observed from Fig. 2 that the study area overall presents low W ground in early spring and high W ground in early summer. The months having the lowest/highest W ground vary by regions. Figure 3a,b shows the spatial distribution of the months having maximum/minimum W ground (M max / M min ), for 1 year for the study period. As discussed earlier, the M max / M min presents large spatial differences over the study area.
Generally, the study region shows M max mainly in late spring and early summer, except for the Northern Quebec and Labrador MDA which has M max mainly in October or November. In contrast to M max , the M min 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 M min in March and the north portion of the Southwestern Hudson Bay MDA, which has M min in fall (October or November). Figure 3c,d shows the interannual variations (σ max and σ min ) of M max and M min , 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 W ground variation (ΔW ground ), averaged over the study period. The ΔW ground presents large spatial differences with continuous increase from west to east in the study region. The large portion in east Canada has ΔW ground 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 ΔW ground around 120 mm. For the western region, the ΔW ground is low with the smallest ΔW ground of <30 mm in the Prairie region.
The interannual variation of ΔW ground among the 14 years represented by the standard deviation (σ ΔWground ), as shown in Fig. 4b, presents a similar spatial pattern to ΔW ground (Fig. 4a).
However, in relative terms, the interannual variations of ΔW ground 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 W ground 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 Lawrence River water level plunged to its lowest point in more than 30 years, which led to a very low W ground 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 W ground after the severe drought. Significant decreasing trends are observed mainly in the western region. The highest decreasing trend of W ground appears in some areas of the Great Slave Lake MDA. The large decreasing trends of W ground in the western regions are likely due to the decrease of precipitation ).

W ground for major drainage areas (MDAs)
The variations of W ground for the MDAs are summarized in Table 2  The ΔW ground varied between 54.8 and 280.5 mm among the MDAs, with the highest and lowest ΔW ground 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 ΔW ground larger than 100 mm. All the four MDAs in the western region have small ΔW ground values around 50-60 mm, which are less than quarter of that for the Maritime Provinces MDA. For the entire landmass, the W ground reaches the high peak in July and the low peak in April (Table 2), with a ΔW ground of 117.7 mm or 699.9 km 3 . It is worth noting that there is large spatial variability of W ground in each MDA especially in the Maritime Provinces MDA (Figs. 2 and 6).
The long-term trends of W ground for the MDAs (

Comparisons of different GRACE solutions
The W ground 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 W ground 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 W ground 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 W ground 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.

Comparisons to groundwater level (W level ) measurements
Due to the difficulty in acquiring groundwater storage observations, this study compares the SH-derived W ground to the groundwater level (W level ) 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 W level 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 W level variations. Figure 8 compares the time series of GRACE-based W ground and W level (well average) over the study period for the two GRACE SH grids. The W ground and W level present similar trends in both grids. For the grid in the Nelson River MDA (Fig. 8a) (Fig. 8b), both W ground and W level do not present notable long-term trends. The correlation coefficients between W ground and W level are 0.3 for the Nelson River MDA, and 0.26 for the St. Lawrence MDA. It is observed that the W ground and W level seasonal variations show phase differences, which partially result in the low correlation coefficients between W ground and W level . The phase differences can be clearly seen from the monthly average values shown in Fig. 9.
Generally, the W ground shows the seasonal pattern 1 month earlier than the W level for the grid in Nelson River MDA. When moving the W level data 1 month forward, the R increases to 0.37 from 0.3. For the grid in the St. Lawrence MDA, the W ground and W level show different phase differences. After moving the W ground time series 1 month forward, the R between W ground and W level increases to 0.51 from 0.26. be seen from the low correlation coefficients among the well measurements as discussed in the preceding. It is noted that the variability of W level among the well measurements is twice that of the seasonal variation range (ΔW level ) for the grid in the Nelson River MDA, and about half of the ΔW level for the grid in the St. Lawrence MDA, suggesting relatively large uncertainty in W level measurements and small variations of the W level 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 W ground and W level discussed in the preceding.

Discussion
In this study, W ground is estimated using the GRACE-observed TWS and LSM-based W soil , W snow , and W surf from the EALCO model. Errors or uncertainties in TWS, W soil , W snow , and W surf estimates would directly propagate to W ground . Mortimer et al. (2020) evaluated nine Northern Hemishphere gridded W s n o w products from three categories-reanalysis models, passive microwave remote sensing combined with surface observations, and standalone 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 W snow ) 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. 2014Wang et al. , 2015a; 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 W ground 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 longterm trend rather than the seasonal patterns.
Integrating all the errors from TWS, W soil , W snow and W surf , Frappart and Ramillien (2018) showed that the errors on trend estimates of W ground 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. W snow is a primary contribution to TWS in winter time. The large uncertainty of LSM-based estimates of W snow might significantly impact accuracy of the GRACE-derived W ground . Directly quantifying the uncertainty of W ground is still difficult due to the lack of independent data and it needs to be addressed in future studies.  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 W ground for the entire landmass does not show a significant trend over 2003-2016. The comparison of W ground from the GRACE spherical harmonic solution and two Mascon solutions shows that they present consistent seasonal patterns and trends. The comparison of GRACE-based W ground with groundwater well measurements shows that they present similar long-term trends but with phase difference (~1 month) in seasonal variations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.