Climate Change, Mean Sea Level and High Tides in the Bay of Fundy

Abstract

Global warming is predicted to result in a rise in sea level which will lead to increased flood risk. Two other factors that will affect high water are the existing trends in mean sea level and changing tides. We illustrate here that in the Bay of Fundy and Gulf of Maine these two are related. An analysis of long-term sea level records shows that, independent of global warming related to climate change, sea level and tidal range have been increasing in this system. Our numerical model investigation indicates that recent changes in sea level, attributed in part to post-glacial rebound, are giving rise to increasing tides. The combined effects of modern sea level rise, global warming induced sea level rise, and the expanded tidal range they induce, will produce a significant increase in the high water level, much greater than that found when considering modern climate-induced sea level changes in isolation. We are predicting a dramatic increase in the risk of flooding at higher high water during the twenty-first century.

RÉSUMÉ [Traduit par la rédaction] Le réchauffement climatique devrait entraîner une hausse du niveau de la mer, d'où l'augmentation des risques d'inondation. Mais deux autres facteurs sont à l'origine de cette hausse, soit les tendances actuelles liées au niveau moyen de la mer et au changement dans l'amplitude des marées. Dans notre étude, nous démontrons que ces deux facteurs sont interreliés dans la baie de Fundy et dans le golfe du Maine. L'analyse des relevés du niveau de la mer à long terme démontre que le niveau de la mer et l'amplitude des marées sont en hausse dans ce système, abstraction faite du réchauffement causé par les changements climatiques. L’étude que nous avons effectuée au moyen d'un modèle numérique fait ressortir que les changements survenus récemment dans le niveau de la mer, attribuables en partie au relèvement post-glaciaire, expliquent l'amplitude croissante des marées. L'effet de l’élévation actuelle du niveau de la mer, conjugué à celui de l’élévation du niveau de la mer provoquée par le réchauffement climatique et de l'amplitude croissante des marées que ces deux phénomènes induisent augmenteront considérablement le niveau des hautes eaux, beaucoup plus que ce qu'indiquent les conclusions des études sur les changements du niveau de la mer induits par le climat actuel pris isolément. Nous prévoyons une augmentation importante des risques d'inondation dans des conditions de pleine mer supérieure au cours du XXI^e siècle.

Keywords:

• sea level rise
• tide
• Fundy
• climate change

1 Introduction

In the Bay of Fundy (Fig. 1) the large tidal range is always a factor in incidents of flooding. Surges not occurring near high water will most likely produce water levels that are within the normal tidal range. The high tides are a result of the near resonance of the semi-diurnal forcing with the natural period of the system. The domain of the resonance encompasses the Bay and the Gulf of Maine to the edge of the continental shelf (Garrett, 1972, 1974; Greenberg, 1979). As the mean water depth rises because of climate change and crustal subsidence, the wave propagation speed will increase and, hence, the resonance frequency will change. Godin's (1992) study of Saint John sea level suggests that the amplitudes of the semi-diurnal tidal constituents are changing and identifies the relative sea level rise due to isostatic processes as a possible cause. Ray (2006) and Müller et al. (2010) noted the same trend in Saint John and other Fundy–Maine tide observations. If this is a Fundy–Maine phenomenon, the response is expected to be similar throughout the Bay of Fundy and Gulf of Maine but different in the adjacent waters as noted in Garrett (1972).

Fig. 1 The Bay of Fundy and Gulf of Maine with locations mentioned in the text indicated. The bathymetry contours are in metres.

Storms combined with high tides have given rise to flooding in the Bay of Fundy. Some, such as the Saxby Gale, 1869 and the Groundhog Day Storm of 1976 (Parkes et al., 1997), have reached legendary status for their destructive histories. This study found that storm surges along Canada's Atlantic Coast with a magnitude greater than 0.6 m occur two to four times per year, but do not always inflict damage. The average duration of such surges is 2.2 hours, but they have been observed to persist for over 12 hours. In the macrotidal Bay of Fundy, serious flooding will only occur when such surges coincide with the one to two hours that the tide is near its maximum. Storm surges arriving at other times in the tidal cycle are much less problematic because the surge-augmented water level is most likely to fall within the normal tidal range. Even at higher high water, the flooding risk is diminished in the Bay of Fundy if the surge coincides with an average or mean tide instead of a large tide which at Saint John can differ by 1.2 m (DFO, 2011). In the outer Bay this difference is less (0.7 m at Yarmouth); at the head of the Bay the difference is significantly greater (2.1 m at Burntcoat Head in Minas Basin).

Investigations of sea level rise in Canada's Maritime Provinces were carried out earlier for Charlottetown (P. Lane and Associates Ltd., 1988) and Saint John (Martec Ltd., 1987). These looked at storm surge statistics and land topography to assess climate change impacts. They took an almost arbitrary one metre rise in sea level and assumed other factors would remain constant. Both studies predicted major residential, industrial and infrastructure impacts. Neither study made any allowances for other sources of sea level change or changing dynamics from shifts in sea level.

The coast of the Bay of Fundy and Gulf of Maine is still adjusting to the ending of the last ice age. The glaciers centred in Hudson Bay (north central North America) depressed the middle of the continent, tilting up the margins out to the edge of the continental shelf (Shaw et al., 1997). The retreating of the ice has reversed the picture with the centre of the continent rising and the coast and continental shelf falling. Tectonic models of this glacial rebound (Peltier, 2002, 2004) show how this process is not uniform in the area. These models compute the interactions of ice sheets, the earth's mantle and the redistribution of the ocean water mass. Computations have indicated sea level changes in the Gulf of Maine and Bay of Fundy that increase out to the edge of the continental shelf.

The determination of changing flood risk over the next century will be more complex than simply adding future contributions determined from climate models. The sea level changes from other processes such as glacial rebound and the changes in tides in the Bay of Fundy and Gulf of Maine also need to be considered. In this study we give an overview of an analysis of long-term sea level data to determine rates of sea level rise and changes in tidal constants; we also use a numerical tidal model to investigate how climate change scenarios might combine with the the already changing local sea level and tides to affect tidal high water. We do not take into account short-term meteorologically driven sea level events. Wind-induced surges, inverse barometer effects, sea surface waves and swell are not included here and will be additive (in some complex manner) to the effects we do consider. The possible changing nature of these events in climate change scenarios definitely needs consideration.

2 Results from the data analysis

There were three objectives for this analysis: first, we wanted to confirm the earlier work of Godin (1992), Ray (2006) and Müller (2011) that predicted increasing tides in the region. Second, consistent with Garrett's (1972) thesis that the Bay of Fundy and Gulf of Maine were part of a single resonant system, we wanted to demonstrate that increased tides would be seen in the Gulf of Maine and Bay of Fundy stations (Boston, Portland, Eastport and Saint John) and that changes within the system would be different from those seen at a station outside the system (Halifax). Third, to the extent possible, we wanted to establish quantitative estimates for the present rates of sea level rise and changes in tidal amplitudes.

The data from the various stations were not of equal quality. There were gaps in the data; some records were short, and some tide gauges had operational problems. We have excluded data from Yarmouth and Bar Harbor for these reasons. The stations in which we have the most confidence in the results are Boston, Portland, Saint John and Halifax (Figs 1 –3). We have made only limited use of the Eastport analysis because the record is somewhat shorter than the others and includes several data gaps, using it in model–data tidal change comparisons but not in the determination of existing mean sea level rise. We understand that because of wharf construction near the gauge, there were siltation problems influencing the Saint John water level data for several years starting around 1980; the gauge was moved in 1999. The intake pipe was modified to prevent damping problems in 2006, and there are possible Bernoulli effects from high currents (Charles O'Reilly and Phillip MacAulay CHS, personal communication, 2004, 2008).

Fig. 2 Time series of the annual mean sea level (Z₀) from analysis. The change in trend pre- and post-1980 was found to be statistically significant, as was a drop in mean sea level. The straight black lines show the trends and the coloured lines indicate the confidence intervals.

Fig. 3 The analyzed annual value of M₂ amplitude showing the increasing trends in the Gulf of Maine and Bay of Fundy and the abrupt shift in value at 1982. The black curve is our fit to the data; the coloured curves are the 95% confidence intervals.

Fig. 4 Sea level rise per century as determined from the data analysis. The pre-break analyses from Boston, Portland, Saint John and Halifax are considered to be the most reliable.

Data for Saint John were obtained from DFO (2010), and all other data records were obtained from Hawaii (2010). Yearly records of hourly sea level measurements were subjected to harmonic regression analysis using the program of Pawlowicz et al. (2002). Only years with greater than 50% non-missing data (between 4380 and 8766 hourly measurements) were considered, and the amplitudes and phases for the fixed set of constituents listed in Table 1 were estimated. Nodal adjustments were carried out for estimating phases but not for amplitudes caused by non-linear effects as discussed by Ku et al. (1985). No constituents were inferred. A summary of selected results for the M₂ and Z₀ tidal constituents are included here and put in the context of the modelling undertaken in this study.

Table1. Tidal constituents analyzed. Frequencies are in cycles per hour.

Constituent	Freqency	Constituent	Frequency
Z0	0	M2	0.0805114
S sa	0.0002282	MKS2	0.0807396
M sm	0.0013098	λ2	0.0818212
M m	0.0015122	L2	0.0820236
M sf	0.0028219	S2	0.0833333
M f	0.0030501	K2	0.0835615
α i	0.0343966	MSN2	0.0848455
2Q1	0.0357064	η2	0.0850736
σ1	0.0359087	MO3	0.1192421
Q1	0.0372185	M3	0.1207671
ρ1	0.0374209	SO3	0.1220640
O1	0.0387307	MK3	0.1222921
τ1	0.0389588	SK3	0.1251141
β1	0.0400404	MN4	0.1595106
NO1	0.0402686	M4	0.1610228
χ1	0.0404710	SN4	0.1623326
P1	0.0415526	MS4	0.1638447
K1	0.0417807	MK4	0.1640729
φ1	0.0420089	S4	0.1666667
θ1	0.0430905	SK4	0.1668948
J1	0.0432929	2MK5	0.2028035
SO1	0.0446027	2SK5	0.2084474
OO1	0.0448308	2MN6	0.2400221
υ1	0.0463430	M6	0.2415342
OQ2	0.0759749	2MS6	0.2443561
ϵ2	0.0761773	2MK6	0.2445843
2N2	0.0774871	2SM6	0.2471781
MU2	0.0776895	MSK6	0.2474062
N2	0.0789992	3MK7	0.2833149
ν2	0.0792016	M8	0.3220456

A feature of the data is an unexplained shift around 1980 that we see in the M₂ amplitude (Fig. 3). This has been noted both by Ray (2006) and Müller (2011). Ray felt that because the amplitude change was evident in the Gulf of Maine and Bay of Fundy stations as well as Halifax, the process giving rise to the shift might not be just local to the Bay of Fundy and Gulf of Maine, although he did not see this in several other nearby American stations and could not explain it. Müller (2011) carried out an extensive analysis of many North Atlantic stations and did see a widespread negative offset or at least a negative trend in many North Atlantic M₂ amplitudes and speculated that this may be associated with global warming. Woodworth (2010) itemized changes in tidal components in many parts of the world but could not ascribe clear processes to these changes.

We carried out a regression analysis of estimated M₂ amplitudes allowing for shifts in slope and intercept beyond 1981. We chose the breakpoint in the M₂ amplitude (before 1 January 1982, after 31 December 1981) by minimizing the combined residual sum of squares over all ports, allowing the location of the breakpoint to vary between 1960 and 1990.

The following model was fit to Y_t , the M₂ amplitude at year t, where t is the year, ω denotes the frequency of the 18.599 year nodal cycle, and Z = 1 for 1982 onwards and 0 otherwise.

The harmonic terms of frequency 2ω are included to allow for possible asymmetry in the nodal variation. Apart from the nodal variation terms, the mean M₂ amplitude at year t is β₀ + β₁ t for years t ≤ 1981 and β₀ + β₂ + (β₁ + β₃)t for years t ≥ 1982. The slope of the amplitude trend is β₁ prior to 1982 and β₁ + β₃ from 1982 onwards, with β₂ representing a shift in mean level after 1981.

For each station, subset models were fitted and a best fitting model was chosen to minimize the Akaike Information Criterion (AIC) (Akaike, 1974). The minimum AIC models included a positive trend β₁ before 1982 at all locations but Halifax where a trend term prior to 1982 was not included. At Boston, Portland and Saint John, a negative level shift after 1981 (β₂ < 0) was included, and this was of particularly large magnitude at Saint John. At Eastport and Halifax there was a reduced trend post-break (β₃  < 0) while this was positive at Saint John. Nodal variation terms cos(ωt) and/or sin(ωt) were included at all locations, and a harmonic term sin(2ωt) was included at Portland, Eastport and Saint John.

Autocorrelation and partial autocorrelation plots of the residual deviations of the best fitting models suggest a variety of structures (white noise at Portland; approximate white noise at Boston and Eastport (first-order autoregressive with coefficients approximately 0.3); second-order autoregressive at Halifax; and perhaps a higher order mixed autoregressive moving average model at Saint John), although we did not attempt formal residual model identification. In general, the residuals appeared to be approximately normally distributed.

The results are summarized in Fig. 3 and Table 2. The 18.6 year nodal cycle is clear, and there is evidence of a reduction in the M₂ amplitude trend after 1981. The point and interval estimates in Table 2 are confidence intervals on the slope in M₂ amplitude prior to 1982 (β₁), from 1982 onwards (β₃), and for a regression model without a trend change at 1981/82, estimated in the full model (Eq. (1)). The relatively wide confidence intervals for the trend after 1981 are caused by the relatively short length of the record after that point.

Table 2. Estimated trends in Z₀ and M₂ tidal constituents with associated confidence intervals. The break in Z₀ trend is at 1979/1980. The break in M₂ amplitude and phase is at 1981/1982. Amplitude trends are in centimetres per century and phase trends are in degrees per century. The triples consist of lower confidence bound, point estimate of trend and upper confidence bound.

	Pre-break	Post-break	Single slope
Z0
Boston	(27.23, 32.34, 37.46)	(20.91, 35.76, 50.62)	(22.95, 26.03, 29.11)
Portland	(17.89, 22.36, 26.83)	(−5.78, 10.64, 27.06)	(14.51, 17.31, 20.10)
Eastport	(25.41, 31.35, 37.30)	(4.13, 17.00, 29.86)	(16.07, 19.51, 22.96)
Saint John	(21.12, 25.84, 30.57)	(−6.53, 13.05, 32.63)	(18.6, 21.7, 24.80)
Halifax	(34.86, 38.81, 42.76)	(12.59, 23.63, 34.66)	(29.12, 31.59, 34.06)
M2 amplitude
Boston	(2.90, 4.32, 5.74)	(0.79, 5.63, 10.46)	(1.95, 2.83, 3.70)
Portland	(5.66, 7.43, 9.21)	(0.59, 8.15, 15.71)	(5.13, 6.24, 7.35)
Eastport	(8.74, 12.4, 16.07)	(0.88, 10.26, 19.64)	(6.02, 8.05, 10.07)
Saint John	(4.80, 7.95, 11.11)	(22.34, 37.12, 51.9)	(2.8, 5.27, 7.73)
Halifax	(−1.00, −0.24, 0.51)	(−3.37, −0.89, 1.59)	(−2.07, −1.57, −1.06)
M2 phase
Boston	(−4.8, −3.5, −2.2)	(−6.1, −1.7, 2.7)	(−2.5, −1.7, −0.8)
Portland	(−3.3, −1.9, −0.4)	(−6.3, −0.2, 5.9)	(−1.8, −0.9, 0.0)
Eastport	(−2.5, −0.4, 1.7)	(−4.5, 0.9, 6.3)	(−1.6, −0.5, 0.5)
Saint John	(−2.3, −1.3, −0.3)	(−6.3, −1.6, 3.0)	(−0.2, 0.6, 1.4)
Halifax	(−1.6, 0.0, 1.5)	(−10.3, −5.3, −0.2)	(0.3, 1.2, 2.2)

The confidence intervals in Fig. 3 and elsewhere are 99% intervals. It follows from Bonferroni's inequality (see Miller, 1981) that this will provide a simultaneous coverage probability of at least 0.95 for the intervals at five stations. We did not attempt to adjust the confidence intervals to account for correlated errors, as in Ray (2009). No such adjustment was necessary for Portland where the residuals appear to be white noise, and any adjustment at Boston and Eastport would be fairly small as the residuals are approximately white noise. It would be challenging to find the appropriate adjustment for Saint John, where the residuals have a complex dependence structure.

There seem to be some peculiarities in the results that are explained by the inclusion of the mean level shift parameter β₂. For example, for the Saint John M₂ amplitude the overall slope is less than either of the pre-break or post-break slopes. This is a result of a large negative estimated shift in mean level beyond the break. As an aside, we believe that the post-break M₂ slope for Saint John may be inflated because of the aforementioned problems beginning with the relocation of the tide gauge.

The trend and other estimates are sensitive to the period of the data being analyzed. For M₂ at Saint John, Godin's (1992) analysis using data from 1947–79 gave an estimated trend of 14.5 cm cy⁻¹, which is substantially higher than our estimate of 7.95 cm cy⁻¹ before 1982. The difference reflects the sensitivity of the analysis to the data window. For example, our estimate of trend was 12.2 cm cy⁻¹ when we fit our model to the 1947–79 data window. The remaining difference (14.5–12.2) is likely caused by different methods of analysis. Where we model the nodal cycle explicitly, Godin fit a linear model to 19-year running means, the latter used to average over the nodal period.

Ray's (2009) estimates of the M₂ amplitude trend using data from 1935–2005 were −2.4 cm cy⁻¹ (Halifax), 7.5 cm cy⁻¹ (Eastport), 6.2 cm cy⁻¹ (Portland) and 2.7 cm cy⁻¹ (Boston), which are in fairly close agreement with our single slope estimates in Table 3.

Table 3. Percentage change per century in amplitude. Triples consist of lower confidence bound, percentage change per century and upper confidence bound. Pre-1982 is percentage change per century pre-1982, relative to the 1970 predicted mean amplitude. Post-1981 is percentage change per century post-1981, relative to the 1982 predicted mean amplitude. Single slope is percentage change per century relative to the 1970 predicted mean amplitude.

	Pre-1982	Post-1981	Single slope
M2
Boston	(2.12, 3.15, 4.18)	(0.56, 4.12, 7.69)	(1.43, 2.07, 2.71)
Portland	(4.18, 5.47, 6.75)	(0.39, 6.00, 11.62)	(3.79, 4.6, 5.41)
Eastport	(3.31, 4.69, 6.06)	(0.32, 3.89, 7.46)	(2.28, 3.05, 3.82)
Saint John	(1.58, 2.61, 3.64)	(7.40, 12.43, 17.45)	(0.93, 1.74, 2.55)
Halifax	(−1.57, −0.38, 0.81)	(−5.38, −1.42, 2.54)	(−3.29, −2.50, −1.7)

For fitting the M₂ phase, the same break point (1981/82) was used as for the M₂ amplitude, with the best models being again selected using the AIC criterion. The selected models for Boston and Portland included a negative trend prior to 1982 (β₁ < 0) with a reduced rate of decrease after 1981 (β₃ > 0). Saint John showed a negative trend with no change in trend after 1981 but with a positive mean level shift after 1981 (β₂ > 0). The model for Eastport included a negative mean level shift but no trend. The model for Halifax included a positive mean level shift with a negative trend after 1981, but no trend before 1982. All stations except Boston had one or more nodal modulation terms included.

Ray's (2009) estimates of the M₂ phase trend are 1.8° cy⁻¹ (Halifax), −0.1° cy⁻¹ (Eastport), 1.6° cy⁻¹ (Portland) and −0.7° cy⁻¹ (Boston). Godin (1992) estimates 0.72° cy ^{− 1} at Saint John. These phase estimates differ considerably from our estimates in Table 2. We suspect that this reflects either the different data windows used or the method of processing.

We also carried out a trend analysis of estimated yearly mean sea level Z₀, using the same procedures as for M₂, and found evidence of a reduced trend in mean sea level beginning in 1980 (Fig. 2, Table 2), with changes in intercept and slope being highly significant for Halifax and Eastport and moderately significant at Portland.

Although the inclusion of trend changes and level shifts were supported by the statistical analysis, there is no a priori reason to expect linear trends. A smoother change in mean sea level or M₂ amplitude might be more appropriate and we are currently investigating this with a more thorough examination of hourly sea level records. Godin (1992), for example, considered using a quadratic model.

The rates of change of annual mean sea level are not spatially uniform (Figs 2, 4 and Table 2). The differences are consistent with Peltier's (2002) results indicating a minimum rate of change in the western Gulf of Maine (Portland) and increasing rates along the coast in both directions, particularly closer to the deep ocean (Halifax).

The analyzed rate of change in M₂ tidal amplitude (Fig. 5, Table 2) is increased for stations with higher tidal amplitude. This is further emphasized in Fig. 6 and Table 3.

Fig. 5 Change per century in M₂ tidal amplitude as determined from the data analysis. The post-break Saint John line is off the scale (see Table 2 and Section 2).

Fig. 6 Percentage change per century in M₂ tidal amplitude. Note the proximity of the pre-break results from the Gulf of Maine and Bay of Fundy and the separation from Halifax. We believe this to be indicative of a phenomenon within the Gulf that is different from the adjacent surrounding ocean. Confidence bounds truncated by axis limits can be found in Table 2.

The relative rate of M₂ amplitude change, expressed as a fraction of the amplitude at year t is given by β₁/(β₀ + β₁ t) for years prior to the break and by (β₁ + β₃)/(β₀ + β₁ t + β₂ + β₃ t) after the break. Table 3 shows these rates of change relative to the 1970 and 1982 mean amplitudes, adjusted to percentage per century. Confidence limits on the relative rates of change were calculated using the delta method (van der Vaart, 2000).

Valid comparisons of rate of change with other published values require equivalent denominators. Müller et al. (2010) calculated the estimated rate of M₂ amplitude change relative to the 1970 mean as −3.0 (Halifax), 0.9 (Saint John), 2.8 (Eastport), 4.4 (Portland) and 1.9 (Boston). These are in reasonably close agreement with the single slope estimates in Table 3. Ray (2009) computes a similar quantity but relative to the 1935–2005 average. Müller's (2011) rates of change are the mean difference in M₂ means for the periods 1930–80 and 1985–95 relative to the 1930–80 mean, so not directly comparable to our estimates.

We note that in the pre-break period, those stations within the system are grouped together at a 2.8%–5.4% increase in amplitude per century. These are well removed from the Halifax amplitude change which is negative but near zero. This suggests that the increase in M₂ tidal amplitudes is a Bay of Fundy–Gulf of Maine phenomenon, not seen in the adjacent waters. The analyzed rate of change in M₂ tidal phase is seen in Fig. 7 and Table 2. The rate of change in pre-1982 M₂ amplitude is low compared to Godin (1992) but can be explained by the different data windows (see discussion in this section). The changes found in phase were near zero, as in Godin (1992), but different in sign.

Fig. 7 Change per century in M₂ tidal phase.

Given the extra uncertainty in the annual estimates beginning around 1982, we are more comfortable using data prior to this break to establish local pre-climate change trends in sea level and M₂ amplitude.

3 The modelling methodology

The model used in this study is referred to as the Toulouse Unstructured Grid Ocean model (T-UGOm; Pairaud et al., 2008). This is a flexible, fully non-linear, three-dimensional, finite-element model in spherical-polar coordinates. We run it in a two-dimensional barotropic mode, with P1P1 elements (velocity components and elevation computation points coincident at the triangle vertices) solving the wave equation as in Lynch et al. (1996) with the wetting and drying of inter-tidal areas as in Greenberg et al. (2005).

Earlier studies (Greenberg, 1975, 1979) have shown that to simulate changes in the tidal regime adequately in this area, there is a need to:

1.	resolve the upper Bay,
2.	include drying inter-tidal areas,
3.	cover the full domain of the Bay of Fundy–Gulf of Maine resonance, and
4.	include the full non-linear representation of the governing equations.

This variable-resolution finite-element model is ideally suited to this task. The upper reaches of the Bay are well resolved. The model mesh covers the full resonant domain plus the adjacent continental shelf and deep sea (Fig. 8). It has 14,070 nodes, 26,251 elements and a maximum, mean and minimum node separation of 53.3 km, 5.4 km and 170 m, respectively.

Fig. 8 The finite element model mesh used in tidal regime computations. There is higher resolution in shallow areas, in areas with steep gradients and in the Upper Bay of Fundy. The bathymetry colour scale is in metres.

In the absence of meteorological forcing, the sea level (mean and tidal) within the Gulf of Maine and Bay of Fundy will be driven predominantly by the sea level at the open boundary. Within this domain the direct influences on tides from gravity, loading and self attraction will be minimal. Climate-induced local steric variations are thought to have little effect on tidal elevation. The different model-specified components of sea level driving are described in Section 5. The sea level changes not associated with modern climate change are prescribed by adjusting the model bathymetry according to our spatial model (Section 5d); the steric and land-ice-generated sea level changes are prescribed as a spatially constant addition to this (Sections 5b and 5c). The same M₂ tidal amplitudes and phases are prescribed on the open boundary for all sea level configurations (see Section 5a).

4 Global average sea level change

Global average sea level rise, from 1990 to 2100, was projected to lie between 0.11 and 0.77 m (Church et al., 2001). These projections were derived from 14 atmosphere-ocean global climate models (AOGCMs) driven with the IS92a emissions scenario (including sulphate aerosols). Thermal expansion of sea water was the main contributor to global sea level rise, with lesser contributions from glaciers and the Greenland and Antarctic ice sheets.

Adoption of the Special Report on Emissions Scenarios (SRES) (Nakicenovic et al., 2000) led to an update of these estimates to between 0.09 and 0.88 m for 1990 to 2100 (Fig. 9), with thermal expansion again being the major contributor. The central value is 0.48 m, which corresponds to an average rate of about two to four times the rate for the twentieth century (Church et al., 2001). Although models agree on the qualitative conclusion that the range of regional variation in sea level change is substantial compared to global average sea level rise, confidence in the regional distribution of sea level change from AOGCMs is low because there is little similarity between models (Church et al., 2001). However, nearly all models project a greater than average rise in the Arctic Ocean and a less than average rise in the Southern Ocean.

Fig. 9 Global average sea level rise 1990 to 2100 for the SRES scenarios. Thermal expansion and land-ice changes were calculated using a simple climate model calibrated separately for each of seven AOGCMs; contributions from changes in permafrost, the effect of sediment deposition and the long-term adjustment of the ice sheets to past climate change were added. Each of the six lines appearing in the key is the average of AOGCMs for one of the six illustrative scenarios. The region in dark shading shows the range of the average of AOGCMs for all thirty-five SRES scenarios. The region in light shading shows the range of all AOGCMs for all thirty-five scenarios. The region delimited by the outermost lines shows the range of all AOGCMs and scenarios including uncertainty in land-ice changes, permafrost changes and sediment deposition. Note that this range does not allow for uncertainty relating to ice-dynamic changes in the West Antarctic ice sheet (Figure 11.12 Church et al. (2001)).

In the most recent Fourth Assessment Report (Meehl et al., 2007), global average sea level is projected to rise from the present (1980–99) to the end of this century (2090–99) under the SRES B1 scenario by 0.18 to 0.38 m, B2 by 0.20 to 0.43 m, A1B by 0.21 to 0.48 m, A1T by 0.20 to 0.45 m, A2 by 0.23 to 0.51 m, and A1FI by 0.26 to 0.59 m (Fig. 10). For each emissions scenario, the midpoint of the range is within 10% of the Intergovernmental Panel on Climate Change (IPCC) Third Assessment Report's (Church et al., 2001) model average for 2090–99. The ranges are narrower than in Church et al. (2001) mainly because of improved information about some uncertainties in the projected contributions. Thermal expansion contributes 70% to 75% of the central estimate in these projections for all scenarios.

Fig. 10 Projections and uncertainties (5% to 95% ranges) of global average sea level rise and its components from 2090 to 2099 (relative to 1980 to 1999) for the six SRES marker scenarios (Figure 10.33 Meehl et al. (2007)).

Because of ocean density and circulation changes, the distribution of sea level rise during the twenty-first century is projected to have substantial geographical variability. Figure 11 illustrates a distribution of local sea level change based on an ensemble mean of 14 AOGCMs. A smaller than average sea level rise can be seen in the Southern Ocean and a larger than average rise in the Arctic, which is in keeping with results reported in Church et al. (2001). A narrow band of pronounced sea level rise stretching across the southern Atlantic and Indian oceans is also apparent. Large deviations among models, however, make estimates of sea level distribution across the Caribbean, Indian and Pacific oceans uncertain.

Fig. 11 Local sea level change (m) caused by ocean density and circulation change relative to the global average (i.e., positive values indicate greater local sea level change than global) during the twenty-first century, calculated as the difference between averages for 2070 to 2099 and 1980 to 1999, as an ensemble mean over 16 AOGCMs forced with the SRES A1B scenario. Stippling denotes regions where the magnitude of the multi-model ensemble mean divided by the multi-model standard deviation exceeds 1.0 (Figure 10.32 Meehl et al. (2007)).

5 Components of sea level change

In this section, we describe how we specify the model inputs corresponding to the different physical processes contributing to sea level change.

a Tide

The amplitudes and phases of the M₂ tide are needed along the open (sea) boundary to drive the tides in the interior of the model. We only consider the M₂ constituent because it dominates the tides in the Bay of Fundy and Gulf of Maine and can be considered to represent mean conditions. The initial tidal boundary specification was from BIO (2008) which was based on Dupont et al. (2005). Calibration adjustments to these boundaries and friction parameters were made to obtain agreement with the present-day tidal regime (see Section 6).

We assume that the tides on this open boundary will not change with rising sea level. A simulation of global tides with a 5 m increase in sea level (Green, 2010) found minimal changes in the deep sea tides adjacent to the Gulf of Maine. Further, small changes in the North Atlantic M₂ tide are somewhat muted in the Gulf of Maine as increased tides are damped by increased friction, and decreased tides have lower self damping (Ku et al., 1985).

In each run the model is started from rest. The open boundary tide is ramped up from zero to full tide over 12 hours. The model is run for 16 tidal cycles to get a steady (repeatable) tidal wave. For each run, the last cycle is analyzed to get the amplitude and phase.

b Steric changes

The expansion of sea water caused by the warming associated with climate change will lead to a rise in sea level. Thus, to get approximate sea level changes for model input, we have obtained multiple results from different coupled atmosphere–ocean global climate models: the Canadian Climate Centre for Climate Modelling and Analysis CGCM1 model and the Hadley Centre HadCM2 and HadCM3 models. Although these runs were made using the IS92a emissions scenario (with and without sulphate aerosols), the global average sea level rise projections from these runs are still in keeping with the most recent projections given in the IPCC Fourth Assessment Report (Meehl et al., 2007).

Average sea level was computed over the area delimited by 38.24°–43.83° N, 61.875°–73.125° W for each AOGCM run. Three time periods were averaged: 1961–90 (the baseline climate), 2040–69 (referred to as 2055) and 2070–99 (referred to as 2085). For the CGCM1 model, the sea level values available were changed with respect to the model's control run (an experiment where no changes in forcing were imposed). For all models, the sea level rise was obtained by subtracting the baseline climate sea levels from the 2055 and 2085 levels. For the CGCM1 model, this assumes that the control run sea level values are similar in each time period. There was a wide range in predicted sea level changes: 0.16–0.28 m for 2055 and 0.27–0.42 m for 2085 (Table 4).

Table 4. Changes applied to model bathymetry (in metres) to account for cryospheric and steric effects on sea level. The ice changes are as derived in Section 5c, with values for 2055 and 2085 being obtained by linear interpolation. The steric sea level changes for 2055 and 2085 are taken from model runs analyzed for this study (Section 5b). The 2100 steric values were extrapolated from those computations.

	2055	2085	2100
Ice
Minimum	−0.02	−0.03	−0.04
Central value	0.09	0.14	0.17
Maximum	0.19	0.29	0.34
Steric
Minimum	0.16	0.27	0.33
Mean	0.21	0.28	0.41
Maximum	0.28	0.42	0.49

c Ice

Projections of sea level change as a result of variations in land ice are made in the IPCC Fourth Assessment Report (Meehl et al., 2007). They have separated components from Greenland and Antarctic ice sheets and the glaciers and ice caps in the rest of the world. Although the melting of all Greenland and Antarctic ice would raise sea level by close to 70 m (Huybrechts et al., 2004), it is expected that in the twenty-first century these two regions will only make a small contribution to changing sea level. The Antarctic ice sheet is predicted to grow in the twenty-first century, diminishing the extent of sea level rise. Initially, the world's non-polar melting ice caps and glaciers will give rise to the major ice-derived sea level effects.

There is, however, little confidence in the predictions of sea level change caused by changing land ice. In Section 10.6.5 of Meehl et al. (2007) they assert: “… current scientific understanding leaves poorly known uncertainties in the methods used to make projections …”. The problem is that flows from polar ice sheets have recently accelerated, and it is possible that this recent trend could continue, stop, or further increase. Meehl et al. (2007) continue: “… we cannot assess the likelihood of any of these three alternatives … The state of understanding prevents a best estimate from being made.”

Given the uncertainties above, we have constructed a mid-value of the land-ice contribution to sea level rise in the twenty-first century by averaging over all scenarios the central values of “Glaciers and ice caps” in Fig. 10 [Figure 10.33 in Meehl et al. (2007)] and adding one-half the contribution from “Scaled-up ice sheet dynamical imbalance” similarly computed. To find an upper bound for our computations, we add the maximum of the range for any scenario of the “Glaciers and ice caps” plotted variable and the maximum of the range for any scenario of the “Scaled-up ice sheet dynamical imbalance.” To compute a lower bound for any scenario we summed the lower limit of the “Glaciers and ice caps” and the Greenland and Antarctic “ice sheet dynamical imbalance.” These extremes, upper and lower bounds, were all found in scenario A1F1. This procedure gives us the lower bound, mid-value and upper bound for sea level change caused by land-ice evolution from 2000 to 2100 of −0.04 m, 0.17 m and 0.34 m, respectively. Values for 2055 and 2085 were linearly interpolated assuming zero for the year 2000 (Table 4).

d Pre-climate change sea level trend

As described in Section 1, the Bay of Fundy and Gulf of Maine are still descending relative to mean sea level because of the glacial melting after the last ice age. This is one component of the existing trend in sea level. If the changes were uniform across the region, we could equivalently either change Z₀ on the model open boundary or increase all depths by the same amount. But because the existing modelled and historically observed sea level changes are spatially varying (Shaw et al., 1997; Peltier, 2002, 2004) as are the trends in our analyses (Figs 2, 4 and Table 2), we need a method to account for the variation.

In the forerunner of this work (Greenberg, 2001) experimentation was carried out with simple depth changes applied uniformly over different areas. Adding a constant 0.25 m to the model bathymetry (an approximation from the Z₀ analysis similar to that in our Table 2) produced changes in tides well below those seen in the observations. Earlier modelling of historical sea level (Scott and Greenberg, 1983) indicated that increasing depths on the outer shelf had a strong effect on tidal amplitudes. Studies of coastal and shelf subsidence (Shaw et al., 1997; Peltier, 2002, 2004) indicate that the tectonic changes are more pronounced offshore, but these models do not give a detailed picture. Considering this, a better match to observed tidal change was obtained by adding 0.5 m to the topography seaward of a line between Yarmouth and Cape Cod, on top of the 0.25 m uniformly applied in the initial run.

In this paper, we construct our own depth adjustment model (Fig. 12) to reproduce the observed tidal amplitude changes. We use points from analyses of long-term sea level observations and vary the offshore depths to get optimal agreement. With that in mind, for shore stations, we used our analyzed values (metres per century) for Boston (0.323), Portland (0.224), Saint John (0.258) and Halifax (0.388) (Table 2). These were supplemented with other values found in the literature for Woods Hole (0.261; Zervas, 2009), Newport (0.258; Zervas, 2009) and Charlottetown (0.32; Forbes et al., 2008). A value of 0.4 m was prescribed along the deep sea open boundary. We adjusted the value along the shelf break to calibrate for good agreement with our analyses of the changing M₂ tides. We found that what we consider to be a large value, 2.5 m, was needed to give a good match with our data analysis. Along the Middle Atlantic Bight cross-shelf boundary, the Cabot Strait boundary and the land points, values were set to be consistent with coastal stations and the specified deep sea and shelf break points. Our aim was to mimic the form of tectonic models in which the sea level change increased to about the shelf edge and decreased in the deep sea, but note that our value at the shelf break is much greater than those seen in tectonic models.

Fig. 12 Our estimated pre-climate change rate of increase in sea level in metres per century. Named ports were specified from our analyses (Boston, Portland, Saint John and Halifax) and other publications (Newport, Woods Hole and Charlottetown). The values along the shelf break were set to 2.5 m during calibration to obtain agreement with the analyzed increase in M₂ amplitude. Values on the deep sea open boundary were set to 0.4 metres per century. Other boundary values, needed to encompass the model domain, were specified to be consistent with nearby interior points. The increase in model bathymetry was linearly interpolated on the triangles of this mesh.

6 Model calibration and experiments

Using the model and grid as described in Section 3, calibration adjustments were made to tidal boundaries and frictional parameters to achieve a good fit to the existing data. The starting point was the open boundary used in the Dupont et al. (2005) outer model. Using trial and error, the model was adjusted by

1.	setting the quadratic friction coefficient to 2.4 × 10−3;
2.	scaling the initially prescribed open boundary tidal elevations so that the average value of the ratio between the modelled and observed amplitudes is 1; and
3.	adjusting the phase to minimize the error vector, defined as the distance in polar coordinates between the modelled and observed amplitude and phase (as in Fig. 13).

Fig. 13 Polar plot comparison between the calibrated model and the observations (red ▴). Radial variations correspond to differences in amplitude. Angular variations correspond to differences in phase.

This process involved scaling the amplitudes by less than 1% and changing the phase by approximately 2.5°. The agreement in tidal amplitude is very good in this macro-tidal environment. The average length of the error vector is 0.117 m. The average difference between modelled and observed amplitude and phase components is +0.006 m and −1.3°, respectively. The average absolute value of the difference between amplitude and phase components is 0.036 m and 3.4°.

Using bathymetry and boundary sea level changes as given in Section 5 (see Table 4 and Figure 12) and applied as described in Section 3, we have carried out several modelling experiments. The results from selected model runs are given in Table 5. Reference is made to the following case studies:

•	pre-climate change: This is the 100 year pre-climate change simulation described in Section 5d and in Fig. 12 where we increased sea level using our interpreted map. The comparison of the model with observed tidal change is shown in Figs 14 and 15. Influences on sea level from steric and melting land ice related to climate change are not included. This is the second stage of calibration, run to determine the pre-climate change depth change at the shelf break that is in good agreement with our analyzed tidal changes.
•	2055: This is the mean sea level rise value from the coupled atmosphere–ocean global climate models: +0.21 m and the linearly interpolated ice and pre-climate change sea level increase.
•	2085: This is the mean sea level rise value from the coupled atmosphere–ocean global climate models: +0.28 m and the linearly interpolated ice and pre-climate change sea level increase.
•	2100: This is the linearly extrapolated mean sea level rise values from the coupled atmosphere–ocean global climate models: +0.41 m and the total ice and pre-climate change increase.

Fig. 14 Polar plot comparison between the calibrated model present state (red ▴) and the model predicted analyzed change (green line) and the existing change determined from the pre-break data analysis (blue line - see also Fig. 15). Radial variations correspond to differences in amplitude in metres. Angular variations correspond to differences in phase.

Fig. 15 Expanded scale for data analysis of the stations in Fig. 14. The stations shown, in order of increasing amplitude (radius), are Halifax, Portland (right), Boston (left), Eastport and Saint John.

Fig. 16 The increase in the height of tidal high water from the sum of all factors as described in the text and the tables. PCC refers to the pre-climate change trends computed in Section 2. The error bar gives the maximum and minimum variation on this based on the literature for the ice contribution, the variability in the climate model results, and the tide computed taking into consideration these extremes.

TABLE 5. The predicted increase in tidal high water in metres at selected locations for the years 2055, 2085 and 2100. The values reflect tectonic influence, ice melt, steric effects and changing M₂ tides. Maximum and minimum values are obtained from tidal runs using our estimated existing sea level trend and the extreme values of the ice and steric inputs.

Time	2055			2085			2100
	min	mid	max	min	mid	max	min	mid	max
Boston	0.35	0.51	0.69	0.57	0.75	1.05	0.68	0.98	1.23
Portland	0.29	0.46	0.63	0.48	0.66	0.96	0.57	0.87	1.13
Saint John	0.33	0.50	0.68	0.54	0.73	1.04	0.64	0.95	1.21
Cobequid Bay	0.41	0.60	0.79	0.66	0.87	1.20	0.79	1.12	1.40
Yarmouth	0.37	0.54	0.71	0.60	0.78	1.08	0.71	1.01	1.27
Halifax	0.36	0.58	0.69	0.57	0.76	1.04	0.68	0.97	1.22

Upper and lower values of change in M₂ tidal elevation are computed for the uncertainties in ice and steric values (minima and maxima in Table 4). The effects on M₂ tidal elevation of uncertainties in our mean sea level analysis and interpretation are not evaluated.

7 Discussion

We have been looking at factors that will change predicted high water in the Bay of Fundy and Gulf of Maine over the coming century. This has included evidence of mean sea level change and change in the amplitude of the tides. The past increase in the tides is seen to be a process confined to this system. Our modelling has indicated that this tidal process was related to increasing sea level which we feel could be, at least in part, tectonic in origin from the post-glacial rebound of the last ice age but probably includes other unknown processes. Whether tectonic in origin or not, we have assumed that the trends that are found in our analysis of pre-break data are independent of modern climate change and will continue at the same rate over the coming century. We are able to reproduce the magnitude of these trends in tidal amplitude reasonably well by adjusting depths in a hydrodynamic model of the area. We then predict future trends in tidal amplitude by adding the predicted sea level changes resulting from steric effects (the expansion of sea water as it warms) and the melting of glaciers and ice caps.

We have presented results exclusively for the M₂ tidal constituent. The changes from analysis of the other principal constituents (not shown) had minimal significance. Within the Gulf of Maine and Bay of Fundy, changes in S₂ amplitude were small and negative (typically around −0.8 cm cy⁻¹ prior to 1973, except at Portland +0.8 cm cy⁻¹ and generally with an increased negative slope after 1972). Changes in N₂ amplitude were small and positive pre-break at 1990/91, except at Halifax where the pre-break slope was negative. Only at Portland was the pre-break slope significantly different from zero. Because of the short length of the post-break records, no slopes post-break were significantly different from zero. We think this fits well with our understanding that the primary resonant period of the system is between 12.5 and 13.3 hours. The dominant M₂ tide has a strong damping effect on other constituents but is minimally affected by them. Thus, M₂ could be damping both S₂ and N₂ amplitudes. With the dominant resonant period being reduced with increasing depth, the resonance with the N₂ (period ≈12.5 hours) will be increased, leading to a small net increase in amplitude. Changes in diurnal constituents O₁ and K₁ were small and not even consistent in sign.

Our estimates of rates of M₂ amplitude change, when using comparable data windows, match well with the results of Godin (1992), Ray (2009) and Müller (2011). Because of issues with the data and/or large error bars in our models for the last two decades of the twentieth century, we have based our work on the analysis of the earlier data.

8 Conclusions

Even in the absence of climate change, the Bay of Fundy–Gulf of Maine region can expect an increase in tidal high water of order 0.3 m over the next century (Fig. 16). Combined with the influence of climate change, high water in the Bay of Fundy is predicted to rise on the order of 0.5 m over the next 50 years, and on the order of 1 m by the end of the century.

We have been looking only at mean sea level and tide to the exclusion of meteorologically driven storm surges, inverse barometer sea level changes, surface waves and swell. The latter will assume added significance in a new regime of higher tidal peaks. In an area such as this, where flooding levels are so close to the present astronomical higher high water (Parkes et al., 1997), any increase in high water will be significant. Given our predicted increases in high water that include existing trends, climate change effects and modified tides as computed in this study, there will be significant implications for planning and land use in Bay of Fundy and Gulf of Maine coastal communities.

Acknowledgements

This project benefited greatly from discussions with Malte Müller, University of Victoria; Bill Crawford and Mike Foreman, Institute of Ocean Sciences, Fisheries and Oceans Canada, Patricia Bay, British Columbia; Keith Thompson, Dalhousie University, Halifax; and Charles O'Reilly and Phillip MacAulay of the Canadian Hydrographic Service, Bedford Institute of Oceanography. Dr. Jason Lowe of the UK Hadley Centre for Climate Prediction and Research kindly sent the raw data for the Hadley Centre model. The initial work on this project was largely funded by Environment Canada's Climate Change Action Fund.