A Wind-driven Nonseasonal Barotropic Fluctuation of the Canadian Inland Seas

A wind-driven, spatially coherent mode of nonsea-sonal, depth-independent variability in the Canadian inland seas (i.e., the collective of Hudson Bay, James Bay, and Foxe Basin) is identified based on Gravity Recovery and Climate Experiment (GRACE) retrievals, a tide-gauge record, and a barotropic model over 2003–2013. This dominant mode of nonseasonal variability is correlated with the North Atlantic Oscillation and is associated with net flows into and out of the Canadian inland seas; the anomalous inflows and outflows, which are reflected in mean sea level and bottom pressure changes, are driven by wind stress anomalies over Hudson Strait, probably related to wind setup, as well as over the northern North Atlantic Ocean, possibly mediated by various wave mechanisms. The mode is also associated with mass redistribution within the Canadian inland seas, reflecting linear response to local wind stress variations under the combined influences of rotation, gravity, and variable bottom topography. Results exemplify the usefulness of GRACE for studying regional ocean circulation and climate.


Introduction
Hudson Bay, James Bay, and Foxe Basin together constitute the Canadian inland seas (CIS; Fig. 1).This set of marginal seas connects to the Labrador Sea and North Atlantic through Hudson Strait and Ungava Bay in the east, and to the Arctic Ocean through Fury and Hecla Strait and the Gulf of Boothia in the north.These seas are shallow, having depths of ∼ 90-150 m, and broad, spanning an area of ∼ 1 × 10 6 km 2 (Mac-Donald and Kuyzk, 2011).The mean hydrography in Hudson Bay in summer and autumn is such that a thinner, shallower layer of fresher, warmer water sits atop a thicker, deeper layer of saltier, cooler water; during winter and spring, the surface waters cool, and the mixed layer reaches deeper down in the water column (Prinsenberg, 1986a(Prinsenberg, , 1987)).Related climatological features include a seasonal cycle in sea ice, which oscillates between complete ice cover in wintertime and icefree conditions in summertime (Markham, 1986), as well as volume input due to runoff (Déry et al., 2005(Déry et al., , 2011)), transport through Fury and Hecla Strait, and flow from Baffin Bay through Hudson Strait along the Baffin Island coast, which is mostly balanced by volume outflow through Hudson Strait along the Québec coast (Straneo and Saucier, 2008a, b).
The CIS play an important role in the ocean general circulation.Numerical simulations suggest that Hudson Strait is one of the most important regions in the world ocean for the dissipation of tidal energy (Egbert and Ray, 2001;Webb, 2014).Barotropic models show how, on synoptic timescales, dynamic response of Hudson Bay to barometric pressure drives flows through Hudson Strait, which generate coastal waves that subsequently affect sea level downstream as they propagate over the continental shelf (Wright et al., 1987;Greatbatch et al., 1996).Direct measurements of the baroclinic boundary current by a moored current array deployed in Hudson Strait reveal that the outflow through Hudson Strait is responsible for a substantial portion of the fresh water supplied to the Labrador Current and the North Atlantic Ocean (Straneo and Saucier, 2008a).
This region is also of interest in the context of changes ongoing in the Arctic system (White et al., 2007).Passive microwave data reveal that concentrations and extents of seasonal sea ice have decreased in Hudson Bay over recent decades, while historical climate data show that surface air temperatures around Hudson Bay have warmed since 1950 (Hochheim and Barber, 2010;Hochheim et al., 2011).Hy-Published by Copernicus Publications on behalf of the European Geosciences Union.
Figure 1.Shading is the logarithm of ocean depth from ETOPO5 5 min gridded elevation data (National Oceanic and Atmospheric Administration 1988) in the Hudson Bay study area.Units are log 10 (m).Color shading saturates at a value equivalent to 1000 m.Gray arrows schematically represent the sense of the mean regional surface circulation after Prinsenberg (1986b, c) and Drinkwater (1986).The black circle near 59 • N, 94 • W is the location of the Churchill tide gauge.Letters are acronyms for major regional features; alphabetically, they are Foxe Basin (FB), Fury and Hecla Strait (FHS), Greenland (GRN), Gulf of Boothia (GB), Hudson Bay (HB), Hudson Strait (HS), James Bay (JB), Labrador Current (LC), Labrador Sea (LS), Prince Charles Island (PCI), and Ungava Bay (UB).Red outlining in the inset displays the study region location relative to the world ocean.
drometric data indicate strong interannual changes in Hudson Bay streamflow along with a marked shift in the seasonality of river discharge (Déry and Wood, 2004;Déry et al., 2011).However, it remains unclear whether the subsurface waters of the CIS have also undergone change.
Despite their relevance to circulation and climate, the CIS have been grossly undersampled: due to their vast expanse, harsh conditions, and remote location, few campaigns have been dedicated to continuously measuring their subsurface waters; even estimates of the bathymetry in this region can show large uncertainties, especially in the more northern reaches of the CIS.Early observational descriptions of circulation patterns and current structures are derived from sparse data (Prinsenberg, 1986b, c;Drinkwater, 1986).More recent data have facilitated more nuanced descriptions, for example, of the spatial structure of the boundary currents and the role of synoptic eddies in transporting fresh water through Hudson Strait (Straneo and Saucier, 2008a;Sutherland et al., 2011) and the mean state and seasonal cycle in the circulation and hydrography in Hudson and James bays (St-Laurent et al., 2012).However, continuous measurements of subsurface properties remain sparse, leaving open basic questions regarding regional ocean behavior on nonseasonal periods longer than a few days.
Concerns over impacts of climate change (Laidler and Gough, 2003) motivate best use of extant data to provide an understanding of anomalous behavior in this region.While a tide gauge situated at Churchill in southwestern Hudson Bay has measured sea level fluctuations since 1940, and the Gravity Recovery and Climate Experiment (GRACE) spacecraft have observed changes in the mass of ocean and ice over the CIS since 2002, only a few studies have made use of these data to understand the nature of variability in the CIS.Based on the Churchill tide gauge and hydrometric data, Gough and Robinson (2000) posit that sea level variations observed at Churchill partly result from local discharge from the Churchill River.Considering GRACE and an atmospheric reanalysis, Piecuch and Ponte (2014) submit that wind setup could effect mass changes in Hudson Bay.However, these hypotheses are based on statistical metrics (correlation coefficients and coefficients of determination), and it remains to test them using a more dynamically rigorous approach.
In this paper, we investigate nonseasonal oceanic behavior in the CIS.We provide an exploration and interpretation of the data from GRACE and the tide-gauge measurements mostly based on a coarse-resolution barotropic model driven by surface wind stress.The remainder of this paper is organized as follows: in Sect.2, we describe and contrast the observational data; in Sect.3, we describe the ocean model, comparing it to the available observations as well as output from a higher-resolution ocean/sea-ice model, and then use it to understand the leading mode of nonseasonal behavior of the CIS; in Sect.4, we summarize and discuss our results.

Satellite gravimetry
Since their launch in March 2002, the twin GRACE satellites have been monitoring the exchange of water mass between the land and the sea (e.g., Boening et al., 2012).We use monthly ocean bottom pressure estimates derived from Release-05 GRACE time-variable gravity coefficients over the period 2003-2013 to study mass variability in the CIS.The data are taken from the GRACE Tellus server (data version RL05.DSTvDPC1409) and are processed at the University of Texas Center for Space Research (Bettadpur, 2012).Postprocessing by Don P. Chambers (University of South Florida) follows methods described by Chambers and Bonin (2012).Relevant for our purposes, the data are smoothed with a 500 km Gaussian filter, which reduces errors with short wavelengths in the satellite recoveries, but which can also attenuate the magnitudes of the oceanic signals.Relative to Chambers and Bonin (2012), updated estimates are used for degree 2 order 0 coefficients (Cheng et al., 2011) and glacial isostatic adjustment (A et al., 2013).Global spatial-mean values are subtracted from the ocean mass estimates at each time step.The values are provided on a regular 1 • × 1 • horizontal grid; at this resolution, Hudson Bay and James Bay are together represented by 75 grid cells, while Foxe Basin is represented by 11 grid cells.Throughout the paper, we quote values in equivalent seawater thickness units.
Given our interest in nonseasonal behavior, we remove a seasonal cycle from all time series, which we compute by averaging together all January entries, February entries, etc., over 2003-2013 into a 12 month time series.To circumvent difficulties of interpreting the gravity data over the ocean in the presence of large rates of glacial isostatic adjustment over Canada and cryospheric mass loss from Greenland (e.g., Tamisiea et al., 2007;Velicogna, 2009;Rignot et al., 2011), linear trends are removed from all time series using least squares.

Tide-gauge data
A tide gauge maintained by the Canadian Hydrographic Service has measured relative sea level at the mouth of the Churchill River in Churchill, Manitoba (Fig. 1), for more than 70 years.Revised local reference monthly data were extracted from the Permanent Service for Mean Sea Level database (Holgate et al., 2013) on 18 August 2014.Data cover 90 % of months between January 1940 and December 2013, with a complete record existing since July 1991.Given the GRACE record, we consider tide-gauge data over 2003-2013.A set of adjustments is applied to the tide-gauge data.As with the gravimetric estimates, a seasonal cycle and a linear trend are removed.To focus on ocean dynamical signals, we also subtract from the tide-gauge record a global mean sea level time series based on altimetry data (Ablain et al., 2009) as well as the inverted barometer response based on monthly Interim European Centre for Medium-Range Weather Forecasts Reanalysis (ERA-Interim) (Dee et al., 2011) mean sea level pressure fields and Eq.(1) of Ponte (2006).We note that removal of the global mean and inverse barometer signals reduces the detrended monthly variance in the tide-gauge sea level time series over 2003-2013 by 40 %.Given our focus on detrended behavior, we do not make any further corrections for vertical land motion, instead assuming that the relevant geophysical processes (e.g., postglacial rebound) are represented by a linear trend over the analysis period (cf.Santamaría-Gómez et al., 2012).

Data comparisons
One concern of using GRACE ocean bottom pressure estimates over the CIS is that they might be contaminated by transient terrestrial water storage from surrounding watersheds.Root mean square values of monthly water storage estimated by a land hydrology model can be 5-10 cm equiv-alent water thickness in parts of the Hudson Bay drainage basin (Landerer and Swenson, 2012).Given the averaging function applied to the gravity data, such land signals could leak into the ocean data (Wahr et al., 1998).
To determine whether the ocean bottom pressure estimates are polluted by leakage of terrestrial water storage, we consider nonseasonal time series of GRACE bottom pressure averaged over the CIS alongside GRACE water storage 1 averaged over the Hudson Bay drainage basin (Fig. 2a). 2 An analogous technique is used by Peralta-Ferriz et al. (2014) to determine whether GRACE data over the Kara and Barents seas are contaminated by land leakage over the respective watersheds that drain into them.The two time series appear visually distinct and their correlation coefficient (0.26) is not statistically significant (Fig. 2a).The two signals do not show meaningful coherence at any frequency (not shown).These results demonstrate that the GRACE data over the CIS are not overwhelmed by leakage of terrestrial water storage at nonseasonal periods.Additional analysis corroborates this conclusion; computing correlation coefficients between the averaged terrestrial water storage curve and bottom pressure time series at individual GRACE ocean grid cells, we find that there are no points within the CIS where the local ocean bottom pressure is significantly correlated with the large-scale land signal on nonseasonal timescales (not shown).
As an additional check on the GRACE ocean data quality, and also to give physical insight, we compare bottom pressure estimates averaged over the CIS to sea level observed at the Churchill tide gauge (Fig. 2b).Notwithstanding the difference in amplitude, which probably partly reflects attenuation of the true ocean bottom pressure signal by the spatial averaging involved in the postprocessing (see Sect. 2.1.1),there is close correspondence between the two curves (Fig. 2b).The overall correlation coefficient between the nonseasonal sea level and bottom pressure (0.58) is statistically significant at the 95 % confidence level.The correspondence between the two time series in Fig. 2b is consistent with our interpretation of the results in Fig. 2a, attesting to the general meaningfulness of the nonseasonal signals in the GRACE gravity data over the ocean.What is more, this result suggests that the nonseasonal behavior at Churchill partly reflects bay-wide, depth-independent (that is, barotropic) variability.
1 GRACE terrestrial water storage estimates were processed at the University of Texas Center for Space Research and postprocessed by Sean Swenson (National Center for Atmospheric Research).The gridded estimates, provided on a 1 • × 1 • grid, were downloaded from the GRACE Tellus server (data version RS05.DSTvSCS1409) and scaled following Landerer and Swenson (2012). 2 The Hudson Bay drainage basin has been defined as the union of the Hudson Bay seaboard and Nelson River basins, determined based on watersheds data provided by the Commission for Environmental Cooperation (http://www.cec.org/Page.asp?PageID=924&ContentID=2866).

Model framework
The main drivers of mass redistribution in a homogeneous ocean are surface loading and wind stress (Hughes, 2008).In the case of synoptic timescales (i.e., periods of a few days), surface loading by barometric pressure can drive a relatively large non-equilibrium sea level response on the continental shelf (standard deviations 10 cm) (Greatbatch et al., 1996), related to the fact that propagation speeds of barotropic gravity waves are reduced in shallow regions.Indeed, Ponte and Vinogradov (2007) suggest that the assumption of an inverted-barometer response is not appropriate at periods of about 1 month and shorter for Hudson Bay (their Fig. 6).However, in the case of sub-synoptic timescales (e.g., periods longer than 1 month), non-isostatic adjustment of sea level on the shelf to pressure loading is thought to be comparatively small (standard deviations 1 cm) (Greatbatch et al., 1996).Thus, we expect that observed sea level variations in Fig. 2 (which are corrected for an inverted barometer response) are driven by wind stress rather than by surface loading.
To assess this expectation, we consider numerical solutions from the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997).The setup solves the primitive equations on a coarse-resolution (0.5 • × 0.5 • ) spherical polar grid having quasi-global spatial coverage, with solid walls imposed at 79 • north and south latitude.(At this resolution, the model represents the meridional breadth of Hudson Strait, which varies from ∼ 100 to 200 km, depending on longitude (Fig. 1), using between 2 and 4 grid cells.)Boundary conditions are in the form of surface fluxes based on monthly means of instantaneous zonal and meridional turbulent surface wind stresses taken from ERA-Interim (Dee et al., 2011).The reanalysis fields, which are provided on a regular 1.5 • × 1.5 • grid, are bilinearly interpolated onto the model grid.Bottom topography is based on 5 min gridded elevations and bathymetry for the world (ETOPO5) data (National Oceanic and Atmospheric Administration, 1988); due to sparsity of measurements, it is likely that this bathymetric data set has large uncertainties in the CIS.The bathymetry is averaged within 0.5 • × 0.5 • bins and then smoothed with a two-dimensional 2 • × 2 • boxcar function.We use a constant value for density (1029 kg m −3 ) and a single layer in the vertical.Variable ocean depths are implemented using partial cells (Adcroft et al., 1997).The model uses a linear free surface boundary condition along with a 900 s time step for the momentum equations.We use a vertical eddy viscosity of 5 × 10 −4 m 2 s −1 and a grid-dependent lateral eddy viscosity varying from about 1 × 10 4 m 2 s −1 at low latitudes to roughly 3×10 3 m 2 s −1 at high latitudes.Simulations are run forward in time from rest for 35 years beginning on 1 January 1979 (the temporal range of ERA-Interim).To be consistent with the observations (Fig. 2), we consider monthly averaged model output over 2003-2013 with the seasonal cycle and a linear trend removed.
This framework is admittedly simple; many effects (e.g., mesoscale eddies, sea ice, river runoff, wind stress over the Arctic Ocean) have been precluded.On account of the coarse grid resolution, we do not resolve the topographic gyres and current separations induced by bathymetry in Hudson Bay discussed by Wang et al. (1994) in the context of a finerresolution model.Due to the lack of ocean stratification, we do not capture the baroclinic boundary current in Hudson and James bays discussed by St-Laurent et al. (2012) among others.Given the lack of an interactive sea-ice model, we do not simulate any role played by sea ice in mediating the transfer of momentum between the atmosphere and the ocean (St-Laurent et al., 2011).However, to the extent that our model agrees with the data of interest, we can conclude that any omitted physics is unimportant in the present context.

Comparing model and data
Before comparing it to the GRACE ocean data, we smooth the model bottom pressure using the same 500 km spatial filter used in the GRACE postprocessing (Chambers and Bonin, 2012) and then average over the CIS, interpolating onto the GRACE ocean grid.The statistically significant correlation coefficient between the model and data bottom pressure is 0.69, with the model curve explaining 47 % of the data curve's variance (Fig. 3a).The standard deviation from the model (1.0 cm) is smaller than the standard deviation from the data (1.5 cm), which could be partly due to residual noise in the data.The gross correspondence between the time series speaks to the realism of the two independent estimates, consistent with our prior assessment of GRACE data quality based on observed water storage and sea level (Fig. 2).It demonstrates that the ocean model suffices for capturing the major features of observed nonseasonal bay-wide fluctuations in ocean mass.
For comparing against the tide-gauge observations, we consider model sea level from the grid cell whose centroid is the closest to the Churchill site (Fig. 3b).The model time series roughly reproduces the gross features of the observed sea level curve.However, the model underestimates the data amplitudes: standard deviations of the observed and simulated signals are 5.7 and 3.0 cm, respectively, possibly reflecting the importance of forcing terms that have been omitted from our model, for example, riverine discharge (Gough and Robinson, 2000), or perhaps indicating that the model underestimates the wind-driven ocean response, for instance, on account of the coarse resolution of the forcing data set.Overall, the correlation coefficient between model and data curves is 0.82, with the model explaining a majority (58 %) of the observed variance.
To gauge the influence of horizontal resolution and missing physics (ocean stratification, sea-ice dynamics, etc.) on the correspondence between the data and our model, we also consider a higher-resolution ocean/sea-ice model.Monthly sea level and bottom pressure from the Estimating the Cir-culation and Climate of the Ocean Phase-II (ECCO2; Menemenlis et al., 2005) cube92 solution were obtained for 2003-2012.This estimate of the ocean/sea-ice state, generated by the MITgcm coupled with a fully interactive sea-ice model, is defined on a global "cubed sphere" topology, with a nominal horizontal resolution of 0.25 • ×0.25 • and 50 vertical levels.Surface forcing is essentially based on unadjusted fields from the Japanese 25 year Re-Analysis (JRA-25; Onogi et al. 2007), except for precipitation (JRA-25 adjusted to remove large resultant drifts in salinity and global sea level in the model solution) and runoff (which, for the CIS and the Arctic Ocean, is derived from monthly mean river discharge from the Arctic Runoff Database).Some internal model parameters were previously adjusted to better fit observations.Due to inclusion of ocean stratification, sea level and bottom pressure are not generally equivalent in the ECCO2 solution.Also, in the presence of sea ice, "sea level" is defined as the physical depression of the sea surface plus the sea-ice load in equivalent water thickness units.
Comparable ECCO2 curves for nonseasonal ocean bottom pressure averaged over the CIS and sea level at the Churchill tide-gauge location are overlaid in Fig. 3a and b, respectively.Perhaps surprisingly, for the common period 2003-2012, our simple barotropic model simulation performs as well as (if not better than) the ECCO2 cube92 solution in reproducing the data.While the ECCO2 solution explains 36 and 54 % of the variances in the GRACE and tide-gauge time series, respectively, the barotropic simulation explains 50 and 58 % of the respective observed variances (Fig. 3).
Results in Fig. 3 show that our simple dynamical framework is sufficient to capture the lowest-order monthly nonseasonal behavior observed in the CIS across a range of spatial scales, and imply that the influences of ocean stratification, sea ice, and surface fluxes of mass and buoyancy are higher order.The realism of the model encourages its further exploration to more fully understand the nonseasonal variability in the CIS.In what follows, we focus on sea level, but note that, since sea level and bottom pressure are equivalent quantities in a barotropic ocean, the results also apply to bottom pressure.

Empirical orthogonal functions
The relationship between observed sea level and bottom pressure changes (Fig. 2b) led us to hypothesize the existence of a bay-wide depth-independent oscillation.To determine more rigorously whether there is in fact such a wind-driven nonseasonal barotropic fluctuation of the CIS, we perform an empirical orthogonal function (EOF) analysis, which boils down to solving for the eigenvalues and eigenvectors of the covariance matrix of a scalar that varies in space and time (von Storch and Zwiers, 1999).
The leading eigenvector of simulated sea level over the CIS shows a single-signed spatial structure (Fig. 4a), but values are larger over the deep interior region and smaller in the shallow boundary area.The mode could be caused by a combination of local and remote effects, perhaps with remote mechanisms forcing water into and out of the domain and local drivers acting to redistribute mass within the domain.This leading empirical mode explains 69.4 % of the total nonseasonal simulated sea level variance in the CIS (Fig. 4a).Local variance explained is highest (> 90 %) along a "ring" around Hudson Bay separating interior and boundary regions; this ring of strong correlation could reflect rapid Kelvin wave propagation around Hudson Bay, with a possible analogy to the coherent sea level fluctuations along the global continental slope observed by Hughes and Meredith (2006).Explained variance is lowest (< 10 %) in shallow regions of Foxe Basin southeast of Prince Charles Island.Over most of Hudson Bay's shallow boundary (e.g., near Churchill) and deep interior, the mode accounts for about two-thirds of the local sea level variance.
The leading expansion coefficients show variability across all accessible timescales, and a dominant period is not visually obvious (Fig. 4b), while an estimate of the associated power spectral density is slightly red in nature (not shown).
We observe that the expansion-coefficient time series is essentially equivalent to the bay-mean sea level signal: the correlation coefficient between the two curves is ∼ 1 (not shown).There appears to be a relationship between the North Atlantic Oscillation (NAO) and the leading expansion coefficients such that anomalous sea level in the CIS is high when the NAO index is low, and vice versa; the statistically significant correlation coefficient between the two time series (−0.59) confirms this out-of-phase relationship.
Looking further afield, we compute correlation coefficients between the expansion coefficients (Fig. 4b) and nonseasonal sea level time series at each model grid cell over the global ocean (not shown).A statistically significant in-phase relationship is apparent between fluctuations in the CIS, Baffin Bay, and the Mediterranean Sea, while a significant outof-phase relationship is evident between the expansion coefficients and variations over the midlatitude North Atlantic and along parts of the North Sea.These relationships are similar to those suggested by Piecuch and Ponte (2014), who computed correlations between the leading mode of nonseasonal bottom pressure variability over the midlatitude North Atlantic and anomalous bottom pressure elsewhere based on GRACE data (e.g., their Fig.3a).
Assuming geostrophy, fluctuations in sea level are proportional to variations in the barotropic stream function, and therefore this mode can be physically interpreted in terms of ocean circulation changes.For example, the domed shape of the eigenvector in Hudson Bay (Fig. 4a) suggests anomalous anticyclonic (cyclonic) circulation when the expansion coefficients are positive (negative).Given the cyclonic sense of Hudson Bay's mean circulation (e.g., St-Laurent et al., 2012), this mode thus corresponds to spin-up and -down of the barotropic component of the mean circulation in Hudson Bay roughly during positive and negative NAO periods, respectively.
Model results in Fig. 4 corroborate our earlier suspicion based on data that there exists a wind-driven barotropic fluctuation of the CIS that explains most of the nonseasonal sea level variance across a range of spatial scales.What is more, this mode of oscillation is correlated with the NAO, implying that some of the anomalous sea level behavior in the CIS is tied to climate variability more broadly over the North Atlantic sector, consistent with suggestions made by Gough and Robinson (2000).It remains to be determined, however, what the important regions of wind forcing are, and what the relevant ocean dynamics are.We turn to these questions in the next section.

Forcing and dynamics
Seeing as the leading expansion coefficients are correlated with the NAO (Fig. 4b), we now consider the structure of nonseasonal wind stress forcing over the northern North Atlantic sector (Fig. 5).Standard deviations of nonseasonal wind stress are on the order of a few hundredths of a N m −2 .Noteworthy are strong variations in zonal wind stress near Cape Farewell, Greenland (Moore and Renfrew, 2005).
To suggest relationships between wind stress over the northern North Atlantic and sea level in the CIS, we compute correlations between the expansion coefficients (Fig. 4b) and nonseasonal zonal and meridional wind stress (Fig. 5).Zonal winds over a broad swath extending from Hudson Bay and the Labrador Sea across the northern North Atlantic Ocean to the North, Norwegian and Barents seas are significantly negatively correlated with sea levels over the CIS (Fig. 6a).Sea levels over the CIS are significantly negatively correlated with meridional winds over eastern Hudson Bay and the western Labrador Sea as well as from the northeastern North Atlantic to the Norwegian Sea, while positive correlations are seen off the southwestern coast of Greenland (Fig. 6b).Given the correlation between the expansion coefficients and the NAO (Fig. 4b), these correlation patterns are consistent with wind stress anomalies associated with the NAO; for example, anomalous westerly winds occur during positive phases of the NAO (e.g., Marshall et al., 2001), when sea levels over the CIS are anomalously negative (Fig. 4b).
The strongest correlations between the expansion coefficients and zonal winds occur over Hudson Strait (Fig. 6a), suggesting that the sea level in the CIS might be influenced by wind stress variations over adjacent regions.For example, wind stress along Hudson Strait directed towards the CIS would force flow into the CIS until dynamic balance is established between the along-strait sea level gradient and wind stress, i.e., a wind setup (e.g., Csanady, 1981); in the absence of additional boundary forcing, the CIS would undergo barotropic adjustment in response to the mass inflow, likely resulting in a horizontally uniform sea level increase.
(An analogous scenario can be entertained for winds over Hudson Strait directed away from the CIS.) To establish what the influence of wind stress over Hudson Strait is, we perform another numerical simulation based on the model framework described previously (Sect.3.1) by setting the surface wind stress to zero everywhere except over Hudson Strait (see Fig. 5), all else (e.g., the time period of integration) being equal.This Hudson Strait winds experiment captures some of the variability in the CIS from the baseline experiment examined in Fig. 4. Namely, winds over Hudson Strait effect bay-wide changes in sea level over the CIS, and the correlation between the leading sea level expansion coefficients from the baseline and Hudson Strait winds experiments is 0.58 (Fig. 7), demonstrating that Hudson Strait winds are important to variability in the region.
But there are also disagreements between behaviors generated by the two experiments.Wind stress over Hudson Strait effects sea level changes over the CIS whose amplitudes are horizontally uniform, implying that this experiment lacks important local effects (e.g., winds over the CIS) responsible for generating the spatially varying amplitudes that are manifested in the baseline experiment.What is more, the baymean sea level changes from the Hudson Strait winds experiment are smaller than those from the baseline experiment (cf.amplitudes in Figs.4a and 7a), indicating that wind driv- ing in remote regions (e.g., over the North Atlantic Ocean) is also important in forcing water into and out of the CIS.
To more completely account for the behavior from the baseline experiment, another simulation is performed over the same time period by driving the model with wind stress only over regions where the correlation coefficient between wind stress and baseline expansion coefficients is statistically significant (see the colored regions in Fig. 6).This correlated winds experiment generates variability in the CIS that is extremely close to the behavior produced by the baseline simulation (Fig. 8), explaining 95 % of the variance in the baymean sea level signal from the latter (Fig. 8b); spatial patterns of the leading modes of variability from the two model experiments are practically identical (Figs.4a and 8a).
The remote surface forcing is potentially communicated to the CIS by a variety of physical mechanisms that can influence sea level in and around Hudson Strait.We speculate that these could include, for example, coastally trapped propagating waves forced over shallow areas (e.g., around Cape Farewell, Greenland), and planetary and topographic Rossby waves forced by wind curl patterns over the deep ocean.Finally, to infer what the relevant local dynamics are producing the spatial structure of the variability inside the bay (Fig. 4a), we run a suite of experiments treating the western entrance of Hudson Strait (surrounding Mill, Nottingham, and Salisbury islands) as an open boundary, with flow into and out of the CIS (i.e., bay-mean sea level changes) from the baseline experiment specified a priori.Simulations are identical in all respects other than that we choose to either variously "turn off" one of the terms in the model's momentum equations (e.g., advection, Coriolis, surface forcing, pressure gradient) or change the topography of the CIS.
Turning off winds over the CIS results in a leading mode of sea level variability whose spatial structure is horizontally uniform, similar to the Hudson Strait winds experiment case (Fig. 7a).In contrast, removal of quadratic bottom drag or nonlinear terms from the model dynamics has no perceptible effect, and the variable spatial structure from the baseline experiment is recovered (cf.Fig. 4a).Finally, setting the depth to a constant value (of about 250 m), or turning off either the Coriolis acceleration term or the pressure-gradient contribution, leads to spatial structures for the leading sea level variability very different from those plotted in Fig. 4a.Thus, based on these experiments, we reason that the spatial structure of the leading mode of sea level variability in the CIS reflects linear response to local winds governed by rotation, gravity, and topography.

Conclusions
Using satellite gravimetry, a tide gauge, and a barotropic model, we identified a wind-driven nonseasonal barotropic fluctuation of the Canadian inland seas (CIS) that is correlated with the NAO (Figs. 2-6).Anomalous inflows and outflows, which are reflected in spatially averaged changes in sea level and bottom pressure over the CIS, are driven by wind stress over Hudson Strait (Fig. 7), probably through a wind setup, and the northern North Atlantic Ocean (Fig. 8), possibly communicated by means of wave mechanisms.Anomalous mass redistribution within the CIS, which relates to changes in the depth-mean circulation, is governed by the linear ocean response to more local wind stress variations under the joint influences of rotation, gravity, and variable bottom topography.We observe that, while it suggests broad regions of wind forcing over the northern North Atlantic potentially contributing to CIS variability (Fig. 6), our analysis does not unambiguously pinpoint which forcing regions are most relevant -in fact, it could be that forcing over just one small area of the northern North Atlantic is the main driver.An adjoint model -capable of quantifying the sensitivity of a particular modeled quantity at a specific place and time to all model inputs and states at preceding times -could be used to shed more light on which regions of wind forcing most influence the CIS, as is done by Fukumori et al. (2007) to elucidate a near-uniform basin-wide sea level fluctuation of the Mediterranean Sea, but such an analysis is beyond our scope and deferred to future study.
Our findings complement previous modeling work on Hudson Bay (Wang et al., 1994;Saucier and Dionne, 1998;Saucier et al., 2004;St-Laurent, 2011, 2012).Whereas past studies tend to regard Hudson Strait as an open boundary, specifying inflow and outflow at the outset, our forcing experiments highlight wind stress changes over adjacent and remote areas responsible for driving mean sea level changes in the CIS (Figs. 7 and 8), emphasizing the need for accurate estimates of atmospheric variability for modeling the circulation in the CIS.Our ability to reproduce qualitatively the data (Fig. 3) without recourse to ice-ocean interactions accords with St-Laurent et al. (2011), who find that ice plays only a small role in mediating seasonal momentum transfer between air and sea, reflecting the loose, mobile nature of sea ice in Hudson Bay.Similar to Wang et al. (1994), we find that the effects of variable bottom topography, which are ignored in the flat-bottomed conceptual model of Hudson and James bays due to St-Laurent et al. (2012), are an important determinant of circulation changes in the CIS.(We note that St-Laurent et al. (2012) recover some of the effects of topo-graphic steering by assuming that there is a strong current in the boundary region.) Investigating all periods from monthly data between 1974 and 1994, and based on a correlation analysis, Gough and Robinson (2000) conclude that 43 % of the variance in the Churchill tide-gauge record can be understood in terms of a response to Churchill River discharge.Considering nonseasonal periods from detrended monthly observations over 2003-2013, and using a barotropic model, we submit that at least 58 % of the sea level variance at Churchill is driven by wind stress anomalies (Fig. 3b).The fact that our emphasis (on remote driving and wind stress) differs from that of Gough and Robinson (2000) (on local forcing and river runoff) could reflect the distinct timescales and periods being considered; for example, during the late 1970s and early 1980s, the tide-gauge record is dominated by decadal decline in sea level (Gough and Robinson, 2000, Fig. 3).These considerations, along with the uniqueness of the tide-gauge record, provide ample motivation for more general future works to reconcile the relative roles of wind stress and river runoff, thus painting a more complete portrait of the sea level behavior at Churchill.
One of the challenges of using GRACE data over the ocean in near-coastal regimes is separating oceanic signals from non-oceanic noise (Chambers and Schröter, 2011).In the case of the CIS, this is an especially difficult task, given the large rates of mass loss from the Greenland ice sheet, ongoing postglacial rebound over Canada, and any terrestrial water storage tied to changes in river discharge.Our results (Figs. 2 and 3) suggest that meaningful estimates of nonseasonal ocean bottom pressure behavior can be derived from GRACE retrievals over the CIS.

Figure 2 .
Figure 2. (a) Nonseasonal GRACE Release-05 2003-2013 ocean bottom pressure averaged over the Canadian inland seas (black) and terrestrial water storage averaged over the Hudson Bay drainage basin (red).(b) Nonseasonal GRACE Release-05 2003-2013 ocean bottom pressure averaged over the Canadian inland seas (black) and Churchill tide-gauge sea level (red).Sea level time series is multiplied by 0.25 to fit within axis limits.All quantities shown in equivalent water thickness with units of cm.

Figure 3 .
Figure 3. (a) Nonseasonal time series of ocean bottom pressure averaged over the Canadian inland seas from GRACE Release-05 (black), the barotropic model (red), and the ECCO2 solution (cyan).(b) Nonseasonal time series of sea level at Churchill measured by a tide gauge (black), the barotropic model (red), and the ECCO2 solution (cyan).Units are cm.

Figure 4 .
Figure 4. Leading empirical orthogonal function (EOF) of nonseasonal sea level determined from the barotropic model in the Canadian inland seas.(a) Shading is the value of the leading eigenvector (cm).Contouring is the local fraction of simulated anomalous sea level variance explained by the leading EOF.(b) Nonseasonal time series of the leading expansion coefficients (black) normalized to have unit variance.Also shown is the NAO with a linear trend and a seasonal cycle removed (red).The NAO time series is taken from the National Oceanic and Atmospheric Administration's Earth System Research Laboratory Physical Sciences Division website.

Figure 5 .
Figure 5.Standard deviations of nonseasonal (a) zonal and (b) meridional wind stress from ERA-Interim used to force the barotropic model.Units are N m −2 .Thick black box drawn around Hudson Strait is the control volume for the numerical experiment with altered surface wind stress forcing described in Sect.3.4.Thin black contours are the model's 200, 1000, and 2000 m isobaths, shown for reference.

Figure 6 .
Figure 6.Correlation coefficients between the expansion coefficients of the leading empirical orthogonal function (Fig. 4b) and nonseasonal (a) zonal and (b) meridional ERA-Interim wind stress.Only values statistically significant at the 95 % confidence level are shown.Thin black contours are as in Fig. 5.

Figure 7 .
Figure 7.As in Fig. 4 but shown for the experiment driven only by wind stress over Hudson Strait.In (b), the nonseasonal time series of the leading expansion coefficients from the baseline experiment is shown in gray for reference.

Figure 8 .
Figure 8.As in Fig. 4 but shown for the experiment driven only by wind stress over regions of statistically significant correlation coefficients shown in Fig. 6.In (b), the nonseasonal time series of the leading expansion coefficients from the baseline experiment is shown in gray for reference.